text:
# Similar to plotTrajectory.jl, but in 3D

using Plots

function plotTrajPos3D_Multiple(ptList, showN = 10)
    plt = plot3d(legend = :outerright, zlabel = "Z (m)")

    if length(ptList) >= showN
        indices = floor.(Int, LinRange(1, length(ptList), showN))
    else
        indices = collect(1:length(ptList))
    end

    println("Plotting the following indices: ")
    println(indices)

    minXY = 0
    maxXY = 0

    for ind in indices
        pt = ptList[ind]
        xyzList = splitDimensions(pt.sList)
        xList = xyzList[1]
        yList = xyzList[2]
        zList = xyzList[3]
        plot3d!(xList, yList, zList, markershape = :circle,
                                label = "Traj $ind")

        minXY = min(minXY, minimum(xList), minimum(yList))
        maxXY = max(maxXY, maximum(xList), maximum(yList))
    end

    xlabel!("X (m)")
    ylabel!("Y (m)")
    # zlabel!("Z (km)")
    title!("Trajectory Convergence")
    # yflip!()

    # Check that the XY plane looks square
    xlims!(minXY, maxXY)
    ylims!(minXY, maxXY)

    return plt
end

function plotSVUTime3D_Multiple(sListStack, vListStack, uListStack)

    # Variables needed for later
    nDim = size(sListStack[1][1], 1)
    sxyzLabels = ["x", "y", "z"]
    vxyzLabels = ["vx", "vy", "vz"]
    uxyzLabels = ["ux", "uy", "uz"]

    seLabels = ["Initial", "Final"]

    markerShapeArr = [:rect, :diamond]
    markerSizeArr = [8; 4]
    markerColorArr = [:blue, :purple, :darkgreen, :orange2, :gold, :pink]

    # Plot the Positions
    plt_s = plot()

    colorSelect = 1
    for (ind, sList) in enumerate(sListStack)
        xyzList = splitDimensions(sList)
        for d in 1:nDim
            thisLabel = sxyzLabels[d] * " for " * seLabels[ind] * " Traj"
            plot!(xyzList[d], markershape = markerShapeArr[ind],
                     markersize = markerSizeArr[ind],
                     markercolor = markerColorArr[colorSelect],
                     linecolor = markerColorArr[colorSelect],
                     label = thisLabel)
            colorSelect += 1
        end
    end
    title!("Change in position over time")
    xlabel!("Time")
    ylabel!("Position")
    display(plt_s)

    # Plot the Velocities
    plt_v = plot()
    colorSelect = 1
    for (ind, vList) in enumerate(vListStack)
        vxyzList = splitDimensions(vList)
        for d in 1:nDim
            thisLabel = vxyzLabels[d] * " for " * seLabels[ind] * " Traj"
            plot!(vxyzList[d], markershape = markerShapeArr[ind],
                     markersize = markerSizeArr[ind],
                     markercolor = markerColorArr[colorSelect],
                     linecolor = markerColorArr[colorSelect],
                     label = thisLabel)
            colorSelect += 1
        end
    end
    title!("Change in velocity over time")
    xlabel!("Time")
    ylabel!("Velocity")
    display(plt_v)

    # Plot the Controls
    colorSelect = 1
    markerColorArrU = [ markerColorArr[1], markerColorArr[4],
                        markerColorArr[2], markerColorArr[5],
                        markerColorArr[3], markerColorArr[6]]

    uxyzListSplit = splitDimensions.(uListStack)

    uPltList = []

    for d in 1:nDim
        pltNew = plot(legend = :outerright)
        for ind in 1:length(uxyzListSplit)
            thisLabel = uxyzLabels[d] * " for " * seLabels[ind] * " Traj"
            uxyzList = uxyzListSplit[ind][d]
            plot!(uxyzList, markershape = markerShapeArr[ind],
                     markersize = markerSizeArr[ind],
                     markercolor = markerColorArrU[colorSelect],
                     linecolor = markerColorArrU[colorSelect],
                     label = thisLabel)
            colorSelect += 1
        end
        push!(uPltList, pltNew)
        # display(pltNew)
    end

    plt_u = plot(uPltList..., layout = (length(uPltList), 1))

    title!("Change in thrust over time")
    xlabel!("Time")
    ylabel!("Contol (Thrust)")
    display(plt_u)

    return  plt_s, plt_v, plt_u
end

function plotSVUTime3D_StartEnd(ptList::Array{parseTrajectory, 1})
    ptStart = ptList[1]
    ptEnd = ptList[end]

    sListStack = [ptStart.sList, ptEnd.sList]
    vListStack = [ptStart.vList, ptEnd.vList]
    uListStack = [ptStart.uList, ptEnd.uList]

    plotSVUTime3D_Multiple(sListStack, vListStack, uListStack)
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
abstract type Optimizer end

# Define general use solver
function optimize!(opt::Optimizer, opts::Options)

    # Dispatch to _solve() with typeof opt
    res = _optimize!(opt, opts)

    #return results
    return res
end

function optimize!(opt::Optimizer)

    # Initialize options with defaults
    opts = Options()

    # Solve
    res = optimize!(opt, opts)

    # return results
    return res
end

ext: jl


text:
using PhomericsTest
using Test

@testset "Tools.jl" begin
    x = collect(0:5)
    y = x.^2

    @test first_derivative(x,y) == [1.00, 3.00, 5.00, 7.00, 9.00]
end


text:
module BNP

import MAT
import HDF5
import Distributions

export
Gaussian1D, Dirichlet,
Gaussian1DGaussian1D, MultinomialDirichlet,
BMM, LDA, DPM, HDP, dHDP,
sample, collapsed_gibbs_sampler, truncated_gibbs_sampler, posterior

include("common.jl")
include("distributions.jl")
include("conjugates.jl")
include("BMM.jl")
include("LDA.jl")
include("DPM.jl")
include("HDP.jl")
include("dHDP.jl")
include("generate_data.jl")

end # module


text:
################################################################################
#
# Run inversion for the Gaussian lens model
#
################################################################################

### Module loading

using LinearAlgebra, JLD2, PyPlot
using Optim
using JUDI.TimeModeling
push!(LOAD_PATH, string(pwd(), "/src/"))
using TWRIdual

### Load synthetic data

@load "./data/GaussLens/GaussLens_data.jld"

### Background model

m0 = 1f0./2f0.^2*ones(R, model_true.n)
n = model_true.n
d = model_true.d
o = model_true.o
model0 = Model(n, d, o, m0)

### Set objective functional

# Time sampling
dt_comp = dat.geometry.dt

# Pre- and post-conditioning
mask = BitArray(undef, model0.n); mask .= false
mask[6:end-5, 6:end-5] .= true
function proj_bounds(m, mmin, mmax)
    m[m[:] .< mmin] .= mmin
    m[m[:] .> mmax] .= mmax
    return m
end
mmin = 1f0/4f0^2
mmax = Inf
preproc(x) = proj_bounds(contr2abs(x, mask, m0), mmin, mmax)
postproc(g) = gradprec_contr2abs(g, mask, m0)

# # [FWI]
# @show inv_name = "FWI"
# fun!(F, G, x) = objFWI!(F, G, preproc(x), n, d, o, fsrc, dat; dt_comp = dt_comp, gradmprec_fun = postproc)

# [TWRIdual]
@show inv_name = "TWRIdual"
ε = Array{Float32, 1}(undef, fsrc.nsrc); ε .= 0.01f0
v_bg = sqrt(1/m0[1])
freq_peak = 0.006f0
δ = 0.1f0*R(sqrt(2)/2)*v_bg/freq_peak
weight_fun_pars = ("srcfocus", δ)
grad_corr = true
fun!(F, G, x) = objTWRIdual!(F, G, preproc(x), n, d, o, fsrc, dat, ε; comp_alpha = true, grad_corr = grad_corr, weight_fun_pars = weight_fun_pars, dt_comp = dt_comp, gradmprec_fun = postproc)

### Optimization

# Starting guess
x0 = zeros(R, length(findall(mask .== true)))

# Options
method = LBFGS()
niter = 20
optimopt = Optim.Options(iterations = niter, store_trace = true, show_trace = true, show_every = 1)

# Run
result = optimize(Optim.only_fg!(fun!), x0, method, optimopt)
x0 .= Optim.minimizer(result)
m_inv = preproc(x0)
loss_log = Optim.f_trace(result)


text:
for N in [Float64, Rational{Int}, Float32]
    # singleton
    S = Singleton(N[1, -2, 3, -4])
    H = Hyperrectangle(zeros(N, 4), N[1, 2, 3, 4])
    @test symmetric_interval_hull(S) == H

    # linear map of singleton
    M = Diagonal(ones(N, 4))
    @test symmetric_interval_hull(M * S) == H

    # Minkowski sum of singletons
    H2 = Hyperrectangle(zeros(N, 4), 2 * N[1, 2, 3, 4])
    @test symmetric_interval_hull(S + S) == H2

    # linear map of hyperrectangle
    M = Diagonal(N[1, -1, -2, 2])
    H1 = Hyperrectangle(5 * ones(N, 4), N[1, 2, 3, 4])
    H2 = Hyperrectangle(zeros(N, 4), N[6, 7, 16, 18])
    @test symmetric_interval_hull(M * H1) == H2
end


text:
using LinearAlgebra
using ToeplitzMatrices

