text:
## Default algorithm

# Allows A === nothing as a stand-in for dense matrix
function defaultalg(A,b)
    if A isa DiffEqArrayOperator
        A = A.A
    end

    # Special case on Arrays: avoid BLAS for RecursiveFactorization.jl when
    # it makes sense according to the benchmarks, which is dependent on
    # whether MKL or OpenBLAS is being used
    if (A === nothing && !isgpu(b)) || A isa Matrix
        if (A === nothing || eltype(A) <: Union{Float32,Float64,ComplexF32,ComplexF64}) &&
                    ArrayInterface.can_setindex(b)
            if length(b) <= 10
                alg = GenericLUFactorization()
            elseif (length(b) <= 100 || (isopenblas() && length(b) <= 500))
                alg = RFLUFactorization()
            else
                alg = LUFactorization()
            end
        else
            alg = LUFactorization()
        end

    # These few cases ensure the choice is optimal without the
    # dynamic dispatching of factorize
    elseif A isa Tridiagonal
        alg = GenericFactorization(;fact_alg=lu!)
    elseif A isa SymTridiagonal
        alg = GenericFactorization(;fact_alg=ldlt!)
    elseif A isa SparseMatrixCSC
        alg = UMFPACKFactorization()

    # This catches the cases where a factorization overload could exist
    # For example, BlockBandedMatrix
    elseif A !== nothing && ArrayInterface.isstructured(A)
        alg = GenericFactorization()

    # This catches the case where A is a CuMatrix
    # Which does not have LU fully defined
    elseif isgpu(A) || isgpu(b)
        alg = QRFactorization(false)

    # Not factorizable operator, default to only using A*x
    else
        alg = KrylovJL_GMRES()
    end
    alg
end

## Other dispatches are to decrease the dispatch cost

function SciMLBase.solve(cache::LinearCache, alg::Nothing,
                         args...; kwargs...)
    @unpack A = cache
    if A isa DiffEqArrayOperator
        A = A.A
    end

    # Special case on Arrays: avoid BLAS for RecursiveFactorization.jl when
    # it makes sense according to the benchmarks, which is dependent on
    # whether MKL or OpenBLAS is being used
    if A isa Matrix
        b = cache.b
        if (A === nothing || eltype(A) <: Union{Float32,Float64,ComplexF32,ComplexF64}) &&
                    ArrayInterface.can_setindex(b)
            if length(b) <= 10
                alg = GenericLUFactorization()
                SciMLBase.solve(cache, alg, args...; kwargs...)
            elseif (length(b) <= 100 || (isopenblas() && length(b) <= 500))
                alg = RFLUFactorization()
                SciMLBase.solve(cache, alg, args...; kwargs...)
            else
                alg = LUFactorization()
                SciMLBase.solve(cache, alg, args...; kwargs...)
            end
        else
            alg = LUFactorization()
            SciMLBase.solve(cache, alg, args...; kwargs...)
        end

    # These few cases ensure the choice is optimal without the
    # dynamic dispatching of factorize
    elseif A isa Tridiagonal
        alg = GenericFactorization(;fact_alg=lu!)
        SciMLBase.solve(cache, alg, args...; kwargs...)
    elseif A isa SymTridiagonal
        alg = GenericFactorization(;fact_alg=ldlt!)
        SciMLBase.solve(cache, alg, args...; kwargs...)
    elseif A isa SparseMatrixCSC
        alg = UMFPACKFactorization()
        SciMLBase.solve(cache, alg, args...; kwargs...)

    # This catches the cases where a factorization overload could exist
    # For example, BlockBandedMatrix
    elseif ArrayInterface.isstructured(A)
        alg = GenericFactorization()
        SciMLBase.solve(cache, alg, args...; kwargs...)

    # This catches the case where A is a CuMatrix
    # Which does not have LU fully defined
    elseif isgpu(A)
        alg = QRFactorization(false)
        SciMLBase.solve(cache, alg, args...; kwargs...)

    # Not factorizable operator, default to only using A*x
    # IterativeSolvers is faster on CPU but not GPU-compatible
    else
        alg = KrylovJL_GMRES()
        SciMLBase.solve(cache, alg, args...; kwargs...)
    end
end

function init_cacheval(alg::Nothing, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    if A isa DiffEqArrayOperator
        A = A.A
    end

    # Special case on Arrays: avoid BLAS for RecursiveFactorization.jl when
    # it makes sense according to the benchmarks, which is dependent on
    # whether MKL or OpenBLAS is being used
    if A isa Matrix
        if (A === nothing || eltype(A) <: Union{Float32,Float64,ComplexF32,ComplexF64}) &&
                    ArrayInterface.can_setindex(b)
            if length(b) <= 10
                alg = GenericLUFactorization()
                init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
            elseif (length(b) <= 100 || (isopenblas() && length(b) <= 500))
                alg = RFLUFactorization()
                init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
            else
                alg = LUFactorization()
                init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
            end
        else
            alg = LUFactorization()
            init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
        end

    # These few cases ensure the choice is optimal without the
    # dynamic dispatching of factorize
    elseif A isa Tridiagonal
        alg = GenericFactorization(;fact_alg=lu!)
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    elseif A isa SymTridiagonal
        alg = GenericFactorization(;fact_alg=ldlt!)
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    elseif A isa SparseMatrixCSC
        alg = UMFPACKFactorization()
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)

    # This catches the cases where a factorization overload could exist
    # For example, BlockBandedMatrix
    elseif ArrayInterface.isstructured(A)
        alg = GenericFactorization()
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)

    # This catches the case where A is a CuMatrix
    # Which does not have LU fully defined
    elseif isgpu(A)
        alg = QRFactorization(false)
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)

    # Not factorizable operator, default to only using A*x
    # IterativeSolvers is faster on CPU but not GPU-compatible
    else
        alg = KrylovJL_GMRES()
        init_cacheval(alg, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose)
    end
end


text:
"""
Integrator Block
"""

export IntegratorBlock

abstract type AbstractIntegratorBlock <: AbstractBlock end

mutable struct IntegratorBlock <: AbstractIntegratorBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}
    state::Symbol
    initialcondition::AbstractSignalLine
    saturationlimits::Vector{AbstractSignalLine}

    function IntegratorBlock(::Type{Tv} = Float64; initialcondition = 0, saturationlimits = []) where Tv
        @createblock new(Vector{InPort}(), Vector{OutPort}(), gensym(),
            signal(initialcondition, Tv),
            [signal(x, Tv) for x = saturationlimits]
            ) 1 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::IntegratorBlock)
    println(io, "IntegratorBlock()")
end

"""
to expr
"""

function _toexpr(blk::IntegratorBlock)
    return blk.state
end

"""
tosystemexpr

This is an expression to build system equations.

"""

function _tosystemexpr(blk::IntegratorBlock)
    saturationlimits = [_toexpr(x) for x = blk.saturationlimits]
    if length(saturationlimits) == 0
        _toexpr(blk.inport[1])
    else
        lower, upper = saturationlimits
        Expr(:if, Expr(:comparison, lower, :<=, blk.state, :<=, uupper), _toexpr(blk.inport[1]), 0)
    end
end

ext: jl


text:
@testset "NoiseFunction" begin

using DiffEqNoiseProcess, DiffEqBase, Test
f = (t) -> exp(t)

W = NoiseFunction(0.0,f)

dt = 0.1
calculate_step!(W,dt)

for i in 1:10
  accept_step!(W,dt)
end

@test W(W.curt + W.dt)[1] - W.curW == W.dW

prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

@test sol(W.curt + W.dt)[1] - W.curW == W.dW

f = (out,t) -> (out.=exp(t))
W = NoiseFunction(0.0,f,noise_prototype=rand(4))
prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

out1 = rand(4)
out2 = nothing
sol(out1,out2,0.5)

@test out1 == sol(0.5)[1]

end


text:
# We never need to permute indices on the CPU.
@inline   permute_index(solver_type, ::CPU, i, j, k, Nx, Ny, Nz) = i, j, k
@inline unpermute_index(solver_type, ::CPU, i, j, k, Nx, Ny, Nz) = i, j, k

# For why we use Base.unsafe_trunc instead of floor below see:
# https://github.com/CliMA/Oceananigans.jl/issues/828
# https://github.com/CliMA/Oceananigans.jl/pull/997

"""
    _permute_index(i, N)

Permute `i` such that, for example, `i ∈ 1:N` becomes

    [1, 2, 3, 4, 5, 6, 7, 8] -> [1, 8, 2, 7, 3, 6, 4, 5]

if `N=8`.
"""
@inline _permute_index(i, N)::Int = ifelse(isodd(i),
                                           Base.unsafe_trunc(Int, i/2) + 1,
                                           N - Base.unsafe_trunc(Int, (i-1)/2))

"""
    _unpermute_index(i, N)

Permute `i` in the opposite manner as `_permute_index`, such that,
for example, `i ∈ 1:N` becomes

   [1, 2, 3, 4, 5, 6, 7, 8] -> [1, 3, 5, 7, 8, 6, 4, 3]

if `N=8`.
"""
@inline _unpermute_index(i, N) = ifelse(i <= ceil(N/2), 2i-1, 2(N-i+1))

@inline   permute_index(::TriplyPeriodic, ::GPU, i, j, k, Nx, Ny, Nz) = i, j, k
@inline unpermute_index(::TriplyPeriodic, ::GPU, i, j, k, Nx, Ny, Nz) = i, j, k

@inline   permute_index(::HorizontallyPeriodic, ::GPU, i, j, k, Nx, Ny, Nz) = i, j,   _permute_index(k, Nz)
@inline unpermute_index(::HorizontallyPeriodic, ::GPU, i, j, k, Nx, Ny, Nz) = i, j, _unpermute_index(k, Nz)

@inline   permute_index(::Channel, ::GPU, i, j, k, Nx, Ny, Nz) = i,   _permute_index(j, Ny),   _permute_index(k, Nz)
@inline unpermute_index(::Channel, ::GPU, i, j, k, Nx, Ny, Nz) = i, _unpermute_index(j, Ny), _unpermute_index(k, Nz)


text:
################################################################################
#  Copyright 2021, Marta Vanin                                                 #
################################################################################
function run_case_study_E(nlsolver)

    season = "winter"
    time_step = 144
    elm = ["load"]
    pfs = [0.95]

    msr_path = joinpath(mktempdir(),"temp.csv")

    solver = _PMD.optimizer_with_attributes(nlsolver...)

    data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(1, 2),data_model=_PMD.ENGINEERING); #EU LV feeder
    _PMDSE.rm_enwl_transformer!(data)
    _PMDSE.reduce_enwl_lines_eng!(data)
    data["settings"]["sbase_default"] = 100.0 # set power_base to convert in per unit
    _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

    data = _PMD.transform_data_model(data);
    _PMDSE.reduce_single_phase_loadbuses!(data, exclude = [])

    pf_results = _PMD.solve_mc_pf(data, _PMD.ACPUPowerModel, solver)

    # Write measurements based on power flow
    v_pu = data["settings"]["vbases_default"]["1"]* data["settings"]["voltage_scale_factor"] # divider [V] to get the voltage in per units.
    v_max_err = 1.15
    σ_v = 1/3*v_max_err/v_pu

    p_pu = data["settings"]["sbase"]
    p_max_err = 0.01
    σ_p = 1/3*p_max_err/p_pu

    # Write measurements based on power flow
    σ_dict = Dict("load" => Dict("load" => σ_p,
                    "bus"  => σ_v),
                    "gen"  => Dict("gen" => σ_p/100,
                    "bus"  => σ_v/100)
                            )

    _PMDSE.write_measurements!(_PMD.ACPUPowerModel, data, pf_results, msr_path, σ = σ_dict)

    # Read-in measurement data and set initial values
    _PMDSE.add_measurements!(data, msr_path, actual_meas = true)

    data["se_settings"] = Dict{String,Any}("criterion" => "rwlav",
                                            "rescaler" => 100)

    _PMDSE.assign_start_to_variables!(data)
    _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

    ###############################################################
    ######                                                    #####
    ######      ADD BAD DATA AND RUN STATE ESTIMATION         #####
    ######                                                    #####
    ###############################################################

    bad_data_lav = deepcopy(data)
    _SEF.add_some_noise!(bad_data_lav)
    bad_data_lav["meas"]["43"]["dst"] = [_DST.Normal(0.019, σ_p)] # This is P_5 in the paper
    bad_data_lav["meas"]["53"]["dst"] = [_DST.Normal(0.01, σ_p)] # This is Q_27 in the paper
    bad_data_lav["meas"]["33"]["dst"] = [_DST.Normal(1.01, σ_v)] #This is U_8 in the paper
    bad_data_chi = deepcopy(bad_data_lav)
    bad_data_lnr = deepcopy(bad_data_lav)

    bad_data_lav["se_settings"] = Dict{String,Any}("criterion" => "rwlav",
    "rescaler" => 1)

    for (m, meas) in bad_data_lav["meas"]
        σ = [1.0, 1.0, 1.0]
        μ = _DST.mean.(meas["dst"])
        meas["dst"] = [_DST.Normal(μ[i], σ[i]) for i in 1:length(meas["dst"])]
    end

    deleted_data_points_lav = []
    for i in 1:4 # runs SE+chi square test+lav 4 times, three of which have a bad data point, the last one has no bad data (see paper)

        bad_data_n = [3, 2, 1, "no"]
        display("Now running SE with Chi-squares and LAV for the $i time. There are now $(bad_data_n[i]) bad data points.")

        ###############################################################
        ######                                                    #####
        ######                Chi-SQUARE ANALYSIS                 #####
        ######                                                    #####
        ###############################################################

        # NB: bad_data_chi and bad_data_lav are identical except for the SE settings: all weights are 1 in bad_data_lav, but that assumption won't work with chi square analysis
        se_result_bd = _PMDSE.solve_acp_red_mc_se(bad_data_chi, solver)

        excds_bd, obj_bd, trsh_bd = _PMDSE.exceeds_chi_squares_threshold(se_result_bd, bad_data_chi, suppress_display=true)

        if excds_bd
            display("The objective of ACP state estimation WITH bad data is $(obj_bd), and exceeds the Chi-square threshold of $trsh_bd. Bad data detected.")
        else
            display("The objective of ACP state estimation WITH bad data is $(obj_bd), and does not exceed the Chi-square threshold of $trsh_bd.")
        end

        ###############################################################
        ######                                                    #####
        ###### LAV ESTIMATION AND HIGHEST RESIDUAL IDENTIFICATION #####
        ######                                                    #####
        ###############################################################

        se_result_bd_lav = _PMDSE.solve_acp_red_mc_se(bad_data_lav, solver)
        residual_tuples = [(m, maximum(meas["res"])[1]) for (m, meas) in se_result_bd_lav["solution"]["meas"]]
        sorted_tuples = sort(residual_tuples, by = last, rev = true)

        delete!(bad_data_lav["meas"], first(sorted_tuples[1]))
        delete!(bad_data_chi["meas"], first(sorted_tuples[1]))
        push!(deleted_data_points_lav, first(sorted_tuples[1]))

        display("The largest residuals and their relative measurement number are the following:")
        display(sorted_tuples[1:4])

    end

    deleted_data_points_lnr = []
    for i in 1:4 # runs LNR 4 times, three of which have a bad data point, the last one has no bad data (see paper)

        bad_data_n = [3, 2, 1, "no"]
        display("Now running SE + LNR for the $i time. There are now $(bad_data_n[i]) bad data points.")
        se_result_bd = _PMDSE.solve_acp_red_mc_se(bad_data_lnr, solver)

        ###############################################################
        ######                                                    #####
        ######           LARGEST NORMALIZED RESIDUALS             #####
        ######                                                    #####
        ###############################################################

        max_meas_idx = maximum([parse(Int64, m) for (m,meas) in data["meas"]])

        _PMDSE.add_zib_virtual_meas!(bad_data_lnr, 1e-15, exclude = [2, 1, 30, 48, 36])
        for (_, load) in bad_data_lnr["load"]
            load["status"] = 1
            if !haskey(load, "connections") load["connections"] = [1,2,3] end
        end

        variable_dict = _PMDSE.build_variable_dictionary(bad_data_lnr)
        h_array = _PMDSE.build_measurement_function_array(bad_data_lnr, variable_dict)
        state_array = _PMDSE.build_state_array(se_result_bd, variable_dict)

        H = _PMDSE.build_H_matrix(h_array, state_array)
        R = _PMDSE.build_R_matrix(bad_data_lnr)
        G = _PMDSE.build_G_matrix(H, R)
        K = _PMDSE.build_K_matrix(H, G, R)
        S = _PMDSE.build_S_matrix(K)
        Ω = _PMDSE.build_omega_matrix(S, R)

        for (m,meas) in bad_data_lnr["meas"]
            if parse(Int64, m) > max_meas_idx delete!(bad_data_lnr["meas"], m) end
        end

        id_val, exc = _PMDSE.normalized_residuals(bad_data_lnr, se_result_bd, Ω)
        lnr_tuples = [(m, maximum(se_result_bd["solution"]["meas"][m]["nr"])[1]) for (m, meas) in bad_data_lnr["meas"]]
        sorted_lnr_tuples = sort(lnr_tuples, by = last, rev = true)

        push!(deleted_data_points_lnr, first(id_val))
        bad_data_lnr = deepcopy(data)
        bad_data_lnr["meas"]["43"]["dst"] = [_DST.Normal(0.019, σ_p)]
        bad_data_lnr["meas"]["53"]["dst"] = [_DST.Normal(0.01, σ_p)]
        bad_data_lnr["meas"]["33"]["dst"] = [_DST.Normal(1.01, σ_v)]
        for m in deleted_data_points_lnr
            delete!(bad_data_lnr["meas"], m)
        end

        display("The largest normalizes residuals and their relative measurement number are the following:")
        display(sorted_lnr_tuples[1:4])
        display("If bigger than 3, the first LNR is deleted.")

    end
