text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module CCAS

using Compat

if Sys.WORD_SIZE == 32
  const LIBCCAS = is_windows() ? Pkg.dir("CCAS/libccas/lib/libccas_x32.dll") : is_linux() ? Pkg.dir("CCAS/libccas/lib/libccas_x32.so") : error("OS not supported")
elseif Sys.WORD_SIZE == 64
  const LIBCCAS = is_windows() ? Pkg.dir("CCAS/libccas/lib/libccas_x64.dll") : is_linux() ? Pkg.dir("CCAS/libccas/lib/libccas_x64.so") : error("OS not supported")
else
  error("CCAS: Architecture must be 32-bit or 64-bit")
end

export Constants, CASShared, reset, reset!, version, error_msg,
            max_intruders, update!, author

using CASInterface
import Base: reset, resize!, size, getindex

get_equipage(::Type{Val{EQUIPAGE_ATCRBS}}) = ccall((:enum_EQUIPAGE_ATCRBS, LIBCCAS),Int32,())
get_equipage(::Type{Val{EQUIPAGE_MODES}}) = ccall((:enum_EQUIPAGE_MODES, LIBCCAS),Int32,())
get_equipage(::Type{Val{EQUIPAGE_TCASTA}}) = ccall((:enum_EQUIPAGE_TCASTA, LIBCCAS),Int32,())
get_equipage(::Type{Val{EQUIPAGE_TCAS}}) = ccall((:enum_EQUIPAGE_TCAS, LIBCCAS),Int32,())

get_equipage(equip::EQUIPAGE) = get_equipage(Val{equip})

type OwnInputRef
  handle::Ptr{Void}
end

type IntruderInputRef
  handle::Ptr{Void}
end

type CollectionIntruderInputRef
  handle::Ptr{Void}
  intruders::Vector{IntruderInputRef}
end
CollectionIntruderInputRef(handle::Ptr{Void}) = CollectionIntruderInputRef(handle, IntruderInputRef[])

type InputRef
  handle::Ptr{Void}
  ownInput::OwnInputRef
  intruder_collection::CollectionIntruderInputRef

  function InputRef(nintruders::Int)
    obj = new()
    obj.handle = ccall((:newCInput, LIBCCAS), Ptr{Void}, ())
    obj.ownInput = OwnInputRef(obj)
    obj.intruder_collection = CollectionIntruderInputRef(obj, nintruders)

    finalizer(obj, obj -> ccall((:delCInput, LIBCCAS), Void, (Ptr{Void},), obj.handle))
    return obj
  end
end

function OwnInputRef(input::InputRef, ownInputVals::Union{OwnInput,Void}=nothing)
  handle = ccall((:getRefCOwnInput, LIBCCAS), Ptr{Void}, (Ptr{Void},), input.handle)
  obj = OwnInputRef(handle)

  if ownInputVals != nothing
    set!(obj, ownInputVals)
  end
  return obj
end

function IntruderInputRef(collection::CollectionIntruderInputRef, index_::Integer,
                       intrVals::Union{IntruderInput,Void}=nothing)
  index = checked_convert(UInt32, index_-1) #C uses 0 indexing
  handle = ccall((:getRefCIntruderInput, LIBCCAS), Ptr{Void},
                 (Ptr{Void}, UInt32), collection.handle, index)
  obj = IntruderInputRef(handle)

  if intrVals != nothing
    set!(obj, intrVals)
  end
  return obj
end

function CollectionIntruderInputRef(input::InputRef, nintruders::Int)
  handle = ccall((:getRefCCollectionIntruderInput, LIBCCAS), Ptr{Void}, (Ptr{Void},), input.handle)
  obj = CollectionIntruderInputRef(handle)
  resize!(obj,nintruders)
  resize!(obj.intruders, nintruders)

  for i = 1:nintruders
    obj.intruders[i] = IntruderInputRef(obj, i)
  end
  return obj
end

function set!(input::OwnInputRef, inputVals::OwnInput)
  ccall((:setCOwnInput, LIBCCAS), Void,
        (Ptr{Void}, Float64, Float64, Float64, Float64, UInt32), input.handle, inputVals.dz, inputVals.z,
        inputVals.psi, inputVals.h, inputVals.modes)
end

function set!(input::IntruderInputRef, iinputVals::IntruderInput)
  equipage = get_equipage(iinputVals.equipage)
  ccall((:setCIntruderInput, LIBCCAS), Void,
        (Ptr{Void}, Bool, UInt32, UInt32, Float64, Float64, Float64,
         UInt8, UInt8, UInt8, Int32, UInt8,
         UInt8, UInt8), input.handle, iinputVals.valid, iinputVals.id,
        iinputVals.modes, iinputVals.sr,
        iinputVals.chi, iinputVals.z,
        iinputVals.cvc, iinputVals.vrc, iinputVals.vsb, equipage, iinputVals.quant,
        iinputVals.sensitivity_index, iinputVals.protection_mode)
end

function set!(input::InputRef, inputVals::Input)
  set!(input.ownInput, inputVals.ownInput)
  for (input_intruder, inputVals_intruder) in zip(input.intruder_collection.intruders,
                                                 inputVals.intruders)
    set!(input_intruder, inputVals_intruder)
  end
end

function resize!(collection::CollectionIntruderInputRef, size_::Integer)
  size = checked_convert(UInt32, size_)
  ccall((:resizeCCollectionIntrInput, LIBCCAS), Void, (Ptr{Void},UInt32), collection.handle, size)
end

function size(collection::CollectionIntruderInputRef)
  csize = ccall((:sizeCCollectionIntrInput, LIBCCAS), UInt32, (Ptr{Void},), collection.handle)
  csize = convert(Int64, csize)
  @assert csize == length(collection.intruders) #make sure we're in sync
  return csize
end

type IntruderOutputRef
  handle::Ptr{Void}
end

type CollectionIntruderOutputRef
  handle::Ptr{Void}
  intruders::Vector{IntruderOutputRef}
end
CollectionIntruderOutputRef(handle::Ptr{Void}) = CollectionIntruderOutputRef(handle, IntruderOutputRef[])

type OutputRef
  handle::Ptr{Void}
  intruder_collection::CollectionIntruderOutputRef

  function OutputRef(nintruders::Int)
    obj = new()
    obj.handle = ccall((:newCOutput, LIBCCAS), Ptr{Void}, ())
    obj.intruder_collection = CollectionIntruderOutputRef(obj, nintruders)

    finalizer(obj,obj->ccall((:delCOutput, LIBCCAS), Void, (Ptr{Void},), obj.handle))
    return obj
  end
end

function IntruderOutputRef(collection::CollectionIntruderOutputRef, index_::Integer,
                        outputVals::Union{IntruderOutput,Void}=nothing)
  index = checked_convert(UInt32, index_ - 1) #C uses 0 indexing
  handle = ccall((:getRefCIntruderOutput, LIBCCAS), Ptr{Void},
                     (Ptr{Void}, UInt32), collection.handle, index)
  obj = IntruderOutputRef(handle)

  if outputVals != nothing
    set!(obj,outputVals)
  end
  return obj
end

function CollectionIntruderOutputRef(output::OutputRef, nintruders::Int)
  handle = ccall((:getRefCCollectionIntruderOutput, LIBCCAS), Ptr{Void},
                             (Ptr{Void},), output.handle)
  obj = CollectionIntruderOutputRef(handle)

  resize!(obj, nintruders)
  resize!(obj.intruders, nintruders)

  for i = 1:nintruders
    obj.intruders[i] = IntruderOutputRef(obj, i)
  end
  return obj
end

function set_id!(intruderOutput::IntruderOutputRef, id_::Integer)
  id = checked_convert(UInt32, id_)
  ccall((:setCIntrOutput_id, LIBCCAS), Void, (Ptr{Void}, UInt32), intruderOutput.handle, id)
end

function get!(output::IntruderOutputRef, outputVals::IntruderOutput)
  outputVals.id = ccall((:getCIntrOutput_id, LIBCCAS), UInt32, (Ptr{Void},), output.handle)
  outputVals.cvc = ccall((:getCIntrOutput_cvc, LIBCCAS), UInt8, (Ptr{Void},), output.handle)
  outputVals.vrc = ccall((:getCIntrOutput_vrc, LIBCCAS), UInt8, (Ptr{Void},), output.handle)
  outputVals.vsb = ccall((:getCIntrOutput_vsb, LIBCCAS), UInt8, (Ptr{Void},), output.handle)
  outputVals.tds = ccall((:getCIntrOutput_tds, LIBCCAS), Float64, (Ptr{Void},), output.handle)
  outputVals.code = ccall((:getCIntrOutput_code, LIBCCAS), UInt8, (Ptr{Void},), output.handle)

  return outputVals
end

function get!(output::OutputRef, outputVals::Output)
  outputVals.cc       = ccall((:getCOutput_cc, LIBCCAS), UInt8, (Ptr{Void},), output.handle)
  outputVals.vc       = ccall((:getCOutput_vc, LIBCCAS), UInt8, (Ptr{Void},), output.handle)
  outputVals.ua       = ccall((:getCOutput_ua, LIBCCAS), UInt8, (Ptr{Void},), output.handle)
  outputVals.da       = ccall((:getCOutput_da, LIBCCAS), UInt8, (Ptr{Void},), output.handle)
  outputVals.target_rate       = ccall((:getCOutput_target_rate, LIBCCAS), Float64, (Ptr{Void},), output.handle)
  outputVals.turn_off_aurals   = ccall((:getCOutput_turn_off_aurals, LIBCCAS), Bool, (Ptr{Void},), output.handle)
  outputVals.crossing          = ccall((:getCOutput_crossing, LIBCCAS), Bool, (Ptr{Void},), output.handle)
  outputVals.alarm    = ccall((:getCOutput_alarm, LIBCCAS), Bool, (Ptr{Void},), output.handle)
  outputVals.alert    = ccall((:getCOutput_alert, LIBCCAS), Bool, (Ptr{Void},), output.handle)
  outputVals.dh_min   = ccall((:getCOutput_dh_min, LIBCCAS), Float64, (Ptr{Void},), output.handle)
  outputVals.dh_max   = ccall((:getCOutput_dh_max, LIBCCAS), Float64, (Ptr{Void},), output.handle)
  outputVals.sensitivity_index = ccall((:getCOutput_sensitivity_index, LIBCCAS), UInt8, (Ptr{Void},), output.handle)
  outputVals.ddh      = ccall((:getCOutput_ddh, LIBCCAS), Float64, (Ptr{Void},), output.handle)

  for (out, outVals) in zip(output.intruder_collection.intruders, outputVals.intruders)
    get!(out, outVals)
  end

  return outputVals
end

function resize!(collection::CollectionIntruderOutputRef, size_::Integer)
  size = checked_convert(UInt32, size_)
  ccall((:resizeCCollectionIntrOutput, LIBCCAS), Void, (Ptr{Void}, UInt32), collection.handle, size)
end

function size(collection::CollectionIntruderOutputRef)
  csize = ccall((:sizeCCollectionIntrOutput, LIBCCAS), UInt32, (Ptr{Void},),
                 collection.handle)
  csize = convert(Int64,csize)
  @assert csize == length(collection.intruders) #make sure we're in sync

  return csize
end

type Constants
  handle::Ptr{Void} #Pointer to C object

  function Constants(quant_::Integer, config_filename::AbstractString, max_intruders_::Integer)
    quant = checked_convert(UInt8, quant_)
    max_intruders = checked_convert(UInt8, max_intruders_)
    handle = ccall((:newCConstants, LIBCCAS), Ptr{Void}, (UInt8, Ptr{UInt8}, UInt32),
                   quant, config_filename, max_intruders)
    obj = new(handle)

    finalizer(obj, obj->ccall((:delCConstants, LIBCCAS),Void, (Ptr{Void},), obj.handle))
    return obj
  end
end

type CASShared
  handle::Ptr{Void} #Pointer to C object
  libcas::AbstractString #path to library
  max_intruders::Int64
  input::InputRef
  output::OutputRef

  function CASShared(libcas::AbstractString, consts::Constants)

    handle = ccall((:newCCASShared, LIBCCAS), Ptr{Void}, (Ptr{Void}, Ptr{UInt8}),
                   consts.handle, libcas)

    obj = new(handle)
    obj.libcas = libcas
    obj.max_intruders = max_intruders(obj)
    obj.input = InputRef(obj.max_intruders)
    obj.output = OutputRef(obj.max_intruders)

    finalizer(obj,obj->ccall((:delCCASShared, LIBCCAS),Void, (Ptr{Void},),obj.handle))
    return obj
  end
