text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@with_kw mutable struct CombinatorialCutsData{T <: AbstractFloat}
    L::T = -1e10
    solve_relaxed::Bool = false
end

@with_kw mutable struct CombinatorialCutsParameters{T <: AbstractFloat}
    lower_bound::T = -1e10
    alternate::Bool = false
    update_L_every::Int = 0
end
"""
    CombinatorialCutsMaster

Master functor object for using weak optimality cuts in an integer L-shaped algorithm. Requires all first-stage decisions to be binary. Create by supplying a [`CombinatorialCuts`](@ref) object through `integer_strategy` in `LShaped.Optimizer` or set the [`IntegerStrategy`](@ref) attribute.

"""
struct CombinatorialCutsMaster{T <: AbstractFloat} <: AbstractIntegerAlgorithm
    parameters::CombinatorialCutsParameters{T}

    function CombinatorialCutsMaster(::Type{T}; kw...) where T <: AbstractFloat
        T_ = promote_type(T, Float32)
        return new{T_}(CombinatorialCutsParameters{T}(; kw...))
    end
end

function initialize_integer_algorithm!(master::CombinatorialCutsMaster, first_stage::JuMP.Model)
    # Sanity check
    if !all(is_binary, all_decision_variables(first_stage, 1))
        error("Combinatorial cuts require all first-stage decisions to be binary.")
    end
    return nothing
end

function handle_integrality!(lshaped::AbstractLShaped, master::CombinatorialCutsMaster)
    if count(active_model_objectives(lshaped)) != num_thetas(lshaped)
        # Only update lower bound if all master variables have been added
        return nothing
    end
    # Ensure all binary decisions are rounded
    for (i,dvar) in enumerate(all_decision_variables(lshaped.structure.first_stage, 1))
        lshaped.x[i] = round(lshaped.x[i])
    end
    if master.parameters.lower_bound <= -1e10 ||
        (master.parameters.update_L_every != 0 &&
         rem(lshaped.data.iterations, master.parameters.update_L_every) == 0)
        # Check if lower bound should be calculated
        if evaluate_first_stage(lshaped, current_decision(lshaped)) > sqrt(eps())
            # Get sense
            sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
            coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
            # Cache current objective
            F = MOI.get(lshaped.master, MOI.ObjectiveFunctionType())
            objective = MOI.get(lshaped.master, MOI.ObjectiveFunction{F}())
            # Replace objective with only model objectives
            F = MOI.ScalarAffineFunction{Float64}
            MOI.set(lshaped.master, MOI.ObjectiveFunction{F}(), zero(F))
            for var in lshaped.master_variables
                MOI.modify(lshaped.master,
                           MOI.ObjectiveFunction{F}(),
                           MOI.ScalarCoefficientChange(var, coeff))
            end
            MOI.optimize!(lshaped.master)
            θs = map(lshaped.master_variables) do vi
                if vi.value == 0
                    -1e10
                else
                    MOI.get(lshaped.master, MOI.VariablePrimal(), vi)
                end
            end
            # Re-add objective
            F = typeof(objective)
            MOI.set(lshaped.master, MOI.ObjectiveFunction{F}(), objective)
        else
            θs = model_objectives(lshaped)
        end
        master.parameters.lower_bound = sum(θs)
        mutate_subproblems!(lshaped) do subproblem
            update_lower_bound!(subproblem.integer_algorithm, θs[subproblem.id])
        end
    end
    return nothing
end
"""
    CombinatorialCutsWorker

Worker functor object for using weak optimality cuts in an integer L-shaped algorithm. Create by supplying a [`CombinatorialCuts`](@ref) object through `integer_strategy` in `LShaped.Optimizer` or set the [`IntegerStrategy`](@ref) attribute.

"""
struct CombinatorialCutsWorker{T <: AbstractFloat} <: AbstractIntegerAlgorithm
    data::CombinatorialCutsData{T}
    parameters::CombinatorialCutsParameters{T}
    integer_variables::Vector{MOI.VariableIndex}

    function CombinatorialCutsWorker(::Type{T}; kw...) where T <: AbstractFloat
        T_ = promote_type(T, Float32)
        worker = new{T_}(CombinatorialCutsData{T_}(),
                         CombinatorialCutsParameters{T}(; kw...),
                         Vector{MOI.VariableIndex}())
        worker.data.L = worker.parameters.lower_bound
        return worker
    end
end

function initialize_integer_algorithm!(worker::CombinatorialCutsWorker, subproblem::SubProblem)
    # Gather integer variables
    append!(worker.integer_variables, gather_integer_variables(subproblem))
    # Sanity check
    if isempty(worker.integer_variables)
        @warn "No integer variables in subproblem $(subproblem.id). Integer strategy is superfluous."
    end
    return nothing
end

function integer_variables(worker::CombinatorialCutsWorker)
    return worker.integer_variables
end

function update_lower_bound!(worker::CombinatorialCutsWorker, L::AbstractFloat)
    worker.data.L = L
    return nothing
end

function solve_subproblem(subproblem::SubProblem,
                          metadata,
                          ::NoFeasibilityAlgorithm,
                          worker::CombinatorialCutsWorker,
                          x::AbstractVector)
    # Check if lower bound has been set
    if worker.data.L <= -1e10 || worker.data.solve_relaxed
        # Relax integrality and compute regular optimality cuts
        unrelax = relax_decision_integrality(subproblem.model)
        cut = solve_subproblem(subproblem, x)
        # Unrelax integrality restrictions again
        unrelax()
        # Check if integer restrictions are satisfied
        set_metadata!(metadata,
                      subproblem.id,
                      :integral_solution,
                      check_integrality_restrictions(subproblem))
        # Solve integer problem next iteration
        worker.data.solve_relaxed = false
        # Return standard optimality cut (if successful)
        return cut
    end
    # Solve relaxed problem next iteration if alternating
    worker.data.solve_relaxed = worker.parameters.alternate
    # Solve subproblem with integer restrictions
    MOI.optimize!(subproblem.optimizer)
    status = MOI.get(subproblem.optimizer, MOI.TerminationStatus())
    if status ∈ AcceptableTermination
        # Integer restrictions are satisfied if optimal
        set_metadata!(metadata,
                      subproblem.id,
                      :integral_solution,
                      true)
        return CombinatorialOptimalityCut(subproblem, x, worker.data.L)
    elseif status == MOI.INFEASIBLE
        return Infeasible(subproblem)
    elseif status == MOI.DUAL_INFEASIBLE
        return Unbounded(subproblem)
    else
        error("Subproblem $(subproblem.id) was not solved properly, returned status code: $status")
    end
end

function solve_subproblem(subproblem::SubProblem,
                          metadata,
                          feasibility_algorithm::FeasibilityCutsWorker,
                          worker::CombinatorialCutsWorker,
                          x::AbstractVector)
    # Prepare auxiliary problem
    model = subproblem.optimizer
    if !prepared(feasibility_algorithm)
        prepare!(model, feasibility_algorithm)
    else
        activate!(model, feasibility_algorithm)
    end
    # Relax integrality
    unrelax = relax_decision_integrality(subproblem.model)
    # Optimize auxiliary problem
    MOI.optimize!(model)
    # Sanity check that aux problem could be solved
    status = MOI.get(subproblem.optimizer, MOI.TerminationStatus())
    if !(status ∈ AcceptableTermination)
        error("Subproblem $(subproblem.id) was not solved properly during feasibility check, returned status code: $status")
    end
    # check objective
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    w = correction * MOI.get(model, MOI.ObjectiveValue())
    # Ensure correction is available in master
    set_metadata!(metadata, subproblem.id, :correction, correction)
    # Check feasibility
    if w > sqrt(eps())
        # Subproblem is infeasible, create feasibility cut
        cut = FeasibilityCut(subproblem, x)
        # Unrelax integrality restrictions again
        unrelax()
        return cut
    end
    # Unrelax integrality restrictions again
    unrelax()
    # Restore subproblem and solve as usual
    deactivate!(model, feasibility_algorithm)
    return solve_subproblem(subproblem, metadata, NoFeasibilityAlgorithm(), worker, x)
end

function CombinatorialOptimalityCut(subproblem::SubProblem{T}, x::AbstractVector, L::T) where T <: AbstractFloat
    π = subproblem.probability
    cols = collect(1:length(x))
    vals = zeros(T, length(x))
    S = 0
    for (i,val) in enumerate(x)
        if isapprox(val, 0., rtol = 1e-6)
            vals[i] = -1.
        else
            # Assume x == 1
            S += 1
            vals[i] = 1.
        end
    end
    # Get sense
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    # Get sense-corrected optimal value
    Q = correction * π * MOI.get(subproblem.optimizer, MOI.ObjectiveValue())
    G = -(Q - L) * sparsevec(cols, vals, length(x))
    q = -(Q - L) * (S - 1) + L
    return OptimalityCut(G, q, subproblem.id)
end

# API
# ------------------------------------------------------------
"""
    CombinatorialCuts

Factory object for [`CombinatorialCuts`](@ref). Pass to `integer_strategy` in `LShaped.Optimizer` or set the [`IntegerStrategy`](@ref) attribute.

...
# Parameters
- `lower_bound::AbstractFloat = -1e10`: Set a lower bound on the second-stage objective, removing the need to approximate it.
- `alternate::Bool = false`: Specify if algorithm should alternate between solving relaxed problems (generating regular optimality cuts) and unrelaxed problems (generating combinatorial cuts)
- `update_L_every::Integer = 0`: Set the frequency at which the lower bound approximation should be updated. Only approximate once if set to zero.
- `optimizer = nothing`: Optionally specify an optimizer used to solve auxilliary problems in the `LiftAndProject` or `CuttingPlaneTree` strategies.
...
"""
struct CombinatorialCuts <: AbstractIntegerStrategy
    parameters::CombinatorialCutsParameters{Float64}
end
CombinatorialCuts(; kw...) = CombinatorialCuts(CombinatorialCutsParameters(; kw...))

function master(wc::CombinatorialCuts, ::Type{T}) where T <: AbstractFloat
    return CombinatorialCutsMaster(T; type2dict(wc.parameters)...)
end

function worker(wc::CombinatorialCuts, ::Type{T}) where T <: AbstractFloat
    return CombinatorialCutsWorker(T; type2dict(wc.parameters)...)
end
function worker_type(::CombinatorialCuts)
    return CombinatorialCutsWorker{Float64}
end


text:
using Plots, CSV, DataFrames, Jedi

plotlyjs()

# set plotting style
Jedi.default_plotlyjs!()

# Read data file
data = CSV.read("script1_results.csv")

# Find unique mutation rates
rho = data.rho |> unique

plot(xlabel="l", ylabel="γ")

for x in rho
    l = data[data.rho .== x, :l]
    gamma = data[data.rho .== x, :gamma]
    plot!(l, gamma, label=x, linewidth=2)
end

ext: jl