end


text:
export A2CLearner

"""
    A2CLearner(;kwargs...)

# Keyword arguments

- `approximator`::[`ActorCritic`](@ref)
- `γ::Float32`, reward discount rate.
- `actor_loss_weight::Float32`
- `critic_loss_weight::Float32`
- `entropy_loss_weight::Float32`
- `update_freq::Int`, usually set to the same with the length of trajectory.
"""
Base.@kwdef mutable struct A2CLearner{A<:ActorCritic} <: AbstractLearner
    approximator::A
    γ::Float32
    max_grad_norm::Union{Nothing,Float32} = nothing
    actor_loss_weight::Float32
    critic_loss_weight::Float32
    entropy_loss_weight::Float32
    update_freq::Int
    update_step::Int = 0
    # for logging
    actor_loss::Float32 = 0.0f0
    critic_loss::Float32 = 0.0f0
    entropy_loss::Float32 = 0.0f0
    loss::Float32 = 0.0f0
    norm::Float32 = 0.0f0
end

Flux.functor(x::A2CLearner) = (app = x.approximator,), y -> @set x.approximator = y.app

function (learner::A2CLearner)(env::MultiThreadEnv)
    learner.approximator.actor(send_to_device(device(learner), state(env))) |> send_to_host
end

function (learner::A2CLearner)(env)
    s = state(env)
    s = Flux.unsqueeze(s, ndims(s) + 1)
    s = send_to_device(device(learner), s)
    learner.approximator.actor(s) |> vec |> send_to_host
end

function RLBase.update!(learner::A2CLearner, t::CircularArraySARTTrajectory)
    length(t) == 0 && return  # in the first update, only state & action is inserted into trajectory
    learner.update_step += 1
    if learner.update_step % learner.update_freq == 0
        _update!(learner, t)
    end
end

function _update!(learner::A2CLearner, t::CircularArraySARTTrajectory)
    n = length(t)

    AC = learner.approximator
    γ = learner.γ
    w₁ = learner.actor_loss_weight
    w₂ = learner.critic_loss_weight
    w₃ = learner.entropy_loss_weight
    to_device = x -> send_to_device(device(AC), x)

    states = select_last_dim(t[:state], 1:n)
    states_flattened = flatten_batch(states) |> to_device # (state_size..., n_thread * update_freq)

    actions = select_last_dim(t[:action], 1:n)
    actions = flatten_batch(actions)
    actions = CartesianIndex.(actions, 1:length(actions))

    next_state_values =
        t[:state] |> select_last_frame |> Array |> to_device |> AC.critic |> send_to_host

    gains =
        discount_rewards(
            t[:reward],
            γ;
            dims = 2,
            init = send_to_host(next_state_values),
            terminal = t[:terminal],
        ) |> to_device

    ps = Flux.params(AC)
    gs = gradient(ps) do
        logits = AC.actor(states_flattened)
        probs = softmax(logits)
        log_probs = logsoftmax(logits)
        log_probs_select = log_probs[actions]
        values = AC.critic(states_flattened)
        advantage = vec(gains) .- vec(values)
        actor_loss = -mean(log_probs_select .* Zygote.dropgrad(advantage))
        critic_loss = mean(advantage .^ 2)
        entropy_loss = -sum(probs .* log_probs) * 1 // size(probs, 2)
        loss = w₁ * actor_loss + w₂ * critic_loss - w₃ * entropy_loss
        ignore() do
            learner.actor_loss = actor_loss
            learner.critic_loss = critic_loss
            learner.entropy_loss = entropy_loss
            learner.loss = loss
        end
        loss
    end
    if !isnothing(learner.max_grad_norm)
        learner.norm = clip_by_global_norm!(gs, ps, learner.max_grad_norm)
    end
    update!(AC, gs)
end

RLCore.check(::QBasedPolicy{<:A2CLearner}, ::MultiThreadEnv) = nothing


text:
module InducingPoints

import Base: rand, show

using AbstractGPs: AbstractGPs
using Clustering: kmeans!
using DataStructures
using DeterminantalPointProcesses: DPP
using Distances
using KernelFunctions
using KernelFunctions: ColVecs, RowVecs, vec_of_vecs
using LinearAlgebra#: Symmetric, Eigen, eigen, eigvals, I, logdet, diag, norm
using Random: rand, bitrand, AbstractRNG, MersenneTwister, GLOBAL_RNG
using StatsBase: Weights, sample

export AbstractInducingPointsSelectionAlg

## Generic functions

export inducingpoints
export initZ
export updateZ!, updateZ

## Offline algorithms
export KmeansAlg
export RandomSubset
export StdDPP, kDPP
export Greedy

## Online algorithms
export OIPS
export SeqDPP
export StreamKmeans
export UniGrid
export Webscale

const jitt = 1e-5

abstract type AbstractInducingPointsSelectionAlg end

const AIPSA = AbstractInducingPointsSelectionAlg

abstract type OfflineInducingPointsSelectionAlg <: AIPSA end

const OffIPSA = OfflineInducingPointsSelectionAlg

abstract type OnlineInducingPointsSelectionAlg <: AIPSA end

const OnIPSA = OnlineInducingPointsSelectionAlg

## Wrapper for matrices
"""
     inducingpoints([rng::AbstractRNG], alg::OffIPSA, X::AbstractVector; kwargs...)
     inducingpoints([rng::AbstractRNG], alg::OffIPSA, X::AbstractMatrix; obsdim=1, kwargs...)

Select inducing points according to the algorithm `alg`.
"""
inducingpoints

function inducingpoints(
    rng::AbstractRNG, alg::AIPSA, X::AbstractMatrix; obsdim=1, kwargs...
)
    return inducingpoints(rng, alg, vec_of_vecs(X; obsdim=obsdim); kwargs...)
end

function inducingpoints(alg::AIPSA, X::AbstractMatrix; obsdim=1, kwargs...)
    return inducingpoints(GLOBAL_RNG, alg, X; obsdim=obsdim, kwargs...)
end

## Wrapper for the RNG
function inducingpoints(alg::AIPSA, X::AbstractVector; kwargs...)
    return inducingpoints(GLOBAL_RNG, alg, X; kwargs...)
end

## Online IP selection functions
"""
     initZ([rng::AbstractRNG], alg::OnIPSA, X::AbstractVector; kwargs...)
     initZ([rng::AbstractRNG], alg::OnIPSA, X::AbstractMatrix; obsdim=1, kwargs...)

Select inducing points according to the algorithm `alg` and return a Vector of Vector.
"""
initZ

initZ(Z::OnIPSA, X::AbstractVector; kwargs...) = initZ(GLOBAL_RNG, Z, X; kwargs...)

function initZ(alg::OnIPSA, X::AbstractMatrix; obsdim=1, kwargs...)
    return initZ(GLOBAL_RNG, alg, X; obsdim=obsdim, kwargs...)
end

function initZ(rng::AbstractRNG, alg::OnIPSA, X::AbstractMatrix; obsdim=1, kwargs...)
    return initZ(rng, vec_of_vecs(X; obsdim=obsdim), X; kwargs...)
end

"""
    updateZ!([rng::AbstractRNG], Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...)

Update inducing points `Z` with data `X` and algorithm `alg`
"""
updateZ!

function updateZ!(Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...)
    return updateZ!(GLOBAL_RNG, Z, alg, X; kwargs...)
end
# Default behavior is to simply add points
function updateZ!(
    rng::AbstractRNG, Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...
)
    return add_point!(rng, Z, alg, X; kwargs...)
end

"""
    updateZ([rng::AbstractRNG], Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...)

Return new vector of inducing points `Z` with data `X` and algorithm `alg` without changing the original one
"""
updateZ

function updateZ(Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...)
    return updateZ(GLOBAL_RNG, Z, alg, X; kwargs...)
end
# Default behavior is to simply add points
function updateZ(
    rng::AbstractRNG, Z::AbstractVector, alg::OnIPSA, X::AbstractVector; kwargs...
)
    return add_point(rng, Z, alg, X; kwargs...)
end

## Offline algorithms
include(joinpath("offline", "kmeans.jl"))
include(joinpath("offline", "randomsubset.jl"))
include(joinpath("offline", "stddpp.jl"))
include(joinpath("offline", "kdpp.jl"))
include(joinpath("offline", "greedyip.jl"))

## Online algorithms
include(joinpath("online", "seqdpp.jl"))
include(joinpath("online", "streamkmeans.jl"))
include(joinpath("online", "webscale.jl"))
include(joinpath("online", "oips.jl"))
include(joinpath("online", "unigrid.jl"))

## Utilities
include("utils.jl")

end


text:
using PersistenceDiagramsBase

using Compat
using DataFrames
using Test

@testset "PersistenceInterval" begin
    int1 = PersistenceInterval(1, 2)
    int2 = PersistenceInterval(1, Inf)
    r = [1, 2, 6, 3]
    int3 = PersistenceInterval(1, 2; birth_simplex=:σ, death_simplex=:τ, representative=r)
    int4 = PersistenceInterval(int1; birth_simplex=:σ, death_simplex=:τ, representative=r)

    @testset "Equality, order" begin
        @test int1 ≠ int2
        @test int1 == int3
        @test int1 < int2
        @test isless(int3, int2)
    end

    @testset "Comparison with tuples" begin
        @test int1 == (1, 2)
        @test int2 == (1, Inf)
        @test (1.0, 2.0) == int4
        @test int1 == int3

        @test PersistenceInterval((1, 2)) == int1
    end

    @testset "Birth, death, iteration" begin
        @test int1[1] == birth(int1) == 1
        @test int1[2] == death(int1) == 2
        @test int3[1] == birth(int1) == 1
        @test int4[2] == death(int1) == 2
        @test isfinite(int1)
        @test !isfinite(int2)
        @test persistence(int1) == 1
        @test persistence(int2) == Inf
        @test midlife(int1) == 1.5
        @test midlife(int2) == Inf

        @test eltype(int1) ≡ Float64
        @test eltype(PersistenceInterval) ≡ Float64
        @test length(int1) == 2
        @test collect(int1) == [1, 2]
        @test tuple(int1...) ≡ (1.0, 2.0)
        @test firstindex(int1) == 1
        @test lastindex(int1) == 2
        @test first(int1) == 1
        @test last(int1) == 2
        @test Base.IteratorSize(int1) == Base.HasLength()
        @test Base.IteratorEltype(int1) == Base.HasEltype()

        @test_throws BoundsError int2[0]
        @test_throws BoundsError int2[3]
    end

    @testset "Metadata access" begin
        @test int1.meta == NamedTuple()
        @test birth_simplex(int3) == :σ
        @test int4.birth_simplex == :σ
        @test death_simplex(int4) == :τ
        @test int3.death_simplex == :τ
        @test representative(int3) == r

        @test propertynames(int1) == (:birth, :death)
        @test propertynames(int1, true) == (:birth, :death, :meta)
        @test propertynames(int3) ==
              (:birth, :death, :birth_simplex, :death_simplex, :representative)
        @test propertynames(int3, true) ==
              (:birth, :death, :birth_simplex, :death_simplex, :representative, :meta)

        @test_throws ErrorException birth_simplex(int1)
        @test_throws ErrorException death_simplex(int1)
        @test_throws ErrorException representative(int1)
        @test_throws ErrorException int3.something
    end

    @testset "Printing" begin
        @test sprint(print, int1) == "[1.0, 2.0)"
        @test sprint(print, int2) == "[1.0, ∞)"
        @test sprint(print, int3) == "[1.0, 2.0)"

        print_text_plain(io, x) = show(io, MIME"text/plain"(), x)
        @test sprint(print_text_plain, int1) == "[1.0, 2.0)"
        @test sprint(print_text_plain, int2) == "[1.0, ∞)"
        @test sprint(print_text_plain, int3) ==
              "[1.0, 2.0) with:\n" *
              " birth_simplex: Symbol\n" *
              " death_simplex: Symbol\n" *
              " representative: 4-element $(typeof(r))"
    end
