text:
# MMA Algorithms

"""
    MMA87

The original method of moving asymptotes (MMA) algorithm from the [1987 paper](https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620240207).
"""
@params struct MMA87 <: AbstractOptimizer
    dualoptimizer
end
function MMA87(; dualoptimizer = Optim.GradientDescent())
    return MMA87(dualoptimizer)
end

const MMA = MMA87

"""
    MMA02

The globally convergent method of moving asymptotes (MMA) algorithm from the [2002 paper](https://epubs.siam.org/doi/abs/10.1137/S1052623499362822).
"""
@params struct MMA02 <: AbstractOptimizer
    dualoptimizer
end
function MMA02(; dualoptimizer = Optim.GradientDescent())
    return MMA02(dualoptimizer)
end

const GCMMA = MMA02

"""
    MMAOptions

A struct that stores all the options of the MMA algorithms. Th following are the fields of `MMAOptions`:
 - `maxiter`: the maximum number of inner iterations. For `MMA87`, there is 1 inner iteration per outer iteration.
 - `outer_maxiter`: the maximum number of outer iterations.
 - `maxinner`: the maximum number of inner iterations per outer iteration of [`MMA02`](@ref). Not applicable for [`MMA87`](@ref).
 - `tol`: a tolerance struct of type [`Tolerance`](@ref).
 - `s_init`: defined in the original [`MMA02`](@ref) paper.
 - `s_incr`: defined in the original [`MMA02`](@ref) paper.
 - `s_decr`: defined in the original [`MMA02`](@ref) paper.
 - `store_trace`: if true, a trace will be stored.
 - `dual_options`: the options passed to the dual optimizer from [`Optim.jl`](https://github.com/JuliaNLSolvers/Optim.jl).
"""
@with_kw mutable struct MMAOptions{T, Ttol <: Tolerance, TSubOptions <: Optim.Options}
    maxiter::Int = 1000
    outer_maxiter::Int = 10^8
    maxinner::Int = 10
    tol::Ttol = Tolerance()
    s_init::T = 0.5
    s_incr::T = 1.2
    s_decr::T = 0.7
    store_trace::Bool = false
    show_trace::Bool = false
    auto_scale::Bool = false
    keep_best::Bool = false
    dual_options::TSubOptions = Optim.Options(
        allow_outer_f_increases = true,
        allow_f_increases = true,
        iterations = 1000,
        outer_iterations=1000,
    )
end

"""
    Trace

A struct that stores the history of solutions.
"""
struct Trace
    tr::Vector
end

"""
    MMAWorkspace

A struct that stores all the intermediate states and memory allocations needed for the optimization. The following are the fields of the struct:
 - `model`: the original instance of [`MMAModel`](@ref) being optimized.
 - `dualmodel`: an instance of [`DualMMAModel`](@ref) which resembles the dual of the MMA approximation of the original model.
 - `x0`: the initial primal solution.
 - `solution`: an instance of [`Solution`](@ref) that resembles the current solution in the optimization algorithm.
 - `σ`: the workspace move limit, explained in [`MMAApprox`](@ref).
 - `ρ`: the `ρ` parameter as explained in [`MMAApprox`](@ref).
 - `tempx`: a temporary vector used to store the 2nd previous primal solution.
 - `options`: an instance of [`MMAOptions`](@ref) that resembles the options of the MMA algorithm.
 - `convcriteria`: an instance of [`ConvergenceCriteria`](@ref) that specifies the convergence criteria of the MMA algorithm.
 - `callback`: a function that is called on `solution` in every iteration of the algorithm. This can be used to store information about the optimization process.
 - `optimizer`: an instance of [`AbstractOptimizer`](@ref) such as `MMA87()` or `MMA02()` that specifies the variant of MMA used to optimize the model.
 - `suboptimizer`: the dual optimization algorithm used to optimize the barrier problem. This should be an [`Optim.jl`](https://github.com/JuliaNLSolvers/Optim.jl) optimizer.
 - `tracing`: a Boolean that when true stores the trace of solutions evaluated during the optimization process.
 - `outer_iter`: the current outer iteration index. [`MMA87`](@ref) has 1 inner iteration per outer iteration. [`MMA02`](@ref) does multiple inner iterations per outer iteration.
 - `iter`: the current inner iteration index.
 - `fcalls`: the current number of objective and constraint function calls.
"""
@params mutable struct MMAWorkspace <: Workspace
    model::VecModel
    dualmodel::MMADualModel
    x0::AbstractVector
    solution::Solution
    σ::AbstractVector
    ρ::AbstractVector
    tempx::AbstractVector
    convcriteria::ConvergenceCriteria
    callback::Function
    optimizer::AbstractOptimizer
    options
    trace::Trace
    outer_iter::Int
    iter::Int
	fcalls::Int
end
function MMAWorkspace(
    model::VecModel,
    optimizer::AbstractOptimizer,
    x0::AbstractVector;
    options = default_options(model, optimizer),
    convcriteria::ConvergenceCriteria = KKTCriteria(),
    plot_trace::Bool = false,
    show_plot::Bool = plot_trace,
    save_plot = nothing,
    callback::Function = plot_trace ? LazyPlottingCallback(; show_plot = show_plot, save_plot = save_plot) : NoCallback(),
    kwargs...,
)
    T = eltype(x0)
    init!(model)
    dualmodel = MMADualModel(MMAApproxModel(model, x0; kwargs...))

    # Convergence
    λ = ones(getdim(getineqconstraints(model)))
    solution = Solution(dualmodel, λ)
    assess_convergence!(solution, model, options.tol, convcriteria)
    correctsolution!(solution, model, options)

    # Trace
    trace = Trace([])

    # Iteraton counter
    fcalls, outer_iter, iter = 1, 0, 0
    σ = similar(x0)
    ρ = similar(λ, length(λ) + 1)
    tempx = similar(x0)

    return MMAWorkspace(
        model,
        dualmodel,
        x0,
        solution,
        σ,
        ρ,
        tempx,
        convcriteria,
        callback,
        optimizer,
        options,
        trace,
        outer_iter,
        iter,
        fcalls,
    )
end

function reset!(w::MMAWorkspace, x0 = nothing)
    @unpack solution = w
    outer_iter, iter, fcalls = 0, 0, 0, 0
    @pack! w = fcalls, iter, outer_iter
    if x0 !== nothing
        w.x0 .= x0
        w.tempx .= solution.prevx
        solution.prevx .= solution.x
        solution.x .= x0
        updateapprox!(w.dualmodel, x0)
    end
    assess_convergence!(w)
    return w
end

default_options(model::VecModel, alg::MMA87) = MMAOptions()
default_options(model::VecModel, alg::MMA02) = MMAOptions()

init!(model::VecModel) = model

function Workspace(model::VecModel, optimizer::Union{MMA87, MMA02}, args...; kwargs...)
    return MMAWorkspace(model, optimizer, args...; kwargs...)
end

function optimize!(workspace::MMAWorkspace)
    @unpack dualmodel, solution, convcriteria = workspace
    @unpack callback, optimizer, options, trace = workspace
    @unpack x0, σ, ρ, outer_iter, iter, fcalls = workspace
    @unpack dualoptimizer = optimizer
    @unpack dual_options, maxiter, outer_maxiter, auto_scale = options
    @unpack prevx, x, g, λ = solution
    best_solution = deepcopy(solution)

    T = eltype(x)
    ni, nj = length(λ), length(x)

    # Original model
    model = dualmodel |> getparent |> getparent

    # Approximate MMA model
    approxmodel = getparent(dualmodel)

    # Dual objective to be maximized
    dualobj = getobjective(dualmodel)

    # Optim-compatible objective
    # Calculates negative the dual objective and its gradient
    optimobj = getoptimobj(dualobj, false)

    # Lower and upper bounds on the dual variables
    λl = zeros(ni) .+ 1e-10
    λu = fill(Inf, ni)
    λ .= 1

    # Callback, e.g. a trace plotting callback
    callback(solution)

    # Initialize the workspace's trust region σ
    initializeσ!(workspace)
    # Initialize the workspace's lift ρ, only used in MMA02
    initializeρ!(workspace)

    while (!hasconverged(solution) || outer_iter == 0) && iter < maxiter && outer_iter < outer_maxiter
        outer_iter += 1

        # Adapt workspace's trust region σ using s_incr and s_decr
        outer_iter > 2 && updateσ!(workspace)
        # Decrease the workspace's lift ρ, only used in MMA02
        decreaseρ!(workspace, optimizer)

        # Update the model's σ
        setσ!(dualmodel, σ)
        # Update the model's ρ, only used in MMA02
        setρ!(dualmodel, ρ)

        # Update the convex approximation at the current x
        updateapprox!(dualmodel)

        # Assume the convex approximation is not an upper bound approximation at x
        upperbounded = false
        prev_iter = iter
        while (!hasconverged(solution) || outer_iter == 1) && !upperbounded && iter < maxiter && iter < prev_iter + options.maxinner
            iter += 1
            # Solve the dual problem by minimizing negative the dual objective value
            if length(λ) > 0
                λ .= 1.0
                λ .= Optim.optimize(
                    Optim.only_fg!(optimobj),
                    λl,
                    λu,
                    λ,
                    Optim.Fminbox(dualoptimizer),
                    dual_options,
                ).minimizer
            end
            if debugging[]
                @show λ
            end

            # Revaluate the optimal primal solution at the optimal dual solution
            dualobj(λ)

            # Update the x vector to the new one and update the previous 2 xs
            updatex!(workspace)

            # Evaluate the approximate objective and constraints at the optimal x
            approxfg = getapproxfg(dualobj)

            # Evaluates the exact objective and constraints and their gradients at the optimal x
            optimalx = getoptimalx(dualmodel)

            fg, ∇fg = value_jacobian(
                approxmodel.objective_ineq_constraints,
                optimalx,
            )
            fcalls += 1

            # Scale the objective appropriately
            # Reduces the chance of overflow or underflow in λs
            # Helps with the convergence
            auto_scale && iter < 5 && scaleobjective!(model, fg, ∇fg, approxfg, λ)

            # Check if the approximation is an upper bound at the current x.
            # If not, increase ρ. Only used in MMA02
            if iter > 1
                upperbounded = increaseρ!(workspace, fg, approxfg, optimizer)
                upperbounded || setρ!(dualmodel, ρ)
            end
            # Update the convex approximation at the current x
            updateapprox!(dualmodel, optimalx, fg, ∇fg)
            # Update the objective and constraint values and gradients in `solution`
            updatefg!(solution, fg, ∇fg)

            # Check if the algorithm has converged
            assess_convergence!(solution, model, options.tol, convcriteria)

            # Callback, e.g. a trace plotting callback
            callback(solution)
        end
        @pack! workspace = outer_iter, iter, fcalls

        if options.keep_best
            best_solution = get_best_solution(solution, best_solution, options)
        else
            best_solution = deepcopy(solution)
        end
        hasconverged(best_solution) && break

        # Print some trace if flag is on
        # @mmatrace()
    end
    @pack! workspace = outer_iter, iter, fcalls
    workspace.solution = best_solution

    # Reset the objective scaling factor to 1
    set_objective_multiple!(model, 1)
    callback(best_solution, update = true)

    results = GenericResult(
        optimizer,
        x0,
        best_solution.x,
        best_solution.f,
        iter,
        iter == options.maxiter,
        options.tol,
        best_solution.convstate,
        fcalls,
    )
    return results