text:
"""
    search_cv(estimator, X, y, splits, params; [verbosity = 1])

Exhaustive search over the specified sets of parameter values for the `estimator` with features
data `X` and label `y`. The iterable `splits` provides vectors of indices for the training dataset.
The remaining indices are used to create the validation dataset.

Return an array with a tuple for each set of parameters value, where the first entry is a set of
parameter values and the second entry the cross-validation outcome of those values. This outcome is
a dictionary with an entry for the validation dataset and, if the parameter `is_training_metric` is
set in the `estimator`, an entry for the training dataset. Each entry of the dictionary is
another dictionary with an entry for each validation metric in the `estimator`. Each of these
entries is an array that holds the validation metric's value for each dataset, at the last valid
iteration.

# Arguments
* `estimator::LGBMEstimator`: the estimator to be fit.
* `X::Matrix{TX<:Real}`: the features data.
* `y::Vector{Ty<:Real}`: the labels.
* `splits`: the iterable providing arrays of indices for the training dataset.
* `params`: the iterable providing dictionaries of pairs of parameters (Symbols) and values to
    configure the `estimator` with.
* `verbosity::Integer`: keyword argument that controls LightGBM's verbosity. `< 0` for fatal logs
    only, `0` includes warning logs, `1` includes info logs, and `> 1` includes debug logs.
"""
function search_cv(estimator::LGBMEstimator, X::Matrix{TX}, y::Vector{Ty},
                                      splits, params; verbosity::Integer = 1) where {TX<:Real,Ty<:Real}
    n_params = length(params)
    results = Array{Tuple{Dict{Symbol,Any},Dict{String,Dict{String,Vector{Float64}}}}}(undef,n_params)
    for (search_idx, search_params) in enumerate(params)
        log_info(verbosity, "\nSearch: ", search_idx, "\n", search_params, "\n")

        search_estimator = deepcopy(estimator)
        foreach(param -> setfield!(search_estimator, param[1], param[2]), search_params)
        search_results = cv(search_estimator, X, y, deepcopy(splits),
                            verbosity = ifelse(verbosity == 1, 0, verbosity))

        results[search_idx] = (search_params, search_results)
        cv_logsummary(search_results, verbosity)
    end
    log_info(verbosity, "\nSearch finished\n")

    return results
end


text:
# =========================================================================== #
# Training
# =========================================================================== #

struct TrainingOptions
    learning_rate::Float64
    batch_size::Int
    max_iter::Int
    max_norm::Float64
end
function TrainingOptions(;
    learning_rate = 0.1,
    batch_size = 2,
    max_iter = 100,
    max_norm = 1.0)
    return TrainingOptions(learning_rate, batch_size, max_iter, max_norm)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
struct Kernel
    knl::PyObject
    shape::Tuple{Vararg{Int}}
    variables::Array{PyObject}
    function Kernel(node::Node)
        variables = findgemvariables(node)
        gemexpr = togem(node)
        shape = gemexpr.shape
        knl = gemtoop2knl(gemexpr, variables)
        new(knl, shape, variables)
    end
end

function execute(knl::Kernel, variables::Dict{String, T}) where T
    _preprocessvariables(variables)
    return executeop2knl(knl.knl, knl.shape, [variables[var.name] for var in knl.variables])
end

function execute(node::Node, variables::Dict{String, T}) where T
    knl = Kernel(node)
    return execute(knl, variables)
end

function execute(expr::Union{Node, Kernel}, variables::Vararg{Pair})
    return execute(expr, Dict{String, Any}(variables))
end

function _findvariables(visited, tensor::VariableTensor)
    return Set{VariableTensor}([tensor])
end

function _findvariables(visited, tensor::Tensor{T}) where T<:Union{Any, Tensor}
    return union(findvariables.(tensor.value)...)
end

function _findvariables(visited, tensor::ConstantTensor{T}) where T<:VariableTensor
    return Set{VariableTensor}([tensor.value])
end

function _findvariables(visited, node::Node)
    return Set{VariableTensor}()
end

function _findvariables(visited, root::RootNode, found::Vararg{Set{VariableTensor}})
    return RootNode(union(found...))
end

function _findvariables(visited, node::Node, found::Vararg{Set{VariableTensor}})
    return union(found...)
end

function findvariables(node::Node)
    return traversal(node, x->x, _findvariables, nothing, nothing, nothing)
end

function findvariables(n::Number)
    return Set{VariableTensor}()
end

function findgemvariables(node::Node)
    return [togem(var) for var in findvariables(node)]
end

function _preprocessvariables(variables::Dict{String, Any})
    for pair in variables
        if variables[first(pair)] isa Number
            variables[first(pair)] = fill(last(pair), ())
        end
    end
end

function _preprocessvariables(variables::Dict{String, A}) where A<:Array
    return variables
end

max_stars_repo_name: Sohva/TensorDSL

max_issues_repo_name: Sohva/TensorDSL

max_forks_repo_name: Sohva/TensorDSL


text:
module GaussNewton

using LinearAlgebra, ForwardDiff, FiniteDiff, Printf

export optimize, optimize!, has_converged, is_fatal

struct GNTrace
    states::Vector{String}
end
GNTrace() = GNTrace(Vector{String}[])

function update(tr::GNTrace, iter::Int, f0::Real, x0::Vector, store_trace::Bool, show_trace::Bool, show_every::Int = 1)
    os = string("History ", iter, " ", f0, " ", x0)
    update(tr, os, store_trace, show_trace, show_every)
end

function update(tr::GNTrace, str::String, store_trace::Bool, show_trace::Bool, show_every::Int = 1)
    os = str
    if store_trace
        push!(tr.states, os)
    end
    if show_trace
        if show_every == 0
            show(os)
        end
    end
    return
end

struct GNResult{TR,TJ}
    ssr::Real #final sum of squares
    minimizer::TR
    best_jacobian::TJ
    info::Int
    iterations::Int
    reltol::Real
    abstol::Real
    stptol::Real
    tr::GNTrace
    f_calls::Int
    g_calls::Int
end

is_fatal(result::GNResult) = result.info >= 128
has_converged(result::GNResult) = has_converged_abstol(result) || has_converged_reltol(result) || has_converged_stptol(result)
has_converged_stptol(result::GNResult) = result.info & 0x1 == 0x1
has_converged_reltol(result::GNResult) = result.info & 0x2 == 0x2
has_converged_abstol(result::GNResult) = result.info & 0x8 == 0x8

function Base.show(io::IO, r::GNResult)
    @printf io "Results of Gauss-Newton Algorithm\n"
    @printf io " * Minimizer: [%s]\n" join(r.minimizer, ",")
    @printf io " * Sum of squares at Minimum: %f\n" r.ssr
    @printf io " * Iterations: %d\n" r.iterations
    @printf io " * Convergence: %s\n" has_converged(r)
    @printf io " * |x - x'| < %.1e: %s\n" r.stptol has_converged_stptol(r)
    @printf io " * |f(x) - f(x')| / |f(x)| < %.1e: %s\n" r.reltol has_converged_reltol(r)
    @printf io " * |f(x)| < %.1e: %s\n" r.abstol has_converged_abstol(r)
    @printf io " * Function Calls: %d\n" r.f_calls
    @printf io " * Gradient Calls: %d\n" r.g_calls
    @printf io " * Trace: %s\n" join(r.tr.states, "\n")
    return
end

#stptol in: step size for relative convergence test.
# ! reltol, abstol in: value relative/absolute convergence test.
# !    Setting accuracies too small wastes evaluations.
# !    Setting accuracies too large may not get the best solution.
# ! derivstp in: the step for derivatives.
# !    Must be large enough to get some change in the function.
# !    Must be small enough for some accuracy in the derivative.
# ! iprint in: degree of printout, if value is:
# !    0=none, 1=final f0 and x0, 2=iterations, 3=each x tried,
# !    4=all except Jacobian and Hessian, 5=all.
# ! ihist in: if > 0, write history file HistFile using unit ihist
# ! limit  in: maximum number of all evaluations allowed, approximate.
# ! tuning constants:
# !    ZLOW,ZHIGH change bounds of trust region (del) (.5,2).
# !    Small changes in the path set by them greatly affect results.
# !    ZCP Cauchy step size (1).
# !    ZCPMIN minimum Cauchy step size (0.1).
# !    ZCPMAX maximum Cauchy step size (1000).
# !    MXBAD number bad steps with rank-1 updates (1).
# ! NewtStep1st=.true. start with Newton step. Use on linear problems.
# ! mu0,del0,r0,H0,Jac0,iHOOKmax out: final values
# ! dfj/dxi=Jac0(j,i)
# ! iHOOKmax is max number of iterations in hook search
# ! iscale=0, no scaling
# ! iscale=1, variable scaling
# ! iscale=2, fixed scaling based on D0, which must be allocated
# !  and filled by user before call to GN

function optimize(fcn, x0::AbstractArray{T};
    show_trace::Bool = false,
    store_trace::Bool = false,
    NewtStep1st::Bool = false,
    iscale::Int = 0,
    limit::Int = 10000,
    abstol::T = 1e-11,
    derivstp::T = 1e-4,
    reltol::T = 3e-7,
    stptol::T = 1e-4,
    MXBAD::Int = 1,
    ZCONSEC::T = 1.005,
    ZHIGH::T = 2.0,
    ZLOW::T = 0.65,
    ZCP::T = 0.6,
    ZCPMAX::T = 1000.0,
    ZCPMIN::T = 0.18,
    D0 = zeros(T, length(x0)),
    autodiff = :single) where {T}
    fcn! = (out, x) -> copyto!(out, fcn(x))
    r = fcn(x0)
    return optimize!(fcn!, nothing, x0, zeros(eltype(r), length(r)),
        NewtStep1st = NewtStep1st,
        iscale = iscale,
        limit = limit,
        abstol = abstol,
        derivstp = derivstp,
        reltol = reltol,
        stptol = stptol,
        MXBAD = MXBAD,
        ZCONSEC = ZCONSEC,
        ZHIGH = ZHIGH,
        ZCP = ZCP,
        ZCPMAX = ZCPMAX,
        ZCPMIN = ZCPMIN,
        D0 = D0,
        autodiff = autodiff,
        show_trace = show_trace,
        store_trace = store_trace
    )

end
function optimize!(fcn!, x0::AbstractArray{T}, r::AbstractArray{TR};
    show_trace::Bool = false,
    store_trace::Bool = false,
    NewtStep1st::Bool = false,
    iscale::Int = 0,
    limit::Int = 10000,
    abstol::T = 1e-11,
    derivstp::T = 1e-4,
    reltol::T = 3e-7,
    stptol::T = 1e-4,
    MXBAD::Int = 1,
    ZCONSEC::T = 1.005,
    ZHIGH::T = 2.0,
    ZLOW::T = 0.65,
    ZCP::T = 0.6,
    ZCPMAX::T = 1000.0,
    ZCPMIN::T = 0.18,
    D0 = zeros(T, length(x0)),
    autodiff = :single) where {T,TR}

    return optimize!(fcn!, nothing, x0, r,
        NewtStep1st = NewtStep1st,
        iscale = iscale,
        limit = limit,
        abstol = abstol,
        derivstp = derivstp,
        reltol = reltol,
        stptol = stptol,
        MXBAD = MXBAD,
        ZCONSEC = ZCONSEC,
        ZHIGH = ZHIGH,
        ZCP = ZCP,
        ZCPMAX = ZCPMAX,
        ZCPMIN = ZCPMIN,
        D0 = D0,
        autodiff = autodiff,
        show_trace = show_trace,
        store_trace = store_trace
    )

