text:
# TODO: maybe try 16 bit?
const Priority = UInt64

@enum PriorityType::UInt8 begin
    DEPTH_FIRST_PRIORITY
    CONSTANT_PRIORITY
    RANDOM_PRIORITY
end

struct DepthFirstScheduler <: Scheduler
    multiq::MultiQueue{Priority,Work}
    workers::Vector{Worker}
    state::SchedulerState
    priority_type::PriorityType
end

function DepthFirstScheduler(; priority_type::PriorityType = DEPTH_FIRST_PRIORITY)
    # TODO: Support and test "oversubscription"; useful for waitable tasks.
    nworkers = Threads.nthreads()
    sch = DepthFirstScheduler(
        MultiQueue{Priority,Work}(),
        Worker[],  # to be filled
        SchedulerState(nworkers),
        priority_type,
    )
    for wid in 1:nworkers
        function workerloop()
            try
                workerloop!(sch, wid)
            catch err
                @error "Worker loop exited" exception = (err, catch_backtrace())
                rethrow()
            end
        end
        task = Task(workerloop)
        task.sticky == true
        ccall(:jl_set_task_tid, Cvoid, (Any, Cint), task, (wid - 1) % nworkers)
        push!(sch.workers, Worker(WorkerState(WKR_WORKING), task))
    end
    for wkr in sch.workers
        schedule(wkr.task)
    end
    return sch
end

function workerloop!(sch::DepthFirstScheduler, wid::Integer)
    multiq = sch.multiq
    wkr = sch.workers[wid].state
    schstate = sch.state
    npolls = 0
    while true
        while true
            item = maybepopmin!(multiq)
            if item === nothing
                npolls += 1
                if npolls < 1_000_000
                    GC.safepoint()
                    ccall(:jl_cpu_pause, Cvoid, ())
                    # Main.@tlc loop_poll
                else
                    break
                end
            else
                run!(something(item::Some{Work}))
                GC.safepoint()
                npolls = 0
            end
        end

        @trace(label = :worker_start_sleep, threadid = Threads.threadid())
        @atomic wkr.state = WKR_WAITING
        @atomic schstate.counter -= 1
        wait()
        @assert (@atomic wkr.state) == WKR_NOTIFIED
        @atomic wkr.state = WKR_WORKING
        @trace(label = :worker_done_sleep, threadid = Threads.threadid())
        # Main.@tlc loop_wakeup
    end
end

const CURRENT_PRIORITY = gensym(:_CURRENT_PRIORITY_)

const PriorityRange = Pair{Priority,Priority}

entirerange(::Type{R}) where {T,R<:Pair{T,T}} = (typemin(T) => typemax(T) - oneunit(T))::R

# TODO: Better storage than `task_local_storage`?
get_current_priority_range() =
    get(task_local_storage(), CURRENT_PRIORITY, entirerange(PriorityRange))::PriorityRange
set_current_priority_range(pr::PriorityRange) = task_local_storage(CURRENT_PRIORITY, pr)

# For debugging:
reset_current_priority_range() = set_current_priority_range(entirerange(PriorityRange))

function halverange(pr::Pair)
    f = first(pr)
    l = last(pr)
    f >= l && return (pr, pr)
    m = (l - f) ÷ 2 + f
    return (f => m, m => l)
end

"""
    rollback_priority(f)

Run `f()` and rollback the priority range. As such, `f` must synchronize all
child tasks.
"""
function rollback_priority(f)
    pr = get_current_priority_range()
    try
        f()
    finally
        set_current_priority_range(pr)
    end
end

function spawn!(f::F, sch::DepthFirstScheduler) where {F}
    pr = get_current_priority_range()
    low, high = halverange(pr)
    lb, ub = entirerange(PriorityRange)
    rnd1 = rand(lb:ub)
    rnd2 = rand(lb:ub)
    if sch.priority_type == CONSTANT_PRIORITY
        low = high = pr
    elseif sch.priority_type == RANDOM_PRIORITY
        low = rnd1 => rnd1
        high = rnd2 => rnd2
    end
    set_current_priority_range(low)  # more urgent
    subrange = high::PriorityRange   # less urgent

    # TODO: move this explicitly to `Work`?
    pdf_wrapper() = task_local_storage(f, CURRENT_PRIORITY, subrange)

    work = Work(pdf_wrapper)
    work.scheduler = sch
    if !trypush!(sch.multiq, first(subrange) => work)
        work.isenqueued = false
        run!(work)  # no deadlock; due to the serial-projection property
        return work
    end
    trywakeupall!(sch)
    return work
end

helpself_until!(isdone, ::DepthFirstScheduler) = false

function helpothers_until!(isdone, sch::DepthFirstScheduler)
    multiq = sch.multiq
    npolls = 0
    backoff = 1
    while true
        isdone() && return
        item = maybepopmin!(multiq)
        if item === nothing
            # TODO: DON'T SPIN!!!!
            # yield()
            nspins = rand(1:backoff)
            npolls += nspins
            for _ in 1:nspins
                GC.safepoint()
                ccall(:jl_cpu_pause, Cvoid, ())
            end
            backoff *= 2
            backoff = min(backoff, 1000)
            # Main.@tlc helpothers_poll
            if npolls > 1_000_000_000
                error("timeout")
            end
        else
            npolls = 0
            run!(something(item::Some{Work}))
            GC.safepoint()
            # Main.@tlc helpothers_helped
        end
    end
end

const DEPTH_FIRST_SCHEDULER = Ref{DepthFirstScheduler}()
const CONSTANT_PRIORITY_SCHEDULER = Ref{DepthFirstScheduler}()
const RANDOM_PRIORITY_SCHEDULER = Ref{DepthFirstScheduler}()

function init_depth_first_scheduler()
    DEPTH_FIRST_SCHEDULER[] = DepthFirstScheduler()
    CONSTANT_PRIORITY_SCHEDULER[] = DepthFirstScheduler(; priority_type = CONSTANT_PRIORITY)
    RANDOM_PRIORITY_SCHEDULER[] = DepthFirstScheduler(; priority_type = RANDOM_PRIORITY)
end

function Base.show(io::IO, sch::DepthFirstScheduler)
    if sch === DEPTH_FIRST_SCHEDULER[]
        print(io, TapirSchedulers, '.', "DEPTH_FIRST_SCHEDULER[]")
    elseif sch === CONSTANT_PRIORITY_SCHEDULER[]
        print(io, TapirSchedulers, '.', "CONSTANT_PRIORITY_SCHEDULER[]")
    elseif sch === RANDOM_PRIORITY_SCHEDULER[]
        print(io, TapirSchedulers, '.', "RANDOM_PRIORITY_SCHEDULER[]")
    else
        print(io, DepthFirstScheduler, "()")
    end
end

function Base.show(io::IO, ::MIME"text/plain", sch::DepthFirstScheduler)
    show(io, sch)
    print(io, " ", length(sch.workers), " workers")
    print(io, " (", sch.state.counter, " active worker(s))")
    println(io)
    print(io, "#items in multi queue:")
    for pq in sch.multiq.queues
        print(io, " ", pq.nitems)
    end
end