end

"""
    scaleobjective!(model::VecModel, fg, ∇fg, approxfg, λ)

Scales the objective of `model` using the ratio of the ∞ norms of the constraint and objective values and gradients. After scaling, the ∞ norms will be the same.
"""
function scaleobjective!(model::VecModel, fg, ∇fg, approxfg, λ)
    @views normratio = norm(∇fg[2:end,:], Inf) / norm(∇fg[1,:], Inf)
    @views normratio = max(normratio, norm(fg[2:end], Inf) / abs(fg[1]))
    normratio = isfinite(normratio) ? normratio : one(normratio)
    set_objective_multiple!(model, get_objective_multiple(model) * normratio)
    fg[1] = fg[1] * normratio
    approxfg[1] = approxfg[1] * normratio
    ∇fg[1,:] .= ∇fg[1,:] .* normratio
    scalequadweight!(model, normratio)
    λ .*= normratio
    return model
end

scalequadweight!(model::VecModel, s::Real) = model

function correctsolution!(solution::Solution, model::VecModel, options::MMAOptions)
    solution.f = solution.f / get_objective_multiple(model)
    return solution
end

function get_best_solution(solution::Solution, best_solution::Solution, options)
    best_infeas = best_solution.convstate.infeas
    if max(solution.convstate.infeas, solution.convstate.kkt_residual) <=
        max(best_solution.convstate.infeas, best_solution.convstate.kkt_residual)
        best_solution = deepcopy(solution)
    end
    return best_solution
end

function updatex!(workspace::MMAWorkspace)
    @unpack tempx, solution, dualmodel = workspace
    @unpack prevx, x = solution
    tempx .= prevx
    prevx .= x
    x .= getoptimalx(dualmodel)
    return workspace
end

function updatefg!(solution::Solution, fg::AbstractVector, ∇fg::AbstractMatrix)
    solution.prevf = solution.f
    solution.f = fg[1]
    solution.g .= fg[2:end]
    solution.∇f .= ∇fg[1, :]
    solution.∇g .= ∇fg[2:end, :]
    return solution
end

function initializeσ!(workspace::MMAWorkspace)
    @unpack σ, optimizer, options, model = workspace
    T = eltype(σ)
    map!(σ, 1:length(σ)) do j
        diff = getmax(model, j) - getmin(model, j)
        if !isfinite(diff)
            diff = 1000 * one(T)
        end
        return options.s_init * diff
    end
    return workspace
end

function initializeρ!(workspace::MMAWorkspace)
    @unpack ρ = workspace
    ρ .= 0
    return workspace
end

function updateσ!(workspace::MMAWorkspace)
    @unpack σ, solution, tempx, optimizer, options, model = workspace
    @unpack x, prevx = solution
    T = eltype(x)
    map!(σ, 1:length(x)) do j
        if x[j] == prevx[j] || prevx[j] == tempx[j]
            d = σ[j]
        elseif xor(x[j] > prevx[j], prevx[j] > tempx[j])
            d = σ[j] * options.s_decr
        else
            d = σ[j] * options.s_incr
        end
        diff = getmax(model, j) - getmin(model, j)
        if !isfinite(diff)
            diff = 1000 * one(T)
        end
        min = diff/100
        max = 10diff
        if d <= min
            return min
        elseif d >= max
            return max
        else
            return d
        end
    end
    return workspace
end

decreaseρ!(workspace::MMAWorkspace, ::MMA87) = workspace
function decreaseρ!(workspace::MMAWorkspace, ::MMA02)
    @unpack ρ = workspace
    ρ .= max.(ρ ./ 10, 1e-5)
    return workspace
end

increaseρ!(workspace::MMAWorkspace, ::AbstractVector, ::AbstractVector, ::MMA87) = true
function increaseρ!(
    workspace::MMAWorkspace,
    fg::AbstractVector,
    approxfg::AbstractVector,
    ::MMA02,
)
    @unpack dualmodel, solution, σ, ρ = workspace
    @unpack x, prevx = solution
    T = eltype(x)
    w = zero(T)
    for j in 1:length(x)
        diff2 = (x[j] - prevx[j])^2
        w += diff2 / (σ[j]^2 - diff2)
    end
    w /= 2

    upperbounded = true
    for i in 1:length(ρ)
        if fg[i] > approxfg[i]
            upperbounded = false
            ρ[i] = min(10ρ[i], 1.1 * (ρ[i] + (fg[i] - approxfg[i]) / w))
        end
    end
    return upperbounded
end

function getoptimobj(obj, minimize = true)
    optimobj(z) = optimobj(1.0, nothing, z)
    function optimobj(F, G, z)
        if G !== nothing
            val, grad = value_gradient(obj, z)
            if minimize
                G[:] .= grad
            else
                G[:] .= .-grad
            end
            if F !== nothing
                if minimize
                    return val
                else
                    return -val
                end
            end
        end
        # No gradient necessary, just return the log joint.
        if F !== nothing
            if minimize
                return obj(z)
            else
                return -obj(z)
            end
        end
        return nothing
    end
    return optimobj
end


text:
using ProbNumODE
using DifferentialEquations
using Plots
using ProgressMeter
using LinearAlgebra

FILEDIR = dirname(@__FILE__)

prob = ProbNumODE.remake_prob_with_jac(van_der_pol(p=[1e6]))

# The EKF1-IWP3 is able to solve this very stiff problem
println("Timed solve of the VdP problem with EKS1-IBM3:")
@time sol = solve(prob, EKF1(), q=3, smooth=true); sol.destats

# Visualizations
# 1. Solution
p1 = plot(sol.t,
          max.(min.(ProbNumODE.stack(sol.u), 7.0), -7.0),
          xlabel="t", ylabel="y(t)", label="",
          );
savefig(p1, joinpath(FILEDIR, "vdp_solution.png"))

# 2. Step sizes
p2 = plot(sol.t[1:end-1], sol.t[2:end]-sol.t[1:end-1],
          yscale=:log10, xlabel="t", ylabel="dt", label="");
savefig(p2, joinpath(FILEDIR, "vdp_stepsizes.png"))

# 3. Errors
println("Timed very-low-tolerance-solve of the VdP problem with RadauIIA5:")
@time appxsol = solve(
    prob, RadauIIA5(), abstol=1e-14, reltol=1e-14, maxiters=1e7); appxsol.destats
errs = sol.u - appxsol.(sol.t)
ProbNumODE.stack(errs) ./ ProbNumODE.stack(appxsol.(sol.t))

p3 = plot(sol.t, replace(abs.(ProbNumODE.stack(errs)), 0=>missing),
          color=[1 2], label="", yscale=:log10, xlabel="t", ylabel="Error");
plot!(p3, sol.t, replace(sqrt.(ProbNumODE.stack(diag.(sol.pu.Σ))), 0=>missing),
      color=[1 2], linestyle=:dash, label="");
savefig(p3, joinpath(FILEDIR, "vdp_errors.png"))

# Additional plot: Phase plot
u = ProbNumODE.stack(sol.u)
p4 = plot(u[:, 1], u[:, 2], xlabel="y_1(t)", ylabel="y_2(t)", label="");
savefig(p4, joinpath(FILEDIR, "vdp_solution_phase.png"))