function chebdif(N, M)

    # The function x,DM = chebdif(N,M) computes the differentiation
    # matrices D1, D2, ..., DM on Chebyshev nodes.
    # It is adapted from the function with the same name in dmsuite
    # developed for Matlab by JAC Weideman and SC Reddy.
    #
    # Adapted by M. Beneitez at KTH Mechanics.
    # beneitez@mech.kth.se
    #
    # Input:
    # N:  Size of the differentiation matrix.
    # M:  Number of derivates required
    # Note: 0 < M <= N-1.
    #
    # Output:
    # DM: DM[1:N,1:N,ell] contains the ell-th derivative matrix, ell=1...M
    #

    eye = Matrix{Float64}(I, N, N);        # Identity matrix.

    n1 = Int(floor(N/2)); n2  = Int(ceil(N/2));     # Indices used for flipping trick.

    k = collect(0:N-1);                        # Compute theta vector.
    th = k.*pi./(N-1)

    x = sin.(pi.*collect(N-1:-2:1-N)./(2*(N-1))); # Compute Chebyshev points.

    T = repeat(th/2,1,N);
    DX = 2*sin.(T'+T).*sin.(T'-T);          # Trigonometric identity.
    DX = [DX[1:n1,:]; -reverse(reverse(DX[1:n2,:], dims = 2), dims=1)];   # Flipping trick.
    DX[eye.==true] = ones(N,1);                      # Put 1's on the main diagonal of DX.

    C = Matrix(Toeplitz((-1).^vec(k),(-1).^vec(k)));            # C is the matrix with
    C[1,:] = C[1,:].*2; C[N,:] = C[N,:].*2;     # entries c[k]/c[j]
    C[:,1] = C[:,1]/2; C[:,N] = C[:,N]/2

    Z = 1 ./DX;                              # Z contains entries 1/(x[k]-x[j])
    Z[eye.==true] = zeros(N,1);              # with zeros on the diagonal.

    D  = eye;                                 # D contains diff(). matrices.
    DM = zeros(N,N,M);

    for ell = 1:M
        D = ell*Z.*(C.*repeat(diag(D),1,N) - D); # Off-diagonals
        D[eye.==true] = -sum(D', dims = 1);                            # Correct main diagonal of D
        DM[:,:,ell] = D;                                   # Store current D in DM
    end

    return x,DM

end


text:
const N_GRADS_2D_EXPONENT = 7
const N_GRADS_2D = 1 << N_GRADS_2D_EXPONENT

const NORMALIZER_2D = 0.05481866495625118

struct Grad2
    x::Float32
    y::Float32
end

const grad2 = Grad2[
    Grad2( 0.38268343236509,   0.923879532511287),
    Grad2( 0.923879532511287,  0.38268343236509),
    Grad2( 0.923879532511287, -0.38268343236509),
    Grad2( 0.38268343236509,  -0.923879532511287),
    Grad2(-0.38268343236509,  -0.923879532511287),
    Grad2(-0.923879532511287, -0.38268343236509),
    Grad2(-0.923879532511287,  0.38268343236509),
    Grad2(-0.38268343236509,   0.923879532511287),
    #-------------------------------------------#
    Grad2( 0.130526192220052,  0.99144486137381),
    Grad2( 0.608761429008721,  0.793353340291235),
    Grad2( 0.793353340291235,  0.608761429008721),
    Grad2( 0.99144486137381,   0.130526192220051),
    Grad2( 0.99144486137381,  -0.130526192220051),
    Grad2( 0.793353340291235, -0.60876142900872),
    Grad2( 0.608761429008721, -0.793353340291235),
    Grad2( 0.130526192220052, -0.99144486137381),
    Grad2(-0.130526192220052, -0.99144486137381),
    Grad2(-0.608761429008721, -0.793353340291235),
    Grad2(-0.793353340291235, -0.608761429008721),
    Grad2(-0.99144486137381,  -0.130526192220052),
    Grad2(-0.99144486137381,   0.130526192220051),
    Grad2(-0.793353340291235,  0.608761429008721),
    Grad2(-0.608761429008721,  0.793353340291235),
    Grad2(-0.130526192220052,  0.99144486137381),
]

const GRADIENTS_2D = Vector{Grad2}(undef, N_GRADS_2D)

function grad(a, seed, xsvp, ysvp, dx, dy)
    hash = xor(seed, xsvp, ysvp) * HASH_MULTIPLIER
    g = GRADIENTS_2D[(xor(hash, hash >> (64 - N_GRADS_2D_EXPONENT)) & (N_GRADS_2D - 1)) + 1]
    Float32(p4(a) * (g.x * dx + g.y * dy) / NORMALIZER_2D)
end

function fill_2d()
    j = 0
    len = length(grad2)
    for i = 1:N_GRADS_2D
        GRADIENTS_2D[i] = grad2[j += 1]
        j == len && (j = 0)
    end
end

fill_2d()


text:
##################################################

const ∑ = sum

function innerproduct(P::Type{<:AbstractDiscreteOrthogonalPolynomial}, f, g)
    dom = domain(P)
    fn = x -> f(x) * g(x) * weight_function(P)(x)
    a, b = first(dom), last(dom)
    if !isinf(a) && !isinf(b)
        return ∑(fn(x)  for x in first(dom):last(dom))
    else
        ## what to do if infinite
    end
end

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

abstract type AbstractDiscreteWeightFunction{T,X,N} <: AbstractDiscreteOrthogonalPolynomial{T,X} end
abstract type DiscreteWeightFunction{T,X,N} <:  AbstractDiscreteWeightFunction{T,X,N} end
export DiscreteWeightFunction

"""
    DiscreteWeightFunction

For a discrete measure, `dλ = ∑ wᵢ δ(x - xᵢ)`, specified through two
vectors, `xs` and `ws`, a collection of monic orthogonal polynomials is
produced through Darboux's formula for `α_n` and `β_n` using the
3-term recurrence defined by `π_{n+1} = (x-α_n)⋅π_n - β_n⋅π_{n-1}` (`An=1`, `Bn=-α_n`, `Cn=β_n`)
and the discrete Stieltjes method [Guatschi §3.1](https://www.cs.purdue.edu/homes/wxg/Madrid.pdf).

# Example

Discrete Chebyshev by its weight function (uniform  on 0,1,…,N-1)

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> const SP = SpecialPolynomials;

julia> N = 9
9

julia> xs, ws = collect(0:N-1), ones(N);   # w(x) = ∑ wⱼ⋅δ(x-xⱼ)

julia> SP.@register0 DWF DiscreteWeightFunction

julia> SP.@register_discrete_weight_function(DWF, xs, ws)

julia> [SP.Bn.(DWF, 0:N-1) SP.Cn.(DWF, 0:N-1)]
9×2 Matrix{Float64}:
 -4.0  9.0
 -4.0  6.66667
 -4.0  5.13333
 -4.0  4.62857
 -4.0  4.12698
 -4.0  3.53535
 -4.0  2.83217
 -4.0  2.01026
 -4.0  1.06667

julia> i,j = 3,4; ## check  that ∫pᵢpⱼdw  = 0    for i,j=3,4

julia> sum(basis(DWF,i)(x) *  basis(DWF,j)(x) * w for  (x,w) in zip(xs, ws))
5.684341886080802e-14

julia> ## Gogin, Hirvensalo (https://doi.org/10.1007/s10958-017-3410-8) characterization
       D(k,N,x) =  sum((-1)^l * binomial(k+l,k) * binomial(N-l,k-l) *  SP.generalized_binomial(x,l) for l in 0:k)
D (generic function with 1 method)

julia> x = variable()
Polynomials.Polynomial(x)

julia> ps,qs = [D(k,N-1,x)  for  k in 0:N-1], [basis(DWF, k)(x) for k  in 0:N-1];

julia> all(qs .* [p[end] for p  in ps] .≈ ps)
true
```

"""
DiscreteWeightFunction

basis_symbol(::Type{<:AbstractDiscreteWeightFunction}) = "W"

xs_ws(::Type{<:AbstractDiscreteWeightFunction}) = throw(ArgumentError("No default method"))

# (Gautschi](https://www.cs.purdue.edu/homes/wxg/Madrid.pdf), section 3.1
# compute α_n = <tπ_n,π_n>/<π_n,π_n>, β_n = <π_n,π_n>/<π_{n-1},π_{n-1}>
# wherer <p,q> = ∑_1^N w_k p(k) q(k)
function discrete_stieltjes(W::Type{<:AbstractDiscreteWeightFunction})

    xs,ws = xs_ws(W)

    N  =  length(xs)
    n = N

    # k = 0 case
    πk_1, πk = zeros(N), ones(N)

    βk = β0 = norm_k = sum(ws)/1 # <π_0, π_0> = ∑ w_k 1 ⋅ 1
    αk = α0 = ∑(ws[k] * xs[k] for k in eachindex(ws)) / norm_k
    αs = [α0]
    βs = [β0]

    for _ ∈ 1:(n-1)

        πk1 = (xs .- αk) .* πk - βk * πk_1 # use just computed αk, βk to find π_{k+1} = (x-αk)⋅π_k - βk⋅π_{k-1}
        norm_k1 =  ∑(ws[k] * πk1[k]^2  for k in eachindex(ws)) # <π_{k+1}, π_{k+1}>

        # Darboux
        # α_{k+1} = <x ⋅ π_{k+1}, π_{k+1}> / <π_{k+1}, π_{k+1}>,
        # β_{k+1} = <π_{k+1}, π_{k+1}> / <π_k, π_k>,
        αk1 = ∑(ws[k] * xs[k] * πk1[k]^2 for k in eachindex(ws)) / norm_k1
        βk1 = norm_k1/norm_k

        push!(αs, αk1)
        push!(βs, βk1)

        αk, βk = αk1, βk1
        πk_1, πk = πk, πk1
        norm_k = norm_k1

    end

    (-αs, βs)
end

An(::Type{W}, n) where {W <: AbstractDiscreteWeightFunction} = one(eltype(W))

Bn(::Type{W}, k::Int) where {W <:AbstractDiscreteWeightFunction} = discrete_stieltjes(W)[1][k+1]
Cn(::Type{W}, k::Int) where {W <:AbstractDiscreteWeightFunction} = discrete_stieltjes(W)[2][k+1]

Polynomials.domain(::Type{<:DiscreteWeightFunction}) = Polynomials.Interval(-Inf, Inf)

innerproduct(W::Type{<:AbstractDiscreteWeightFunction}, f, g) = ∑(wk * f(xk) * g(xk) for (xk, wk) in zip(xs_ws(W)...))


text:
using Optim
using LinearAlgebra: I
using LineSearches
import NLopt

function loss_fn(X, W)
    n, d = size(X)
    M = X * W
    R = X - M
    loss = 0.5 / n * sum(R .^ 2)
    G_loss = - 1.0 / n * transpose(X) * R
    return loss, G_loss
end

function h_fn(W)
    d, _ = size(W)
    eye = 1 * Matrix(I, d, d)
    M = eye + W .* W / d
    # FIXME E = np.linalg.matrix_power(M, d - 1)
    E = M ^ (d-1)
    h = sum(transpose(E) .* M) - d
    G_h = transpose(E) .* W * 2
    return h, G_h
end

function adj_fn(v)
    # FIXME length or size(v,1)?
    d = convert(Int, sqrt(length(v)/2))
    # NOTE: positive and negative part
    reshape(v[1:d*d] - v[d*d+1:end], d, d)
end

function f(v, X, rho, alpha)
    lambda1 = 0.1
    n, d = size(X)

    W = adj_fn(v)
    # X use 1
    loss, G_loss = loss_fn(X, W)

    # X use 2
    h, G_h = h_fn(W)

    # rho and alpha uses
    obj = loss + 0.5 * rho * h * h + alpha * h + lambda1 * sum(v)
    G_smooth = G_loss + (rho * h + alpha) * G_h

    g_obj = hcat(G_smooth .+ lambda1, - G_smooth .+ lambda1)
    return obj, g_obj
end

# this works
function my_opt_NLopt(_f, _g, lower, upper, w_est)
    d = convert(Int, sqrt(size(w_est, 1) / 2))
    function obj_fn(v, G)
        ret = _f(v)
        g = _g(v)
        copyto!(G, g)
        ret
    end

    opt = NLopt.Opt(:LD_LBFGS, 2*d*d)
    opt.lower_bounds = lower
    opt.upper_bounds = upper
    opt.min_objective = obj_fn

    res = NLopt.optimize(opt, w_est)
    res = res[2]
end

# NOT WORKING
function my_opt_Optim(_f, _g, lower, upper, w_est)
    sol = optimize(_f,
                   _g,
                   # (v)->f(v, X, rho, alpha)[2],
                   # (G, v) -> copyto!(G, f(v, X, rho, alpha)[2]),
                   lower, upper,
                   w_est,
                   # FIXME this does not throw errors, but gives NaN
                   # Fminbox(LBFGS(linesearch=LineSearches.BackTracking())),
                   # FIXME this throw errors
                   Fminbox(LBFGS()),
                   # LBFGS(),
                   inplace=false,
                   Optim.Options())
    Optim.minimizer(sol)
end

# FIXME should I check the exact number matching?
function notears(X, fix_upper=false)
    # max_iter=100
    max_iter=10
    h_tol=1e-8
    rho_max=1e+16
    w_threshold=0.3

    n, d = size(X)
    w_est = zeros(2*d*d)
    rho = 1.0
    alpha = 0.0
    h = Inf

    # bounds
    lower = zeros(2*d*d)
    # DEBUG
    if fix_upper
        upper = [if (i==j || (i==2 && j==1)) 0.1 else Inf end for _ in 1:2 for i in 1:d for j in 1:d]
    else
        upper = [if i==j 0. else Inf end for _ in 1:2 for i in 1:d for j in 1:d]
    end
    # DEBUG setting the direction
    for i in 1:max_iter
        w_new = nothing
        h_new = nothing
        @info "iter $(i) .."
        while rho < rho_max
            # @info "running optim .."
            _f = (v)->f(v, X, rho, alpha)[1]
            _g = (v)->f(v, X, rho, alpha)[2]

            # w_new = my_opt_Optim(_f, _g, lower, upper, w_est)
            w_new = my_opt_NLopt(_f, _g, lower, upper, w_est)
            h_new, _ = h_fn(adj_fn(w_new))

            # VERBOSE
            # @show loss_fn(X, adj_fn(w_new))[1]
            @show h_new

            if h_new > 0.25 * h
                rho *= 10
            else
                break
            end
        end
        w_est, h = w_new, h_new
        alpha += rho * h
        if h <= h_tol || rho >= rho_max
            break
        end
    end
    res = adj_fn(w_est)
    # DEBUG threshold or not, here or out
    # res[abs.(res) .< w_threshold] .= 0
    res
end

function graphical_metrics(W, Ŵ)
    # assuming W is dxd matrix
    d = size(W,1)

    g = DiGraph(W)
    ĝ = DiGraph(Ŵ)

    # this will make it 1 or 1
    mg = adjacency_matrix(g)
    mĝ = adjacency_matrix(ĝ)

    # nnz
    nnz = sum(mĝ)
    # shd: structural hamming distance
    shd = sum(mg .!= mĝ)

    if shd > 0
        @warn "Not same"
        @info "true graph:"
        display(g)
        @info "estimated graph:"
        display(ĝ)
    end

    tp = sum(mĝ[mg .== mĝ] .== 1)
    fp = sum(mĝ[mĝ .== 1] .!= mg[mĝ .== 1])
    tt = sum(mg .== 1)
    ff = sum(mg .== 0)

    prec = tp / (tp + fp)
    recall = tp / tt

    # tpr: true positive rate
    # the intersect of edges / total edges in ĝ
    tpr = tp / tt
    # fpr: false positive rate
    fpr = fp / ff
    # fdr: false discovery rate
    fdr = fp / sum(mĝ .== 1)

    # (:nnz=>nnz, :shd=>shd, :prec=>prec, :recall=>recall, :tpr=>tpr, :fpr=>fpr, :fdr=>fdr)
    (:shd=>shd, :prec=>prec, :recall=>recall)
end

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
using ArchGDAL, GDAL; AG = ArchGDAL
using Test

AG.read("data/utmsmall.tif") do ds_small
    @testset "GDAL Error" begin
        @test_throws GDAL.GDALError AG.gdalinfo(ds_small, ["-novalidoption"])
        @test_throws GDAL.GDALError AG.unsafe_gdaltranslate(ds_small, ["-novalidoption"])
        @test_throws GDAL.GDALError AG.unsafe_gdalbuildvrt([ds_small], ["-novalidoption"])
        @test_throws GDAL.GDALError AG.unsafe_gdaldem(ds_small, "hillshade", ["-novalidoption"])
        @test_throws GDAL.GDALError AG.unsafe_gdalnearblack(ds_small, ["-novalidoption"])
        @test_throws GDAL.GDALError AG.unsafe_gdalwarp([ds_small], ["-novalidoption"])
    end

    @testset "GDAL Info" begin
        infostr = AG.gdalinfo(ds_small, ["-checksum"])
        @test occursin("Checksum=50054", infostr)
        info_default = AG.gdalinfo(ds_small)
        @test occursin("Driver: GTiff/GeoTIFF", info_default)
    end

    AG.gdaltranslate(ds_small, # resample to a 5×5 ascii grid
        ["-of","AAIGrid","-r","cubic","-tr","1200","1200"]
    ) do ds_tiny
        @testset "GDAL Translate" begin
            @test AG.read(ds_tiny, 1) == [128  171  127   93   83;
                                          126  164  148  114  101;
                                          161  175  177  164  140;
                                          185  206  205  172  128;
                                          193  205  209  181  122]
        end

        @testset "GDAL Build VRT" begin
            AG.gdalbuildvrt([ds_tiny]) do ds_vrt
                @test AG.read(ds_vrt, 1) == [128  171  127   93   83;
                                             126  164  148  114  101;
                                             161  175  177  164  140;
                                             185  206  205  172  128;
                                             193  205  209  181  122]
            end
        end

        @testset "GDAL DEM Processing" begin
            AG.gdaldem(ds_tiny, "hillshade", ["-of","AAIGrid"]) do ds_dempr
                @test AG.read(ds_dempr, 1) == [ 0    0    0    0  0;
                                                0  183  180  181  0;
                                                0  184  182  181  0;
                                                0  183  181  177  0;
                                                0    0    0    0  0]
            end
        end

        @testset "GDAL Near Black" begin
            AG.gdalnearblack(ds_tiny, ["-of","GTiff","-color","0"]) do ds_nearblack
                @test AG.read(ds_nearblack, 1) == [ 0  0    0  0  0;
                                                    0  0    0  0  0;
                                                    0  0  177  0  0;
                                                    0  0    0  0  0;
                                                    0  0    0  0  0]
            end
        end
    end

    # cannot reproject file on AppVeyor yet
    # GDALError (CE_Failure, code 4):
    #       Unable to open EPSG support file gcs.csv.  Try setting the
    #       GDAL_DATA environment variable to point to the directory
    #       containing EPSG csv files.
    # @testset "GDAL Warp" begin
    #     AG.gdalwarp([ds_small], ["-of","MEM","-t_srs","EPSG:4326"]) do ds_warped
    #         @test AG.width(ds_small) == 100
    #         @test AG.height(ds_small) == 100
    #         @test AG.width(ds_warped) == 109
    #         @test AG.height(ds_warped) == 91
    #     end
    # end
    @testset "GDAL Warp" begin
        AG.gdalwarp([ds_small], ["-of","MEM"]) do ds_warped
            @test AG.width(ds_small) == 100
            @test AG.height(ds_small) == 100
            @test AG.shortname(AG.getdriver(ds_small)) == "GTiff"
            @test AG.width(ds_warped) == 100
            @test AG.height(ds_warped) == 100
            @test AG.shortname(AG.getdriver(ds_warped)) == "MEM"
        end
    end
end

AG.read("data/point.geojson") do ds_point
    @testset "GDAL Grid" begin
        AG.gdalgrid(ds_point, ["-of","MEM","-outsize","3",
            "10","-txe","100","100.3","-tye","0","0.1"]) do ds_grid
            @test AG.getgeotransform(ds_grid) ≈ [100.0,0.1,0.0,0.0,0.0,0.01]
        end
    end

    @testset "GDAL Rasterize" begin
        AG.gdalrasterize(ds_point, ["-of","MEM","-tr","0.05","0.05"]) do ds_rasterize
            @test AG.getgeotransform(ds_rasterize) ≈ [99.975,0.05,0.0,0.1143,0.0,-0.05]
        end
    end

    @testset "GDAL Vector Translate" begin
        AG.gdalvectortranslate([ds_point], ["-f","CSV","-lco",
            "GEOMETRY=AS_XY"], dest = "data/point.csv") do ds_csv
        end
        @test replace(read("data/point.csv", String), "\r" => "") == """
        X,Y,FID,pointname
        100,0,2,point-a
        100.2785,0.0893,3,point-b
        100,0,0,a
        100.2785,0.0893,3,b
        """
        rm("data/point.csv")
    end
end


text:
"""
    ParametricRepresentation{S, P} <: AbstractFilterRepresentation{S}

Abstract type for the representation of the conditional distribution over the hidden state by a parameter vector containing elements of type `P`.
"""
abstract type ParametricRepresentation{S, P} <: AbstractFilterRepresentation{S} end

"""
    parameter_vector(rep::ParametricRepresentation)

Return a parameter vector describing the representation `rep`.
In the language of differential (information) geometry, the function `parameter_vector` provides a chart on a parametric family of distributions.
"""
function parameter_vector(rep::ParametricRepresentation) end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
@testset "examples" begin
    @test TimeDag.output_type(/, Int64, Int64) == Float64
    @test TimeDag.output_type(/, Missing, Int64) == Missing
    @test TimeDag.output_type(/, Union{Missing,Int64}, Int64) == Union{Missing,Float64}
end


text:
using MLJ

X_iris, y_iris = @load_iris;
schema(X_iris)

y_iris[1:4]

levels(y_iris)

DecisionTree = @load DecisionTreeClassifier pkg=DecisionTree # model type
model = DecisionTree(min_samples_split=5)                    # model instance

mach = machine(model, X_iris, y_iris)

train_rows = vcat(1:60, 91:150); # some row indices (observations are rows not columns)
fit!(mach, rows=train_rows)
fitted_params(mach)

mach.model.min_samples_split  = 10
fit!(mach, rows=train_rows) # re-train with new hyper-parameter

predict(mach, rows=71:73)

Xnew = (sepal_length = [5.1, 6.3],
        sepal_width = [3.0, 2.5],
        petal_length = [1.4, 4.9],
        petal_width = [0.3, 1.5])
yhat = predict(mach, Xnew)

pdf.(yhat, "virginica")

yhat[2]

import DataFrames

data = OpenML.load(42178) # data set from OpenML.org
df0 = DataFrames.DataFrame(data)
first(df0, 4)

scitype(["cat", "mouse", "dog"])

schema(df0) |> DataFrames.DataFrame  # converted to DataFrame for better display

fix_blanks(v) = map(v) do x
    if x == " "
        return "0.0"
    else
        return x
    end
end

df0.TotalCharges = fix_blanks(df0.TotalCharges);

coerce!(df0, :TotalCharges => Continuous);

coerce!(df0, Textual => Multiclass);

coerce!(df0, :Churn => OrderedFactor)
levels(df0.Churn) # to check order

schema(df0) |> DataFrames.DataFrame

df, df_test, df_dumped = partition(df0, 0.07, 0.03, # in ratios 7:3:90
                                   stratify=df0.Churn,
                                   rng=123);

y, X = unpack(df, ==(:Churn), !=(:customerID));
schema(X).names

intersect([:Churn, :customerID], schema(X).names)

ytest, Xtest = unpack(df_test, ==(:Churn), !=(:customerID));

Booster = @load EvoTreeClassifier pkg=EvoTrees

booster = Booster()

scitype(y) <: target_scitype(booster)

scitype(X) <: input_scitype(booster)

pipe = ContinuousEncoder() |> booster

pipe.evo_tree_classifier.max_depth

measures("Brier")

mach_pipe = machine(pipe, X, y)

train, validation = partition(1:length(y), 0.7)
fit!(mach_pipe, rows=train)

fp = fitted_params(mach_pipe);
keys(fp)

Set(fp.continuous_encoder.features_to_keep) == Set(schema(X).names)

rpt = report(mach_pipe)
keys(rpt.evo_tree_classifier)

fi = rpt.evo_tree_classifier.feature_importances
feature_importance_table =
    (feature=Symbol.(first.(fi)), importance=last.(fi)) |> DataFrames.DataFrame

ŷ = predict(mach_pipe, rows=validation);
print(
    "Measurements:\n",
    "  brier loss: ", brier_loss(ŷ, y[validation]) |> mean, "\n",
    "  auc:        ", auc(ŷ, y[validation]),                "\n",
    "  accuracy:   ", accuracy(mode.(ŷ), y[validation])
)

confmat(mode.(ŷ), y[validation])

using Plots

roc_curve = roc(ŷ, y[validation])
plt = scatter(roc_curve, legend=false)
plot!(plt, xlab="false positive rate", ylab="true positive rate")
plot!([0, 1], [0, 1], linewidth=2, linestyle=:dash, color=:black)

e_pipe = evaluate(pipe, X, y,
                  resampling=StratifiedCV(nfolds=6, rng=123),
                  measures=[brier_loss, auc, accuracy],
                  repeats=3,
                  acceleration=CPUThreads())

using Measurements

function confidence_intervals(e)
    factor = 2.0 # to get level of 95%
    measure = e.measure
    nfolds = length(e.per_fold[1])
    measurement = [e.measurement[j] ± factor*std(e.per_fold[j])/sqrt(nfolds - 1)
                   for j in eachindex(measure)]
    table = (measure=measure, measurement=measurement)
    return DataFrames.DataFrame(table)
end

confidence_intervals_basic_model = confidence_intervals(e_pipe)

unimportant_features = filter(:importance => <(0.005), feature_importance_table).feature

pipe2 = ContinuousEncoder() |>
    FeatureSelector(features=unimportant_features, ignore=true) |> booster

controls = [
    Step(1),              # to increment iteration parameter (`pipe.nrounds`)
    NumberSinceBest(4),   # main stopping criterion
    TimeLimit(2/3600),    # never train more than 2 sec
    InvalidValue()        # stop if NaN or ±Inf encountered
]

iterated_pipe = IteratedModel(model=pipe2,
                              controls=controls,
                              measure=brier_loss,
                              resampling=Holdout(fraction_train=0.7))

mach_iterated_pipe = machine(iterated_pipe, X, y)
fit!(mach_iterated_pipe);

show(iterated_pipe, 2)

p1 = :(model.evo_tree_classifier.η)
p2 = :(model.evo_tree_classifier.max_depth)

r1 = range(iterated_pipe, p1, lower=-2, upper=-0.5, scale=x->10^x)
r2 = range(iterated_pipe, p2, lower=2, upper=6)

tuning = RandomSearch(rng=123)

tuned_iterated_pipe = TunedModel(model=iterated_pipe,
                                 range=[r1, r2],
                                 tuning=tuning,
                                 measures=[brier_loss, auc, accuracy],
                                 resampling=StratifiedCV(nfolds=6, rng=123),
                                 acceleration=CPUThreads(),
                                 n=40)

mach_tuned_iterated_pipe = machine(tuned_iterated_pipe, X, y)
fit!(mach_tuned_iterated_pipe)

rpt2 = report(mach_tuned_iterated_pipe);
best_booster = rpt2.best_model.model.evo_tree_classifier

print(
    "Optimal hyper-parameters: \n",
    "  max_depth: ", best_booster.max_depth, "\n",
    "  η:         ", best_booster.η
)

e_best = rpt2.best_history_entry
confidence_intervals(e_best)

rpt2.best_report.controls |> collect

plot(mach_tuned_iterated_pipe, size=(600,450))

MLJ.save("tuned_iterated_pipe.jls", mach_tuned_iterated_pipe)

e_tuned_iterated_pipe = evaluate(tuned_iterated_pipe, X, y,
                                 resampling=StratifiedCV(nfolds=6, rng=123),
                                 measures=[brier_loss, auc, accuracy])

confidence_intervals(e_tuned_iterated_pipe)

confidence_intervals_basic_model

mach_restored = machine("tuned_iterated_pipe.jls")

ŷ_tuned = predict(mach_restored, Xtest);
ŷ_tuned[1]

print(
    "Tuned model measurements on test:\n",
    "  brier loss: ", brier_loss(ŷ_tuned, ytest) |> mean, "\n",
    "  auc:        ", auc(ŷ_tuned, ytest),                "\n",
    "  accuracy:   ", accuracy(mode.(ŷ_tuned), ytest)
)

mach_basic = machine(pipe, X, y)
fit!(mach_basic, verbosity=0)

ŷ_basic = predict(mach_basic, Xtest);

print(
    "Basic model measurements on test set:\n",
    "  brier loss: ", brier_loss(ŷ_basic, ytest) |> mean, "\n",
    "  auc:        ", auc(ŷ_basic, ytest),                "\n",
    "  accuracy:   ", accuracy(mode.(ŷ_basic), ytest)
)

# 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:
function widget end

abstract type AbstractWidget{T, S} <: AbstractObservable{S}; end

mutable struct Widget{T, S} <: AbstractWidget{T, S}
    components::OrderedDict{Symbol, Any}
    output::AbstractObservable{S}
    scope
    layout::Function
    function Widget{T}(components::OrderedDict{Symbol, Any};
        output = Observable{Any}(nothing),
        scope = nothing,
        layout = defaultlayout(get_backend())) where {T, S}

        output_obs = to_abstractobservable(output)
        new{T, eltype(output_obs)}(components, to_abstractobservable(output), scope, layout)
    end
end

Widget{T}(components; kwargs...) where {T} = Widget{T}(OrderedDict{Symbol, Any}(Symbol(key) => val for (key, val) in components); kwargs...)

Widget{T}(; components = OrderedDict{Symbol, Any}(), kwargs...) where {T} = Widget{T}(components; kwargs...)

function Widget{T}(w::Widget; kwargs...) where {T}
    dict = Dict{Symbol, Any}(kwargs)
    for field in fieldnames(Widget)
        get!(dict, field, getfield(w, field))
    end
    Widget{T}(; dict...)
end

Widget(w::Widget{T}; kwargs...) where {T} = Widget{T}(w; kwargs...)
Widget(args...; kwargs...) = Widget{:default}(args...; kwargs...)

function widget(f::Function, args...; init = f(map(Observable._val, args)...), kwargs...)
    Widget{:output}(; output = map(f, args...; init = init), kwargs...)
end

widget(f::Function; kwargs...) = w -> widget(f, w; kwargs...)

widgettype(::AbstractWidget{T}) where {T} = T

"""
`scope(w::Widget)`

Return primary scope for widget `w` if it exists, `nothing` otherwise.
"""
scope(w::Widget) = w.scope

"""
`scope!(w::Widget, sc)`

sets up a primary scope `sc` for widget `w`
"""
function scope!(w::Widget, sc)
    w.scope = sc
    w
end

component(x, u) = getindex(x, u)
component(x, args...) = foldl(component, args, init = x)

components(w::Widget) = w.components

observe(u::Widget, args...) = observe(component(u, args...))
observe(u::Widget) = u.output

_getindex(ui::Widget, i::Symbol) = get(components(ui), i, nothing)

function Base.getindex(ui::Widget, i::Symbol)
    val = _getindex(ui, i)
    val === nothing || return val

    for (key, el) in components(ui)
        if el isa Widget
            val = getindex(el, i)
            val === nothing || return val
        end
    end
    return nothing
end

Base.getindex(ui::Widget, i::AbstractString) = getindex(ui, Symbol(i))
Base.setindex!(ui::Widget, val, i::Symbol) = setindex!(components(ui), val, i)
Base.setindex!(ui::Widget, val, i::AbstractString) = setindex!(ui, val, Symbol(i))

"""
`@auto(exprs...)`

Macro to automatize widget creation. Transform `x = rhs` into `x = widget(rhs, label = "x")`.
Return an `OrderedDict` of widgets, which can be used as `components` in a `Widget` object.

## Examples

```julia
julia> wdgs = Widgets.@auto a = 1:100 b = 12

julia> Widget(wdgs, output = map(+, a, b))
```
"""
macro auto(args...)
    esc(auto_helper!(args...))
end

function _block2args(args...)
    length(args) == 1 || return args
    arg = args[1]
    isa(arg, Expr) && (arg.head == :block) || return (arg,)
    return Tuple(Iterators.filter(t -> !isa(t, LineNumberNode), arg.args))
end

function auto_helper!(args...)
    exprs = _block2args(args...)
    res = Any[]
    dict = gensym()
    push!(res, :($dict = Widgets.OrderedDict{Symbol, Any}()))
    for expr in exprs
        @assert expr.head == :(=)
        var = expr.args[1]
        label = name2string(var)
        expr.args[2] = Expr(:call, :(Widgets.widget), expr.args[2], Expr(:kw, :label, label))
        push!(res, expr)
        push!(res, :($dict[$(Expr(:quote, var))] = $var))
    end
    push!(res, dict)
    Expr(:block, res...)
end

# Placeholder for the input function, to define input widgets.
function input end

render(w::Widget) = w.layout(w)

Base.show(io::IO, x::Widget) = show(io, render(x))
Base.show(io::IO, m::MIME"text/html", x::Widget) = show(io, m, render(x))
Base.display(w::Widget) = display(render(w))


text:
@testset "IntegralForms" begin

  @testset "Discriminant" begin
    A = [i for i in -12:13]
    C = [1,1,0,0,1,1,0,0,1,1,0,0,0,1,0,0,1,1,0,0,1,1,0,0,1,1]
    @test is_discriminant.(A) == C
  end

  @testset "Conductor" begin
    A = [5, 8, 12, 13, 17, 20, 21, 24, 28, 29, 32, 33, 37, 40, 41, 44, 45, 48,
         52, 53, 56, 57, 60, 61, 65, 68, 69, 72, 73, 76, 77, 80, 84, 85, 88,
         89, 92, 93, 96, 97, 101, 104, 105, 108, 109, 112, 113, 116, 117, 120,
         124, 125, 128, 129, 132, 133, 136, 137, 140, 141, 145, 148, 149, 152,
         153, 156, 157, 160, 161, 164, 165, 168, 172, 173, 176, 177, 180, 181,
         184, 185, 188, 189, 192, 193, 197, 200, 201, 204, 205, 208, 209, 212,
         213, 216, 217, 220, 221, 224, 228, 229, 232, 233, 236, 237, 240, 241,
         244, 245, 248, 249, 252, 253, 257, 260, 261, 264, 265, 268, 269, 272,
         273, 276, 277, 280, 281, 284, 285, 288, 292, 293, 296, 297, 300, 301,
         304, 305, 308, 309, 312, 313, 316, 317, 320, 321, 325, 328, 329, 332,
         333, 336, 337, 340, 341, 344, 345, 348, 349, 352, 353, 356, 357, 360,
         364, 365, 368, 369, 372, 373, 376, 377, 380, 381, 384, 385, 388, 389,
         392, 393, 396, 397, 401, 404, 405, 408, 409, 412, 413, 416, 417, 420,
         421, 424, 425, 428, 429, 432, 433, 436, 437, 440, 444, 445, 448, 449,
         452, 453, 456, 457, 460, 461, 464, 465, 468, 469, 472, 473, 476, 477,
         480, 481, 485, 488, 489, 492, 493, 496, 497, 500, 501, 504, 505, 508,
         509, 512, 513, 516, 517, 520, 521, 524, 525, 528, 532, 533, 536, 537,
         540, 541, 544, 545, 548, 549, 552, 553, 556, 557, 560, 561, 564, 565,
         568, 569, 572, 573, 577, 580, 581, 584, 585, 588, 589, 592, 593, 596,
         597, 600, 601, 604, 605, 608, 609, 612, 613, 616, 617, 620, 621, 624,
         628, 629, 632, 633, 636, 637, 640, 641, 644, 645, 648, 649, 652, 653,
         656, 657, 660, 661, 664, 665, 668, 669, 672, 673, 677, 680, 681, 684,
         685, 688, 689, 692, 693, 696, 697, 700, 701, 704, 705, 708, 709, 712,
         713, 716, 717, 720, 721, 724, 725, 728, 732, 733, 736, 737, 740, 741,
         744, 745, 748, 749, 752, 753, 756, 757, 760, 761, 764, 765, 768, 769,
         772, 773, 776, 777, 780, 781, 785, 788, 789, 792, 793, 796, 797, 800,
         801, 804, 805, 808, 809, 812, 813, 816, 817, 820, 821, 824, 825, 828,
         829, 832, 833, 836, 837, 840, 844, 845, 848, 849, 852, 853, 856, 857,
         860, 861, 864, 865, 868, 869, 872, 873, 876, 877, 880, 881, 884, 885,
         888, 889, 892, 893, 896, 897, 901, 904, 905, 908, 909, 912, 913, 916,
         917, 920, 921, 924, 925, 928, 929, 932, 933, 936, 937, 940, 941, 944,
         945, 948, 949, 952, 953, 956, 957, 960, 964, 965, 968, 969, 972, 973,
         976, 977, 980, 981, 984, 985, 988, 989, 992, 993, 996, 997, 1000]
    Acond = [1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 3, 2, 2, 1, 1, 1,
             1, 1, 1, 2, 1, 3, 1, 1, 1, 4, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 3,
             1, 2, 1, 2, 3, 1, 1, 5, 4, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1, 3, 1,
             1, 2, 1, 2, 1, 1, 1, 1, 2, 1, 6, 1, 1, 1, 1, 3, 4, 1, 1, 5, 1, 1,
             1, 4, 1, 2, 1, 3, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 2, 1, 2, 7, 1, 1,
             3, 1, 1, 2, 3, 1, 1, 1, 1, 4, 1, 2, 1, 1, 1, 1, 1, 6, 2, 1, 1, 3,
             5, 1, 2, 1, 2, 1, 1, 1, 1, 1, 8, 1, 5, 1, 1, 1, 3, 4, 1, 2, 1, 1,
             1, 1, 1, 2, 1, 2, 1, 3, 1, 1, 2, 3, 2, 1, 1, 1, 1, 1, 4, 1, 2, 1,
             7, 1, 3, 1, 1, 2, 9, 1, 1, 1, 1, 2, 1, 2, 1, 1, 5, 1, 1, 6, 1, 2,
             1, 1, 1, 1, 4, 1, 2, 1, 1, 1, 1, 1, 4, 1, 6, 1, 1, 1, 1, 3, 2, 1,
             1, 1, 1, 1, 1, 2, 1, 10, 1, 3, 1, 1, 1, 8, 3, 2, 1, 1, 1, 1, 5, 4,
             2, 1, 1, 1, 3, 1, 2, 1, 2, 3, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1,
             1, 2, 1, 1, 3, 7, 1, 4, 1, 2, 1, 5, 1, 1, 11, 2, 1, 6, 1, 1, 1, 1,
             3, 2, 2, 1, 1, 1, 1, 7, 4, 1, 2, 1, 9, 1, 1, 1, 4, 3, 2, 1, 1, 1,
             1, 1, 2, 1, 1, 1, 1, 3, 1, 2, 1, 2, 3, 1, 1, 5, 1, 4, 1, 2, 1, 1,
             1, 1, 1, 12, 1, 2, 5, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 2, 1, 6, 1,
             1, 1, 1, 3, 8, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 3, 1, 1, 1, 10, 3, 2,
             1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 5, 3, 1, 8, 7, 2, 3, 1, 1, 13, 4, 1,
             2, 1, 1, 1, 1, 1, 6, 1, 2, 1, 1, 3, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1,
             4, 1, 1, 1, 1, 1, 3, 4, 1, 2, 1, 1, 1, 1, 5, 2, 1, 2, 1, 3, 1, 1,
             1, 2, 3, 2, 1, 1, 1, 1, 1, 4, 2, 1, 11, 1, 9, 1, 4, 1, 14, 3, 1,
             1, 1, 1, 2, 1, 2, 1, 5]
    @test conductor.(A) == Acond
  end

  @testset "FundamentalDiscriminant" begin
    @test is_fundamental_discriminant(12) == true
    @test is_fundamental_discriminant(-12) == false
    D = [1, 5, 8, 12, 13, 17, 21, 24, 28, 29, 33, 37, 40, 41, 44, 53, 56, 57,
         60, 61, 65, 69, 73, 76, 77, 85, 88, 89, 92, 93, 97, 101, 104, 105,
         109, 113, 120, 124, 129, 133, 136, 137]
    @test all(is_fundamental_discriminant, D)
    D = map(x -> -x, [3, 4, 7, 8, 11, 15, 19, 20, 23, 24, 31, 35, 39, 40, 43,
                      47, 51, 52, 55, 56, 59, 67, 68, 71, 79, 83, 84, 87, 88,
                      91, 95, 103])
    @test all(is_fundamental_discriminant, D)

    @test fundamental_discriminant(20) == 5
    @test fundamental_discriminant(-20) == -20
  end

  @testset "DefiniteForms" begin
    @test is_negative_definite(binary_quadratic_form(-2, 3, -2)) == true
    @test is_positive_definite(binary_quadratic_form(2, -3, 2)) == true
    @test is_indefinite(binary_quadratic_form(1, 3, 1)) == true
    @test is_positive_definite(binary_quadratic_form(1, 3, 1)) == false
    @test is_negative_definite(binary_quadratic_form(1, 3, 1)) == false
  end

  @testset "Reduction" begin
    @test is_reduced(binary_quadratic_form(1, 2, 3)) == false
    @test is_reduced(binary_quadratic_form(2, 1, 3)) == true
    @test is_reduced(binary_quadratic_form(1, -1, 1)) == false
    @test is_reduced(binary_quadratic_form(1, 1, 1)) == true
    @test is_reduced(binary_quadratic_form(-1, 2, 2)) == true
    # indefinite
    @test is_reduced(binary_quadratic_form(1, 9, 4)) == false
    @test is_reduced(binary_quadratic_form(1, 5, -1)) == true
    #@test reduce(binary_quadratic_form(195751, 37615, 1807)) == binary_quadratic_form(1, 1, 1)
    #@test reduce(binary_quadratic_form(33, 11, 5)) == binary_quadratic_form(5, -1, 27)
    #@test reduce(binary_quadratic_form(15, 0, 15)) == binary_quadratic_form(15, 0, 15)

    q1 = binary_quadratic_form(15, -7, 0)
    q1red = Hecke.reduction(q1)
    @test is_reduced(q1red)

    q2 = binary_quadratic_form(15, 7, 0)
    q2red = Hecke.reduction(q2)
    @test is_reduced(q2red)
    @test is_equivalent(q1, q2)
  end

#   @testset "Composition" begin
#   @test reduce(compose(binary_quadratic_form(2 ,2, 1), binary_quadratic_form(5, 4, 1))) == binary_quadratic_form(1, 0, 1)
#   @test reduce(compose(binary_quadratic_form(1, 1, 6), binary_quadratic_form(1, 1, 6))) == binary_quadratic_form(1, 1, 6)
#   @test reduce(compose(binary_quadratic_form(2, -1, 3), binary_quadratic_form(2, -1, 3))) == binary_quadratic_form(2, 1, 3)
#   end

  @testset "Cycle" begin

    A1 = [binary_quadratic_form(1, 7, -6), binary_quadratic_form(6, 5, -2), binary_quadratic_form(2, 7, -3), binary_quadratic_form(3, 5, -4),
    binary_quadratic_form(4, 3, -4), binary_quadratic_form(4, 5, -3), binary_quadratic_form(3, 7, -2), binary_quadratic_form(2, 5, -6), binary_quadratic_form(6, 7, -1)]

    @test cycle(binary_quadratic_form(1,7,-6)) == A1

    A2 = [binary_quadratic_form(1, 8, -3), binary_quadratic_form(3, 4, -5), binary_quadratic_form(5, 6, -2), binary_quadratic_form(2, 6, -5),
    binary_quadratic_form(5, 4, -3), binary_quadratic_form(3, 8, -1)]

    @test cycle(binary_quadratic_form(1,8,-3)) == A2

    A3 = [binary_quadratic_form(14, 17, -2), binary_quadratic_form(2, 19, -5), binary_quadratic_form(5, 11, -14)]
    @test cycle(binary_quadratic_form(14, 17, -2)) == A3
  end

  @testset "CanSolve" begin
    @test can_solve(binary_quadratic_form(3,2,2), 28) == (true, (2, 2))
    @test can_solve(binary_quadratic_form(2,1,3), 3) == (true, (0, 1))
    @test can_solve(binary_quadratic_form(2,1,3), 5) == (false, (0, 0))
    @test can_solve(binary_quadratic_form(2,1,3), 6) == (true, (1, 1))
    @test can_solve(binary_quadratic_form(2,-1,3), 3) == (true, (0, 1))
    @test can_solve(binary_quadratic_form(2,-1,3), 5) == (false, (0, 0))
    @test can_solve(binary_quadratic_form(2,-1,3), 6) == (true, (-1, 1))
    @test can_solve(binary_quadratic_form(1,1,6), 3) == (false, (0, 0))
    @test can_solve(binary_quadratic_form(1,1,6), 5) == (false, (0, 0))
    @test can_solve(binary_quadratic_form(1,1,6), 6) == (true, (0, 1))
  end

  @testset "PrimeForm" begin
    f = prime_form(ZZ(17), ZZ(1237))
    @test discriminant(f) == 17 && f[1] == 1237
    f = prime_form(ZZ(12), ZZ(743))
    @test discriminant(f) == 12 && f[1] == 743
  end

  @testset "Equivalent" begin
    f = binary_quadratic_form(fmpz(-1), fmpz(0), fmpz(3))
    g = Hecke.reduction(f)
    @test g[1] == -1 && g[2] == 2 && g[3] == 2
    g, T = Hecke.reduction_with_transformation(f)
    @test g[1] == -1 && g[2] == 2 && g[3] == 2
    @test T == matrix(FlintZZ, 2, 2, [-1, 1, 0, -1])
    @test Hecke._buchmann_vollmer_action(f, T) == g

    # TODO: We do not have equivalence for indefinite types :(
    # f = binary_quadratic_form(4, 4, 15)
    # g = binary_quadratic_form(4, -4, 15)
    # is_equivalent(f, g)

    f = binary_quadratic_form(33, 11, 5)
    #g = reduction(f)
    #@test g == binary_quadratic_form(5, -1, 27)
    #@test is_equivalent(f, g)
    #@test !is_equivalent(f, binary_quadratic_form(3, 4, 5))

    f = binary_quadratic_form(9, 8, -7)
    g = binary_quadratic_form(9, -8, -7)
    @test Hecke.is_locally_equivalent(f, g)
    @test is_equivalent(f, g, proper = false)
    @test !is_equivalent(f, g, proper = true)

    f = binary_quadratic_form(0, 4, 2)
    g = binary_quadratic_form(2, 4, 0)
    @test Hecke.is_locally_equivalent(f, g)
    @test is_equivalent(f, g, proper = false)

    f = binary_quadratic_form(fmpz(3), fmpz(4), fmpz(-2))
    g = binary_quadratic_form(fmpz(-2), fmpz(4), fmpz(3))
    @test is_equivalent(f, g)
    @test is_equivalent(f, g, proper = false)
    @test Hecke.is_locally_equivalent(f, g)

    f = binary_quadratic_form(2, -120, 1785)
    g = binary_quadratic_form(10, -120, 357)
    @test !is_equivalent(f, g)
    @test !Hecke.is_locally_equivalent(f, g)
    @test !is_equivalent(f, g, proper = false)

  end

  @testset "Automormorphism group" begin
    g = binary_quadratic_form(3, 1, -3)
    gens = automorphism_group_generators(g)
    @test all(T -> Hecke._action(g, T) == g, gens)
    @assert all(T -> abs(det(T)) == 1, gens)
    @assert any(T -> det(T) == -1, gens) # g is ambiguous

    g = binary_quadratic_form(0, 3, 0)
    gens = automorphism_group_generators(g)
    @test gens == [ZZ[-1 0;0 -1], ZZ[0 1; 1 0]]

    g = binary_quadratic_form(-2, 1, 0)
    gens = automorphism_group_generators(g)
    @test gens == [ZZ[-1 0; 0 -1], ZZ[-2 1; -3 2]]

    g = binary_quadratic_form(-4, 3, 0)
    gens = automorphism_group_generators(g)
    @test gens == [ZZ[-1 0;0 -1]]

    g = binary_quadratic_form(1, 2, 0)
    gens = automorphism_group_generators(g)
    @test length(gens) == 2

    g = binary_quadratic_form(0, 2, 1)
    gens = automorphism_group_generators(g)
    @test length(gens) == 2

    f = binary_quadratic_form(0, -2, 4)
    gens = automorphism_group_generators(f)
    @test length(gens) == 2

    f = binary_quadratic_form(-7, 15, 2)
    gens = automorphism_group_generators(f)
    @test length(gens) == 3
    @test all(g -> Hecke._action(f, g) == f, gens)

    f = binary_quadratic_form(-10, 8, 7)
    gens = automorphism_group_generators(f)
    @test length(gens) == 3
    @test all(g -> Hecke._action(f, g) == f, gens)

    f = binary_quadratic_form(-9, 9, 5)
    gens = automorphism_group_generators(f)
    @test length(gens) == 3
    @test all(g -> Hecke._action(f, g) == f, gens)

  end

  @testset "Misc" begin
    f = binary_quadratic_form(1, 8, -3)
    g, T = Hecke._tau(f)
    @test g == binary_quadratic_form(-1, 8, 3)
    @test g == Hecke._buchmann_vollmer_action(f, T)

    f = binary_quadratic_form(1, 8, -3)
    g, T = Hecke._rhotau(f)
    @test g == binary_quadratic_form(3, 4, -5)
    @test g == Hecke._buchmann_vollmer_action(f, T)

    f = binary_quadratic_form(1, 8, -3)
    g, T = Hecke._rho(f)
    @test g == binary_quadratic_form(-3, 4, 5)
    @test g == Hecke._buchmann_vollmer_action(f, T)
  end

  @testset "Reducible forms" begin
    f = binary_quadratic_form(1, 13, 0)
    g = binary_quadratic_form(52142, 972209, 4531809)
    fl, T = Hecke._isequivalent_reducible(f, g)
    @test fl
    @test Hecke._action(f, T) == g
    f = binary_quadratic_form(2, 13, 0)
    g = binary_quadratic_form(7, 13, 0)
    @test !is_equivalent(f, g)
    @test is_equivalent(f, g, proper = false)
    @test is_equivalent(2 * f, 2 * f)
    @test !is_equivalent(2 * f, 3 * f)
    h = binary_quadratic_form(-2, 1, 0)
    @test is_reduced(Hecke.reduction(h))

    f = binary_quadratic_form(0, -3, 1)
    @test Hecke.reduction(f) == binary_quadratic_form(1, 3, 0)
  end

  @testset "Representatives" begin
    # (-4)
    # [x^2 + y^2]
    #
    # (-163)
    # [x^2 + x*y + 41*y^2]
    #
    # (-12)
    # [x^2 + 3*y^2, 2*x^2 + 2*x*y + 2*y^2]
    #
    # (-16)
    # [x^2 + 4*y^2, 2*x^2 + 2*y^2]
    #
    # (-63)
    # [x^2 + x*y + 16*y^2, 2*x^2 - x*y + 8*y^2, 2*x^2 + x*y + 8*y^2, 3*x^2 + 3*x*y + 6*y^2, 4*x^2 + x*y + 4*y^2]
    #
    # (-23*9)
    # [x^2 + x*y + 52*y^2,
    # 2*x^2 - x*y + 26*y^2,
    # 2*x^2 + x*y + 26*y^2,
    # 3*x^2 + 3*x*y + 18*y^2,
    # 4*x^2 - x*y + 13*y^2,
    # 4*x^2 + x*y + 13*y^2,
    # 6*x^2 - 3*x*y + 9*y^2,
    # 6*x^2 + 3*x*y + 9*y^2,
    # 8*x^2 + 7*x*y + 8*y^2]
    #
    # (-23*9, primitive = true)
    # [x^2 + x*y + 52*y^2,
    # 2*x^2 - x*y + 26*y^2,
    # 2*x^2 + x*y + 26*y^2,
    # 4*x^2 - x*y + 13*y^2,
    # 4*x^2 + x*y + 13*y^2,
    # 8*x^2 + 7*x*y + 8*y^2]
    Zx, (x, y) = ZZ["x", "y"]
    poly_to_form(f) = binary_quadratic_form(coeff(f, x^2), coeff(f, x*y), coeff(f, y^2))
    d = fmpz(73)
    @test length(Hecke.binary_quadratic_form_representatives(d)) ==
      length([4*x^2 + 3*x*y - 4*y^2])
    d = fmpz(76)
    @test length(Hecke.binary_quadratic_form_representatives(d, primitive = true)) ==
      length([-3*x^2 + 4*x*y + 5*y^2,
            3*x^2 + 4*x*y - 5*y^2])
    d = fmpz(136)
    @test length(Hecke.binary_quadratic_form_representatives(d)) ==
      length([-5*x^2 + 4*x*y + 6*y^2,
            -2*x^2 + 8*x*y + 9*y^2,
            2*x^2 + 8*x*y - 9*y^2,
            5*x^2 + 4*x*y - 6*y^2])
    @test length(Hecke.binary_quadratic_form_representatives(d, proper = false)) ==
      length([-2*x^2 + 8*x*y + 9*y^2, 2*x^2 + 8*x*y - 9*y^2, 5*x^2 + 4*x*y - 6*y^2])

    d = fmpz(148)
    @test length(Hecke.binary_quadratic_form_representatives(d, proper = false, primitive = false)) ==
      length([x^2 + 12*x*y - y^2, 4*x^2 + 6*x*y - 7*y^2, 6*x^2 + 2*x*y - 6*y^2])
    @test length(Hecke.binary_quadratic_form_representatives(d, proper = false, primitive = true)) ==
      length([x^2 + 12*x*y - y^2, 4*x^2 + 6*x*y - 7*y^2])
    @test length(Hecke.binary_quadratic_form_representatives(d, proper = true, primitive = true)) ==
      length([-7*x^2 + 6*x*y + 4*y^2, x^2 + 12*x*y - y^2, 4*x^2 + 6*x*y - 7*y^2])
    @test length(Hecke.binary_quadratic_form_representatives(d, proper = true, primitive = false)) ==
      length([-7*x^2 + 6*x*y + 4*y^2,
            x^2 + 12*x*y - y^2,
            4*x^2 + 6*x*y - 7*y^2,
            6*x^2 + 2*x*y - 6*y^2])
    # (10^2, proper = false, primitive = false)
    #        [-4*x^2 + 10*x*y,
    #         -3*x^2 + 10*x*y,
    #         -2*x^2 + 10*x*y,
    #         -x^2 + 10*x*y,
    #         10*x*y,
    #         x^2 + 10*x*y,
    #         2*x^2 + 10*x*y,
    #         5*x^2 + 10*x*y]
    #
    # (10^2, proper = false, primitive = true)
    #        [-3*x^2 + 10*x*y, -x^2 + 10*x*y, x^2 + 10*x*y]
    # (10^2, proper = true, primitive = true)
    #        [-3*x^2 + 10*x*y, -x^2 + 10*x*y, x^2 + 10*x*y, 3*x^2 + 10*x*y]
    # (10^2, proper = true, primitive = false)
    #        [-4*x^2 + 10*x*y,
    #         -3*x^2 + 10*x*y,
    #         -2*x^2 + 10*x*y,
    #         -x^2 + 10*x*y,
    #         10*x*y,
    #         x^2 + 10*x*y,
    #         2*x^2 + 10*x*y,
    #         3*x^2 + 10*x*y,
    #         4*x^2 + 10*x*y,
    #         5*x^2 + 10*x*y]
  end

  @testset "Printing" begin
    q = binary_quadratic_form(1, 2, 3)
    K, a = cyclotomic_field(4)
    q1 = binary_quadratic_form(K, K(3), a, K(0))
    q2 = binary_quadratic_form(0, 1, 0)
    @test sprint(show, q) isa String
    @test sprint(show,"text/plain" , q) isa String
    @test sprint(show, q1) isa String
    @test sprint(show,"text/plain" , q1) isa String
    @test sprint(show, q2) isa String
    @test sprint(show,"text/plain" , q2) isa String
  end
end


text:
abstract type Univariate end

struct UnivariateCalibration <: Univariate
    Slope::Float64
    Offset::Float64
    Slope_Uncertainty::Float64
    Offset_Uncertainty::Float64
    Rsq::Float64
    Fstatistic::Float64
    DOF::Int
end

"""
    UnivariateCalibration( Y, X )
    Performs a univariate least squares regression with a bias/offset term(`y = mx + b`).
Returns a Univariate Calibration object containing uncertainty in parameters, and other statistics.
"""
function UnivariateCalibration( Y, X )
    len = length( X )
    @assert(len == length(Y))
    μ_x, μ_y    = sum( X ) / len                , sum( Y ) / len
    SSxx, SSyy  = sum( (X .^ 2) .- (μ_x ^ 2) )  , sum( (Y .^ 2) .- (μ_y ^ 2) )
    SSxy = sum( (X .* Y) .- (μ_x * μ_y) )

    m = SSxy / SSxx
    b = μ_y - ( m * μ_x )

    Y_hat = ( m .* X ) .+ b
    SSe = sum( (Y .- Y_hat) .^ 2 )
    SSr = SSyy - SSe

    Syx = sqrt(SSe / ( len - 2 ))
    Sm = Syx / sqrt( sum( (X .- μ_x) .^ 2 ))
    Sb = Syx * sqrt( sum( X .^ 2 ) / ( len * SSxx ) )

    Rsq = SSr / SSyy
    Fstat = ( SSxx - SSe ) / Syx

    return UnivariateCalibration( m, b, Sm, Sb, Rsq, Fstat, len - 2 )
end

(M::UnivariateCalibration)(X) = (X*M.Slope) + M.Offset

struct StandardAddition
    Unknown::Float64
    Slope::Float64
    Offset::Float64
    Unknown_Uncertainty::Float64
    Slope_Uncertainty::Float64
    Offset_Uncertainty::Float64
    Rsq::Float64
    DOF::Int
end

"""
    StandardAddition( Signal, Spike )

    Performs a univariate standard addition calibration. Returns a StandardAddition
object.

"""
function StandardAddition( Signal, Spike )
    len = length( Signal )
    DOF = len - 2
    @assert(len == length(Spike))
    mu_x, mu_y = sum( Spike ) / len             , sum( Signal ) / len
    SSxx, SSyy = sum( ( Spike .- mu_x ) .^ 2 )  , sum( ( Signal .- mu_y ) .^ 2 )
    SSxy = sum( ( Spike .- mu_x ) .* ( Signal .- mu_y ) )

    m = SSxy / SSxx
    #y = mx + b -> y - mx =  b
    b = mu_y - ( m * mu_x )
    Y_hat = ( m * Spike ) .+ b
    #Find X - intercept...
    x0 = -b / m #y = 0 = mx + b => -b / m = x

    Sy = sqrt( sum( (Signal .- Y_hat) .^ 2 ) / DOF )
    Sx = (Sy / abs(m)) * sqrt( (1/len) * ( mu_y^2 / ( m^2 * sum( ( Spike .- mu_x ) .^ 2 )  ) ) )

    SSe = sum( (Signal .- Y_hat) .^ 2 )
    SSr = SSyy - SSe

    Syx = SSe / DOF
    Sm = Syx / SSxx
    Sb = Syx * ( ( 1 / len ) + ( ( mu_x ^ 2 ) / SSxx ) )

    Rsq = SSr / SSyy
    return StandardAddition( x0, m, b, Sy, Sm, Sb, Rsq, DOF )
end

"""
    Confidence_Offset( UC::Univariate; Significance = 0.05 )
Returns a tuple of (-, +) of the estimated offset/bias's confidence interval from a `UC` Univariate type object.
"""
function Confidence_Offset( UC::Univariate; Significance = 0.05 )
    if UC.DOF < 31
        extent = UC.Offset_Uncertainty * -quantile( TDist( UC.DOF ), Significance / 2 )
    else
        extent = UC.Offset_Uncertainty * -quantile( Normal(), Significance / 2 )
    end
    return ( UC.Offset - extent, UC.Offset + extent )
end

"""
    Confidence_Slope( UC::Univariate; Significance = 0.05 )
Returns a tuple of (-, +) of the estimated slope's confidence interval from a `UC` Univariate type object.
"""
function Confidence_Slope( UC::Univariate; Significance = 0.05 )
    if UC.DOF < 31
        extent = UC.Slope_Uncertainty * -quantile( TDist( UC.DOF ), Significance / 2 )
    else
        extent = UC.Slope_Uncertainty * -quantile( Normal(), Significance / 2 )
    end
    return ( UC.Slope - extent, UC.Slope + extent )
end

max_issues_repo_name: caseykneale/ChemometricsTools

max_forks_repo_name: caseykneale/ChemometricsTools


text:
# Ref: https://github.com/JuliaSmoothOptimizers/Percival.jl/blob/master/src/AugLagModel.jl
mutable struct PenaltyEvaluator{T} <: AbstractNLPEvaluator
    inner::AbstractNLPEvaluator
    cons::AbstractVector{T}
    infeasibility::AbstractVector{T}
    penalties::AbstractVector{T}
    # Scaling
    scaler::AbstractScaler
end

function PenaltyEvaluator(nlp::AbstractNLPEvaluator, u0;
                          scale=false,
                          penalties=Float64[],
                          c₀=0.1)

    if n_constraints(nlp) == 0
        @warn("Original model has no inequality constraint")
    end
    if length(penalties) != length(nlp.constraints)
        penalties = Float64[c₀=c₀ for cons in nlp.constraints]
    end

    cons = similar(nlp.g_min)
    cx = similar(nlp.g_min)

    if scale
        scaler = MaxScaler(nlp, u0)
    else
        scaler = MaxScaler(nlp.g_min, nlp.g_max)
    end

    return PenaltyEvaluator(nlp, cons, cx, penalties, scaler)
end

function update!(pen::PenaltyEvaluator, u)
    conv = update!(pen.inner, u)
    # Update constraints
    constraint!(pen.inner, pen.cons, u)
    # Rescale
    pen.cons .*= pen.scaler.scale_cons
    # Update infeasibility error
    g♭ = pen.scaler.g_min
    g♯ = pen.scaler.g_max
    pen.infeasibility .= max.(0, pen.cons .- g♯) .+ min.(0, pen.cons .- g♭)
    return conv
end

function update_penalty!(pen::PenaltyEvaluator; η=10.0)
    fr_ = 0
    pen.penalties = min.(η * pen.penalties, 10e12)
end

function objective(pen::PenaltyEvaluator, u)
    base_nlp = pen.inner
    # Internal objective
    obj = pen.scaler.scale_obj * objective(base_nlp, u)
    # Add penalty terms
    fr_ = 0
    for (πp, cons) in zip(pen.penalties, base_nlp.constraints)
        n = size_constraint(base_nlp.model, cons)
        mask = fr_+1:fr_+n
        cx = pen.infeasibility[mask]
        obj += 0.5 * πp * dot(cx, cx)
        fr_ += n
    end
    return obj
end

function gradient!(pen::PenaltyEvaluator, grad, u)
    base_nlp = pen.inner
    model = base_nlp.model
    scaler = pen.scaler
    λ = base_nlp.λ
    # Import buffer (has been updated previously in update!)
    buffer = base_nlp.buffer
    # Import AutoDiff objects
    ∇gᵤ = jacobian(model, base_nlp.autodiff.Jgᵤ, buffer)
    ∂obj = base_nlp.autodiff.∇f
    jvx = ∂obj.jvₓ ; fill!(jvx, 0)
    jvu = ∂obj.jvᵤ ; fill!(jvu, 0)

    # compute gradient of objective
    ∂cost(model, ∂obj, buffer)
    jvx .+= scaler.scale_obj .* ∂obj.∇fₓ
    jvu .+= scaler.scale_obj .* ∂obj.∇fᵤ

    # compute gradient of penalties
    fr_ = 0
    for (πp, cons) in zip(pen.penalties, base_nlp.constraints)
        n = size_constraint(base_nlp.model, cons)
        mask = fr_+1:fr_+n
        cx = @view pen.infeasibility[mask]
        v = cx .* πp .* scaler.scale_cons[mask]
        jtprod(model, cons, ∂obj, buffer, v)
        jvx .+= ∂obj.∇fₓ
        jvu .+= ∂obj.∇fᵤ
        fr_ += n
    end
    # evaluate reduced gradient
    LinearSolvers.ldiv!(base_nlp.linear_solver, λ, base_nlp.∇gᵗ, jvx)
    grad .= jvu
    mul!(grad, transpose(∇gᵤ), λ, -1.0, 1.0)
end

primal_infeasibility!(pen::PenaltyEvaluator, cons, u) = primal_infeasibility!(pen.inner, cons, u)
primal_infeasibility(pen::PenaltyEvaluator, u) = primal_infeasibility(pen.inner, u)

function reset!(pen::PenaltyEvaluator)
    reset!(pen.inner)
    fill!(pen.cons, 0)
    fill!(pen.infeasibility, 0)
    fill!(pen.penalties, 0)
end


text:
module CURAND

import CUDAdrv: CUDAdrv, CuContext, CuPtr
import CUDAapi

using ..CuArrays
using ..CuArrays: libcurand, active_context

using GPUArrays

using Random

export rand_logn!, rand_poisson!

include("libcurand_types.jl")
include("error.jl")

const _generators = Dict{CuContext,RNG}()
const _generator = Ref{Union{Nothing,RNG}}(nothing)

function generator()
    if _generator[] == nothing
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _generator[] = get!(_generators, active_context[]) do
            context = active_context[]
            generator = create_generator()
            # FIXME: crashes
            #atexit(()->CUDAdrv.isvalid(context) && destroy_generator(generator))
            generator
        end
    end

    return _generator[]::RNG
end

include("libcurand.jl")
include("highlevel.jl")

version() = VersionNumber(curandGetProperty(CUDAapi.MAJOR_VERSION),
                          curandGetProperty(CUDAapi.MINOR_VERSION),
                          curandGetProperty(CUDAapi.PATCH_LEVEL))

end

const rand = CURAND.rand
const randn = CURAND.randn
const rand_logn = CURAND.rand_logn
const rand_poisson = CURAND.rand_poisson

@deprecate curand CuArrays.rand
@deprecate curandn CuArrays.randn


text:
import LinearAlgebra: Symmetric
@explicit_intercepts Symmetric Tuple{∇Array}
∇(::Type{Symmetric}, ::Type{Arg{1}}, p, Y::∇Array, Ȳ::∇Array, X::∇Array) =
    UpperTriangular(Ȳ) + LowerTriangular(Ȳ)' - Diagonal(Ȳ)


text:
module TestFunctionCompilation

using Test
using DataFramesMeta

@testset "function_compilation" begin
    @eval begin
        df = DataFrame(a = [1], b = [2])

        testfun(x, y) = x .* y
        testdotfun(x, y) = x * y
        testnt(x) = (c = x,)
    end

    # Lazy way of making sure all functions are pre-compiled.
    # @eval prevents julia from caching the intermediate anonymous functions.
    for _ in 1:2
        @eval begin
            @test @select(df, :c =  :a + :b) == DataFrame(c = [3])

            fasttime = @timed @select(df, :c =  :a + :b)
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> a + b) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :c =  begin :a + :b end) == DataFrame(c = [3])

            fasttime = @timed @select(df, :c =  begin :a + :b end)
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> a + b) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:c) = :a + :b) == DataFrame(c = [3])

            fasttime = @timed @select(df, cols(:c) = :a + :b)
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> a + b) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :c =  :a .+ :b) == DataFrame(c = [3])

            fasttime = @timed @select(df, :c =  :a .+ :b)
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> a .+ b) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:c) = :a .+ :b) == DataFrame(c = [3])

            fasttime = @timed @select(df, cols(:c) = :a .+ :b)
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> a .+ b) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :c =  cols(:a) + cols(:b)) == DataFrame(c = [3])

            fasttime = @timed @select(df, :c =  cols(:a) + cols(:b))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> a + b) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:c) = cols(:a) + cols(:b)) == DataFrame(c = [3])

            fasttime = @timed @select(df, cols(:c) = cols(:a) + cols(:b))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> a + b) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :c =  cols(:a) .+ cols(:b)) == DataFrame(c = [3])

            fasttime = @timed @select(df, :c =  cols(:a) .+ cols(:b))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> a .+ b) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:c) = cols(:a) .+ cols(:b)) == DataFrame(c = [3])

            fasttime = @timed @select(df, cols(:c) = cols(:a) .+ cols(:b))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> a .+ b) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :c =  :a) == DataFrame(c = [1])

            fasttime = @timed @select(df, :c =  :a)
            slowtime = @timed select(df, [:a] => (a -> identity(a)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:c) = :a) == DataFrame(c = [1])

            fasttime = @timed @select(df, cols(:c) = :a)
            slowtime = @timed select(df, [:a] => (a -> identity(a)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :c =  cols(:a)) == DataFrame(c = [1])

            fasttime = @timed @select(df, :c =  cols(:a))
            slowtime = @timed select(df, [:a] => (a -> identity(a)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:c) = cols(:a)) == DataFrame(c = [1])

            fasttime = @timed @select(df, cols(:c) = cols(:a))
            slowtime = @timed select(df, [:a] => (a -> identity(a)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :a) == df[:, [:a]]

            fasttime = @timed @select(df, :a)
            slowtime = @timed select(df, [:a] => (a -> identity(a)) => :a)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:a)) == df[:, [:a]]

            fasttime = @timed @select(df, cols(:a))
            slowtime = @timed select(df, [:a] => (a -> identity(a)) => :a)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :c =  testfun(:a, :b)) == DataFrame(c = [2])

            fasttime = @timed @select(df, :c =  testfun(:a, :b))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> testfun(a, b)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:c) = testfun(:a, :b)) == DataFrame(c = [2])

            fasttime = @timed @select(df, cols(:c) = testfun(:a, :b))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> testfun(a, b)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :c =  testfun(cols("a"), cols("b"))) == DataFrame(c = [2])

            fasttime = @timed @select(df, :c =  testfun(cols("a"), cols("b")))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> testfun(a, b)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:c) = testfun(cols("a"), cols("b"))) == DataFrame(c = [2])

            fasttime = @timed @select(df, cols(:c) = testfun(cols("a"), cols("b")))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> testfun(a, b)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :c =  testdotfun.(:a, :b)) == DataFrame(c = [2])

            fasttime = @timed @select(df, :c =  testdotfun.(:a, :b))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> testdotfun.(a, b)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:c) = testdotfun.(:a, :b)) == DataFrame(c = [2])

            fasttime = @timed @select(df, cols(:c) = testdotfun.(:a, :b))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> testdotfun.(a, b)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :c =  testdotfun.(cols("a"), cols("b"))) == DataFrame(c = [2])

            fasttime = @timed @select(df, :c =  testdotfun.(cols("a"), cols("b")))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> testdotfun.(a, b)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, cols(:c) = testdotfun.(cols("a"), cols("b"))) == DataFrame(c = [2])

            fasttime = @timed @select(df, cols(:c) = testdotfun.(cols("a"), cols("b")))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> testdotfun.(a, b)) => :c)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            gd = groupby(df, :a)

            @test @combine(gd, cols(AsTable) = testnt(:b)) == DataFrame(a = [1], c =  [2])

            fasttime = @timed @combine(gd, cols(AsTable) = testnt(:b))
            slowtime = @timed combine(gd, :b => (b -> testnt(b)) => AsTable)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @with df (:a + :b) == [3]

            fasttime = @timed @with df (:a + :b)
            slowtime = @timed @with df ((a, b) -> a + b)(df.a, df.b)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @with df (:a .* :b) == [2]

            @with df (:a .* :b)
            fasttime = @timed @with df (:a .* :b)
            slowtime = @timed @with df ((a, b) -> a .* b)(df.a, df.b)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @with df testfun(:a, :b) == [2]

            fasttime = @timed @with df testfun(:a, :b)
            slowtime = @timed @with df ((a, b) -> testfun(a, b))(df.a, df.b)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @with df testdotfun.(:a, :b) == [2]

            fasttime = @timed @with df testdotfun.(:a, :b)
            slowtime = @timed @with df ((a, b) -> testdotfun.(a, b))(df.a, df.b)
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")
        end
    end