text:
"""
	solar(cmd0::String="", arg1=[]; kwargs...)

Calculate and plot the day-night terminator and the civil, nautical and astronomical twilights.

Full option list at [`pssolar`](http://gmt.soest.hawaii.edu/doc/latest/pssolar.html)

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_B)
- **C** : **vectors** : -- Str --

    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/pssolar.html#c)
- **G** : **fill** : -- Number or Str --

    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/pssolar.html#g)
- **I** : **inquire** : -- Bool or [] --

	[`-I`](http://gmt.soest.hawaii.edu/doc/latest/pssolar.html#i)
- $(GMT.opt_P)
- **M** : **alpha** : -- Str or [] --

	[`-M`](http://gmt.soest.hawaii.edu/doc/latest/pssolar.html#M)
- $(GMT.opt_R)
- **N** : **radius** : -- Bool or [] --

	[`-N`](http://gmt.soest.hawaii.edu/doc/latest/pssolar.html#n)
- **T** : -- Bool or [] --

	[`-T`](http://gmt.soest.hawaii.edu/doc/latest/pssolar.html#t)
- **W** : **pen** : -- Str or tuple --

	[`-W`](http://gmt.soest.hawaii.edu/doc/latest/pssolar.html#w)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bo)
- $(GMT.opt_h)
- $(GMT.opt_o)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function solar(cmd0::String="", arg1=[]; K=false, O=false, first=true, kwargs...)

	arg2 = []		# May be needed if GMTcpt type is sent in via C
	N_args = isempty_(arg1) ? 0 : 1

	length(kwargs) == 0 && N_args == 0 && isempty(data) && return monolitic("pssolar", cmd0, arg1)	# Speedy mode

    d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

    cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12cd/0d")
	cmd = parse_UVXY(cmd, d)
	cmd, = parse_bo(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_o(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd = parse_params(cmd, d)

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	cmd = add_opt(cmd, 'C', d, [:C :])
    cmd = add_opt(cmd, 'G', d, [:G :fill])
    cmd = add_opt(cmd, 'I', d, [:I :sun])
	cmd = add_opt(cmd, 'M', d, [:M :dump])
	cmd = add_opt(cmd, 'N', d, [:M :invert])
	cmd = add_opt(cmd, 'T', d, [:T :terminators])
	cmd = cmd * opt_pen(d, 'W', [:W :pen])

    opt_extra = ""
    if (occursin( "-I", cmd) || occursin("-I", cmd0))
        output = "";    opt_extra = "-I"
    end
	cmd = finish_PS(d, cmd, output, K, O)
    return finish_PS_module(d, cmd, opt_extra, output, fname_ext, opt_T, K, "pssolar", arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
solar!(cmd0::String="", arg1=[]; K=true, O=true,  first=false, kw...) =
	solar(cmd0, arg1; K=K, O=O,  first=first, kw...)

pssolar  = solar				# Alias
pssolar! = solar!				# Alias

ext: jl


text:
# Returns index corresponding to the name in NamedTupleDist.
# If the name corresponds to a multivariate distribution,
# an array containig all indices of this distribution is returned.
function asindex(vs::NamedTupleShape, name::Symbol)
    name = string(name)
    if name in (rn = repeatednames(vs))
        indxs = findall(x-> x == string(name), rn)
        length(indxs) > 1 ? (return indxs) : (return indxs[1])
    else
        throw(ArgumentError("Symbol :$name found not in $vs."))
    end
end

# Returns index corresponding to the name in NamedTupleDist.
function asindex(vs::NamedTupleShape, name::Expr)
    name = string(name)
    if name in (an = all_active_names(vs))
        return findfirst(x-> x == name, an)[1]
    else
        throw(ArgumentError("Symbol :$name found not in $vs."))
    end
end

# for prior
function asindex(ntd::NamedTupleDist, name::Union{Expr, Symbol})
    return asindex(varshape(ntd), name)
end

# Return the name corresponding to the index as Symbol (for univariate)
# or Expr for (multivariate) distributions
function getname(vs::NamedTupleShape, idx::Integer)
    names = all_active_names(vs)(vs)
    return Meta.parse(names[idx])
end

# Return the name corresponding to the index as String
function getstring(vs::NamedTupleShape, idx::Integer)
    names = all_active_names(vs)
    return names[idx]
end

function getstring(prior::NamedTupleDist, idx::Integer)
    vs = varshape(prior)
    names = all_active_names(vs)
    return names[idx]
end

# Return array of strings with the names of all indices.
# For a multivariate distribution, names for each dimension are created by appending "[i]" to the name.
function allnames(vs::NamedTupleShape)
    subindxs = [collect(1:getproperty(vs, k).len) for k in keys(vs)]
    names = [length(subindxs[i]) > 1 ? subname.(string(keys(vs)[i]), subindxs[i]) : [string(keys(vs)[i])] for i in 1:length(vs)]

    return reduce(vcat, names)
end

function all_active_names(vs::NamedTupleShape)
    subindxs = [collect(1:getproperty(vs, k).len) for k in keys(vs)]
    names = [length(subindxs[i]) > 1 ? subname.(string(keys(vs)[i]), subindxs[i]) : [string(keys(vs)[i])] for i in 1:length(vs) if getproperty(vs, keys(vs)[i]).len>0]

    return reduce(vcat, names)
end

function all_active_names(vs::ScalarShape)
    return "v"
end

function all_active_names(vs::ArrayShape)
    ndims = totalndof(vs)
    return ["v_$i" for i in 1:ndims]
end

# Return array of strings with the names of all indices.
# For a multivariate distribution, the name is repeated for all its dimensions.
function repeatednames(vs::NamedTupleShape)
    names = [[string(k) for i in 1:getproperty(vs, k).len] for k in keys(vs) if getproperty(vs, k).len>0]
    return reduce(vcat, names)
end

# Generate sub-name for one dimension of a multivariate distribution
function subname(name::String, indx::Integer)
    return name*"["*string(indx)*"]"
end

get_fixed_names(vs::AbstractValueShape) = Vector{String}()

function get_fixed_names(vs::NamedTupleShape)
    active_names = all_active_names(vs)
    all_names = allnames(vs)
    fixed_names = [n for n in all_names if !in(n, active_names)]
end


text:
#const WireInStimRegAddr = 0x08;
#const WireInStimRegWord = 0x08;

function update_digital_output(fpga::FPGA,d::DigOut)

    stream = 16

    timestep_us = 1.0e6 / fpga.sampleRate

    #Trigger Paramters
    value = (d.triggerEnabled ? (1 << 7) : 0) + (d.triggerOnLow ? (1 << 6) : 0) + (d.edgeTriggered ? (1 << 5) : 0) + d.triggerSource;
    programStimReg(fpga,stream, d.channel, 0, value)

    #Number of Pulses
    value = (d.negStimFirst ? (1 << 10) : 0) + (d.shapeInt << 8) + (d.numPulses - 1);
    programStimReg(fpga,stream, d.channel, 1, value);

    postTriggerDelay = round(Int32,(d.postTriggerDelay / timestep_us + 0.5))
    firstPhaseDuration = round(Int32,(d.firstPhaseDuration / timestep_us + 0.5))
    refractoryPeriod = round(Int32,(d.refractoryPeriod / timestep_us + 0.5))
    pulseTrainPeriod = round(Int32,(d.pulseTrainPeriod / timestep_us + 0.5))

    eventStartStim = postTriggerDelay;
    eventEndStim = eventStartStim + firstPhaseDuration;
    eventEnd = eventEndStim + refractoryPeriod;

    eventRepeatStim = 0;

    if (d.pulseOrTrain == 1)
        eventRepeatStim = eventStartStim + pulseTrainPeriod;
    else
        eventRepeatStim = 4294967295;
    end

    #EventStartStim (post Trigger Delay) #These times are in units of clock ticks
    programStimReg(fpga, stream, d.channel, 4, eventStartStim);

    #eventEndStim (post Trigger Delay + Phase duration)
    programStimReg(fpga,stream, d.channel, 7, eventEndStim);

    #eventRepeatStim
    programStimReg(fpga,stream, d.channel, 8, eventRepeatStim);

    #event End
    programStimReg(fpga,stream, d.channel, 13, eventEnd);

end

function programStimReg(fpga,stream, channel, reg, value)

    SetWireInValue(fpga,WireInStimRegAddr, (stream << 8) + (channel << 4) + reg)
    SetWireInValue(fpga,WireInStimRegWord, value)
    UpdateWireIns(fpga)
    ActivateTriggerIn(fpga,TrigInRamAddrReset, 1)
end


text:
"""
Plotting spin statistics and Hamiltonian weights. Statistics of local fields and
w_i^k are plotted rather ugly. Contains also functions for plotting
wavefunctions, and learning rates and likelihoods as a function of QBM
Witeration.
"""

function plot_weights(Weights, modelname, n)
    all_weights = []
    for k in fieldnames(WEIGHTS)
        weights_k = transpose(getfield(Weights, k))
        weights_k[abs.(weights_k) .< 1e-6] .= 0
        push!(all_weights, weights_k)
    end
    posneg_col = [:pu_or; (-absmax(Weights), +absmax(Weights))]
    pos_col = [:inferno; (0, +absmax(Weights))]
    check_col(weights) = any(weights .< 0) ? posneg_col : pos_col
    heatmap_plots = []
    aspect_ratio = [repeat([n], 3); repeat([1], 3)]

    for i=1:6
        p = plot( heatmap(all_weights[i], aspect_ratio=aspect_ratio[i],
            framestyle=:none, color=check_col(all_weights[i])[1]),
            size=(600,600),
            clims= check_col(all_weights[i])[2] )
        push!(heatmap_plots, p)
    end

    p_k = plot(heatmap_plots[1:3]..., layout=(1,3))
    p_kk = plot(heatmap_plots[4:6]..., layout=(1,3))

    plot(p_k, p_kk, layout=(2,1))
    savefig("$modelname-$n-weights-all.pdf")
end

function plot_statistics(Statistics, namebase, n)
    all_stats = []
    for k in fieldnames(STATISTICS)
        modelstats = transpose(getfield(Statistics, k))
        modelstats[abs.(modelstats) .< 1e-6] .= 0.0
        push!(all_stats, modelstats)
    end
    posneg_col = [:pu_or; (-1.0, 1.0)]
    pos_col = [:inferno; (0.0, 1.0)]
    check_col(stats) = any(stats .< 0) ? posneg_col : pos_col
    heatmap_plots = []
    aspect_ratio = [repeat([0.5 * n], 3); repeat([1], 3)]

    for i=1:6
        p = plot( heatmap(all_stats[i], aspect_ratio=aspect_ratio[i],
            framestyle=:none, color=check_col(all_stats[i])[1]),
            size=(600,600),
            clims= check_col(all_stats[i])[2] )
        push!(heatmap_plots, p)
    end

    p_k = plot(heatmap_plots[1:3]..., layout=(1,3))
    p_kk = plot(heatmap_plots[4:6]..., layout=(1,3))

    plot(p_k, p_kk, layout=(2,1))
    savefig("$namebase-stats-all.pdf")
end


text:
module PushDirectBenchmarkTest
using PushDirectBenchmark
using Test

@testset "PushDirectBenchmark.jl" begin
    n = 1000
    tc_res1 = run_test_case(1, n)
    for i = 2:get_test_cases_count()
        @test tc_res1 == run_test_case(i, n)
    end
    # @test f1(n) == f2(n)
    # @test f1(n) == f3(n)
    # @test_throws BoundsError f4(n)
    # @test f1(n) != f5(n)
    # @test f1(n) == f6(n)
    # @test f1(n) == f7(n)
end

end # module


text:
struct FINITO_DFlbfgs_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg, TH} <: CIAO_iterable
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                    # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    γ::Maybe{Union{Array{R},R}}  # stepsizes
    sweeping::Int8             # to only use one stepsize γ
    batch::Int                # batch size
    α::R                    # in (0, 1), e.g.: 0.99
    H::TH                   # for lbfgs (lbfgs type)
    inseq::Int              # extra 'memory length' needed = length of inner seq (excluding z)   # exm: 1231231234234234... as memory = 3
    rep::Int                # number of reps for each index   # how many times 123 is repeated
end

mutable struct FINITO_DFlbfgs_state{R<:Real,Tx, TH}
    γ::Array{R}             # stepsize parameter
    hat_γ::R                # average γ
    av::Tx                  # the running average
    ind::Array{Array{Int}}  # running index set
    d::Int                  # number of batches
    H::TH                   # Hessian approx
    z_M::Vector{Tx}         # palceholder for z vectors
    # some extra placeholders
    z::Tx
    ∇f_temp::Tx             # placeholder for gradients
    zbar::Tx                # bar z
    zbar_prev::Maybe{Tx}    # bar z previous
    res_zbar::Tx            # v
    res_zbar_prev::Maybe{Tx}# v
    dir::Tx                 # direction
    ∇f_sum::Tx              # for linesearch
    z_trial::Tx             # linesearch candidate
    inds::Array{Int}        # needed for shuffled only!
    τ::Float64              # number of epochs
end

function FINITO_DFlbfgs_state(γ::Array{R}, hat_γ::R, av::Tx, ind, d, H::TH, z_M) where {R,Tx,TH}
    return FINITO_DFlbfgs_state{R,Tx,TH}(
        γ,
        hat_γ,
        av,
        ind,
        d,
        H,
        z_M,
        copy(av),
        copy(av),
        copy(av),
        nothing,
        copy(av),
        nothing,
        copy(av),
        copy(av),
        copy(av),
        collect(1:d),
        1.0
        )
end

function Base.iterate(iter::FINITO_DFlbfgs_iterable{R}) where {R}
    N = iter.N
    r = iter.batch # batch size
    # create index sets
    ind = Vector{Vector{Int}}(undef, 0)
    d = Int(floor(N / r))
    for i = 1:d
        push!(ind, collect(r*(i-1)+1:i*r))
    end
    r * d < N && push!(ind, collect(r*d+1:N))
    # updating the stepsize
    if iter.γ === nothing
        if iter.L === nothing
            @warn "--> smoothness parameter absent"
            return nothing
        else
            γ = zeros(R, N)
            for i = 1:N
                isa(iter.L, R) ? (γ = fill(iter.α * R(iter.N) / iter.L, (N,))) :
                (γ[i] = iter.α * R(N) / (iter.L[i]))
            end
        end
    else
        isa(iter.γ, R) ? (γ = fill(iter.γ, (N,))) : (γ = iter.γ) # provided γ
    end
    #initializing the vectors
    hat_γ = 1 / sum(1 ./ γ)
    av = copy(iter.x0)
    for i = 1:N
        ∇f, ~ = gradient(iter.F[i], iter.x0)
        ∇f .*= hat_γ / N
        av .-= ∇f
    end
    z_M = [copy(av) for i = 1:(iter.inseq+1)]

    state = FINITO_DFlbfgs_state(γ, hat_γ, av, ind, cld(N, r), iter.H, z_M)

    return state, state
end

function Base.iterate(
    iter::FINITO_DFlbfgs_iterable{R},
    state::FINITO_DFlbfgs_state{R},
) where {R}

    if state.zbar_prev === nothing
        state.zbar_prev = zero(state.z)
        state.res_zbar_prev = zero(state.z)
    end

    # full update
    state.zbar, ~ = prox(iter.g, state.av, state.hat_γ)
    envVal = 0.0
    state.av .= state.zbar
    state.∇f_sum .= zero(state.av)
    for i = 1:iter.N
        state.∇f_temp, fi_z = gradient(iter.F[i], state.zbar) # update the gradient
        state.av .-= (state.hat_γ / iter.N) .* state.∇f_temp
        envVal += fi_z / iter.N
        state.∇f_sum .+= state.∇f_temp
    end
    state.res_zbar, gz = prox(iter.g, state.av, state.hat_γ)
    envVal += gz
    # state.count += 1

    state.res_zbar .-= state.zbar # \bar v- \bar z

    if state.zbar_prev !== nothing # for lbfgs
        # update metric ##### bug prone  v = -res not res
        update!(state.H, state.zbar - state.zbar_prev, -state.res_zbar +  state.res_zbar_prev)
        # store vectors for next update
        copyto!(state.zbar_prev, state.zbar)
        copyto!(state.res_zbar_prev, state.res_zbar)
    end

    mul!(state.dir, state.H, state.res_zbar)

    envVal += real(dot(state.∇f_sum, state.res_zbar)) / iter.N
    envVal += norm(state.res_zbar)^2 / (2 *  state.hat_γ)

    # Compute_direction!(iter,state)
    # state.dir .=  state.res_zbar # sanity check

    # println("envelope(bar z) outside is $(envVal)")
    state.τ = 1.0
    # while true  ######## z_line should be removed using some other placeholder
    for i=1:5
        state.z_trial .=  state.zbar .+ (1- state.τ) .* state.res_zbar + state.τ * state.dir

        # compute varphi(z_trial)
        state.av .= state.z_trial
        state.∇f_sum = zero(state.av)
        envVal_trial = 0
        for i = 1:iter.N
            state.∇f_temp, fi_z = gradient(iter.F[i], state.z_trial) # update the gradient
            state.av .-= (state.hat_γ / iter.N) .* state.∇f_temp
            envVal_trial += fi_z / iter.N
            state.∇f_sum .+= state.∇f_temp
        end
        # state.count += 1
        state.z, gz = prox(iter.g, state.av, state.hat_γ)
        envVal_trial += gz
        state.z .-= state.z_trial # \bar v- \bar z

        envVal_trial += real(dot(state.∇f_sum, state.z)) / iter.N
        envVal_trial += norm(state.z)^2 / (2 *  state.hat_γ)

        envVal_trial <= envVal + eps(R) && break
        # println("ls backtracked, tau was $(state.τ)")
        state.τ /= 50       ##### bug prone: change in reporting if you change this! ######
    end
    state.zbar .= state.z_trial # of line search

    ###### at the momoent I am computing the next prox twice (it is already computed in the ls)
    iter.sweeping == 3 && (state.inds = randperm(state.d)) # shuffled

    # println("rep is $(iter.rep)")
    Ind = Index_iterable(1, iter.inseq, iter.rep, iter.N) # in CIAO utilities

    for (it, ind) in Iterators.enumerate(Ind) # here we don't have batches
        # println(it)
        if mod(it, iter.inseq * iter.rep) == 1
            # println("------new cycle starts-----")
            for ell in 1:(iter.inseq+1)       # for initializing the inner loop
                copyto!(state.z_M[ell], state.zbar) # resetting with zbar: x^full
            end
        end
        println(ind)
        prox!(state.z_M[ind.m_n], iter.g, state.av, state.hat_γ) # always m_o is 1 larger than m_n

        gradient!(state.∇f_temp, iter.F[state.inds[ind.i]], state.z_M[ind.m_o]) # update the gradient - ind.i counts 123123123456...
        state.av .+= (state.hat_γ / iter.N) .* state.∇f_temp
        gradient!(state.∇f_temp, iter.F[state.inds[ind.i]], state.z_M[ind.m_n]) # update the gradient
        state.av .-= (state.hat_γ / iter.N) .* state.∇f_temp
        state.av .+= (state.hat_γ / state.γ[state.inds[ind.i]]) .* (state.z_M[ind.m_n] .- state.z_M[ind.m_o])
    end

    return state, state
end

solution(state::FINITO_DFlbfgs_state) = state.z_M[end]
epoch_count(state::FINITO_DFlbfgs_state) = state.τ   # number of epochs is 2+ 1/tau , where 1/tau is from ls

###### there is some bug when rep inseq =(1,1) it is not working......check after adding adaptive stepsizes!!!!


text:
D2hs = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("sigmah", [1 0 0; 0 1 0; 0 0 -1]),
    Symmetry.Symel("i", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2^1", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("sigmav_1", [-1.0 0.0 -0.0; 0.0 1.0 0.0; -0.0 0.0 1.0]),
    Symmetry.Symel("sigmad_1", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])]

D3hs = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("sigmah", [1 0 0; 0 1 0; 0 0 -1]),
    Symmetry.Symel("C_3^1", [-0.4999999999999998 -0.8660254037844387 0.0; 0.8660254037844387 -0.4999999999999998 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^2", [-0.5000000000000003 0.8660254037844384 0.0; -0.8660254037844384 -0.5000000000000003 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("S_3^1", [-0.4999999999999998 -0.8660254037844387 -0.0; 0.8660254037844387 -0.4999999999999998 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_3^5", [-0.5000000000000003 0.8660254037844384 0.0; -0.8660254037844384 -0.5000000000000003 -0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("sigmav_2", [0.5000000000000002 0.8660254037844387 0.0; 0.8660254037844387 -0.5000000000000009 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmav_3", [0.4999999999999991 -0.8660254037844394 0.0; -0.8660254037844394 -0.4999999999999998 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("sigmav_1", [-1.0 0.0 -0.0; 0.0 1.0 0.0; -0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-0.5000000000000002 -0.8660254037844387 0.0; -0.8660254037844387 0.5000000000000009 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(3)", [-0.4999999999999991 0.8660254037844394 0.0; 0.8660254037844394 0.4999999999999998 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])]

D4hs = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("sigmah", [1 0 0; 0 1 0; 0 0 -1]),
    Symmetry.Symel("i", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_4^1", [0.0 -1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2^1", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_4^3", [0.0 1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("S_4^1", [0.0 -1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_4^3", [0.0 1.0 0.0; -1.0 0.0 -0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("sigmav_2", [-1.0 0.0 0.0; 0.0 1.0 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmav_1", [-1.0 0.0 -0.0; 0.0 1.0 0.0; -0.0 0.0 1.0]),
    Symmetry.Symel("sigmad_2", [0.0 1.0 0.0; 1.0 0.0 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_1", [0.0 -1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(3)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(4)", [-1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])]

D5hs = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("sigmah", [1 0 0; 0 1 0; 0 0 -1]),
    Symmetry.Symel("C_5^1", [0.30901699437494745 -0.9510565162951535 0.0; 0.9510565162951535 0.30901699437494745 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^2", [-0.8090169943749473 -0.5877852522924731 0.0; 0.5877852522924731 -0.8090169943749473 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^3", [-0.8090169943749475 0.587785252292473 0.0; -0.587785252292473 -0.8090169943749475 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^4", [0.30901699437494734 0.9510565162951535 0.0; -0.9510565162951535 0.30901699437494734 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("S_5^1", [0.30901699437494745 -0.9510565162951535 0.0; 0.9510565162951535 0.30901699437494745 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_5^7", [-0.8090169943749473 -0.5877852522924731 -0.0; 0.5877852522924731 -0.8090169943749473 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_5^3", [-0.8090169943749475 0.587785252292473 0.0; -0.587785252292473 -0.8090169943749475 -0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_5^9", [0.30901699437494734 0.9510565162951535 0.0; -0.9510565162951535 0.30901699437494734 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("sigmav_2", [0.8090169943749473 -0.5877852522924734 -0.0; -0.5877852522924734 -0.8090169943749481 -0.0; -0.0 -0.0 1.0]),
    Symmetry.Symel("sigmav_3", [-0.3090169943749479 0.951056516295154 0.0; 0.951056516295154 0.3090169943749471 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmav_4", [-0.30901699437494834 -0.9510565162951539 0.0; -0.9510565162951539 0.30901699437494745 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("sigmav_5", [0.8090169943749475 0.5877852522924731 -0.0; 0.5877852522924731 -0.8090169943749481 0.0; -0.0 0.0 1.0]),
    Symmetry.Symel("sigmav_1", [-1.0 0.0 -0.0; 0.0 1.0 0.0; -0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-0.8090169943749473 0.5877852522924734 0.0; 0.5877852522924734 0.8090169943749481 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(3)", [0.3090169943749479 -0.951056516295154 0.0; -0.951056516295154 -0.3090169943749471 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(4)", [0.30901699437494834 0.9510565162951539 0.0; 0.9510565162951539 -0.30901699437494745 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(5)", [-0.8090169943749475 -0.5877852522924731 0.0; -0.5877852522924731 0.8090169943749481 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])]

D6hs = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("sigmah", [1 0 0; 0 1 0; 0 0 -1]),
    Symmetry.Symel("i", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_6^1", [0.5000000000000001 -0.8660254037844386 0.0; 0.8660254037844386 0.5000000000000001 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^1", [-0.4999999999999998 -0.8660254037844388 0.0; 0.8660254037844388 -0.4999999999999998 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2^1", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^2", [-0.5000000000000006 0.8660254037844385 0.0; -0.8660254037844385 -0.5000000000000006 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_6^5", [0.49999999999999944 0.8660254037844392 0.0; -0.8660254037844392 0.49999999999999944 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("S_6^1", [0.5000000000000001 -0.8660254037844386 0.0; 0.8660254037844386 0.5000000000000001 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_3^1", [-0.4999999999999998 -0.8660254037844388 -0.0; 0.8660254037844388 -0.4999999999999998 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_3^8", [-0.5000000000000006 0.8660254037844385 0.0; -0.8660254037844385 -0.5000000000000006 -0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_6^5", [0.49999999999999944 0.8660254037844392 0.0; -0.8660254037844392 0.49999999999999944 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("sigmav_2", [0.5000000000000002 0.8660254037844387 0.0; 0.8660254037844387 -0.5000000000000009 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmav_3", [0.4999999999999991 -0.8660254037844394 0.0; -0.8660254037844394 -0.4999999999999998 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("sigmav_1", [-1.0 0.0 -0.0; 0.0 1.0 0.0; -0.0 0.0 1.0]),
    Symmetry.Symel("sigmad_2", [1.0 0.0 -0.0; 0.0 -1.0 -0.0; -0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_3", [-0.4999999999999998 0.8660254037844392 0.0; 0.8660254037844392 0.49999999999999933 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_1", [-0.5000000000000009 -0.8660254037844383 0.0; -0.8660254037844383 0.5000000000000008 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-0.4999999999999998 0.8660254037844388 0.0; 0.8660254037844388 0.4999999999999998 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(3)", [-0.5000000000000004 -0.8660254037844385 0.0; -0.8660254037844385 0.5000000000000007 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(4)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(5)", [-0.4999999999999989 0.8660254037844395 0.0; 0.8660254037844395 0.49999999999999933 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(6)", [-0.5000000000000013 -0.8660254037844378 0.0; -0.8660254037844378 0.5000000000000009 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])]

D2hcn = ["Ag", "B1g", "B2g", "B3g", "Au", "B1u", "B2u", "B3u"]
D2hct = [1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0;
         1.0  1.0 -1.0 -1.0  1.0  1.0 -1.0 -1.0;
         1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0 -1.0;
         1.0 -1.0 -1.0  1.0  1.0 -1.0 -1.0  1.0;
         1.0  1.0  1.0  1.0 -1.0 -1.0 -1.0 -1.0;
         1.0  1.0 -1.0 -1.0 -1.0 -1.0  1.0  1.0;
         1.0 -1.0  1.0 -1.0 -1.0  1.0 -1.0  1.0;
         1.0 -1.0 -1.0  1.0 -1.0  1.0  1.0 -1.0]
D3hcn = ["A1′", "A2′", "E′", "A1″", "A2″", "E″"]
D3hct = D3dct
D4hcn = ["A1g", "A2g", "B1g", "B2g", "Eg", "A1u", "A2u", "B1u", "B2u", "Eu"]
D4hct = [1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0;
         1.0  1.0  1.0 -1.0 -1.0  1.0  1.0  1.0 -1.0 -1.0;
         1.0 -1.0  1.0  1.0 -1.0  1.0 -1.0  1.0  1.0 -1.0;
         1.0 -1.0  1.0 -1.0  1.0  1.0 -1.0  1.0 -1.0  1.0;
         2.0  0.0 -2.0  0.0  0.0  2.0  0.0 -2.0  0.0  0.0;
         1.0  1.0  1.0  1.0  1.0 -1.0 -1.0 -1.0 -1.0 -1.0;
         1.0  1.0  1.0 -1.0 -1.0 -1.0 -1.0 -1.0  1.0  1.0;
         1.0 -1.0  1.0  1.0 -1.0 -1.0  1.0 -1.0 -1.0  1.0;
         1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0 -1.0;
         2.0  0.0 -2.0  0.0  0.0 -2.0  0.0  2.0  0.0  0.0]
D5hcn = ["A1′", "A2′", "E1′", "E2′", "A1″", "A2″", "E1″", "E2″"]
D5hct = D5dct
D6hcn = ["A1g", "A2g", "B1g", "B2g", "E1g", "E2g", "A1u", "A2u", "B1u", "B2u", "E1u", "E2u"]
D6hct = [1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0;
         1.0  1.0  1.0  1.0 -1.0 -1.0  1.0  1.0  1.0  1.0 -1.0 -1.0;
         1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0 -1.0;
         1.0 -1.0  1.0 -1.0 -1.0  1.0  1.0 -1.0  1.0 -1.0 -1.0  1.0;
         2.0  1.0 -1.0 -2.0  0.0  0.0  2.0  1.0 -1.0 -2.0  0.0  0.0;
         2.0 -1.0 -1.0  2.0  0.0  0.0  2.0 -1.0 -1.0  2.0  0.0  0.0;
         1.0  1.0  1.0  1.0  1.0  1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0;
         1.0  1.0  1.0  1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0  1.0  1.0;
         1.0 -1.0  1.0 -1.0  1.0 -1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0;
         1.0 -1.0  1.0 -1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0  1.0 -1.0;
         2.0  1.0 -1.0 -2.0  0.0  0.0 -2.0 -1.0  1.0  2.0  0.0  0.0;
         2.0 -1.0 -1.0  2.0  0.0  0.0 -2.0  1.0  1.0 -2.0  0.0  0.0]


text:
module Memcache

using Sockets
using Serialization

import Base: close, touch, get, replace

export set, cas, add, replace, append, prepend, get, touch, incr, decr, delete
export MemcacheClient, stats, version, flush_all, close, slabs_reassign, slabs_automove, quit
export MemcacheClients
# not exporting add_client as it could unintentionally change hash distribution calculation in the middle of execution
#export add_client

include("client.jl")
include("multi.jl")

end # module


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	generate_model(setup::Dict,inputs::Dict,OPTIMIZER::MOI.OptimizerWithAttributes,modeloutput = nothing)

This function sets up and solves a constrained optimization model of electricity system capacity expansion and operation problem and extracts solution variables for later processing.

In addition to calling a number of other modules to create constraints for specific resources, policies, and transmission assets, this function initializes two key expressions that are successively expanded in each of the resource-specific modules: (1) the objective function; and (2) the zonal power balance expression. These two expressions are the only expressions which link together individual modules (e.g. resources, transmission assets, policies), which otherwise are self-contained in defining relevant variables, expressions, and constraints.

**Objective Function**

The objective function of GenX minimizes total annual electricity system costs over the following six components shown in the equation below:

```math
\begin{aligned}
	&\sum_{y \in \mathcal{G} } \sum_{z \in \mathcal{Z}}
	\left( (\pi^{INVEST}_{y,z} \times \overline{\Omega}^{size}_{y,z} \times  \Omega_{y,z})
	+ (\pi^{FOM}_{y,z} \times \overline{\Omega}^{size}_{y,z} \times  \Delta^{total}_{y,z})\right) + \notag \\
	&\sum_{y \in \mathcal{O} } \sum_{z \in \mathcal{Z}}
	\left( (\pi^{INVEST,energy}_{y,z} \times    \Omega^{energy}_{y,z})
	+ (\pi^{FOM,energy}_{y,z} \times  \Delta^{total,energy}_{y,z})\right) + \notag \\
	&\sum_{y \in \mathcal{O}^{asym} } \sum_{z \in \mathcal{Z}}
	\left( (\pi^{INVEST,charge}_{y,z} \times    \Omega^{charge}_{y,z})
	+ (\pi^{FOM,charge}_{y,z} \times  \Delta^{total,charge}_{y,z})\right) + \notag \\
	& \sum_{y \in \mathcal{G} } \sum_{z \in \mathcal{Z}} \sum_{t \in \mathcal{T}} \left( \omega_{t}\times(\pi^{VOM}_{y,z} + \pi^{FUEL}_{y,z})\times \Theta_{y,z,t}\right) + \sum_{y \in \mathcal{O \cup DF} } \sum_{z \in \mathcal{Z}} \sum_{t \in \mathcal{T}} \left( \omega_{t}\times\pi^{VOM,charge}_{y,z} \times \Pi_{y,z,t}\right) +\notag \\
	&\sum_{s \in \mathcal{S} } \sum_{z \in \mathcal{Z}} \sum_{t \in \mathcal{T}}\left(\omega_{t} \times n_{s}^{slope} \times \Lambda_{s,z,t}\right) + \sum_{t \in \mathcal{T} } \left(\omega_{t} \times \pi^{unmet}_{rsv} \times r^{unmet}_{t}\right) \notag \\
	&\sum_{y \in \mathcal{H} } \sum_{z \in \mathcal{Z}} \sum_{t \in \mathcal{T}}\left(\omega_{t} \times \pi^{START}_{y,z} \times \chi_{s,z,t}\right) + \notag \\
	& \sum_{l \in \mathcal{L}}\left(\pi^{TCAP}_{l} \times \bigtriangleup\varphi^{max}_{l}\right)
\end{aligned}
```

The first summation represents the fixed costs of generation/discharge over all zones and technologies, which refects the sum of the annualized capital cost, $\pi^{INVEST}_{y,z}$, times the total new capacity added (if any),  plus the Fixed O&M cost, $\pi^{FOM}_{y,z}$, times the net installed generation capacity, $\overline{\Omega}^{size}_{y,z} \times \Delta^{total}_{y,z}$ (e.g., existing capacity less retirements plus additions).

The second summation corresponds to the fixed cost of installed energy storage capacity and is summed over only the storage resources. This term includes the sum of the annualized energy capital cost, $\pi^{INVEST,energy}_{y,z}$, times the total new energy capacity added (if any), plus the Fixed O&M cost, $\pi^{FOM, energy}_{y,z}$, times the net installed energy storage capacity, $\Delta^{total}_{y,z}$ (e.g., existing capacity less retirements plus additions).

The third summation corresponds to the fixed cost of installed charging power capacity and is summed over only over storage resources with independent/asymmetric charge and discharge power components ($\mathcal{O}^{asym}$). This term includes the sum of the annualized charging power capital cost, $\pi^{INVEST,charge}_{y,z}$, times the total new charging power capacity added (if any), plus the Fixed O&M cost, $\pi^{FOM, energy}_{y,z}$, times the net installed charging power capacity, $\Delta^{total}_{y,z}$ (e.g., existing capacity less retirements plus additions).

The fourth and fifth summations corresponds to the operational cost across all zones, technologies, and time steps. The fourth summation represents the sum of fuel cost, $\pi^{FUEL}_{y,z}$ (if any), plus variable O&M cost, $\pi^{VOM}_{y,z}$ times the energy generation/discharge by generation or storage resources (or demand satisfied via flexible demand resources, $y\in\mathcal{DF}$) in time step $t$, $\Theta_{y,z,t}$, and the weight of each time step $t$, $\omega_t$, where $\omega_t$ is equal to 1 when modeling grid operations over the entire year (8760 hours), but otherwise is equal to the number of hours in the year represented by the representative time step, $t$ such that the sum of $\omega_t \forall t \in T = 8760$, approximating annual operating costs. The fifth summation represents the variable charging O&M cost, $\pi^{VOM,charge}_{y,z}$ times the energy withdrawn for charging by storage resources (or demand deferred by flexible demand resources) in time step $t$ , $\Pi_{y,z,t}$ and the annual weight of time step $t$,$\omega_t$.

The sixth summation represents the total cost of unserved demand across all segments $s$ of a segment-wise price-elastic demand curve, equal to the marginal value of consumption (or cost of non-served energy), $n_{s}^{slope}$, times the
amount of non-served energy, $\Lambda_{y,z,t}$, for each segment on each zone during each time step (weighted by $\omega_t$).

The seventh summation represents the total cost of not meeting hourly operating reserve requirements, where $\pi^{unmet}_{rsv}$ is the cost penalty per unit of non-served reserve requirement, and $r^{unmet}_t$ is the amount of non-served reserve requirement in each time step (weighted by $\omega_t$).

The eighth summation corresponds to the startup costs incurred by technologies to which unit commitment decisions apply (e.g. $y \in \mathcal{UC}$), equal to the cost of start-up, $\pi^{START}_{y,z}$, times the number of startup events, $\chi_{y,z,t}$, for the cluster of units in each zone and time step (weighted by $\omega_t$).

The last term corresponds to the transmission reinforcement or construction costs, for each transmission line in the model. Transmission reinforcement costs are equal to the sum across all lines of the product between the transmission reinforcement/construction cost, $pi^{TCAP}_{l}$, times the additional transmission capacity variable, $\bigtriangleup\varphi^{max}_{l}$. Note that fixed O\&M and replacement capital costs (depreciation) for existing transmission capacity is treated as a sunk cost and not included explicitly in the GenX objective function.

In summary, the objective function can be understood as the minimization of costs associated with five sets of different decisions: (1) where and how to invest on capacity, (2) how to dispatch or operate that capacity, (3) which consumer demand segments to serve or curtail, (4) how to cycle and commit thermal units subject to unit commitment decisions, (5) and where and how to invest in additional transmission network capacity to increase power transfer capacity between zones. Note however that each of these components are considered jointly and the optimization is performed over the whole problem at once as a monolithic co-optimization problem.

**Power Balance**

The power balance constraint of the model ensures that electricity demand is met at every time step in each zone. As shown in the constraint, electricity demand, $D_{t,z}$, at each time step and for each zone must be strictly equal to the sum of generation, $\Theta_{y,z,t}$, from thermal technologies ($\mathcal{H}$), curtailable VRE ($\mathcal{VRE}$), must-run resources ($\mathcal{MR}$), and hydro resources ($\mathcal{W}$). At the same time, energy storage devices ($\mathcal{O}$) can discharge energy, $\Theta_{y,z,t}$ to help satisfy demand, while when these devices are charging, $\Pi_{y,z,t}$, they increase demand. For the case of flexible demand resources ($\mathcal{DF}$), delaying demand (equivalent to charging virtual storage), $\Pi_{y,z,t}$, decreases demand while satisfying delayed demand (equivalent to discharging virtual demand), $\Theta_{y,z,t}$, increases demand. Price-responsive demand curtailment, $\Lambda_{s,z,t}$, also reduces demand. Finally, power flows, $\Phi_{l,t}$, on each line $l$ into or out of a zone (defined by the network map $\varphi^{map}_{l,z}$), are considered in the demand balance equation for each zone. By definition, power flows leaving their reference zone are positive, thus the minus sign in the below constraint. At the same time losses due to power flows increase demand, and one-half of losses across a line linking two zones are attributed to each connected zone. The losses function $\beta_{l,t}(\cdot)$ will depend on the configuration used to model losses (see Transmission section).

```math
\begin{aligned}
	& \sum_{y\in \mathcal{H}}{\Theta_{y,z,t}} +\sum_{y\in \mathcal{VRE}}{\Theta_{y,z,t}} +\sum_{y\in \mathcal{MR}}{\Theta_{y,z,t}} + \sum_{y\in \mathcal{O}}{(\Theta_{y,z,t}-\Pi_{y,z,t})} + \notag\\
	& \sum_{y\in \mathcal{DF}}{(-\Theta_{y,z,t}+\Pi_{y,z,t})} +\sum_{y\in \mathcal{W}}{\Theta_{y,z,t}}+ \notag\\
	&+ \sum_{s\in \mathcal{S}}{\Lambda_{s,z,t}}  - \sum_{l\in \mathcal{L}}{(\varphi^{map}_{l,z} \times \Phi_{l,t})} -\frac{1}{2} \sum_{l\in \mathcal{L}}{(\varphi^{map}_{l,z} \times \beta_{l,t}(\cdot))} = D_{z,t}
	\forall z\in \mathcal{Z},  t \in \mathcal{T}
\end{aligned}
```

"""