println("Plots saved to $FILEDIR")

ext: jl

max_stars_repo_path: experiments/1 Stiff Van der Pol/main.jl

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_path: experiments/1 Stiff Van der Pol/main.jl

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_path: experiments/1 Stiff Van der Pol/main.jl

max_forks_repo_name: nathanaelbosch/capos


text:
"""
    StuTLik <: Likelihood

Student-t likelihood (a.k.a. non-standardized [Student's
t-distribution](https://en.wikipedia.org/wiki/Student's_t-distribution))
```math
p(y | f, σ) = Γ((ν + 1)/2)/[Γ(ν/2)√(πν)σ](1 + 1/ν((y - f)/σ)²)^{-(ν + 1)/2}
```
with degrees of freedom ``ν ∈ N₀``, where scale ``σ`` is a non-fixed hyperparameter and
``f`` is the latent Gaussian process.
"""
mutable struct StuTLik <: Likelihood
    "Degrees of freedom"
    ν::Int
    "Scale"
    σ::Float64
    "Priors for likelihood parameters"
    priors::Array

    """
        StuTLik(ν::Int, lσ::Float64)

    Create a `StuTLik` with degrees of freedom `ν` and scale `exp(lσ)`.
    """
    StuTLik(ν::Int, lσ::Float64) = new(ν, exp(lσ), [])
end

#log of probability density
function log_dens(studentT::StuTLik, f::VecF64, y::VecF64)
    ν = studentT.ν
    σ = studentT.σ
    c = lgamma(0.5*(ν+1)) - lgamma(0.5*ν) - 0.5*log(pi*ν) - log(σ)
    return [c - (0.5*(ν+1))*log(1+(1/ν)*((yi-fi)/σ)^2) for (fi,yi) in zip(f,y)]
end

#derivative of log pdf wrt latent function
function dlog_dens_df(studentT::StuTLik, f::VecF64, y::VecF64)
    ν = studentT.ν
    σ = studentT.σ
    return [(ν+1)*(yi-fi)/(ν*σ^2 + (yi-fi)^2) for (fi,yi) in zip(f,y)]
end

#derivative of log pdf wrt to parameters
function dlog_dens_dθ(studentT::StuTLik, f::VecF64, y::VecF64)
    ν = studentT.ν
    σ = studentT.σ
    return σ*[-1/σ+(ν+1)*(yi-fi)^2/(ν*σ^3 + σ*(yi-fi)^2) for (fi,yi) in zip(f,y)]
end

#mean and variance under likelihood
mean_lik(studentT::StuTLik, f::VecF64) = f
var_lik(studentT::StuTLik, f::VecF64) = studentT.σ^2*(studentT.ν/(studentT.ν-2.0))*ones(length(f))

function set_params!(studentT::StuTLik, hyp::VecF64)
    length(hyp) == 1 || throw(ArgumentError("Student-t likelihood has only one free parameter"))
    studentT.σ = exp(hyp[])
end

get_params(studentT::StuTLik) = Float64[log(studentT.σ)]
num_params(studentT::StuTLik) = 1


text:
# promotion.jl
#

# # Promotion
function promote(a::TaylorModelN{N,T,S}, b::R) where {N, T<:Real, S<:Real, R<:Real}
    orderTMN = get_order(a[0])
    apol, bb = promote(a.pol, b)
    return (TaylorModelN(apol, a.rem, a.x0, a.dom), TaylorModelN(bb, zero(a.rem), a.x0, a.dom))
end
promote(b::R, a::TaylorModelN{N,T,S}) where {N, T<:Real, S<:Real, R<:Real} =
    reverse(promote(a,b))

#
function promote(a::TaylorModelN{N,T,S}, b::TaylorN{R}) where {N, T, S, R}
    RR = promote_type(T,R)
    aa = TaylorModelN(convert(TaylorN{RR},a.pol), a.rem, a.x0, a.dom)
    bb = TaylorModelN(convert(TaylorN{RR},b), 0..0, a.x0, a.dom)
    (aa, bb)
end
promote(b::TaylorN{R}, a::TaylorModelN{N,T,S}) where {N, T, S, R} = reverse( promote(a, b) )

function promote(a::Taylor1{TaylorModelN{N,Interval{T},S}}, b::Taylor1{Interval{T}}) where
        {N, T<:Real, S<:Real}

    orderTMN = get_order(a[0])
    bTN = Array{TaylorModelN{N,Interval{T},S}}(undef, get_order(a)+1)
    unoTM = TaylorModelN(Interval(T(1), T(1)), orderTMN, a[0].x0, a[0].dom)
    @inbounds for ord = 1:length(bTN)
        bTN[ord] = b.coeffs[ord] * unoTM
    end
    (a, Taylor1(bTN))
end
promote(b::Taylor1{Interval{T}}, a::Taylor1{TaylorModelN{N,Interval{T},S}}) where
    {N, T<:Real, S<:Real} = reverse( promote(a, b) )

function promote(a::Taylor1{TaylorModelN{N,Interval{T},S}}, b::T) where
        {N, T<:Real, S<:Real}

    orderTMN = get_order(a[0])
    bTN = Array{TaylorModelN{N,Interval{T},S}}(undef, get_order(a)+1)
    unoTM = TaylorModelN(Interval(one(T), one(T)), orderTMN, a[0].x0, a[0].dom)
    bTN[1] = b * unoTM
    @inbounds for ord = 2:length(bTN)
        bTN[ord] = zero(b)*unoTM
    end
    (a, Taylor1(bTN))
end
promote(b::T, a::Taylor1{TaylorModelN{N,Interval{T},S}}) where
        {N, T<:Real, S<:Real} = reverse( promote(a, b) )

function promote(a::Taylor1{TaylorModelN{N,T,S}}, b::R) where
        {N, T<:Real, S<:Real, R<:Real}

    orderTMN = get_order(a[0])
    TT = promote_type(T,R)
    aTN = Array{TaylorModelN{N,TT,S}}(undef, get_order(a)+1)
    bTN = Array{TaylorModelN{N,TT,S}}(undef, get_order(a)+1)
    unoTT = one(TT)
    zeroTT = zero(TT)
    @inbounds for ord = 1:length(aTN)
        aTN[ord] = unoTT * a.coeffs[ord]
        bTN[ord] = zeroTT * a.coeffs[ord]
    end
    bTN[1] += b
    (Taylor1(aTN), Taylor1(bTN))
end
promote(b::R, a::Taylor1{TaylorModelN{N,T,S}}) where
        {N, T<:Real, S<:Real, R<:Real} = reverse( promote(a, b) )

# function promote(a::Taylor1{TaylorModelN{N,T,S}}, b::Taylor1{S}) where {N,T<:Real,S<:Real}
#     (a, Taylor1(TaylorModelN(interval(b), get_order(a), a.x0, a.dom), b.order))
# end
# promote(b::Taylor1{Interval{T}}, a::Taylor1{TaylorModelN{N,Interval{T},S}}) where
#     {N, T<:Real, S<:Real} = promote(a, b)

for TM in tupleTMs
    @eval promote(a::$TM{T,S}, b::T) where {T,S} =
        (a, $(TM)(Taylor1([b], get_order(a)), zero(a.rem), a.x0, a.dom))
    @eval promote(b::T, a::$TM{T,S}) where {T,S} = reverse( promote(a,b) )
    #
    @eval promote(a::$TM{T,S}, b::S) where {T,S} =
        (a, $TM(Taylor1([convert(T, b)], get_order(a)), zero(a.rem), a.x0, a.dom))
    @eval promote(b::S, a::$TM{T,S}) where {T,S} = reverse( promote(a,b) )
    #
    @eval promote(a::$TM{T,S}, b::R) where {T,S,R} = promote(a, convert(S, b))
    @eval promote(b::R, a::$TM{T,S}) where {T,S,R} = reverse( promote(a,b) )
    #
    @eval function promote(a::$TM{TaylorModelN{N,T,S},S}, b::T) where {N,T,S}
        tmN = TaylorModelN(b, get_order(a.pol[0]), a.pol[0].x0, a.pol[0].dom)
        (a, $TM(Taylor1([tmN], get_order(a)), zero(a.rem), a.x0, a.dom))
    end
    @eval promote(b::T, a::$TM{TaylorModelN{N,T,S},S}) where {N,T,S} = reverse( promote(a,b) )
end


text:
@testset "PositivePoisson" begin
using GLM: linkinv, linkfun, mueta, inverselink

λ0=3.4
pp = PositivePoisson(λ0)
@test params(pp) == (λ0,)
@test rate(pp) == λ0
@test mean(pp) == λ0 / (1-exp(-λ0))

xs = 1:10
@test pdf.(pp,xs) ≈ exp.(logpdf.(pp,xs)) atol=1e-4

xs = 1:100:10000
lp1=[logpdf(Poisson(λ0),x)::Float64 for x=xs]
lp2=[HurdleDMR.logpdf_approx(Poisson(λ0),x)::Float64 for x=xs]
@test lp1 ≈ lp2

lp1=[logpdf(PositivePoisson(λ0),x)::Float64 for x=xs]
lp2=[logpdf_exact(PositivePoisson(λ0),x)::Float64 for x=xs]
@test lp1 ≈ lp2