end
function optimize!(fcn!, fcnDer!, x0::AbstractArray{T}, r::AbstractArray{TR};
    show_trace::Bool = false,
    store_trace::Bool = false,
    NewtStep1st::Bool = false,
    iscale::Int = 0,
    limit::Int = 10000,
    abstol::T = 1e-11,
    derivstp::T = 1e-4,
    reltol::T = 3e-7,
    stptol::T = 1e-4,
    MXBAD::Int = 1,
    ZCONSEC::T = 1.005,
    ZHIGH::T = 2.0,
    ZLOW::T = 0.65,
    ZCP::T = 0.6,
    ZCPMAX::T = 1000.0,
    ZCPMIN::T = 0.18,
    D0 = zeros(T, length(x0)),
    autodiff = :central) where {T,TR}
    #a port of fortran90 Netlib.org/misc/gn code
    # AUTHORS: Kenneth Klare (kklare@gmail.com)
    #    !      and Guthrie Miller (guthriemiller@gmail.com)
    #    ! OPEN SOURCE: may be freely used, with author achnowledgment,
    #    !    for any peaceful, beneficial purpose
    #    ! PURPOSE: Minimize sum of squares of residuals r using
    #    !    augmented Gauss-Newton step and Levenberg-Marquardt trust region.
    #    !    Uses finite-difference derivatives and 1-D Jacobian updates.
    #    ! KEYWORDS: Nonlinear Least squares, nonlinear minimization, robust,
    #    !           optimization, finite-difference derivatives,
    #    !           augmented Gauss-Newton, Levenberg-Marquardt.
    #    ! ARGUMENTS:
    #    ! fcn    in: A user subroutine to compute r(1:m) from x(1:n).
    #    !    It is invoked: call FCN(r,x), where
    #    !    x      in. n parameters to evaluate residuals.
    #    !    r      out. m residuals to be reduced. DO NOT SQUARE THEM.
    #    !    It must not change m, n, nor x.
    #    !    Scale the output residuals:
    #    !       r(i)=(fit(x,i)-y(i))/yerror(i), where y is the data.
    #    !    To reject the given x, set the r very large, like sqrt(largest).
    #    !    The starting point (x0) must valid and not be rejected.
    #    !    Declare FCN EXTERNAL so the compiler knows it.
    #    !    An iteration has function and derivative evaluations.
    #    !    Updates use function evaluations only.
    #    !    The program is designed to minimize evaluation count.
    #    ! m      in: number of equations or fitting residuals, the size of r0.
    #    ! n      in: number of variables, the size of x0.
    #    !    Program handles overdetermined, m>n, and underdetermined, m<n,
    #    !    as well as linear and nonlinear equations (m=n).
    #    ! info   out: status sum of bits (LT 128) or single error code (GE 128)
    #    !    Notation: pT q=transpose of p times q, ||q||=norm(q)=sqrt(qT q).
    #    ! CONVERGENCES
    #    !     1 step less than tolerance (stptol).
    #    !     2 rT r relative improvement less than tolerance (reltol).
    #    !     8 rT r sum of squares less than absolute tolerance (abstol).
    #    ! FATAL ERRORS
    #    !    128 input error.
    #    !    129 Jacobian is zero--cannot proceed with unchanging function.
    #    !    130 evaluation limit exceeded without termination.
    #    !    133 hook search error: could not find a suitable mu in 30 tries.
    #    ! x0     in/out n-vector: initial guess and best position found.
    #    ! f0     out scalar: final sum of squares. Think chi-squared.
    #    ! ACCESSIBLE in GN_MOD:
    #    ! r0     out m-vector: final residuals at x0
    #    ! iter   out scalar: actual number of iterations taken
    #    ! nfcn   out scalar: actual number of function evaluations used
    #    ! mu0    out scalar: the "best" L-M parameter used
    #    ! del0   out scalar: the "best" trust region size
    #    ! D0     out n-vector: scale factor, input for iscale=2
    #    ! H0     out n-by-n matrix: the unaugmented Hessian
    #    ! Jac0   out m-by-n array: the "best" Jacobian
    #    ! limitATIONS:
    #    !    Strict descent limits use to reproducible data.
    #    !    x0, r0, and others must be 1-D vectors, not arrays.
    #    !    Large m and n will require large storage and round-off errors.
    #    !    Your residuals should not exceed sqrt(HUGE/m) to avoid overflows.
    #    ! RECOMMENDATIONS:
    #    !    Use double-precision calculation, REAL*8, if possible.
    #    !    Single precision can do 3 to 4 digits in x, at best.
    #    !    Give a good guess or you may find a secondary minimum.
    #    !    Scale parameters and residuals, so chi-square and relative
    #    !    errors are meaningful.
    #    ! CREDITS:
    #    !     * Dennis and Schnabel, "Numerical Methods for Unconstrained
    #    !     Optimization and Nonlinear Equations" (1983),
    #    !     for many ideas and the hook calculation, i.e., mu from phi.
    #    !     * IMSL ZXSSQ for the Jacobian update attributed in
    #    !     Fletcher, "Practical Methods of Optimization" (1980)
    #    !     to Barnes (1965) and Broyden (1965).
    #    !     * x(+) = x-(JacT Jac+mu I)**-1 JacT r augmented Gauss-Newton.
    #    !     * See this for allocations and vector operators in F90.
    #    !     <http://www.stanford.edu/class/me200c/tutorial_90/07_arrays.html>
    #    !     * See "Fortran90 for Fortran77 Programmers" by Clive Page
    #    !     <http://www.star.le.ac.uk/~cgp/f90course/f90.html> 11/26/2001
    #    ! DYNAMICALLY ALLOCATED:
    #    ! Jac  scratch, m-by-n array: Jacobian of derivatives w.r.t. variables.
    #    ! H    scratch, n-by-n array: JacT Jac Hessian without 2nd order derivs.
    #    ! L    scratch, n-by-n array: lower-triangle factor, L LT = H+mu*I.
    #    ! r    scratch, m vector of residuals. ChiSq=Sum(rT r).
    #    ! D,g,sN,s,x   scratch, n vectors.
    #    !    They are scale factor, gradient, Newton step, step tried, position.
    #    !    You must scale by scaling factors D.
    #    !    Use dfj/dxi=Jac(j,i)/D(i) H(i,j)/D(i)/D(j) g(j)/D(j) sN(j)*D(j).
    local hook = false
    local take = false
    local iHOOK = 0
    local nbad = 0
    local nconsec = 0
    local info = 0
    local del = zero(T)
    local f::T
    local f0 = zero(T)
    local sfp::T
    local gnorm::T
    local mu::T
    local mulow::T
    local muhigh::T
    local phi::T
    local phip::T
    local rr::T
    local scal::T
    local scal0::T
    local snewt::T
    local snorm::T
    local temp::T

    m = length(r)
    n = length(x0)

    local tracing = show_trace || store_trace
    local tr = GNTrace()
    r0 = zeros(T, m)
    Jac0 = zeros(T, (m, n))
    H0 = zeros(T, (n, n))
    local Jac = zeros(T, (m, n))
    local H = zeros(T, (n, n))
    local L = zeros(T, (n, n))
    local D = ones(T, n)
    local g = zeros(T, n)
    local s = zeros(T, n)
    local sN = zeros(T, n)
    local epsn12 = sqrt(eps(T) * n)
    local x = copy(x0)
    local newFcnDer! = fcnDer!
    if iscale == 2
        if length(D0) != n
            throw(error("the size of D0 must match n"))
        elseif norm(D0) == 0
            throw(error("the norm of D0 can not be zero"))
        else
            D[1:end] = D0
        end
    end

    if typeof(fcnDer!) == Nothing
        if autodiff == :central
            central_cache = FiniteDiff.JacobianCache(similar(x), similar(r), similar(r))
            newFcnDer! = (J::Matrix, xp::Vector) -> FiniteDiff.finite_difference_jacobian!(J, fcn!, xp, central_cache)
        elseif autodiff == :forward
            jac_cfg = ForwardDiff.JacobianConfig(fcn!, r, x, ForwardDiff.Chunk(x))
            ForwardDiff.checktag(jac_cfg, fcn!, x)
            newFcnDer! = (J::Matrix, xp::Vector) -> ForwardDiff.jacobian!(J, fcn!, deepcopy(r), xp, jac_cfg, Val{false}())
        elseif autodiff == :single
            newFcnDer! = nothing
        else
            throw(DomainError(autodiff, "Invalid automatic differentiation method."))
        end
    end

    nfcn = 0
    nfcnDer = 0
    nrank1 = n
    hook = false
    take = false
    nconsec = 0
    nconsecMax = 0
    iHOOKmax = 0

    hnorm = zero(T)
    js = zeros(T, m)
    itersave = 0
    f_calls = 0
    g_calls = 0
    for iter = 1:limit
        itersave = iter
        info = 0
        nfcn += 1
        f_calls += 1
        fcn!(r, x)
        f = norm(r)^2
        if iter > 1
            take = f < f0
            if take
                nconsec += 1
                nconsec = max(nconsecMax, nconsec)
            else
                nconsec = 0
            end
        end
        if iter == 1
            take = true
        elseif iter == 2 && NewtStep1st && !take
            # ! CAUCHY STEP is distance to Cauchy point in Newton direction.
            # ! Here if Newton step didn't work.
            hook = false
            del = cauchyDist(n, L, g, gnorm, ZCP = ZCP, ZCPMIN = ZCPMIN, ZCPMAX = ZCPMAX)
        else
            #       ! UPDATE TRUST REGION.
            # ! Strict descent requires f<f0.
            # ! Parabolic interpolation/extrapolation using f0, f and slope at x0.
            # ! Update del to be distance to minimum, new del is delFac times old.
            # lastdel = del
            if abs(f - f0) <= max(reltol * max(f, f0), abstol)
                info = 2
            end
            # ! sfp is the step length s times the slope fp of f in the direction s.
            # ! it would be the change in f for negligible curvature.
            # ! sfp is negative because the step is in a decreasing direction.
            sfp = 2 * dot(g, s)
            rr = (f - f0) / sfp
            # ! rr is bounded above by 1, when f < f0 and f behaves linearly,
            # ! and unbounded below when f > f0
            delFac = 1 / (2 * max(1 - rr, eps(T)))
            del = snorm * min(max(delFac, ZLOW), ZHIGH)
            # ! we also increase trust radius based on number of consecutive takes.
            # ! Test problem #36 illustrates the need for this.
            del = del * ZCONSEC^nconsec
        end
        if take
            # ! SAVE THE BEST.
            #    if(iprint.ge.4) write(*,*) 'saving best'\
            tracing && update(tr, "saving best", store_trace, show_trace)
            nbad = 0
            x0[1:end] = x
            f0 = f
            # mu0 = mu
            # del0 = del
            D0[1:end] = D
            H0[1:end, 1:end] = H
            Jac0[1:end, 1:end] = Jac

            # if(ihist.gt.0) write(ihist,*) iter,f0,x0
            tracing && update(tr, iter, f0, x0, store_trace, show_trace)
        elseif iter != 2 || !NewtStep1st
            # ! Overshot the mark. Try, then get a new Jacobian.
            hook = true
            if nrank1 > 0
                nbad += 1
            end
            #if(iprint.ge.4) write(*,*)'Overshot nrank1,nbad',nrank1,nbad
            tracing && update(tr, string("Overshot ", nrank1, " ", nbad), store_trace, show_trace)
        end
        # ! Some convergence tests.
        if iter > 1
            if snorm <= stptol
                info = info + 1
            end
        end
        if f < abstol
            info = info + 8
        end
        done = (info != 0 && nrank1 == 0) || info >= 8
        if !done && nfcn >= limit
            info = 130
        end
        if iter > 1 && (done || info >= 128)
            break
        end
        # ! ---Jacobian dfi/dxj, gradient g=JacT r, Hessian H=JacT Jac+2nd order
        # ! Jacobian update when not stale nor final, else a full Jacobian.
        if nbad <= MXBAD && info == 0 && nrank1 < n
            if !take
                @goto L40
            end
            # ! Rank-1 update to Jacobian. Jac(new) = Jac+((r-r0-Jac s) sT)/(sT s).
            #     if(iprint.ge.4) write(*,*)
            # 1     'Rank=1 Jacobian update: nrank1,nbad,snorm',nrank1,nbad,snorm
            tracing && update(tr, string("Rank=1 Jacobian update: ", nrank1, " ", nbad, " ", snorm), store_trace, show_trace)
            nrank1 = nrank1 + 1
            js = Jac * s
            #Jac = Jac +matmul(reshape(r-r0-matmul(Jac,s),(m,1)),               reshape(s/snorm**2,(1,n))) //!outer product
            for u = 1:m
                for v = 1:n
                    Jac[u, v] += (r[u] - r0[u] - js[u]) * (s[v] / (snorm * snorm))
                end
            end
            if take
                r0[1:end] = r
            end
        else

            # ! Full Jacobian.
            # ! Step away from zero to avoid crossing it.
            #          if(iprint.ge.4) write(*,*)'Full Jacobian: nrank1,nbad,take'
            #      1                                            ,nrank1,nbad,take
            tracing && update(tr, string("Full Jacobian: ", nrank1, " ", nbad, " ", take), store_trace, show_trace)
            nrank1 = 0
            nbad = 0
            if take
                r0[1:end] = r
            end
            nfcn += n
            g_calls += 1
            if iscale == 1
                # ! variable scale
                scal0 = dot(s, D)^2
                @. D = max((abs(x0) + D) / 2, stptol)

                # if(iprint.ge.4) write(*,*) 'variable scale D',D
                tracing && update(tr, string("variable scale ", D), store_trace, show_trace)
                scal = dot(s, D)^2
                if scal > zero(T)
                    del = del * sqrt(scal0 / scal)
                end
            end
            if typeof(newFcnDer!) == Nothing
                computeGNJacobian(fcn!, x0, D, r0, r, Jac, derivstp = derivstp)
            else
                newFcnDer!(Jac, x0)
                for j = 1:n
                    dj = D[j]
                    Jac[1:end, j] .*= dj
                end
            end

        end

        # ! Gradient and Hessian.
        g = Jac' * r0
        H = Jac' * Jac
        #  if(iprint.ge.5) then
        #     do i=1,m
        #        write(*,*) 'J',Jac(i,:)
        #     enddo
        #     do j=1,n
        #        write(*,*) 'H',H(j,:)
        #     enddo
        #  endif
        gnorm = norm(g)
        #       if(iprint.ge.4) write(*,*) 'gnorm,g',gnorm,g
        tracing && update(tr, string(gnorm, " ", g), store_trace, show_trace)
        # ! L1 norm (max of row sums of abs) of H(i,j).
        # ! H=JacT Jac symmetric.
        hnorm = zero(T)
        for j = 1:n
            lsum = sum(x -> abs(x), H[1:end, j])
            hnorm = max(hnorm, lsum)
        end
        # ! Get a small number for further augmentation, check underflow.
        hnorm = hnorm * epsn12
        if hnorm == 0
            info = 129
            break #MainLoop
        end
        # ! Find bad rows and ignore them.
        for j = 1:n
            if H[j, j] <= 0
                H[j, j] = one(T)
                H[j, :] .= zero(T)
                H[:, j] .= zero(T)
            end
        end

        # ! --- GET NEWTON STEP, H sN = -g ---
        mu = zero(T)
        # ! solve (H + mu I)s = g for step s, possibly augmenting diagonal of H
        GN_CHOL(n, mu, H, g, L, sN, hnorm)
        snewt = norm(sN)
        if snewt <= del
            hook = false
        end
        # if(iprint.ge.4) write(*,*) 'snewt,mu,sN',snewt,mu,sN
        tracing && update(tr, string(snewt, " ", mu, " ", sN), store_trace, show_trace)
        if iter == 1
            if NewtStep1st
                # ! Try Newton step first.
                del = snewt
                s[1:end] = sN
                snorm = snewt
                # if(iprint.ge.4) write(*,*)'Taking Newton step first'
                tracing && update(tr, "Taking Newton step first", store_trace, show_trace)
                @goto L100
            else
                # ! more conservative approach, try Cauchy step first
                hook = false
                del = cauchyDist(n, L, g, gnorm, ZCP = ZCP, ZCPMIN = ZCPMIN, ZCPMAX = ZCPMAX)
            end
        end
        @label L40  # continue ! --- CHOOSE NEWTON OR HOOK STEPS -------------------
        if !hook
            # ! Allow Newton step to be up to the current trust radius del
            temp = 1
            if snewt > del
                temp = del / snewt
            end
            s = sN .* temp
            snorm = snewt * temp
            #     if(iprint.ge.4) write(*,*) 'Step in Newton direction:'
            # 1    ,' snewt,del,snorm',snewt,del,snorm
            tracing && update(tr, string("Step in Newton direction: ", snewt, " ", del, " ", snorm), store_trace, show_trace)
        else
            # ! --- Hook search ---
            # ! Find step of length del by finding mu
            # ! that gives snorm = ||s|| = ||(H + mu I)**-1 g|| = del.
            # ! Because del = ||(H + mu I)**-1 g|| and H is positive,
            # ! del is less than ||g||/mu.
            muhigh = gnorm / del
            mulow = zero(T)
            # ! mulow <= mu <= muhigh
            # muStart = mu
            # ! REPEAT UNTIL abs(del-||s||)<.05 or mulow>muhigh.
            #   HookLoop: do iHOOK=1,30
            for iHOOK = 1:30
                GN_CHOL(n, mu, H, g, L, s, hnorm)
                mulow = max(mulow, mu)
                snorm = norm(s)
                if abs(snorm - del) <= 0.05 * del || mulow >= muhigh
                    break
                end
                # 1       exit HookLoop
                phi = snorm - del
                # ! phi<0 indicates mu is too large, use this mu to update muhigh
                if phi < 0
                    muhigh = min(mu, muhigh)
                end
                if phi > 0
                    mulow = max(mu, mulow)
                end
                # ! Want sT (H + mu I)**-1 s = ||L**-1 s||**2.
                # ! Start with x = L**-1 s.
                GN_LSOLVE(n, L, x, s)

                phip = zero(T)
                if snorm > zero(T)
                    sx = norm(x)^2
                    phip = -sx / snorm
                end
                #    ! As mu increases the step size decreases, so phip must be negative.
                if phip < zero(T)
                    mu = mu - (snorm / del) * (phi / phip)
                end
                mu = max(min(mu, muhigh), mulow)
            end
            #! End of REPEAT.

            #     if(iprint.ge.4.or.iHOOK.ge.30) then
            #        write(*,*) 'iHOOK,iHOOKmax,muStart,mu,lastdel,del,snorm'
            # 1                 ,iHOOK,iHOOKmax,muStart,mu,lastdel,del,snorm
            tracing && update(tr, string(iHOOK, " ", iHOOKmax, " ", mu, " ", del, " ", snorm), store_trace, show_trace)
            if iHOOK >= 30
                info = 133
                if iHOOK > 30
                    break # MainLoop
                end
            end
            iHOOKmax = max(iHOOKmax, iHOOK)
        end

        @label L100 #            continue ! TAKE THE STEP.
        # if(iprint.ge.4) write(*,*) 'Taking step s',s
        tracing && update(tr, string("Taking step ", s), store_trace, show_trace)
        @. x = x0 + s * D
    end
    # abs(f - f0) <= max(reltol * max(f, f0), abstol)
    return f0, GNResult(f0, x0, Jac0, info, itersave, reltol, abstol, stptol, tr, f_calls, g_calls)