end

@testset "PersistenceDiagram" begin
    diagram1 = PersistenceDiagram(
        [
            PersistenceInterval(1, 3; a=1),
            PersistenceInterval(3, 4; a=2),
            PersistenceInterval(3, Inf; a=3),
        ];
        dim=1,
    )
    diagram2 = PersistenceDiagram([(1, 3), (3, 4), (3, Inf)]; threshold=0.3)
    diagram3 = PersistenceDiagram(
        view(
            [
                PersistenceInterval(1, 2),
                PersistenceInterval(3, 4),
                PersistenceInterval(3, Inf),
            ],
            1:3,
        );
        a=1,
    )

    @testset "A persistence diagram is an array" begin
        @test diagram1[1] == (1, 3)
        @test diagram2[2] == (3, 4)
        @test diagram3[3] == (3, Inf)

        @test_throws BoundsError diagram1[0]
        @test_throws BoundsError diagram2[4]

        @test diagram1 == diagram2
        @test diagram1 == [(1, 3), (3, 4), (3, Inf)]

        @test length(diagram1) == 3
        @test firstindex(diagram2) == 1
        @test lastindex(diagram3) == 3
        @test length(diagram1) == 3
        @test size(diagram2) == (3,)

        @test first(diagram1) == (1, 3)
        @test last(diagram2) == (3, Inf)

        @test copy(diagram1) == diagram1
        @test copy(diagram2).threshold == 0.3

        @test similar(diagram1) isa typeof(diagram1)
        @test similar(diagram1).dim == 1
        @test similar(diagram3, (Base.OneTo(2),)) isa typeof(diagram3)
        @test similar(diagram3, (Base.OneTo(2),)).a == 1

        @test sort(diagram3) isa typeof(diagram3)
        @test sort(diagram2; by=death, rev=true) == [(3, Inf), (3, 4), (1, 3)]
        @test sort(diagram3; by=death, rev=true).a == 1
    end

    @testset "Metadata access" begin
        @test dim(diagram1) == diagram1.dim == 1
        @test threshold(diagram2) == diagram2.threshold == 0.3
        @test diagram3.a == 1

        @test propertynames(diagram1) == (:intervals, :dim)
        @test propertynames(diagram1, true) == (:intervals, :meta, :dim)

        @test_throws ErrorException diagram1.threshold
        @test_throws ErrorException diagram2.dim

        @test diagram1[1].a == 1
        @test diagram1[2].a == 2
        @test diagram1[3].a == 3
    end

    if VERSION ≥ v"1.5.0"
        @testset "Printing" begin
            print_text_plain(io, x) = show(io, MIME"text/plain"(), x)

            @test sprint(print, diagram1) == "3-element 1-dimensional PersistenceDiagram"
            @test sprint(print_text_plain, diagram1) ==
                  "3-element 1-dimensional PersistenceDiagram:\n" *
                  " [1.0, 3.0)\n" *
                  " [3.0, 4.0)\n" *
                  " [3.0, ∞)"

            @test sprint(print, diagram2) == "3-element PersistenceDiagram"
            @test sprint(print_text_plain, diagram2) ==
                  "3-element PersistenceDiagram:\n" *
                  " [1.0, 3.0)\n" *
                  " [3.0, 4.0)\n" *
                  " [3.0, ∞)"
        end
    end
end

@testset "Tables.jl interface" begin
    diag1 = PersistenceDiagram(
        [PersistenceInterval(1, 2), PersistenceInterval(1, 3)]; dim=0, threshold=4
    )

    df = DataFrame(diag1)
    @test names(df) == ["birth", "death", "dim", "threshold"]
    @test nrow(df) == 2
    @test PersistenceDiagram(df) == diag1

    diag2 = PersistenceDiagram(
        [
            PersistenceInterval(1, 2; a=nothing),
            PersistenceInterval(1, 3; a=nothing),
            PersistenceInterval(1, 4; a=1),
        ];
        dim=1,
        b=2,
    )

    df = DataFrame(diag2)
    @test names(df) == ["birth", "death", "dim", "threshold"]
    @test nrow(df) == 3
    @test all(ismissing, df.threshold)

    table = Tables.columntable((threshold=[1, 1], birth=[0, 0], death=[0, 0]))
    diagram = PersistenceDiagram(table)
    @test ismissing(dim(diagram))
    @test threshold(diagram) == 1

    table = Tables.columntable((dim=[1, 1], birth=[0, 0], death=[0, 0]))
    diagram = PersistenceDiagram(table)
    @test ismissing(threshold(diagram))
    @test dim(diagram) == 1

    table = Tables.columntable((birth=[], death=[]))
    @test isempty(PersistenceDiagram(table))

    @test_throws ErrorException PersistenceDiagram(Tables.columntable((
        birth=[1, 1], death=[2, 2], threshold=[1, 2]
    )))
    @test_throws ErrorException PersistenceDiagram(Tables.columntable((
        birth=[1, 1], death=[2, 2], dim=[1, 2]
    )))
    @test_throws ErrorException PersistenceDiagram(Tables.columntable((
        death=[2, 2], dim=[1, 1]
    )))
    @test_throws ErrorException PersistenceDiagram(Tables.columntable((
        birth=[2, 2], dim=[1, 1]
    )))
end


text:
"""
    SerialExecution

Functor object for using serial execution in a progressive-hedging algorithm. Create by supplying a [`Serial`](@ref) object through `execution` in the `ProgressiveHedgingSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct SerialExecution{T <: AbstractFloat,
                       A <: AbstractVector,
                       PT <: AbstractPenaltyterm} <: AbstractProgressiveHedgingExecution
    subproblems::Vector{SubProblem{T,A,PT}}

    function SerialExecution(::Type{T}, ::Type{A}, ::Type{PT}) where {T <: AbstractFloat, A <: AbstractVector, PT <: AbstractPenaltyterm}
        return new{T,A,PT}(Vector{SubProblem{T,A,PT}}())
    end
end

function initialize_subproblems!(ph::AbstractProgressiveHedging,
                                 execution::SerialExecution{T},
                                 scenarioproblems::ScenarioProblems,
                                 penaltyterm::AbstractPenaltyterm) where T <: AbstractFloat
    # Create subproblems
    for i = 1:num_subproblems(scenarioproblems)
        push!(execution.subproblems, SubProblem(
            subproblem(scenarioproblems, i),
            i,
            T(probability(scenarioproblems, i)),
            copy(penaltyterm)))
    end
    # Initial reductions
    update_iterate!(ph)
    update_dual_gap!(ph)
    return nothing
end

function finish_initilization!(execution::SerialExecution, penalty::AbstractFloat)
    for subproblem in execution.subproblems
        initialize!(subproblem, penalty)
    end
    return nothing
end

function restore_subproblems!(ph::AbstractProgressiveHedging, execution::SerialExecution)
    for subproblem in execution.subproblems
        restore_subproblem!(subproblem)
    end
    return nothing
end

function resolve_subproblems!(ph::AbstractProgressiveHedging, execution::SerialExecution)
    Qs = Vector{SubproblemSolution}(undef, length(execution.subproblems))
    # Reformulate and solve sub problems
    for (i, subproblem) in enumerate(execution.subproblems)
        reformulate_subproblem!(subproblem, ph.ξ, penalty(ph))
        Qs[i] = subproblem(ph.ξ)
    end
    # Return current objective value
    return sum(Qs)
end

function update_iterate!(ph::AbstractProgressiveHedging, execution::SerialExecution)
    # Update the estimate
    ξ_prev = copy(ph.ξ)
    ph.ξ .= mapreduce(+, execution.subproblems, init = zero(ph.ξ)) do subproblem
        π = subproblem.probability
        x = subproblem.x
        π * x
    end
    # Update δ₁
    ph.data.δ₁ = norm(ph.ξ - ξ_prev, 2) ^ 2
    return nothing
end

function update_subproblems!(ph::AbstractProgressiveHedging, execution::SerialExecution)
    # Update dual prices
    update_subproblems!(execution.subproblems, ph.ξ, penalty(ph))
    return nothing
end

function update_dual_gap!(ph::AbstractProgressiveHedging, execution::SerialExecution{T}) where T <: AbstractFloat
    # Update δ₂
    ph.data.δ₂ = mapreduce(+, execution.subproblems, init = zero(T)) do subproblem
        π = subproblem.probability
        x = subproblem.x
        π * norm(x - ph.ξ, 2) ^ 2
    end
    return nothing
end

function calculate_objective_value(ph::AbstractProgressiveHedging, execution::SerialExecution{T}) where T <: AbstractFloat
    return mapreduce(+, execution.subproblems, init = zero(T)) do subproblem
        _objective_value(subproblem)
    end
end

function scalar_subproblem_reduction(value::Function, execution::SerialExecution{T}) where T <: AbstractFloat
    return mapreduce(+, execution.subproblems, init = zero(T)) do subproblem
        π = subproblem.probability
        return π * value(subproblem)
    end
end

function vector_subproblem_reduction(value::Function, execution::SerialExecution{T}, n::Integer) where T <: AbstractFloat
    return mapreduce(+, execution.subproblems, init = zero(T, n)) do subproblem
        π = subproblem.probability
        return π * value(subproblem)
    end
end

# API
# ------------------------------------------------------------
function (execution::Serial)(::Type{T}, ::Type{A}, ::Type{PT}) where {T <: AbstractFloat, A <: AbstractVector, PT <: AbstractPenaltyterm}
    return SerialExecution(T, A, PT)
end

function str(::Serial)
    return ""
end


text:
using FastTransforms, Makie, Random

include("SLPDE.jl")
include("etdrk4.jl")
include("evaluate_lambda.jl")
include("sphcesaro.jl")

Random.seed!(0)

n = 64
U0 = zeros(4n, 8n-1);
U0[1:n, 1:2n-1] = sphrandn(Float64, n, 2n-1)/n
n = size(U0, 1)

θ = [0;(0.5:n-0.5)/n;1]
φ = [(0:2n-2)*2/(2n-1);2]
x = Float32[cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
y = Float32[sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
z = Float32[cospi(θ) for θ in θ, φ in φ]

P = plan_sph2fourier(U0)
PS = plan_sph_synthesis(U0)
F = deepcopy(U0)
Ft = zeros(Float32, size(F, 1)+2, size(F, 2)+1)
lmul!(P, F)
lmul!(PS, F)
fillF!(Ft, F)

scene = Scene(resolution = (1200, 1200));
surf = surface!(scene, x, y, z, color = Ft, colormap = :viridis, colorrange = (-1.0, 1.0))[end];
update_cam!(scene, Vec3f0(2.5), Vec3f0(0), Vec3f0(0, 0, 1))
scene.center = false
scene

T = 8.0
n = round(Int, 6*T)

SLPDE = SemiLinearPDE(LaplaceBeltrami(0.02), NonlinearOperator(u->u-u^3), U0)
io = VideoStream(scene);
ETDRK4(SLPDE, T, n, surf, io)
Makie.save("plots/LAC.gif", io)

SLPDE = SemiLinearPDE(NonlocalLaplaceBeltrami(0.02, 0.0, 0.2), NonlinearOperator(u->u-u^3), U0)
io = VideoStream(scene);
ETDRK4(SLPDE, T, n, surf, io)
Makie.save("plots/NAC.gif", io)


text:
module SpecUtils

using StatsBase, PyPlot, Statistics

export combine, binning

"""
    `combine(x::Array{Array{T,1},1}, y::Array{Array{T,N},1}) where {T, N}`

Combines and sorts datasets. x values of each set is stored in `x` y values
are stored in `y`.
"""
function combine(x::Array{Array{T,1},1}, y::Array{Array{T,N},1}) where {T, N}
    xc = vcat(x...)
    yc = vcat(y...)
    xyc = hcat(xc, yc)
    xys = sortslices(xyc, dims=1)
    xc .= xys[:,1]
    N == 2 && (yc .= xys[:,2:end])
    N == 1 && (yc .= xys[:,2])
    return xc, yc
end

"""
    `binning(x::Array{T,1}, y::Array{T,1}, n::Int; mode=:in, debug=false) where T <: Number`

Averages the array `y` appropriately over a certain amount of 'x' values
to return an array with 'n' data points.

Keyword Arguments:
    * `binning` = `:in` or `:out` controls the position of the first
        and last bin. `:out` centers the first and last data point in the first and
        last bin respectively. `:in` places the first bins left edge to the first
        data point and the last bins right edge to the last data point.

    * `debug` = `true` or `false`: Generate a plot for debugging.