le1 = HurdleDMR.logexpm1.(xs)
le1big = HurdleDMR.logexpm1.(big.(xs))
le1big2 = broadcast(x->log(exp(x)-one(x)),big.(xs))
@test le1 ≈ le1big
@test le1 ≈ le1big2

ηs=-10:1.0:10
μs=broadcast(η->linkinv(LogProductLogLink(),η),ηs)
ηscheck=broadcast(μ->linkfun(LogProductLogLink(),μ),μs)
@test ηs ≈ ηscheck

μs=broadcast(η->mueta(LogProductLogLink(),η),ηs)
μscheck=broadcast(η->inverselink(LogProductLogLink(),η)[2],ηs)
@test μs ≈ μscheck

loglik(y, μ) = GLM.loglik_obs(PositivePoisson(λ0), y, μ, one(y), 0)

μs=1.01:10.0:1000.0
ηs=broadcast(μ->linkfun(LogProductLogLink(),μ),μs)
μscheck=broadcast(η->linkinv(LogProductLogLink(),η),ηs)
@test μs ≈ μscheck
#verify works for large μ
ys = round.(Int,μs) .+ 1.0
devresids = devresid.(PositivePoisson(λ0), ys, μs)
@test all(isfinite,devresids)
logliks = loglik.(ys, μs)
@test all(isfinite,logliks)

@test isinf(devresid(PositivePoisson(λ0), 3.0, 0.1))
@test iszero(devresid(PositivePoisson(λ0), 0.3, 0.1))

ys = fill(1.0, size(μs))
devresids1 = devresid.(PositivePoisson(λ0), ys, μs)
@test all(isfinite,devresids1)
logliks1 = loglik.(ys, μs)
@test all(isfinite,logliks1)

μsbig = big.(μs)
ηs=broadcast(μ->linkfun(LogProductLogLink(),μ),μsbig)
μscheck=broadcast(η->linkinv(LogProductLogLink(),η),ηs)
@test μsbig ≈ μscheck
#verify works for large μ
ysbig = round.(BigInt,μsbig) .+ 1.0
devresidsbig = devresid.(PositivePoisson(λ0), ysbig, μsbig)
@test all(isfinite,devresidsbig)
@test devresids ≈ Float64.(devresidsbig)

logliksbig = loglik.(ysbig, μsbig)
@test all(isfinite,logliksbig)
@test logliks ≈ Float64.(logliksbig)

ysbig = fill(big"1.0", size(μsbig))
devresidsbig1 = devresid.(PositivePoisson(λ0), ysbig, μsbig)
@test all(isfinite,devresidsbig1)
@test devresids1 ≈ Float64.(devresidsbig1)
logliksbig1 = loglik.(ysbig, μsbig)
@test all(isfinite,logliksbig1)

# μs close to 1.0
μs=big.(collect(1.0+1e-10:1e-10:1.0+1000*1e-10))
ηs=broadcast(μ->linkfun(LogProductLogLink(),μ),μs)
μscheck=broadcast(η->linkinv(LogProductLogLink(),η),ηs)
@test μs ≈ μscheck
#verify works for large μ
ys = round.(BigInt,μs) .+ big"1.0"
devresidsbig = devresid.(PositivePoisson(λ0), ys, μs)
@test all(isfinite,devresidsbig)
logliksbig = loglik.(ys, μs)
@test all(isfinite,logliksbig)

ys = fill(big"1.0", size(μs))
devresidsbig = devresid.(PositivePoisson(λ0), ys, μs)
@test all(isfinite,devresidsbig)
logliks1 = loglik.(ys, μs)
@test all(isfinite,logliks1)

# R benchmark
seed=12
nn=1000
b0=1.0
b1=-2.0
coefs0=[b0,b1]

# # uncomment code to generate R benchmark and save it to csv files
# using RCall
# R"if(!require(VGAM)){install.packages(\"VGAM\");library(VGAM)}"
# R"library(VGAM)"
# R"set.seed($seed)"
# R"pdata <- data.frame(x2 = runif(nn <- $nn))"
# R"pdata <- transform(pdata, lambda = exp($b0 + $b1 * x2))"
# R"pdata <- transform(pdata, y1 = rpospois(nn, lambda))"
# R"with(pdata, table(y1))"
# R"fit <- vglm(y1 ~ x2, pospoisson, data = pdata, trace = TRUE, crit = \"coef\")"
# print(R"summary(fit)")
# pdata=rcopy("pdata")
# coefsR=vec(rcopy(R"coef(fit, matrix = TRUE)"))
# coefsRdf = DataFrame(intercept=[coefsR[1]],x2=[coefsR[2]])
# writetable(joinpath(testdir,"data","positive_poisson_pdata.csv"),pdata)
# writetable(joinpath(testdir,"data","positive_poisson_coefsR.csv"),coefsRdf)

# load saved R benchmark
import CSV
pdata=CSV.read(joinpath(testdir,"data","positive_poisson_pdata.csv"))
coefsR=vec(convert(Matrix{Float64},CSV.read(joinpath(testdir,"data","positive_poisson_coefsR.csv"))))

X=convert(Array{Float64,2}, pdata[!,[:x2]])
Xwconst=[ones(size(X,1)) X]
y=convert(Array{Float64,1},pdata.y1)

glmfit = fit(GeneralizedLinearModel,Xwconst,y,PositivePoisson(),LogProductLogLink())

# glmfit = fit(GeneralizedLinearModel,Xwconst,y,PositivePoisson(),LogProductLogLink();convTol=1e-2)
coefsGLM = coef(glmfit)
@test coefsGLM ≈ coefsR rtol=1e-7
@test coefsGLM ≈ coefs0 rtol=0.05

# GammaLassoPath without actualy regularization
glpfit = fit(GammaLassoPath,X,y,PositivePoisson(),LogProductLogLink(); λ=[0.0])
coefsGLP = vec(coef(glpfit))
@test coefsGLP ≈ coefsGLM

# GammaLassoPath doing Lasso
lassofit = fit(GammaLassoPath,X,y,PositivePoisson(),LogProductLogLink(); γ=0.0)
coefsLasso = vec(coef(lassofit;select=MinAICc()))
@test coefsLasso ≈ coefsGLM rtol=0.02
@test coefsLasso ≈ coefs0 rtol=0.05

# GammaLassoPath doing concave regularization
glpfit = fit(GammaLassoPath,X,y,PositivePoisson(),LogProductLogLink(); γ=10.0)
coefsGLP = vec(coef(glpfit;select=MinAICc()))
@test coefsGLP ≈ coefsLasso rtol=0.0002
@test coefsGLP ≈ coefs0 rtol=0.05