## generate_model(setup::Dict,inputs::Dict,OPTIMIZER::MOI.OptimizerWithAttributes,modeloutput = nothing)
################################################################################
##
## description: Sets up and solves constrained optimization model of electricity
## system capacity expansion and operation problem and extracts solution variables
## for later processing
##
## returns: Model EP object containing the entire optimization problem model to be solved by SolveModel.jl
##
################################################################################
function generate_model(setup::Dict,inputs::Dict,OPTIMIZER::MOI.OptimizerWithAttributes,modeloutput = nothing)

	println("Generating model")

	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones

	## Start pre-solve timer
	presolver_start_time = time()

	# Generate Energy Portfolio (EP) Model
	EP = Model(OPTIMIZER)

	# Introduce dummy variable fixed to zero to ensure that expressions like eTotalCap,
	# eTotalCapCharge, eTotalCapEnergy and eAvail_Trans_Cap all have a JuMP variable
	@variable(EP, vZERO == 0);

	# Initialize Power Balance Expression
	# Expression for "baseline" power balance constraint
	@expression(EP, ePowerBalance[t=1:T, z=1:Z], 0)

	# Initialize Objective Function Expression
	@expression(EP, eObj, 0)

	# Initialize Capacity Reserve Margin Expression
	if setup["CapacityReserveMargin"] > 0
		@expression(EP, eCapResMarBalance[res=1:inputs["NCapacityReserveMargin"], t=1:T], 0)
	end

	# Energy Share Requirement
	if setup["EnergyShareRequirement"] >= 1
		@expression(EP, eESR[ESR=1:inputs["nESR"]], 0)
	end

	if (setup["MinCapReq"] == 1)
		@expression(EP, eMinCapRes[mincap = 1:inputs["NumberOfMinCapReqs"]], 0)
	end

	#@expression(EP, :eCO2Cap[cap=1:inputs["NCO2Cap"]], 0)
	@expression(EP, eGenerationByZone[z=1:Z, t=1:T], 0)

	# Infrastructure
	EP = discharge(EP, inputs, setup["EnergyShareRequirement"])

	EP = non_served_energy(EP, inputs, setup["CapacityReserveMargin"])

	EP = investment_discharge(EP, inputs, setup["MinCapReq"])

	if setup["UCommit"] > 0
		EP = ucommit(EP, inputs, setup["UCommit"])
	end

	EP = emissions(EP, inputs)

	if setup["Reserves"] > 0
		EP = reserves(EP, inputs, setup["UCommit"])
	end

	if Z > 1
		EP = transmission(EP, inputs, setup["UCommit"], setup["NetworkExpansion"], setup["CapacityReserveMargin"])
	end

	# Technologies
	# Model constraints, variables, expression related to dispatchable renewable resources

	if !isempty(inputs["VRE"])
		EP = curtailable_variable_renewable(EP, inputs, setup["Reserves"], setup["CapacityReserveMargin"])
	end

	# Model constraints, variables, expression related to non-dispatchable renewable resources
	if !isempty(inputs["MUST_RUN"])
		EP = must_run(EP, inputs, setup["CapacityReserveMargin"])
	end

	# Model constraints, variables, expression related to energy storage modeling
	if !isempty(inputs["STOR_ALL"])
		EP = storage(EP, inputs, setup["Reserves"], setup["OperationWrapping"], setup["LongDurationStorage"], setup["EnergyShareRequirement"], setup["CapacityReserveMargin"], setup["StorageLosses"])
	end

	# Model constraints, variables, expression related to reservoir hydropower resources
	if !isempty(inputs["HYDRO_RES"])
		EP = hydro_res(EP, inputs, setup["Reserves"], setup["CapacityReserveMargin"])
	end

	# Model constraints, variables, expression related to reservoir hydropower resources with long duration storage
	if setup["OperationWrapping"] == 1 && !isempty(inputs["STOR_HYDRO_LONG_DURATION"])
		EP = hydro_inter_period_linkage(EP, inputs)
	end

	# Model constraints, variables, expression related to demand flexibility resources
	if !isempty(inputs["FLEX"])
		EP = flexible_demand(EP, inputs, setup["CapacityReserveMargin"])
	end
	# Model constraints, variables, expression related to thermal resource technologies
	if !isempty(inputs["THERM_ALL"])
		EP = thermal(EP, inputs, setup["UCommit"], setup["Reserves"], setup["CapacityReserveMargin"])
	end

	# Policies
	# CO2 emissions limits
	EP = co2_cap(EP, inputs, setup)

	# Energy Share Requirement
	if setup["EnergyShareRequirement"] >= 1
		EP = energy_share_requirement(EP, inputs, setup)
	end

	#Capacity Reserve Margin
	if setup["CapacityReserveMargin"] > 0
		EP = cap_reserve_margin(EP, inputs, setup)
	end

	if (setup["MinCapReq"] == 1)
		EP = minimum_capacity_requirement(EP, inputs, setup)
	end

	## Define the objective function
	@objective(EP,Min,EP[:eObj])

	## Power balance constraints
	# demand = generation + storage discharge - storage charge - demand deferral + deferred demand satisfaction - demand curtailment (NSE)
	#          + incoming power flows - outgoing power flows - flow losses - charge of heat storage + generation from NACC
	@constraint(EP, cPowerBalance[t=1:T, z=1:Z], EP[:ePowerBalance][t,z] == inputs["pD"][t,z])

	## Record pre-solver time
	presolver_time = time() - presolver_start_time
    	#### Question - What do we do with this time now that we've split this function into 2?
	if setup["PrintModel"] == 1
		if modeloutput === nothing
			filepath = joinpath(pwd(), "YourModel.lp")
			JuMP.write_to_file(EP, filepath)
		else
			filepath = joinpath(modeloutput, "YourModel.lp")
			JuMP.write_to_file(EP, filepath)
		end
		println("Model Printed")
    end

    return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using BroadcastableStructs