"""
function binning(x::Array{T,1}, y::Array{T,1}, n::Int; mode=:in, debug=false) where T <: Number
    @assert length(x) == length(y)
    @assert issorted(x)

    xr = Array{Float64,1}(undef, n)
    yr = Array{Float64,1}(undef, n)

    # We want to place the edges in such a way that the data points lie exactly
    # in the middle of the bins when the number of resampled points matches the
    # number of native points. Therefore edges_tight has to be adjusted.
    weight = ProbabilityWeights(ones(T, length(x)))

    if mode == :out
        edges = range(x[1], stop=x[end], length=n)
        edges_start = x[1]   - 1/2 * step(edges)
        edges_end   = x[end] + 1/2 * step(edges)
        edges = range(edges_start, edges_end, length=n+1)
    elseif mode == :in
        edges = range(x[1] - eps(Float64), stop=x[end], length=n+1)
    else
        error("Mode $mode not defined. Choose :in or :out.")
    end

    hist = fit(Histogram, x, weight, edges, closed=:right)

    i = 1
    j = 1
    while j <= n
        w = trunc(Int, i + hist.weights[j] - 1)
        yr[j] = mean(y[i:w])
        xr[j] = (edges[j] + edges[j+1]) / 2
        i += trunc(Int, hist.weights[j])
        j += 1
    end

    if debug == true
        figure()
        title("Binning Debug")
        scatter(x, y)
        plot(xr, yr, "r.-")
        vlines(edges, ylim()..., alpha=0.5, color="k")
        for i = 1:length(hist.weights)
            isodd(i) ? (p = 0.85) : (p = 0.81)
            text(xr[i], p * ylim()[2] - ylim()[1], trunc.(Int, hist.weights[i]),
                horizontalalignment="center",
                color="k")
        end
        show()
    end

    xr, yr
end

end # module

max_stars_repo_name: MLackner/SpecUtils

max_issues_repo_name: MLackner/SpecUtils

max_forks_repo_name: MLackner/SpecUtils


text:
"""
`Basics.tools()`  ... select different tools from a menu for which no further input is required.
"""
function Basics.tools(x::Int64)
    println("JAC.Tools.select(): A menu if nothing comes in")
    #
    t1 = "Simple tools JAC without input: "
    b1 = dropdown(["Start task", "Grid calculator", "Another task"])
    b2 = slider(1:100, label = "To what extent?", value = 33)
    update = button("Update")
    ui = vbox( hbox( pad(1em, t1) ),
               hbox( pad(1em, b1), pad(1em, b2), pad(1em, update) )
             )
    Interact.display(ui)
    output = Interact.@map (&update;  (observe(b1)[], observe(b2)[]) )
    return( output )
end

"""
`Basics.tools(dict::Dict)`
    ... select different tools from a menu if proper results::dict are given in terms of a dictionary as obtained
        usually from an Atomic.Computation.
"""
function Basics.tools(dict::Dict)
    println("JAC.Tools.select(): A menu if a dict comes in; $dict ")
end


text:
function reciprocal(obj::FieldElement)
    return jcall(obj, "reciprocal", Object, ())
end


text:
module DelayedErrors

export push_delayed_error, pop_delayed_errors

struct DelayedError
    message::T1 where T1 <: AbstractString
    dict::T2 where T2 <: AbstractDict
    bt::T3 where T3
    st::T4 where T4
end

@noinline function DelayedError(
        message::T1,
        dict::T2,
        )::DelayedError where
        T1 <: AbstractString where
        T2 <: AbstractDict
    bt = backtrace()
    st = stacktrace(bt)
    x = DelayedError(
        message,
        dict,
        bt,
        st,
        )
    return x
end

@noinline function __init__()::Nothing
    global delayed_error_list = Vector{DelayedError}()
    return nothing
end

"""
"""
function push_delayed_error end

"""
"""
function pop_delayed_errors end

@noinline function push_delayed_error(
        message::S;
        kwargs...
        )::Nothing where S <: AbstractString
    x = DelayedError(message, Dict(kwargs...),)
    global delayed_error_list
    push!(delayed_error_list, x,)
    @error("Pushed delayed error: $(x.message)", x.dict...)
    Base.show_backtrace(stderr, x.bt,)
    println(stderr)
    println(stderr)
    return nothing
end

@noinline function push_delayed_error(
        message::Vararg{Any,N};
        kwargs...
        )::Nothing where {N}
    push_delayed_error(Main.Base.string(message); kwargs...)
    return nothing
end

@noinline function pop_delayed_errors()::Nothing
    global delayed_error_list
    if isempty(delayed_error_list)
        @debug("There were no delayed errors to pop.")
    else
        while !isempty(delayed_error_list)
            x = popfirst!(delayed_error_list)
            @error("Popped delayed error: $(x.message)", x.dict...)
            Base.show_backtrace(stderr, x.bt,)
            println(stderr)
            println(stderr)
        end
        error("There were one or more delayed errors.")
    end
    return nothing
end

end # end module DelayedErrors


text:
module MosaicViews

using ImageCore
using PaddedViews

export

    MosaicView,
    mosaicview

"""
    MosaicView(A::AbstractArray)

Create a two dimensional "view" of the three or four dimensional
array `A`. The resulting `MosaicView` will display the data in
`A` such that it emulates using `vcat` for all elements in the
third dimension of `A`, and `hcat` for all elements in the fourth
dimension of `A`.

For example, if `size(A)` is `(2,3,4)`, then the resulting
`MosaicView` will have the size `(2*4,3)` which is `(8,3)`.
Alternatively, if `size(A)` is `(2,3,4,5)`, then the resulting
size will be `(2*4,3*5)` which is `(8,15)`.

Another way to think about this is that `MosaicView` creates a
mosaic of all the individual matrices enumerated in the third
(and optionally fourth) dimension of the given 3D or 4D array
`A`. This can be especially useful for creating a single
composite image from a set of equally sized images.

```@jldoctest
julia> using MosaicViews

julia> A = [(k+1)*l-1 for i in 1:2, j in 1:3, k in 1:2, l in 1:2]
2×3×2×2 Array{Int64,4}:
[:, :, 1, 1] =
 1  1  1
 1  1  1

[:, :, 2, 1] =
 2  2  2
 2  2  2

[:, :, 1, 2] =
 3  3  3
 3  3  3

[:, :, 2, 2] =
 5  5  5
 5  5  5

julia> MosaicView(A)
4×6 MosaicViews.MosaicView{Int64,4,Array{Int64,4}}:
 1  1  1  3  3  3
 1  1  1  3  3  3
 2  2  2  5  5  5
 2  2  2  5  5  5
```
"""
struct MosaicView{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,2}
    parent::A
    dims::Tuple{Int,Int}
    pdims::NTuple{N,Int}

    function MosaicView{T,N}(A::AbstractArray{T,N}, dims) where {T,N}
        3 <= N <= 4 || throw(ArgumentError("The given array must have dimensionality N=3 or N=4"))
        new{T,N,typeof(A)}(A, dims, size(A))
    end
end

function MosaicView(A::AbstractArray{T,N}) where {T,N}
    dims = (size(A,1) * size(A,3), size(A,2) * size(A,4))
    MosaicView{T,N}(A, dims)
end

Base.parent(mv::MosaicView) = mv.parent
Base.size(mv::MosaicView) = mv.dims

@inline function Base.getindex(mv::MosaicView{T,3,A}, i::Int, j::Int) where {T,A}
    @boundscheck checkbounds(mv, i, j)
    pdims = mv.pdims
    parent = mv.parent
    idx1 = (i-1) % pdims[1] + 1
    idx2 = (j-1) % pdims[2] + 1
    idx3 = (i-1) ÷ pdims[1] + 1
    @inbounds res = parent[idx1, idx2, idx3]
    res
end

@inline function Base.getindex(mv::MosaicView{T,4,A}, i::Int, j::Int) where {T,A}
    @boundscheck checkbounds(mv, i, j)
    pdims = mv.pdims
    parent = mv.parent
    idx1 = (i-1) % pdims[1] + 1
    idx2 = (j-1) % pdims[2] + 1
    idx3 = (i-1) ÷ pdims[1] + 1
    idx4 = (j-1) ÷ pdims[2] + 1
    @inbounds res = parent[idx1, idx2, idx3, idx4]
    res
end

"""
    mosaicview(A::AbstractArray, [fill]; [npad=0], [nrow], [ncol], [rowmajor=false]) -> MosaicView

Create a two dimensional "view" of the higher dimensional array
`A`. The resulting [`MosaicView`](@ref) will display all the
matrix slices of the first two dimensions of `A` arranged as a
single large mosaic (in the form of a matrix).

In contrast to using the constructor of [`MosaicView`](@ref)
directly, the function `mosaicview` also allows for a couple of
convenience keywords. Note that as a consequence the function is
not type stable and should only be used if performance is not a
priority. A typical use case would be to create an image mosaic
from a set of equally sized input images.

- The optional positional parameter `fill` defines the value that
  that should be used for empty space. This can be padding caused
  by `npad`, or empty mosaic tiles in case the number of matrix
  slices in `A` is smaller than `nrow*ncol`.

- The parameter `npad` defines the empty padding space between
  adjacent mosaic tiles. This can be especially useful if the
  individual tiles (i.e. matrix slices in `A`) are images that
  should be visually separated by some grid lines.

- The parameters `nrow` and `ncol` can be used to choose the
  number of rows and/or columns the mosaic should be arranged in.
  Note that it suffices to specify one of the two parameters, as
  the other one can be inferred accordingly. The default in case
  none of the two are specified is `nrow = size(A,3)`.

- If `rowmajor` is set to `true`, then the slices will be
  arranged left-to-right-top-to-bottom, instead of
  top-to-bottom-left-to-right (default).

# Examples