# # problematic case 1
# X = [6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 6.693264063357201 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 5.354611250685761 4.015958438014321 4.015958438014321 4.015958438014321 4.015958438014321 4.015958438014321 4.015958438014321 4.015958438014321 4.015958438014321 4.015958438014321 2.6773056253428806 2.6773056253428806 2.6773056253428806 2.6773056253428806 1.3386528126714403 1.3386528126714403 1.3386528126714403 1.3386528126714403]'
# y = [1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]
# yalt = [3.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]
# y = [1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]
# y=ones(size(X,1))
# y[end]=1
# problem1fit = fit(LassoPath,X,y,Binomial();verbose=true,naivealgorithm=true)
# problem1fit = fit(GeneralizedLinearModel,X,y,Poisson();verbose=true)
#
# problem1fit = fit(GammaLassoPath,X,y,Poisson();γ=0.0,verbose=true,naivealgorithm=false)
# problem1fit.m.pp
# m=[]
# problem1fit=[]
# (obj,objold,path,f,newcoef,oldcoef,b0,oldb0,b0diff,coefdiff,scratchmu,cd,r,α,curλ)=(0,0,[],0,[],[],0,0,0,[],[],[],[],0,0)
# try
#   problem1fit = fit(GammaLassoPath,X,y,PositivePoisson();γ=0.0,verbose=true,naivealgorithm=false)
# catch e
#   if typeof(e) == Lasso.ConvergenceException
#     (obj,objold,path,f,newcoef,oldcoef,b0,oldb0,b0diff,coefdiff,scratchmu,cd,r,α,curλ) = e.debugvars
#     m = path.m
#     error(e.msg)
#     nothing
#   else
#     e
#   end
# end
# coef(problem1fit)
# newcoef==oldcoef
# b0diff
# coefdiff
# newcoef
# oldcoef
# newcoef-oldcoef
# scratchmu
# m.rr.devresid
# sum(m.rr.devresid)
# println(m.rr.var)
# cd
# m.rr.mu
# m.pp.oldy
# deviance(m)
# fieldnames(path.m.pp)
# curλ
# b0=-3.67186085318965
# newcoef.coef[1]=-0.08566741748972534
#
# rnew = deepcopy(m.rr)
# pnew = deepcopy(m.pp)
# newmu = Lasso.linpred!(scratchmu, pnew, newcoef, b0)
# updateμ!(rnew, newmu)
# devnew = deviance(rnew)
# curλ*Lasso.P(α, newcoef, pnew.ω)
# objnew = devnew/2 + curλ*Lasso.P(α, newcoef, pnew.ω)
#
# rold = deepcopy(m.rr)
# pold = deepcopy(m.pp)
# oldmu = Lasso.linpred!(scratchmu, pold, oldcoef, oldb0)
# updateμ!(rold, oldmu)
# devold = deviance(rold)
# objold = devold/2 + curλ*Lasso.P(α, oldcoef, pold.ω)
#
# devnew==devold
# curλ*Lasso.P(α, oldcoef, pold.ω) > curλ*Lasso.P(α, newcoef, pnew.ω)
# objnew > objold
# # step-halving failure zoomin
# function devf(b0,newcoef,m,cd,f)
#   T = eltype(y)
#   m=deepcopy(m)
#   newcoef=deepcopy(newcoef)
#   b0=deepcopy(b0)
#   cd=deepcopy(cd)
#   r=m.rr
#   p=m.pp
#   for icoef = 1:nnz(newcoef)
#       oldcoefval = icoef > nnz(oldcoef) ? zero(T) : oldcoef.coef[icoef]
#       newcoef.coef[icoef] = oldcoefval+f*(coefdiff.coef[icoef])
#   end
#   b0 = oldb0+f*b0diff
#   updateμ!(r, Lasso.linpred!(scratchmu, cd, newcoef, b0))
#   dev = deviance(r)
#   curλ*Lasso.P(α, newcoef, cd.ω)
#   obj = dev/2 + curλ*Lasso.P(α, newcoef, cd.ω)
# end
# using Gadfly
# Δf=f*1
# fs=-100f:Δf:100f
# devfs=broadcast(f->devf(b0,newcoef,m,cd,f),fs)
# plot(layer(x=fs,y=devfs,Geom.line),layer(x=[0],y=[objold],Geom.point))
# r.var
# #
#
# fieldnames(newcoef)
# size(newcoef)
#
# # problematic case 2
#
# using Lasso, StatsBase
# X=[3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 3.461994189597978 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.7695953516783822 2.0771965137587864 2.0771965137587864 2.0771965137587864 2.0771965137587864 2.0771965137587864 2.0771965137587864 2.0771965137587864 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 1.3847976758391911 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956 0.6923988379195956]'
# Xwconst=[ones(X) X]
# y=[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,1.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]
# sum(y.!=1.0)
# summarystats(y)
# std(y)
# problem2fit = fit(GeneralizedLinearModel,Xwconst,y,Poisson();verbose=true)
# problem2fit = fit(GeneralizedLinearModel,ones(X),y,Poisson();verbose=true)
# problem2fit = fit(LassoPath,X,y,Poisson();verbose=true,λ=[0.0])
# problem2fit = fit(LassoPath,X,y,Poisson();verbose=true,λ=[0.0],dofit=false)
# problem2fit = fit(LassoPath,X,y,Poisson();verbose=true)
# intercept=true
# n=length(y)
# T=eltype(X)
# d=Poisson()
# l=canonicallink(d)
# wts=ones(T, length(y))
# wts .*= convert(T, 1/sum(wts))
# offset=similar(y, 0)
# irls_tol =1e-7
# nullmodel = fit(GeneralizedLinearModel, ones(T, n, ifelse(intercept, 1, 0)), y, d, l;
#                     wts=wts, offset=offset, convTol=irls_tol)
# problem2fit = fit(GeneralizedLinearModel,ones(X),y,Poisson();verbose=true,convTol=irls_tol)
# problem2fit = fit(GeneralizedLinearModel,ones(X),y,Poisson();verbose=true,convTol=irls_tol,dofit=false)
#
# problem2fit = fit(GammaLassoPath,X,y,Poisson();γ=0.0,verbose=true)
# problem2fit = fit(LassoPath,ones(X),y,Poisson();verbose=true,intercept=false)

end


text:
export SimulationDraws

# functions to define simulations

# convert θρ ∈ R ↦ ρ ∈ (0,1)
@inline _ρ(θρ) = logistic(θρ)
@inline _dρdθρ(θρ) = (z = logistic(θρ); z*(1-z) )
@inline _ρsq(θρ) = _ρ(θρ)^2
@deprecate _dρdσ(θρ) _dρdθρ(θρ)
# @deprecate _ρ2(θρ) _ρsq(θρ)

# go from iid normal u,v to correlated shocks
@inline _ψ1(u::Number, v::Number, ρ::Number) = ρ*u + sqrt(1-ρ^2)*v
@inline _ψ2(u::Number, v::Number, ρ::Number) = _ψ2(u,v)
@inline _ψ2(u::Number, v::Number) = u

# derivative of ψ1 wrt ρ and θρ
@inline _dψ1dρ( u::Number, v::Number, ρ::Number) = u - ρ/sqrt(1-ρ^2)*v
@inline _dψ1dθρ(u::Number, v::Number, ρ::Number, θρ::Number) = _dψ1dρ(u,v,ρ)*_dρdθρ(θρ)

# z(ψ2|ψ1)
@inline cond_z(x2::Number, x1::Number, Δ::Number, ρ::Number) = (x2 - ρ*x1 + Δ)/sqrt(1-ρ^2)
# @deprecate  _z(x2::Number, x1::Number, Δ::Number, ρ::Number) conditional_z(x2, x1, Δ, ρ)

# derivatives
@inline dcond_zdρ(x2::Number, x1::Number, ρ::Number, z::Number) = -x1/sqrt(1-ρ^2) + ρ*z/(1-ρ^2)
# @deprecate  _dzdρ(x2::Number, x1::Number, ρ::Number, z::Number) dcond_zdρ(x2,x1,ρ,z)

@inline _dcond_probdρ(x2::Number, x1::Number, Δ::Number, ρ::Number) = (z = cond_z(x2,x1,Δ,ρ);  normpdf(z) * dcond_zdρ(x2,x1,ρ,z))
# @deprecate      _dπdρ(x2::Number, x1::Number, Δ::Number, ρ::Number) _dcond_probdρ(x2,x1,Δ,ρ)

# finite difference versions
@inline _ρfd(x::Number, h::Number) = _ρ(x+h)
@inline cond_zfd(x2::Number, x1::Number, Δ::Number, ρ::Number, h::Number) = cond_z(x2, x1+h, Δ, ρ)
# @deprecate    _z(x2::Number, x1::Number, Δ::Number, ρ::Number, h::Number) cond_zfd(x2,x1,Δ,ρ,h)

# data strucutres
#---------------------------

struct SimulationDraw{T,V<:AbstractVector{T}}
    psi1::T
    u::T
    dpsidrho::T
    drillgradm::V
end

eltype(::SimulationDraw{T}) where {T} = T

function SimulationDraw(u::Real, v::Real, θρ::T) where {T}
    ρ = _ρ(θρ)
    psi1 = _ψ1(u,v,ρ)
    @assert u == _ψ2(u,v,ρ)
    dpsidrho = _dψ1dθρ(u,v,ρ,θρ)
    zerovec = Vector{T}(undef,0)
    return SimulationDraw(psi1, u, dpsidrho, zerovec)
end

function SimulationDrawFD(s::SimulationDraw, x::Real)
    SimulationDraw(_ψ1(s)+x, _u(s)+x, _dψ1dθρ(s), _drillgradm(s))
end

SimulationDraw(θρ) = SimulationDraw(randn(), randn(), θρ)

struct SimulationDraws{T, N, A<:AbstractRealArray{T,N}}
    u::A
    v::A
    psi1::A
    dpsidrho::A
    qm::Vector{T} # holds log L_{im}
    am::Vector{T} # pdf(η₁) / cdf(η₁)
    bm::Vector{T} # pdf(η₂) / cdf(η₂)
    cm::Vector{T} # (pdf(η₂) - pdf(η₁)) / (cdf(η₂) - cdf(η₁))
    drillgradm::Matrix{T}

    function SimulationDraws(
        u::A, v::A, psi1::A, dpsidrho::A, qm::V, am::V, bm::V, cm::V,
        drillgradm::Matrix{T}
    ) where {T,N,A<:AbstractRealArray{T,N},V<:Vector{T}}
        size(dpsidrho) == size(psi1) == size(u) == size(v)  || throw(DimensionMismatch())
        length(qm) == length(am) == length(bm) == length(cm) == size(u,1) || throw(DimensionMismatch())
        length(qm) == size(drillgradm,2) || throw(DimensionMismatch())
        return new{T,N,A}(u,v,psi1,dpsidrho,qm,am,bm,cm,drillgradm)
    end
end

const SimulationDrawsVector{T,A} = SimulationDraws{T,1,A} where {T,A}
const SimulationDrawsMatrix{T,A} = SimulationDraws{T,2,A} where {T,A}
const SimulationDrawOrDraws = Union{<:SimulationDraws,<:SimulationDraw}

# generate simulations using Halton
#---------------------------

function SimulationDraws(M::Integer, N::Integer, K::Integer=0; bases::NTuple{2,Int}=(2,3))
    uvs = map((b) -> Matrix{Float64}(undef, M,N), bases)
    HaltonDraws!.(uvs, bases; skip=5000, distr=Normal())
    psi1 = similar(uvs[1])
    dpsidrho = similar(uvs[1])
    qm = Vector{Float64}(undef, M)
    am = similar(qm)
    bm = similar(qm)
    cm = similar(qm)
    drillgradm = Matrix{Float64}(undef,K,M)
    return SimulationDraws(uvs..., psi1, dpsidrho, qm, am, bm, cm, drillgradm)
end

function SimulationDraws(M, data::AbstractDataSet; kwargs...)
    return SimulationDraws(M, length(data), _nparm(_model(data)); kwargs...)