end

reset(cas::CASShared) = ccall((:reset, LIBCCAS), Void, (Ptr{Void},), cas.handle)

version(cas::CASShared) = unsafe_string(ccall((:version, LIBCCAS), Ptr{UInt8}, (Ptr{Void},), cas.handle))

function error_msg(cas::CASShared)
  err = ccall((:error, LIBCCAS), Ptr{UInt8},(Ptr{Void},), cas.handle)

  return err == C_NULL ? nothing : unsafe_string(err)
end

max_intruders(cas::CASShared) = ccall((:max_intruders, LIBCCAS), Int64, (Ptr{Void},), cas.handle)

function update(cas::CASShared, input::InputRef, output::OutputRef)
  ccall((:update, LIBCCAS), Void,(Ptr{Void}, Ptr{Void}, Ptr{Void}), cas.handle, input.handle, output.handle)
end

function update!(cas::CASShared, inputVals::Input, outputVals::Output)

  #id in output isn't populated by libcas for some reason, so we'll do it on our side
  for i = 1:cas.max_intruders
    set_id!(cas.output.intruder_collection.intruders[i], inputVals.intruders[i].id)
  end

  set!(cas.input, inputVals)
  update(cas, cas.input,cas.output)
  get!(cas.output, outputVals)

  return outputVals #for convenience, since we are directly modifying outputVals input arg
end

function checked_convert{T <: Unsigned}(::Type{T},n::Integer)
  return (typemin(T) <= n <= typemax(T)) ?
    convert(T,n) : error("Cannot convert numeric type: Out of range")
end

author() = unsafe_string(ccall((:author, LIBCCAS), Ptr{UInt8},()))

end

max_stars_repo_name: sisl/CCAS

max_issues_repo_name: sisl/CCAS

max_forks_repo_name: sisl/CCAS


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base.Test

# test the basic floating point functions

# flipsign

for elty in (Float32,Float64)
    x = convert(elty,-2.0)
    x = flipsign(x,-1.0)
    @test flipsign(x,big(-1.0)) == convert(elty,-2.0)
end

# maxintfloat

@test maxintfloat(Float16) === Float16(2048f0)
for elty in (Float16,Float32,Float64)
    @test maxintfloat(rand(elty)) === maxintfloat(elty)
end
@test maxintfloat() === maxintfloat(Float64)
@test maxintfloat(Float64, Int32) === 2147483647.0
@test maxintfloat(Float32, Int32) === maxintfloat(Float32)
@test maxintfloat(Float64, Int16) === 32767.0
@test maxintfloat(Float64, Int64) === maxintfloat(Float64)

# isinteger
for elty in (Float16,Float32,Float64)
    @test !isinteger(elty(1.2))
    @test isinteger(elty(12))
    @test isinteger(zero(elty))
    @test isinteger(-zero(elty))
    @test !isinteger(nextfloat(zero(elty)))
    @test !isinteger(prevfloat(zero(elty)))
    @test isinteger(maxintfloat(elty))
    @test isinteger(-maxintfloat(elty))
    @test !isinteger(elty(Inf))
    @test !isinteger(-elty(Inf))
    @test !isinteger(elty(NaN))
end

# num2hex, hex2num
for elty in (Float16,Float32,Float64), _ = 1:10
    x = rand(elty)
    @test hex2num(num2hex(x)) ≈ x
end

# round
for elty in (Float32,Float64)
    x = rand(elty)
    A = fill(x,(10,10))
    @test round.(A,RoundToZero) == fill(trunc(x),(10,10))
    @test round.(A,RoundUp) == fill(ceil(x),(10,10))
    @test round.(A,RoundDown) == fill(floor(x),(10,10))
    A = fill(x,(10,10,10))
    @test round.(A,RoundToZero) == fill(trunc(x),(10,10,10))
    @test round.(A,RoundUp) == fill(ceil(x),(10,10,10))
    @test round.(A,RoundDown) == fill(floor(x),(10,10,10))
    for elty2 in (Int32,Int64)
        A = fill(x,(10,))
        @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,))
        @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,))
        @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,))
        A = fill(x,(10,10))
        @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,10))
        @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,10))
        @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,10))
        A = fill(x,(10,10,10))
        @test round.(elty2,A,RoundToZero) == fill(trunc(elty2,x),(10,10,10))
        @test round.(elty2,A,RoundUp) == fill(ceil(elty2,x),(10,10,10))
        @test round.(elty2,A,RoundDown) == fill(floor(elty2,x),(10,10,10))
        @test round.(elty2,A) == fill(round(elty2,x),(10,10,10))
    end
end

@testset "Types" begin
    for x in (Int16(0), 1, 2f0, pi, 3//4, big(5//6), 7.8, big(9), big(e))
        @test float(typeof(x)) == typeof(float(x))
        @test float(typeof(complex(x, x))) == typeof(float(complex(x, x)))
    end
end

ext: jl

max_stars_repo_name: jdlangs/julia

max_issues_repo_name: jdlangs/julia

max_forks_repo_name: jdlangs/julia


text:
function max_space_ids(ec::NestedEinsum{OMEinsum.DynamicEinCode{Int}}, ids_dic)
    OMEinsum.isleaf(ec) && return Int[]
    max_ids = Int[]
    max_sc = -Inf
    for arg in ec.args
        new_max_ids = max_space_ids(arg, ids_dic)
        new_sc = space_complexity(new_max_ids, ids_dic)
        if new_sc > max_sc
            max_ids = new_max_ids
            max_sc = new_sc
        end
    end
    for ids_x in ec.eins.ixs
        if space_complexity(ids_x, ids_dic) > max_sc
            max_ids = ids_x
            max_sc = space_complexity(ids_x, ids_dic)
        end
    end
    if space_complexity(ec.eins.iy, ids_dic) > max_sc
        max_ids = ec.eins.iy
        max_sc = space_complexity(ec.eins.iy, ids_dic)
    end
    return copy(max_ids)
end
space_complexity(ids::Vector{Int}, ids_dic) = isempty(ids) ? -Inf : sum(log2(ids_dic[i]) for i in ids)


text:
# KDE plot

# w/ Gadfly
# Pkg.add("KernelDensityEstimatePlotting")
# Pkg.add("KernelDensityEstimate")
# Pkg.add("Gadfly")
# Pkg.add("Plots")
# Pkg.add("PyPlot")

using Plots
using Gadfly

using KernelDensityEstimate, KernelDensityEstimatePlotting

p100 = kde!([randn(50);10.0.+2*randn(50)])
p2 = kde!([0.0;10.0],[1.0])
p75 = resample(p2,75)

Gadfly.plot([p100;p2;p75],c=["red";"green";"blue"])

# w/ Makie
Pkg.add("CairoMakie")
using CairoMakie
let
    μs = [0.25, 2, -1]
    σs = [1.25, 1, 2]
    colors = ["#E69F00", "#56B4E9", "#009E73"]
    x = LinRange(-10, 10, 200)
    fg(x, μ, σ) = exp.(.-(x .- μ) .^ 2 ./ (2σ^2)) ./ (σ * √(2π))
    fig = Figure(resolution = (600, 400), font = "CMU Serif")
    ax = Axis(fig[1, 1], xlabel = L"x", ylabel = L"y", xlabelsize = 22, ylabelsize = 22)
    for (idx, μ) in enumerate(μs)
        lines!(x, fg(x, μ, σs[idx]); color = colors[idx],
            label = L"\mu = %$(μ), \sigma = %$(σs[idx])")
        band!(x, fill(0, length(x)), fg(x, μ, σs[idx]); color = (colors[idx], 0.1),
            label = L"\mu = %$(μ), \sigma = %$(σs[idx])")
    end
    axislegend(; merge = true, framecolor = :red)
    display(fig)
end;

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
module ProgramConstantsModule

const MAX_DEGREE = 2
const CONST_TYPE = Float32
const BATCH_DIM = 2
const FEATURE_DIM = 1
const RecordType = Dict{String,Any}

"""Enum for concurrency type (to get function specialization)"""
abstract type SRConcurrency end
struct SRSerial <: SRConcurrency end
struct SRThreaded <: SRConcurrency end
struct SRDistributed <: SRConcurrency end

end


text:
module TestCRQ

using Base.Experimental: @sync
using ConcurrentCollections
using ConcurrentCollections.Implementations:
    CRQSlot, IndirectConcurrentRingQueueNode, trypush!, isclosed
using ProgressLogging: @logprogress, @withprogress
using Test
using ..Utils: ⊏

function test_CRQSlot()
    for index in [111, 222],
        safe in [false, true],
        storage in UInt32[0xaaa, 0xbbb]

        @test NamedTuple(CRQSlot(; index, safe, storage)) == (; index, safe, storage)
    end
end

function test_print_crqslot()
    slot = CRQSlot(; index = 111, safe = false, storage = UInt32(0xaaa))
    str = sprint(show, "text/plain", slot)
    @test "CRQSlot" ⊏ str
    @test r"index *= *111" ⊏ str
    @test r"safe *= *false" ⊏ str
    @test r"storage *= *0x0+aaa" ⊏ str
end

function unfair_sleep(seconds::Real)
    t0 = time_ns()
    ns = seconds * 1e9
    while time_ns() - t0 < ns
        GC.safepoint()
    end
end

function concurrent_denqueue!(
    crq::IndirectConcurrentRingQueueNode,
    nitems::Integer,
    nsend::Integer,
    nrecv::Integer,
)
    @assert nsend > 0
    @assert nrecv > 0
    use_yield = nsend + nrecv > Threads.nthreads()
    # TODO: If `!use_yield`, make sure the tasks are spawned in different
    # threads and they are sticky.
    function poll()
        if use_yield
            yield()
        else
            GC.safepoint()
        end
        rand() < 0.1 && unfair_sleep(rand() * 1e-6)
    end

    received = Vector{Int}[]
    senders = Task[]
    receivers = Task[]
    global SENDERS = senders
    global RECEIVERS = receivers
    ref = Threads.Atomic{Int}(0)
    activesenders = Threads.Atomic{Int}(nsend)
    @sync begin
        for offset in 1:nsend
            t = Threads.@spawn try
                local y = nothing
                local i = 0
                local allpushed = true
                for outer i in offset:nsend:nitems
                    local x = eltype(crq)(i)
                    if !trypush!(crq, x)
                        allpushed = false
                        break
                    end
                    local s = Threads.atomic_add!(ref, 1)
                    if s > nsend
                        while ref[] > 1
                            poll()
                        end
                    end
                end
                return (; y, i, allpushed)
            finally
                Threads.atomic_sub!(activesenders, 1)
            end
            push!(senders, t)
        end
        for _ in 1:nrecv
            ys = Int[]
            push!(received, ys)
            t = Threads.@spawn begin
                local ys_nb = Int[]
                local ys_b = Int[]
                while true
                    local y = maybepopfirst!(crq)
                    if y === nothing
                        if activesenders[] == 0
                            y = maybepopfirst!(crq)
                            if y === nothing  # Confirm that CRQ is empty
                                return y
                            end
                            # Reaching here means that there were some enqueues
                            # between our `maybepopfirst!(crq) === nothing` and
                            # `activesenders[] == 0`.
                        else
                            poll()
                            continue
                        end
                    end
                    local i = something(y)
                    push!(ys, i)

                    local s = Threads.atomic_sub!(ref, 1)
                    if s < -nrecv
                        while ref[] < -1
                            poll()
                        end
                    end
                end
            end
            push!(receivers, t)
        end
    end
    return received, senders, receivers
end

function check_consecutive(xs)
    notfound = Int[]
    dups = Int[]
    pre = xs[begin] - 1
    for x in xs
        e = pre + 1
        append!(notfound, e:x-1)
        append!(dups, x:e-1)
        pre = x
    end
    return (; notfound, dups)
end

function test_concurrent_push_pop(ntrials = 100)
    @withprogress name = "concurrent push-pop" begin
        @testset for trial in 1:ntrials
            @logprogress (trial - 1) / ntrials
            check_concurrent_push_pop()
        end
    end
end

function check_concurrent_push_pop()
    global received, notfound, dups, allreceived
    nsend = cld(Threads.nthreads(), 2)
    nrecv = max(1, Threads.nthreads() - nsend)
    crq = IndirectConcurrentRingQueueNode{Int16}(32)
    global CRQ = crq
    nitems = 2^20
    nitems = typemax(Int16)
    received, senders, receivers = concurrent_denqueue!(crq, nitems, nsend, nrecv)
    allreceived = reduce(vcat, received)

    if isclosed(crq)
        @info "CRQ closed. Skipping the tests..."
        return
    end

    @test [fetch(t).allpushed for t in senders] == fill(true, nsend)

    sort!(allreceived)
    (; notfound, dups) = check_consecutive(allreceived)
    @test length(allreceived) == nitems
    @test notfound == []
    @test dups == []
    @test allreceived == 1:nitems
end

function test_print()
    crq = IndirectConcurrentRingQueueNode{Int16}(32)
    str = sprint(show, "text/plain", crq)
    @test "CRQ: " ⊏ str
end

end  # module


text:
models = (
    LogisticClassifier,
    LogisticCVClassifier,
    PassiveAggressiveClassifier,
    PerceptronClassifier,
    RidgeClassifier,
    RidgeCVClassifier,
    SGDClassifier,
    ProbabilisticSGDClassifier,
    )

fparams = (
    LogisticClassifier=(:classes, :coef, :intercept),
    LogisticCVClassifier=(:classes, :coef, :intercept, :Cs, :l1_ratios, :coefs_paths, :scores, :C, :l1_ratio),
    PassiveAggressiveClassifier=(:coef, :intercept),
    PerceptronClassifier=(:coef, :intercept),
    RidgeClassifier=(:coef, :intercept),
    RidgeCVClassifier=(:coef, :intercept),
    SGDClassifier=(:coef, :intercept),
    ProbabilisticSGDClassifier=(:coef, :intercept)
    )

@testset "Fit/Predict" begin
    X, y = simple_binaryclf()
    for mod in models
        m = mod()
        acc, f = test_clf(m, X, y)
        @test acc >= 0.75
        @test keys(MB.fitted_params(m, f)) == getproperty(fparams, Symbol(mod))
    end
end

@testset "Docstrings" begin
    for mod in models
        m = mod()
        @test !isempty(MB.docstring(m))
    end
end


text:
using JuMP, EAGO, Gurobi, CSV, DataFrames, Ipopt

m = Model(with_optimizer(
    EAGO.Optimizer,
    relaxed_optimizer = Gurobi.Optimizer(OutputFlag = 0),
))

# ----- Variables ----- #
x_Idx = Any[1, 2]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[1], -10.0)
JuMP.set_upper_bound(x[1], 10.0)
JuMP.set_lower_bound(x[2], -10.0)
JuMP.set_upper_bound(x[2], 10.0)