using Test

bcapp(f, args...) = f.(args...)

struct WeightedAdd{A, B} <: BroadcastableCallable
    a::A
    b::B
end

(f::WeightedAdd)(x, y) = f.a * x + f.b * y

struct ManyAdd{T1, T2, T3, T4, T5} <: BroadcastableCallable
    x1::T1
    x2::T2
    x3::T3
    x4::T4
    x5::T5
end

ManyAdd() = ManyAdd(1:5...)

(f::ManyAdd)(y1, y2, y3, y4, y5) =
    f.x1 * y1 + f.x2 * y2 + f.x3 * y3 + f.x4 * y4 + f.x5 * y5

struct Poly3{T0, T1, T2, T3} <: BroadcastableCallable
    c0::T0
    c1::T1
    c2::T2
    c3::T3
end

poly(c0=false, c1=false, c2=false, c3=false) = Poly3(c0, c1, c2, c3)

(p::Poly3)(x) = p.c0 + p.c1 * x + p.c2 * x^2 + p.c3 * x^3

struct AddCall{F, G} <: BroadcastableCallable
    f::F
    g::G
end

(c::AddCall)(x) = c.f(x) + c.g(x)

macro test_inferred(ex)
    esc(:($Test.@test (($Test.@inferred $ex); true)))