```julia-repl
julia> using MosaicViews

julia> A = [k for i in 1:2, j in 1:3, k in 1:5]
2×3×5 Array{Int64,3}:
[:, :, 1] =
 1  1  1
 1  1  1

[:, :, 2] =
 2  2  2
 2  2  2

[:, :, 3] =
 3  3  3
 3  3  3

[:, :, 4] =
 4  4  4
 4  4  4

[:, :, 5] =
 5  5  5
 5  5  5

julia> mosaicview(A, ncol=2)
6×6 MosaicViews.MosaicView{Int64,4,...}:
 1  1  1  4  4  4
 1  1  1  4  4  4
 2  2  2  5  5  5
 2  2  2  5  5  5
 3  3  3  0  0  0
 3  3  3  0  0  0

julia> mosaicview(A, nrow=2)
4×9 MosaicViews.MosaicView{Int64,4,...}:
 1  1  1  3  3  3  5  5  5
 1  1  1  3  3  3  5  5  5
 2  2  2  4  4  4  0  0  0
 2  2  2  4  4  4  0  0  0

julia> mosaicview(A, nrow=2, rowmajor=true)
4×9 MosaicViews.MosaicView{Int64,4,...}:
 1  1  1  2  2  2  3  3  3
 1  1  1  2  2  2  3  3  3
 4  4  4  5  5  5  0  0  0
 4  4  4  5  5  5  0  0  0

julia> mosaicview(A, nrow=2, npad=1, rowmajor=true)
5×11 MosaicViews.MosaicView{Int64,4,...}:
 1  1  1  0  2  2  2  0  3  3  3
 1  1  1  0  2  2  2  0  3  3  3
 0  0  0  0  0  0  0  0  0  0  0
 4  4  4  0  5  5  5  0  0  0  0
 4  4  4  0  5  5  5  0  0  0  0

julia> mosaicview(A, -1, nrow=2, npad=1, rowmajor=true)
5×11 MosaicViews.MosaicView{Int64,4,...}:
  1   1   1  -1   2   2   2  -1   3   3   3
  1   1   1  -1   2   2   2  -1   3   3   3
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
  4   4   4  -1   5   5   5  -1  -1  -1  -1
  4   4   4  -1   5   5   5  -1  -1  -1  -1
```
"""
function mosaicview(A::AbstractArray{T,3},
                    fill = zero(T);
                    npad = 0,
                    nrow = -1,
                    ncol = -1,
                    rowmajor = false) where T
    nrow == -1 || nrow > 0 || throw(ArgumentError("The parameter \"nrow\" must be greater than 0"))
    ncol == -1 || ncol > 0 || throw(ArgumentError("The parameter \"ncol\" must be greater than 0"))
    npad >= 0 || throw(ArgumentError("The parameter \"npad\" must be greater than or equal to 0"))
    ntile = size(A,3)
    ntile_ceil = ntile # ntile need not be integer divideable
    if nrow == -1 && ncol == -1
        # automatically choose nrow to reflect what MosaicView does
        nrow = ntile
        ncol = 1
    elseif nrow == -1
        # compute nrow based on ncol
        nrow = ceil(Int, ntile / ncol)
        ntile_ceil = nrow * ncol
    elseif ncol == -1
        # compute ncol based on nrow
        ncol = ceil(Int, ntile / nrow)
        ntile_ceil = nrow * ncol
    else
        # accept nrow and ncol as is if it covers at least all
        # existing tiles
        ntile_ceil = nrow * ncol
        ntile_ceil < ntile && throw(ArgumentError("The product of the parameters \"ncol\" (value: $ncol) and \"nrow\" (value: $nrow) must be greater than or equal to $ntile"))
    end
    # we pad size(A,3) to nrow*ncol. we also pad the first two
    # dimensions according to npad. think of this as border
    # between tiles (useful for images)
    pad_dims = (size(A,1) + npad, size(A,2) + npad, ntile_ceil)
    A_pad = PaddedView(T(fill), A, pad_dims)
    # next we reshape the image such that it reflects the
    # specified nrow and ncol
    A_new = if !rowmajor
        res_dims = (size(A_pad,1), size(A_pad,2), nrow, ncol)
        reshape(A_pad, res_dims)
    else
        # same as above but we additionally permute dimensions
        # to mimic row first layout for the tiles. this is useful
        # for images since user often reads these from left-to-right
        # before top-to-bottom. (note the swap of "ncol" and "nrow")
        res_dims = (size(A_pad,1), size(A_pad,2), ncol, nrow)
        A_tp = reshape(A_pad, res_dims)
        permuteddimsview(A_tp, (1, 2, 4, 3))
    end
    # decrease size of the resulting MosaicView by npad to not have
    # a border on the right side and bottom side of the final mosaic.
    dims = (size(A_new,1) * size(A_new,3) - npad, size(A_new,2) * size(A_new,4) - npad)
    MosaicView{T,4}(A_new, dims)
end

function mosaicview(A::AbstractArray{T,N},
                    args...;
                    npad = 0,
                    nrow = -1,
                    ncol = -1,
                    rowmajor = false) where {T,N}
    3 <= N || throw(ArgumentError("The given array must have dimensionality of N=3 or higher"))
    # if neither nrow nor ncol is provided then automatically choose
    # nrow and ncol to reflect what MosaicView does (i.e. use size)
    if nrow == -1 && ncol == -1
        nrow = size(A, 3)
        # ncol = size(A, 4)
    end
    mosaicview(reshape(A, (size(A,1), size(A,2), :)), args...;
               npad=npad, nrow=nrow, ncol=ncol, rowmajor=rowmajor)
end

end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

module CoDa

using CSV
using Tables
using DataFrames
using StatsBase
using StaticArrays
using UnicodePlots

import Base: +, -, *, ==
import Base: names, getproperty
import LinearAlgebra: norm, dot, ⋅

include("composition.jl")
include("transforms.jl")
include("utils.jl")

export
  # composition
  Composition,
  parts, components,
  norm, dot, ⋅,
  distance,

  # transforms
  alr, alrinv,
  clr, clrinv,
  ilr, ilrinv,

  # utils
  readcoda,
  compose

end


text:
using LinearAlgebra
using ProximalOperators
using Random
using Test

Random.seed!(0)

@testset "IndHyperslab" begin

for R in [Float16, Float32, Float64]
    for shape in [(5,), (3, 5), (3, 4, 5)]
        c = randn(R, shape)
        x = randn(R, shape)
        cx = dot(c, x)

        for (low, upp) in [(cx-R(1), cx+R(1)), (cx-R(2), cx-R(1)), (cx+R(1), cx+R(2))]
            f = IndHyperslab(low, c, upp)
            predicates_test(f)
            call_test(f, x)
            prox_test(f, x, R(0.5)+rand(R))
        end
    end
end

end


text:
"FailUnsat on condition"
struct FailUnsatAlg <: SamplingAlgorithm end

softhard(::Type{FailUnsatAlg}) = IsHard{FailUnsatAlg}()

const FailUnsat = FailUnsatAlg()

isapproximate(::FailUnsatAlg) = false

"`n` samples from `x`"
function Base.rand(rng::AbstractRNG,
                   ΩT::Type{OT},
                   pred::RandVar,
                   n::Integer,
                   alg::FailUnsatAlg) where {OT <: Ω}
  ωsamples = ΩT[]
  ωsamples = Vector{ΩT}(undef, n)
  accepted = 0
  for i = 1:n
    ω = ΩT()
    issat = apl(pred, Omega.Space.tagrng(ω, rng))
    !issat && error("Condition unsatisfied. Use appropriate infrence alg.")
    @inbounds ωsamples[i] = ω
    lens(Loop, (ω = ω, accepted = accepted, p = float(issat), i = i))
  end
  ωsamples
end

function Base.rand(rng::AbstractRNG,
                   x::RandVar,
                   n::Integer,
                   alg::FailUnsatAlg;
                   ΩT::Type{OT} = defΩ(alg)) where {OT <: Ω}
  pred = Omega.indomain(x)
  ωsamples = rand(rng, ΩT, pred, n, alg)
  map(ω -> apl(x, ω), ωsamples)
end


text:
#-----------------------------------------------------------------------
#   VECTOR AND MATRIX DATATYPES
#-----------------------------------------------------------------------

const RealVec{T<:Real} = Array{T,1}      # Defining vector datatype
const RealMtx{T<:Real} = Array{T,2}      # Defining matrix datatype

#-----------------------------------------------------------------------


text:
module NBodyLeapFrog

using LinearAlgebra, NiLang

export Bodies, leapfrog
export Body, V3
export LeapFrog, PEFRL, FR

include("V3.jl")
include("starData.jl")
include("common.jl")
include("leapfrog.jl")
include("reversible_leapfrog.jl")
include("pefrl.jl")
include("reversible_pefrl.jl")

end


text:
# Before running this, please make sure to activate and instantiate the
# tutorial-specific package environment, using this
# [`Project.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/ISL-lab-9/Project.toml) and
# [this `Manifest.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/ISL-lab-9/Manifest.toml), or by following
# [these](https://juliaai.github.io/DataScienceTutorials.jl/#learning_by_doing) detailed instructions.

# ## Getting started

using MLJ
import RDatasets: dataset
using PrettyPrinting
using Random

# We start by generating a 2D cloud of points

Random.seed!(3203)
X = randn(20, 2)
y = vcat(-ones(10), ones(10))

# which we can visualise

using PyPlot
figure(figsize=(8,6))

ym1 = y .== -1
ym2 = .!ym1
plot(X[ym1, 1], X[ym1, 2], ls="none", marker="o")
plot(X[ym2, 1], X[ym2, 2], ls="none", marker="x")

# \figalt{Toy points}{ISL-lab-9-g1.svg}

# let's wrap the data as a table:

X = MLJ.table(X)
y = categorical(y);

# and fit a SVM classifier

SVC = @load SVC pkg=LIBSVM

svc_mdl = SVC()
svc = machine(svc_mdl, X, y)

fit!(svc);

# As usual we can check how it performs

ypred = MLJ.predict(svc, X)
misclassification_rate(ypred, y)

# Not bad.

# ### Basic tuning

# As usual we could tune the model, for instance the penalty encoding the tradeoff between margin width and misclassification:

rc = range(svc_mdl, :cost, lower=0.1, upper=5)
tm = TunedModel(model=svc_mdl, ranges=[rc], tuning=Grid(resolution=10),
                resampling=CV(nfolds=3, rng=33), measure=misclassification_rate)
mtm = machine(tm, X, y)

fit!(mtm)

ypred = MLJ.predict(mtm, X)
misclassification_rate(ypred, y)

# You could also change the kernel etc.

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"""
```
type AltPolicy
```

Type defining an alternative policy rule.

### Fields

- `key::Symbol`: alternative policy identifier

- `eqcond::Function`: a version of `DSGE.eqcond` which computes the equilibrium
  condition matrices under the alternative policy. Like `eqcond`, it should take
  in one argument of type `AbstractModel` and return the `Γ0`, `Γ1`, `C`, `Ψ`,
  and `Π` matrices.

- `solve::Function`: a version of `DSGE.solve` which solves the model under the
  alternative policy. Like `DSGE.solve`, it should take in one argument of type
  `AbstractModel` and return the `TTT`, `RRR`, and `CCC` matrices.

- `forecast_init::Function`: a function that initializes forecasts under the
  alternative policy rule. Specifically, it accepts a model, an `nshocks` x
  `n_forecast_periods` matrix of shocks to be applied in the forecast, and a
  vector of initial states for the forecast. It must return a new matrix of
  shocks and a new initial state vector. If no adjustments to shocks or initial
  state vectors are necessary under the policy rule, this field may be omitted.

- `color::Colorant`: color to plot this alternative policy in. Defaults to blue.

- `linestyle::Symbol`: line style for forecast plots under this alternative
  policy. See options from `Plots.jl`. Defaults to `:solid`.

"""
type AltPolicy
    key::Symbol
    eqcond::Function
    solve::Function
    forecast_init::Function
    color::Colorant
    linestyle::Symbol
end

function AltPolicy(key::Symbol, eqcond_fcn::Function, solve_fcn::Function;
                   forecast_init::Function = identity,
                   color::Colorant = RGB(0., 0., 1.),
                   linestyle::Symbol = :solid)

    AltPolicy(key, eqcond_fcn, solve_fcn, forecast_init, color, linestyle)
end

Base.string(a::AltPolicy) = string(a.key)


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractNeighborhood

A neighborhood of points in a spatial object.
"""
abstract type AbstractNeighborhood end

"""
    isneighbor(neighborhood, xₒ, x)

Tells whether or not the coordinates `x` are in the `neighborhood`
centered at coordinates `xₒ`.
"""
isneighbor(neigh::AbstractNeighborhood,
           xₒ::AbstractVector, x::AbstractVector) = @error "not implemented"

#------------------
# IMPLEMENTATIONS
#------------------
include("neighborhoods/ball.jl")
include("neighborhoods/cylinder.jl")


text:
const KN_var_type_map = Dict(
    KN_VARTYPE_CONTINUOUS => :Cont,
    KN_VARTYPE_INTEGER => :Int,
    KN_VARTYPE_BINARY => :Bin
)

const KN_rev_var_type_map = Dict(
    :Cont => KN_VARTYPE_CONTINUOUS,
    :Int => KN_VARTYPE_INTEGER,
    :Bin => KN_VARTYPE_BINARY
)