# ----- Objective ----- #
@NLobjective(
    m,
    Min,
    ((-5 + sqrt((x[1])^2)) * (-5 + sqrt((x[1])^2)) +
     (-5 + sqrt((x[2])^2)) * (-5 + sqrt((x[2])^2)))
)

JuMP.optimize!(m)
val = objective_value(m)
bnd = objective_bound(m)
gap = bnd - val
backend_opt = backend(m).optimizer.model.optimizer
last_relaxed_opt = backend_opt.relaxed_optimizer
run_time = backend_opt._run_time
println("run time: $run_time")


text:
###############################################################################
#
# The data used through the algorithmic process in the Stopping framework
# are stored in a State.
# We illustrate here the NLPAtX which is a specialization of the State for
# non-linear programming.
#
###############################################################################
#using Test, NLPModels, Stopping

include("../test-stopping/rosenbrock.jl")
#Formulate the problem with NLPModels
x0 = ones(6)
y0 = ones(1)
c(x) = [x[1] - x[2]]
lcon = [0.0]
ucon = [0.0]

#We can create a NLPAtX for constrained optimization.
#Here we provide y0 = [1.0]
#Note that the default value is [0.0]
nlp = ADNLPModel(x->rosenbrock(x), x0, y0 = y0,
                 c=c, lcon=lcon, ucon=ucon,
                 lvar=zeros(6), uvar = Inf * ones(6))
#We can create a NLPAtX for bounds-constrained optimization:
nlp2 = ADNLPModel(x->rosenbrock(x), x0,
                 lvar=zeros(6), uvar = Inf * ones(6))
#We can create a NLPAtX for unconstrained optimization:
nlp3 = ADNLPModel(x->rosenbrock(x), x0)

###############################################################################
#I. Initialize a NLPAtX:
#
#There are two main constructor for the States.
#The unconstrained:
state_unc = NLPAtX(x0)
#The constrained:
state_con = NLPAtX(x0, y0)

#By default, all the values in the State are set to nothing except x and lambda
#In the unconstrained case lambda is a vector of length 0
@test !(state_unc.lambda == nothing)
#From the default initialization, all the other entries are void:
@test state_unc.mu == nothing && state_con.mu == nothing
@test state_unc.fx == nothing && state_con.fx == nothing
#exception is the counters which is initialized as a default Counters:
@test (sum_counters(state_unc.evals) + sum_counters(state_con.evals)) == 0

#Note that the constructor proceeds to a size checking on gx, Hx, mu, cx, Jx.
#It returns an error if this test fails.
try
  NLPAtX(x0, Jx = ones(1,1))
  @test false
catch
  #printstyled("NLPAtX(x0, Jx = ones(1,1)) is invalid as length(lambda)=0\n")
  @test true
end

###############################################################################
#II. Update the entries
#
#At the creation of a NLPAtX, keyword arguments populate the state:
state_bnd = NLPAtX(x0, mu = zeros(6))
@test state_bnd.mu == zeros(6) #initialize multipliers with bounds constraints

#The NLPAtX has two functions: update! and reinit!
#The update! has the same behavior as in the GenericState:
update!(state_bnd, fx = 1.0, blah = 1) #update! ignores unnecessary keywords
@test state_bnd.mu == zeros(6) && state_bnd.fx == 1.0 && state_bnd.x == x0

#reinit! by default reuse x and lambda and reset all the entries at their
#default values (void or empty Counters):
reinit!(state_bnd, mu = ones(6))
@test state_bnd.mu == ones(6) && state_bnd.fx == nothing
@test state_bnd.x == x0 && state_bnd.lambda == zeros(0)
#Trying to inherit reinit!(AbstractState, Vector) would not work here as
#lambda is a mandatory entry.
try
 reinit!(state_bnd, 2 * ones(6))
 @test false
catch
 @test true
end
#However, we can specify both entries
reinit!(state_bnd, 2 * ones(6), zeros(0))
@test state_bnd.x == 2*ones(6) && state_bnd.lambda == zeros(0)
@test state_bnd.mu == nothing && sum_counters(state_bnd.evals) == 0
#Giving a new Counters update as well:
test = Counters(); setfield!(test, :neval_obj, 102)
reinit!(state_bnd, evals = test)
@test getfield(state_bnd.evals, :neval_obj) == 102
@test sum_counters(state_bnd.evals) - 102 == 0

###############################################################################
#III. Domain Error
#Similar to the GenericState we can use _domain_check to verify there are no NaN
@test Stopping._domain_check(state_bnd) == false
update!(state_bnd, fx = NaN)
@test Stopping._domain_check(state_bnd) == true

###############################################################################
#IV. Use the NLPAtX
#
#For algorithmic use, it might be conveninent to fill in all the entries of then
#State. In this case, we can use the Stopping:
stop = NLPStopping(nlp, (x,y) -> unconstrained_check(x,y), state_unc)
#Note that the fill_in! can receive known informations via keywords.
#If we don't want to store the hessian matrix, we turn the keyword
#matrix_info as false.
fill_in!(stop, x0, matrix_info = false)

#printstyled("Hx has not been updated: ",stop.current_state.Hx == nothing,"\n")
@test stop.current_state.Hx == nothing

# We can now use the updated step in the algorithmic procedure
@test start!(stop) #return true