end

macro test_broken_inferred(ex)
    esc(:($Test.@test_broken (($Test.@inferred $ex); true)))
end


text:
# Structured results of esimation

# y:            Time series
# θ:            Estimates
# pars:         Estimates in parameter struct
# λ:            Fitted values
# residuals:    Residuals (for INARMA using most probable sequence of R)
# LL:           Maximum value of log-likelihood
# LLs:          Single log-likelihood contributions
# nPar:         Number of parameters
# nObs:         Number of Observations
# se:           Standard errors
# CI:           Confidence intervals (normal approximation)
# model:        Model specification
# converged:    Success of optimization routine

"""
    INGARCHResults(y, θ, pars, λ, residuals,
          LL, LLs, nPar, nObs, se, CI,
          model, converged, MLEControl)
Structure for estimation results.

* `y`: Time series
* `θ`: Estimates (vector)
* `pars`: Estimates (parameter)
* `λ`: Conditional means
* `residuals`: Residuals (y - λ)
* `LL`: Maximum of likelihood
* `LLs`: Likelihood contributions
* `nPar`: Number of parameters
* `nObs`: Number of observations
* `se`: Standard errors
* `CI`: Confidence intervals
* `model`: Model specification
* `converged`: Indicator, convergence of optimization routine?
* `MLEControl`: Estimation settings used
"""
mutable struct INGARCHresults
    y::Vector{T} where T<:Integer
    θ::Vector{T} where T<: AbstractFloat
    pars::parameter
    λ::Vector{T} where T<:AbstractFloat
    residuals::Vector{T} where T<:AbstractFloat
    LL::T where T<: AbstractFloat
    LLs::Array{T, 1} where T<:AbstractFloat
    nPar::T where T<:Integer
    nObs::T where T<: Integer
    se::Array{T, 1} where T<:AbstractFloat
    CI::Array{T, 2} where T<:AbstractFloat
    model::T where T<:INGARCH
    converged::Bool
    MLEControl::MLEControl
end

"""
    INARMAResults(y, θ, pars,
      LL, LLs, nPar, nObs, se, CI,
      model, converged, MLEControl)

* `y`: Time series
* `θ`: Estimates (vector)
* `pars`: Estimates (parameter)
* `LL`: Maximum of likelihood
* `LLs`: Likelihood contributions
* `nPar`: Number of parameters
* `nObs`: Number of observations
* `se`: Standard errors
* `CI`: Confidence intervals
* `model`: Model specification
* `converged`: Indicator, convergence of optimization routine?
* `MLEControl`: Estimation settings used
"""
mutable struct INARMAresults
    y::Vector{T} where T<:Integer
    θ::Vector{T} where T<: AbstractFloat
    pars::parameter
    LL::T where T<: AbstractFloat
    LLs::Array{T, 1} where T<:AbstractFloat
    nPar::T where T<:Integer
    nObs::T where T<: Integer
    se::Array{T, 1} where T<:AbstractFloat
    CI::Array{T, 2} where T<:AbstractFloat
    model::T where T<:INARMA
    converged::Bool
    MLEControl::MLEControl
end


text:
struct QuitEvent <: T.Event end

Base.show(io::IO, ::QuitEvent) = Base.print(io, "QuitEvent")

struct App
    pipeline::Vector{Channel}
end

get_event_queue(app::App) = app.pipeline[end]

function init_pipeline(size=Inf)
    sequence_queue = Channel{String}(size, spawn=true) do ch
        while true
            sequence = T.read_stream()
            put!(ch, sequence)
        end
    end
    event_queue = Channel{T.Event}(size, spawn=true) do ch
        while true
            sequence = take!(sequence_queue)
            put!(ch, T.parse_sequence(sequence))
        end
    end

    return [sequence_queue, event_queue]
end

function initial_app()
    pipeline = init_pipeline()
    app = App(pipeline)

    return app
end

emit_quit_event(app::App) = put!(get_event_queue(app), QuitEvent())

function handle_quit(app::App)
    keep_running = false
    foreach(close, app.pipeline)
    T.println("Shutted down...")

    return keep_running
end

function handle_event(app::App)
    T.raw!(true)
    is_running = true
    while is_running
        e = take!(get_event_queue(app))
        @show e

        # sleep(1) # previous time-consuming calculation
        if e === QuitEvent()
            is_running = handle_quit(app)
        elseif e === T.KeyPressedEvent(T.ESC)
            emit_quit_event(app)
        end
        # sleep(1) # next time-consuming calculation
    end
    T.raw!(false)
end

function debug_manually()
    app = initial_app()
    handle_event(app)
end

debug_manually()


text:
#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic
using IntervalContractors

#Preamble
setprecision(53)
setprecision(Interval, Float64)

# Using approximate checks for validaty update later?
import Base.isapprox
isapprox(x::Interval,y::Interval) = isapprox(x.lo,y.lo,atol=1E-4) && isapprox(x.hi,y.hi,atol=1E-4)

@testset "exp2_rev_test" begin
    @test exp2_rev(∅, entireinterval(Float64))[2] == ∅
    @test exp2_rev(Interval(0.0, 1.0), ∅)[2] == ∅
    @test exp2_rev(∅, Interval(0.0, 1.0))[2] == ∅
    @test exp2_rev(Interval(-2.0, -1.0), entireinterval(Float64))[2] == ∅
    @test exp2_rev(Interval(1.0, 1.0), entireinterval(Float64))[2] == Interval(0.0, 0.0)
    @test exp2_rev(entireinterval(Float64), entireinterval(Float64))[2] == Interval(-∞, ∞)
    @test exp2_rev(Interval(-Inf, 0.0), entireinterval(Float64))[2] == ∅
end

@testset "exp10_rev_test" begin
    @test exp10_rev(∅, entireinterval(Float64))[2] == ∅
    @test exp10_rev(Interval(0.0, 1.0), ∅)[2] == ∅
    @test exp10_rev(∅, Interval(0.0, 1.0))[2] == ∅
    @test exp10_rev(Interval(-2.0, -1.0), entireinterval(Float64))[2] == ∅
    @test exp10_rev(Interval(1.0, 1.0), entireinterval(Float64))[2] == Interval(0.0, 0.0)
    @test exp10_rev(entireinterval(Float64), entireinterval(Float64))[2] == Interval(-∞, ∞)
    @test exp10_rev(Interval(-Inf, 0.0), entireinterval(Float64))[2] == ∅
end

@testset "expm1_rev_test" begin
    @test expm1_rev(∅, entireinterval(Float64))[2] == ∅
    @test expm1_rev(Interval(0.0, 1.0), ∅)[2] == ∅
    @test expm1_rev(∅, Interval(0.0, 1.0))[2] == ∅
    @test expm1_rev(Interval(-2.0, -1.0), entireinterval(Float64))[2] == ∅
    @test isapprox(expm1_rev(Interval(1.0, 1.0), entireinterval(Float64))[2],Interval(0.693147, 0.693148))
    @test expm1_rev(entireinterval(Float64), entireinterval(Float64))[2] == Interval(-∞, ∞)
    @test expm1_rev(Interval(-Inf, -1.0), entireinterval(Float64))[2] == ∅
end

@testset "log2_rev_test" begin
    @test log2_rev(∅, entireinterval(Float64))[2] == ∅
    @test log2_rev(Interval(0.0, 1.0), ∅)[2] == ∅
    @test log2_rev(∅, Interval(0.0, 1.0))[2] == ∅
    @test log2_rev(Interval(-2.0, -1.0), entireinterval(Float64))[2] == Interval(0.25, 0.5)
    @test log2_rev(Interval(1.0, 1.0), entireinterval(Float64))[2] == Interval(2.0, 2.0)
    @test log2_rev(entireinterval(Float64), entireinterval(Float64))[2] == Interval(0.0, Inf)
    @test log2_rev(Interval(-Inf, 0.0), entireinterval(Float64))[2] == Interval(0.0, 1.0)
end

@testset "log10_rev_test" begin
    @test log10_rev(∅, entireinterval(Float64))[2] == ∅
    @test log10_rev(Interval(0.0, 1.0), ∅)[2] == ∅
    @test log10_rev(∅, Interval(0.0, 1.0))[2] == ∅
    @test isapprox(log10_rev(Interval(-2.0, -1.0), entireinterval(Float64))[2], Interval(0.00999999, 0.100001))
    @test log10_rev(Interval(1.0, 1.0), entireinterval(Float64))[2] == Interval(10.0, 10.0)
    @test log10_rev(entireinterval(Float64), entireinterval(Float64))[2] == Interval(0.0, Inf)
    @test log10_rev(Interval(-Inf, 0.0), entireinterval(Float64))[2] == Interval(0.0, 1.0)