# This list was obtained through AWK with Knitro 11.1 with the following command:
# grep "#define" knitro.h | grep "KN_PARAM_" | awk '{ printf("\"%s\" => Cint(%s),\n",$2,$3) }'
const KN_paramName2Indx = Dict(
    "KN_PARAM_NEWPOINT" => Cint(1001),
    "KN_PARAM_HONORBNDS" => Cint(1002),
    "KN_PARAM_ALGORITHM" => Cint(1003),
    "KN_PARAM_ALG" => Cint(1003),
    "KN_PARAM_BAR_MURULE" => Cint(1004),
    "KN_PARAM_BAR_FEASIBLE" => Cint(1006),
    "KN_PARAM_GRADOPT" => Cint(1007),
    "KN_PARAM_HESSOPT" => Cint(1008),
    "KN_PARAM_BAR_INITPT" => Cint(1009),
    "KN_PARAM_ACT_LPSOLVER" => Cint(1012),
    "KN_PARAM_CG_MAXIT" => Cint(1013),
    "KN_PARAM_MAXIT" => Cint(1014),
    "KN_PARAM_OUTLEV" => Cint(1015),
    "KN_PARAM_OUTMODE" => Cint(1016),
    "KN_PARAM_SCALE" => Cint(1017),
    "KN_PARAM_SOC" => Cint(1019),
    "KN_PARAM_DELTA" => Cint(1020),
    "KN_PARAM_BAR_FEASMODETOL" => Cint(1021),
    "KN_PARAM_FEASTOL" => Cint(1022),
    "KN_PARAM_FEASTOLABS" => Cint(1023),
    "KN_PARAM_MAXTIMECPU" => Cint(1024),
    "KN_PARAM_BAR_INITMU" => Cint(1025),
    "KN_PARAM_OBJRANGE" => Cint(1026),
    "KN_PARAM_OPTTOL" => Cint(1027),
    "KN_PARAM_OPTTOLABS" => Cint(1028),
    "KN_PARAM_LINSOLVER_PIVOTTOL" => Cint(1029),
    "KN_PARAM_XTOL" => Cint(1030),
    "KN_PARAM_DEBUG" => Cint(1031),
    "KN_PARAM_MULTISTART" => Cint(1033),
    "KN_PARAM_MSENABLE" => Cint(1033),
    "KN_PARAM_MSMAXSOLVES" => Cint(1034),
    "KN_PARAM_MSMAXBNDRANGE" => Cint(1035),
    "KN_PARAM_MSMAXTIMECPU" => Cint(1036),
    "KN_PARAM_MSMAXTIMEREAL" => Cint(1037),
    "KN_PARAM_LMSIZE" => Cint(1038),
    "KN_PARAM_BAR_MAXCROSSIT" => Cint(1039),
    "KN_PARAM_MAXTIMEREAL" => Cint(1040),
    "KN_PARAM_CG_PRECOND" => Cint(1041),
    "KN_PARAM_BLASOPTION" => Cint(1042),
    "KN_PARAM_BAR_MAXREFACTOR" => Cint(1043),
    "KN_PARAM_LINESEARCH_MAXTRIALS" => Cint(1044),
    "KN_PARAM_BLASOPTIONLIB" => Cint(1045),
    "KN_PARAM_OUTAPPEND" => Cint(1046),
    "KN_PARAM_OUTDIR" => Cint(1047),
    "KN_PARAM_CPLEXLIB" => Cint(1048),
    "KN_PARAM_BAR_PENRULE" => Cint(1049),
    "KN_PARAM_BAR_PENCONS" => Cint(1050),
    "KN_PARAM_MSNUMTOSAVE" => Cint(1051),
    "KN_PARAM_MSSAVETOL" => Cint(1052),
    "KN_PARAM_PRESOLVEDEBUG" => Cint(1053),
    "KN_PARAM_MSTERMINATE" => Cint(1054),
    "KN_PARAM_MSSTARTPTRANGE" => Cint(1055),
    "KN_PARAM_INFEASTOL" => Cint(1056),
    "KN_PARAM_LINSOLVER" => Cint(1057),
    "KN_PARAM_BAR_DIRECTINTERVAL" => Cint(1058),
    "KN_PARAM_PRESOLVE" => Cint(1059),
    "KN_PARAM_PRESOLVE_TOL" => Cint(1060),
    "KN_PARAM_BAR_SWITCHRULE" => Cint(1061),
    "KN_PARAM_HESSIAN_NO_F" => Cint(1062),
    "KN_PARAM_MA_TERMINATE" => Cint(1063),
    "KN_PARAM_MA_MAXTIMECPU" => Cint(1064),
    "KN_PARAM_MA_MAXTIMEREAL" => Cint(1065),
    "KN_PARAM_MSSEED" => Cint(1066),
    "KN_PARAM_MA_OUTSUB" => Cint(1067),
    "KN_PARAM_MS_OUTSUB" => Cint(1068),
    "KN_PARAM_XPRESSLIB" => Cint(1069),
    "KN_PARAM_TUNER" => Cint(1070),
    "KN_PARAM_TUNER_OPTIONSFILE" => Cint(1071),
    "KN_PARAM_TUNER_MAXTIMECPU" => Cint(1072),
    "KN_PARAM_TUNER_MAXTIMEREAL" => Cint(1073),
    "KN_PARAM_TUNER_OUTSUB" => Cint(1074),
    "KN_PARAM_TUNER_TERMINATE" => Cint(1075),
    "KN_PARAM_LINSOLVER_OOC" => Cint(1076),
    "KN_PARAM_BAR_RELAXCONS" => Cint(1077),
    "KN_PARAM_MSDETERMINISTIC" => Cint(1078),
    "KN_PARAM_BAR_REFINEMENT" => Cint(1079),
    "KN_PARAM_DERIVCHECK" => Cint(1080),
    "KN_PARAM_DERIVCHECK_TYPE" => Cint(1081),
    "KN_PARAM_DERIVCHECK_TOL" => Cint(1082),
    "KN_PARAM_LINSOLVER_INEXACT" => Cint(1083),
    "KN_PARAM_LINSOLVER_INEXACTTOL" => Cint(1084),
    "KN_PARAM_MAXFEVALS" => Cint(1085),
    "KN_PARAM_FSTOPVAL" => Cint(1086),
    "KN_PARAM_DATACHECK" => Cint(1087),
    "KN_PARAM_DERIVCHECK_TERMINATE" => Cint(1088),
    "KN_PARAM_BAR_WATCHDOG" => Cint(1089),
    "KN_PARAM_FTOL" => Cint(1090),
    "KN_PARAM_FTOL_ITERS" => Cint(1091),
    "KN_PARAM_ACT_QPALG" => Cint(1092),
    "KN_PARAM_BAR_INITPI_MPEC" => Cint(1093),
    "KN_PARAM_XTOL_ITERS" => Cint(1094),
    "KN_PARAM_LINESEARCH" => Cint(1095),
    "KN_PARAM_OUT_CSVINFO" => Cint(1096),
    "KN_PARAM_INITPENALTY" => Cint(1097),
    "KN_PARAM_ACT_LPFEASTOL" => Cint(1098),
    "KN_PARAM_CG_STOPTOL" => Cint(1099),
    "KN_PARAM_RESTARTS" => Cint(1100),
    "KN_PARAM_RESTARTS_MAXIT" => Cint(1101),
    "KN_PARAM_BAR_SLACKBOUNDPUSH" => Cint(1102),
    "KN_PARAM_CG_PMEM" => Cint(1103),
    "KN_PARAM_BAR_SWITCHOBJ" => Cint(1104),
    "KN_PARAM_OUTNAME" => Cint(1105),
    "KN_PARAM_OUT_CSVNAME" => Cint(1106),
    "KN_PARAM_ACT_PARAMETRIC" => Cint(1107),
    "KN_PARAM_ACT_LPDUMPMPS" => Cint(1108),
    "KN_PARAM_ACT_LPALG" => Cint(1109),
    "KN_PARAM_ACT_LPPRESOLVE" => Cint(1110),
    "KN_PARAM_ACT_LPPENALTY" => Cint(1111),
    "KN_PARAM_BNDRANGE" => Cint(1112),
    "KN_PARAM_BAR_CONIC_ENABLE" => Cint(1113),
    "KN_PARAM_CONVEX" => Cint(1114),
    "KN_PARAM_OUT_HINTS" => Cint(1115),
    "KN_PARAM_EVAL_FCGA" => Cint(1116),
    "KN_PARAM_BAR_MAXCORRECTORS" => Cint(1117),
    "KN_PARAM_STRAT_WARM_START" => Cint(1118),
    "KN_PARAM_FINDIFF_TERMINATE" => Cint(1119),
    "KN_PARAM_CPUPLATFORM" => Cint(1120),
    "KN_PARAM_PRESOLVE_PASSES" => Cint(1121),
    "KN_PARAM_PRESOLVE_LEVEL" => Cint(1122),
    "KN_PARAM_FINDIFF_RELSTEPSIZE" => Cint(1123),
    "KN_PARAM_INFEASTOL_ITERS" => Cint(1124),
    "KN_PARAM_PRESOLVEOP_TIGHTEN" => Cint(1125),
    "KN_PARAM_BAR_LINSYS" => Cint(1126),
    "KN_PARAM_PRESOLVE_INITPT" => Cint(1127),
    "KN_PARAM_ACT_QPPENALTY" => Cint(1128),
    "KN_PARAM_BAR_LINSYS_STORAGE" => Cint(1129),
    "KN_PARAM_LINSOLVER_MAXITREF" => Cint(1130),
    "KN_PARAM_BFGS_SCALING" => Cint(1131),
    "KN_PARAM_MIP_METHOD" => Cint(2001),
    "KN_PARAM_MIP_BRANCHRULE" => Cint(2002),
    "KN_PARAM_MIP_SELECTRULE" => Cint(2003),
    "KN_PARAM_MIP_INTGAPABS" => Cint(2004),
    "KN_PARAM_MIP_INTGAPREL" => Cint(2005),
    "KN_PARAM_MIP_MAXTIMECPU" => Cint(2006),
    "KN_PARAM_MIP_MAXTIMEREAL" => Cint(2007),
    "KN_PARAM_MIP_MAXSOLVES" => Cint(2008),
    "KN_PARAM_MIP_INTEGERTOL" => Cint(2009),
    "KN_PARAM_MIP_OUTLEVEL" => Cint(2010),
    "KN_PARAM_MIP_OUTINTERVAL" => Cint(2011),
    "KN_PARAM_MIP_OUTSUB" => Cint(2012),
    "KN_PARAM_MIP_DEBUG" => Cint(2013),
    "KN_PARAM_MIP_IMPLICATNS" => Cint(2014),
    "KN_PARAM_MIP_GUB_BRANCH" => Cint(2015),
    "KN_PARAM_MIP_KNAPSACK" => Cint(2016),
    "KN_PARAM_MIP_ROUNDING" => Cint(2017),
    "KN_PARAM_MIP_ROOTALG" => Cint(2018),
    "KN_PARAM_MIP_LPALG" => Cint(2019),
    "KN_PARAM_MIP_TERMINATE" => Cint(2020),
    "KN_PARAM_MIP_MAXNODES" => Cint(2021),
    "KN_PARAM_MIP_HEURISTIC" => Cint(2022),
    "KN_PARAM_MIP_HEUR_MAXIT" => Cint(2023),
    "KN_PARAM_MIP_HEUR_MAXTIMECPU" => Cint(2024),
    "KN_PARAM_MIP_HEUR_MAXTIMEREAL" => Cint(2025),
    "KN_PARAM_MIP_PSEUDOINIT" => Cint(2026),
    "KN_PARAM_MIP_STRONG_MAXIT" => Cint(2027),
    "KN_PARAM_MIP_STRONG_CANDLIM" => Cint(2028),
    "KN_PARAM_MIP_STRONG_LEVEL" => Cint(2029),
    "KN_PARAM_MIP_INTVAR_STRATEGY" => Cint(2030),
    "KN_PARAM_MIP_RELAXABLE" => Cint(2031),
    "KN_PARAM_MIP_NODEALG" => Cint(2032),
    "KN_PARAM_MIP_HEUR_TERMINATE" => Cint(2033),
    "KN_PARAM_MIP_SELECTDIR" => Cint(2034),
    "KN_PARAM_MIP_CUTFACTOR" => Cint(2035),
    "KN_PARAM_MIP_ZEROHALF" => Cint(2036),
    "KN_PARAM_MIP_MIR" => Cint(2037),
    "KN_PARAM_MIP_CLIQUE" => Cint(2038),
    "KN_PARAM_PAR_NUMTHREADS" => Cint(3001),
    "KN_PARAM_PAR_CONCURRENT_EVALS" => Cint(3002),
    "KN_PARAM_PAR_BLASNUMTHREADS" => Cint(3003),
    "KN_PARAM_PAR_LSNUMTHREADS" => Cint(3004),
    "KN_PARAM_PAR_MSNUMTHREADS" => Cint(3005),
    "KN_PARAM_PAR_CONICNUMTHREADS" => Cint(3006),
)

const KNITRO_OPTIONS = String[
    "newpoint",                  # KN_PARAM_NEWPOINT=            #
    "honorbnds",                 # KN_PARAM_HONORBNDS            #
    "algorithm",                 # KN_PARAM_ALGORITHM            #
    "bar_murule",                # KN_PARAM_BAR_MURULE           #
    "bar_feasible",              # KN_PARAM_BAR_FEASIBLE         #
    "gradopt",                   # KN_PARAM_GRADOPT              #
    "hessopt",                   # KN_PARAM_HESSOPT              #
    "bar_initpt",                # KN_PARAM_BAR_INITPT           #
    "act_lpsolver",              # KN_PARAM_ACT_LPSOLVER         #
    "cg_maxit",                  # KN_PARAM_CG_MAXIT             #
    "maxit",                     # KN_PARAM_MAXIT                #
    "outlev",                    # KN_PARAM_OUTLEV               #
    "outmode",                   # KN_PARAM_OUTMODE              #
    "scale",                     # KN_PARAM_SCALE                #
    "soc",                       # KN_PARAM_SOC                  #
    "delta",                     # KN_PARAM_DELTA                #
    "bar_feasmodetol",           # KN_PARAM_BAR_FEASMODETOL      #
    "feastol",                   # KN_PARAM_FEASTOL              #
    "feastolabs",                # KN_PARAM_FEASTOLABS           #
    "maxtimecpu",                # KN_PARAM_MAXTIMECPU           #
    "bar_initmu",                # KN_PARAM_BAR_INITMU           #
    "objrange",                  # KN_PARAM_OBJRANGE             #
    "opttol",                    # KN_PARAM_OPTTOL               #
    "opttolabs",                 # KN_PARAM_OPTTOLABS            #
    "linsolver_pivottol",        # KN_PARAM_LINSOLVER_PIVOTTOL   #
    "xtol",                      # KN_PARAM_XTOL                 #
    "debug",                     # KN_PARAM_DEBUG                #
    "ms_enable",                 # KN_PARAM_MULTISTART           #
    "ms_maxsolves",              # KN_PARAM_MSMAXSOLVES          #
    "ms_maxbndrange",            # KN_PARAM_MSMAXBNDRANGE        #
    "ms_maxtime_cpu",            # KN_PARAM_MSMAXTIMECPU         #
    "ms_maxtime_real",           # KN_PARAM_MSMAXTIMEREAL        #
    "lmsize",                    # KN_PARAM_LMSIZE               #
    "bar_maxcrossit",            # KN_PARAM_BAR_MAXCROSSIT       #
    "maxtime_real",              # KN_PARAM_MAXTIMEREAL          #
    "cg_precond",                # KN_PARAM_CG_PRECOND           #
    "blasoption",                # KN_PARAM_BLASOPTION           #
    "bar_maxrefactor",           # KN_PARAM_BAR_MAXREFACTOR      #
    "linesearch_maxtrials",      # KN_PARAM_LINESEARCH_MAXTRIALS #
    "blasoptionlib",             # KN_PARAM_BLASOPTIONLIB        #
    "outappend",                 # KN_PARAM_OUTAPPEND            #
    "outdir",                    # KN_PARAM_OUTDIR               #
    "cplexlibname",              # KN_PARAM_CPLEXLIB             #
    "bar_penaltyrule",           # KN_PARAM_BAR_PENRULE          #
    "bar_penaltycons",           # KN_PARAM_BAR_PENCONS          #
    "ms_num_to_save",            # KN_PARAM_MSNUMTOSAVE          #
    "ms_savetol",                # KN_PARAM_MSSAVETOL            #
    "ms_terminate",              # KN_PARAM_MSTERMINATE          #
    "ms_startptrange",           # KN_PARAM_MSSTARTPTRANGE       #
    "infeastol",                 # KN_PARAM_INFEASTOL            #
    "linsolver",                 # KN_PARAM_LINSOLVER            #
    "bar_directinterval",        # KN_PARAM_BAR_DIRECTINTERVAL   #
    "presolve",                  # KN_PARAM_PRESOLVE             #
    "presolve_tol",              # KN_PARAM_PRESOLVE_TOL         #
    "bar_switchrule",            # KN_PARAM_BAR_SWITCHRULE       #
    "hessian_no_f",              # KN_PARAM_HESSIAN_NO_F         #
    "ma_terminate",              # KN_PARAM_MA_TERMINATE         #
    "ma_maxtime_cpu",            # KN_PARAM_MA_MAXTIMECPU        #
    "ma_maxtime_real",           # KN_PARAM_MA_MAXTIMEREAL       #
    "ms_seed",                   # KN_PARAM_MSSEED               #
    "ma_outsub",                 # KN_PARAM_MA_OUTSUB            #
    "ms_outsub",                 # KN_PARAM_MS_OUTSUB            #
    "xpresslibname",             # KN_PARAM_XPRESSLIB            #
    "tuner",                     # KN_PARAM_TUNER                #
    "tuner_optionsfile",         # KN_PARAM_TUNER_OPTIONSFILE    #
    "tuner_maxtime_cpu",         # KN_PARAM_TUNER_MAXTIMECPU     #
    "tuner_maxtime_real",        # KN_PARAM_TUNER_MAXTIMEREAL    #
    "tuner_outsub",              # KN_PARAM_TUNER_OUTSUB         #
    "tuner_terminate",           # KN_PARAM_TUNER_TERMINATE      #
    "linsolver_ooc",             # KN_PARAM_LINSOLVER_OOC        #
    "bar_relaxcons",             # KN_PARAM_BAR_RELAXCONS        #
    "ms_deterministic",          # KN_PARAM_MSDETERMINISTIC      #
    "bar_refinement",            # KN_PARAM_BAR_REFINEMENT       #
    "derivcheck",                # KN_PARAM_DERIVCHECK           #
    "derivcheck_type",           # KN_PARAM_DERIVCHECK_TYPE      #
    "derivcheck_tol",            # KN_PARAM_DERIVCHECK_TOL       #
    "maxfevals",                 # KN_PARAM_MAXFEVALS            #
    "fstopval",                  # KN_PARAM_FSTOPVAL             #
    "datacheck",                 # KN_PARAM_DATACHECK            #
    "derivcheck_terminate",      # KN_PARAM_DERIVCHECK_TERMINATE #
    "bar_watchdog",              # KN_PARAM_BAR_WATCHDOG         #
    "ftol",                      # KN_PARAM_FTOL                 #
    "ftol_iters",                # KN_PARAM_FTOL_ITERS           #
    "act_qpalg",                 # KN_PARAM_ACT_QPALG            #
    "bar_initpi_mpec",           # KN_PARAM_BAR_INITPI_MPEC      #
    "xtol_iters",                # KN_PARAM_XTOL_ITERS           #
    "linesearch",                # KN_PARAM_LINESEARCH           #
    "out_csvinfo",               # KN_PARAM_OUT_CSVINFO          #
    "initpenalty",               # KN_PARAM_INITPENALTY          #
    "act_lpfeastol",             # KN_PARAM_ACT_LPFEASTOL        #
    "cg_stoptol",                # KN_PARAM_CG_STOPTOL           #
    "restarts",                  # KN_PARAM_RESTARTS             #
    "restarts_maxit",            # KN_PARAM_RESTARTS_MAXIT       #
    "bar_slackboundpush",        # KN_PARAM_BAR_SLACKBOUNDPUSH   #
    "cg_pmem",                   # KN_PARAM_CG_PMEM              #
    "bar_switchobj",             # KN_PARAM_BAR_SWITCHOBJ        #
    "outname",                   # KN_PARAM_OUTNAME              #
    "out_csvname",               # KN_PARAM_OUT_CSVNAME          #
    "act_parametric",            # KN_PARAM_ACT_PARAMETRIC       #
    "act_lpdumpmps",             # KN_PARAM_ACT_LPDUMPMPS        #
    "act_lpalg",                 # KN_PARAM_ACT_LPALG            #
    "act_lppresolve",            # KN_PARAM_ACT_LPPRESOLVE       #
    "act_lppenalty",             # KN_PARAM_ACT_LPPENALTY        #
    "bndrange",                  # KN_PARAM_BNDRANGE             #
    "bar_conic_enable",          # KN_PARAM_BAR_CONIC_ENABLE     #
    "convex",                    # KN_PARAM_CONVEX               #
    "out_hints",                 # KN_PARAM_OUT_HINTS            #
    "eval_fcga",                 # KN_PARAM_EVAL_FCGA            #
    "bar_maxcorrectors",         # KN_PARAM_BAR_MAXCORRECTORS    #
    "strat_warm_start",          # KN_PARAM_STRAT_WARM_START     #
    "findiff_terminate",         # KN_PARAM_FINDIFF_TERMINATE    #
    "cpuplatform",               # KN_PARAM_CPUPLATFORM          #
    "presolve_passes",           # KN_PARAM_PRESOLVE_PASSES      #
    "presolve_level",            # KN_PARAM_PRESOLVE_LEVEL       #
    "findiff_relstepsize",       # KN_PARAM_FINDIFF_RELSTEPSIZE  #
    "infeastol_iters",           # KN_PARAM_INFEASTOL_ITERS      #
    "mip_method",                # KN_PARAM_MIP_METHOD           #
    "mip_branchrule",            # KN_PARAM_MIP_BRANCHRULE       #
    "mip_selectrule",            # KN_PARAM_MIP_SELECTRULE       #
    "mip_integral_gap_abs",      # KN_PARAM_MIP_INTGAPABS        #
    "mip_integral_gap_rel",      # KN_PARAM_MIP_INTGAPREL        #
    "mip_maxtimecpu",            # KN_PARAM_MIP_MAXTIMECPU       #
    "mip_maxtimereal",           # KN_PARAM_MIP_MAXTIMEREAL      #
    "mip_maxsolves",             # KN_PARAM_MIP_MAXSOLVES        #
    "mip_integer_tol",           # KN_PARAM_MIP_INTEGERTOL       #
    "mip_outlevel",              # KN_PARAM_MIP_OUTLEVEL         #
    "mip_outinterval",           # KN_PARAM_MIP_OUTINTERVAL      #
    "mip_outsub",                # KN_PARAM_MIP_OUTSUB           #
    "mip_debug",                 # KN_PARAM_MIP_DEBUG            #
    "mip_implications",          # KN_PARAM_MIP_IMPLICATNS       #
    "mip_gub_branch",            # KN_PARAM_MIP_GUB_BRANCH       #
    "mip_knapsack",              # KN_PARAM_MIP_KNAPSACK         #
    "mip_rounding",              # KN_PARAM_MIP_ROUNDING         #
    "mip_rootalg",               # KN_PARAM_MIP_ROOTALG          #
    "mip_lpalg",                 # KN_PARAM_MIP_LPALG            #
    "mip_terminate",             # KN_PARAM_MIP_TERMINATE        #
    "mip_maxnodes",              # KN_PARAM_MIP_MAXNODES         #
    "mip_heuristic",             # KN_PARAM_MIP_HEURISTIC        #
    "mip_heuristic_maxit",       # KN_PARAM_MIP_HEUR_MAXIT       #
    "mip_heuristic_maxtimecpu",  # KN_PARAM_MIP_HEUR_MAXTIMECPU  #
    "mip_heuristic_maxtimereal", # KN_PARAM_MIP_HEUR_MAXTIMEREAL #
    "mip_pseudoinit",            # KN_PARAM_MIP_PSEUDOINIT       #
    "mip_strong_maxit",          # KN_PARAM_MIP_STRONG_MAXIT     #
    "mip_strong_candlim",        # KN_PARAM_MIP_STRONG_CANDLIM   #
    "mip_strong_level",          # KN_PARAM_MIP_STRONG_LEVEL     #
    "mip_intvar_strategy",       # KN_PARAM_MIP_INTVAR_STRATEGY  #
    "mip_relaxable",             # KN_PARAM_MIP_RELAXABLE        #
    "mip_nodealg",               # KN_PARAM_MIP_NODEALG          #
    "mip_heuristic_terminate",   # KN_PARAM_MIP_HEUR_TERMINATE   #
    "mip_selectdir",             # KN_PARAM_MIP_SELECTDIR        #
    "mip_cutfactor",             # KN_PARAM_MIP_CUTFACTOR        #
    "mip_zerohalf",              # KN_PARAM_MIP_ZEROHALF         #
    "mip_mir",                   # KN_PARAM_MIP_MIR              #
    "mip_clique",                # KN_PARAM_MIP_CLIQUE           #
    "par_numthreads",            # KN_PARAM_PAR_NUMTHREADS       #
    "par_concurrent_evals",      # KN_PARAM_PAR_CONCURRENT_EVALS #
    "par_blasnumthreads",        # KN_PARAM_PAR_BLASNUMTHREADS   #
    "par_lsnumthreads",          # KN_PARAM_PAR_LSNUMTHREADS     #
    "par_msnumthreads",          # KN_PARAM_PAR_MSNUMTHREADS     #
    "par_conicnumthreads",       # KN_PARAM_PAR_CONICNUMTHREADS  #
    "findiff_relstepsize",       # KN_PARAM_FINDIFF_RELSTEPSIZE  #
    "infeastol_iters",           # KN_PARAM_INFEASTOL_ITERS      #
    "presolveop_tighten",        # KN_PARAM_PRESOLVEOP_TIGHTEN   #
    "pre_redundancylevel",       # KN_CNT_REDUNDANCY_DETECTION   #
    "pre_improvecoefficients",   # KN_IMPROVE_COEFFICIENTS       #
    "bar_linsys",                # KN_PARAM_BAR_LINSYS           #
    "presolve_initpt",           # KN_PARAM_PRESOLVE_INITPT      #
    "act_qppenalty",             # KN_PARAM_ACT_QPPENALTY        #
    "bar_linsys_storage",        # KN_PARAM_BAR_LINSYS_STORAGE   #
    "linsolver_maxitref",        # KN_PARAM_LINSOLVER_MAXITREF   #
    "bfgs_scaling",              # KN_PARAM_BFGS_SCALING         #
    "option_file",
    "tuner_file",
]