text:
"""
	gmtconnect(cmd0::String="", arg1=nothing, kwargs...)

Connect individual lines whose end points match within tolerance

Full option list at [`gmtconnect`]($(GMTdoc)gmtconnect.html)

Parameters
----------

- **C** | **closed** :: [Type => Str | []]        `Arg = [closed]`

    Write all the closed polygons to closed [gmtgmtconnect_closed.txt] and return all other
    segments as they are. No gmtconnection takes place.
    ($(GMTdoc)gmtconnect.html#c)
- **D** | **dump** :: [Type => Str | []]   `Arg = [template]`

    For multiple segment data, dump each segment to a separate output file
    ($(GMTdoc)gmtconnect.html#d)
- **L** | **linkfile** :: [Type => Str | []]      `Arg = [linkfile]`

    Writes the link information to the specified file [gmtgmtconnect_link.txt].
    ($(GMTdoc)gmtconnect.html#l)
- **Q** | **list_file** :: [Type => Str | []]      `Arg =  [listfile]`

    Used with **D** to write a list file with the names of the individual output files.
    ($(GMTdoc)gmtconnect.html#q)
- **T** | **tolerance** :: [Type => Str | List]    `Arg = [cutoff[unit][/nn_dist]]`

    Specifies the separation tolerance in the data coordinate units [0]; append distance unit.
    If two lines has end-points that are closer than this cutoff they will be joined.
    ($(GMTdoc)gmtconnect.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function gmtconnect(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("gmtconnect", cmd0, arg1)

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:V_params :b :d :e :f :g :h :i :o :yx])
	cmd  = parse_these_opts(cmd, d, [[:C :closed], [:D :dump], [:L :linkfile], [:Q :list_file], [:T :tolerance]])

	common_grd(d, cmd0, cmd, "gmtconnect ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtconnect(arg1, arg2=nothing, cmd0::String=""; kw...) = gmtconnect(cmd0, arg1, arg2; kw...)


text:
using Base.REPLCompletions

module CompletionFoo
    module CompletionFoo2

    end
    const bar = 1
    foo() = bar
end

test_complete(s) = completions(s,endof(s))
test_scomplete(s) = shell_completions(s,endof(s))

s = ""
c,r = test_complete(s)
@test in("CompletionFoo",c)
@test isempty(r)
@test s[r] == ""

s = "Comp"
c,r = test_complete(s)
@test in("CompletionFoo",c)
@test r == 1:4
@test s[r] == "Comp"

s = "Main.Comp"
c,r = test_complete(s)
@test in("CompletionFoo",c)
@test r == 6:9
@test s[r] == "Comp"

s = "Main.CompletionFoo."
c,r = test_complete(s)
@test in("bar",c)
@test r == 20:19
@test s[r] == ""

s = "Main.CompletionFoo.f"
c,r = test_complete(s)
@test in("foo",c)
@test r == 20:20
@test s[r] == "f"

# issue #6333
s = "Base.return_types(getin"
c,r = test_complete(s)
@test in("getindex",c)
@test r == 19:23
@test s[r] == "getin"

## Test completion of packages
#mkp(p) = ((@assert !isdir(p)); mkdir(p))
#temp_pkg_dir() do
#    mkp(Pkg.dir("MyAwesomePackage"))
#    mkp(Pkg.dir("CompletionFooPackage"))
#
#    s = "using MyAwesome"
#    c,r = test_complete(s)
#    @test in("MyAwesomePackage",c)
#    @test s[r] == "MyAwesome"
#
#    s = "using Completion"
#    c,r = test_complete(s)
#    @test in("CompletionFoo",c) #The module
#    @test in("CompletionFooPackage",c) #The package
#    @test s[r] == "Completion"
#
#    s = "using CompletionFoo.Completion"
#    c,r = test_complete(s)
#    @test in("CompletionFoo2",c)
#    @test s[r] == "Completion"
#
#    rmdir(Pkg.dir("MyAwesomePackage"))
#    rmdir(Pkg.dir("CompletionFooPackage"))
#end

@unix_only begin
    #Assume that we can rely on the existence and accessibility of /tmp
    s = "/t"
    c,r = test_scomplete("/t")
    @test in("tmp/",c)
    @test r == 2:2
    @test s[r] == "t"

    s = "/tmp"
    c,r = test_scomplete(s)
    @test in("tmp/",c)
    @test r == 2:4
    @test s[r] == "tmp"

    # This should match things that are inside the tmp directory
    if !isdir("/tmp/tmp")
        s = "/tmp/"
        c,r = test_scomplete(s)
        @test !in("tmp/",c)
        @test r == 6:5
        @test s[r] == ""
    end

    s = "cd \$(Pk"
    c,r = test_scomplete(s)
    @test in("Pkg",c)
    @test r == 6:7
    @test s[r] == "Pk"
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
#
# Test a sparse rosenbrock, i.e. we have many factors but only a few are
# actually active in determining the response from the function.
#
function rosenbrock(x)
  n = length(x)
  return sum( 100*( x[2:n] - x[1:(n-1)].^2 ).^2 + ( x[1:(n-1)] - 1 ).^2 )
end
# FIXME this should be a more efficient implementation, but switching to it would reset all benchmarks
#rosenbrock(x) = sum(i -> 100*abs2(x[i+1] - x[i]^2) + abs2(x[i] - 1), Base.OneTo(length(x)-1))

function sphere(x)
  sum(x.^2)
end
# FIXME this should be a more efficient implementation, but switching to it would reset all benchmarks
#sphere(x) = sum(abs2, x)

mutable struct SparseShiftedProblem
  active_factors::Array{Int, 1}
  shift::Array{Float64, 2}
  func::Function

  function SparseShiftedProblem(n = 1000, k = 20, baseProblem = rosenbrock)
    afs = shuffle(collect(1:n))[1:k]
    new(afs, randn(k,1), baseProblem)
  end
end

function calc(sr::SparseShiftedProblem, x)
  sr.func(x[sr.active_factors] - sr.shift)
end

function sample_points_and_calc_response(n, s, func)
  xs = randn(n, s);
  y = zeros(s)
  [(y[i] = func(xs[:,i])) for i in 1:s]
  return xs, y
end

# Correlation selection is based on correlating each covariate with the response
# and keeping the covariates with the highest correlations with the response.
function correlation_selection(xs, y, maxSelected)
  correlations = y' * xs'
  find_top_largest_correlations = sort(correlations[:], rev=true)[1:maxSelected]
  index_of_largest_correlations = findall(in(find_top_largest_correlations), correlations)
  index_of_largest_correlations
end

# Absolute Correlation selection is based on correlating each covariate with the response
# and keeping the covariates with the highest absolute correlations with the response.
function abs_correlation_selection(xs, y, maxSelected)
  correlations = abs(y' * xs')
  find_top_largest_correlations = sort(correlations[:], rev=true)[1:maxSelected]
  index_of_largest_correlations = findall(in(find_top_largest_correlations), correlations)
  index_of_largest_correlations
end

# Stepwise filtering abs correlation selection. This is based on the idea that if
# there are complex interactive effects among the covariates the covariates that
# are part of an interaction should also have some individual effect. This is similar
# to the assumption of "effect heredity". Thus we do a stepwise filtering of the possible
# covariates to include by first eliminating half of them, then redoing the filtering
# on the interactions between the first selection, and so on.

# Random selection.
function random_selection(xs, y, maxSelected)
  n = size(xs, 1)
  shuffle(collect(1:n))[1:maxSelected]
end

# Call everything once to ensure they have been compiled.
a = rosenbrock([10.0, 4.5])
b = sphere([10.0, 4.5])
p = SparseShiftedProblem(100, 5, rosenbrock)
xs, y = sample_points_and_calc_response(100, 10, (x) -> calc(p, x))
s = random_selection(xs, y, 5)
s2 = correlation_selection(xs, y, 5)

# A problem is characterized mainly by:
#   1. How complex the interrelations are among the active covariates, basically determined by the objective func itself
#   2. Sparseness degree, i.e. how many of all available covariates that actually affects the response (0.0-1.0)
#   3. Dimensionality, the total number of covariates (100-10000)
#
# A selection scheme then saves an incomplete "trail" of evaluations of the objective function
# and uses them to predict which covariates are important. The size of this set is also
# important, as well as how it evolves over time.
#
# So we want to study how well a covariate subset selection (CSS) scheme predicts
# which covariates are important as we vary the four main parameters above.
# Dependent variables:
#   O1. Percentage of missed active covariates over NR repetitions
#   O2. Variation in missed covariates over NR repetitions
#   O3. Average time for a prediction
# Independent variables:
Problems = [rosenbrock, sphere]
Dimensionality = [1000, 10000, 100000]
ActualSparsenessDegree = [0.001, 0.01, 0.10, 0.25]
ArchiveSizeDegree = [0.1, 0.50, 1.00]
GuessedSparsenessDegree = [0.001, 0.01, 0.05, 0.10]

NumRepeats = 10
CovariateSelectionSchemes = [abs_correlation_selection, correlation_selection]
ncss = length(CovariateSelectionSchemes)

csvfile = open("comparing_covariate_selection_schemes.csv", "w")

using HypothesisTests

# Return the name of the best method of two methods given a Wilcoxon test.
# By default we are minimizing so better means smaller average, set minimizing
# to false if maximizing.
function best_of_if_better(x, y, names = ["1", "2"], minimizing = true, pvalue_treshold = 0.05)
  if pvalue(MannWhitneyUTest(x, y)) < pvalue_treshold
    op = minimizing ? < : >
    op(mean(x), mean(y)) ? names[1] : names[2]
  else
    "no_diff"
  end
end

println(csvfile, "CovariateSelectionScheme,Problem,Dimensionality,ActualSparseness,GuessedSparseness,ArchiveSize,MeanTime,StdMissed,MeanMissed")
for problem in Problems
  for dim in Dimensionality
    for asd in ActualSparsenessDegree
      for sizedegree in ArchiveSizeDegree
        for gsd in GuessedSparsenessDegree

          k = int(asd * dim)

          num_selected = int(gsd * dim)
          percent_missed = zeros(NumRepeats, ncss)
          times = zeros(NumRepeats, ncss)
          size = int(sizedegree * dim)

          for ir in 1:NumRepeats
            ssp = SparseShiftedProblem(dim, k, problem)
            xs, y = sample_points_and_calc_response(dim, size, (x) -> calc(ssp, x))

            for cssi in 1:ncss
              css = CovariateSelectionSchemes[cssi]
              tic()
              selected_factors = css(xs, y, num_selected)
              times[ir, cssi] = toq()
              missed = setdiff(Set(ssp.active_factors...), Set(selected_factors...))
              percent_missed[ir, cssi] = length(missed)/k*100.0
            end
          end

          for cssi in 1:ncss
            css = CovariateSelectionSchemes[cssi]
            println(csvfile, css, ",", problem, ",", dim, ",", asd, ",", gsd, ",", size, ",",
              mean(times[:,cssi]), ",",
              std(percent_missed[:,cssi]), ",", mean(percent_missed[:,cssi]))
            flush(stdout)
          end
          ms = mean(percent_missed, 1) # ms = mean(percent_missed, 1)
          best_method = best_of_if_better(percent_missed[:,1], percent_missed[:,2], CovariateSelectionSchemes)
          println(join([problem, dim, asd, gsd, size, minimum(ms), ms[2]-ms[1], best_method], ","))
        end
      end
    end
  end
end

# The conclusions from all these runs are:
#  For both problems:
#    1. Absolute correlation selection is almost always better than correlation selection (as expected)
#  For Sphere:
#    1. Only if the actual sparseness degree is very low (i.e. with < 30 actual active covariates) are the methods very good
#    2. And even so one must guess at least 10-50 times more active covariates
#  For Rosenbrock:
#    1. 1000: Only with an intermediate actual sparseness (0.01) are the methods very good
#    2. 10000: Only with very low actual sparseness (0.001) are the methods very good
#    3. 10000: Increasing the number of samples has quite a large effect i.e. lowers the missed percentage
#
# What is not clear is how good the method needs to be to actually give an advantage since they
# are applied in rounds. Maybe a small improvement can have quite a large effect in speeding up convergence.
#
# OTOH, we already have information about the correlation of different covariates in the diagonal
# of the covar matrix so maybe these additional selection algs are not needed...


text:
"""
    function addUpdateLocal_!(p::Program_,platform::String)

Generate the functions related with Global Updates.
"""
function addUpdateGlobal_!(p::Program_,platform::String)

    if "UpdateGlobal" in keys(p.agent.declaredUpdates)

        #Construct functions
        f = simpleFirstLoopWrapInFunction_(platform,:globStep_!,
                        :(begin
                            if ic1_ == 1
                                $(p.agent.declaredUpdates["UpdateGlobal"])
                            end
                        end)
                        )
        f = vectorize_(p.agent,f,p)

        push!(p.declareF.args,
            f)

        push!(p.execInloop.args,
                :(@platformAdapt globStep_!(ARGS_))
            )
    end

    return nothing
end


text:
"""
An aggregator interface for SSA-like algorithms.