end

@testset "composed compilation" begin
    @eval begin
        df = DataFrame(a = [1], b = [2])

        f(x) = identity(x)
        g(x, y) = x + y

        df_wide = DataFrame(rand(10, 1000), :auto)
    end

    for _ in 1:2
        @eval begin
            @test @select(df, :y = (f ∘ g)(:a, :b)).y == [3]

            fasttime = @timed @select(df, :y = (f ∘ g)(:a, :b))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> (f ∘ g)(a, b)) => :y )
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :y = (f ∘ g).(:a, :b)).y == [3]

            fasttime = @timed @select(df, :y = (f ∘ g).(:a, :b))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> (f ∘ g).(a, b)) => :y )
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            fasttime = @timed @rselect df_wide :y = (sum ∘ skipmissing)(AsTable(:))
            slowtime = @timed select(df_wide, AsTable(:) => ByRow(t -> (sum ∘ skipmissing)(t)) => :y)

            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            @test @select(df, :y = f(g(:a, :b))).y == [3]

            fasttime = @timed @select(df, :y = f(g(:a, :b)))
            slowtime = @timed select(df, [:a, :b] => ((a, b) -> f(g(a, b))) => :y )
            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")

            fasttime = @timed @rselect df_wide :y = sum(skipmissing(AsTable(:)))
            slowtime = @timed select(df_wide, AsTable(:) => ByRow(t -> sum(skipmissing(t))) => :y)

            (slowtime[2] > fasttime[2]) || @warn("Slow compilation")
        end
    end

    # Tests for correctness
    t = @select df :y = sum(skipmissing(AsTable([:a, :b])))
    @test t == DataFrame(y = [3])

    t = @select df :y = f(g(:a, :b))
    @test t == DataFrame(y = [3])

    a_str = "a"
    b_str = "b"

    t = @select df :y = f(g($a_str, $(b_str)))
    @test t == DataFrame(y = [3])

    t = @select df :y = f(g(:a, $b_str))
    @test t == DataFrame(y = [3])

    t = @select df :y = sum(skipmissing(AsTable(Cols(:))))
    @test t == DataFrame(y = [3])

    t = @select df :y = sum(skipmissing(AsTable(Cols(:))))
    @test t == DataFrame(y = [3])

    t = DataFramesMeta.@col :y = (identity ∘ first)(:x)
    @test t == ([:x] => (identity ∘ first => :y))

    t = DataFramesMeta.@col :y = identity(first(:x))
    @test t == ([:x] => (identity ∘ first => :y))

    t = DataFramesMeta.@col :y = identity(first(AsTable(Cols(:))))

    @test t == (AsTable(Cols(:)) => ((identity ∘ first) => :y))