end

function computeGNJacobian(fcn, x0::AbstractArray{T}, D::AbstractArray{T}, r0::AbstractArray{T}, r::AbstractArray{T}, Jac::AbstractMatrix{T}; derivstp::T = 1e-4) where {T}
    #start jacobian
    for j = 1:length(x0)
        hold = x0[j]
        step = copysign(derivstp, hold)
        x0[j] += step * D[j]
        fcn(r, x0) #jacevaluation
        Jac[:, j] = (r - r0) / step
        x0[j] = hold
    end
    #end jacobian => could be analytical
end

# !-----------------------------------------------------------------------
# !+Cauchy distance is distance along gradient to minimum.
# ! Cauchy step is Cauchy distance in the Newton direction.
# ! Cauchy distance = -||g||**2/(gT H g) ||g||, gT H g = ||LT g||**2.
# ! Use ZCP * Cauchy distance.
# ! ||v|| = sqrt(vT v) is L2 or Euclidean norm of vector v.
function cauchyDist(n::Int, L::AbstractArray{T}, g::AbstractArray{T}, gnorm::T; ZCP::T, ZCPMIN::T, ZCPMAX::T)::T where {T}
    # use GN_MOD, only: ZCP,ZCPMIN,ZCPMAX,iprint
    # implicit none
    # intrinsic max,min,sum
    # real*8 L(n,n),g(n),gnorm,del,temp
    if gnorm == zero(T)
        return ZCPMIN
    end
    temp = zero(T)
    # ! calculate temp = gT H g/||g||**2 = ||LT g||**2/||g||**2
    for j = 1:n
        sum = dot(L[j:end, j], g[j:end])
        temp += (sum / gnorm)^2
    end
    temp = ZCP * gnorm / temp
    del = max(ZCPMIN, min(temp, n * ZCPMAX))
    #  if(iprint.ge.4) write(*,*)
    # 1 'Cauchy step,del,stpmax',temp,del,n*ZCPMAX
    return del
end #function GN_CauchyDist

# !-----------------------------------------------------------------------
# !+Cholesky decomposition and solution, (H + (mu+add)*I) s = -g.
# ! Decomposition to lower triangle L without addition if it works.
# ! Break as soon as OK.
function GN_CHOL(n::Int, mu::T, H::AbstractMatrix{T}, g::AbstractArray{T}, L::AbstractMatrix{T}, s::AbstractArray{T}, add::T) where {T}
    # implicit none
    # intrinsic sqrt,min,sum
    # integer n,i,iadd,j
    # real*8 add,mu,tmp,H(n,n),g(n),L(n,n),s(n)
    #loop1:
    isCycle = false
    for iadd = 1:n
        for j = 1:n
            for i = j:n
                sum = dot(L[i, 1:j-1], L[j, 1:j-1])
                L[i, j] = H[j, i] - sum
            end
            tmp = L[j, j] + mu
            if tmp <= 0
                mu = mu + add
                isCycle = true # loop1
                break
            end
            tmp = sqrt(tmp)
            L[j, j] = tmp
            L[j+1:end, j] ./= tmp
        end
        if !isCycle
            break
        else
            isCycle = false
        end
    end
    # ! forward row reduction and backward substitution
    for i = 1:n
        sum = dot(L[i, 1:i-1], s[1:i-1])
        s[i] = (-g[i] - sum) / L[i, i]
    end
    for i = n:-1:1
        sum = dot(L[i+1:end, i], s[i+1:end])
        s[i] = (s[i] - sum) / L[i, i]
    end
end

# !-----------------------------------------------------------------------
# !+solve L x = y. x and y may be the same vector, L lower triangle.
function GN_LSOLVE(n::Int, L::AbstractMatrix{T}, x::AbstractArray{T}, y::AbstractArray{T}) where {T}
    # implicit none
    # intrinsic sum
    # integer n,i
    # real*8 L(n,n),x(n),y(n)
    for i = 1:n
        sum = dot(L[i, 1:i-1], x[1:i-1])
        x[i] = (y[i] - sum) / L[i, i]
    end
end #subroutine GN_LSOLVE

end # module