end

@testset "log1p_rev_test" begin
    @test log1p_rev(∅, entireinterval(Float64))[2] == ∅
    @test log1p_rev(Interval(0.0, 1.0), ∅)[2] == ∅
    @test log1p_rev(∅, Interval(0.0, 1.0))[2] == ∅
    @test isapprox(log1p_rev(Interval(-2.0, -1.0), entireinterval(Float64))[2], Interval(-0.864665, -0.63212))
    @test isapprox(log1p_rev(Interval(1.0, 1.0), entireinterval(Float64))[2], Interval(1.71828, 1.71829))
    @test log1p_rev(entireinterval(Float64), entireinterval(Float64))[2] == Interval(-1.0, Inf)
    @test log1p_rev(Interval(-Inf, 0.0), entireinterval(Float64))[2] == Interval(-1.0, 0.0)
end


text:
struct KnownBounds
    min::Float32
    max::Float32
end

# Holds the min-max values of the tree.
mutable struct MinMaxStats
    minimum::Float32
    maximum::Float32
end

MinMaxStats(known_bounds::KnownBounds) = MinMaxStats(known_bounds.min, known_bounds.max)
MinMaxStats() = MinMaxStats(typemin(Float32), typemax(Float32))
function MinMaxStats(known_bounds::KnownBounds, k)
    [MinMaxStats(known_bounds) for _ in 1:k]
end

function update!(x::MinMaxStats, value::Float32)
    x.minimum = min(x.minimum, value)
    x.maximum = max(x.maximum, value)
end

function normalize(x::MinMaxStats, value::Float32)
    if x.maximum > x.minimum
        # We normalize only when we have set the maximum and minimum values.
        return (value - x.minimum) / (x.maximum - x.minimum)
    end
    value
end


text:
# Cartestian geometry functions

module CartesianGeometry

using GeometricalPredicates
using MatrixHelper

import
	Base.isless,
	Base.==,
	#GeometricalPredicates.area,
	GeometricalPredicates.getx,
	GeometricalPredicates.gety

export  # Export functions
	#area,
	find_bounding_polygon,
	generate_triangles,
	getx,
	gety,
	intersection,
	#isbetween,
	isless,
	line_isequal,
	line2D_xsolve,
	line2D_ysolve,
	rand_line2D,
	segment_intersection,
	simple_convex_polygon_area,
	triangle_area,

	# types
	Point2D_Unbounded,

	# constants
	unbounded_max_coord,
	unbounded_min_coord

# unscaled 2D point for use with GeometricalPredicates
const unbounded_max_coord = Float64(2^24)
const unbounded_min_coord = -unbounded_max_coord

immutable Point2D_Unbounded <: AbstractPoint2D
	_x::Float64
	_y::Float64
	_min_x::Float64
	_max_x::Float64
	_min_y::Float64
	_max_y::Float64
	Point2D_Unbounded(x::Float64,y::Float64,min_x::Float64,max_x::Float64,min_y::Float64,max_y::Float64) = new(x,y,min_x,max_x,min_y,max_y)
end
Point2D_Unbounded(x::Float64, y::Float64) = Point2D_Unbounded(x,y,unbounded_min_coord,unbounded_max_coord,unbounded_min_coord,unbounded_max_coord)

function rescale_coordinate(c::Float64, cmin::Float64, cmax::Float64)
	if c < cmin # Error?
		c = cmin
	end
	if c > cmax # Error?
		c = cmax
	end
	if c == cmax
		return max_coord
	end
	return min_coord + (c-cmin)/(cmax - cmin)
end

GeometricalPredicates.getx(p::Point2D_Unbounded) = rescale_coordinate(p._x, p._min_x, p._max_x)
GeometricalPredicates.gety(p::Point2D_Unbounded) = rescale_coordinate(p._y, p._min_y, p._max_y)

function isless{T<:AbstractPoint2D}(p1::T, p2::T)
	(getx(p1) < getx(p2)) || (getx(p1) == getx(p2) && gety(p1) < gety(p2))
end

# Line functions
function line_isequal(ax, ay, bx, by, cx, cy, dx, dy)
	return ((cy-ay)*(bx-ax) == (by-ay)*(cx-ax)) && ((dy-ay)*(bx-ax) == (by-ay)*(dx-ax))
end

function =={T<:AbstractPoint2D}(lineA::Line2D{T}, lineB::Line2D{T})
	return (orientation(lineA, geta(lineB)) == 0 && orientation(lineA, getb(lineB)) == 0)
end

function line2D_ysolve{T<:AbstractPoint2D}(line::Line2D{T}, x)
	return (gety(getb(line))-gety(geta(line)))*(x-getx(geta(line)))/(getx(getb(line))-getx(geta(line)))+gety(geta(line))
end

function line2D_xsolve{T<:AbstractPoint2D}(line::Line2D{T}, y)
	return (y-gety(geta(line)))*(getx(getb(line))-getx(geta(line)))/(gety(getb(line))-gety(geta(line)))+getx(geta(line))
end

function rand_line2D()
	coords = sortrows(rand(2,2))
	return Line2D(Point(coords[1],coords[3]), Point(coords[2],coords[4]))
end

function intersection{T<:AbstractPoint2D}(line::Line2D{T}, point::T)
	if orientation(line, point) == 0
		return point
	else
		return AbstractPoint2D
	end
end
intersection{T<:AbstractPoint2D}(point::T, line::Line2D{T}) = intersection(line, point)

function intersection{T<:AbstractPoint2D}(lineA::Line2D{T}, lineB::Line2D{T})
	# lines are the same
	if lineA == lineB
		return lineA
	end

	lineA_dx = getx(getb(lineA)) - getx(geta(lineA))
	lineB_dx = getx(getb(lineB)) - getx(geta(lineB))

	# Vertical lines
	if (lineA_dx == 0. && lineB_dx == 0.) # Both lines vertical & not identical (test above)
		return AbstractPoint2D
	elseif lineA_dx == 0. # lineA is vertical
		return Point(getx(geta(lineA)), line2D_ysolve(lineB, getx(geta(lineA))))
	elseif lineB_dx == 0. # lineB is vertical
		return Point(getx(geta(lineB)), line2D_ysolve(lineA, getx(geta(lineB))))
	end

	lineA_slope = (gety(getb(lineA)) - gety(geta(lineA))) / lineA_dx
	lineA_int = gety(geta(lineA)) - lineA_slope*getx(geta(lineA))

	lineB_slope = (gety(getb(lineB)) - gety(geta(lineB))) / lineB_dx
	lineB_int = gety(geta(lineB)) - lineB_slope*getx(geta(lineB))

	# Lines are parallel & not identical (test above)
	if lineA_slope == lineB_slope
		return AbstractPoint2D
	end

	return Point((lineB_int-lineA_int)/(lineA_slope-lineB_slope),
			lineA_slope*(lineB_int-lineA_int)/(lineA_slope-lineB_slope)+lineA_int)
end

function isbetween(x, a, b)
	if a >= b
		return b <= x <= a
	else
		return a <= x <= b
	end
end

function segment_intersection{T<:AbstractPoint2D}(line::Line2D{T}, point::T)
	pt = intersection(line, point)

	if (pt == point &&
		isbetween(getx(pt), getx(geta(line)), getx(getb(line))) &&
		isbetween(gety(pt), gety(geta(line)), gety(getb(line))) )
		return pt
	else
		return AbstractPoint2D
	end
end
segment_intersection{T<:AbstractPoint2D}(point::T, line::Line2D{T}) = segment_intersection(line,point)

function segment_intersection{T<:AbstractPoint2D}(lineA::Line2D{T}, lineB::Line2D{T})
	pt = intersection(lineA, lineB)

	if pt == segment_intersection(lineA, pt) && pt == segment_intersection(lineB, pt)
		return pt
	else
		return AbstractPoint2D
	end
end

# Polygon functions

function intersection{T<:AbstractPoint2D}(polygon::Polygon2D{T}, line::Line2D{T})
	intersection_points = AbstractPoint2D[]

	for l in getlines(polygon)
		if orientation(line, geta(l)) != orientation(line, getb(l))
			line_int = intersection(line, l)
			if line_int != AbstractPoint2D
				push!(intersection_points, line_int)
			end
		end
	end

	if length(intersection_points) > 0
		return intersection_points
	else
		return AbstractPoint2D
	end
end
intersection{T<:AbstractPoint2D}(line::Line2D{T}, polygon::Polygon2D{T}) = intersection(polygon, line)

function generate_triangles{T}(pointarray::Vector{T})
	startpt = selectperm(pointarray,1)

	vertextups = [ map(x->(x<size(pointarray,1)+1?x:mod(x,size(pointarray,1))), (startpt, startpt+i, startpt + i + 1) ) for i in 1:(size(pointarray,1)-2) ]

	return [ pointarray[union(x),:] for x in vertextups ]
end

function generate_triangles{T<:AbstractMatrix}(pointarray::T)
	startpt = selectpermrows(pointarray,1)

	vertextups = [ map(x->(x<size(pointarray,1)+1?x:mod(x,size(pointarray,1))), (startpt, startpt+i, startpt + i + 1) ) for i in 1:(size(pointarray,1)-2) ]

	return [ pointarray[union(x),:] for x in vertextups ]
end

function triangle_area{T<:Number}(verticies::Array{T,2})
	v = broadcast(-, verticies, verticies[1,:])

	return abs(v[2,1]*v[3,2] - v[2,2]*v[3,1])/2.
end

function triangle_area{T<:AbstractPoint2D}(m::Matrix{T})
	return triangle_area(vec(m))
end

function triangle_area{T<:AbstractPoint2D}(v::Vector{T})
	if length(v) != 3
		error("Not a triangle")
	end

	orig_x = getx(v[1])
	orig_y = gety(v[1])

	return abs((getx(v[2])-orig_x)*(gety(v[3])-orig_y) - (gety(v[2])-orig_y)*(getx(v[3])-orig_x))*0.5
end

function simple_convex_polygon_area(ptarray)
	trianglelist = generate_triangles(ptarray)

	return mapreduce(triangle_area, +, trianglelist)
end

end # Module


text:
# ---
# title: Rotate
# cover: rotate.gif
# description: rotate image anticlockwise
# ---

# The type [`Rotate`](@ref) defines a generic anticlockwise rotation operation around the center
# of the image. It is also possible to pass some abstract vector to the constructor, in which case
# Augmentor will randomly sample one of its elements every time the operation is applied.

using Augmentor
using ImageShow, ImageCore
using Random
Random.seed!(0)

img_in = testpattern(RGB, ratio=0.5)
mosaicview(
    img_in,

    ## deterministic rotation
    augment(img_in, Rotate(45)),

    ## random rotation
    augment(img_in, Rotate(-45:45));
    fillvalue=colorant"white", nrow=1, npad=10
)

# Note that the output image size will be changed after rotation, [`CropNative`](@ref) can be particalually
# useful to preserve the image size.

mosaicview(
    augment(img_in, Rotate(45)),
    augment(img_in, Rotate(45) |> CropNative(axes(img_in)));
    nrow=1, npad=10
)

# Rotation by some special degree (e.g.,90, 180 and 270) can be handled more efficiently without interpolation.
# Compared to `Rotate(90)`, it is recommended to use [`Rotate90`](@ref) when possible. [`Rotate180`](@ref) and
# [`Rotate270`](@ref) are available, too.

# ## References

#md # ```@docs
#md # Rotate
#md # Rotate90
#md # Rotate180
#md # Rotate270
#md # ```

## save covers #src
using ImageMagick #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(testpattern(RGB, ratio=0.5), Rotate(-20:20), 5) #src
ImageMagick.save("rotate.gif", cover; fps=1) #src


text:
@enum Boundary_type Dirichlet Neumann

struct Boundary
    mom_bot_type::Boundary_type
    mom_top_type::Boundary_type
    s_bot_type::Boundary_type
    s_top_type::Boundary_type
end

function string_to_type(s::String)
    if s == "Dirichlet"
        return Dirichlet::Boundary_type
    elseif s == "Neumann"
        return Neumann::Boundary_type
    end
end