end

function SimulationDraws(M, data::AbstractDataSetofSets; kwargs...)
    drillmod = _model(drill(data))
    return SimulationDraws(M, num_i(data), _nparm(drillmod); kwargs...)
end

# interface
#---------------------------

# access fields
_v(    s::SimulationDraws) = s.v
_u(    s::SimulationDrawOrDraws) = s.u
_ψ1(   s::SimulationDrawOrDraws) = s.psi1
_dψ1dθρ(s::SimulationDrawOrDraws) = s.dpsidrho
_qm(   s::SimulationDrawOrDraws) = s.qm
_drillgradm(s::SimulationDrawOrDraws) = s.drillgradm
_am(   s::SimulationDraws) = s.am
_bm(   s::SimulationDraws) = s.bm
_cm(   s::SimulationDraws) = s.cm

_ψ2(   s::SimulationDrawOrDraws) = _u(s)
_psi1( s::SimulationDrawOrDraws) = _ψ1(s)
_psi2( s::SimulationDrawOrDraws) = _ψ2(s)
_llm(  s::SimulationDrawOrDraws) = _qm(s)

# @deprecate _LLm(s::SimulationDrawOrDraws) _llm(s)
# @deprecate _dψ1dρ(s::SimulationDrawOrDraws) _dψ1dθρ(s)

_num_sim(s::SimulationDraws) = size(_u(s),1)
num_i(s::SimulationDrawsMatrix) = size(_u(s),2)
num_i(s::SimulationDrawsVector) = 1

_num_sim(M::Integer) = M

# manipulate like array
tup(s::SimulationDraws) = _u(s), _v(s), _ψ1(s), _dψ1dθρ(s)
size(s::SimulationDraws) = size(_u(s))
_nparm(s::SimulationDraws) = size(_drillgradm(s),1)

function view(s::SimulationDrawsMatrix, i::Integer)
    psiviews = view.(tup(s), :, i)
    SimulationDraws(psiviews..., _qm(s), _am(s), _bm(s), _cm(s), _drillgradm(s))
end

function getindex(s::SimulationDrawsVector, m)
    u, v, ψ1, dψ1dρ = getindex.(tup(s), m)
    gradmvw = view(_drillgradm(s), :, m)
    return SimulationDraw(ψ1, u, dψ1dρ, gradmvw)
end

# doing mean(psi) and mean(psi^2)
function psi2_wtd_sum_and_sumsq(s::SimulationDrawsVector{T}) where {T}
    qm = _qm(s)
    ψ = _ψ2(s)
    ψbar = dot(qm,ψ)
    ψ2bar = sumprod3(qm,ψ,ψ)
    return ψbar, ψ2bar
end

# so we can pre-calculate shocks
function update_ψ1!(ψ1::AA, u::AA, v::AA, θρ::Real) where {AA<:AbstractArray}
    size(ψ1) == size(u) == size(v) || throw(DimensionMismatch())
    ρ = _ρ(θρ)
    0 <= ρ <= 1 || @warn "0 <= ρ=$ρ <= 1 is false. Have θρ = $θρ"
    ψ1 .= _ψ1.(u, v, ρ)
    return nothing
end

function update_dψ1dθρ!(dψ1dθρ::A, u::A, v::A, θρ::Real) where {A<:AbstractArray}
    size(dψ1dθρ) == size(u) == size(v) || throw(DimensionMismatch())
    ρ = _ρ(θρ)
    0 <= ρ <= 1 || @warn "0 <= ρ=$ρ <= 1 is false. Have θρ = $θρ"
    dψ1dθρ .= _dψ1dθρ.(u,v,ρ,θρ)
    return nothing
end

update_ψ1!(    s::SimulationDraws, θρ) = update_ψ1!(     _ψ1(s),    _u(s), _v(s), θρ)
update_dψ1dθρ!(s::SimulationDraws, θρ) = update_dψ1dθρ!(_dψ1dθρ(s), _u(s), _v(s), θρ)

function update!(sim::SimulationDraws,θρ::Real)
    update_ψ1!(sim, θρ)
    update_dψ1dθρ!(sim, θρ)
end


text:
using XLSX, CSV, DataFrames

function main()
  faltantes = CSV.read("faltantes.csv", DataFrame)
  n = size(faltantes, 1)
  faltantes = [faltantes DataFrame(horas = zeros(Int, n))]
  xlsx = XLSX.readxlsx("presença.xlsx")
  pesos = [3, 2, 3, 2, 4, 2, 4, 2, 2, 4]
  for i in 1:XLSX.sheetcount(xlsx)
    sh = xlsx[i]
    emails = unique(sort(lowercase.(sh["B"][2:end])))
    I = indexin(emails, faltantes[!, "email"])
    I = I[findall(I .!== nothing)]
    faltantes[I, "horas"] .+= pesos[i]
  end
  CSV.write("faltantes-horas.csv", faltantes)
end

main()

max_stars_repo_name: CiDAMO/scripts-certificados

max_issues_repo_name: CiDAMO/scripts-certificados

max_forks_repo_name: CiDAMO/scripts-certificados


text:
@testset "TimeSimulation" begin
      # Time response from a single particle
      include("../docs/src/example/time_response_single_particle/time_response_single_particle.jl")
      freq_result, time_result = run_time_response_single_particle()
      # Spike at start and then two diffracted spikes later, then
      # not much else

      # incident plane wave
      @test abs( maximum(field(time_result)) - 1.0) < 1e-2

      # first diffraction
      i1 = findmin(abs.(time_result.t .- 198.5))[2]
      m, i = findmax(field(time_result)[i1-50:i1])
      t = time_result.t[i1-50:i1][i]
      @test abs(m - 0.046) < 2e-2
      @test abs(t - 198) < 1e-1

      # second diffraction
      i2 = findmin(abs.(time_result.t .- 210.0))[2]
      m, i = findmax(field(time_result)[i1+5:i2])
      t = time_result.t[i1+5:i2][i]
      @test abs(m - 0.05) < 2e-2
      @test abs(t - 206) < 2e-1

end


text:
export getvars, getvar, cutdata, subvolume, subsurface

"""
	cutdata(data, var; plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1)

Get 2D plane cut in orientation `dir` for `var` out of 3D box `data` within `plotrange`.
The returned 2D data lies in the `sequence` plane from - to + in `dir`.
"""
function cutdata(data::Data, var::AbstractString;
   plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1)

   x, w = data.x, data.w
   VarIndex_ = findfirst(x->x==lowercase(var), lowercase.(data.head.wnames))
   isempty(VarIndex_) && error("$(var) not found in header variables!")

   X = @view x[:,:,:,1]
   Y = @view x[:,:,:,2]
   Z = @view x[:,:,:,3]

   W = w[:,:,:,VarIndex_]

   if dir == "x"
      cut1 = @view X[sequence,:,:]
      cut2 = @view Y[sequence,:,:]
      W    = @view W[sequence,:,:]
   elseif dir == "y"
      cut1 = @view X[:,sequence,:]
      cut2 = @view Z[:,sequence,:]
      W    = @view W[:,sequence,:]
   elseif dir == "z"
      cut1 = @view X[:,:,sequence]
      cut2 = @view Y[:,:,sequence]
      W    = @view W[:,:,sequence]
   end

   if !all(isinf.(plotrange))
      cut1, cut2, W = subsurface(cut1, cut2, W, plotrange)
   end

   cut1, cut2, W
end

@inline function checkvalidlimits(limits, dim=2)
   if dim == 2
      if length(limits) != 4
         throw(ArgumentError("Reduction range $limits should be [xmin xmax ymin ymax]!"))
      end

      if limits[1] > limits[2] || limits[3] > limits[4]
         throw(DomainError(limits, "Invalid reduction range!"))
      end
   elseif dim == 3
      if length(limits) != 6
         throw(ArgumentError(
            "Reduction range $limits should be [xmin xmax ymin ymax zmin max]!"))
      end

      if limits[1] > limits[2] || limits[3] > limits[4] || limits[5] > limits[6]
         throw(DomainError(limits, "Invalid reduction range!"))
      end
   end
end

@inline function findindexes(x, y, limits)
   hx = @view x[:,1]
   hy = @view y[1,:]

   if isinf(limits[1]) limits[1] = minimum(hx) end
   if isinf(limits[3]) limits[3] = minimum(hy) end
   if isinf(limits[2]) limits[2] = maximum(hx) end
   if isinf(limits[4]) limits[4] = maximum(hy) end

   xind = findall(limits[1] .≤ hx .≤ limits[2])
   yind = findall(limits[3] .≤ hy .≤ limits[4])
   xind, yind
end

@inline function findindexes(x, y, z, limits)
   hx = @view x[:,1,1]
   hy = @view y[1,:,1]
   hz = @view z[1,1,:]

   if isinf(limits[1]) limits[1] = minimum(hx) end
   if isinf(limits[3]) limits[3] = minimum(hy) end
   if isinf(limits[5]) limits[5] = minimum(hz) end
   if isinf(limits[2]) limits[2] = maximum(hx) end
   if isinf(limits[4]) limits[4] = maximum(hy) end
   if isinf(limits[6]) limits[6] = maximum(hz) end

   xind = findall(limits[1] .≤ hx .≤ limits[2])
   yind = findall(limits[3] .≤ hy .≤ limits[4])
   zind = findall(limits[5] .≤ hy .≤ limits[6])
   xind, yind, zind