### Required Fields
- `next_jump`
- `next_jump_time`
- `end_time`
- `cur_rates`
- `sum_rate`
- `ma_jumps`
- `rates`
- `affects!`
- `save_positions`
- `rng`
"""
abstract type AbstractSSAJumpAggregator <: AbstractJumpAggregator end

DiscreteCallback(c::AbstractSSAJumpAggregator) = DiscreteCallback(c, c, initialize = c, save_positions = c.save_positions)

########### The following routines should be templates for all SSAs ###########

# # condition for jump to occur
# @inline function (p::SSAJumpAggregator)(u, t, integrator)
#   p.next_jump_time == t
# end

# # executing jump at the next jump time
# function (p::SSAJumpAggregator)(integrator)
#   execute_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
#   generate_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
#   register_next_jump_time!(integrator, p, integrator.t)
#   nothing
# end

# # setting up a new simulation
# function (p::SSAJumpAggregator)(dj, u, t, integrator) # initialize
#   initialize!(p, integrator, u, integrator.p, t)
#   register_next_jump_time!(integrator, p, integrator.t)
#   nothing
# end

############################## Generic Routines ###############################

@inline function register_next_jump_time!(integrator, p::AbstractSSAJumpAggregator, t)
    if p.next_jump_time < p.end_time
        add_tstop!(integrator, p.next_jump_time)
    end
    nothing
end

# helper routine for setting up standard fields of SSA jump aggregations
function build_jump_aggregation(jump_agg_type, u, p, t, end_time, ma_jumps, rates,
                                affects!, save_positions, rng; kwargs...)

  # mass action jumps
    majumps = ma_jumps
    if majumps === nothing
        majumps = MassActionJump(Vector{typeof(t)}(),
                             Vector{Vector{Pair{Int,eltype(u)}}}(),
                             Vector{Vector{Pair{Int,eltype(u)}}}())
    end

  # current jump rates, allows mass action rates and constant jumps
    cur_rates = Vector{typeof(t)}(undef, get_num_majumps(majumps) + length(rates))

    sum_rate = zero(typeof(t))
    next_jump = 0
    next_jump_time = typemax(typeof(t))
    jump_agg_type(next_jump, next_jump_time, end_time, cur_rates, sum_rate,
                majumps, rates, affects!, save_positions, rng; kwargs...)
end

# reevaluate all rates and total rate
function fill_rates_and_sum!(p::AbstractSSAJumpAggregator, u, params, t)
    sum_rate = zero(typeof(p.sum_rate))

  # mass action jumps
    majumps   = p.ma_jumps
    cur_rates = p.cur_rates
    @inbounds for i in 1:get_num_majumps(majumps)
        cur_rates[i] = evalrxrate(u, i, majumps)
        sum_rate    += cur_rates[i]
    end

  # constant rates
    rates = p.rates
    idx   = get_num_majumps(majumps) + 1
    @inbounds for rate in rates
        cur_rates[idx] = rate(u, params, t)
        sum_rate += cur_rates[idx]
        idx += 1
    end

    p.sum_rate = sum_rate
    nothing
end

# recalculate jump rates for jumps that depend on the just executed jump
# requires dependency graph
function update_dependent_rates!(p::AbstractSSAJumpAggregator, u, params, t)
    @inbounds dep_rxs = p.dep_gr[p.next_jump]
    cur_rates   = p.cur_rates
    sum_rate    = p.sum_rate
    @inbounds for rx in dep_rxs
        sum_rate -= cur_rates[rx]
        @inbounds cur_rates[rx] = calculate_jump_rate(p.ma_jumps, p.rates, u,params,t,rx)
        sum_rate += cur_rates[rx]
    end

    p.sum_rate = sum_rate
    nothing
end

# Update state based on the p.next_jump
@inline function update_state!(p::AbstractSSAJumpAggregator, integrator, u)
    @unpack ma_jumps, next_jump = p
    num_ma_rates = get_num_majumps(ma_jumps)
    if next_jump <= num_ma_rates # is next jump a mass action jump
        if u isa SVector
            integrator.u = executerx(u, next_jump, ma_jumps)
        else
            @inbounds executerx!(u, next_jump, ma_jumps)
        end
    else
        idx = next_jump - num_ma_rates
        @inbounds p.affects![idx](integrator)
    end
    return integrator.u
end

"check if the total rate is 0 and if it is, make the next jump time Inf"
@inline function nomorejumps!(p, sum_rate) :: Bool
    if sum_rate < eps(typeof(sum_rate))
        p.next_jump = 0
        p.next_jump_time = convert(typeof(sum_rate), Inf)
        return true
    end
    return false
end

"perform linear search of r over array. Output element j s.t. array[j-1] < r <= array[j]. Will error if no such r exists"
@inline function linear_search(array, r)
    jidx = 1
    @inbounds parsum = array[jidx]
    while parsum < r
        jidx   += 1
        @inbounds parsum += array[jidx]
    end
    return jidx
end

"perform rejection sampling test"
@inline function rejectrx(ma_jumps, rates, cur_rate_high, cur_rate_low, rng, u, jidx, params, t)
    # rejection test
    @inbounds r2     = rand(rng) * cur_rate_high[jidx]
    @inbounds crlow  = cur_rate_low[jidx]

    if crlow > zero(crlow) && r2 <= crlow
        return false
    else
        # calculate actual propensity, split up for type stability
        crate = calculate_jump_rate(ma_jumps, rates, u, params, t, jidx)
        if crate > zero(crate) && r2 <= crate
            return false
        end
    end
    return true
end

"update the jump rate, assuming p.rates is a vector of functions"
@inline function calculate_jump_rate(ma_jumps, rates, u, params, t, rx)
    num_majumps = get_num_majumps(ma_jumps)
    if rx <= num_majumps
        return evalrxrate(u, rx, ma_jumps)
    else
        @inbounds return rates[rx - num_majumps](u, params, t)
    end
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ b17843a0-05ae-11ec-1825-8b2b03bd5811
begin
    using Pkg
    Pkg.activate("..")
    using Revise
    using Javis
    using JUtils
    using Colors
    using Random
    using ROCAnalysis
    using GLM
    using DataFrames
    using StatsPlots
    # import Luxor
    import MLJ: make_regression
end

# ╔═╡ 5e3a6b7d-f15d-4541-9b48-7c3c3ce5ddd4
begin
    X, y = make_regression(
        100,
        1,
        noise = 0.0,
        sparse = 0.0,
        outliers = 0.0,
        binary = true,
        as_table = false,
    )
    df = DataFrame(x = (X[:, 1]), y = y)
    log = glm(@formula(y ~ x), df, Binomial())
    preds = predict(log)
    tar = preds[df.y .== 1.0]
    nontar = preds[df.y .== 0.0]
    myroc = roc(tar, nontar)

    myplot = plot(myroc.pmiss, 1 .- myroc.pfa, label = "")
    plot!(myplot, [0, 1], [0, 1], label = "")
end

# ╔═╡ e0bd9a55-5e4b-4d13-adf6-24ccd11a0384
@draw begin
    draw_roc(Point(100, 0), tar, nontar, 0.05)
end 700 500

# ╔═╡ f9061057-61fd-492f-9624-4aa44db95cdd
begin
    function ground(args...)
        background("white")
        sethue("black")
    end
end

# ╔═╡ 61a37185-92fe-4b31-b609-43167b9c3f76
begin
    scatter(sort(nontar), bins = 100)
    scatter!(sort(tar), bins = 100)
end

# ╔═╡ 54de45e1-f204-4521-a31c-3ce5e2e76a7d
begin
    n_frames = 200
    n_steps = 20
    step_size = n_frames ÷ n_steps

    size = 300
    global_shift = Point(100, 50)

    my_video = Video(800, 800)

    Background(1:n_frames, ground)

    max_pred = maximum(preds)

    for (frame_n, p) in zip(1:step_size:n_frames, range(0.0, 1.0, length = n_steps))
        Object(
            frame_n:(frame_n + step_size - 1),
            @JShape begin
                draw_roc(global_shift, tar, nontar, p)
            end
        )
    end

    margin = 50

    Object(
        1:n_frames,
        @JShape begin
            rocpoints = Point.(myroc.pmiss .* size, myroc.pfa .* size)
            Javis.translate(-midpoint(rocpoints[1], rocpoints[end]))
            bottomleft = Point(rocpoints[1].x, rocpoints[end].y - margin)
            topright = Point(rocpoints[end].x, bottomleft.y - margin - max_pred * size / 2)
            Javis.translate(global_shift)
            newidxs = sortperm(preds)
            barchart(
                preds[newidxs],
                boundingbox = Luxor.BoundingBox([bottomleft, topright]),
                margin = 0,
                barfunction = (values, i, lowpos, highpos, barwidth, scaledvalue) -> begin
                    @layer begin
                        height = (topright.y - bottomleft.y)
                        y[newidxs][i] ? sethue("blue") : sethue("red")
                        circle(Point(lowpos.x, lowpos.y + values[i] * height), 3, :fill)
                    end
                end,
            )
        end
    )

    for (frame_n, p) in zip(1:step_size:n_frames, range(0.0, 1.0, length = n_steps))
        Object(
            frame_n:(frame_n + step_size - 1),
            @JShape begin
                rocpoints = Point.(myroc.pmiss .* size, myroc.pfa .* size)
                Javis.translate(-midpoint(rocpoints[1], rocpoints[end]))
                Javis.translate(global_shift)
                bottomleft = Point(rocpoints[1].x, rocpoints[end].y - margin)
                topright =
                    Point(rocpoints[end].x, bottomleft.y - max_pred * size / 2 - margin)

                bottomright = Point(topright.x, bottomleft.y)

                line(bottomleft, Point(bottomleft.x, topright.y - margin))
                line(bottomleft, Point(topright.x + margin, bottomleft.y))

                height = (topright.y - bottomleft.y)
                shift = p * height
                line(bottomleft + (0, shift), bottomright + (0, shift), :stroke)
                setcolor(sethue("blue")..., 0.4)
                box(bottomleft, bottomright + (0, shift), :fill)
                setcolor(sethue("red")..., 0.4)
                box(bottomleft + (0, shift), bottomright + (0, height), :fill)
            end
        )
    end

    render(my_video, pathname = "../output/roc_animation.gif")
end

# ╔═╡ Cell order:
# ╠═b17843a0-05ae-11ec-1825-8b2b03bd5811
# ╠═5e3a6b7d-f15d-4541-9b48-7c3c3ce5ddd4
# ╠═e0bd9a55-5e4b-4d13-adf6-24ccd11a0384
# ╟─f9061057-61fd-492f-9624-4aa44db95cdd
# ╠═61a37185-92fe-4b31-b609-43167b9c3f76
# ╠═54de45e1-f204-4521-a31c-3ce5e2e76a7d


text:
module MHDSolver
# ----------
# Navier–Stokes Solver for 3D Magnetohydrodynamics Problem
# ----------

export
	UᵢUpdate!,
	BᵢUpdate!,
	MHDcalcN_advection!,
	MHDupdatevars!

using
  CUDA,
  Reexport,
  DocStringExtensions

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum

# δ function
δ(a::Int,b::Int) = ( a == b ? 1 : 0 );

function UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")

  if direction == "x"

  	# a = {1,2,3} -> {x,y,z} direction
  	a    = 1;
  	kₐ   = grid.kr;
  	k⁻²  = grid.invKrsq;
  	uᵢh  = vars.uxh;
  	∂uᵢh∂t = @view N[:,:,:,params.ux_ind];

  elseif direction == "y"

  	a    = 2;
  	kₐ   = grid.l;
  	k⁻²  = grid.invKrsq;
  	uᵢh  = vars.uyh;
  	∂uᵢh∂t = @view N[:,:,:,params.uy_ind];

  elseif direction == "z"

  	a    = 3;
  	kₐ   = grid.m;
  	k⁻²  = grid.invKrsq;
  	uᵢh  = vars.uzh;
  	∂uᵢh∂t = @view N[:,:,:,params.uz_ind];

  else

  	@warn "Warning : Unknown direction is declerad"

  end

  @. ∂uᵢh∂t*= 0;

  for (bᵢ,uᵢ,kᵢ) ∈ zip([vars.bx,vars.by,vars.bz],[vars.ux,vars.uy,vars.uz],[grid.kr,grid.l,grid.m])
        for (bⱼ,uⱼ,kⱼ,j) ∈ zip([vars.bx,vars.by,vars.bz],[vars.ux,vars.uy,vars.uz],[grid.kr,grid.l,grid.m],[1, 2, 3])

          #We split the computation into two parts to save memory
          @. vars.nonlin1  *= 0;
          @. vars.nonlinh1 *= 0;
          uᵢuⱼ  = vars.nonlin1;
          uᵢuⱼh = vars.nonlinh1;
          @. uᵢuⱼ = uᵢ*uⱼ;
          mul!(uᵢuⱼh, grid.rfftplan, uᵢuⱼ);

          # Perform the actual calculation
          @. ∂uᵢh∂t += -im*kᵢ*(δ(a,j)-kₐ*kⱼ*k⁻²)*uᵢuⱼh;

          @. vars.nonlin1  *= 0;
          @. vars.nonlinh1 *= 0;
          bᵢbⱼ  = vars.nonlin1;
          bᵢbⱼh = vars.nonlinh1;
          @. bᵢbⱼ = bᵢ*bⱼ;;
          mul!(bᵢbⱼh, grid.rfftplan, bᵢbⱼ);

          # Perform the actual calculation
          @. ∂uᵢh∂t += im*kᵢ*(δ(a,j)-kₐ*kⱼ*k⁻²)*bᵢbⱼh;

        end
    end

    #Compute the diffusion term  - νk^2 u_i
    @. ∂uᵢh∂t += -grid.Krsq*params.ν*uᵢh;
    return nothing

end

# B function
function BᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")

	#To Update B_i, we have two terms to compute:
	# ∂B_i/∂t = im ∑_j k_j*(b_iu_j - u_ib_j)  - η k^2 B_i
	#We split it into two part for sparating the computation.

  # declare the var u_i, b_i for computation
	if direction == "x"

		uᵢ  = vars.ux;
		bᵢ  = vars.bx;
		bᵢh = vars.bxh;
		∂Bᵢh∂t = @view N[:,:,:,params.bx_ind];

	elseif direction == "y"

		uᵢ  = vars.uy;
		bᵢ  = vars.by;
		bᵢh = vars.byh;
		∂Bᵢh∂t = @view N[:,:,:,params.by_ind];

	elseif direction == "z"

		uᵢ  = vars.uz;
		bᵢ  = vars.bz;
		bᵢh = vars.bzh;
		∂Bᵢh∂t = @view N[:,:,:,params.bz_ind];

	else

		@warn "Warning : Unknown direction is declerad"

	end

    @. ∂Bᵢh∂t*= 0;

    #Compute the first term, im ∑_j k_j*(b_iu_j - u_ib_j)
    for (bⱼ,uⱼ,kⱼ) ∈ zip([vars.bx,vars.by,vars.bz],[vars.ux,vars.uy,vars.uz],[grid.kr,grid.l,grid.m])

        # Initialization
        @. vars.nonlin1  *= 0;
        @. vars.nonlinh1 *= 0;
        uᵢbⱼ  = vars.nonlin1;
        uᵢbⱼh = vars.nonlinh1;
        # Perform Computation in Real space
        @. uᵢbⱼ = uᵢ*bⱼ;
        mul!(uᵢbⱼh, grid.rfftplan, uᵢbⱼ);
        # Actual Advection
        @. ∂Bᵢh∂t += im*kⱼ*uᵢbⱼh;

        @. vars.nonlin1  *= 0;
        @. vars.nonlinh1 *= 0;
        bᵢuⱼ  = vars.nonlin1;
        bᵢuⱼh = vars.nonlinh1;
        @. bᵢuⱼ = bᵢ*uⱼ;
        mul!(bᵢuⱼh, grid.rfftplan, bᵢuⱼ);
        @. ∂Bᵢh∂t += -im*kⱼ*bᵢuⱼh;

    end

    #Compute the diffusion term  - ηk^2 B_i
    @. ∂Bᵢh∂t += -grid.Krsq*params.η*bᵢh;

    return nothing

end

function MHDcalcN_advection!(N, sol, t, clock, vars, params, grid)

  #Update V + B Fourier Conponment
  @. vars.uxh = sol[:, :, :, params.ux_ind];
  @. vars.uyh = sol[:, :, :, params.uy_ind];
  @. vars.uzh = sol[:, :, :, params.uz_ind];

  @. vars.bxh = sol[:, :, :, params.bx_ind];
  @. vars.byh = sol[:, :, :, params.by_ind];
  @. vars.bzh = sol[:, :, :, params.bz_ind];

  #Update V + B Real Conponment
  ldiv!(vars.ux, grid.rfftplan, deepcopy(vars.uxh))
  ldiv!(vars.uy, grid.rfftplan, deepcopy(vars.uyh))
  ldiv!(vars.uz, grid.rfftplan, deepcopy(vars.uzh))
  ldiv!(vars.bx, grid.rfftplan, deepcopy(vars.bxh))
  ldiv!(vars.by, grid.rfftplan, deepcopy(vars.byh))
  ldiv!(vars.bz, grid.rfftplan, deepcopy(vars.bzh))

  #Update V Advection
  UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")
  UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="y")
  UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="z")

  #Update B Advection
  BᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")
  BᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="y")
  BᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="z")

  return nothing
end

function MHDupdatevars!(prob)
  vars, grid, sol, params = prob.vars, prob.grid, prob.sol, prob.params

  dealias!(sol, grid)

  #Update V + B Fourier Conponment
  @. vars.uxh = sol[:, :, :, params.ux_ind];
  @. vars.uyh = sol[:, :, :, params.uy_ind];
  @. vars.uzh = sol[:, :, :, params.uz_ind];

  @. vars.bxh = sol[:, :, :, params.bx_ind];
  @. vars.byh = sol[:, :, :, params.by_ind];
  @. vars.bzh = sol[:, :, :, params.bz_ind];

  #Update V + B Real Conponment
  ldiv!(vars.ux, grid.rfftplan, deepcopy(vars.uxh)) # deepcopy() since inverse real-fft destroys its input
  ldiv!(vars.uy, grid.rfftplan, deepcopy(vars.uyh)) # deepcopy() since inverse real-fft destroys its input
  ldiv!(vars.uz, grid.rfftplan, deepcopy(vars.uzh)) # deepcopy() since inverse real-fft destroys its input
  ldiv!(vars.bx, grid.rfftplan, deepcopy(vars.bxh)) # deepcopy() since inverse real-fft destroys its input
  ldiv!(vars.by, grid.rfftplan, deepcopy(vars.byh)) # deepcopy() since inverse real-fft destroys its input
  ldiv!(vars.bz, grid.rfftplan, deepcopy(vars.bzh)) # deepcopy() since inverse real-fft destroys its input

  return nothing
end

end

max_stars_repo_name: MHDFlows/MHDFlows

max_issues_repo_name: MHDFlows/MHDFlows

max_forks_repo_name: MHDFlows/MHDFlows


text:
using TensorOperations

include("RCCSDHelper.jl")

function RCCSD(moints::IntegralHelper{T,E,O}, newT1::AbstractArray{T,2}, newT2::AbstractArray{T,4},
                    alg::RCCSDa) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    # Print intro
    cc_header()
    Eref = moints.orbitals.sd_energy
    ndocc = moints.molecule.Nα
    nvir = size(moints.orbitals.C,1) - ndocc
    core = Options.get("drop_occ")
    inac = Options.get("drop_vir")

    if core ≥ ndocc
        throw(FermiException("invalid number of frozen orbitals ($core) for $ndocc doubly occupied orbitals"))
    end
    if inac ≥ nvir
        throw(FermiException("invalid number of inactive virtual orbitals ($inac) for $nvir total virtual orbitals"))
    end

    # Get relevant Options
    cc_max_iter = Options.get("cc_max_iter")
    cc_e_conv = Options.get("cc_e_conv")
    cc_max_rms = Options.get("cc_max_rms")
    precision_override = Options.get("precision_override")

    dp = Options.get("cc_damp_ratio")
    do_diis = Options.get("cc_diis")
    ndiis = Options.get("cc_ndiis")
    diis_start = Options.get("diis_start")
    diis_relax = Options.get("cc_diis_relax")

    diis_prec = Options.get("diis_prec")
    diis_prec = diis_prec == "double" ? Float64 : diis_prec == "single" ? Float32 : Float16
    do_diis ? DM_T1 = Fermi.DIIS.DIISManager{T, diis_prec}(size=ndiis) : nothing
    do_diis ? DM_T2 = Fermi.DIIS.DIISManager{T, diis_prec}(size=ndiis) : nothing

    output("\tDropped Occupied Orbitals →  {:3.0f}", core)
    output("\tDropped Virtual Orbitals  →  {:3.0f}", inac)
    output("\n"*repeat("-",80))
    output("\tOptions:")
    output("\tPrecision? {}", T)
    output("\tDIIS? {}", do_diis)
    output("\tDIIS Vectors? {}", (do_diis ? ndiis : 0))
    output("\tDIIS Precision? {}", diis_prec)
    output("\tDamping percentage ? {}", dp)
    output("\t\tcc_max_iter →  {:3.0d}", cc_max_iter)
    output("\t\tcc_e_conv   →  {:2.0e}", cc_e_conv)
    output("\t\tcc_max_rms  →  {:2.0e}", cc_max_rms)
    if (cc_e_conv < eps(T) || cc_max_rms < eps(T)) && !(precision_override)
        output("\t⚠ WARNING ⚠  cc_e_conv set to less than ϵ ({}) for precision {}", eps(T), T)
        output("\t             CCSD is unlikely to converge to your standards.")
        output("\t             OVERRIDING set convergence criterion.")
        output("\t             Use @set `precision_override` true to perform the computation as entered.")
        cc_e_conv = eps(T)
        cc_max_rms = eps(T)
    end

    # Initialize Loop parameters
    r1 = 1
    r2 = 1
    dE = 1
    rms = 1
    ite = 1
    T1 = deepcopy(newT1)
    T2 = deepcopy(newT2)

    output(repeat("-", 80))

    # Compute Guess Energy
    Ecc = cc_update_energy(newT1, newT2, moints, alg)
    Eguess = Ecc + Eref

    output("\tGuess Correlation Energy:   {:15.10f}", Ecc)
    output("\tGuess Total Energy:         {:15.10f}\n", Eguess)
    output("    Starting CC Iterations\n")

    dE = 1
    rms = 1
    ite = 1

    if do_diis
        DM_T1 = Fermi.DIIS.DIISManager{T, diis_prec}(size=ndiis)
        DM_T2 = Fermi.DIIS.DIISManager{T, diis_prec}(size=ndiis)
    end

    main_time = 0
    output("{:10s}    {: 15s}    {: 12s}    {:12s}    {:10s}    {:s5}","Iteration","CC Energy","ΔE","Max RMS","time (s)", "diis")
    while (abs(dE) > cc_e_conv || rms > cc_max_rms)
        if ite > cc_max_iter
            output("\n⚠  CC Equations did not converge in {:1.0d} iterations.", cc_max_iter)
            break
        end
        t = @elapsed begin

            T1 .= newT1
            T2 .= newT2
            oldE = Ecc

            update_amp!(newT1, newT2, T1, T2, moints, alg)

            # Compute residues
            r1 = sqrt(sum((newT1 .- T1).^2)/length(T1))
            r2 = sqrt(sum((newT2 .- T2).^2)/length(T2))

            diis = false
            if do_diis
                # Update DIIS vector
                e1 = (newT1 - T1)
                e2 = (newT2 - T2)
                push!(DM_T1,newT1,e1)
                push!(DM_T2,newT2,e2)
                # Using DIIS at every iteration is BAD! So do it every three steps
                if ite ≥ diis_start && ite % diis_relax == 0
                    # Get new extrapolated amplitudes from DIIS
                    newT2 = Fermi.DIIS.extrapolate(DM_T2;add_res=true)
                    newT1 = Fermi.DIIS.extrapolate(DM_T1;add_res=true)
                    diis = true
                end
            end

            # Apply damping
            if dp > 0
                newT1 .= (1-dp)*newT1 .+ dp*T1
                newT2 .= (1-dp)*newT2 .+ dp*T2
            end

            Ecc = cc_update_energy(newT1, newT2, moints, alg)
        end
        rms = max(r1,r2)
        dE = Ecc - oldE
        main_time += t
        output("    {:<5.0d}    {:<15.10f}    {:>+12.10f}    {:<12.10f}    {:<10.5f}    {:s5}", ite, Ecc, dE, rms, t, diis)
        ite += 1
    end
    output("\nMain CCSD iterations done in {:5.5f} s", main_time)
    output("Average time per iteration {:5.5f}", main_time/(ite-1))

    # Converged?
    conv = false
    if abs(dE) < cc_e_conv && rms < cc_max_rms
        output("\n 🍾 Equations Converged!")
        conv = true
    end

    output(" @Final CCSD Correlation Energy:     {:15.10f}", Ecc)
    output(" @Final CCSD Energy:                 {:15.10f}", Ecc+Eref)
    output(repeat("-",80))

    return RCCSD(Eguess, Ecc+Eref, Ecc, newT1, newT2, dE, rms)
end


text:
export Editor, editor, setbars, barson, barsoff, setcursor

type Editor
  file
  w::Window
end

filetitle(e) = e.file == nothing ? "Julia" : basename(e.file)

Blink.msg(e::Editor, args...) = Blink.msg(e.w, args...)
Blink.handlers(e::Editor) = Blink.handlers(e.w)

# I know, I know, this code is a mess.

function loadeditor(p::Page; value = "", ver = "5.0.0")
  for f in ["codemirror.min.js"
            "codemirror.min.css"
            "addon/display/rulers.min.js"
            "addon/selection/active-line.min.js"
            "addon/search/searchcursor.min.js"
            "keymap/sublime.min.js"]
    load!(p, "https://cdnjs.cloudflare.com/ajax/libs/codemirror/$ver/$f")
  end

  for f in ["julia.js", "editor.css", "june.css", "bars.js", "bars.css"]
    load!(p, Pkg.dir("DevTools", "res", f))
  end

  @js_ p CodeMirror.keyMap.blink = d("fallthrough"=>"sublime")

  body!(p, "", fade = false)
  @js_ p cm = CodeMirror(document.body,
                         $(d(:value => value,
                             :mode=>"julia2",
                             :theme=>"june",
                             :keyMap=>"blink",
                             :lineNumbers=>true,
                             :styleActiveLine=>true,
                             :rulers=>[80],
                             :showCursorWhenSelecting=>true)))
  @js_ p Bars.hook(cm)
  handle_eval(p)
end

function Editor(value = ""; file = nothing)
  w = Window(d(:title=>file == nothing ? "Julia" : basename(file)))
  loadeditor(w.content, value = value)
  ed = Editor(file, w)
  handle_dirty(ed)
  handle_save(ed)
  return ed
end

editor(f) = Editor(read(f, String), file = f)

setbars(e, ls) = @js_ e Bars.set(cm, $ls)
barson(e) = @js_ e Bars.on(cm)
barsoff(e) = @js_ e Bars.off(cm)

function centrecursor(ed)
  @js_ ed begin
    @var l = cm.getCursor().line
    @var y = cm.charCoords(d("line"=>l, "ch"=>0), "local").top
    @var height = cm.getScrollerElement().offsetHeight
    cm.scrollTo(0, y - height/2 + 55)
  end
end

function setcursor(ed, line, ch = 0)
  @js_ ed cm.setCursor($(line-1), $ch)
  centrecursor(ed)
end

norm(key::AbstractString) = replace(key, r"\bc\b"i, Sys.is_apple() ? "Cmd" : "Ctrl")

function keymap(ed, key, res)
  @js_ ed begin
    @var map = CodeMirror.keyMap.blink
    map[$(norm(key))] = $res
#     CodeMirror.normalizeKeyMap(map)
  end
end

function handle_dirty(e::Editor)
  @js_ e cm.on("changes", () -> Blink.msg("change", d("clean"=>cm.isClean())))

  handle(e, "change") do data
    t = filetitle(e)
    data["clean"] || (t *= "*")
    title(e.w, t)
  end
end

function handle_save(ed::Editor)
  keymap(ed, "C-S", :(cm -> Blink.msg("save", d("code"=>cm.getValue()))))

  handle(ed, "save") do data
    if ed.file != nothing && isfile(ed.file)
      @js_ ed cm.markClean()
      title(ed.w, filetitle(ed))
      open(io -> write(io, data["code"]), ed.file, "w")
    end
  end
end

function handle_eval(ed)
  keymap(ed, "Shift-Enter", :(cm -> Blink.msg("eval", d("code"=>cm.getSelection()))))

  handle(ed, "eval") do data
    try
      result = eval(Main, parse(data["code"]))
      @js_ ed console.log($(sprint(show, result)))
    catch e
      @js_ ed console.log($(sprint(showerror, e, catch_backtrace())))
    end
  end
end


text:
@cache struct SSPRK22Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end

struct SSPRK22ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::SSPRK22,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  SSPRK22Cache(u,uprev,k,fsalfirst,alg.stage_limiter!,alg.step_limiter!)
end

alg_cache(alg::SSPRK22,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = SSPRK22ConstantCache()

@cache struct SSPRK33Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end

struct SSPRK33ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::SSPRK33,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  SSPRK33Cache(u,uprev,k,fsalfirst,alg.stage_limiter!,alg.step_limiter!)
end

alg_cache(alg::SSPRK33,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = SSPRK33ConstantCache()

@cache struct KYKSSPRK42Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
   u::uType
   uprev::uType
   k::rateType
   tmp::uType
   fsalfirst::rateType
   tab::TabType
 end

struct KYKSSPRK42ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache
  α20::T
  α21::T
  α30::T
  α32::T
  α40::T
  α43::T
  β10::T
  β21::T
  β30::T
  β32::T
  β40::T
  β43::T
  c1::T2
  c2::T2
  c3::T2
end

function KYKSSPRK42ConstantCache(T, T2)
  α20 = T(0.394806441339829)
  α21 = T(0.605193558660171)
  α30 = T(0.002797307087390)
  α32 = T(0.997202692912610)
  α40 = T(0.252860909354373)
  α43 = T(0.747139090645627)
  β10 = T(0.406584463657504)
  β21 = T(0.246062298456822)
  β30 = T(0.013637216641451)
  β32 = T(0.405447122055692)
  β40 = T(0.016453567333598)
  β43 = T(0.303775146447707)
  c1 = T2(0.406584463657504)
  c2 = T2(0.4921245969136438)
  c3 = T2(0.9098323119879613)
  KYKSSPRK42ConstantCache(α20, α21, α30, α32, α40, α43, β10, β21, β30, β32, β40, β43, c1, c2, c3)
end

function alg_cache(alg::KYKSSPRK42,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = KYKSSPRK42ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  KYKSSPRK42Cache(u,uprev,k,tmp, fsalfirst,tab)
end

function alg_cache(alg::KYKSSPRK42,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  KYKSSPRK42ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
end

@cache struct SSPRK53Cache{uType,rateType,StageLimiter,StepLimiter,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  tmp::uType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  tab::TabType
end

struct SSPRK53ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α30::T
  α32::T
  α40::T
  α43::T
  α52::T
  α54::T
  β10::T
  β21::T
  β32::T
  β43::T
  β54::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2

  function SSPRK53ConstantCache(T, T2)
    α30 = T(0.355909775063327)
    α32 = T(0.644090224936674)
    α40 = T(0.367933791638137)
    α43 = T(0.632066208361863)
    α52 = T(0.237593836598569)
    α54 = T(0.762406163401431)
    β10 = T(0.377268915331368)
    β21 = T(0.377268915331368)
    β32 = T(0.242995220537396)
    β43 = T(0.238458932846290)
    β54 = T(0.287632146308408)
    c1 = T2(0.377268915331368)
    c2 = T2(0.754537830662736)
    c3 = T2(0.728985661612188)
    c4 = T2(0.699226135931670)

    new{T,T2}(α30, α32, α40, α43, α52, α54, β10, β21, β32, β43, β54, c1, c2, c3, c4)
  end
end

function alg_cache(alg::SSPRK53,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u)
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  tab = SSPRK53ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SSPRK53Cache(u,uprev,k,fsalfirst,tmp,alg.stage_limiter!,alg.step_limiter!,tab)
end

function alg_cache(alg::SSPRK53,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SSPRK53ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
end

@cache struct SHLDDRK52Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  tab::TabType
end

struct SHLDDRK52ConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache
  α2::T1
  α3::T1
  α4::T1
  α5::T1
  β1::T1
  β2::T1
  β3::T1
  β4::T1
  β5::T1
  c2::T2
  c3::T2
  c4::T2
  c5::T2
end

function SHLDDRK52ConstantCache(T1,T2)
  α2 = T1(-0.6913065)
  α3 = T1(-2.655155)
  α4 = T1(-0.8147688)
  α5 = T1(-0.6686587)
  β1 = T1(0.1)
  β2 = T1(0.75)
  β3 = T1(0.7)
  β4 = T1(0.479313)
  β5 = T1(0.310392)
  c2 = T2(0.1)
  c3 = T2(0.3315201)
  c4 = T2(0.4577796)
  c5 = T2(0.8666528)
  SHLDDRK52ConstantCache(α2, α3, α4, α5, β1, β2, β3, β4, β5, c2, c3, c4, c5)
end

function alg_cache(alg::SHLDDRK52,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SHLDDRK52ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
end

function alg_cache(alg::SHLDDRK52,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = SHLDDRK52ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SHLDDRK52Cache(u,uprev,k,tmp,fsalfirst,tab)
end

@cache mutable struct SHLDDRK_2NCache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  tab::TabType
  step::Int
end

mutable struct SHLDDRK_2NConstantCache{T1,T2} <: OrdinaryDiffEqConstantCache
  α21::T1
  α31::T1
  α41::T1
  α51::T1
  β11::T1
  β21::T1
  β31::T1
  β41::T1
  β51::T1
  c21::T2
  c31::T2
  c41::T2
  c51::T2

  α22::T1
  α32::T1
  α42::T1
  α52::T1
  α62::T1
  β12::T1
  β22::T1
  β32::T1
  β42::T1
  β52::T1
  β62::T1
  c22::T2
  c32::T2
  c42::T2
  c52::T2
  c62::T2

  step::Int
end

function SHLDDRK_2NConstantCache(T1,T2)
  α21 = T1(-0.6051226)
  α31 = T1(-2.0437564)
  α41 = T1(-0.7406999)
  α51 = T1(-4.4231765)
  β11 = T1(0.2687454)
  β21 = T1(0.8014706)
  β31 = T1(0.5051570)
  β41 = T1(0.5623568)
  β51 = T1(0.0590065)
  c21 = T2(0.2687454)
  c31 = T2(0.5852280)
  c41 = T2(0.6827066)
  c51 = T2(1.1646854)

  α22 = T1(-0.4412737)
  α32 = T1(-1.0739820)
  α42 = T1(-1.7063570)
  α52 = T1(-2.7979293)
  α62 = T1(-4.0913537)
  β12 = T1(0.1158488)
  β22 = T1(0.3728769)
  β32 = T1(0.7379536)
  β42 = T1(0.5798110)
  β52 = T1(1.0312849)
  β62 = T1(0.15)
  c22 = T2(0.1158485)
  c32 = T2(0.3241850)
  c42 = T2(0.6193208)
  c52 = T2(0.8034472)
  c62 = T2(0.9184166)
  SHLDDRK_2NConstantCache(α21, α31, α41, α51, β11, β21, β31, β41, β51, c21, c31, c41, c51, α22, α32, α42, α52, α62, β12, β22, β32, β42, β52, β62, c22, c32, c42, c52, c62, 1)
end

function alg_cache(alg::SHLDDRK_2N,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SHLDDRK_2NConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
end

function alg_cache(alg::SHLDDRK_2N,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = SHLDDRK_2NConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SHLDDRK_2NCache(u,uprev,k,tmp,fsalfirst,tab,1)
end

@cache struct SSPRK53_2N1Cache{uType,rateType,StageLimiter,StepLimiter,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  tab::TabType
end

struct SSPRK53_2N1ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α40::T
  α43::T
  β10::T
  β21::T
  β32::T
  β43::T
  β54::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2

  function SSPRK53_2N1ConstantCache(T, T2)
    α40 = T(0.571403511494104)
    α43 = T(0.428596488505896)
    β10 = T(0.443568244942995)
    β21 = T(0.291111420073766)
    β32 = T(0.270612601278217)
    β43 = T(0.110577759392786)
    β54 = T(0.458557505351052)
    c1 = T2(0.443568244942995)
    c2 = T2(0.734679665016762)
    c3 = T2(1.005292266294979)
    c4 = T2(0.541442494648948)

    new{T,T2}(α40, α43,β10, β21, β32, β43, β54, c1, c2, c3, c4)
  end
end

function alg_cache(alg::SSPRK53_2N1,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  tab = SSPRK53_2N1ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SSPRK53_2N1Cache(u,uprev,k,fsalfirst,alg.stage_limiter!,alg.step_limiter!,tab)
end

function alg_cache(alg::SSPRK53_2N1,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SSPRK53_2N1ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
end

@cache struct SSPRK53_2N2Cache{uType,rateType,StageLimiter,StepLimiter,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  tab::TabType
end

struct SSPRK53_2N2ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α30::T
  α32::T
  α50::T
  α54::T
  β10::T
  β21::T
  β32::T
  β43::T
  β54::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2

  function SSPRK53_2N2ConstantCache(T, T2)
    α30 = T(0.682342861037239)
    α32 = T(0.317657138962761)
    α50 = T(0.045230974482400)
    α54 = T(0.954769025517600)
    β10 = T(0.465388589249323)
    β21 = T(0.465388589249323)
    β32 = T(0.124745797313998)
    β43 = T(0.465388589249323)
    β54 = T(0.154263303748666)
    c1 = T2(0.465388589249323)
    c2 = T2(0.930777178498646)
    c3 = T2(0.420413812847710)
    c4 = T2(0.885802402097033)

    new{T,T2}(α30, α32,α50,α54,β10, β21, β32, β43, β54, c1, c2, c3, c4)
  end
end

function alg_cache(alg::SSPRK53_2N2,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  tab = SSPRK53_2N2ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SSPRK53_2N2Cache(u,uprev,k,fsalfirst,alg.stage_limiter!,alg.step_limiter!,tab)
end

function alg_cache(alg::SSPRK53_2N2,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SSPRK53_2N2ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
end

@cache struct SSPRK53_HCache{uType,rateType,StageLimiter,StepLimiter,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  tmp::uType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  tab::TabType
end

struct SSPRK53_HConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α30::T
  α32::T
  α40::T
  α41::T
  α43::T
  β10::T
  β21::T
  β32::T
  β43::T
  β54::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2

  function SSPRK53_HConstantCache(T,T2)
    α30 = T(0.308684154602513)
    α32 = T(0.691315845397487)
    α40 = T(0.280514990468574)
    α41 = T(0.270513101776498)
    α43 = T(0.448971907754928)
    β10 = T(0.377268915331368)
    β21 = T(0.377268915331368)
    β32 = T(0.260811979144498)
    β43 = T(0.169383144652957)
    β54 = T(0.377268915331368)
    c1 = T2(0.377268915331368)
    c2 = T2(0.754537830662737)
    c3 = T2(0.782435937433493)
    c4 = T2(0.622731084668631)

    new{T,T2}(α30, α32, α40, α41, α43, β10, β21, β32, β43, β54, c1, c2, c3, c4)
  end
end

function alg_cache(alg::SSPRK53_H,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u)
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  tab = SSPRK53_HConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SSPRK53_HCache(u,uprev,k,fsalfirst,tmp,alg.stage_limiter!,alg.step_limiter!,tab)
end

function alg_cache(alg::SSPRK53_H,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SSPRK53_HConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
end

@cache struct SSPRK63Cache{uType,rateType,StageLimiter,StepLimiter,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  tmp::uType
  u₂::uType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  tab::TabType
end

struct SSPRK63ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α40::T
  α41::T
  α43::T
  α62::T
  α65::T
  β10::T
  β21::T
  β32::T
  β43::T
  β54::T
  β65::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2

  function SSPRK63ConstantCache(T, T2)
    α40 = T(0.476769811285196)
    α41 = T(0.098511733286064)
    α43 = T(0.424718455428740)
    α62 = T(0.155221702560091)
    α65 = T(0.844778297439909)
    β10 = T(0.284220721334261)
    β21 = T(0.284220721334261)
    β32 = T(0.284220721334261)
    β43 = T(0.120713785765930)
    β54 = T(0.284220721334261)
    β65 = T(0.240103497065900)
    c1 = T2(0.284220721334261)
    c2 = T2(0.568441442668522)
    c3 = T2(0.852662164002783)
    c4 = T2(0.510854218958172)
    c5 = T2(0.795074940292433)

    new{T,T2}(α40, α41, α43, α62, α65, β10, β21, β32, β43, β54, β65, c1, c2, c3, c4, c5)
  end
end

function alg_cache(alg::SSPRK63,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u)
  u₂ = zero(u)
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  tab = SSPRK63ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SSPRK63Cache(u,uprev,k,fsalfirst,tmp,u₂,alg.stage_limiter!,alg.step_limiter!,tab)
end

function alg_cache(alg::SSPRK63,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SSPRK63ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
end

@cache struct SSPRK73Cache{uType,rateType,StageLimiter,StepLimiter,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  tmp::uType
  u₁::uType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  tab::TabType
end

struct SSPRK73ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α40::T
  α43::T
  α50::T
  α51::T
  α54::T
  α73::T
  α76::T
  β10::T
  β21::T
  β32::T
  β43::T
  β54::T
  β65::T
  β76::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2

  function SSPRK73ConstantCache(T, T2)
    α40 = T(0.184962588071072)
    α43 = T(0.815037411928928)
    α50 = T(0.180718656570380)
    α51 = T(0.314831034403793)
    α54 = T(0.504450309025826)
    α73 = T(0.120199000000000)
    α76 = T(0.879801000000000)
    β10 = T(0.233213863663009)
    β21 = T(0.233213863663009)
    β32 = T(0.233213863663009)
    β43 = T(0.190078023865845)
    β54 = T(0.117644805593912)
    β65 = T(0.233213863663009)
    β76 = T(0.205181790464579)
    c1 = T2(0.233213863663009)
    c2 = T2(0.466427727326018)
    c3 = T2(0.699641590989027)
    c4 = T2(0.760312095463379)
    c5 = T2(0.574607439040817)
    c6 = T2(0.807821302703826)

    new{T,T2}(α40, α43, α50, α51, α54, α73, α76, β10, β21, β32, β43, β54, β65, β76, c1, c2, c3, c4, c5, c6)
  end
end

function alg_cache(alg::SSPRK73,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u)
  u₁ = zero(u)
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  tab = SSPRK73ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SSPRK73Cache(u,uprev,k,fsalfirst,tmp,u₁,alg.stage_limiter!,alg.step_limiter!,tab)
end

function alg_cache(alg::SSPRK73,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SSPRK73ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
end

@cache struct SSPRK83Cache{uType,rateType,StageLimiter,StepLimiter,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  tmp::uType
  u₂::uType
  u₃::uType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  tab::TabType
end

struct SSPRK83ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α50::T
  α51::T
  α54::T
  α61::T
  α65::T
  α72::T
  α73::T
  α76::T
  β10::T
  β21::T
  β32::T
  β43::T
  β54::T
  β65::T
  β76::T
  β87::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  c7::T2

  function SSPRK83ConstantCache(T, T2)
    α50 = T(0.421366967085359)
    α51 = T(0.005949401107575)
    α54 = T(0.572683631807067)
    α61 = T(0.004254010666365)
    α65 = T(0.995745989333635)
    α72 = T(0.104380143093325)
    α73 = T(0.243265240906726)
    α76 = T(0.652354615999950)
    β10 = T(0.195804015330143)
    β21 = T(0.195804015330143)
    β32 = T(0.195804015330143)
    β43 = T(0.195804015330143)
    β54 = T(0.112133754621673)
    β65 = T(0.194971062960412)
    β76 = T(0.127733653231944)
    β87 = T(0.195804015330143)
    c1 = T2(0.195804015330143)
    c2 = T2(0.391608030660286)
    c3 = T2(0.587412045990429)
    c4 = T2(0.783216061320572)
    c5 = T2(0.561833689734037)
    c6 = T2(0.755247658555329)
    c7 = T2(0.804195984669857)

    new{T,T2}(α50, α51, α54, α61, α65, α72, α73, α76, β10, β21, β32, β43, β54, β65, β76, β87, c1, c2, c3, c4, c5, c6, c7)
  end
end

function alg_cache(alg::SSPRK83,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u)
  u₂ = zero(u)
  u₃ = zero(u)
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  tab = SSPRK83ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SSPRK83Cache(u,uprev,k,fsalfirst,tmp,u₂,u₃,alg.stage_limiter!,alg.step_limiter!,tab)
end

function alg_cache(alg::SSPRK83,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SSPRK83ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
end

@cache struct SSPRK43Cache{uType,rateType,uNoUnitsType,StageLimiter,StepLimiter,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  utilde::uType
  atmp::uNoUnitsType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  tab::TabType
end

struct SSPRK43ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  one_third_u::T
  two_thirds_u::T
  half_u::T
  half_t::T2

  function SSPRK43ConstantCache(T, T2)
    one_third_u = inv(T(3))
    two_thirds_u = 2 * one_third_u
    half_u = T(0.5)
    half_t = T2(0.5)

    new{T,T2}(one_third_u, two_thirds_u, half_u, half_t)
  end
end

function alg_cache(alg::SSPRK43,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  utilde = zero(u)
  atmp = similar(u,uEltypeNoUnits)
  tab = SSPRK43ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SSPRK43Cache(u,uprev,k,fsalfirst,utilde,atmp,alg.stage_limiter!,alg.step_limiter!,tab)
end

function alg_cache(alg::SSPRK43,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SSPRK43ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
end

@cache struct SSPRK432Cache{uType,rateType,uNoUnitsType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  utilde::uType
  atmp::uNoUnitsType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end

struct SSPRK432ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::SSPRK432,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  utilde = zero(u)
  atmp = similar(u,uEltypeNoUnits)
  SSPRK432Cache(u,uprev,k,fsalfirst,utilde,atmp,alg.stage_limiter!,alg.step_limiter!)
end

alg_cache(alg::SSPRK432,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = SSPRK432ConstantCache()

@cache mutable struct SSPRKMSVS32Cache{uType,rateType,dtArrayType,dtType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  u_2::uType
  u_1::uType
  k::rateType
  tmp::uType
  dts::dtArrayType
  dtf::dtArrayType
  μ::dtType
  v_n::Float64
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  step::Int
end

@cache mutable struct SSPRKMSVS32ConstantCache{uType,dtArrayType,dtType} <: OrdinaryDiffEqConstantCache
  u_2::uType
  u_1::uType
  dts::dtArrayType
  dtf::dtArrayType
  μ::dtType
  v_n::Float64
  step::Int
end

function alg_cache(alg::SSPRKMSVS32,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  fsalfirst = zero(rate_prototype)
  dts = fill(zero(dt),3)
  dtf = fill(zero(dt),2)
  μ = zero(dt)
  u_2 = zero(u)
  u_1 = zero(u)
  k  = zero(rate_prototype)
  tmp = zero(u)
  SSPRKMSVS32Cache(u,uprev,fsalfirst,u_2,u_1,k,tmp,dts,dtf,μ,0.5,alg.stage_limiter!,alg.step_limiter!,1)
end

function alg_cache(alg::SSPRKMSVS32,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  dts = fill(zero(dt),3)
  dtf = fill(zero(dt),2)
  μ = zero(dt)
  u_2 = u
  u_1 = u
  SSPRKMSVS32ConstantCache(u_2,u_1,dts,dtf,μ,0.5,1)
end

@cache mutable struct SSPRKMSVS43Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  u_3::uType
  u_2::uType
  u_1::uType
  k::rateType
  k1::rateType
  k2::rateType
  k3::rateType
  tmp::uType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  step::Int
end

@cache mutable struct SSPRKMSVS43ConstantCache{uType,rateType} <: OrdinaryDiffEqConstantCache
  u_3::uType
  u_2::uType
  u_1::uType
  k1::rateType
  k2::rateType
  k3::rateType
  step::Int
end

function alg_cache(alg::SSPRKMSVS43,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  fsalfirst = zero(rate_prototype)
  u_3 = zero(u)
  u_2 = zero(u)
  u_1 = zero(u)
  k  = zero(rate_prototype)
  k1  = zero(rate_prototype)
  k2  = zero(rate_prototype)
  k3  = zero(rate_prototype)
  tmp = zero(u)
  SSPRKMSVS43Cache(u,uprev,fsalfirst,u_3,u_2,u_1,k,k1,k2,k3,tmp,alg.stage_limiter!,alg.step_limiter!,1)
end

function alg_cache(alg::SSPRKMSVS43,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  u_3 = u
  u_2 = u
  u_1 = u
  k1  = rate_prototype
  k2  = rate_prototype
  k3  = rate_prototype
  SSPRKMSVS43ConstantCache(u_3,u_2,u_1,k1,k2,k3,1)
end

@cache struct SSPRK932Cache{uType,rateType,uNoUnitsType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  utilde::uType
  atmp::uNoUnitsType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end

struct SSPRK932ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::SSPRK932,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  utilde = zero(u)
  atmp = similar(u,uEltypeNoUnits)
  SSPRK932Cache(u,uprev,k,fsalfirst,utilde,atmp,alg.stage_limiter!,alg.step_limiter!)
end

alg_cache(alg::SSPRK932,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = SSPRK932ConstantCache()

@cache struct SSPRK54Cache{uType,rateType,StageLimiter,StepLimiter,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  k₃::rateType
  u₂::uType
  u₃::uType
  tmp::uType # should be u₄, but tmp is needed for callbacks
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  tab::TabType
end

struct SSPRK54ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  β10::T
  α20::T
  α21::T
  β21::T
  α30::T
  α32::T
  β32::T
  α40::T
  α43::T
  β43::T
  α52::T
  α53::T
  β53::T
  α54::T
  β54::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2

  function SSPRK54ConstantCache(T, T2)
    β10 = T(0.391752226571890)
    α20 = T(0.444370493651235)
    α21 = T(0.555629506348765)
    β21 = T(0.368410593050371)
    α30 = T(0.620101851488403)
    α32 = T(0.379898148511597)
    β32 = T(0.251891774271694)
    α40 = T(0.178079954393132)
    α43 = T(0.821920045606868)
    β43 = T(0.544974750228521)
    α52 = T(0.517231671970585)
    α53 = T(0.096059710526147)
    β53 = T(0.063692468666290)
    α54 = T(0.386708617503269)
    β54 = T(0.226007483236906)
    c1 = T2(0.391752226571890)
    c2 = T2(0.586079689311540)
    c3 = T2(0.474542363121400)
    c4 = T2(0.935010630967653)

    new{T,T2}(β10, α20, α21, β21, α30, α32, β32, α40, α43, β43, α52, α53, β53, α54, β54, c1, c2, c3, c4)
  end
end

function alg_cache(alg::SSPRK54,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  u₂ = zero(u)
  u₃ = zero(u)
  tmp = zero(u)
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  k₃ = zero(rate_prototype)
  tab = SSPRK54ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
  SSPRK54Cache(u,uprev,k,fsalfirst,k₃,u₂,u₃,tmp,alg.stage_limiter!,alg.step_limiter!,tab)
end

function alg_cache(alg::SSPRK54,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  SSPRK54ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
end

@cache struct SSPRK104Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  k₄::rateType
  tmp::uType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end

struct SSPRK104ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::SSPRK104,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  tmp = zero(u)
  k = zero(rate_prototype)
  if calck
    fsalfirst = similar(k)
  else
    fsalfirst = k
  end
  k₄ = zero(rate_prototype)
  SSPRK104Cache(u,uprev,k,fsalfirst,k₄,tmp,alg.stage_limiter!,alg.step_limiter!)
end

alg_cache(alg::SSPRK104,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = SSPRK104ConstantCache()