text:
"""
    legendre(x, p::Integer)

Evaluate the Legendre polynomial of degree `p` at `x` using the three term
recursion [Kopriva, Implementing Spectral Methods for PDEs, Algorithm 20].
"""
function legendre(x, p::Integer)
    a::typeof(x) = one(x)
    b::typeof(x) = x

    if p <= 0
        return a
    elseif p == 1
        return b
    end

    for j in 2:p
        a, b = b, ( (2j-1)*x*b - (j-1)*a ) / j
    end

    b
end

"""
    legendre_and_derivative(x, p::Integer)

Evaluate the Legendre polynomial of degree `p` and its derivative at `x` using
the three term recursion [Kopriva, Implementing Spectral Methods for PDEs,
Algorithm 22].
"""
function legendre_and_derivative(x, p::Integer)
    # coefficients for the polynomial...
    a::typeof(x) = one(x)
    b::typeof(x) = x
    # ... and for the derivative
    aa::typeof(x) = zero(x)
    bb::typeof(x) = one(x)

    if p <= 0
        return a, aa
    elseif p == 1
        return b, bb
    end

    for j in 2:p
        a, b = b, ( (2j-1)*x*b - (j-1)*a ) / j
        aa, bb = bb, aa + (2j-1)*a
    end

    b, bb
end

"""
    gauss_legendre_nodes_and_weights(p, T=Float64::Type, tol=4*eps(T), maxit=100)

Compute the Gauss-Legendre nodes and weights for polynomials of degree `p`
using the scalar type `T`, tolerance `tol` and maximal number of Newton iterations
`maxit` [Kopriva, Implementing Spectral Methods for PDEs, Algorithm 23].
"""
function gauss_legendre_nodes_and_weights(p, T=Float64::Type, tol=4*eps(T), maxit=100)
    if p <= 0
        return T[0], T[2]
    elseif p == 1
        x0 = -1 / sqrt(T(3))
        return [x0, -x0], T[1, 1]
    end

    nodes = Vector{T}(undef, p+1)
    weights = Vector{T}(undef, p+1)
    for j in 0:((p+1)÷2-1)
        x = -cospi(T(2j+1)/(2p+2))
        for k in 1:maxit
            pol, der = legendre_and_derivative(x, p+1)
            Δ = -pol / der
            x = x + Δ
            abs(Δ) <= tol*abs(x) && break
        end
        pol, der = legendre_and_derivative(x, p+1)
        nodes[j+1] = x
        nodes[p-j+1] = -x
        weights[j+1] = weights[p-j+1] = 2 / ( (1-x^2)*der^2 )
    end

    if mod(p,2) == 0
        pol, der = legendre_and_derivative(zero(T), p+1)
        nodes[p÷2+1]   = 0
        weights[p÷2+1] = 2 / der^2
    end

    return nodes, weights
end

# helper function [Kopriva, Implementing Spectral Methods for PDEs, Algorithm 24]
function q_and_L_evaluation(x, p::Integer)
    a = one(x)
    b = x
    aa = zero(x)
    bb = one(x)

    for j in 2:p
        a, b = b, ( (2j-1)*x*b - (j-1)*a ) / j
        aa, bb = bb, aa + (2j-1)*a
    end
    pol = ( (2p+1)*x*b - p*a ) / (p+1)
    der = aa + (2p+1)*b

    pol-a, der-aa, b
end