end

end # module


text:
using JuMP
using Xpress

function reopt_model(MAXTIME)

    REopt = direct_model(Xpress.Optimizer(MAXTIME=-MAXTIME))
	return REopt

end


text: abstract type Storage <: Injection end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: lex_runtime_v2
using AWS.Compat
using AWS.UUIDs

"""
    delete_session(bot_alias_id, bot_id, locale_id, session_id)
    delete_session(bot_alias_id, bot_id, locale_id, session_id, params::Dict{String,<:Any})

Removes session information for a specified bot, alias, and user ID.  You can use this
operation to restart a conversation with a bot. When you remove a session, the entire
history of the session is removed so that you can start again. You don't need to delete a
session. Sessions have a time limit and will expire. Set the session time limit when you
create the bot. The default is 5 minutes, but you can specify anything between 1 minute and
24 hours. If you specify a bot or alias ID that doesn't exist, you receive a
BadRequestException.  If the locale doesn't exist in the bot, or if the locale hasn't been
enables for the alias, you receive a BadRequestException.

# Arguments
- `bot_alias_id`: The alias identifier in use for the bot that contains the session data.
- `bot_id`: The identifier of the bot that contains the session data.
- `locale_id`: The locale where the session is in use.
- `session_id`: The identifier of the session to delete.

"""
delete_session(botAliasId, botId, localeId, sessionId; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_v2("DELETE", "/bots/$(botId)/botAliases/$(botAliasId)/botLocales/$(localeId)/sessions/$(sessionId)"; aws_config=aws_config)
delete_session(botAliasId, botId, localeId, sessionId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_v2("DELETE", "/bots/$(botId)/botAliases/$(botAliasId)/botLocales/$(localeId)/sessions/$(sessionId)", params; aws_config=aws_config)

"""
    get_session(bot_alias_id, bot_id, locale_id, session_id)
    get_session(bot_alias_id, bot_id, locale_id, session_id, params::Dict{String,<:Any})

Returns session information for a specified bot, alias, and user. For example, you can use
this operation to retrieve session information for a user that has left a long-running
session in use. If the bot, alias, or session identifier doesn't exist, Amazon Lex returns
a BadRequestException. If the locale doesn't exist or is not enabled for the alias, you
receive a BadRequestException.

# Arguments
- `bot_alias_id`: The alias identifier in use for the bot that contains the session data.
- `bot_id`: The identifier of the bot that contains the session data.
- `locale_id`: The locale where the session is in use.
- `session_id`: The identifier of the session to return.

"""
get_session(botAliasId, botId, localeId, sessionId; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_v2("GET", "/bots/$(botId)/botAliases/$(botAliasId)/botLocales/$(localeId)/sessions/$(sessionId)"; aws_config=aws_config)
get_session(botAliasId, botId, localeId, sessionId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_v2("GET", "/bots/$(botId)/botAliases/$(botAliasId)/botLocales/$(localeId)/sessions/$(sessionId)", params; aws_config=aws_config)

"""
    put_session(bot_alias_id, bot_id, locale_id, session_id, session_state)
    put_session(bot_alias_id, bot_id, locale_id, session_id, session_state, params::Dict{String,<:Any})

Creates a new session or modifies an existing session with an Amazon Lex bot. Use this
operation to enable your application to set the state of the bot.

# Arguments
- `bot_alias_id`: The alias identifier of the bot that receives the session data.
- `bot_id`: The identifier of the bot that receives the session data.
- `locale_id`: The locale where the session is in use.
- `session_id`: The identifier of the session that receives the session data.
- `session_state`: Sets the state of the session with the user. You can use this to set the
  current intent, attributes, context, and dialog action. Use the dialog action to determine
  the next step that Amazon Lex should use in the conversation with the user.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ResponseContentType"`: The message that Amazon Lex returns in the response can be
  either text or speech depending on the value of this parameter.    If the value is
  text/plain; charset=utf-8, Amazon Lex returns text in the response.
- `"messages"`: A list of messages to send to the user. Messages are sent in the order that
  they are defined in the list.
- `"requestAttributes"`: Request-specific information passed between Amazon Lex and the
  client application. The namespace x-amz-lex: is reserved for special attributes. Don't
  create any request attributes with the prefix x-amz-lex:.
"""
put_session(botAliasId, botId, localeId, sessionId, sessionState; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_v2("POST", "/bots/$(botId)/botAliases/$(botAliasId)/botLocales/$(localeId)/sessions/$(sessionId)", Dict{String, Any}("sessionState"=>sessionState); aws_config=aws_config)
put_session(botAliasId, botId, localeId, sessionId, sessionState, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_v2("POST", "/bots/$(botId)/botAliases/$(botAliasId)/botLocales/$(localeId)/sessions/$(sessionId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("sessionState"=>sessionState), params)); aws_config=aws_config)

"""
    recognize_text(bot_alias_id, bot_id, locale_id, session_id, text)
    recognize_text(bot_alias_id, bot_id, locale_id, session_id, text, params::Dict{String,<:Any})

Sends user input to Amazon Lex. Client applications use this API to send requests to Amazon
Lex at runtime. Amazon Lex then interprets the user input using the machine learning model
that it build for the bot. In response, Amazon Lex returns the next message to convey to
the user and an optional response card to display.

# Arguments
- `bot_alias_id`: The alias identifier in use for the bot that processes the request.
- `bot_id`: The identifier of the bot that processes the request.
- `locale_id`: The locale where the session is in use.
- `session_id`: The identifier of the user session that is having the conversation.
- `text`: The text that the user entered. Amazon Lex interprets this text.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"requestAttributes"`: Request-specific information passed between the client application
  and Amazon Lex  The namespace x-amz-lex: is reserved for special attributes. Don't create
  any request attributes with the prefix x-amz-lex:.
- `"sessionState"`: The current state of the dialog between the user and the bot.
"""
recognize_text(botAliasId, botId, localeId, sessionId, text; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_v2("POST", "/bots/$(botId)/botAliases/$(botAliasId)/botLocales/$(localeId)/sessions/$(sessionId)/text", Dict{String, Any}("text"=>text); aws_config=aws_config)
recognize_text(botAliasId, botId, localeId, sessionId, text, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_v2("POST", "/bots/$(botId)/botAliases/$(botAliasId)/botLocales/$(localeId)/sessions/$(sessionId)/text", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("text"=>text), params)); aws_config=aws_config)

"""
    recognize_utterance(content-_type, bot_alias_id, bot_id, locale_id, session_id)
    recognize_utterance(content-_type, bot_alias_id, bot_id, locale_id, session_id, params::Dict{String,<:Any})

Sends user input to Amazon Lex. You can send text or speech. Clients use this API to send
text and audio requests to Amazon Lex at runtime. Amazon Lex interprets the user input
using the machine learning model built for the bot.

# Arguments
- `content-_type`: Indicates the format for audio input or that the content is text. The
  header must start with one of the following prefixes:   PCM format, audio data must be in
  little-endian byte order.   audio/l16; rate=16000; channels=1   audio/x-l16;
  sample-rate=16000; channel-count=1   audio/lpcm; sample-rate=8000; sample-size-bits=16;
  channel-count=1; is-big-endian=false     Opus format
  audio/x-cbr-opus-with-preamble;preamble-size=0;bit-rate=256000;frame-size-milliseconds=4
   Text format   text/plain; charset=utf-8
- `bot_alias_id`: The alias identifier in use for the bot that should receive the request.
- `bot_id`: The identifier of the bot that should receive the request.
- `locale_id`: The locale where the session is in use.
- `session_id`: The identifier of the session in use.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Response-Content-Type"`: The message that Amazon Lex returns in the response can be
  either text or speech based on the responseContentType value.   If the value is
  text/plain;charset=utf-8, Amazon Lex returns text in the response.   If the value begins
  with audio/, Amazon Lex returns speech in the response. Amazon Lex uses Amazon Polly to
  generate the speech using the configuration that you specified in the requestContentType
  parameter. For example, if you specify audio/mpeg as the value, Amazon Lex returns speech
  in the MPEG format.   If the value is audio/pcm, the speech returned is audio/pcm at 16 KHz
  in 16-bit, little-endian format.   The following are the accepted values:   audio/mpeg
  audio/ogg   audio/pcm (16 KHz)   audio/* (defaults to mpeg)   text/plain; charset=utf-8
- `"inputStream"`: User input in PCM or Opus audio format or text format as described in
  the requestContentType parameter.
- `"x-amz-lex-request-attributes"`: Request-specific information passed between the client
  application and Amazon Lex  The namespace x-amz-lex: is reserved for special attributes.
  Don't create any request attributes for prefix x-amz-lex:.
- `"x-amz-lex-session-state"`: Sets the state of the session with the user. You can use
  this to set the current intent, attributes, context, and dialog action. Use the dialog
  action to determine the next step that Amazon Lex should use in the conversation with the
  user.
"""
recognize_utterance(Content_Type, botAliasId, botId, localeId, sessionId; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_v2("POST", "/bots/$(botId)/botAliases/$(botAliasId)/botLocales/$(localeId)/sessions/$(sessionId)/utterance", Dict{String, Any}("headers"=>Dict{String, Any}("Content-Type"=>Content_Type)); aws_config=aws_config)
recognize_utterance(Content_Type, botAliasId, botId, localeId, sessionId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = lex_runtime_v2("POST", "/bots/$(botId)/botAliases/$(botAliasId)/botLocales/$(localeId)/sessions/$(sessionId)/utterance", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("headers"=>Dict{String, Any}("Content-Type"=>Content_Type)), params)); aws_config=aws_config)


text:
function jadd_vertex!(g::JGraph; starting_frame=1)

    add_vertex!(g.data.graph)
    positions = _adapt_to_lims(g.data.layout(g.data.graph), g.data.width, g.data.height)
    push!(
        g.jnodes,
        Object(starting_frame:g.data.frames[end], g.data.node_templates, positions[end])
    )
    jgraph_morph(g, g.data.layout, (g.data.width, g.data.height), frames=starting_frame:starting_frame)
end


text:
using RipQP, LinearOperators, LLSModels, QPSReader, QuadraticModels, Quadmath
using DelimitedFiles, LinearAlgebra, MatrixMarket, SparseArrays, Test

qps1 = readqps("QAFIRO.SIF") #lower bounds
qps2 = readqps("HS21.SIF") # low/upp bounds
qps3 = readqps("HS52.SIF") # free bounds
qps4 = readqps("AFIRO.SIF") # LP
qps5 = readqps("HS35MOD.SIF") # fixed variables

include("solvers/test_augmented.jl")
include("solvers/test_newton.jl")
include("solvers/test_normal.jl")

include("test_multi.jl")
include("test_alternative_methods.jl")
include("test_alternative_problems.jl")


text:
include("rw_stat.jl")

dim = length(ARGS)>0 ? parse(ARGS[1]) : 1
num = length(ARGS)>1 ? parse(ARGS[2]) : 1000
nsteps = length(ARGS)>2 ? parse(Int, ARGS[3]) : 1000
filename = length(ARGS)>3 ? ARGS[4] : ""

if dim == 1
    P = Point1D
elseif dim == 2
    P = Point2D
else
    error("dimension error!")
end

if isempty(filename)
    io = stdout
else
    io = open(filename, "w")
end

means, vars = randomwalk(P, num, nsteps)

for (i, (m, v)) in enumerate(zip(means, vars))
    println(io, i-1, " ", m, " ", v)
end

max_stars_repo_name: isuke/JuliaBook-Samples

max_issues_repo_name: isuke/JuliaBook-Samples

max_forks_repo_name: isuke/JuliaBook-Samples


text:
using StatsBase
using JLD2

function daily!(f, counts::AbstractVector{T} where T<:Integer, arr)
  fill!(counts, 0)
  max_day = length(counts)
  for a in arr
    t = f(a)
    if ismissing(t) || !isfinite(t)
      continue
    end
    day = floor(Int, t) + 1
    if day <= max_day
      counts[day] += 1
    end
  end
  counts
end

daily!(counts, arr) = daily!(identity, counts, arr)
daily(f::Function, times) = daily!(f, fill(0, ceil(Int, maximum(f, times))), times)
daily(times, max_time=ceil(Int, maximum(times))) = daily!(fill(0, max_time), times)

struct DailyTrajectories <: Output
  file::JLD2.JLDFile
end

DailyTrajectories(fname::AbstractString) = DailyTrajectories(JLD2.jldopen(fname, "w+", compress=true))
DailyTrajectories(fname::AbstractString, ::Integer) = DailyTrajectories(fname)

function pushtrajectory!(d::DailyTrajectories, trajectory_id::Integer, writelock::Base.AbstractLock, state::MocosSim.SimState, params::MocosSim.SimParams, cb::DetectionCallback)
  try
    lock(writelock)
    trajectory_group = JLD2.Group(d.file, string(trajectory_id))
    save_daily_trajectories(trajectory_group, state, params, cb)
  finally
    unlock(writelock)
  end
  nothing
end

aftertrajectories(d::DailyTrajectories) = close(d.file)

function save_daily_trajectories(dict, state::MocosSim.SimState, params::MocosSim.SimParams, cb::DetectionCallback)
  max_days = MocosSim.time(state) |> floor |> Int
  num_individuals = MocosSim.numindividuals(state)

  infection_times = Vector{OptTimePoint}(missing, num_individuals)
  contact_kinds = Vector{MocosSim.ContactKind}(undef, num_individuals)

  for i in 1:num_individuals
    event = MocosSim.backwardinfection(state, i)
    kind = contactkind(event)
    contact_kinds[i] = kind
    infection_times[i] = ifelse(kind == MocosSim.NoContact, missing, time(event))
  end

  hospitalization_progressions = getproperty.(state.progressions, :severe_symptoms_time)
  recovery_progressions = getproperty.(state.progressions, :recovery_time)
  death_progressions = getproperty.(state.progressions, :death_time)
  release_progressions = coalesce.(recovery_progressions, death_progressions)
  hospital_release_progressions = ifelse.(ismissing.(hospitalization_progressions), missing,  release_progressions)

  dict["daily_infections"] = daily(filter(!ismissing, infection_times), max_days)
  dict["daily_detections"] = daily(filter(!ismissing, cb.detection_times), max_days)
  dict["daily_deaths"] = daily(filter(!ismissing, infection_times.+death_progressions), max_days)
  dict["daily_hospitalizations"] = daily(filter(!ismissing, infection_times.+hospitalization_progressions), max_days)
  dict["daily_hospital_releases"] = daily(filter(!ismissing, infection_times.+hospital_release_progressions), max_days)
  for kind in instances(ContactKind)
    if kind != NoContact
      dict["daily_" * lowercase(string(kind))] = daily(infection_times[contact_kinds.==Int(kind)], max_days)
    end
  end
end


text:
export SeisPha

# ===========================================================================
# SeisPha
"""
    SeisPha()

IRIS-style seismic phase and pick container

Field     | Type    | Meaning         | SeisIO conventions/behavior
--------: |:------- |:--------------  | :----------
 amp      | Float64 | amplitude       | uses units of data source
 d        | Float64 | distance        | no unit conversion; can be m, km, or °
 ia       | Float64 | incidence angle | uses units of data source
 res      | Float64 | pick residual   |
 rp       | Float64 | ray parameter   |
 ta       | Float64 | takeoff angle   |
 tt       | Float64 | travel time     |
 unc      | Float64 | uncertainty     |
 pol      | Char    | polarity        |
 qual     | Char    | pick quality    | not (re)calculated

"""
mutable struct SeisPha
  amp ::Float64   # amplitude
  d   ::Float64   # distance
  ia  ::Float64   # incidence angle
  res ::Float64   # residual
  rp  ::Float64   # ray parameter
  ta  ::Float64   # takeoff angle
  tt  ::Float64   # travel time
  unc ::Float64   # uncertainty
  pol ::Char      # polarity
  qual::Char      # quality

  function SeisPha(
                    amp ::Float64 , # amplitude
                    d   ::Float64 , # distance
                    ia  ::Float64 , # incidence angle
                    res ::Float64 , # residual
                    rp  ::Float64 , # ray parameter
                    ta  ::Float64 , # takeoff angle
                    tt  ::Float64 , # travel time
                    unc ::Float64 , # uncertainty
                    pol ::Char    , # polarity
                    qual::Char      # quality
                    )
    return new(amp, d, ia, res, rp, ta, tt, unc, pol, qual)
  end
end

SeisPha( ;
        amp ::Float64   = zero(Float64),
        d   ::Float64   = zero(Float64),
        ia  ::Float64   = zero(Float64),
        res ::Float64   = zero(Float64),
        rp  ::Float64   = zero(Float64),
        ta  ::Float64   = zero(Float64),
        tt  ::Float64   = zero(Float64),
        unc ::Float64   = zero(Float64),
        pol ::Char      = ' ',
        qual::Char      = ' '
        ) = SeisPha(amp, d, ia, res, rp, ta, tt, unc, pol, qual)

function write(io::IO, Pha::SeisPha)
  write(io, Pha.amp)
  write(io, Pha.d)
  write(io, Pha.ia)
  write(io, Pha.res)
  write(io, Pha.rp)
  write(io, Pha.ta)
  write(io, Pha.tt)
  write(io, Pha.unc)
  write(io, Pha.pol)
  write(io, Pha.qual)
  return nothing
end

read(io::IO, ::Type{SeisPha}) =
  SeisPha(read(io, Float64),
          read(io, Float64),
          read(io, Float64),
          read(io, Float64),
          read(io, Float64),
          read(io, Float64),
          read(io, Float64),
          read(io, Float64),
          read(io, Char),
          read(io, Char)
          )

function isempty(Pha::SeisPha)
  q::Bool = min(getfield(Pha, :pol) == ' ', getfield(Pha, :qual) == ' ')
  if q == true
    for f in (:amp, :d, :ia, :res, :rp, :ta, :tt, :unc)
      q = min(q, getfield(Pha, f) == zero(Float64))
    end
  end
  return q
end

function isequal(S::SeisPha, U::SeisPha)
  q::Bool = isequal(getfield(S, :pol), getfield(U, :pol))
  if q == true
    for f in (:amp, :d, :ia, :res, :rp, :ta, :tt, :unc, :pol, :qual)
      q = min(q, getfield(S,f) == getfield(U,f))
    end
  end
  return q
end
==(S::SeisPha, U::SeisPha) = isequal(S, U)
sizeof(P::SeisPha) = 146


text:
using LinearAlgebra
using Random

export rfi
export rfi_b

"""
    rfi(nep,nept,[λ=0,][errmeasure,][tol=eps()*100,][maxit=100,][v=randn,][u=randn,][logger=0,][linsolvecreator=DefaultLinSolverCreator(),])

This is an implementation of the two-sided Rayleigh functional Iteration (RFI) to compute an eigentriplet of the problem specified by `nep`.
This method requires the transpose of the NEP, specified in `nept`.
`λ`, `u` and `v` are initial guesses for the eigenvalue, the right eigenvector and the left eigenvector respectively.
See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> nep=nep_gallery("dep0");
julia> nept=DEP([nep.A[1]',nep.A[2]'])
julia> λ,v,u=rfi(nep,nept)
julia> compute_resnorm(nep,λ,v) # v is a right eigenvector
3.0171586304599647e-16
julia> compute_resnorm(nept,λ,u) # u is a right eigenvector
7.145081514857076e-17
```
# Reference
*  Algorithm 4 in  Schreiber, Nonlinear Eigenvalue Problems: Newton-type Methods and Nonlinear Rayleigh Functionals, PhD thesis, TU Berlin, 2008.
"""
rfi(nep::NEP, nept::NEP; kwargs...) = rfi(ComplexF64,nep, nept,;kwargs...)
function rfi(::Type{T},
            nep::NEP,
            nept::NEP;
            errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
            tol = eps(real(T))*1000,
            maxit=100,
            λ::Number = zero(T),
            v::Vector = randn(T,size(nep,1)),
            u::Vector = randn(T,size(nep,1)),
            linsolvercreator=BackslashLinSolverCreator(), # Note: Not default. It's better to use backslash
            inner_solver=@default_compute_rf_inner_solver(nep),
            logger=0) where {T <: Number}

        @parse_logger_param!(logger)

        err = Inf

        #Normalize v and u
        normalize!(v)
        normalize!(u)

        for k=1:maxit
            err = estimate_error(errmeasure,λ,u)

            if(err < tol)
                return λ,u,v
            end

            push_iteration_info!(logger,k,err=err,λ=λ,v=v, continues=true);
            push_info!(logger," u=$u");

            local linsolver=create_linsolver(linsolvercreator,nep,λ)
            local linsolver_t=create_linsolver(linsolvercreator,nept,λ)

            x = lin_solve(linsolver,compute_Mlincomb(nep,λ,u,[T(1)],1),tol = tol)
            u[:] = normalize(x)

            y = lin_solve(linsolver_t,compute_Mlincomb(nept,λ,v,[T(1)],1),tol = tol)
            v[:] = normalize(y)

            λ_vec = compute_rf(T,nep, u, inner_solver, y=v)
            λ = closest_to(λ_vec,  λ)
        end

        msg="Number of iterations exceeded. maxit=$(maxit)."
        throw(NoConvergenceException(λ,u,err,msg))
end

"""
    rfi_b(nep,nept,[λ=0,][errmeasure,][tol=eps()*100,][maxit=100,][v=randn,][u=randn,][logger=0,][linsolvecreator=DefaultLinSolverCreator(),])

This is an implementation of the two-sided Rayleigh functional Iteration(RFI)-Bordered version to compute an eigentriplet of the problem specified by `nep`.
This method requires the transpose of the NEP, specified in `nept`.
`λ`, `u` and `v` are initial guesses for the eigenvalue, the right eigenvector and the left eigenvector respectively.
See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> nep=nep_gallery("dep0");
julia> nept=DEP([nep.A[1]',nep.A[2]'])
julia> λ,v,u=rfi_b(nep,nept,v=ones(ComplexF64,size(nep,1)))
julia> compute_resnorm(nep,λ,v) # v is a right eigenvector
5.343670589284583e-15
```

# Reference
*  Algorithm 5 in  Schreiber, Nonlinear Eigenvalue Problems: Newton-type Methods and Nonlinear Rayleigh Functionals, PhD thesis, TU Berlin, 2008.

"""
rfi_b(nep::NEP, nept::NEP; kwargs...) = rfi_b(ComplexF64,nep, nept,;kwargs...)
function rfi_b(::Type{T},
            nep::NEP,
            nept::NEP;
            errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
            tol = eps(real(T))*1000,
            maxit=100,
            λ::Number = zero(T),
            v::Vector = randn(T,size(nep,1)),
            u::Vector = randn(T,size(nep,1)),
            inner_solver=@default_compute_rf_inner_solver(nep),
            logger=0) where {T <: Number}

        @parse_logger_param!(logger)

        err = Inf
        #Normalize v and u
        normalize!(v)
        normalize!(u)

        for k=1:maxit
            err = estimate_error(errmeasure,λ,u)

            if(err < tol)
                return λ,u,v
            end

            push_iteration_info!(logger,k,err=err,λ=λ,v=v, continues=true);
            push_info!(logger," u=$u");

            #Construct C_k
            C = [compute_Mder(nep,λ,0) compute_Mlincomb(nep,λ,u,[T(1)],1);v'*compute_Mder(nep,λ,1) T(0.0)]

            l1 = C\-[compute_Mlincomb(nep,λ,u,[T(1)],0);0]
            s = l1[1:end-1]
            u[:] = normalize(u+s)

            l2 = C\-[compute_Mlincomb(nept,λ,v,[T(1)],0);0]
            t = l2[1:end-1]
            v[:] = normalize(v+t)

            λ_vec = compute_rf(T,nep, u, inner_solver, y=v)
            λ = closest_to(λ_vec,  λ)
        end

        msg="Number of iterations exceeded. maxit=$(maxit)."
        throw(NoConvergenceException(λ,u,err,msg))
end