function Boundary(d::Dict)
    mom_bot_type = string_to_type(d["mom_bot_type"])
    mom_top_type = string_to_type(d["mom_top_type"])
    s_bot_type = string_to_type(d["s_bot_type"])
    s_top_type = string_to_type(d["s_top_type"])

    Boundary(
        mom_bot_type,
        mom_top_type,
        s_bot_type,
        s_top_type)
end

function boundary_cyclic_kernel!(
    a, is, ie, js, je, igc, jgc)
    # East-west BCs
    @tturbo for k in 1:size(a, 3)
        for j in 1:size(a, 2)
            for i in 1:igc
                a[i, j, k] = a[ie-i+1, j, k]
                a[ie+i, j, k] = a[is+i-1, j, k]
            end
        end
    end

    # North-south BCs
    @tturbo for k in 1:size(a, 3)
        for j in 1:jgc
            for i in 1:size(a, 1)
                a[i, j, k] = a[i, je-j+1, k]
                a[i, je+j, k] = a[i, js+j-1, k]
            end
        end
    end
end

function boundary_cyclic_kernel_2d!(
    a, is, ie, js, je, igc, jgc)
    # East-west BCs
    @tturbo for j in 1:size(a, 2)
        for i in 1:igc
            a[i, j] = a[ie-i+1, j]
            a[ie+i, j] = a[is+i-1, j]
        end
    end

    # North-south BCs
    @tturbo for j in 1:jgc
        for i in 1:size(a, 1)
            a[i, j] = a[i, je-j+1]
            a[i, je+j] = a[i, js+j-1]
        end
    end
end

function set_ghost_cells_bot_kernel!(
    a, a_bot, a_gradbot, dzh, ks, bot_type::Boundary_type)

    if bot_type == Dirichlet::Boundary_type
        dzhi = 1 / dzh[ks]
        @tturbo for j in 1:size(a, 2)
            for i in 1:size(a, 1)
                a[i, j, ks-1] = 2a_bot[i, j] - a[i, j, ks]
                a_gradbot[i, j] = (a[i, j, ks] - a[i, j, ks-1]) * dzhi
            end
        end

    elseif bot_type == Neumann::Boundary_type
        @tturbo for j in 1:size(a, 2)
            for i in 1:size(a, 1)
                a[i, j, ks-1] = -a_gradbot[i, j]*dzh[ks] + a[i, j, ks]
                a_bot[i, j] = 1//2 * (a[i, j, ks-1] + a[i, j, ks])
            end
        end
    end
end

function set_ghost_cells_top_kernel!(
    a, a_top, a_gradtop, dzh, ke, bot_type::Boundary_type)

    if bot_type == Dirichlet::Boundary_type
        dzhi = 1 / dzh[ke+1]
        @tturbo for j in 1:size(a, 2)
            for i in 1:size(a, 1)
                a[i, j, ke+1] = 2a_top[i, j] - a[i, j, ke]
                a_gradtop[i, j] = (a[i, j, ke+1] - a[i, j, ke]) * dzhi
            end
        end

    elseif bot_type == Neumann::Boundary_type
        @tturbo for j in 1:size(a, 2)
            for i in 1:size(a, 1)
                a[i, j, ke+1] = a_gradtop[i, j]*dzh[ke+1] + a[i, j, ke]
                a_top[i, j] = 1//2 * (a[i, j, ke] + a[i, j, ke+1])
            end
        end
    end
end