text:
using FixedPointAcceleration
func(x) = cos.(x)
Inputs = 1.1
# We should not have any slowdown from starting in simple for two. Then SEA for 3. Then Simple. Then VEA until the end.
# The reason is that SEA/VEA does nothing of interest until 7 by default anyway.
fp_chain      = fixed_point(func, Inputs; Algorithm = :Simple, MaxIter = Integer(2))
fp_chain      = fixed_point(func, fp_chain; Algorithm = :SEA, MaxIter = Integer(3))
fp_chain      = fixed_point(func, fp_chain; Algorithm = :Simple, MaxIter = Integer(1))
fp_chain      = fixed_point(func, fp_chain; Algorithm = :VEA, MaxIter = Integer(100))

fp_nochain = fixed_point(func, Inputs; Algorithm = :VEA, MaxIter = Integer(100))
fp_chain.Iterations_ == fp_nochain.Iterations_
all(abs.(fp_nochain.Inputs_ .- fp_chain.Inputs_) .< 1e-14)

func(x) = [0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2]]
Inputs = [1.1,2.2]

fp_chain      = fixed_point(func, Inputs; Algorithm = :Simple, MaxIter = Integer(2))
fp_chain      = fixed_point(func, fp_chain; Algorithm = :MPE, MaxIter = Integer(3))
fp_chain      = fixed_point(func, fp_chain; Algorithm = :Simple, MaxIter = Integer(1))
fp_chain      = fixed_point(func, fp_chain; Algorithm = :RRE, MaxIter = Integer(100))

fp_nochain = fixed_point(func, Inputs; Algorithm = :RRE, MaxIter = Integer(100))
fp_chain.Iterations_ == fp_nochain.Iterations_
all(abs.(fp_nochain.Inputs_ .- fp_chain.Inputs_) .< 1e-14)


text:
#! format: off

abstract type AbstractStorageFormulation <: AbstractDeviceFormulation end
abstract type AbstractEnergyManagement  <: AbstractStorageFormulation end
struct BookKeeping <: AbstractStorageFormulation end
struct BookKeepingwReservation <: AbstractStorageFormulation end
struct BatteryAncillaryServices <: AbstractStorageFormulation end
struct EnergyTarget <: AbstractEnergyManagement end