end

"""
	subsurface(x, y, data, limits)
	subsurface(x, y, u, v, limits)

Extract subset of 2D surface dataset. See also: [`subvolume`](@ref).
"""
function subsurface(x, y, data, limits)

   checkvalidlimits(limits)

   xind, yind = findindexes(x, y, limits)

   newdata = subdata(data, xind, yind, size(data))

   newx = x[xind, yind]
   newy = y[xind, yind]

   newx, newy, newdata
end

function subsurface(x, y, u, v, limits)

   checkvalidlimits(limits)

   xind, yind = findindexes(x, y, limits)

   sz = size(u)
   newu = subdata(u, xind, yind, sz)
   newv = subdata(v, xind, yind, sz)

   newx = x[xind, yind]
   newy = y[xind, yind]

   newx, newy, newu, newv
end

"""
	subvolume(x, y, z, data, limits)
	subvolume(x, y, z, u, v, w, limits)

Extract subset of 3D dataset in ndgrid format. See also: [`subsurface`](@ref).
"""
function subvolume(x, y, z, data, limits)

   checkvalidlimits(limits, 3)

   xind, yind, zind = findindexes(x, y, z, limits)

   newdata = subdata(data, xind, yind, zind, size(data))

   newx = x[xind,yind,zind]
   newy = y[xind,yind,zind]
   newz = z[xind,yind,zind]

   newx, newy, newz, newdata
end

function subvolume(x, y, z, u, v, w, limits)

   checkvalidlimits(limits, 3)

   sz = size(u)

   xind, yind, zind = findindexes(x, y, z, limits)

   newu = subdata(u, xind, yind, zind, sz)
   newv = subdata(v, xind, yind, zind, sz)
   neww = subdata(w, xind, yind, zind, sz)

   newx = x[xind,yind,zind]
   newy = y[xind,yind,zind]
   newz = z[xind,yind,zind]

   newx, newy, newz, newu, newv, neww
end

"""
    subdata(data, xind, yind, sz)
    subdata(data, xind, yind, zind, sz)

Return the sliced data based on indexes `xind` and `yind` of size `sz`.
"""
function subdata(data, xind, yind, sz)
   newdata = data[xind,yind]
   newsz = size(newdata)

   if length(sz) > 2
      newdata = reshape(newdata, (newsz[1:2]..., sz[3:end]))
   end

   newdata
end

function subdata(data, xind, yind, zind, sz)

   newdata = data[xind,yind,zind]
   newsz = size(newdata)

   if length(sz) > 3
      newdata = reshape(newdata, (newsz[1:3]..., sz[4:end]))
   end

   newdata
end

"Return variable data from string `var`."
function getvar(data::Data, var)
   VarIndex_ = findfirst(x->x==lowercase(var), lowercase.(data.head.wnames))
   isnothing(VarIndex_) && error("$(var) not found in file header variables!")

   ndim = data.head.ndim
   if ndim == 1
      w = data.w[:,VarIndex_]
   elseif ndim == 2
      w = data.w[:,:,VarIndex_]
   elseif ndim == 3
      w = data.w[:,:,:,VarIndex_]
   end
   w
end

"""
    getvars(data::Data, Names::Vector) -> Dict

Return variables' data as a dictionary from string vector.
See also: [`getvar`](@ref).
"""
function getvars(data::Data, Names::Vector{T}) where T<:AbstractString

   dict = Dict()
   for name in Names
      dict[name] = getvar(data, name)
   end

   dict
end

const variables_predefined = Dict(
   "B" => data -> sqrt.(getvar(data, "Bx").^2 .+ getvar(data, "By").^2 .+ getvar(data, "Bz").^2),
   "E" => data -> sqrt.(getvar(data, "Ex").^2 .+ getvar(data, "Ey").^2 .+ getvar(data, "Ez").^2),
   "U" => data -> sqrt.(getvar(data, "Ux").^2 .+ getvar(data, "Uy").^2 .+ getvar(data, "Uz").^2),
   #"beta" => data -> getvar(data, "P") ./ getvar(data, "B").^2 * 2μ,
)


text:
__precompile__(true)
"""
Main module for `Psychometrics.jl` -- A Julia package that provides tools for psychometric data analysis.

# Exports

    TruncatedInverseGaussian
    PolyaGamma
    TruncatedGamma
    rand
    AbstractExaminee
    Examinee
    AbstractLatent
    Latent1D
    LatentND
    AbstractItem
    Item
    AbstractParameters
    AbstractParametersBinary
    Parameters1PL
    Parameters2PL
    Parameters3PL
    ParametersNPL
    AbstractResponse
    Response
    probability
    log_likelihood
    likelihood
    latent_information
    item_observed_information
    item_expected_information
    add_prior!
    get_item_by_id
    get_parameters
    empty_chain!
    add_response!
    answer
    get_responses_by_item_id
    get_responses_by_item_idx
    get_responses_by_examinee_id
    get_response_matrix
    get_design_matrix
    get_responses
    get_examinee_by_id
    get_items_idx_answered_by_examinee
    get_examinees_idx_who_answered_item
    generate_w
    set_val_from_chain!
    set_val_from_posterior!
    chain_append!
    mcmc_iter_pg!
    update_estimate!
    posterior
    update_posterior!
    get_latents_vals
    get_latents
    find_best_item
    find_best_examinee
    calibrate_item!
    assess!
    joint_estimation!
"""
module Psychometrics
import Distributed
   using NLopt

    import Distributions
    import SharedArrays
    import Dates
    import LinearAlgebra
    import Interpolations
    import Random
    import SpecialFunctions
    import Base.copy
    using Requires

    include("utils/math/utils_math.jl")
    include("distributions/distributions.jl")
    include("examinee/examinee.jl")
    include("examinee/examinees.jl")
    include("item/item.jl")
    include("item/items.jl")
    include("response.jl")
    include("probability.jl")
    include("likelihood.jl")
    include("utils/polyagamma/utils_pg.jl")
    include("posterior/posterior.jl")
    include("information/information.jl")
    include("online/online.jl")
    include("utils/mmle/utils_mmle.jl")
    include("item/calibration/calibration.jl")
    include("examinee/assessment/assessment.jl")
    include("joint_estimation/joint_estimation.jl")
    include("bootstrap/bootstrap.jl")

export
    TruncatedInverseGaussian,
    TruncatedNormal,
    PolyaGamma,
    TruncatedGamma,
    rand,
    AbstractExaminee,
    Examinee,
    AbstractLatent,
    Latent1D,
    LatentND,
    AbstractItem,
    Item,
    AbstractParameters,
    AbstractParametersBinary,
    Parameters1PL,
    Parameters2PL,
    Parameters3PL,
    ParametersNPL,
    AbstractResponse,
    Response,
    probability,
    log_likelihood,
    likelihood,
    latent_information,
    item_observed_information,
    item_expected_information,
    latent_information_3PL,
    add_prior!,
    get_item_by_id,
    get_parameters,
    empty_chain!,
    add_response!,
    answer,
    get_responses_by_item_id,
    get_responses_by_item_idx,
    get_responses_by_examinee_id,
    get_response_matrix,
    get_design_matrix,
    get_responses,
    get_examinee_by_id,
    get_items_idx_answered_by_examinee,
    get_examinees_idx_who_answered_item,
    generate_w,
    set_val_from_chain!,
    set_val_from_posterior!,
    chain_append!,
    get_latents,
    find_best_item,
    find_best_examinee,
    calibrate_item!,
    assess_examinee!,
    joint_estimate!,
    rescale!,
    bootstrap!
end # module


text:
# buidling connection matrices

# Utility functions
# for connectivity

"""
        lognorm_reparametrize(m,std) -> distr::LogNormal
# parameters
  + `m`   sample mean
  + `std` sample std
"""
function lognorm_reparametrize(m::Real,std::Real)
    vm2= (std/m)^2
    μ = log(m / sqrt(1. + vm2))
    σ = sqrt(log( 1. + vm2))
    return LogNormal(μ,σ)
end

function no_autapses!(mat::SparseMatrixCSC)
  for i in 1: min(size(mat)...)
    mat[i,i]=0.0
  end
  return dropzeros!(mat)
end

# scale rows so that the sum of the nonzero elements is exactly x
# for simplicity ignores the sign of x
function scale_rows!(x::Real,mat::SparseMatrixCSC)
  x=abs(x) # there are better ways to change sign!
  n = size(mat,1)
  vals=nonzeros(mat)
  rows=rowvals(mat)
  for i in 1:n
    idxs = findall(r->r==i,rows)
    if !isempty(idxs) # is it control useful? or less efficient?
      s = sum(vals[idxs])
      vals[idxs] .*= x/s
    end
  end
  return mat
end

function _sparse_wmat(npost::Integer,npre::Integer,p::Real,make_weights!::Function ;
    scal::Float64=1.0, # final scaling (e.g. -1 for inhibitory neurons)
    noself::Bool=true, # no autapses ?
    rowsum::Union{Nothing,Float64}=nothing) #mean of each row (before scal)
  wmat = sprand(Float64,npost,npre,p)
  if noself
    no_autapses!(wmat)
  end
  make_weights!(wmat)
  if !isnothing(rowsum) # scale so that each row has mean exactly rowsum mean
    scale_rows!(rowsum,wmat)
  end
  if scal != 1.0
    lmul!(scal,wmat)
  end
  return wmat