function set_boundary!(fields::Fields, grid::Grid, boundary::Boundary)
    # Cyclic BC.
    boundary_cyclic_kernel!(
        fields.u, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel!(
        fields.v, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel!(
        fields.w, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel!(
        fields.s, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)

    # Cyclic BCs of boundary fields.
    boundary_cyclic_kernel_2d!(
        fields.u_bot, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel_2d!(
        fields.u_gradbot, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel_2d!(
        fields.u_top, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel_2d!(
        fields.u_gradtop, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)

    boundary_cyclic_kernel_2d!(
        fields.v_bot, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel_2d!(
        fields.v_gradbot, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel_2d!(
        fields.v_top, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel_2d!(
        fields.v_gradtop, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)

    boundary_cyclic_kernel_2d!(
        fields.s_bot, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel_2d!(
        fields.s_gradbot, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel_2d!(
        fields.s_top, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)
    boundary_cyclic_kernel_2d!(
        fields.s_gradtop, grid.is, grid.ie, grid.js, grid.je, grid.igc, grid.jgc)

    # Bottom BC.
    set_ghost_cells_bot_kernel!(
        fields.u, fields.u_bot, fields.u_gradbot, grid.dzh, grid.ks, boundary.mom_bot_type)
    set_ghost_cells_bot_kernel!(
        fields.v, fields.v_bot, fields.v_gradbot, grid.dzh, grid.ks, boundary.mom_bot_type)
    set_ghost_cells_bot_kernel!(
        fields.s, fields.s_bot, fields.s_gradbot, grid.dzh, grid.ks, boundary.s_bot_type)

    # Top BC.
    set_ghost_cells_top_kernel!(
        fields.u, fields.u_top, fields.u_gradtop, grid.dzh, grid.ke, boundary.mom_top_type)
    set_ghost_cells_top_kernel!(
        fields.v, fields.v_top, fields.v_gradtop, grid.dzh, grid.ke, boundary.mom_top_type)
    set_ghost_cells_top_kernel!(
        fields.s, fields.s_top, fields.s_gradtop, grid.dzh, grid.ke, boundary.s_top_type)
end


text:
using Optimisers
using ChainRulesCore, Functors, StaticArrays, Zygote
using LinearAlgebra, Statistics, Test, Random

Random.seed!(1)

RULES = [
  # All the rules at default settings:
  Descent(), ADAM(), Momentum(), Nesterov(), RMSProp(),
  ADAGrad(), AdaMax(), ADADelta(), AMSGrad(), NADAM(),
  ADAMW(), RADAM(), OADAM(), AdaBelief(),
  # A few chained combinations:
  OptimiserChain(WeightDecay(), ADAM(0.001)),
  OptimiserChain(ClipNorm(), ADAM(0.001)),
  OptimiserChain(ClipGrad(0.5), Momentum()),
  OptimiserChain(WeightDecay(), OADAM(), ClipGrad(1)),
]

name(o) = typeof(o).name.name  # just for printing testset headings
name(o::OptimiserChain) = join(name.(o.opts), " → ")

LOG = Dict()  # for debugging these testsets, this makes it easy to plot each optimiser's loss

loggradient(o) = (f, xs...) -> begin
  y, dxs = Zygote.withgradient(f, xs...)
  push!(get!(() -> Float32[], LOG, name(o)), y)
  dxs  # save the loss, return the gradient
end

@testset "independence" begin
  empty!(LOG)
  @testset "$(name(o))" for o in RULES
    w = randn(10, 10)
    w′ = randn(10, 10)
    iloss(x, w, w′) = mean((w*x .- w′*x) .^ 2)
    @test iloss(rand(10, 10), w, w′) > 1
    st = Optimisers.setup(o, w)
    for t = 1:10^5
      x = rand(10)
      gs = loggradient(o)(w -> iloss(x, w, w′), w)
      st, w = Optimisers.update!(st, w, gs...)
    end
    @test iloss(rand(10, 10), w, w′) < 0.01
  end
end

@testset "simple sum" begin
  empty!(LOG)
  @testset "$(name(o))" for o in RULES
    m = shuffle!(reshape(1:64, 8, 8) .+ 0.0)
    s = Optimisers.setup(o, m)
    for _ in 1:10^5
      g = loggradient(o)(x -> sum(abs2, x + x'), m)[1]
      s, m = Optimisers.update!(s, m, g)
    end
    @test sum(m) < sum(1:64)
    if sum(m) < 1
      @test sum(m) < 1
    else
      @show name(o) sum(m)/sum(1:64)
      @test_broken sum(m) < 1
    end
  end
end

@testset "original" begin
  empty!(LOG)
  @testset "$(name(o))" for o in RULES
    w′ = (α = rand(3, 3), β = rand(3, 3))
    w = (α = 5rand(3, 3), β = rand(3, 3))
    st = Optimisers.setup(o, w)
    loss(x, y) = mean((x.α .* x.β .- y.α .* y.β) .^ 2)
    @test loss(w, w′) > 1
    for i = 1:10^4
      gs = loggradient(o)(x -> loss(x, w′), w)
      st, w = Optimisers.update(st, w, gs...)
    end
    @test loss(w, w′) < 0.001
  end
end

@testset "StaticArrays" begin
  empty!(LOG)
  @testset "$(name(o))" for o in RULES
    W1 = @SMatrix randn(10, 10)
    b1 = @SVector randn(10)
    W2 = @SMatrix randn(10, 10)
    model = (; W1, b1, W2, tanh)
    s_loss(m, x, y) = sum(abs2, m.W2 * (m.tanh).(m.W1*x .+ m.b1) .- y)
    # x = @SMatrix randn(10, 10)
    # y = @SMatrix randn(10, 10)  # gives an error from sum(; dims=())
    x = @SVector randn(10)
    y = @SVector randn(10)
    @test s_loss(model, x, y) > 10
    state = Optimisers.setup(o, model)
    for t = 1:10^3
      g = loggradient(o)(m -> s_loss(m, x, y), model)[1]
      state, model = Optimisers.update!(state, model, g)
    end
    if o isa Descent
      @show name(o) s_loss(model, x, y)
      @test_broken s_loss(model, x, y) < 1
    else
      @test s_loss(model, x, y) < 1
    end
  end
end

@testset "element types" begin
  @testset "$(name(o))" for o in RULES
    marray = (Float32[1,2], Float64[3,4], Float16[5,6])
    types = map(eltype, marray)

    # This is a weak test, as it copies & then does `update!`
    uparray = Optimisers.update(Optimisers.setup(o, marray), marray, marray)[2]
    @test map(eltype, uparray) == types

    # Static version is truly out-of-place:
    mstatic = (SA{Float32}[1,2], SA{Float64}[3,4]) # , SA{Float16}[5,6])  with Float16, all fail
    upstatic = Optimisers.update(Optimisers.setup(o, mstatic), mstatic, mstatic)[2]
    @test map(eltype, upstatic) == types[1:2]
    @test upstatic[1] isa SVector

    # With ordinary Array gradient, what happens? Not so important!
    upstatic2 = Optimisers.update(Optimisers.setup(o, mstatic), mstatic, marray[1:2])[2]
    # @test map(eltype, upstatic2) == types[1:2]  # same information
    if upstatic2[1] isa SVector
      @test upstatic2[1] isa SVector
    else
      @test_broken upstatic2[1] isa SVector
    end
  end
end

@testset "gradient types" begin
  @testset "$(name(o))" for o in RULES
    x = (a = ones(2,2), b = transpose(ones(2,2)))
    s = Optimisers.setup(o, x)

    _, x1 = Optimisers.update(s, x, (a = [1 2; 3 4], b = nothing))
    @test x1.a != ones(2,2)
    @test x1.b == ones(2,2)

    _, xfill = Optimisers.update(s, x, (a = Zygote.Fill(2.0,2,2), b = Zygote.Fill(true,2,2)))
    @test xfill.a != ones(2,2)
    @test xfill.b != ones(2,2)

    bc = Optimisers.@lazy 1 + log([2 3; 4 5]) / 6
    _, xbc = Optimisers.update(s, x, (a = bc, b = bc))
    @test xbc.a != ones(2,2)
    @test xbc.b != ones(2,2)

    th = ChainRulesCore.@thunk @. 1 + log([2 3; 4 5]) / 6
    _, xth = Optimisers.update(s, x, (a = bc, b = bc))
    @test xth.a != ones(2,2)
    @test xth.b != ones(2,2)
  end
end

@testset "mutation check" begin
  # If @lazy captures a matrix which is later mutated, the results won't agree here:
  @testset "$(name(o))" for o in RULES
    model = Float64.(rand(Int8, 8))
    s_model = SVector{8}(model)
    grads = [Float64.(rand(Int8, 8)) for t in 1:13]
    s_grads = [SVector{8}(x) for x in grads]
    state = Optimisers.setup(o, model)
    s_state = Optimisers.setup(o, s_model)
    for t in 1:13
      state, model = Optimisers.update!(state, model, grads[t])
      s_state, s_model = Optimisers.update!(s_state, s_model, s_grads[t])
    end
    @test model == s_model
  end
end

@testset "with complex numbers: Flux#1776" begin
  empty!(LOG)
  @testset "$(name(opt))" for opt in [
              # The Flux PR had 1e-2 for all. But ADADelta(ρ) needs ρ≈0.9 not small. And it helps to make ε not too small too:
              ADAM(1e-2), RMSProp(1e-2), RADAM(1e-2), OADAM(1e-2), ADAGrad(1e-2), ADADelta(0.9, 1e-5), NADAM(1e-2), AdaBelief(1e-2),
              # These weren't in Flux PR:
              Descent(1e-2), Momentum(1e-2), Nesterov(1e-2), ADAMW(1e-2),
              ]
    # Our "model" is just a complex number
    model = (w = zeros(ComplexF64, 1),)

    # Our model attempts to learn `f(x) = conj(x)` where `f(x) = w*x`
    function loss(m)
      # Deterministic training data is the best training data
      x = ones(1, 1) + 1im*ones(1, 1)
      # Manually implement `mse()` to allow demonstration of brokenness
      # on older Flux builds that don't have a fixed `mse()`
      return sum(abs2.(m.w * x .- conj(x)))
    end
    @test loss(model) ≈ 2.0

    state = Optimisers.setup(opt, model)

    # Train for 10 iterations, enforcing that loss is monotonically decreasing
    last_loss = Inf
    for idx in 1:10
      grads = loggradient(opt)(loss, model)
      state, model = Optimisers.update!(state, model, grads...)
      opt isa Union{Momentum, Nesterov} && idx > 8 && continue  # these are very flat at the end
      @test loss(model) < last_loss
      last_loss = loss(model)
    end
    @test loss(model) < 1.9

    # Repeat with StaticArrays
    static_model = (w = SA[0.0 + 0im],)
    static_state = Optimisers.setup(opt, static_model)
    function static_loss(m)
      x = hcat(SA[1.0 + im])
      sum(abs2.(m.w * x .- conj(x)))
    end
    last_loss = Inf
    for idx in 1:10
      grads = gradient(static_loss, static_model)
      static_state, static_model = Optimisers.update!(static_state, static_model, grads...)
      opt isa Union{Momentum, Nesterov} && idx > 8 && continue
      @test static_loss(static_model) < last_loss
      last_loss = static_loss(static_model)
    end
    @test static_loss(static_model) < 1.9
  end
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

tests for MOI wrapper
=#

using Test
using LinearAlgebra
import SparseArrays
import Random
import MathOptInterface
const MOI = MathOptInterface
import Hypatia
import Hypatia.Cones

function test_moi_cones(T::Type{<:Real})
    # MOI predefined cones

    @testset "PositiveSemidefiniteConeTriangle" begin
        moi_cone = MOI.PositiveSemidefiniteConeTriangle(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.PosSemidefTri{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 6
        @test !Cones.use_dual_barrier(hyp_cone)
    end

    @testset "NormInfinityCone" begin
        moi_cone = MOI.NormInfinityCone(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiNormInf{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test !Cones.use_dual_barrier(hyp_cone)
    end

    @testset "NormOneCone" begin
        moi_cone = MOI.NormOneCone(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiNormInf{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test Cones.use_dual_barrier(hyp_cone)
    end

    @testset "SecondOrderCone" begin
        moi_cone = MOI.SecondOrderCone(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiNormEucl{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test !Cones.use_dual_barrier(hyp_cone)
    end

    @testset "RotatedSecondOrderCone" begin
        moi_cone = MOI.RotatedSecondOrderCone(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiPerSquare{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test !Cones.use_dual_barrier(hyp_cone)
    end

    @testset "NormSpectralCone" begin
        moi_cone = MOI.NormSpectralCone(2, 3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiNormSpectral{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7
        @test !Cones.use_dual_barrier(hyp_cone)
    end

    @testset "NormNuclearCone" begin
        moi_cone = MOI.NormNuclearCone(2, 3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiNormSpectral{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7
        @test Cones.use_dual_barrier(hyp_cone)
    end

    @testset "PowerCone" begin
        iT5 = inv(T(5))
        moi_cone = MOI.PowerCone{T}(iT5)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.GeneralizedPower{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test !Cones.use_dual_barrier(hyp_cone)
        @test hyp_cone.α == T[iT5, 1 - iT5]
    end

    @testset "DualPowerCone" begin
        iT5 = inv(T(5))
        moi_cone = MOI.DualPowerCone{T}(iT5)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.GeneralizedPower{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test Cones.use_dual_barrier(hyp_cone)
        @test hyp_cone.α == T[iT5, 1 - iT5]
    end

    @testset "GeometricMeanCone" begin
        moi_cone = MOI.GeometricMeanCone(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoGeoMean{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test !Cones.use_dual_barrier(hyp_cone)
    end

    @testset "RootDetConeTriangle" begin
        moi_cone = MOI.RootDetConeTriangle(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoRootdetTri{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7
        @test !Cones.use_dual_barrier(hyp_cone)
    end

    @testset "ExponentialCone" begin
        moi_cone = MOI.ExponentialCone()
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoPerLog{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test !Cones.use_dual_barrier(hyp_cone)
    end

    @testset "DualExponentialCone" begin
        moi_cone = MOI.DualExponentialCone()
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoPerLog{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test Cones.use_dual_barrier(hyp_cone)
    end

    @testset "LogDetConeTriangle" begin
        moi_cone = MOI.LogDetConeTriangle(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoPerLogdetTri{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 8
        @test !Cones.use_dual_barrier(hyp_cone)
    end

    @testset "RelativeEntropyCone" begin
        moi_cone = MOI.RelativeEntropyCone(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiRelEntropy{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test !Cones.use_dual_barrier(hyp_cone)
    end

    # Hypatia predefined cones

    @testset "Nonnegative" begin
        moi_cone = Hypatia.NonnegativeCone{T}(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.Nonnegative{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
    end

    @testset "PosSemidefTri" begin
        moi_cone = Hypatia.PosSemidefTriCone{T, T}(6)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.PosSemidefTri{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 6

        moi_cone = Hypatia.PosSemidefTriCone{T, Complex{T}}(9)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.PosSemidefTri{T, Complex{T}}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 9
    end

    @testset "DoublyNonnegativeTri" begin
        moi_cone = Hypatia.DoublyNonnegativeTriCone{T}(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.DoublyNonnegativeTri{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
    end

    @testset "PosSemidefTriSparse" begin
        for (psdsparseimpl, realtypes) in Cones.PSDSparseImplList
            (T <: realtypes) || continue

            Random.seed!(1)
            side = 5
            (row_idxs, col_idxs, _) = SparseArrays.findnz(
                tril!(SparseArrays.sprand(Bool, side, side, 0.3)) + I)

            moi_cone = Hypatia.PosSemidefTriSparseCone{psdsparseimpl, T, T}(
                side, row_idxs, col_idxs)
            hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
            @test hyp_cone isa Cones.PosSemidefTriSparse{psdsparseimpl, T, T}
            @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) ==
                length(row_idxs)

            moi_cone = Hypatia.PosSemidefTriSparseCone{psdsparseimpl, T, Complex{T}}(
                side, row_idxs, col_idxs)
            hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
            @test hyp_cone isa
                Cones.PosSemidefTriSparse{psdsparseimpl, T, Complex{T}}
            @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) ==
                2 * length(row_idxs) - side
        end
    end

    @testset "LinMatrixIneq" begin
        As = [Symmetric(Matrix(one(T) * I, 2, 2)),
            Hermitian(Complex{T}[1 0; 0 -1])]
        moi_cone = Hypatia.LinMatrixIneqCone{T}(As)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.LinMatrixIneq{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 2
        @test hyp_cone.As == As
    end

    @testset "EpiNormInf" begin
        moi_cone = Hypatia.EpiNormInfCone{T, T}(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiNormInf{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3

        moi_cone = Hypatia.EpiNormInfCone{T, Complex{T}}(5)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiNormInf{T, Complex{T}}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 5
    end

    @testset "EpiNormEucl" begin
        moi_cone = Hypatia.EpiNormEuclCone{T}(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiNormEucl{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
    end

    @testset "EpiPerSquare" begin
        moi_cone = Hypatia.EpiPerSquareCone{T}(4)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiPerSquare{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 4
    end

    @testset "EpiNormSpectral" begin
        moi_cone = Hypatia.EpiNormSpectralCone{T, T}(2, 3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiNormSpectral{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7

        moi_cone = Hypatia.EpiNormSpectralCone{T, Complex{T}}(2, 3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiNormSpectral{T, Complex{T}}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 13
    end

    @testset "MatrixEpiPerSquare" begin
        moi_cone = Hypatia.MatrixEpiPerSquareCone{T, T}(2, 3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.MatrixEpiPerSquare{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 10

        moi_cone = Hypatia.MatrixEpiPerSquareCone{T, Complex{T}}(2, 3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.MatrixEpiPerSquare{T, Complex{T}}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 17
    end

    @testset "GeneralizedPower" begin
        α = rand(T, 2)
        α ./= sum(α)
        moi_cone = Hypatia.GeneralizedPowerCone{T}(α, 3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.GeneralizedPower{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 5
        @test hyp_cone.α == α
    end

    @testset "HypoPowerMean" begin
        α = rand(T, 2)
        α ./= sum(α)
        moi_cone = Hypatia.HypoPowerMeanCone{T}(α)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoPowerMean{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test hyp_cone.α == α
    end

    @testset "HypoGeoMean" begin
        moi_cone = Hypatia.HypoGeoMeanCone{T}(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoGeoMean{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
    end

    @testset "HypoRootdetTri" begin
        moi_cone = Hypatia.HypoRootdetTriCone{T, T}(7)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoRootdetTri{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 7

        moi_cone = Hypatia.HypoRootdetTriCone{T, Complex{T}}(10)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoRootdetTri{T, Complex{T}}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 10
    end

    @testset "HypoPerLog" begin
        moi_cone = Hypatia.HypoPerLogCone{T}(4)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoPerLog{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 4
    end

    @testset "HypoPerLogdetTri" begin
        moi_cone = Hypatia.HypoPerLogdetTriCone{T, T}(8)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoPerLogdetTri{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 8

        moi_cone = Hypatia.HypoPerLogdetTriCone{T, Complex{T}}(11)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.HypoPerLogdetTri{T, Complex{T}}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 11
    end

    @testset "EpiPerSepSpectral" begin
        h_fun = Cones.NegEntropySSF()
        QTs = [
            Cones.VectorCSqr{T},
            Cones.MatrixCSqr{T, T},
            Cones.MatrixCSqr{T, Complex{T}},
            ]
        for QT in QTs
            moi_cone = Hypatia.EpiPerSepSpectralCone{T}(h_fun, QT, 3)
            hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
            @test hyp_cone isa Cones.EpiPerSepSpectral{QT, T}
            @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone)
        end
    end

    @testset "EpiRelEntropy" begin
        moi_cone = Hypatia.EpiRelEntropyCone{T}(5)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiRelEntropy{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 5
    end

    @testset "EpiTrRelEntropyTriCone" begin
        moi_cone = Hypatia.EpiTrRelEntropyTriCone{T}(3)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.EpiTrRelEntropyTri{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
    end

    @testset "WSOSInterpNonnegative" begin
        Ps = [rand(T, 3, 2), rand(T, 3, 1)]
        moi_cone = Hypatia.WSOSInterpNonnegativeCone{T, T}(3, Ps)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.WSOSInterpNonnegative{T, T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
        @test hyp_cone.Ps == Ps

        Ps = [rand(Complex{T}, 4, 3), rand(Complex{T}, 4, 2)]
        moi_cone = Hypatia.WSOSInterpNonnegativeCone{T, Complex{T}}(4, Ps)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.WSOSInterpNonnegative{T, Complex{T}}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 4
        @test hyp_cone.Ps == Ps
    end

    @testset "WSOSInterpPosSemidefTri" begin
        Ps = [rand(T, 3, 2), rand(T, 3, 1)]
        moi_cone = Hypatia.WSOSInterpPosSemidefTriCone{T}(2, 3, Ps)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.WSOSInterpPosSemidefTri{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 9
        @test hyp_cone.Ps == Ps
    end

    @testset "WSOSInterpEpiNormOne" begin
        Ps = [rand(T, 3, 2), rand(T, 3, 1)]
        moi_cone = Hypatia.WSOSInterpEpiNormOneCone{T}(2, 3, Ps)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.WSOSInterpEpiNormOne{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 6
        @test hyp_cone.Ps == Ps
    end

    @testset "WSOSInterpEpiNormEucl" begin
        Ps = [rand(T, 3, 2), rand(T, 3, 1)]
        moi_cone = Hypatia.WSOSInterpEpiNormEuclCone{T}(2, 3, Ps)
        hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
        @test hyp_cone isa Cones.WSOSInterpEpiNormEucl{T}
        @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 6
        @test hyp_cone.Ps == Ps
    end

    return
end