get_variable_sign(_, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = NaN
########################### ActivePowerInVariable, Storage #################################

get_variable_binary(::ActivePowerInVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_expression_name(::ActivePowerInVariable, ::Type{<:PSY.Storage}) = :nodal_balance_active

get_variable_lower_bound(::ActivePowerInVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_upper_bound(::ActivePowerInVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_input_active_power_limits(d).max
get_variable_sign(::ActivePowerInVariable, d::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = -1.0

########################### ActivePowerOutVariable, Storage #################################

get_variable_binary(::ActivePowerOutVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_expression_name(::ActivePowerOutVariable, ::Type{<:PSY.Storage}) = :nodal_balance_active

get_variable_lower_bound(::ActivePowerOutVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_upper_bound(::ActivePowerOutVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_output_active_power_limits(d).max
get_variable_sign(::ActivePowerOutVariable, d::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = 1.0

############## ReactivePowerVariable, Storage ####################
get_variable_sign(::PowerSimulations.ReactivePowerVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = 1.0
get_variable_binary(::ReactivePowerVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_expression_name(::ReactivePowerVariable, ::Type{<:PSY.Storage}) = :nodal_balance_reactive

############## EnergyVariable, Storage ####################

get_variable_binary(::EnergyVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_upper_bound(::EnergyVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_state_of_charge_limits(d).max
get_variable_lower_bound(::EnergyVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_state_of_charge_limits(d).min
get_variable_initial_value(::EnergyVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_initial_energy(d)

############## ReserveVariable, Storage ####################

get_variable_binary(::ReserveVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = true

get_efficiency(v::T, var::Type{<:InitialConditionType}) where T <: PSY.Storage = PSY.get_efficiency(v)

############## EnergyShortageVariable, Storage ####################

get_variable_binary(::EnergyShortageVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_lower_bound(::EnergyShortageVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_upper_bound(::EnergyShortageVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_rating(d)

############## EnergySlackDown, Storage ####################

get_variable_binary(::EnergySurplusVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_upper_bound(::EnergySurplusVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_lower_bound(::EnergySurplusVariable, d::PSY.Storage, ::AbstractStorageFormulation) = - PSY.get_rating(d)
#! format: on

################################## output power constraints#################################

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerOutVariable},
    ::Type{T},
    ::Type{<:BookKeeping},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.Storage}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerOutVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerOutVariable, T),
            limits_func = x -> PSY.get_output_active_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerInVariable},
    ::Type{T},
    ::Type{<:BookKeeping},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.Storage}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerInVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerInVariable, T),
            limits_func = x -> PSY.get_input_active_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerOutVariable},
    ::Type{T},
    ::Type{<:AbstractStorageFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.Storage}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerOutVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerOutVariable, T),
            bin_variable_names = [make_variable_name(ReserveVariable, T)],
            limits_func = x -> PSY.get_output_active_power_limits(x),
            constraint_func = device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerInVariable},
    ::Type{T},
    ::Type{<:AbstractStorageFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.Storage}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerInVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerInVariable, T),
            bin_variable_names = [make_variable_name(ReserveVariable, T)],
            limits_func = x -> PSY.get_input_active_power_limits(x),
            constraint_func = reserve_device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function adds the reactive  power limits of generators when there are CommitmentVariables
"""
function add_constraints!(
    optimization_container::OptimizationContainer,
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    devices::IS.FlattenIteratorWrapper{St},
    model::DeviceModel{St, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    constraint_infos = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_reactive_power_limits(d)
        constraint_infos[ix] = DeviceRangeConstraintInfo(name, limits)
    end

    device_range!(
        optimization_container,
        RangeConstraintSpecInternal(
            constraint_infos,
            make_constraint_name(RangeConstraint, ReactivePowerVariable, St),
            make_variable_name(ReactivePowerVariable, St),
        ),
    )
    return
end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{St},
    formulation::AbstractStorageFormulation,
) where {St <: PSY.Storage}
    storage_energy_initial_condition!(optimization_container, devices, formulation)
    return
end

######################### Initialize Functions for Storage #################################
# TODO: This IC needs a cache for Simulation over long periods of tim
function storage_energy_initial_condition!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::D,
) where {T <: PSY.Storage, D <: AbstractStorageFormulation}
    key = ICKey(InitialEnergyLevel, T)
    _make_initial_conditions!(
        optimization_container,
        devices,
        D(),
        EnergyVariable(),
        key,
        _make_initial_condition_energy,
        _get_variable_initial_value,
        StoredEnergy,
    )

    return
end

############################ Energy Capacity Constraints####################################
function energy_capacity_constraints!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{St},
    model::DeviceModel{St, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    constraint_infos = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_state_of_charge_limits(d)
        constraint_info = DeviceRangeConstraintInfo(name, limits)
        add_device_services!(constraint_info, d, model)
        constraint_infos[ix] = constraint_info
    end

    device_range!(
        optimization_container,
        RangeConstraintSpecInternal(
            constraint_infos,
            make_constraint_name(ENERGY_CAPACITY, St),
            make_variable_name(ENERGY, St),
        ),
    )
    return
end

############################ book keeping constraints ######################################

function DeviceEnergyBalanceConstraintSpec(
    ::Type{<:EnergyBalanceConstraint},
    ::Type{EnergyVariable},
    ::Type{St},
    ::Type{<:AbstractStorageFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {St <: PSY.Storage}
    return DeviceEnergyBalanceConstraintSpec(;
        constraint_name = make_constraint_name(ENERGY_LIMIT, St),
        energy_variable = make_variable_name(ENERGY, St),
        initial_condition = InitialEnergyLevel,
        pin_variable_names = [make_variable_name(ACTIVE_POWER_IN, St)],
        pout_variable_names = [make_variable_name(ACTIVE_POWER_OUT, St)],
        constraint_func = energy_balance!,
    )
end

############################ reserve constraints ######################################

function reserve_contribution_constraint!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{<:PM.AbstractPowerModel},
    ::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.Storage, D <: AbstractStorageFormulation}
    constraint_infos_up = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    constraint_infos_dn = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    constraint_infos_energy = Vector{ReserveRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        up_info = DeviceRangeConstraintInfo(name, PSY.get_output_active_power_limits(d))
        down_info = DeviceRangeConstraintInfo(name, PSY.get_input_active_power_limits(d))
        energy_info = ReserveRangeConstraintInfo(
            name,
            PSY.get_state_of_charge_limits(d),
            PSY.get_efficiency(d),
        )
        add_device_services!(up_info, down_info, d, model)
        add_device_services!(energy_info, d, model)
        constraint_infos_energy[ix] = energy_info
        constraint_infos_up[ix] = up_info
        constraint_infos_dn[ix] = down_info
    end

    reserve_power_ub!(
        optimization_container,
        constraint_infos_up,
        constraint_infos_dn,
        make_constraint_name(RESERVE_POWER, T),
        (make_variable_name(ACTIVE_POWER_IN, T), make_variable_name(ACTIVE_POWER_OUT, T)),
    )

    reserve_energy_ub!(
        optimization_container,
        constraint_infos_energy,
        make_constraint_name(RESERVE_ENERGY, T),
        make_variable_name(ENERGY, T),
    )

    return
end

############################ Energy Management constraints ######################################
function energy_target_constraint!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, EnergyTarget},
    ::Type{<:PM.AbstractPowerModel},
    ::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.Storage}
    parameters = model_has_parameters(optimization_container)
    use_forecast_data = model_uses_forecasts(optimization_container)
    time_steps = model_time_steps(optimization_container)
    target_forecast_label = "storage_target"
    constraint_infos_target = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    if use_forecast_data
        for (ix, d) in enumerate(devices)
            ts_vector_target =
                get_time_series(optimization_container, d, target_forecast_label)
            constraint_info_target =
                DeviceTimeSeriesConstraintInfo(d, x -> PSY.get_rating(x), ts_vector_target)
            constraint_infos_target[ix] = constraint_info_target
        end
    else
        for (ix, d) in enumerate(devices)
            ts_vector_target =
                length(time_steps) == 1 ? [PSY.get_storage_target(d)] :
                vcat(zeros(time_steps[end - 1]), PSY.get_storage_target(d))
            constraint_info_target =
                DeviceTimeSeriesConstraintInfo(d, x -> PSY.get_rating(x), ts_vector_target)
            constraint_infos_target[ix] = constraint_info_target
        end
    end

    if parameters
        energy_target_param!(
            optimization_container,
            constraint_infos_target,
            make_constraint_name(ENERGY_TARGET, T),
            (
                make_variable_name(ENERGY, T),
                make_variable_name(ENERGY_SHORTAGE, T),
                make_variable_name(ENERGY_SURPLUS, T),
            ),
            UpdateRef{T}(TARGET, target_forecast_label),
        )
    else
        energy_target!(
            optimization_container,
            constraint_infos_target,
            make_constraint_name(ENERGY_TARGET, T),
            (
                make_variable_name(ENERGY, T),
                make_variable_name(ENERGY_SHORTAGE, T),
                make_variable_name(ENERGY_SURPLUS, T),
            ),
        )
    end

    constraint_infos = Vector{DeviceRangeConstraintInfo}()
    for d in devices
        op_cost = PSY.get_operation_cost(d)
        if PSY.get_energy_shortage_cost(op_cost) == 0.0
            dev_name = PSY.get_name(d)
            limits = (min = 0.0, max = 0.0)
            constraint_info = DeviceRangeConstraintInfo(dev_name, limits)
            push!(constraint_infos, constraint_info)
            @warn(
                "Device $dev_name has energy shortage cost set to 0.0, as a result the model will turnoff the EnergyShortageVariable to avoid infeasible/unbounded problem."
            )
        end
    end
    if !isempty(constraint_infos)
        device_range!(
            optimization_container,
            RangeConstraintSpecInternal(
                constraint_infos,
                make_constraint_name(RangeConstraint, EnergyShortageVariable, T),
                make_variable_name(EnergyShortageVariable, T),
                Vector{Symbol}(),
            ),
        )
    end

    return
end

function AddCostSpec(
    ::Type{PSY.BatteryEMS},
    ::Type{EnergyTarget},
    optimization_container::OptimizationContainer,
)
    return AddCostSpec(;
        variable_type = ActivePowerOutVariable,
        component_type = PSY.BatteryEMS,
        variable_cost = PSY.get_variable,
        multiplier = OBJECTIVE_FUNCTION_POSITIVE,
    )
end


text:
module BenchForeachSeqDouble

using BenchmarkTools
using Referenceables
using ThreadsX

@inline function manual_double!(ys)
    @simd ivdep for i in eachindex(ys)
        @inbounds ys[i] *= 2
    end
end

@inline function manual_nested_double!(ys::AbstractMatrix)
    for j in axes(ys, 2)
        @simd ivdep for i in axes(ys, 1)
            @inbounds ys[i, j] *= 2
        end
    end
end

double!(y) = y[] *= 2

suite = BenchmarkGroup()

let s0 = suite["linear"] = BenchmarkGroup()
    ys = randn(2^10)
    s0["man"] = @benchmarkable manual_double!(ys) setup = (ys = copy($ys))
    s1 = s0["tx"] = BenchmarkGroup()
    for simd in [false, true, :ivdep]
        s1[:simd=>simd] = @benchmarkable(
            ThreadsX.Implementations.foreach_linear_seq(
                double!,
                referenceable(ys),
                $(Val(simd)),
            ),
            setup = (ys = copy($ys)),
        )
    end
end

let s0 = suite["cartesian"] = BenchmarkGroup()
    # `size(ys, 1)` seems to be large enough to make this benchmark
    # somewhat stable
    ys = randn(2^8, 2^5)
    s0["man"] = @benchmarkable manual_nested_double!(ys') setup = (ys = copy($ys))
    s1 = s0["tx"] = BenchmarkGroup()
    for simd in [false, true, :ivdep]
        s1[:simd=>simd] = @benchmarkable(
            ThreadsX.Implementations.foreach_cartesian_seq(
                double!,
                referenceable(ys'),
                $(Val(simd)),
            ),
            setup = (ys = copy($ys)),
        )
    end
end

end  # module
BenchForeachSeqDouble.suite


text:
using SoftPredicates
using Test

@testset "SoftPredicates.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
############################################################################################
#     a
#    ___
# b |   | c
#   |_d_|
#   |   |
# e |___| f
#     g
############################################################################################
const zero = Set("abcefg")
const one = Set("cf")
const two = Set("acdeg")
const three = Set("acdfg")
const four = Set("bcdf")
const five = Set("abdfg")
const six = Set("abdefg")
const seven = Set("acf")
const eight = Set("abcdefg")
const nine = Set("abcdfg")
const numbers = Dict(zero => '0', one => '1', two => '2', three => '3', four => '4',
                     five => '5', six => '6', seven => '7', eight => '8', nine => '9')

isOne(s::String) :: Bool = (length(s) == 2)
isFour(s::String) :: Bool = (length(s) == 4)
isSeven(s::String) :: Bool = (length(s) == 3)
isEight(s::String) :: Bool = (length(s) == 7)
haslengthsix(s::String) :: Bool = (length(s) == 6)
haslengthfive(s::String) :: Bool = (length(s) == 5)

struct Combination
  digits :: Vector{String}
  message :: Vector{String}
end

function getInitialMap()
  return Dict('a' => Set("abcdefg"), 'b' => Set("abcdefg"),
              'c' => Set("abcdefg"), 'd' => Set("abcdefg"),
              'e' => Set("abcdefg"), 'f' => Set("abcdefg"),
              'g' => Set("abcdefg"))
end

# Chacks if map is fully solved
function solved(d::Dict{Char,Set{Char}}) :: Bool
  for entry in d
    if length(entry.second) != 1
      return false
    end
  end
  return true
end

# Eliminates the letters we know from other lists.
function clean!(d::Dict{Char,Set{Char}})
  for key in d
    if length(key.second) == 1
      clean!(d,first(key.second))
    end
  end
end

# Remove a letter from all except the ones which have been found.
function clean!(d::Dict{Char,Set{Char}}, c::Char)
  for key in d
    if length(key.second) > 1
      delete!(key.second, c)
    end
  end
end

# Adds known combination to the map
function add!(d::Dict{Char,Set{Char}},key::Char,vals::Set{Char})
  d[key] = intersect(d[key],vals)
end

function add!(d::Dict{Char,Set{Char}},keys::Set{Char},vals::Set{Char})
  for key in keys
    add!(d,key,vals)
  end
end

# Find the map that decifers the combination
function findMap(comb::Combination) :: Dict{Char,Char}
  digits  = comb.digits
  posible = getInitialMap()

  # Initial knowledge:
  add!(posible, one  , Set(digits[findfirst(isOne, digits)]))   # One
  add!(posible, four , Set(digits[findfirst(isFour, digits)]))  # Four
  add!(posible, seven, Set(digits[findfirst(isSeven, digits)])) # Seven
  add!(posible, eight, Set(digits[findfirst(isEight, digits)])) # Eight
  add!(posible, intersect(two,three,five), Set(intersect(digits[findall(haslengthfive, digits)]...)))
  add!(posible, intersect(zero,six,nine), Set(intersect(digits[findall(haslengthsix, digits)]...)))

  # Iterate until we decifer all:
  it = 0
  while !solved(posible) && it < 10
    clean!(posible)
    it += 1
  end

  # Create output (inverse map)
  res :: Dict{Char,Char} = Dict()
  for key in posible
    res[first(key.second)] = key.first
  end

  return res
end

# Takes a coded digit and decifers it (returns digit as char)
function decifer(d::Dict{Char,Char}, s::String) :: Char
  res = Set(map(c -> d[c], collect(s)))
  return numbers[res]
end

############################################################################################

function parseWords(s::String) :: Vector{String}
  return split(s, " ", keepempty=false)
end

function parseCombination(s::String) :: Combination
  parts = split(s, "|", keepempty=false)
  digits = parseWords(String(parts[1]))
  message = parseWords(String(parts[2]))
  return Combination(digits,message)
end

function readInput() :: Vector{Combination}
  f = open("./inputs/day8.txt")

  combinations :: Vector{Combination} = []
  for comb in eachline(f)
    push!(combinations, parseCombination(comb))
  end
  return combinations
end

############################################################################################

function problem1(combs::Vector{Combination}) :: Int
  n = 0
  for comb in combs
    for word in comb.message
      if (isOne(word) || isFour(word) || isSeven(word) || isEight(word))
        n += 1
      end
    end
  end
  return n
end

function problem2(combs::Vector{Combination}) :: Int
  n = 0
  for comb in combs
    solvemap = findMap(comb)
    n += parse(Int,join(map(word -> decifer(solvemap,word), comb.message)))
  end
  return n
end

combinations = readInput()
println("Problem 1: ", problem1(combinations))
println("Problem 1: ", problem2(combinations))


text:
export project_l1_Duchi!

"""
    project_l1_Duchi!(v::Union{Vector{TF},Vector{Complex{TF}}}, b::TF)

Projects ector onto L1 ball of specified radius.

w = ProjectOntoL1Ball(v, b) returns the vector w which is the solution
  to the following constrained minimization problem:

   min   ||w - v||_2
   s.t.  ||w||_1 <= b.

  That is, performs Euclidean projection of v to the 1-norm ball of radius
  b.

Author: John Duchi (jduchi@cs.berkeley.edu)
Translated (with some modification) to Julia 1.1 by Bas Peters
"""

function project_l1_Duchi!(v::Union{Vector{TF},Vector{Complex{TF}}}, b::TF) where {TF<:Real}
  b <= TF(0) && error("Radius of L1 ball is negative")
  norm(v, 1) <= b && return v

  lv = length(v)
  u  = similar(v)
  sv = Vector{TF}(undef, lv)

  #use RadixSort for Float32 (short keywords)
  copyto!(u, v)
  u .= abs.(u)
  u  = convert(Vector{TF},u)
  if TF==Float32
    sort!(u, rev=true, alg=RadixSort)
  else
    u = sort!(u, rev=true, alg=QuickSort)
  end

  # if TF==Float32
  #   u = sort!(abs.(u), rev=true, alg=RadixSort)
  # else
  #   u = sort!(abs.(u), rev=true, alg=QuickSort)
  # end

  cumsum!(sv, u)

  # Thresholding level
  temp  = TF(1.0):TF(1.0):TF(lv)
  rho   = max(1, min(lv, findlast(u .> ((sv.-b) ./ temp ) ) ))::Int
  theta = max.(TF(0) , (sv[rho] .- b) ./ rho)::TF

  # Projection as soft thresholding
  v .= sign.(v) .* max.(abs.(v) .- theta, TF(0))

  return v
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
#=
This file contains functions taken and modified from Julia/base/show.jl
license of show.jl:
Copyright (c) 2009-2016: Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors:
https://github.com/JuliaLang/julia/contributors
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

import Base: indent_width, uni_ops, expr_infix_wide
import Base: all_ops, expr_calls, show_block
import Base: show_list, show_enclosed_list, operator_precedence
import Base: is_quoted, is_expr, expr_infix_any
import Base: show_call, show_unquoted, show

import Sugar: ASTIO, LazyMethod, typename, functionname, _typename, show_name
import Sugar: supports_overloading, show_type, show_function, is_native_type, isintrinsic

@compat abstract type CIO <: ASTIO end
immutable EmptyCIO <: CIO
end

immutable EmptyStruct
    # Emtpy structs are not supported in OpenCL, which is why we emit a struct
    # with one floating point field
    x::Int32
    EmptyStruct() = new()
end

# helper function to fake a return type to type inference for intrinsic function stabs
@noinline function ret{T}(::Type{T})::T
    unsafe_load(Ptr{T}(C_NULL))
end
# Number types
# Abstract types
# for now we use Int, more accurate would be Int32. But to make things simpler
# we rewrite Int to Int32 implicitely like this!
const int = Int32
# same goes for float
const float = Float32
const uint = UInt32
const uchar = UInt8

const ints = (UInt64, UInt32, UInt8, Int64, Int32, Int8)
const floats = (Float32, Float64)
const numbers = (ints..., floats..., Bool)

const Ints = Union{ints...}
const Floats = Union{floats...}
const Numbers = Union{numbers...}

const vector_lengths = (2, 3, 4, 8, 16)
# Maybe we should remove this and leave it to `is_fixedsize_array`, which is extensible
_vecs = []
for i in vector_lengths, T in numbers
    T == Bool && continue
    push!(_vecs, NTuple{i, T})
    push!(_vecs, SVector{i, T})
end
const vecs = (_vecs...)
const Vecs = Union{vecs...}

##################################################
# Intrinsics that don't have a julia counterpart:

pow{T <: Numbers}(a::T, b::T) = a ^ b
cl_pow{T1, T2}(a::T1, b::T2) = pow(promote(a, b)...)
cl_pow{T <: Numbers}(a::T, b::T) = pow(a, b)

"""
smoothstep performs smooth Hermite interpolation between 0 and 1 when edge0 < x < edge1. This is useful in cases where a threshold function with a smooth transition is desired. smoothstep is equivalent to:
```
    t = clamp((x - edge0) / (edge1 - edge0), 0.0, 1.0);
    return t * t * (3.0 - 2.0 * t);
```
Results are undefined if edge0 ≥ edge1.
"""
function smoothstep{T}(edge0, edge1, x::T)
    t = clamp.((x .- edge0) ./ (edge1 .- edge0), T(0), T(1))
    return t * t * (T(3) - T(2) * t)
end

"""
mix performs a linear interpolation between x and y using a to weight between them. The return value is computed as
`x .* (T(1) .- a) .+ y .* a`
"""
mix{T}(x, y, a::T) = x .* (T(1) .- a) .+ y .* a

fract(x) = x - floor(x)
fabs(x::AbstractFloat) = abs(x)

#######################################
# shared common intrinsic functions

const functions = (
    +, -, *, /, ^, <=, .<=, !, <, >, ==, !=, |, &, %,
    <<, >>,
    sqrt, fract, log,
    round, floor, ceil, trunc,
    sin, sinpi, sinh, asin, asinh,
    cos, cospi, cosh, acos, acosh,
    tan, tanh, atan, atanh, atan2,# atanpi, atan2pi, <- julia doesnt have those?!
    max, min,
    abs, pow, normalize, cross, dot, smoothstep, mix,
    exp2, exp10, expm1,
    log, log2, log10, log1p,
    length, clamp, fma, fabs, isinf, isnan, sign,
    cbrt, copysign, signbit
    # fast math
)

function fixed_array_length(T)
    N = if T <: Tuple
        length(T.parameters)
    else
        length(T)
    end
end
is_ntuple(x) = false
is_ntuple{N, T}(x::Type{NTuple{N, T}}) = true

const vector_lengths = (2, 3, 4, 8, 16)

# This IO thing gets annoying! TODO, remove IO as dispatch object and use exclusively LazyMethods

is_fixedsize_array{T}(::Type{T}) = is_fixedsize_array(nothing, T)
function is_fixedsize_array{T}(m, ::Type{T})
    (T <: StaticVector || is_ntuple(T)) &&
    fixed_array_length(T) in vector_lengths &&
    eltype(T) <: Numbers &&
    eltype(T) != Bool
end

# Functionality to remove unsupported characters from the source
global replace_unsupported, empty_replace_cache!
let _unsupported_id = 0
    const unsupported_replace_dict = Dict{Char, String}()
    function empty_replace_cache!()
        _unsupported_id = 0
        empty!(unsupported_replace_dict)
        return
    end
    """
    Creates a unique replacement for some character
    """
    function replace_unsupported(char::Char)
        get!(unsupported_replace_dict, char) do
            _unsupported_id += 1
            string(_unsupported_id)
        end
    end
end

is_supported_char(io::IO, char) = true
function is_supported_char(io::CIO, char)
    # Lets just assume for simplicity, that only ascii non operators are supported
    # in a name
    isascii(char) &&
    !Base.isoperator(Symbol(char)) &&
    !(char in ('.', '#', '(', ')', ',', '{', '}', ' '))  # some ascii codes are not allowed
end

"""
Makes sure a symbol is well formed in target source
"""
function symbol_hygiene(io::IO, sym)
    # TODO figure out what other things are not allowed
    # TODO startswith gl_, but allow variables that are actually valid inbuilds
    res_io = IOBuffer()
    for (i, char) in enumerate(string(sym))
        res = if is_supported_char(io, char)
            print(res_io, char)
        else
            i == 1 && print(res_io, 'x') # can't start with number
            print(res_io, replace_unsupported(char)) # get a
        end
    end
    String(take!(res_io))
end

typename(io::CIO, x) = Symbol(symbol_hygiene(io, _typename(io, x)))

_typename(io::IO, T::QuoteNode) = _typename(io, T.value)
julia_name(x::Type{Type{T}}) where T = string(T)

function _typename(io::IO, x)
    str = if isa(x, Expr) && x.head == :curly
        string(x, "_", join(x.args, "_"))
    elseif isa(x, Symbol)
        string(x)
    elseif isa(x, UnionAll)
        string(x)
    elseif isa(x, DataType)
        T = x
        if T <: Tuple{X} where X <: Numbers
            # We map Tuple{1, X} to the intrinsic types of it for better performance (usually this would need to be a struct)
            # TODO benchmark if this terrible idea is worth it
            typename(io, eltype(T))
        elseif T <: Type
            # make names unique when it was a type of Type{X}
            string(T)
        elseif is_fixedsize_array(io, T)
            Sugar.vecname(io, T)
        elseif T <: Tuple
            str = if isempty(T.parameters)
                "EmptyTuple_"
            elseif T == Tuple
                "EmptyTuple"
            else
                str = "Tuple_"
                tstr = map(T.parameters) do t
                    if t <: Tuple{X} where X
                        string("Tupl_", typename(io, t))
                    else
                        typename(io, t)
                    end
                end
                str *= join(tstr, "_")
            end
            str
        else
            str = string(T.name.name)
            if !isempty(T.parameters)
                tstr = map(T.parameters) do t
                    if isa(t, DataType)
                        if t <: Tuple{X} where X
                            string("Tupl_", typename(io, t))
                        else
                            typename(io, t)
                        end
                    else
                        string(t)
                    end
                end
                str *= string("_", join(tstr, "_"))
            end
            str
        end
    else
        error("Not transpilable: $x")
    end
    return str
end

_typename(io::CIO, x::Union{AbstractString, Symbol}) = x

_typename{T <: Numbers}(io::IO, x::Type{Tuple{T}}) = _typename(io, T)
_typename(io::IO, x::Type{Float64}) = "double"
_typename(io::IO, x::Type{Float32}) = "float"
_typename(io::IO, x::Type{Int64}) = "long"
_typename(io::IO, x::Type{Int32}) = "int"
_typename(io::IO, x::Type{UInt32}) = "uint"
_typename(io::IO, x::Type{UInt64}) = "ulong"
_typename(io::IO, x::Type{UInt8}) = "uchar"
_typename(io::IO, x::Type{Bool}) = "JLBool"
_typename{T}(io::IO, x::Type{Ptr{T}}) = "$(typename(io, T)) *"

_typename{F <: Function}(io::CIO, f::F) = _typename(io, F.name.mt.name)
_typename{F <: Function}(io::CIO, f::Type{F}) = string(F)

global signature_hash, empty_hash_dict!
let hash_dict = Dict{Any, Int}(), counter = 0
    function empty_hash_dict!()
        empty!(hash_dict)
        counter = 0
        return
    end
    """
    Returns a unique ID for a type signature, which is as small as possible!
    """
    function signature_hash(types)
        get!(hash_dict, Sugar.to_tuple(types)) do
            counter += 1
            counter
        end
    end
end

function Sugar.functionname(io::CIO, method::Type{T}) where T
    string('(', typename(io, T), ')')
end
function Sugar.functionname(io::CIO, method::LazyMethod)
    func = if istype(method)
        # This should only happen, if the function is actually a type constructor
        method.signature
    else
        Sugar.getfunction(method)
    end
    if func == %
        # % seems to be an operator that is printed as rem (?!)
        # TODO, are there more? This is only important for operators that are also intrinsics
        return :(%)
    end
    if func == cli.intrinsic_select
        return :select
    end
    f_sym = if isa(func, Type)
        functionname(io, func)
    else
        Base.function_name(func)
    end
    if Sugar.isintrinsic(method)
        return f_sym # intrinsic operators don't need hygiene!
    end
    str = if isfunction(method) && !supports_overloading(io)
        string(f_sym, '_', signature_hash(method.signature[2]))
    elseif istype(method)
        return string(f_sym)
    end
    if isa(io, Sugar.ASTIO)
        symbol_hygiene(io, str)
    else
        str
    end
end

function show_name(io::CIO, x)
    print(io, symbol_hygiene(io, x))
end
function show_name(io::CIO, x::Union{Slot, SSAValue})
    name = Sugar.slotname(io.method, x)
    show_name(io, name)
end

function show_type(io::CIO, x)
    print(io, typename(io, x))
end

function Base.show_unquoted(io::CIO, slot::Slot, ::Int, ::Int)
    show_name(io, slot)
end

function c_fieldname(m::LazyMethod, T, i::Integer)
    str = if isleaftype(T)
        name = if is_fixedsize_array(m, T)
            fixed_size_array_fieldname(m, T, i)
        else
            Base.fieldname(T, i)
        end
        if isa(name, Integer) # for types without fieldnames (Tuple)
            "field$name"
        else
            symbol_hygiene(EmptyCIO(), name)
        end
    else
        error("Found abstract type: $T")
    end
    Symbol(str)
end

function typed_type_fields(T)
    nf = nfields(T)
    fields = []
    if nf == 0 # structs can't be empty
        # we use bool as a short placeholder type.
        # TODO, are there cases where bool is no good?
        push!(fields, :(emtpy::Float32))
    else
        for i in 1:nf
            FT = fieldtype(T, i)
            tname = typename(EmptyCIO(), FT)
            fname = c_fieldname(T, i)
            push!(fields, :($fname::$tname))
        end
    end
    fields
end
Sugar.show_comment(io::CIO, comment) = println(io, "// ", comment)

show_linenumber(io::CIO, line) = show_comment(io, " line $line:")
show_linenumber(io::CIO, line, file) = show_comment(io, "$file $line $line")

function show_unquoted(io::CIO, x::Type, ::Int, ::Int)
    print(io, "TYP_INST_", typename(io, Type{x}))
end
function show_unquoted(io::CIO, sym::Symbol, ::Int, ::Int)
    print(io, Symbol(symbol_hygiene(io, sym)))
end
function show_unquoted(io::CIO, ssa::SSAValue, ::Int, ::Int)
    print(io, Sugar.ssavalue_name(ssa))
end
function show_unquoted(io::CIO, f::F, ::Int, ::Int) where F <: Sugar.AllFuncs
    print(io, "FUNC_INST_", typename(io, F))
end

const expr_calls_extended = merge(expr_calls, Dict(:constructor => ('{', '}')))
# show a normal (non-operator) function call, e.g. f(x,y) or A[z]
function show_call(io::CIO, head, func, func_args, indent)
    if head == :curly && Sugar.expr_type(func) <: Type# typeconstructors
        print(io, "TYP_INST_") # use the global constant instance of the type
        show_unquoted(io, func, indent)
        return
    end
    op, cl = expr_calls_extended[head]
    # print(io, '(')
    if head == :ref
        FT = Sugar.expr_type(func)
        if is_fixedsize_array(FT) && length(func_args) == 1
            # Special case for small vectors an an integer variable as index
            IDXT = Sugar.expr_type(func_args[1])
            if IDXT <: Integer && !isa(func_args[1], Integer)
                ET = eltype(FT)
                print(io, "(($(typename(io, ET))*)&") # convert to pointer
                show_unquoted(io, func, indent)
                print(io, ")[")
                show_unquoted(io, func_args[1])
                print(io, "]")
                return
            end
        end
        print(io, '(')
        show_unquoted(io, func, indent)
        print(io, ')')
        if FT <: Tuple{T} where T <: cli.Numbers
            # we Tuple{<: Numbers} is treated as scalar, so we don't print the index expression
            return
        end
    else
        print(io, functionname(io, func))
    end
    if head == :(.)
        print(io, '.')
    end
    addition_tracked_object = []
    m = io.method
    for elem in func_args
        T = Sugar.expr_type(m, elem)
        contains, fields = Sugar.contains_tracked_type(m, T)
        if contains && haskey(m.cache, :tracked_types)
            fields = m.cache[:tracked_types][elem]
            for field in fields
                push!(addition_tracked_object, last(field))
            end
        end
    end
    show_enclosed_list(io, op, vcat(func_args, addition_tracked_object), ", ", cl, indent)
    return
end

# show a block, e g if/for/etc
function show_block(io::CIO, head, args::Vector, body, indent::Int)
    if isempty(args) && isa(body, Expr) && isempty(body.args)
        return # everything empty, let's not make a fool of ourselves and print something
    end
    if isempty(args)
        print(io, head, '{')
    else
        print(io, head, '(')
        show_list(io, args, ", ", indent)
        print(io, "){")
    end
    ind = head === :module || head === :baremodule ? indent : indent + indent_width
    exs = (is_expr(body, :block) || is_expr(body, :body)) ? body.args : Any[body]
    for (i, ex) in enumerate(exs)
        sep = i == 1 ? "" : ";"
        print(io, sep, '\n', " "^ind)
        show_unquoted(io, ex, ind, -1)
    end
    print(io, ";\n", " "^indent)
end

function show_unquoted(io::CIO, ex::Expr, indent::Int, prec::Int)
    line_number = 0 # TODO include line numbers
    head, args, nargs = ex.head, ex.args, length(ex.args)
    # dot (i.e. "x.y"), but not compact broadcast exps
    if head == :(.) && !is_expr(args[2], :tuple)
        show_unquoted(io, args[1], indent + indent_width)
        print(io, '.')
        if is_quoted(args[2])
            show_unquoted(io, unquoted(args[2]), indent + indent_width)
        else
            print(io, '(')
            show_unquoted(io, args[2], indent + indent_width)
            print(io, ')')
        end

    # variable declaration TODO might this occure in normal code?
    elseif head == :(::) && nargs == 2
        show_type(io, args[2])
        print(io, ' ')
        show_name(io, args[1])
    # infix (i.e. "x<:y" or "x = y")

    elseif (head in expr_infix_any && (nargs == 2) || (head == :(:)) && nargs == 3)
        func_prec = operator_precedence(head)
        head_ = head in expr_infix_wide ? " $head " : head
        if func_prec <= prec
            show_enclosed_list(io, '(', args, head_, ')', indent, func_prec, true)
        else
            show_list(io, args, head_, indent, func_prec, true)
        end

    # function call
    elseif head === :call && nargs >= 1
        f = first(args)
        func_args = args[2:end]
        fname = Symbol(functionname(io, f))
        if fname == :getfield && nargs == 3
            accessed, fieldname = args[2], args[3]
            m = io.method
            tt = Sugar.expr_type(ex)
            show_unquoted(io, accessed, indent) # type to be accessed
            print(io, '.')
            fieldname = isa(fieldname, QuoteNode) ? fieldname.value : fieldname
            print(io, fieldname)
        else
            func_prec = operator_precedence(fname)
            # scalar multiplication (i.e. "100x")
            if (
                    fname === :* &&
                    length(func_args) == 2 && isa(func_args[1], Real) &&
                    isa(func_args[2], Symbol)
                )
                if func_prec <= prec
                    show_enclosed_list(io, '(', func_args, "", ')', indent, func_prec)
                else
                    show_list(io, func_args, "", indent, func_prec)
                end

            # unary operator (i.e. "!z")
            elseif isa(fname, Symbol) && fname in uni_ops && length(func_args) == 1
                print(io, fname)
                if isa(func_args[1], Expr) || func_args[1] in all_ops
                    show_enclosed_list(io, '(', func_args, ",", ')', indent, func_prec)
                else
                    show_unquoted(io, func_args[1])
                end

            # binary operator (i.e. "x + y")
            elseif func_prec > 0 # is a binary operator
                na = length(func_args)
                if (na == 2 || (na > 2 && fname in (:+, :++, :*))) && all(!isa(a, Expr) || a.head !== :... for a in func_args)
                    sep = " $fname "
                    if func_prec <= prec
                        show_enclosed_list(io, '(', func_args, sep, ')', indent, func_prec, true)
                    else
                        show_list(io, func_args, sep, indent, func_prec, true)
                    end
                elseif na == 1
                    # 1-argument call to normally-binary operator
                    op, cl = expr_calls[head]
                    show_unquoted(io, f, indent)
                    show_enclosed_list(io, op, func_args, ",", cl, indent)
                else
                    show_call(io, head, f, func_args, indent)
                end

            # normal function (i.e. "f(x,y)")
            else
                show_call(io, head, f, func_args, indent)
            end
        end
    # other call-like expressions ("A[1,2]", "T{X,Y}", "f.(X,Y)")
    elseif haskey(expr_calls, head) && nargs >= 1  # :ref/:curly/:calldecl/:(.)
        funcargslike = head == :(.) ? ex.args[2].args : ex.args[2:end]
        show_call(io, head, ex.args[1], funcargslike, indent)
    # comparison (i.e. "x < y < z")
    elseif (head == :comparison) && nargs >= 3 && (nargs & 1==1)
        comp_prec = minimum(operator_precedence, args[2:2:end])
        if comp_prec <= prec
            show_enclosed_list(io, '(', args, " ", ')', indent, comp_prec)
        else
            show_list(io, args, " ", indent, comp_prec)
        end
    # block with argument
    elseif head in (:while, :function, :if) && nargs == 2
        show_block(io, head, args[1], args[2], indent)
        print(io, "}")
    elseif head == :for
        forheader = args[1]
        forheader.head == :(=) || error("Unsupported for: $ex")
        i, range = forheader.args
        range.head == :(:) || error("Unsupported for: $ex")
        from, to = range.args
        print(io, "for(")
        show_unquoted(io, i)
        print(io, " = ")
        show_unquoted(io, from)
        print(io, "; ")
        show_unquoted(io, i)
        print(io, " <= ")
        show_unquoted(io, to)
        print(io, "; ")
        show_unquoted(io, i)
        print(io, "++)")
        show_block(io, "", [], args[2], indent)
        print(io, "}")
    elseif (head == :module) && nargs==3 && isa(args[1],Bool)
        show_block(io, args[1] ? :module : :baremodule, args[2], args[3], indent)
        print(io, "}")

    # type annotation (i.e. "::Int")
    elseif head == :(::) && nargs == 1
        print(io, ' ')
        print(io, typename(io, args[1]))

    elseif (nargs == 0 && head in (:break, :continue))
        print(io, head)

    elseif (nargs == 1 && head in (:abstract, :const)) ||
                          head in (:local,  :global, :export)
        print(io, head, ' ')
        show_list(io, args, ", ", indent)

    elseif (head === :macrocall) && nargs >= 1
        # expand macros
        show_unquoted(io, expand(ex), indent)

    elseif (head === :line) && 1 <= nargs <= 2
        show_linenumber(io, args...)

    elseif (head === :if) && nargs == 3     # if/else
        show_block(io, "if",   args[1], args[2], indent)
        show_block(io, "} else", args[3], indent)
        print(io, "}")

    elseif (head === :block) || (head === :body)
        show_block(io, "", ex, indent); print(io, "}")

    elseif (head === :new)
        constr_args = args[2:end]
        if isempty(constr_args)
            push!(constr_args, 0f0)
        end
        show_call(io, :constructor, args[1], constr_args, indent)
    elseif head === :return
        if length(args) == 1
            if args[1] != nothing
                if !(isa(args[1], Expr) && (args[1].typ == Void))
                    # if returns void, we need to omit the return statement
                    print(io, "return ")
                end
                show_unquoted(io, args[1])
            end
        else
            # ignore if empty, otherwise, LOL? What's a return with multiple args?
            if isempty(args)
                print(io, "return")
            else
                error("Unknown return Expr: $ex")
            end
        end
    elseif (head in (:meta, :inbounds, :boundscheck))
        # TODO, just ignore this? Log this? We definitely don't need it in GLSL
    elseif head == :local_memory_init
        # e.g. __local float test[100]
        T, N, name = args
        print(io, "__local ")
        show_type(io, T)
        print(io, ' ')
        show_name(io, name)
        print(io, "[$N]")
    elseif head == :local_memory
        # take the above initialized name and returns a pointer (local_memory_init & local_memory
        # are guaranteed to be emitted together)
        # and now just return &test;
        T, name = args
        print(io, "( __local ")
        show_type(io, T)
        print(io, " *)(&")
        show_name(io, name)
        print(io, ')')
    else
        println(ex)
        unsupported_expr(string(ex), line_number)
    end
    nothing
end

function show_typed_list(io, list, seperator, intent)
    for (i, elem) in enumerate(list)
        i != 1 && print(io, seperator)
        name, T = elem.args
        print(io, "    "^(intent))
        show_type(io, T)
        print(io, ' ')
        show_name(io, name)
    end
end

function show_returntype(io, method)
    T = Sugar.returntype(method)
    # void can only be used for return types in C, but in Julia it's used for all sorts of things
    # so we special case the return type and otherwise treat Void not as an inbuild.
    if T == Void
        print(io, "void")
    else
        show_type(io, T)
    end
    return
end
function Sugar.getfuncheader!(x::Union{LazyMethod{:CL}, LazyMethod{:GL}})
    if !isdefined(x, :funcheader)
        x.funcheader = if Sugar.isfunction(x)
            prototype = sprint() do io
                args = Sugar.getfuncargs(x)
                glio = CIO(io, x)
                show_returntype(glio, x)
                print(glio, ' ')
                show_unquoted(glio, x)
                print(glio, '(')
                for (i, elem) in enumerate(args)
                    if i != 1
                        print(glio, ", ")
                    end
                    name, T = elem.args
                    show_type(glio, T)
                    print(glio, ' ')
                    show_name(glio, name)
                end
                print(glio, ')')
            end
            # first is the function prototype!
            header = prototype
            prototype * ";\n" * header
        else
            ""
        end
    end
    x.funcheader
end

function Sugar.getfuncsource(x::Union{LazyMethod{:CL}, LazyMethod{:GL}})
    # TODO make this lazy as well?
    sprint() do io
        show_unquoted(CIO(io, x), Sugar.getast!(x), 0, 0)
    end
end

function Sugar.getfuncargs(x::Union{LazyMethod{:CL}, LazyMethod{:GL}})
    functype = x.signature[1]
    calltypes, slots = Sugar.to_tuple(x.signature[2]), Sugar.getslots!(x)
    n = Sugar.method_nargs(x)
    start = ifelse(Sugar.isclosure(functype), 1, 2)
    unpacked_pointers = []
    args = map(start:n) do i
        argtype, name = slots[i]
        # Slot types might be less specific, e.g. when the variable is unused it might end up as Any.
        # but generally the slot type is the correct one, especially in the context of varargs.
        if !isleaftype(argtype) && length(calltypes) <= i
            argtype = calltypes[i - 1]
        end
        if Sugar.contains_tracked_type(x, argtype)[1] && haskey(x.cache, :tracked_types)
            pointers = map(x.cache[:tracked_types][TypedSlot(i, argtype)]) do field
                ptr_typ = Sugar.get_fields_type(argtype, field[1:end-1])
                :($(last(field))::$ptr_typ)
            end
            append!(unpacked_pointers, pointers)
        end
        expr = :($(name)::$(argtype))
        expr.typ = argtype
        expr
    end
    vcat(args, unpacked_pointers)
end

include("rewriting.jl")