end

function sparse_wmat_rowfun(make_a_row::Function,
    npost::Integer,npre::Integer,p::Real;no_autapses::Bool=false)
  wmat = sprand(Float64,npost,npre,p)
  if no_autapses
    no_autapses!(wmat)
  end
  rw = rowvals(wmat)
  wvec = nonzeros(wmat)
  for r in unique(rw)
    ridx = r .== rw
    n = count(ridx)
    newrow = make_a_row(n)
    wvec[ridx] .= newrow
  end
  return wmat
end

function sparse_constant_wmat(npost::Integer,npre::Integer,p::Real,j_val::Real ;
    scal::Float64=1.0, # final scaling
    no_autapses::Bool=true, # no autapses ?
    rowsum::Union{Nothing,Float64}=nothing)
  make_weights=function(n::Integer)
    ret=fill(j_val,n)
    if !isnothing(rowsum)
      ret .*= rowsum / sum(ret)
    end
    if scal != 1
      lmul!(scal,ret)
    end
    return ret
  end
  return sparse_wmat_rowfun(make_weights,npost,npre,p;no_autapses=no_autapses)
end

function sparse_wmat_lognorm(npost::Integer,npre::Integer,
    p::Real,μ::Real,σ::Real;
    noself::Bool=true, # no autapses ?
    rowsum::Union{Nothing,Float64}=nothing)
  μp = abs(μ)
  distr = lognorm_reparametrize(μp,σ)
  make_weights=function (smat)
    vals=nonzeros(smat)
    Distributions.rand!(distr,vals)
  end
  scal = μ < 0 ? -1.0 : 1.0
  return _sparse_wmat(npost,npre,p,make_weights;
    scal=scal,noself=noself,rowsum=rowsum)
end

function make_sparse_weights(npost::Integer,npre::Integer,p::Real,weigth_value::Real ;
    scal::Float64=1.0, # final scaling
    noself::Bool=true, # no autapses ?
    rowsum::Union{Nothing,Float64}=nothing)
  make_weights=function (smat)
    vals=nonzeros(smat)
    fill!(vals,weigth_value)
  end
  return _sparse_wmat(npost,npre,p,make_weights;
    scal=scal,noself=noself,rowsum=rowsum)
end

# not needed with new Poisson input types!
# """
#     wmat_train_assemblies_protocol(Nneus::Integer,p::Real; scal::Float64=1.0)
#       -> (weights::SparseMatrixCSC,assembly_dictionary::Dict)

#   # Arguments
#   + Npost::Integer : number of postsynaptic neurons
#   + Nassemblies::Integer : number of assemblies to train
#   + p_assembly::Real : probability for a neuron to be in any assembly
#   + scal::Real=1.0 : connection weights

#   # Outputs
#   + weights : sparse weight matrix
#   + assembly_list::Vector{NamedTuple{(:neupost,:neupre),Vector{Int64,Vector{Int64}}}}
#     element `i` contains the pre and post neurons associated to assembly `i`
# """
# function wmat_train_assemblies_protocol(Npost::Integer,Nassemblies::Integer,
#     p_assembly::Real;scal::Float64=1.0)
#   # unoptimized... but called only once
#   post_neus = [findall(rand(Npost) .< p_assembly) for _ in 1:Nassemblies]
#   c = 0
#   pre_neus = map(post_neus) do _po
#     _ln = length(_po)
#     ret = c .+ collect(1:_ln)
#     c += _ln
#     return ret
#   end
#   _pre_neus_all = vcat(pre_neus...)
#   wmat = sparse(vcat(post_neus...),_pre_neus_all,scal,
#     Npost,length(_pre_neus_all)) # include total matrix size
#   assembly_list = map(zip(post_neus,pre_neus)) do (_post,_pre)
#     (neupost=_post,neupre=_pre)
#   end
#   return wmat,assembly_list
# end

# to order assemblies (no idea where to place this function)
function order_by_pattern_idxs(pattern_idxs::Vector{Vector{Int64}},
    nneus::Int64;rev::Bool=false)
  _tosort = map(1:nneus) do neu
    pa = findall(patt->neu in patt,pattern_idxs)
    isnothing(pa) && return 10_000_000
    (length(pa) == 1) && return pa[1]
    (length(pa)>1) && return sample(pa)
  end
  return sortperm(_tosort;rev=rev)
end


text:
# ---
# title: 782. Transform to Chessboard
# id: problem782
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Math
# link: <https://leetcode.com/problems/transform-to-chessboard/description/>
# hidden: true
# ---
#
# An N x N `board` contains only `0`s and `1`s. In each move, you can swap any 2
# rows with each other, or any 2 columns with each other.
#
# What is the minimum number of moves to transform the board into a "chessboard"
# \- a board where no `0`s and no `1`s are 4-directionally adjacent? If the task
# is impossible, return -1.
#
#
#
#     **Examples:**
#     Input: board = [[0,1,1,0],[0,1,1,0],[1,0,0,1],[1,0,0,1]]
#     Output: 2
#     Explanation:
#     One potential sequence of moves is shown below, from left to right:
#
#     0110     1010     1010
#     0110 --> 1010 --> 0101
#     1001     0101     1010
#     1001     0101     0101
#
#     The first move swaps the first and second column.
#     The second move swaps the second and third row.
#
#
#     Input: board = [[0, 1], [1, 0]]
#     Output: 0
#     Explanation:
#     Also note that the board with 0 in the top left corner,
#     01
#     10
#
#     is also a valid chessboard.
#
#     Input: board = [[1, 0], [1, 0]]
#     Output: -1
#     Explanation:
#     No matter what sequence of moves you make, you cannot end with a valid chessboard.
#
#
# **Note:**
#
#   * `board` will have the same number of rows and columns, a number in the range `[2, 30]`.
#   * `board[i][j]` will be only `0`s or `1`s.
#
#
## @lc code=start
using LeetCode

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


text:
type SGDRegressor
    weights::Array{Float64,2}
    bias::Float64
    params::Dict
end

function fit(::Type{SGDRegressor}, X::Matrix{Float64}, y::Array,
             shuffle=true, lr=0.01, n_iter=500)
    params = Dict("shuffle"=>shuffle, "lr"=>lr, "n_iter"=>n_iter)
    lm = SGDRegressor(randn(size(X)[1])', 0.0, params)
    inds = 1:endof(y)
    if shuffle
        inds = collect(inds)
    end
    for i = 1:n_iter
        if shuffle
            shuffle!(inds)
        end
        for j = inds
            pred = predict(lm, X[:, j:j])[1]
            lm.weights -= lr* (pred - y[j]) * X[:, j]'
            lm.bias -= lr * (pred - y[j])
        end
    end
    return lm
end

function predict(lm::SGDRegressor, X::Matrix{Float64})
    return lm.weights * X + lm.bias
end

X = randn(10,1000)
w = randn(10)'
y = w * X + 5


text:
using ArgParse
using WindMixing

s = ArgParseSettings()
@add_arg_table s begin
    "--input"
        help = "location of input file in training_output directory"
        required = true
    "--output"
        help = "destination of output file in extracted_training_output directory"
        required = true
    "--type"
        help = "training type: NDE or NN"
        default = "NDE"
end

arg_parse = parse_args(s)

FILE_PATH = joinpath(pwd(), "training_output", arg_parse["input"])
OUTPUT_PATH = joinpath(pwd(), "extracted_training_output", arg_parse["output"])
# OUTPUT_PATH = joinpath("D:\\University Matters\\MIT\\CLiMA Project\\OceanParameterizations.jl\\training_output", arg_parse["output"])

type = arg_parse["type"]

extract_NN(FILE_PATH, OUTPUT_PATH, type)


text:
module InventoryModels
#using DataStructures : reset! as dsreset!
using Distributions, MacroTools, Base.Iterators, Reexport

@reexport using ReinforcementLearningBase
using Requires

const NumDist = Union{Number, Distribution}
const State = Union{NumDist, AbstractVector{<:NumDist}}
#const RLBase = ReinforcementLearningBase

abstract type AbstractItem end

include("elementary blocks/inventory.jl")
export Inventory, LinearHoldingCost
include("elementary blocks/leadtime.jl")
export LeadTime
include("elementary blocks/policies.jl")
export sSPolicy, RQPolicy, QPolicy, BQPolicy
include("elementary blocks/supplier.jl")
export Supplier, LinearOrderCost, FixedLinearOrderCost
include("elementary blocks/market.jl")
export Market, LinearStockoutCost
include("elementary blocks/assembly.jl")
export Assembly
include("items/endproduct.jl")
export EndProduct
include("items/item.jl")
export Item
include("cvnormal.jl")
export CVNormal, cv
include("constraints.jl")
export RessourceConstraint
include("inventory_system.jl")
export InventorySystem, state_size, action_size, print_state, print_action
include("utils.jl")
export test_policy
include("zoo.jl")
export sl_sip
include("Scarf/Scarf.jl")
export Scarf
include("Scarf/interface.jl")
export ISLogger
include("logger.jl")

#state less objects are unchanged when reseted
reset!(::Any) = nothing

#=export dashboard, draw_graph
include("dashboard/dashboard.jl")=#
end