"""
    lobatto_legendre_nodes_and_weights(p, T=Float64::Type, tol=4*eps(T), maxit=100)

Compute the Lobatto-Legendre nodes and weights for polynomials of degree `p`
using the scalar type `T`, tolerance `tol` and maximal number of Newton iterations
`maxit` [Kopriva, Implementing Spectral Methods for PDEs, Algorithm 25].
"""
function lobatto_legendre_nodes_and_weights(p, T=Float64::Type, tol=4*eps(T), maxit=100)
    if p <= 0
        return T[0], T[2]
    elseif p == 1
        return T[-1, 1], T[1, 1]
    end

    nodes = Vector{T}(undef, p+1)
    weights = Vector{T}(undef, p+1)

    nodes[1] = -1
    nodes[end] = 1
    weights[1] = weights[end] = T(2) / (p*(p+1))

    for j in 1:((p+1)÷2-1)
        x = -cos( T(j+1//4)*π/p - 3/(8p*T(j+1//4)*π) )
        for k in 1:maxit
            pol, der, leg = q_and_L_evaluation(x, p)
            Δ = -pol / der
            x = x + Δ
            abs(Δ) <= tol*abs(x) && break
        end
        pol, der, leg = q_and_L_evaluation(x, p)
        nodes[j+1] = x
        nodes[p-j+1] = -x
        weights[j+1] = weights[p-j+1] = 2 / ( p*(p+1)*leg^2 )
    end

    if mod(p,2) == 0
        pol, der, leg = q_and_L_evaluation(zero(T), p)
        nodes[p÷2+1]   = 0
        weights[p÷2+1] = 2 / ( p*(p+1)*leg^2 )
    end

    return nodes, weights
end

"""
    legendre_vandermonde(nodes)

Computes the Vandermonde matrix with respect to the Legendre polynomials and
the nodal basis on `nodes`.
The Vandermonde matrix `V` is the transformation matrix from the modal Legendre
basis to the nodal Lagrange basis associated with `nodes`.
"""
function legendre_vandermonde(nodes::AbstractVector)
    T = eltype(nodes)
    pp1 = length(nodes)
    V = Array{T}(undef, pp1, pp1)
    for j in 1:pp1, (i,x) in enumerate(nodes)
      V[i, j] = legendre(x, j-1)
    end
    V
end

legendre_vandermonde(basis::NodalBasis{Line}) = legendre_vandermonde(basis.nodes)

"""
    legendre_M(p, T=Float64)

Computes the diagonal mass matrix in the modal Legendre basis up to degree `p`
using the scalar type `T`.
"""
function legendre_M(p, T=Float64)
  Diagonal( T[2//(2n+1) for n in 0:p] )
end

"""
    legendre_D(p, T=Float64)

Computes the derivative matrix in the modal Legendre basis up to degree `p`
using the scalar type `T`.
"""
function legendre_D(p, T=Float64)
  D = fill(zero(T), p+1, p+1)
  if p >= 1
    D[1, 2] = 1
  end
  if p >= 2
    D[2, 3] = 3
  end
  for col in 4:p+1
    D[:, col] = D[:, col-2]
    D[col-1, col] = 2*col-3
  end
  D
end


text:
module PlotFunctions
	import  ..MasslistFunctions
	using PyPlot

	export massDefectPlot, bananaPlot, addClickToggle

	function massDefectPlot(masses, compositions, concentrations, colors, plotTitle, colorCodeTitle; dotSize = 10, maxMass = 450, maxDefect = 0.25, minConc = 0.02, sumformulas = false)
	  figure()

	  h2o = MasslistFunctions.createCompound(H=2,O=1, Hplus=0)
	  o = MasslistFunctions.createCompound(O=1, Hplus=0)

	  for i=1:length(masses)
	    m = masses[i]
	    adduct = h2o
	    if (inCompositions(compositions[:,i] + adduct[4], compositions))
	      m1 = m + adduct[1]
	      plot([m, m1], [m-round(m), m1 - round(m1)], color="lightblue", zorder=1)
	    end
	    adduct = o
	    if (inCompositions(compositions[:,i] + adduct[4], compositions))
	      m1 = m + adduct[1]
	      plot([m, m1], [m-round(m), m1 - round(m1)], color="red", zorder=1)
	    end
	    if sumformulas text(masses[i],masses[i]-round(masses[i]),sumFormulaStringFromCompositionArray(compositions[:,i]), color="grey", clip_on=true, verticalalignment="center", size=10, zorder=100) end
	  end

	  scatter(masses, masses-round.(masses),dotSize*log.(concentrations./minConc), colors, zorder=10, linewidths=0.5)
	  xlim(0,maxMass)
	  ylim(0,maxDefect)
	  cb=colorbar()
	  cb["ax"]["set_ylabel"](colorCodeTitle)
	  xlabel("Mass [amu]")
	  ylabel("Kendrick Mass Defect")
	  title(plotTitle)
	  grid("on")
	end

	function bananaPlot(xbins,ybins,meshdataXY;subplotAx=0)
	    println("Not implemented yet!!!")
	    #if subplotAx == 0
	    #    figure()
	    #    subplotAx = subplot(1,1,1)
	    #end
	    pcolormesh([DateTime(2018,05,05),DateTime(2018,05,06),DateTime(2018,05,07)],[1,7,9],[1 2 3; 4 5 6; 7 8 9][1:end-1,1:end-1])
	    #yscale("log")
	end

	function addClickToggle(plotAxis)
	    f = plotAxis.get_figure()
	    plotLines = plotAxis.get_lines()
	    legLines = plotAxis.get_legend().get_lines()
	    for ll in legLines
	        ll.set_picker(5)
	    end
	    linesDict = Dict(zip(legLines,plotLines))
	    function onLegendPick(event)
	        legLine = event.artist
	        plotLine = linesDict[legLine]
	        isVisible = plotLine.get_visible()
	        plotLine.set_visible(! isVisible)
	        if ! isVisible
	            legLine.set_alpha(1.0)
	        else
	            legLine.set_alpha(0.1)
	        end
	        f.canvas.draw()
	    end
	    f.canvas.mpl_connect("pick_event",onLegendPick)
	end
	function getVisibleTraces(axis)
		return [ axis.get_lines()[i].get_visible() for i in 1:length(axis.get_lines())]
	end

end

module CalibrationFunctions
	import ..InterpolationFunctions
	import PyPlot

	export generateCalibFactorTrace, generateBgTraces, interpolateBgTraces
	function getMeanOfQuantile(samples,quant)
	    q = quantile(samples,quant)
	    return mean(samples[samples.<q])
	end
	function generateBgTraces(times, traces; slices=10, quant=0.05)
	    dt = false
	    if typeof(times[1]) == DateTime
		times = Dates.datetime2unix.(times)
		dt = true
	    end
	    l=size(traces,1)
	    w=size(traces,2)
	    m=slices
	    bgtraces = Array{Float64}(size(traces,1), size(traces,2))
	    bgtimes = [mean(times[Int(floor((n-1)*l/m+1)):Int(floor(n*l/m)),1]) for n = 1:m]
	    for i=1:w
		bgtraces[:,i] = InterpolationFunctions.interpolate(times, bgtimes, [getMeanOfQuantile(traces[Int(floor((n-1)*l/m+1)):Int(floor(n*l/m)),i],0.08) for n = 1:m])
	    end
	    if dt
		bgtimes = Dates.unix2datetime.(bgtimes)
	    end
	    return bgtraces
	end

	function interpolateBgTraces(times, bgtimes, bgvalues)
	    if typeof(times[1]) == DateTime
		times = Dates.datetime2unix.(times)
	    end
	    if typeof(bgtimes[1]) == DateTime
		bgtimes = Dates.datetime2unix.(bgtimes)
	    end

	    nMasses = size(bgvalues,2)
	    bgtraces = Array{Float64}(size(times,1), nMasses)

	    for i=1:nMasses
		bgtraces[:,i] = InterpolationFunctions.interpolate(times, bgtimes, bgvalues[:,i])
	    end

	    return bgtraces
	end

	function generateCalibFactorTrace(traceTimes, calibTimes, calibValues, transitionTimes)
	    if typeof(traceTimes[1]) == DateTime
		println("Convertung traceTimes to unixtime")
		traceTimes = Dates.datetime2unix(traceTimes)
	    end
	    if typeof(calibTimes[1]) == DateTime
		println("Convertung calibTimes to unixtime")
		calibTimes = Dates.datetime2unix(calibTimes)
	    end
	    if typeof(transitionTimes[1]) == DateTime
		println("Convertung transitionTimes to unixtime")
		transitionTimes = Dates.datetime2unix(transitionTimes)
	    end
	  # Stuff might be unsorted
	  traceTimesSorted = sort(traceTimes)
	  sortedCalibIndices = sortperm(calibTimes)
	  calibTimesSorted = calibTimes[sortedCalibIndices]
	  calibValuesSorted = calibValues[sortedCalibIndices]
	  transitionTimesSorted = sort(transitionTimes)

	  piecewiseTimes = Array{Float64}(0)
	  piecewiseValues = Array{Float64}(0)
	  currentValue = 0
	  currentCalibIndex = 1
	  currentTransitionIndex = 1
	  currentPiecewiseIndex = 1

	  # Drop transitions before start of trace
	  while (transitionTimesSorted[currentTransitionIndex] < traceTimes[1]) && (currentTransitionIndex < length(transitionTimesSorted))
	    currentTransitionIndex += 1
	  end
	  # Drop transitions before first calibration
	  while (transitionTimesSorted[currentTransitionIndex] < calibTimes[1]) && (currentTransitionIndex < length(transitionTimesSorted))
	    currentTransitionIndex += 1
	  end
	  println("First transition within data range: $currentTransitionIndex - $(transitionTimesSorted[currentTransitionIndex])")

	  # If no calib was done before start of trace, insert copy of first calib
	  println("No Calib before start of data range found, copying first successive calib")
	  push!(piecewiseTimes, minimum([traceTimes[1] calibTimes[1]]))
	  push!(piecewiseValues, calibValuesSorted[1])
	  currentPiecewiseIndex += 1

	  # Iterate over calibs and transitions
	  while (currentCalibIndex <= length(calibTimesSorted))
	    if currentTransitionIndex < length(transitionTimesSorted)
	      println("CurrCalib: $currentCalibIndex ($(calibTimesSorted[currentCalibIndex])), CurrTransition: $currentTransitionIndex ($(transitionTimesSorted[currentTransitionIndex]))")
	    end
	    # define piecewise values depending on what comes next, calib or transition
	    if ((currentTransitionIndex > length(transitionTimesSorted)) || (calibTimesSorted[currentCalibIndex] < transitionTimesSorted[currentTransitionIndex])) && (currentCalibIndex <= length(calibTimesSorted))
	      # add one point for a calib
	      println("Adding calib point  $(calibValuesSorted[currentCalibIndex]) at $(calibTimesSorted[currentCalibIndex])")
	      push!(piecewiseTimes, calibTimesSorted[currentCalibIndex])
	      push!(piecewiseValues, calibValuesSorted[currentCalibIndex])
	      currentPiecewiseIndex += 1
	      currentCalibIndex += 1
	    else
	      # skip all following transitions that come before a calibration, since there are no calibs in between
	      while (currentTransitionIndex < length(transitionTimesSorted)) && (transitionTimesSorted[currentTransitionIndex+1] < calibTimesSorted[currentCalibIndex])
		println("Skipping transition $currentTransitionIndex, there are no calibs in between!")
		currentTransitionIndex +=1
	      end
	      # add one point with last value before transition
	      println("Adding pre-transition point $(piecewiseValues[end]) at $(transitionTimesSorted[currentTransitionIndex])")
	      push!(piecewiseTimes, transitionTimesSorted[currentTransitionIndex] - 1e-99)
	      push!(piecewiseValues, piecewiseValues[end])
	      currentPiecewiseIndex += 1

	      # add one point with next value after transition
	      println("Adding post-transition point $(calibValuesSorted[currentCalibIndex]) at $(transitionTimesSorted[currentTransitionIndex])")
	      push!(piecewiseTimes, transitionTimesSorted[currentTransitionIndex] + 1e-99)
	      push!(piecewiseValues, calibValuesSorted[currentCalibIndex])
	      currentPiecewiseIndex += 1

	      currentTransitionIndex += 1
	    end
	  end
	  PyPlot.figure()
	  PyPlot.plot(piecewiseTimes, piecewiseValues, "x-")

	  return InterpolationFunctions.interpolate(traceTimes, piecewiseTimes, piecewiseValues)
	end
end

module ExportFunctions
	using ..InterpolationFunctions
	using Dates
	import ..MasslistFunctions

	export exportTracesCSV, exportTracesCSVLossCorr, toMatlabTime, fromMatlabTime

	function exportTracesCSV(saveFolderPath, elementNames, compositions, times, traces; average=0)
	  sumformulas = MasslistFunctions.sumFormulaStringListFromCompositionArrayList(compositions)
	  f = open("$saveFolderPath/ptr3compositions.txt", "w")
	  writedlm(f, hcat(["Mass" "SumFormula"],reshape(elementNames,(1,length(elementNames)))))
	  writedlm(f, hcat(MasslistFunctions.massFromCompositionArrayList(compositions),sumformulas , compositions'))
	  close(f)
	  f = open("$saveFolderPath/ptr3traces.csv", "w")
	  writedlm(f, hcat(["Time"], reshape(sumformulas,(1,length(sumformulas)))))
	  if (average==0)
	    writedlm(f, hcat(times ,traces))
	  else
	    writedlm(f, hcat(averageSamples(times,average) ,averageSamples(traces,average)))
	  end
	  close(f)
	end

	function exportTracesCSVLossCorr(saveFolderPath, elementNames, compositions, times, traces, lossfactor, lossfactorerr, corrfactor, corrfactorerr, corrnotes; average=0)
	  sumformulas = MasslistFunctions.sumFormulaStringListFromCompositionArrayList(compositions)
	  f = open("$saveFolderPath/ptr3compositions.txt", "w")
	  writedlm(f, hcat(["Mass"	"SumFormula"],reshape(elementNames,(1,length(elementNames))),["LossFactor"	"LossFactorError"	"CorrFactor"	"CorrFactorErr"	"CorrNotes"]))
	  writedlm(f, hcat(MasslistFunctions.massFromCompositionArrayList(compositions),sumformulas , compositions', lossfactor, lossfactorerr, corrfactor, corrfactorerr, corrnotes))
	  close(f)
	  f = open("$saveFolderPath/ptr3tracesInletLossCorr.csv", "w")
	  writedlm(f, hcat(["Time"], reshape(sumformulas,(1,length(sumformulas)))))
	  if (average==0)
	    writedlm(f, hcat(times , (corrfactor.*traces' )' ))
	  else
	    writedlm(f, hcat(averageSamples(times,average) ,(corrfactor.*(averageSamples(traces,average))' )' ))
	  end
	  close(f)
	end

	################ EXAMPLE from Matlab: 737551.673515479 should be 06-May-2019 16:09:51 #############
	function toMatlabTime(t::Dates.DateTime)
	    timespan = ((t+Dates.Day(1)) - Dates.DateTime(0,1,1,0,0,0))
	    millis::Float64 = Dates.value(timespan)
	    return millis/24/3600/1000
	end

	function fromMatlabTime(timestamp::Number)
	    days=Int(floor(timestamp))
	    millisecondsRemainder = Int(round((timestamp-days)*24*3600*1000))
	    return Dates.DateTime(0,1,1,0,0,0)+Dates.Day(days)+Dates.Millisecond(millisecondsRemainder)-Dates.Day(1)
	end

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -11.981958290664823 <= q <= 14.886364381032921)

                     add_NL_constraint(m, :((0.5346214578242647 + 0.7328398696189291*$(x[1]) + 0.5275606331718761*$(x[2]) + -0.05282758308498137*$(x[3]) + -0.6801317780174543*$(x[4]) + 0.6177845475163379*$(x[5]))*(1 + erf((0.5346214578242647 + 0.7328398696189291*$(x[1]) + 0.5275606331718761*$(x[2]) + -0.05282758308498137*$(x[3]) + -0.6801317780174543*$(x[4]) + 0.6177845475163379*$(x[5]))/sqrt(2)))/2 + (-0.36399156729020854 + 0.9651821996491203*$(x[1]) + -0.1221090890096228*$(x[2]) + -0.0037390408554469623*$(x[3]) + -0.25471778237955167*$(x[4]) + 0.004904328924779211*$(x[5]))*(1 + erf((-0.36399156729020854 + 0.9651821996491203*$(x[1]) + -0.1221090890096228*$(x[2]) + -0.0037390408554469623*$(x[3]) + -0.25471778237955167*$(x[4]) + 0.004904328924779211*$(x[5]))/sqrt(2)))/2 + (0.6792969956775154 + -0.9847763250624375*$(x[1]) + -0.7653951229931786*$(x[2]) + 0.7304491787436271*$(x[3]) + -0.81259603694894*$(x[4]) + -0.9953733056123499*$(x[5]))*(1 + erf((0.6792969956775154 + -0.9847763250624375*$(x[1]) + -0.7653951229931786*$(x[2]) + 0.7304491787436271*$(x[3]) + -0.81259603694894*$(x[4]) + -0.9953733056123499*$(x[5]))/sqrt(2)))/2 + (-0.022503806788797487 + -0.877492362258455*$(x[1]) + 0.9093714868360583*$(x[2]) + -0.36771946270068234*$(x[3]) + 0.3633871862886231*$(x[4]) + -0.6803298907372901*$(x[5]))*(1 + erf((-0.022503806788797487 + -0.877492362258455*$(x[1]) + 0.9093714868360583*$(x[2]) + -0.36771946270068234*$(x[3]) + 0.3633871862886231*$(x[4]) + -0.6803298907372901*$(x[5]))/sqrt(2)))/2 + (0.6247799657612747 + 0.2803258134848612*$(x[1]) + 0.15247610330367944*$(x[2]) + -0.6290763102374144*$(x[3]) + -0.40592776579339773*$(x[4]) + -0.5176681326197756*$(x[5]))*(1 + erf((0.6247799657612747 + 0.2803258134848612*$(x[1]) + 0.15247610330367944*$(x[2]) + -0.6290763102374144*$(x[3]) + -0.40592776579339773*$(x[4]) + -0.5176681326197756*$(x[5]))/sqrt(2)))/2 - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions