text:
using HDF5
using SparseArrays
using Random

file = h5open(tempname(), "w")
Random.seed!(42)
spmatrix = sprand(20, 20, 0.3)
matrix = rand(20, 20)

for (mat, T, name, issparse) in (
    (spmatrix, Muon.SparseDataset, "sparse", true),
    (matrix, Muon.TransposedDataset, "transposed", false),
)
    mat_t = mat'

    Muon.write_attr(file, "mat", mat)
    Muon.write_attr(file, "mat_t", mat_t)

    @testset "$name matrix read" begin
        @test Muon.read_matrix(file["mat"]) == mat
        @test Muon.read_matrix(file["mat_t"]) == mat_t
    end

    @testset "$name dataset" begin
        dset = T(file["mat"])
        dset_t = T(file["mat_t"])

        @testset "full read" begin
            @test read(dset) == mat
            @test read(dset_t) == mat_t
        end

        # double nesting to have more compact output in case someting fails
        @testset "single element access" begin
            @testset "index: $i,$j" for i in 1:size(mat, 1), j in 1:size(mat, 2)
                @test dset[i, j] == mat[i, j]
                @test dset_t[i, j] == mat_t[i, j]
            end
        end
        @testset "row access" begin
            @testset "row: $i" for i in 1:size(mat, 1)
                @test dset[i, :] == mat[i, :]
                @test dset_t[i, :] == mat_t[i, :]
            end
        end
        @testset "column access" begin
            @testset "column: $j" for j in 1:size(mat, 2)
                @test dset[:, j] == mat[:, j]
                @test dset_t[:, j] == mat_t[:, j]
            end
        end

        @testset "range access" begin
            @testset "range: $rowstart:$(rowstart + 4), $colstart:$(colstart + 4)" for (
                rowstart,
                colstart,
            ) in zip(
                1:(size(mat, 1) - 4),
                1:(size(mat, 2) - 4),
            )
                @test dset[rowstart:(rowstart + 4), colstart:(colstart + 4)] ==
                      mat[rowstart:(rowstart + 4), colstart:(colstart + 4)]
                @test dset_t[rowstart:(rowstart + 4), colstart:(colstart + 4)] ==
                      mat_t[rowstart:(rowstart + 4), colstart:(colstart + 4)]
            end
        end

        @testset "row range access" begin
            @testset "row: $i, column range: $colstart:$(colstart + 4)" for i in 1:size(mat, 1),
                colstart in 1:(size(mat, 2) - 4)

                @test dset[i, colstart:(colstart + 4)] == mat[i, colstart:(colstart + 4)]
                @test dset_t[i, colstart:(colstart + 4)] == mat_t[i, colstart:(colstart + 4)]
            end
        end
        @testset "column range access" begin
            @testset "column: $j, row range: $rowstart:$(rowstart + 4)" for j in 1:size(mat, 2),
                rowstart in 1:(size(mat, 1) - 4)

                @test dset[rowstart:(rowstart + 4), j] == mat[rowstart:(rowstart + 4), j]
                @test dset_t[rowstart:(rowstart + 4), j] == mat_t[rowstart:(rowstart + 4), j]
            end
        end

        @testset "single element modification" begin
            if issparse
                I, J, V = findnz(mat)
                ci = 1
                for j in 1:size(mat, 2)
                    while j > J[ci] && ci < length(I)
                        ci += 1
                    end
                    for i in 1:size(mat, 1)
                        while i > I[ci] && j >= J[ci] && ci < length(I)
                            ci += 1
                        end

                        @testset "index: $i, $j" begin
                            if i == I[ci] && j == J[ci]
                                dset[i, j] = mat[i, j] = dset_t[j, i] = mat_t[j, i] = π
                                @test dset == mat
                                @test dset_t == mat_t
                            else
                                @test_throws KeyError dset[i, j] = π
                                @test_throws KeyError dset_t[j, i] = π
                            end
                        end
                    end
                end
            else
                @testset "index: $i,$j" for i in 1:size(mat, 1), j in 1:size(mat, 2)
                    dset[i, j] = mat[i, j] = dset_t[j, i] = mat_t[j, i] = π
                    @test dset == mat
                    @test dset_t == mat_t
                end
            end
        end

        @testset "range modification" begin
            @testset "range: $rowstart:$(rowstart + 4), $colstart:$(colstart + 4)" for (
                rowstart,
                colstart,
            ) in zip(
                1:(size(mat, 1) - 4),
                1:(size(mat, 2) - 4),
            )
                rI, rJ = rowstart:(rowstart + 4), colstart:(colstart + 4)
                submat = mat[rI, rJ]
                issparse ? nonzeros(submat) .= rand(eltype(submat), nnz(submat)) : submat .= rand(eltype(submat), size(submat)...)

                dset[rI, rJ] = mat[rI, rJ] = submat
                dset_t[rJ, rI] = mat_t[rJ, rI] = submat'

                @test dset == mat
                @test dset_t == mat_t

                if issparse
                    @test_throws KeyError dset[rI, rJ] =
                        rand(eltype(dset), (length(rI), length(rJ)))
                    @test_throws KeyError dset_t[rJ, rI] =
                        rand(eltype(dset_t), (length(rJ), length(rI)))
                end
            end
        end
    end
end
close(file)


text:
using ValSplit, Test

## Test basic functionality

# Define soundof(::Val) methods for :dog and :cat
soundof(animal::Val{:dog}) = "woof"
soundof(animal::Val{:cat}) = "nyan"

# Test first version of macro
@valsplit function soundof(Val(animal::Symbol))
    error("Sound not defined for animal: $animal")
end

@test soundof(:dog) == "woof"
@test soundof(:cat) == "nyan"

# soundof(::Symbol) should error for values without corresponding methods
@test_throws ErrorException soundof(:human) == "meh"

# soundof(::Symbol) should recompile after new methods are defined
soundof(animal::Val{:human}) = "meh"
@test soundof(:human) == "meh"

# Test macro with index argument
@valsplit 1 function soundof(animal::Symbol)
    return "???"
end

@test soundof(:unknown) == "???"

## Test splitting on multiple arguments

# Define soundof(animal::Val, lang::Val) methods
soundof(animal::Val{:cat}, lang::Val{:japanese}) = "nyan"
soundof(animal::Val{:cat}, lang::Val{:indonesian}) = "meong"

soundof(animal::Val{:frog}, lang::Val{:korean}) = "gaegul"
soundof(animal::Val{:frog}, lang::Val{:hindi}) = "tarr"
soundof(animal::Val{:frog}, lang::Val{:english}) = "ribbit"

@valsplit function soundof(Val(animal::Symbol), Val(lang::Symbol))
    error("Sound not defined.")
end

@test soundof(:cat, :japanese) == "nyan"
@test soundof(:cat, :indonesian) == "meong"

@test soundof(:frog, :korean) == "gaegul"
@test soundof(:frog, :hindi) == "tarr"
@test soundof(:frog, :english) == "ribbit"

@test_throws ErrorException soundof(:cat, :english) == "meow"
@test_throws ErrorException soundof(:dog, :english) == "woof"

soundof(::Val{:cat}, ::Val{:english}) = "meow"
soundof(::Val{:dog}, ::Val{:english}) = "woof"

@test soundof(:cat, :english) == "meow"
@test soundof(:dog, :english) == "woof"

## Test non-symbol method signatures

@valsplit function spelling(Val(n::T)) where {T <: Real}
    error("Unknown number.")
end

spelling(::Val{1}) = "one"
spelling(::Val{2}) = "two"
spelling(::Val{1.0}) = "one point zero"

@test spelling(1) == "one"
@test spelling(2) == "two"
@test spelling(1.0) == "one point zero"

@test_throws ErrorException spelling(π) == "pi"
@test_throws ErrorException spelling(ℯ) == "euler's number"

spelling(::Val{π}) = "pi"
spelling(::Val{ℯ}) = "euler's number"

@test spelling(π) == "pi"
@test spelling(ℯ) == "euler's number"

# Test that we can't split Vararg arguments
@test_throws LoadError @macroexpand @valsplit function spelling(Val(n::Real...))
    string(n)
end

## Test closures

struct Exponentiator{T <: Real}
    base::T
end

(ex::Exponentiator)(exponent::Val{1}) = ex.base
(ex::Exponentiator)(exponent::Val{2}) = ex.base * ex.base

@valsplit (ex::Exponentiator)(Val(exponent::Int)) =
    error("Unsupported exponent: $exponent")

@test Exponentiator(3)(1) == 3
@test Exponentiator(3)(2) == 9
@test_throws ErrorException Exponentiator(3)(3) == 27

(ex::Exponentiator)(exponent::Val{3}) = ex.base * ex.base * ex.base

@test Exponentiator(3)(3) == 27

## Test that we can extend functions defined in another module

module Animals
    using ValSplit
    soundof(animal::Val{:dog}) = "woof"
    soundof(animal::Val{:cat}) = "nyan"
    @valsplit function soundof(Val(animal::Symbol))
        error("Sound not defined for animal: $animal")
    end
end

using .Animals

@test Animals.soundof(:dog) == "woof"
@test Animals.soundof(:cat) == "nyan"
@test_throws ErrorException Animals.soundof(:human) == "meh"

Animals.soundof(animal::Val{:human}) = "meh"

@test Animals.soundof(:human) == "meh"

## Test utilities

# Test valarg_params
@test Set(valarg_params(soundof, Tuple{Any}, 1, Symbol)) ==
    Set((:dog, :cat, :human))
@test Set(valarg_params(soundof, Tuple{Val{:frog}, Any}, 2, Symbol)) ==
    Set((:korean, :hindi, :english))
@test Set(valarg_params(soundof, Tuple{Any, Val{:english}}, (1,2), Symbol)) ==
    Set(((:frog, :english), (:cat, :english), (:dog, :english)))
@test Set(valarg_params(soundof, Tuple{Any}, (1,), Tuple{Symbol})) ==
    Set(((:dog,), (:cat,), (:human,)))

# Test valarg_has_param
@test valarg_has_param(π, spelling, 1, Irrational)
@test valarg_has_param(:japanese, soundof, Tuple{Val{:cat}, Any}, 2)
@test valarg_has_param((:frog, :english), soundof, Tuple{Any, Any}, (1, 2))

ext: jl


text:
import LinearAlgebra.isposdef
import Base.in
### Positive semidefinite cone constraint

# TODO: Terrible documentation. Please fix.
struct SDPConstraint <: Constraint
    head::Symbol
    id_hash::UInt64
    child::AbstractExpr
    size::Tuple{Int, Int}
    dual::ValueOrNothing

    function SDPConstraint(child::AbstractExpr)
        sz = child.size
        if sz[1] != sz[2]
            error("Positive semidefinite expressions must be square")
        end
        id_hash = hash((child, :sdp))
        return new(:sdp, id_hash, child, sz, nothing)
    end
end

function vexity(c::SDPConstraint)
    vex = vexity(c.child)
    if vex == AffineVexity() || vex == ConstVexity()
        return AffineVexity()
    else
        return NotDcp()
    end
end

# users specify SDPs as `A in :SDP` where A is an n x n square matrix
# solvers (Mosek and SCS) specify the *lower triangular part of A* is in the SDP cone
# so we need the lower triangular part (n*(n+1)/2 entries) of c.child to be in the SDP cone
# and we need the corresponding upper elements to match the lower elements,
# which we enforce via equality constraints
function conic_form!(c::SDPConstraint, unique_conic_forms::UniqueConicForms)
    # TODO:
        # 1) propagate dual values
        # 2) presolve to eliminate (many) variables --- variables on upper triangular part often don't matter at all
    if !has_conic_form(unique_conic_forms, c)
        n = c.size[1]
        # construct linear indices to pick out the lower triangular part (including diagonal),
        # the upper triangular part (not including diagonal)
        # and the corresponding entries in the lower triangular part, so
        # symmetry => c.child[upperpart]
        rescale = triu(ones(n,n))
        rescale[diagind(n, n)] .= 1.0
        diagandupperpart = findall(!iszero, vec(rescale))
        lowerpart = Vector{Int}(undef, div(n*(n-1),2))
        upperpart = Vector{Int}(undef, div(n*(n-1),2))
        kupper = 0
        # diagandupperpart in column-major order:
        # consider using find(triu(ones(3,3)))
        @inbounds for j = 1:n
            for i = j+1:n
                kupper += 1
                upperpart[kupper] = n*(i-1) + j # (j,i)th element
                lowerpart[kupper] =n*(j-1) + i # (i,j)th element
            end
        end
        objective = conic_form!((broadcast(*,rescale,c.child))[diagandupperpart], unique_conic_forms)
        sdp_constraint = ConicConstr([objective], :SDP, [div(n*(n+1),2)])
        cache_conic_form!(unique_conic_forms, c, sdp_constraint)
        # make sure upper and lower triangular part match in the solution
        # note that this introduces all-zero rows into the constraint matrix
        # if the matrix we require to be PSD has already been constructed to be symmetric
        equality_constraint = conic_form!(c.child[lowerpart] == c.child[upperpart], unique_conic_forms)

        # for computing duals --- won't work b/c
            # 1) sizes are wrong (will be n(n+1)/2, expects n^2), and
            # 2) the "correct" dual is the sum of the sdp dual and the equality constraint dual (so yes, dual of sdp constraint is *not* symmetric)
        # conic_constr_to_constr[sdp_constraint] = c

    end
    return get_conic_form(unique_conic_forms, c)
end

# TODO: Remove isposdef, change tests to use in. Update documentation and notebooks
function isposdef(x::AbstractExpr)
    if sign(x) == ComplexSign()
        SDPConstraint([real(x) -imag(x);imag(x) real(x)])
    else
        SDPConstraint(x)
    end
end

# TODO: Throw error if symbol is invalid.
function in(x::AbstractExpr, y::Symbol)
    if y == :semidefinite || y == :SDP
        if sign(x) == ComplexSign()
            SDPConstraint([real(x) -imag(x);imag(x) real(x)])
        else
            SDPConstraint(x)
        end
    end
end

function ⪰(x::AbstractExpr, y::AbstractExpr)
    if sign(x) == ComplexSign() || sign(y) == ComplexSign()
        SDPConstraint([real(x-y) -imag(x-y);imag(x-y) real(x-y)])
    else
        SDPConstraint(x-y)
    end
end

function ⪯(x::AbstractExpr, y::AbstractExpr)
    if sign(x) == ComplexSign() || sign(y) == ComplexSign()
        SDPConstraint([real(y-x) -imag(y-x);imag(y-x) real(y-x)])
    else
        SDPConstraint(y-x)
    end
end

function ⪰(x::AbstractExpr, y::Value)
    if sign(x) == ComplexSign() || !isreal(y)
        all(y .== 0) ? SDPConstraint([real(x) -imag(x);imag(x) real(x)]) : SDPConstraint([real(x-Constant(y)) -imag(x-Constant(y));imag(x-Constant(y)) real(x-Constant(y))])
    else
        all(y .== 0) ? SDPConstraint(x) : SDPConstraint(x - Constant(y))
    end

end

function ⪰(x::Value, y::AbstractExpr)
    if sign(y) == ComplexSign() || !isreal(x)
        all(x .== 0) ? SDPConstraint([real(-y) -imag(-y);imag(-y) real(-y)]) : SDPConstraint([real(Constant(x)-y) -imag(Constant(x)-y);imag(Constant(x)-y) real(Constant(x)-y)])
    else
        all(x .== 0) ? SDPConstraint(-y) : SDPConstraint(Constant(x) - y)
    end
end

function ⪯(x::Value, y::AbstractExpr)
    if sign(y) == ComplexSign() || !isreal(x)
        all(x .== 0) ? SDPConstraint([real(y) -imag(y);imag(y) real(y)]) : SDPConstraint([real(y-Constant(x)) -imag(y-Constant(x));imag(y-Constant(x)) real(y-Constant(x))])
    else
        all(x .== 0) ? SDPConstraint(y) : SDPConstraint(y - Constant(x))
    end
end

function ⪯(x::AbstractExpr, y::Value)
    if sign(x) == ComplexSign() || !isreal(y)
        all(y .== 0) ? SDPConstraint([real(-x) -imag(-x);imag(-x) real(-x)]) : SDPConstraint([real(Constant(y)-x) -imag(Constant(y)-x);imag(Constant(y)-x) real(Constant(y)-x)])
    else
        all(y .== 0) ? SDPConstraint(-x) : SDPConstraint(Constant(y) - x)
    end
end


text:
sample = "00100
11110
10110
10111
10101
01111
00111
11100
10000
11001
00010
01010"

function flip(bin)
    s = ""
    for i in bin
        if i == '1'
            s = s * "0"
        else
            s = s * "1"
        end
    end
    return s
end

function sol1(sample)
    bit = ""
    for i in 1:sizeof(sample[1])
        most = ""

        for j in 1:size(sample)[1]
            most = most * sample[j][i]
        end

        if count(i->(i=='0'), most) >= size(sample)[1]/2
            bit = bit * "0"
        else
            bit = bit * "1"
        end
    end

    return parse(Int, bit, base=2) * parse(Int, flip(bit), base=2)
end

println("Part 1: ", sol1(split(sample, "\n")))


text:
using Statistics:mean, std
using SpecialFunctions:erf

"""
normalize(scoresTrain,
          scoresTest)

Transform an array of scores into a range between [0,1] using min-max scaling.

Parameters
----------
      scores_train::AbstractVector{<:Real}
A vector of training scores, typically the result of [`fit`](@ref) with a detector.

    scores_test::AbstractVector{<:Real}
A vector of test scores, typically the result of [`score`](@ref) using a previously fitted detector.

Returns
----------
normalized_scores::Tuple{AbstractVector{<:Real}, AbstractVector{<:Real}}
The normalized train and test scores.

Examples
----------
scores_train, scores_test = ([1, 2, 3], [4, 3, 2, 1, 0])
scale_minmax(scores_train, scores_test) # ([0.0, 0.5, 1.0], [1.0, 1.0, 0.5, 0.0, 0.0])
scale_minmax(scores_train) # [0.0, 0.5, 1.0]
"""
function scale_minmax(scores_train::Scores, scores_test::Scores)::Tuple{Scores, Scores}
    minTrain, maxTrain = extrema(scores_train)
    @assert minTrain < maxTrain "Cannot normalize scores if they are all equal"
    f = scores -> clamp.((scores .- minTrain) ./ (maxTrain - minTrain), 0, 1)
    f(scores_train), f(scores_test)
end

"""
unify(scores_train,
      scores_test)

Transform an array of scores into a range between [0,1] using unifying scores as described in [1].

Parameters
----------
    scores_train::AbstractVector{<:Real}
A vector of training scores, typically the result of [`fit`](@ref) with a detector.

    scores_test::AbstractVector{<:Real}
A vector of test scores, typically the result of [`score`](@ref) using a previously fitted detector.

Returns
----------
unified_scores::Tuple{AbstractVector{<:Real}, AbstractVector{<:Real}}
The unified train and test scores.

Examples
----------
scores_train, scores_test = ([1, 2, 3], [4, 3, 2, 1, 0])
unify(scores_train, scores_test) # ([0.0, 0.0, 0.68..], [0.95.., 0.68.., 0.0, 0.0, 0.0])
unify(scores_train) # [0.0, 0.0, 0.68..]

References
----------
Kriegel, Hans-Peter; Kroger, Peer; Schubert, Erich; Zimek, Arthur (2011): Interpreting and Unifying Outlier Scores.
"""
function scale_unify(scores_train::Scores, scores_test::Scores)::Tuple{Scores, Scores}
    μ, σ = mean(scores_train), std(scores_train)
    @assert σ > 0 "Cannot normalize scores if they are all equal"
    f = scores -> clamp.(erf.((scores .- μ) ./ (σ * √2)), 0, 1)
    f(scores_train), f(scores_test)
end


text:
export HessianPreconditioner,getSSORCGRegularizationPreconditioner,getSSORRegularizationPreconditioner,getExactSolveRegularizationPreconditioner
export getEmptyRegularizationPreconditioner
############ General Hessian Preconditioner:
mutable struct HessianPreconditioner
	param
	applyPrec::Function # a function with parameters: Hs::function, d2R::SparseMatrixCSC,v::Vector,param::Any
	setupPrec::Function # a function with parameters: Hs::function, d2R::SparseMatrixCSC,v::Vector,param::Any
end
######################################################################################################################
############ SSORCG: A preconditioner that inverts the regularization matrix by SSOR - CG ############################
######################################################################################################################

mutable struct SSORCGParam
	diagonal	::Vector
	auxVec		::Vector
	omega		::Float64
	tol			::Float64
	maxCGIter	::Int64
end

function getSSORCGRegularizationPreconditioner(omega::Float64=1.0,tol::Float64=1e-2,maxCGIter::Int64=100)
	return HessianPreconditioner(SSORCGParam([],[],omega,tol,maxCGIter),applySSORCG,setupSSORCG);
end

function applySSORCG(Hs::Function, d2R::SparseMatrixCSC,v::Vector,param)
	aux = param.auxVec;
	SSOR(r) = (aux[:].=0.0; xt=ssorPrecTrans!(d2R,aux,r,param.diagonal); return aux);
	x = KrylovMethods.cg(d2R,v,tol=param.tol,maxIter=param.maxCGIter,M=SSOR,out=-1)[1]
	return x;
end

function setupSSORCG(Hs::Function, d2R::SparseMatrixCSC,param)
	param.diagonal = param.omega./diag(d2R);
	param.auxVec   = zeros(size(d2R,2));
	return;
end

######################################################################################################################
############ SSOR: A preconditioner that inverts the regularization matrix by SSOR only ##############################
######################################################################################################################

mutable struct SSORParam
	omega		::Float64
	tol			::Float64
	maxIter	::Int64
end

function getSSORRegularizationPreconditioner(omega::Float64=1.0,tol::Float64=1e-2,maxIter::Int64=100)
	return HessianPreconditioner(SSORParam(omega,tol,maxIter),applySSOR,setupSSOR);
end

function applySSOR(Hs::Function, d2R::SparseMatrixCSC,v::Vector,param)
	x = KrylovMethods.ssor(d2R,copy(v);x=[],tol=param.tol,maxIter=param.maxIter,omega=param.omega,out=-1,storeInterm=false)[1];
	return x;
end

function setupSSOR(Hs::Function, d2R::SparseMatrixCSC,param)
	return;
end

######################################################################################################################
############ ExactSolveRegularization: A preconditioner that inverts the regularization matrix by a backslash operator
######################################################################################################################

function getExactSolveRegularizationPreconditioner()
	return HessianPreconditioner([],applyExactRegSolve,setupExactRegSolve);
end

function applyExactRegSolve(Hs::Function, d2R::SparseMatrixCSC,v::Vector,param)
	return d2R\v;
end

function setupExactRegSolve(Hs::Function, d2R::SparseMatrixCSC,param)
	return;
end

########################################################################################
function getEmptyRegularizationPreconditioner()
	return HessianPreconditioner([],applyEmptyRegSolve,setupEmptyRegSolve);
end

function applyEmptyRegSolve(Hs::Function, d2R::SparseMatrixCSC,v::Vector,param)
	return copy(v);
end

function setupEmptyRegSolve(Hs::Function, d2R::SparseMatrixCSC,param)
	return;
end


text:
mutable struct AB3Cache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  ralk2::rateType
  k::rateType
  tmp::uType
  step::Int
end

u_cache(c::AB3Cache) = ()
du_cache(c::AB3Cache) = (c.fsalfirst,c.k2,c.k3,c.ralk2,c.k)

mutable struct AB3ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
  k2::rateType
  k3::rateType
  step::Int
end

function alg_cache(alg::AB3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  fsalfirst = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  ralk2 = zero(rate_prototype)
  k  = zero(rate_prototype)
  tmp = similar(u)
  AB3Cache(u,uprev,fsalfirst,k2,k3,ralk2,k,tmp,1)
end

function alg_cache(alg::AB3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  k2 = rate_prototype
  k3 = rate_prototype
  AB3ConstantCache(k2,k3,1)
end

mutable struct ABM32Cache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  ralk2::rateType
  k::rateType
  tmp::uType
  step::Int
end

u_cache(c::ABM32Cache) = ()
du_cache(c::ABM32Cache) = (c.fsalfirst,c.k2,c.k3,c.ralk2,c.k)

mutable struct ABM32ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
  k2::rateType
  k3::rateType
  step::Int
end

function alg_cache(alg::ABM32,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  fsalfirst = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  ralk2 = zero(rate_prototype)
  k  = zero(rate_prototype)
  tmp = similar(u)
  ABM32Cache(u,uprev,fsalfirst,k2,k3,ralk2,k,tmp,1)
end

function alg_cache(alg::ABM32,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  k2 = rate_prototype
  k3 = rate_prototype
  ABM32ConstantCache(k2,k3,1)
end

mutable struct AB4Cache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  ralk2::rateType
  k::rateType
  tmp::uType
  t2::rateType
  t3::rateType
  t4::rateType
  step::Int
end

u_cache(c::AB4Cache) = ()
du_cache(c::AB4Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.ralk2,c.k)

mutable struct AB4ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
  k2::rateType
  k3::rateType
  k4::rateType
  step::Int
end

function alg_cache(alg::AB4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  fsalfirst = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  ralk2 = zero(rate_prototype)
  k  = zero(rate_prototype)
  tmp = similar(u)
  t2 = zero(rate_prototype)
  t3 = zero(rate_prototype)
  t4 = zero(rate_prototype)
  AB4Cache(u,uprev,fsalfirst,k2,k3,k4,ralk2,k,tmp,t2,t3,t4,1)
end

function alg_cache(alg::AB4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  k2 = rate_prototype
  k3 = rate_prototype
  k4 = rate_prototype
  AB4ConstantCache(k2,k3,k4,1)
end

mutable struct ABM43Cache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  ralk2::rateType
  k::rateType
  tmp::uType
  t2::rateType
  t3::rateType
  t4::rateType
  t5::rateType
  t6::rateType
  t7::rateType
  step::Int
end

u_cache(c::ABM43Cache) = ()
du_cache(c::ABM43Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.ralk2,c.k)

mutable struct ABM43ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
  k2::rateType
  k3::rateType
  k4::rateType
  step::Int
end

function alg_cache(alg::ABM43,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  fsalfirst = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  ralk2 = zero(rate_prototype)
  k  = zero(rate_prototype)
  tmp = similar(u)
  t2 = zero(rate_prototype)
  t3 = zero(rate_prototype)
  t4 = zero(rate_prototype)
  t5 = zero(rate_prototype)
  t6 = zero(rate_prototype)
  t7 = zero(rate_prototype)
  ABM43Cache(u,uprev,fsalfirst,k2,k3,k4,ralk2,k,tmp,t2,t3,t4,t5,t6,t7,1)
end

function alg_cache(alg::ABM43,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  k2 = rate_prototype
  k3 = rate_prototype
  k4 = rate_prototype
  ABM43ConstantCache(k2,k3,k4,1)
end

mutable struct AB5Cache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k::rateType
  tmp::uType
  t2::rateType
  t3::rateType
  t4::rateType
  step::Int
end

u_cache(c::AB5Cache) = ()
du_cache(c::AB5Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k5,c.k)

mutable struct AB5ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  step::Int
end

function alg_cache(alg::AB5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  fsalfirst = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k  = zero(rate_prototype)
  tmp = similar(u)
  t2 = zero(rate_prototype)
  t3 = zero(rate_prototype)
  t4 = zero(rate_prototype)
  AB5Cache(u,uprev,fsalfirst,k2,k3,k4,k5,k,tmp,t2,t3,t4,1)
end

function alg_cache(alg::AB5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  k2 = rate_prototype
  k3 = rate_prototype
  k4 = rate_prototype
  k5 = rate_prototype
  AB5ConstantCache(k2,k3,k4,k5,1)
end

mutable struct ABM54Cache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k::rateType
  tmp::uType
  t2::rateType
  t3::rateType
  t4::rateType
  t5::rateType
  t6::rateType
  t7::rateType
  t8::rateType
  step::Int
end

u_cache(c::ABM54Cache) = ()
du_cache(c::ABM54Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k5,c.k)

mutable struct ABM54ConstantCache{rateType} <: OrdinaryDiffEqConstantCache
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  step::Int
end

function alg_cache(alg::ABM54,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  fsalfirst = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k  = zero(rate_prototype)
  tmp = similar(u)
  t2 = zero(rate_prototype)
  t3 = zero(rate_prototype)
  t4 = zero(rate_prototype)
  t5 = zero(rate_prototype)
  t6 = zero(rate_prototype)
  t7 = zero(rate_prototype)
  t8 = zero(rate_prototype)
  ABM54Cache(u,uprev,fsalfirst,k2,k3,k4,k5,k,tmp,t2,t3,t4,t5,t6,t7,t8,1)
end

function alg_cache(alg::ABM54,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  k2 = rate_prototype
  k3 = rate_prototype
  k4 = rate_prototype
  k5 = rate_prototype
  ABM54ConstantCache(k2,k3,k4,k5,1)
end

mutable struct VCAB3ConstantCache{TabType,tArrayType,rArrayType,cArrayType,dtArrayType} <: OrdinaryDiffEqConstantCache
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::rArrayType
  ϕstar_nm1::rArrayType
  ϕstar_n::rArrayType
  β::tArrayType
  order::Int
  tab::TabType
  step::Int
end

mutable struct VCAB3Cache{uType,rateType,TabType,uArrayType,bs3Type,tArrayType,cArrayType,uEltypeNoUnits,coefType,dtArrayType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  bs3cache::bs3Type
  k4::rateType
  ϕstar_nm1::coefType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::coefType
  ϕstar_n::coefType
  β::tArrayType
  order::Int
  atmp::uEltypeNoUnits
  tmp::uType
  utilde::uArrayType
  tab::TabType
  step::Int
end

u_cache(c::VCAB3Cache) = ()
du_cache(c::VCAB3Cache) = ()

function alg_cache(alg::VCAB3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  dts = fill(zero(typeof(dt)),3)
  c = fill(zero(typeof(t)), 3, 3)
  g = fill(zero(typeof(t)), 3)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 3)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 3)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 3)
  for i in 1:3
    ϕ_n[i] = copy(rate_prototype)
    ϕstar_nm1[i] = copy(rate_prototype)
    ϕstar_n[i] = copy(rate_prototype)
  end
  β = fill(zero(typeof(t)),3)
  order = 3
  tab = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  VCAB3ConstantCache(dts,c,g,ϕ_n,ϕstar_nm1,ϕstar_n,β,order,tab,1)
end

function alg_cache(alg::VCAB3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  bk1 = zero(rate_prototype)
  bk2 = zero(rate_prototype)
  bk3 = zero(rate_prototype)
  bk4 = zero(rate_prototype)
  butilde = similar(u,axes(u))
  batmp = similar(u,uEltypeNoUnits)
  btmp = similar(u)
  bs3cache = BS3Cache(u,uprev,bk1,bk2,bk3,bk4,butilde,btmp,batmp,tab)
  fsalfirst = zero(rate_prototype)
  k4 = zero(rate_prototype)
  dts = fill(zero(typeof(dt)),3)
  c = fill(zero(typeof(t)), 3, 3)
  g = fill(zero(typeof(t)), 3)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 3)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 3)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 3)
  for i in 1:3
    ϕ_n[i] = zero(rate_prototype)
    ϕstar_nm1[i] = zero(rate_prototype)
    ϕstar_n[i] = zero(rate_prototype)
  end
  β = fill(zero(typeof(t)),3)
  order = 3
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  utilde = similar(u,axes(u))
  VCAB3Cache(u,uprev,fsalfirst,bs3cache,k4,ϕstar_nm1,dts,c,g,ϕ_n,ϕstar_n,β,order,atmp,tmp,utilde,tab,1)
end

mutable struct VCAB4ConstantCache{rk4constcache,tArrayType,rArrayType,cArrayType,dtArrayType} <: OrdinaryDiffEqConstantCache
  ϕstar_nm1::rArrayType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::rArrayType
  ϕstar_n::rArrayType
  β::tArrayType
  order::Int
  rk4constcache::rk4constcache
  step::Int
end

mutable struct VCAB4Cache{uType,rateType,uArrayType,rk4cacheType,tArrayType,cArrayType,uEltypeNoUnits,coefType,dtArrayType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  rk4cache::rk4cacheType
  k4::rateType
  ϕstar_nm1::coefType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::coefType
  ϕstar_n::coefType
  β::tArrayType
  order::Int
  atmp::uEltypeNoUnits
  tmp::uType
  utilde::uArrayType
  step::Int
end

u_cache(c::VCAB4Cache) = ()
du_cache(c::VCAB4Cache) = ()

function alg_cache(alg::VCAB4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  dts = fill(zero(typeof(dt)),4)
  c = fill(zero(typeof(t)), 4, 4)
  g = fill(zero(typeof(t)), 4)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 4)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 4)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 4)
  for i in 1:4
    ϕ_n[i] = copy(rate_prototype)
    ϕstar_nm1[i] = copy(rate_prototype)
    ϕstar_n[i] = copy(rate_prototype)
  end
  β = fill(zero(typeof(t)),4)
  order = 4
  rk4constcache = RK4ConstantCache()
  VCAB4ConstantCache(ϕstar_nm1,dts,c,g,ϕ_n,ϕstar_n,β,order,rk4constcache,1)
end

function alg_cache(alg::VCAB4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  rk1 = zero(rate_prototype)
  rk2 = zero(rate_prototype)
  rk3 = zero(rate_prototype)
  rk4 = zero(rate_prototype)
  rk  = zero(rate_prototype)
  rtmp = similar(u); ratmp = similar(u, uEltypeNoUnits)
  rk4cache = RK4Cache(u,uprev,rk1,rk2,rk3,rk4,rk,rtmp,ratmp)
  fsalfirst = zero(rate_prototype)
  k4 = zero(rate_prototype)
  dts = fill(zero(typeof(dt)),4)
  c = fill(zero(typeof(t)), 4, 4)
  g = fill(zero(typeof(t)), 4)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 4)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 4)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 4)
  for i in 1:4
    ϕ_n[i] = zero(rate_prototype)
    ϕstar_nm1[i] = zero(rate_prototype)
    ϕstar_n[i] = zero(rate_prototype)
  end
  β = fill(zero(typeof(t)),4)
  order = 4
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  utilde = similar(u,axes(u))
  VCAB4Cache(u,uprev,fsalfirst,rk4cache,k4,ϕstar_nm1,dts,c,g,ϕ_n,ϕstar_n,β,order,atmp,tmp,utilde,1)
end

# VCAB5

mutable struct VCAB5ConstantCache{rk4constcache,tArrayType,rArrayType,cArrayType,dtArrayType} <: OrdinaryDiffEqConstantCache
  ϕstar_nm1::rArrayType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::rArrayType
  ϕstar_n::rArrayType
  β::tArrayType
  order::Int
  rk4constcache::rk4constcache
  step::Int
end

mutable struct VCAB5Cache{uType,rateType,uArrayType,rk4cacheType,tArrayType,cArrayType,uEltypeNoUnits,coefType,dtArrayType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  rk4cache::rk4cacheType
  k4::rateType
  ϕstar_nm1::coefType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::coefType
  ϕstar_n::coefType
  β::tArrayType
  order::Int
  atmp::uEltypeNoUnits
  tmp::uType
  utilde::uArrayType
  step::Int
end

u_cache(c::VCAB5Cache) = ()
du_cache(c::VCAB5Cache) = ()

function alg_cache(alg::VCAB5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  dts = fill(zero(typeof(dt)),5)
  c = fill(zero(typeof(t)), 5, 5)
  g = fill(zero(typeof(t)), 5)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 5)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 5)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 5)
  for i in 1:5
    ϕ_n[i] = copy(rate_prototype)
    ϕstar_nm1[i] = copy(rate_prototype)
    ϕstar_n[i] = copy(rate_prototype)
  end
  β = fill(zero(typeof(t)),5)
  order = 5
  rk4constcache = RK4ConstantCache()
  VCAB5ConstantCache(ϕstar_nm1,dts,c,g,ϕ_n,ϕstar_n,β,order,rk4constcache,1)
end

function alg_cache(alg::VCAB5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  rk1 = zero(rate_prototype)
  rk2 = zero(rate_prototype)
  rk3 = zero(rate_prototype)
  rk4 = zero(rate_prototype)
  rk  = zero(rate_prototype)
  rtmp = similar(u); ratmp = similar(u, uEltypeNoUnits)
  rk4cache = RK4Cache(u,uprev,rk1,rk2,rk3,rk4,rk,rtmp,ratmp)
  fsalfirst = zero(rate_prototype)
  k4 = zero(rate_prototype)
  dts = fill(zero(typeof(dt)),5)
  c = fill(zero(typeof(t)), 5, 5)
  g = fill(zero(typeof(t)), 5)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 5)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 5)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 5)
  for i in 1:5
    ϕ_n[i] = zero(rate_prototype)
    ϕstar_nm1[i] = zero(rate_prototype)
    ϕstar_n[i] = zero(rate_prototype)
  end
  β = fill(zero(typeof(t)),5)
  order = 5
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  utilde = similar(u,axes(u))
  VCAB5Cache(u,uprev,fsalfirst,rk4cache,k4,ϕstar_nm1,dts,c,g,ϕ_n,ϕstar_n,β,order,atmp,tmp,utilde,1)
end

# VCABM3

mutable struct VCABM3ConstantCache{TabType,tArrayType,rArrayType,cArrayType,dtArrayType} <: OrdinaryDiffEqConstantCache
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::rArrayType
  ϕ_np1::rArrayType
  ϕstar_nm1::rArrayType
  ϕstar_n::rArrayType
  β::tArrayType
  order::Int
  tab::TabType
  step::Int
end

mutable struct VCABM3Cache{uType,rateType,TabType,uArrayType,bs3Type,tArrayType,cArrayType,uEltypeNoUnits,coefType,dtArrayType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  bs3cache::bs3Type
  k4::rateType
  ϕstar_nm1::coefType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::coefType
  ϕ_np1::coefType
  ϕstar_n::coefType
  β::tArrayType
  order::Int
  atmp::uEltypeNoUnits
  tmp::uType
  utilde::uArrayType
  tab::TabType
  step::Int
end

u_cache(c::VCABM3Cache) = ()
du_cache(c::VCABM3Cache) = ()

function alg_cache(alg::VCABM3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  dts = fill(zero(typeof(dt)),3)
  c = fill(zero(typeof(t)), 4, 4)
  g = fill(zero(typeof(t)), 4)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 3)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 3)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 3)
  ϕ_np1 = Vector{typeof(rate_prototype)}(undef, 4)
  for i in 1:3
    ϕ_n[i] = copy(rate_prototype)
    ϕstar_nm1[i] = copy(rate_prototype)
    ϕstar_n[i] = copy(rate_prototype)
  end
  β = fill(zero(typeof(t)),3)
  order = 3
  tab = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  VCABM3ConstantCache(dts,c,g,ϕ_n,ϕ_np1,ϕstar_nm1,ϕstar_n,β,order,tab,1)
end

function alg_cache(alg::VCABM3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  bk1 = zero(rate_prototype)
  bk2 = zero(rate_prototype)
  bk3 = zero(rate_prototype)
  bk4 = zero(rate_prototype)
  butilde = similar(u,axes(u))
  batmp = similar(u,uEltypeNoUnits)
  btmp = similar(u)
  bs3cache = BS3Cache(u,uprev,bk1,bk2,bk3,bk4,butilde,btmp,batmp,tab)
  fsalfirst = zero(rate_prototype)
  k4 = zero(rate_prototype)
  dts = fill(zero(typeof(dt)),3)
  c = fill(zero(typeof(t)), 4, 4)
  g = fill(zero(typeof(t)), 4)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 3)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 3)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 3)
  ϕ_np1 = Vector{typeof(rate_prototype)}(undef, 4)
  for i in 1:3
    ϕ_n[i] = zero(rate_prototype)
    ϕstar_nm1[i] = zero(rate_prototype)
    ϕstar_n[i] = zero(rate_prototype)
  end
  for i in 1:4
    ϕ_np1[i] = zero(rate_prototype)
  end
  β = fill(zero(typeof(t)),3)
  order = 3
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  utilde = similar(u,axes(u))
  VCABM3Cache(u,uprev,fsalfirst,bs3cache,k4,ϕstar_nm1,dts,c,g,ϕ_n,ϕ_np1,ϕstar_n,β,order,atmp,tmp,utilde,tab,1)
end

# VCABM4

mutable struct VCABM4ConstantCache{rk4constcache,tArrayType,rArrayType,cArrayType,dtArrayType} <: OrdinaryDiffEqConstantCache
  ϕstar_nm1::rArrayType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::rArrayType
  ϕ_np1::rArrayType
  ϕstar_n::rArrayType
  β::tArrayType
  order::Int
  rk4constcache::rk4constcache
  step::Int
end

mutable struct VCABM4Cache{uType,rateType,uArrayType,rk4cacheType,tArrayType,cArrayType,uEltypeNoUnits,coefType,dtArrayType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  rk4cache::rk4cacheType
  k4::rateType
  ϕstar_nm1::coefType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::coefType
  ϕ_np1::coefType
  ϕstar_n::coefType
  β::tArrayType
  order::Int
  atmp::uEltypeNoUnits
  tmp::uType
  utilde::uArrayType
  step::Int
end

u_cache(c::VCABM4Cache) = ()
du_cache(c::VCABM4Cache) = ()

function alg_cache(alg::VCABM4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  dts = fill(zero(typeof(dt)),4)
  c = fill(zero(typeof(t)), 5, 5)
  g = fill(zero(typeof(t)), 5)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 4)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 4)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 4)
  ϕ_np1 = Vector{typeof(rate_prototype)}(undef, 5)
  for i in 1:4
    ϕ_n[i] = copy(rate_prototype)
    ϕstar_nm1[i] = copy(rate_prototype)
    ϕstar_n[i] = copy(rate_prototype)
  end
  β = fill(zero(typeof(t)),4)
  order = 4
  rk4constcache = RK4ConstantCache()
  VCABM4ConstantCache(ϕstar_nm1,dts,c,g,ϕ_n,ϕ_np1,ϕstar_n,β,order,rk4constcache,1)
end

function alg_cache(alg::VCABM4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  rk1 = zero(rate_prototype)
  rk2 = zero(rate_prototype)
  rk3 = zero(rate_prototype)
  rk4 = zero(rate_prototype)
  rk  = zero(rate_prototype)
  rtmp = similar(u); ratmp = similar(u, uEltypeNoUnits)
  rk4cache = RK4Cache(u,uprev,rk1,rk2,rk3,rk4,rk,rtmp,ratmp)
  fsalfirst = zero(rate_prototype)
  k4 = zero(rate_prototype)
  dts = fill(zero(typeof(dt)),4)
  c = fill(zero(typeof(t)), 5, 5)
  g = fill(zero(typeof(t)), 5)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 4)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 4)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 4)
  ϕ_np1 = Vector{typeof(rate_prototype)}(undef, 5)
  for i in 1:4
    ϕ_n[i] = zero(rate_prototype)
    ϕstar_nm1[i] = zero(rate_prototype)
    ϕstar_n[i] = zero(rate_prototype)
  end
  for i in 1:5
    ϕ_np1[i] = zero(rate_prototype)
  end
  β = fill(zero(typeof(t)),4)
  order = 4
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  utilde = similar(u,axes(u))
  VCABM4Cache(u,uprev,fsalfirst,rk4cache,k4,ϕstar_nm1,dts,c,g,ϕ_n,ϕ_np1,ϕstar_n,β,order,atmp,tmp,utilde,1)
end

# VCABM5

mutable struct VCABM5ConstantCache{rk4constcache,tArrayType,rArrayType,cArrayType,dtArrayType} <: OrdinaryDiffEqConstantCache
  ϕstar_nm1::rArrayType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::rArrayType
  ϕ_np1::rArrayType
  ϕstar_n::rArrayType
  β::tArrayType
  order::Int
  rk4constcache::rk4constcache
  step::Int
end

mutable struct VCABM5Cache{uType,rateType,uArrayType,rk4cacheType,tArrayType,cArrayType,uEltypeNoUnits,coefType,dtArrayType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  rk4cache::rk4cacheType
  k4::rateType
  ϕstar_nm1::coefType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::coefType
  ϕ_np1::coefType
  ϕstar_n::coefType
  β::tArrayType
  order::Int
  atmp::uEltypeNoUnits
  tmp::uType
  utilde::uArrayType
  step::Int
end

u_cache(c::VCABM5Cache) = ()
du_cache(c::VCABM5Cache) = ()

function alg_cache(alg::VCABM5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  dts = fill(zero(typeof(dt)),5)
  c = fill(zero(typeof(t)), 6, 6)
  g = fill(zero(typeof(t)), 6)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 5)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 5)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 5)
  ϕ_np1 = Vector{typeof(rate_prototype)}(undef, 6)
  for i in 1:5
    ϕ_n[i] = copy(rate_prototype)
    ϕstar_nm1[i] = copy(rate_prototype)
    ϕstar_n[i] = copy(rate_prototype)
  end
  β = fill(zero(typeof(t)),5)
  order = 5
  rk4constcache = RK4ConstantCache()
  VCABM5ConstantCache(ϕstar_nm1,dts,c,g,ϕ_n,ϕ_np1,ϕstar_n,β,order,rk4constcache,1)
end

function alg_cache(alg::VCABM5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  rk1 = zero(rate_prototype)
  rk2 = zero(rate_prototype)
  rk3 = zero(rate_prototype)
  rk4 = zero(rate_prototype)
  rk  = zero(rate_prototype)
  rtmp = similar(u); ratmp = similar(u, uEltypeNoUnits)
  rk4cache = RK4Cache(u,uprev,rk1,rk2,rk3,rk4,rk,rtmp,ratmp)
  fsalfirst = zero(rate_prototype)
  k4 = zero(rate_prototype)
  dts = fill(zero(typeof(dt)),5)
  c = fill(zero(typeof(t)), 6, 6)
  g = fill(zero(typeof(t)), 6)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 5)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 5)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 5)
  ϕ_np1 = Vector{typeof(rate_prototype)}(undef, 6)
  for i in 1:5
    ϕ_n[i] = zero(rate_prototype)
    ϕstar_nm1[i] = zero(rate_prototype)
    ϕstar_n[i] = zero(rate_prototype)
  end
  for i in 1:6
    ϕ_np1[i] = zero(rate_prototype)
  end
  β = fill(zero(typeof(t)),5)
  order = 5
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  utilde = similar(u,axes(u))
  VCABM5Cache(u,uprev,fsalfirst,rk4cache,k4,ϕstar_nm1,dts,c,g,ϕ_n,ϕ_np1,ϕstar_n,β,order,atmp,tmp,utilde,1)
end

# VCABM

mutable struct VCABMConstantCache{tArrayType,rArrayType,cArrayType,dtType,dtArrayType} <: OrdinaryDiffEqConstantCache
  ϕstar_nm1::rArrayType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::rArrayType
  ϕ_np1::rArrayType
  ϕstar_n::rArrayType
  β::tArrayType
  ξ::dtType
  ξ0::dtType
  order::Int
  max_order::Int
  step::Int
end

mutable struct VCABMCache{uType,rateType,uArrayType,dtType,tArrayType,cArrayType,uEltypeNoUnits,coefType,dtArrayType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k4::rateType
  ϕstar_nm1::coefType
  dts::dtArrayType
  c::cArrayType
  g::tArrayType
  ϕ_n::coefType
  ϕ_np1::coefType
  ϕstar_n::coefType
  β::tArrayType
  order::Int
  max_order::Int
  atmp::uEltypeNoUnits
  tmp::uType
  ξ::dtType
  ξ0::dtType
  utilde::uArrayType
  utildem1::uArrayType
  utildem2::uArrayType
  utildep1::uArrayType
  atmpm1::uEltypeNoUnits
  atmpm2::uEltypeNoUnits
  atmpp1::uEltypeNoUnits
  step::Int
end

u_cache(c::VCABMCache) = ()
du_cache(c::VCABMCache) = ()

function alg_cache(alg::VCABM,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  dts = fill(zero(typeof(dt)),13)
  c = fill(zero(typeof(t)), 13, 13)
  g = fill(zero(typeof(t)), 13)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 13)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 13)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 13)
  ϕ_np1 = Vector{typeof(rate_prototype)}(undef, 14)
  for i in 1:13
    ϕ_n[i] = copy(rate_prototype)
    ϕstar_nm1[i] = copy(rate_prototype)
    ϕstar_n[i] = copy(rate_prototype)
  end
  β = fill(zero(typeof(t)),13)
  ξ = zero(dt)
  ξ0 = zero(dt)
  order = 1
  max_order = 12
  VCABMConstantCache(ϕstar_nm1,dts,c,g,ϕ_n,ϕ_np1,ϕstar_n,β,ξ,ξ0,order,max_order,1)
end

function alg_cache(alg::VCABM,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  fsalfirst = zero(rate_prototype)
  k4 = zero(rate_prototype)
  dts = fill(zero(typeof(dt)),13)
  c = fill(zero(typeof(t)), 13, 13)
  g = fill(zero(typeof(t)), 13)
  ϕ_n = Vector{typeof(rate_prototype)}(undef, 13)
  ϕstar_nm1 = Vector{typeof(rate_prototype)}(undef, 13)
  ϕstar_n = Vector{typeof(rate_prototype)}(undef, 13)
  ϕ_np1 = Vector{typeof(rate_prototype)}(undef, 14)
  for i in 1:13
    ϕ_n[i] = zero(rate_prototype)
    ϕstar_nm1[i] = zero(rate_prototype)
    ϕstar_n[i] = zero(rate_prototype)
  end
  for i in 1:14
    ϕ_np1[i] = zero(rate_prototype)
  end
  β = fill(zero(typeof(t)), 13)
  order = 1
  max_order = 12
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  ξ = zero(dt)
  ξ0 = zero(dt)
  utilde = similar(u,axes(u))
  utildem2 = similar(u,axes(u))
  utildem1 = similar(u,axes(u))
  utildep1 = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits)
  atmpm1 = similar(u,uEltypeNoUnits)
  atmpm2 = similar(u,uEltypeNoUnits)
  atmpp1 = similar(u,uEltypeNoUnits)
  VCABMCache(u,uprev,fsalfirst,k4,ϕstar_nm1,dts,c,g,ϕ_n,ϕ_np1,ϕstar_n,β,order,max_order,atmp,tmp,ξ,ξ0,utilde,utildem1,utildem2,utildep1,atmpm1,atmpm2,atmpp1,1)
end

# IMEX Multistep methods

# CNAB2

mutable struct CNAB2ConstantCache{rateType,F,N,uType,tType} <: OrdinaryDiffEqConstantCache
  k2::rateType
  uf::F
  nlsolve::N
  uprev3::uType
  tprev2::tType
end

mutable struct CNAB2Cache{uType,rateType,uNoUnitsType,J,W,UF,JC,N,tType,F} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  fsalfirst::rateType
  k::rateType
  k1::rateType
  k2::rateType
  du₁::rateType
  du1::rateType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::J
  W::W
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  uprev3::uType
  tprev2::tType
end

u_cache(c::CNAB2Cache)    = ()
du_cache(c::CNAB2Cache)   = ()

function alg_cache(alg::CNAB2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  k2 = rate_prototype
  uprev3 = u
  tprev2 = t

  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,1//2,1,ηold,z₊,dz,tmp,b,k))
  CNAB2ConstantCache(k2,uf,nlsolve,uprev3,tprev2)
end

function alg_cache(alg::CNAB2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  atmp = similar(u,uEltypeNoUnits,axes(u))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  du₁ = zero(rate_prototype)

  uprev3 = similar(u)
  tprev2 = t

  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,1//2,1,ηold,z₊,dz,tmp,b,k))
  CNAB2Cache(u,uprev,uprev2,fsalfirst,k,k1,k2,du₁,du1,z,dz,b,tmp,atmp,J,W,uf,jac_config,linsolve,nlsolve,uprev3,tprev2)
end

# CNLF2

mutable struct CNLF2ConstantCache{rateType,F,N,uType,tType} <: OrdinaryDiffEqConstantCache
  k2::rateType
  uf::F
  nlsolve::N
  uprev2::uType
  uprev3::uType
  tprev2::tType
end

mutable struct CNLF2Cache{uType,rateType,uNoUnitsType,J,W,UF,JC,N,tType,F} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  fsalfirst::rateType
  k::rateType
  k1::rateType
  k2::rateType
  du₁::rateType
  du1::rateType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::J
  W::W
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  uprev3::uType
  tprev2::tType
end

u_cache(c::CNLF2Cache)    = ()
du_cache(c::CNLF2Cache)   = ()

function alg_cache(alg::CNLF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  k2 = rate_prototype
  uprev2 = u
  uprev3 = u
  tprev2 = t

  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,1//1,1,ηold,z₊,dz,tmp,b,k))
  CNLF2ConstantCache(k2,uf,nlsolve,uprev2,uprev3,tprev2)
end

function alg_cache(alg::CNLF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  atmp = similar(u,uEltypeNoUnits,axes(u))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  du₁ = zero(rate_prototype)

  uprev2 = similar(u)
  uprev3 = similar(u)
  tprev2 = t

  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,1//1,1,ηold,z₊,dz,tmp,b,k))
  CNLF2Cache(u,uprev,uprev2,fsalfirst,k,k1,k2,du₁,du1,z,dz,b,tmp,atmp,J,W,uf,jac_config,linsolve,nlsolve,uprev3,tprev2)
end


text:
using StanBase

ProjDir = @__DIR__

bernoulli = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
  real empty[0];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

sm = HelpModel( "bernoulli", bernoulli;
  method = StanBase.Help(:sample))

# Other options for method are:
# :optimize, :diagnose or :variational
# see ?CmdStanModels

rc = stan_help(sm; n_chains=1)

if success(rc)
	run(`cat $(sm.log_file[1])`)
end


text:
function exportstockdata(project)
    stockdata = project.dat["master"]
    timestamp = replace(split(string(now()),".")[1],":"=>"")
    CSV.write("StockExtract_$(timestamp).csv",stockdata)

    sum(ismissing.(stockdata[:,4]))
    prep_df(stockdata,col) = stockdata[findall(isnothing.(stockdata[:,col])),col] .= missing
    prep_df.((stockdata,),collect(4:24))
end


text:
#
# this is a trivial first write-up of NW, Algorithm 17.4
# with the new `update!` mechanism, I'd like to incorporate the
# Lagrange-multiplier and penalty updates into the minimisation
# procedure. But maybe this is not a good idea, to be tested.
#

function optimize( F::DifferentiableFunction,
                   C::EqualityConstraint,
                   x0::AbstractVector;
                   iterations = 5000,
                   c_tol = 1e-6,
                   g_tol = 1e-6,
                   verbose = 1,
                   Optimizer = Optim.ConjugateGradient )

   # initialise
   al = AugmentedLagrangian(F, C, x0)
   eta = al.mu^(-0.1)
   # TODO: this needs to distinguish whether F, C are DifferentiableFunction
   # or TwiceDifferentiableFunction!!!!
   ALobj = DifferentiableFunction( x_ -> evaluate(x_, al),
                                   (x_,g_) -> gradient!(g_, x_, al),
                                   (x_,g_) -> eval_and_grad!(g_, x_, al) )

   x = copy(x0)
   iteration = 0
   if verbose >= 1
      @printf("  it  |     |c|_∞        |g_AL|_∞      μ         η  \n")
      @printf("------|----------------------------------------------\n")
   end
   while iteration < iterations

      # solve the sub-problem
      options = OptimizationOptions(g_tol=g_tol, iterations=iterations - iteration, store_trace = true)
      result = Optim.optimize(ALobj, x, Optimizer(), options)
      # TODO: unclear what to do if this step fails? it could still be ok after
      # we update the lagrange multipliers? For now just print a warning and continue.
      if !Optim.converged(result)
         warn("an inner AL iteration has not converged")
         # return minimizer(result)
      end

      # check for convergence, perform update for multipliers and tolerances
      iteration += Optim.iterations(result)
      x = Optim.minimizer(result)
      nrm_g = Optim.g_norm_trace(result)[end]

      c = C.c(x)    # TODO: this is a superfluous call to C (but #282 would fix this)

      if verbose >= 1
         @printf(" %4d |  %1.4e   %1.4e   %1.2e   %1.2e\n",
                  iteration, norm(x, Inf), nrm_g, al.mu, eta)
      end

      if norm(c, Inf) < c_tol
         return x, al
      elseif norm(c, Inf) <= eta
         al.lambda = al.lambda - al.mu * c
         eta /= al.mu^0.9  # tighten constraint tolerance
      else
         al.mu *= 100.0
         eta = al.mu^(-0.1)
      end
      if al.mu > 1e10
         warn("something horrible is happening")
         return x, al
      end

   end

   # iteration > iterations
   warn("too many iterations in `ConstrainedOptim.optimize` ")
   return x, al
end

max_stars_repo_name: cortner/AugmentedLagrangian

max_issues_repo_name: cortner/AugmentedLagrangian

max_forks_repo_name: cortner/AugmentedLagrangian


text:
@doc raw"""
    cyclic_proximal_point(M, F, proxes, x)

perform a cyclic proximal point algorithm.

# Input

* `M` – a manifold ``\mathcal M``
* `F` – a cost function ``F:\mathcal M→ℝ`` to minimize
* `proxes` – an Array of proximal maps (`Function`s) `(λ,x) -> y` for the summands of ``F``
* `x` – an initial value ``x ∈ \mathcal M``

# Optional
the default values are given in brackets
* `evaluation` – ([`AllocatingEvaluation`](@ref)) specify whether the proximal maps work by allocation (default) form `prox(M, λ, x)`
  or [`MutatingEvaluation`](@ref) in place, i.e. is of the form `prox!(M, y, λ, x)`.
* `evaluation_order` – (`:Linear`) – whether
  to use a randomly permuted sequence (`:FixedRandom`), a per
  cycle permuted sequence (`Random`) or the default linear one.
* `λ` – ( `iter -> 1/iter` ) a function returning the (square summable but not
  summable) sequence of λi
* `stopping_criterion` – ([`StopWhenAny`](@ref)`(`[`StopAfterIteration`](@ref)`(5000),`[`StopWhenChangeLess`](@ref)`(10.0^-8))`) a [`StoppingCriterion`](@ref).
* `return_options` – (`false`) – if activated, the extended result, i.e. the
  complete [`Options`](@ref) are returned. This can be used to access recorded values.
  If set to false (default) just the optimal value `x_opt` if returned
and the ones that are passed to [`decorate_options`](@ref) for decorators.

# Output
* `x_opt` – the resulting (approximately critical) point of gradientDescent
OR
* `options` - the options returned by the solver (see `return_options`)
"""
function cyclic_proximal_point(
    M::AbstractManifold, F::TF, proxes::Union{Tuple,AbstractVector}, x0; kwargs...
) where {TF}
    x_res = copy(M, x0)
    return cyclic_proximal_point!(M, F, proxes, x_res; kwargs...)
end

@doc raw"""
    cyclic_proximal_point!(M, F, proxes, x)

perform a cyclic proximal point algorithm in place of `x`.

# Input

* `M` – a manifold ``\mathcal M``
* `F` – a cost function ``F:\mathcal M→ℝ`` to minimize
* `proxes` – an Array of proximal maps (`Function`s) `(λ,x) -> y` for the summands of ``F``
* `x` – an initial value ``x ∈ \mathcal M``

for all options, see [`cyclic_proximal_point`](@ref).
"""
function cyclic_proximal_point!(
    M::AbstractManifold,
    F::TF,
    proxes::Union{Tuple,AbstractVector},
    x0;
    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
    evaluation_order::Symbol=:Linear,
    stopping_criterion::StoppingCriterion=StopWhenAny(
        StopAfterIteration(5000), StopWhenChangeLess(10.0^-12)
    ),
    λ=i -> 1 / i,
    return_options=false,
    kwargs..., #decorator options
) where {TF}
    p = ProximalProblem(M, F, proxes; evaluation=evaluation)
    o = CyclicProximalPointOptions(x0, stopping_criterion, λ, evaluation_order)

    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end
function initialize_solver!(p::ProximalProblem, o::CyclicProximalPointOptions)
    c = length(p.proximal_maps!!)
    o.order = collect(1:c)
    (o.order_type == :FixedRandom) && shuffle!(o.order)
    return o
end
function step_solver!(p::ProximalProblem, o::CyclicProximalPointOptions, iter)
    c = length(p.proximal_maps!!)
    λi = o.λ(iter)
    for k in o.order
        get_proximal_map!(p, o.x, λi, o.x, k)
    end
    (o.order_type == :Random) && shuffle(o.order)
    return o
end
get_solver_result(o::CyclicProximalPointOptions) = o.x


text:
module AircraftSystems

# import packages
import CCBlade
CC = CCBlade
import Dates
import DelimitedFiles
DF = DelimitedFiles
import FLOWMath
FM = FLOWMath
# import GeometricTools
# GT = GeometricTools
using LaTeXStrings
import LinearAlgebra
LA = LinearAlgebra
import PyPlot
plt = PyPlot
import Splines
import StaticArrays
import Statistics
using Test
import VortexLattice
VL = VortexLattice
import Xfoil
XF = Xfoil

# set path
const topdirectory = normpath(joinpath(@__DIR__, ".."))

# set date
const TODAY = replace(string(Dates.today()),"-" => "")

# load files
subdirectories = [
    "system/",
    "action/fundamental/",
    "action/compound/",
    "postaction/",
    "simulation/",
    "template/",
]
for directory in subdirectories
    for (root, dirs, files) in walkdir(joinpath(topdirectory,"src",directory))
        for file in files
            if splitext(file)[end] == ".jl" #&& file != "TrajOpt2020.jl" && !(file in excluded)
            # println("file = ",joinpath(root, file))
                include(joinpath(root,file))
            end
        end
    end
end

# export packages
export CCBlade, FLOWMath, LaTeXStrings, VortexLattice, Test, Xfoil

# export functions
export runsimulation!

end # module

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
"""
Optimize by randomly generating the candidates.
"""
mutable struct RandomSearcher{S<:SearchSpace} <: AskTellOptimizer
    name::String
    search_space::S
    best_fitness          # FIXME fitness type should be known
    best::Individual

    RandomSearcher(searchSpace::S) where {S<:SearchSpace} =
        new{S}("Random Search", searchSpace, nothing)
end

function ask(rs::RandomSearcher)
    # Just randomly generate a new individual and return it with a dummy index
    # (since we do not have a population there are no indices).
    [Candidate{Float64}(rand_individual(rs.search_space), 1)]
end

function tell!(rs::RandomSearcher, rankedCandidates::Vector{Candidate{F}}) where F
    candidate = rankedCandidates[1]
    if rs.best_fitness == nothing || candidate.fitness < rs.best_fitness
        rs.best = candidate.params
        rs.best_fitness = candidate.fitness
        return 1
    else
        return 0
    end
end

random_search(problem::OptimizationProblem, parameters::Parameters) =
    RandomSearcher(search_space(problem))

random_search(ss::SearchSpace) =
    RandomSearcher(ss)


text:
using Documenter, Stheno

makedocs(
	modules = [Stheno],
    format = :html,
    sitename = "Stheno.jl",
    pages = [
        "Home" => "index.md",
        "BlockArrays extensions" => "block_arrays_ext.md",
        "Interfaces" => "mean_and_kernel_interfaces.md",
    ],
)


text:
module TestFilterCircuit2

using ModiaLang
using Unitful

using DifferentialEquations
@usingModiaPlot
using Test

setLogMerge(false)

include("../models/Electric.jl")

FilterCircuit2 = Model(
    R  = Resistor  | Map(R=100u"Ω"),
    Ri = Resistor  | Map(R=10u"Ω"),
    C  = Capacitor | Map(C=2.5e-3u"F", v=Var(init=0.0u"V")),
    V  = ConstantVoltage | Map(V=10.0u"V"),
    ground = Ground,
    connect = :[
      (V.p, Ri.n)
      (Ri.p, R.p)
      (R.n, C.p)
      (C.n, V.n, ground.p)
    ]
)

filterCircuit2 = @instantiateModel(FilterCircuit2, unitless=true, logCode=true)

simulate!(filterCircuit2, Tsit5(), stopTime = 1.0, requiredFinalStates =[9.736520087687524])
plot(filterCircuit2, ("V.v", "Ri.v", "R.v", "C.v"), figure=1)

end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
function setup_basic()
    # Simple QP problem
    problem = Dict()
    problem[:P] = sparse([11. 0.; 0. 0.])
    problem[:q] = [3.; 4]
    problem[:A] = sparse([-1 0; 0 -1; -1 -3; 2 5; 3 4])
    problem[:u] = [0.; 0.; -15; 100; 80]
    problem[:l] = -Inf * ones(length(problem[:u]))
    problem[:n] = size(problem[:P], 1)
    problem[:m] = size(problem[:A], 1)
    options = Dict(:verbose => false,
                   :eps_abs => 1e-09,
                   :eps_rel => 1e-09,
                   :check_termination => 1,
                   :polish => false,
                   :max_iter => 4000,
                   :rho => 0.1,
                   :adaptive_rho => false,
                   :warm_start => true)
    return problem, options
end

tol = 1e-5

@testset "basic" begin

    @testset "basic_QP" begin
        problem, options = setup_basic()

        model = OSQP.Model()
        OSQP.setup!(model; P=problem[:P], q=problem[:q],
                    A=problem[:A], l=problem[:l], u=problem[:u], options...)
        results = OSQP.solve!(model)

        @test isapprox(norm(results.x - [0.; 5.]), 0., atol=tol)
        @test isapprox(norm(results.y - [1.666666666666; 0.; 1.3333333; 0.; 0.]), 0., atol=tol)
        @test isapprox(results.info.obj_val, 20., atol=tol)

    end

    @testset "update_q" begin
        problem, options = setup_basic()

        model = OSQP.Model()
        OSQP.setup!(model; P=problem[:P], q=problem[:q],
                    A=problem[:A], l=problem[:l], u=problem[:u], options...)

        OSQP.update!(model, q=[10.; 20.])
        results = OSQP.solve!(model)

        @test isapprox(norm(results.x - [0.; 5.]), 0., atol=tol)
        @test isapprox(norm(results.y - [3.33333333; 0.; 6.66666666; 0.; 0.]), 0., atol=tol)
        @test isapprox(results.info.obj_val, 100., atol=tol)

    end

    @testset "update_l" begin
        problem, options = setup_basic()

        model = OSQP.Model()
        OSQP.setup!(model; P=problem[:P], q=problem[:q],
                    A=problem[:A], l=problem[:l], u=problem[:u], options...)

        OSQP.update!(model, l=-100 * ones(problem[:m]))
        results = OSQP.solve!(model)

        @test isapprox(norm(results.x - [0.; 5.]), 0., atol=tol)
        @test isapprox(norm(results.y - [1.6666666666; 0.; 1.333333333333; 0.; 0.]), 0., atol=tol)
        @test isapprox(results.info.obj_val, 20., atol=tol)

    end

    @testset "update_u" begin
        problem, options = setup_basic()

        model = OSQP.Model()
        OSQP.setup!(model; P=problem[:P], q=problem[:q],
                    A=problem[:A], l=problem[:l], u=problem[:u], options...)

        OSQP.update!(model, u=1000 * ones(problem[:m]))
        results = OSQP.solve!(model)

        @test isapprox(norm(results.x - [-1.51515152e-01, -3.33282828e+02]), 0., atol=tol)
        @test isapprox(norm(results.y - [0.; 0.; 1.333333333333; 0.; 0.]), 0., atol=tol)
        @test isapprox(results.info.obj_val, -1333.459595961, atol=tol)

    end

    @testset "update_max_iter" begin
        problem, options = setup_basic()

        model = OSQP.Model()
        OSQP.setup!(model; P=problem[:P], q=problem[:q],
                    A=problem[:A], l=problem[:l], u=problem[:u], options...)

        OSQP.update_settings!(model, max_iter=80)
        results = OSQP.solve!(model)

        @test results.info.status == :Max_iter_reached
    end

    @testset "update_check_termination" begin
        problem, options = setup_basic()

        model = OSQP.Model()
        OSQP.setup!(model; P=problem[:P], q=problem[:q],
                    A=problem[:A], l=problem[:l], u=problem[:u], options...)

        OSQP.update_settings!(model, check_termination=false)
        results = OSQP.solve!(model)

        @test results.info.iter == options[:max_iter]
    end

    @testset "update_rho" begin
        problem, options = setup_basic()

        # Setup default problem
        model = OSQP.Model()
        OSQP.setup!(model; P=problem[:P], q=problem[:q],
                    A=problem[:A], l=problem[:l], u=problem[:u], options...)

        results_default = OSQP.solve!(model)

        # Setup different rho and update to same rho
        new_opts = copy(options)
        new_opts[:rho] = 0.7
        model = OSQP.Model()
        OSQP.setup!(model; P=problem[:P], q=problem[:q],
                    A=problem[:A], l=problem[:l], u=problem[:u], new_opts...)
        OSQP.update_settings!(model, rho=options[:rho])
        results_new_rho = OSQP.solve!(model)

        @test results_default.info.iter == results_new_rho.info.iter
    end

    @testset "time_limit" begin

        problem, options = setup_basic()

        model = OSQP.Model()
        OSQP.setup!(model; P=problem[:P], q=problem[:q],
                    A=problem[:A], l=problem[:l], u=problem[:u],  options...)
        results = OSQP.solve!(model)

        @test results.info.status == :Solved

        # Ensure solver will time out
        OSQP.update_settings!(model, time_limit=1e-6, max_iter=1000000, check_termination=0)

        results_time_limit = OSQP.solve!(model)

        @test results_time_limit.info.status == :Time_limit_reached

    end

end


text:
module AnalysisFTIR

using MultiLayerNFRHT, Plots
using CSV, DataFrames, LsqFit
using DSP
using LinearAlgebra, Statistics

export load_data,convertto_freq,convertto_wavel,
       unnormalize,mynormalize,plot_data,reflectivity_data,
       emissivity_data,comparison,emissivity_model,
       save_data, fit_data,
       window_average, savitsky_golay

" Fitting options "
function option(opt :: Symbol)
    return options[opt]
end

" Load the data and put into an array "
function load_data(filename)
    df    = DataFrame(CSV.File(filename))
    mdata = convert(Array,df)
    return mdata
end

" Convert cm^{-1} to frequency in rad/s "
convertto_freq(wavenum) = 1e2*2.0*pi*c0*wavenum

" Convert cm^{-1} to wavelength in microns "
convertto_wavel(wavenum) = 1e4/wavenum

" Remove the reference from the experimental data"
unnormalize(data,refdata) = data[:,2].*refdata[:,2]

" Normalize the experimental data"
mynormalize(data,normdata) = data[:,2]./normdata[:,2]

" Plot reflectivity/emissivity as a function of frequency"
function plot_data(data1,data2)
    plot(data1[:,1],data1[:,2]; xscale= :log10,ylim=(0,1))
    plot!(data2[:,1],data2[:,2]; xscale= :log10,ylim=(0,1))
end

"""
Extract reflectivity with a given reference where srcdata, initrefdata and refdata are filenames.
initrefdata is the data that was defined as reference data during the measurements.
All measurements are normalized by this dataset. All source data have to be multiplied by this set.
"""
function reflectivity_data(srcdata,initrefdata,refdata)
    data1      = load_data(srcdata)
    initrefdat = load_data(initrefdata)
    refdat     = load_data(refdata)

    data2      = similar(refdat)
    data2[:,1] = convertto_freq.(data1[:,1])
    data2[:,2] = mynormalize(data1,refdat) # since both datasets are devided by the same signal, unnormalize is not necessary
    return data2
end

" Get the emissivity from the reflectivity measurement "
function emissivity_data(srcdata,initrefdata,refdata)
    data      = reflectivity_data(srcdata,initrefdata,refdata)
    data[:,2].= 1.0 .- data[:,2]
    return data
end

" Compare simulation and measurement "
function comparison(srcdata,initrefdata,refdata,substrate :: MultiLayer; save = false, outputfilename = "")
    emdata   = emissivity_data(srcdata,initrefdata,refdata)
    emw      = copy(emdata)
    emw[:,2] = emissivity_kx_w.([substrate],0.0, emdata[:,1]).*0.5
    if save == true
        save_data(emdata,emw,outputfilename)
    end
    plot_data(emdata,emw)
end

" Fit model to experimental data "
function fit_data(emissivity_model, xdata, ydata, p; save = false, outputfilename = "" ,opt = Dict("parameter 1 (unit)"=> 1))
    fit     = curve_fit(emissivity_model,xdata,ydata,p)
    fit_dat = similar(xdata)
    fit_dat = emissivity_model(xdata, fit.param)
    df_fit = value_of_fit(fit,opt)
    if save == true
        save_data([xdata  ydata] ,[xdata fit_dat], outputfilename)
        CSV.write("fit_parameters_"*outputfilename,df_fit)
    end
    data = [xdata  ydata  fit_dat]
    println(df_fit)
    plot_data([xdata ydata],[xdata fit_dat])
end

function value_of_fit(fit, dict)
    df = DataFrame()
    df[!,:Parameters] = [collect(keys(dict)) ; "Sum residuals"]
    df[!,:Values]     = [[fit.param[i] for i in collect(values(dict))] ; sum(fit.resid.^2)]
    return df
end

" Save experimental and simulated emissivity to a *.dat file"
function save_data(data_exp,data_sim,outputfilename)
    df        = DataFrame()
    df[!,:freq] = data_exp[:,1]
    df[!,:Data] = data_exp[:,2]
    df[!,:Fit]  = data_sim[:,2]
    CSV.write(outputfilename, df)
end

" Sliding window average "
function window_average(data :: AbstractArray, window :: Integer)
    data_av = zeros(length(data))
    for i=1:length(data)
        data_av[i] = mean(data[i:min(i+window-1,end)])
    end
    return data_av
end

" Savitsky-Golay filter for smoothing data "
function savitsky_golay(x::Vector, windowSize::Integer, polyOrder::Integer; deriv::Integer=0)

    #Some error checking
    @assert isodd(windowSize) "Window size must be an odd integer."
    @assert polyOrder < windowSize "Polynomial order must me less than window size."

    halfWindow = Int((windowSize-1)/2)

    #Setup the S matrix of basis vectors.
    S = zeros(windowSize, polyOrder+1)
    for ct = 0:polyOrder
	    S[:,ct+1] = collect(-halfWindow:halfWindow).^(ct)
    end

    #Compute the filter coefficients for all orders
    #From the scipy code it seems pinv(S) and taking rows should be enough
    G = S*pinv(S'*S)

    #Slice out the derivative order we want
    filterCoeffs = G[:,deriv+1] * factorial(deriv)

    #Pad the signal with the endpoints and convolve with filter
    paddedX = [x[1]*ones(halfWindow); x ; x[end]*ones(halfWindow)]
    y = conv(filterCoeffs[end:-1:1], paddedX)

    #Return the valid midsection
    return y[2*halfWindow+1:end-2*halfWindow]

end

end # module


text:
#=##############################################################################################
Filename: _plot.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: after finishing a simulation, a list of ploting functions is executed
=###############################################################################################

"""
Describes a function used for postprocessing in `runsimulation!` with:

Inputs:

* `aircraft::Aircraft` : aircraft system struct
* `parameters <: Parameters` : inherits from the `Parameters` type; object containing data required by `<: Action` functions
* `step_range::AbstractArray` : range of times defining simulation
* `step_symbol::String` : defines the step, e.g. `alpha` or `time`

Modifies:

* `aircraft::Aircraft`
* `parameters <: Parameters`

Outputs:

* `flag::Bool` : true if an action experiences errors

Dispatch is as follows:

```
function mypostaction!(aircraft, parameters, step_range, step_symbol)
    ...
    return flag
end
```
"""
abstract type PostAction end


text:
function createMesh(L::Float64, dx::Float64)
    N = Int(round(L/dx))
    Nc, Ni = 4, (N-1)^2
    Nm, Ns = 2*(N-1), 2*(N-1)
    Nel = N^2
    NDOF = 2*(N+1)*(N+1)
    nodes = zeros(Int64, N+1, N+1)

    # Corners
    nodes[end,1], nodes[end,end] = (1,2)
    nodes[1,end], nodes[1,1] = (3,4)
    # Internal
    nodes[2:N,2:N] = reshape(5 : Ni+4, (N-1, N-1))
    # Left
    nodes[2:end-1,1] = Ni+5 : Ni+N+3
    # Bottom
    nodes[end,2:end-1] = Ni+N+4 : Ni+2*(N+1)
    # Right
    nodes[2:end-1,end] = Ni+2*N+3 : Ni+3*N+1
    # Top
    nodes[1,2:end-1] = Ni+3*N+2 : Ni+4*N

    # Generate connectivity of each element based on their node numbers.
    NC = zeros(Int64, Nel, 4)
    for j in 1:N
        for i in 1:N
            NC[i+(j-1)*N, :] = [nodes[end-j+1,i], nodes[end-j+1,i+1], nodes[end-j,i+1], nodes[end-j,i]]
        end
    end
    # Create square mesh over the domain.
    x, y = (0:dx:L, 0:dx:L)
    XY = zeros(Float64, (N+1)^2, 2)
    for j in 1:N+1
        for i in 1:N+1
            XY[i+(j-1)*(N+1), :] = [x[i], y[j]]
        end
    end

    elCoords = zeros(Float64, Nel, 4, 2)
    for j in 1:N
        for i in 1:N
            k = i+(j-1)*(N+1)
            ind = [k, k+1, k+N+2, k+N+1]
            elCoords[i+(j-1)*N, :, :] = XY[ind,:]
        end
    end
    DOFs = zeros(Int64, Nel, 4, 2)
    for i in 1:Nel
        DOFs[i, :, :] = reshape(hcat(2*NC[i,:]' .- 1, 2*NC[i,:]'), (4,2))
    end

    # Dependency matrix
    rowInd = 1 : 2*Nm
    colInd = vcat(repeat([3,4], N-1), repeat([7,8], N-1))
    data = ones(Int64, 2*Nm)
    D = sparse(rowInd, colInd, data)
    Nci = Nc + Ni
    Nim = Ni + Nm
    Ncim = Nc + Ni + Nm
    Nims = Ni + Nm + Ns
    Ncims = Nc + Nims
    Ns = [Nc, Ni, Nm, Ns, Nel, NDOF, Nci, Nim, Ncim, Nims, Ncims]

   return N, NC, XY, elCoords, DOFs, D, Ns;
end

max_stars_repo_name: likemaestro/Nonlinear-Homogenization

max_issues_repo_name: likemaestro/Nonlinear-Homogenization

max_forks_repo_name: likemaestro/Nonlinear-Homogenization


text:
# BL solution
# plotting the results of water-flooding oil recovery
# optimizing the rel-perm parametrs using the recovery data
using Roots, Dierckx, Plots, JLD
# include("rel_perms.jl")

# plot the fractional flow functions for a water-wet and an oil-wet system
include("rel_perms_real.jl")
# Water-wet Mixed-wet Oil-wet
# Swi 0.2 0.2 0.2
# Sor 0.1 0.2 0.3
# Ew=Krw(Sor) 0.6 0.8 0.8
# Eo=Kro(Swi) 0.8 0.8 0.5
# Nw 3 2 3
# No 2 2 3
krw0_ww = 0.3
krw0_ow = 1.0
kro0_ww = 0.6
kro0_ow = 0.76
nw_ww = 2.4
nw_ow= 2.4
no_ww = 2.0
no_ow= 2.0
sor_ww=0.15
sor_ow=0.2
sor_mw=0.08
swc_ww=0.09
swc_ow=0.09
teta_ow=deg2rad(180)
teta_ww=deg2rad(0)
teta_mw=deg2rad(90)
teta_prime=teta_mw-π/2 # see my report
a=(sor_mw-0.5*(sor_ww-sor_ow)*cos(teta_mw)-0.5*(sor_ww+sor_ow))/(cos(teta_mw)^2-1)
b=0.5*(sor_ww-sor_ow)
c=0.5*(sor_ww+sor_ow)-a

SF=0.0 # 1.0 is water-wet, 0.0 is oil-wet
sw_plot_ow=linspace(swc_ow,1-sor_ow,100)
sw_plot_ww=linspace(swc_ww,1-sor_ww,100)

plot(size=(500,400), xtickfont = font(12, "Courier"), ytickfont=font(12, "Courier"),
  ylabel="krw, kro", xlabel="Water saturation", legendfont=font(12, "Courier"),
  guidefont=font(12, "Courier"))
plot!(sw_plot_ow, krw.(sw_plot_ow, krw0_ow, sor_ow, swc_ow, nw_ow),
  xlims=(0.0,1.0), linewidth=3, linestyle=:dash, linecolor=:red,
  label="krw, Oil-wet")
plot!(sw_plot_ow, kro.(sw_plot_ow, kro0_ow, sor_ow, swc_ow, no_ow),
  xlims=(0.0,1.0), linewidth=3, linestyle=:dash, linecolor=:red,
  label="kro, Oil-wet")
plot!(sw_plot_ww, krw.(sw_plot_ww, krw0_ww, sor_ww, swc_ww, nw_ww),
  xlims=(0.0,1.0), linewidth=1, linestyle=:solid, linecolor=:blue,
  label="krw, Water-wet")
plot!(sw_plot_ww, kro.(sw_plot_ww, kro0_ww, sor_ww, swc_ww, no_ww),
  xlims=(0.0,1.0), linewidth=1, linestyle=:solid, linecolor=:blue,
  label="kro, Water-wet")
savefig("two_rel_perms.png")

# plot the recovery factors for a water-wet, a mixed-wet, and an oil-wet system
include("frac_flow_funcs.jl")
L_res=1000 # m
perm_res=0.001e-12 # Darcy
inj_pv=5.0 # pore volumes
poros_res=0.35 # reservoir porosity
inj_vel=1.0/(3600*24) # m/s
mu_oil=2e-3 # oil viscosity
mu_water=1e-3 # water viscosity

# water-wet reservoir
(xt_shock_ww, sw_shock_ww, xt_prf_ww, sw_prf_ww, t_ww, p_inj_ww, R_oil_ww)=
frac_flow_wf(muw=mu_water, muo=mu_oil, ut=inj_vel, phi=poros_res,
  k=perm_res, swc=swc_ww, sor=sor_ww, kro0=kro0_ww, no=no_ww, krw0=krw0_ww,
  nw=nw_ww, sw0=swc_ww, sw_inj=1.0, L=L_res, pv_inj=inj_pv)
# oil-wet reservoir
(xt_shock_ow, sw_shock_ow, xt_prf_ow, sw_prf_ow, t_ow, p_inj_ow, R_oil_ow)=
frac_flow_wf(muw=mu_water, muo=mu_oil, ut=inj_vel, phi=poros_res,
  k=perm_res, swc=swc_ow, sor=sor_ow, kro0=kro0_ow, no=no_ow, krw0=krw0_ow,
  nw=nw_ow, sw0=swc_ow, sw_inj=1.0, L=L_res, pv_inj=inj_pv)

plot(t_ow/(24*3600), R_oil_ow, size=(400,300), linecolor=:black, linewidth=2,
  xtickfont = font(10, "Courier"), ytickfont=font(10, "Courier"),
  xlabel="time [day]", ylabel="Recovery factor", label="Oil-wet",
  legendfont=font(10, "Courier"), guidefont=font(12, "Courier"), ylims=(0.0,1.0))
# savefig("recovery_time_ow.png")
plot!(t_ww/(24*3600), R_oil_ww, linecolor=:blue, label="Water-wet", linewidth=2,
  linestyle=:dash)
# savefig("recovery_time_ww_ow.png")

# BL for different wettability
plot(size=(500,400), xtickfont = font(10, "Courier"), ytickfont=font(10, "Courier"),
  xlabel="t [day]", ylabel="Recovery factor", legendfont=font(10, "Courier"),
  guidefont=font(12, "Courier"))
for SF in 0.0:0.25:1.0
  krw0=krw0_ww*SF+krw0_ow*(1-SF)
  kro0=kro0_ww*SF+kro0_ow*(1-SF)
  #sor=sor_ww*SF+sor_ow*(1-SF)
  swc=swc_ww*SF+swc_ow*(1-SF)
  no= no_ww*SF+no_ow*(1-SF)
  nw= nw_ww*SF+nw_ow*(1-SF)
  teta=teta_ww*SF+teta_ow*(1-SF)
  # sor=(max(sor_ow, sor_ww)-sor_mw)*cos(teta-teta_prime)^2+sor_mw
  sor=a*cos(teta).^2+b*cos(teta)+c
  (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R_oil)=
  frac_flow_wf(muw=mu_water, muo=mu_oil, ut=inj_vel, phi=poros_res,
    k=perm_res, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0,
    nw=nw, sw0=swc, sw_inj=1.0, L=L_res, pv_inj=inj_pv)
    plot!(t/(24*3600), R_oil, size=(500,400), linewidth=2,
      xtickfont = font(10, "Courier"), ytickfont=font(10, "Courier"),
      xlabel="time [day]", ylabel="Recovery factor", label="$SF Sal. Factor",
      legendfont=font(10, "Courier"), guidefont=font(12, "Courier"), ylims=(0.0,1.0))
    savefig("recovery_forced_imb_$SF-WW.png")
end
plot!()

# plot(t, R_oil)
# plot the imbibition Pc-Sw curve for different SF
plot(size=(500,400), xtickfont = font(10, "Courier"), ytickfont=font(10, "Courier"),
  xlabel="Sw", ylabel="Pc [Pa]", legendfont=font(10, "Courier"),
  guidefont=font(12, "Courier"))
for SF in 0.1:0.2:0.9
  krw0=krw0_ww*SF+krw0_ow*(1-SF)
  kro0=kro0_ww*SF+kro0_ow*(1-SF)
  #sor=sor_ww*SF+sor_ow*(1-SF)
  swc=swc_ww*SF+swc_ow*(1-SF)
  no= no_ww*SF+no_ow*(1-SF)
  nw= nw_ww*SF+nw_ow*(1-SF)
  teta=teta_ww*SF+teta_ow*(1-SF)
  # sor=(max(sor_ow, sor_ww)-sor_mw)*cos(teta-teta_prime)^2+sor_mw
  sor=a*cos(teta).^2+b*cos(teta)+c
  pce= 5e2
  sw_plot=linspace(0.0,1.0,100)

  pc=pc_imb.(sw_plot, pce, swc, sor, teta=teta)
  plot!(sw_plot, pc, xlims=(0.0,1.0), ylims=(-2000, 2000), linewidth=2,
  label="$SF Sal. Factor")
  savefig("pc_sw_$SF-WW.png")
end
plot!()

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SphericalVariogram(range=r, sill=s, nugget=n)
    SphericalVariogram(ball; sill=s, nugget=n)

A spherical variogram with range `r`, sill `s` and nugget `n`.
Optionally, use a custom metric `ball`.
"""
struct SphericalVariogram{V,B} <: Variogram
  sill::V
  nugget::V
  ball::B
end

SphericalVariogram(ball; sill=1.0, nugget=zero(typeof(sill))) =
  SphericalVariogram(sill, nugget, ball)

SphericalVariogram(; range=1.0, sill=1.0, nugget=zero(typeof(sill))) =
  SphericalVariogram(sill, nugget, MetricBall(range))

function (γ::SphericalVariogram)(h::T) where {T}
  r = radius(γ.ball)
  s = γ.sill
  n = γ.nugget

  # constants
  c1 = T(3) / T(2)
  c2 = T(1) / T(2)
  s1 = c1*(h/r) - c2*(h/r)^3
  s2 = T(1)

  (h < r) * (s - n) * s1 +
  (h ≥ r) * (s - n) * s2 +
  (h > 0) * n
end

isstationary(::Type{<:SphericalVariogram}) = true


text:
#=
This file handles fetching data from the StaffJoy API, parsing it, running a calculation, and sending it back to StaffJoy.
=#
import StaffJoy

function get_task()
    results = Requests.post(build_uri(config, "/api/v1/tasking/"))
    if results.status == 404
        return nothing
    elseif results.status != 200
        # Bad request? Internal error? This is a pageable offense!
        code = results.status
        Logging.warn("Tasking API Down - code $code")
        return nothing
    end

    # Return task
    task = JSON.parse(Requests.text(results))
    Logging.info("Task Received - $task")
    # Task Fields: schedule_id, role_id, location_id, organization_id
    return task
end

function mark_task_done(task)
    path = string("/api/v1/tasking/", task["schedule_id"])
    results = Requests.delete(build_uri(config, path))
    if results.status != 204  # 204 is "success / empty response"
        # Bad request? Internal error? This is a pageable offense!
        code = results.status
        schedule = task["schedule_id"]
        Logging.err( "Failed to mark schedule $schedule as done. Code $code")
        error()
    end

    Logging.info( "Task Completed - $task")
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
##### Generic methods #####

"""
    size(d::MatrixDistribution)

Return the size of each sample from distribution `d`.
"""
size(d::MatrixDistribution)

"""
    length(d::MatrixDistribution)

The length (*i.e* number of elements) of each sample from the distribution `d`.
"""
Base.length(d::MatrixDistribution)

"""
    mean(d::MatrixDistribution)

Return the mean matrix of `d`.
"""
mean(d::MatrixDistribution)

# sampling

rand!(d::MatrixDistribution, A::AbstractArray{M}) where {M<:Matrix} = _rand!(sampler(d), A)

"""
    rand(d::MatrixDistribution, n)

Draw a sample matrix from the distribution `d`.
"""
rand(d::MatrixDistribution, n::Int) = _rand!(sampler(d), Vector{Matrix{eltype(d)}}(n))

# pdf & logpdf

_pdf(d::MatrixDistribution, x::AbstractMatrix{T}) where {T<:Real} = exp(_logpdf(d, x))

"""
    logpdf(d::MatrixDistribution, AbstractMatrix)

Compute the logarithm of the probability density at the input matrix `x`.
"""
function logpdf(d::MatrixDistribution, x::AbstractMatrix{T}) where T<:Real
    size(x) == size(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _logpdf(d, x)
end

"""
    pdf(d::MatrixDistribution, x::AbstractArray)

Compute the probability density at the input matrix `x`.
"""
function pdf(d::MatrixDistribution, x::AbstractMatrix{T}) where T<:Real
    size(x) == size(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _pdf(d, x)
end

function _logpdf!(r::AbstractArray, d::MatrixDistribution, X::AbstractArray{M}) where M<:Matrix
    for i = 1:length(X)
        r[i] = logpdf(d, X[i])
    end
    return r
end

function _pdf!(r::AbstractArray, d::MatrixDistribution, X::AbstractArray{M}) where M<:Matrix
    for i = 1:length(X)
        r[i] = pdf(d, X[i])
    end
    return r
end

function logpdf!(r::AbstractArray, d::MatrixDistribution, X::AbstractArray{M}) where M<:Matrix
    length(X) == length(r) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _logpdf!(r, d, X)
end

function pdf!(r::AbstractArray, d::MatrixDistribution, X::AbstractArray{M}) where M<:Matrix
    length(X) == length(r) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _pdf!(r, d, X)
end

function logpdf(d::MatrixDistribution, X::AbstractArray{M}) where M<:Matrix
    T = promote_type(partype(d), eltype(M))
    _logpdf!(Array{T}(size(X)), d, X)
end

function pdf(d::MatrixDistribution, X::AbstractArray{M}) where M<:Matrix
    T = promote_type(partype(d), eltype(M))
    _pdf!(Array{T}(size(X)), d, X)
end

"""
    _logpdf(d::MatrixDistribution, x::AbstractArray)

Evaluate logarithm of pdf value for a given sample `x`. This function need not perform dimension checking.
"""
_logpdf(d::MatrixDistribution, x::AbstractArray)

##### Specific distributions #####

for fname in ["wishart.jl", "inversewishart.jl"]
    include(joinpath("matrix", fname))
end


text:
abstract type AbstractVsStrategy end

"""Abstract type for sources (e.g., a game)."""
abstract type AbstractVsSource end

"""Abstract type for scenes."""
abstract type AbstractVsScene end

"""Abstract type for streams. Should support the common functions."""
abstract type AbstractVsStream end

Base.@kwdef struct VsStream <: AbstractVsStream
    video::VideoIO.VideoReader
end
@forward VsStream.video Base.close, Base.eof, Base.seek, Base.seekstart, Base.seekend,
    VideoIO.skipframe, VideoIO.skipframes, VideoIO.framerate

"""Frame."""
Base.@kwdef struct VsFrame{T, TM <: AbstractMatrix{T}} <: AbstractMatrix{T}
    image::TM
    time::Rational{Int64} = 0 // 1
end
image(frame::VsFrame) = frame.image
Base.time(frame::VsFrame) = frame.time
@forward VsFrame.image Base.length, Base.size, Base.getindex, Base.setindex!

Base.@kwdef struct VsContextData
    dict::Dict{Symbol, Any} = Dict()
end
Base.getproperty(data::VsContextData, key::Symbol) = get(getfield(data, :dict), key, nothing)
Base.setproperty!(data::VsContextData, key::Symbol, value) = getfield(data, :dict)[key] = value
Base.propertynames(data::VsContextData) = keys(getfield(data, :dict))

Base.@kwdef struct VsConfig{
    TStrategy <: AbstractVsStrategy,
    TSource <: AbstractVsSource,
}
    # in seconds
    process_interval::Float64 = 1.0
    num_skip_frames::Int = 59
    use_gray_image::Bool = true
    ocr_language::String = "eng"
    strategy::TStrategy
    source::TSource
end

"""Context."""
Base.@kwdef mutable struct VsContext{
    TStrategy <: AbstractVsStrategy,
    TSource <: AbstractVsSource
}
    config::VsConfig{TStrategy, TSource}
    stream::Nullable{AbstractVsStream} = nothing
    # Language => TessInst
    ocr_instances::OrderedDict{String, TessInst}
    current_frame::Nullable{VsFrame} = nothing
    current_scene::Nullable{AbstractVsScene} = nothing
    current_time_skippable::Float64 = 0.0
    data::VsContextData = VsContextData()
end

"""
    vs_setup(source_type; kwargs...)

Create a `VsConfig` with given args.
"""
vs_setup(::Type{T}; kwargs...) where T <: AbstractVsSource =
    error("Not implemented: vs_setup($T)")

"""
    vs_init!(context)

Initialize.
"""
vs_init!(::T) where T <: VsContext =
    error("Not implemented: vs_init!($T)")

"""
    vs_parse_frame!(context, frame)

Parse a frame. Returns a `Union{AbstractVsScene, Nothing}`.
"""
vs_parse_frame!(::T, ::VsFrame) where T <: VsContext =
    error("Not implemented: vs_parse_frame!($T, VsFrame)")

"""
    vs_update!(context, scene)

Store `scene` in `context`.
"""
vs_update!(::T, ::TS) where {T <: VsContext, TS <: AbstractVsScene} =
    error("Not implemented: vs_parse_frame!($T, $TS)")

"""
    vs_result(context)

Fetch the result.
"""
vs_result(::T) where T <: VsContext = error("Not implemented: vs_result($T)")

vs_tryparse_scene(
    ::Type{<:AbstractVsScene},
    ::VsFrame,
    ::Union{VsContext, Nothing} = nothing;
    force::Bool = false
) = nothing


text:
# load required packages and codes -
include("Include.jl")

# set the path to the test model file -
path_to_model_file = joinpath(_PATH_TO_MODEL, "Test.vff")

# build the default model -
model = build_default_model(path_to_model_file)

# the order of the species will be in the species_symbol_array -> E, E_S, P and S
# update the default initial conditions (default: zeros)
#all xₒ units are nM
xₒ = model["initial_conditions_array"]
#xₒ[1] = 3.4e3
xₒ[2] = 1.4e3
xₒ[8] = 9e1
xₒ[13] = 1e-3
#xₒ[16] = 1.0
xₒ[14] = 2.5
xₒ[23] = 2e1
xₒ[24] = 1e1
xₒ[25] = 0.7
xₒ[29] = 1e-1
xₒ[31] = 1.6e2

# set up rate constants - order of rates in the reaction_symbol_array -
#κ have been converted to nM-1s-1 where there is a second order reaction
κ = model["κ"]
κ[1] = 3.1e-3
κ[2] = 3.1e-3
κ[3] = 2.5e-2
κ[4] = 2e-3
κ[5] = 9e-2 #Anand* - indicates where i got the value from, * indicates that i'm not convinced with the values because it wasnt the exact same reaction
κ[6] = 5e-3 #Anand*
κ[7] = 0.5 #Anand*
κ[8] = 0.1 #Varner et al.
κ[9] = 1. #Varner et al.
κ[10] = 0.5 #Varner et al.
κ[11] = 0.1 #Varner et al.
κ[12] = 10. #Varner et al.
κ[13] = 0.5 #Varner et al.
κ[14] = 0.1 #Varner et al.
κ[15] = 5.5 #Varner et al.
κ[16] = 19. #Varner et al.
κ[17] = 2.2e-2
κ[18] = 19
κ[19] = 0.01
κ[20] = 2.4
κ[21] = 1.8
κ[22] = 12e-6 #Varner et al.
κ[23] = 1e-9 #Varner et al.
κ[24] = 30e-4 #Varner et al.
κ[25] = 0.1 #Varner et al.
κ[26] = 15. #Varner et al.
κ[27] = 0.9 #Varner et al.
κ[28] = 1e7*1e-9
κ[29] = 5e-3
κ[30] = 1e8*1e-9
κ[31] = 1e-3
κ[32] = 8.2
κ[33] = 6e-3
κ[34] = 2.2e-5
κ[35] = 1e-3
κ[36] = 1e-3
κ[37] = 0.1 #Varner et al.
κ[38] = 7.2 #Varner et al.
κ[39] = 4.2 #Varner et al.
κ[40] = 4e-1
κ[41] = 0.2
κ[42] = 1e-1
κ[43] = 103
κ[44] = 75.
κ[45] = 0.01 #Varner et al.*
κ[46] = 1.1 #Varner et al.*
κ[47] = 0.15 #Varner et al.*
κ[48] = 0.3
κ[49] = 3.6e-2
κ[50] = 0.32
κ[51] = 1.1e-3
κ[52] = 5e-2
κ[53] = 10e-6
κ[54] = 7.1e-6
κ[55] = 4.9e-7
κ[56] = 7e-6
κ[57] = 3e-7

# solve the model -
tspan = (0.0,1200.0)
Δt = 0.1
prob = ODEProblem(balances, xₒ, tspan, model; saveat = Δt)
soln = solve(prob)

# get the results from the solver -
T = soln.t
Xₘ = hcat(soln.u...)
Thrombin = Xₘ[3,:]+1.2*Xₘ[39,:] #FIIa is the 3rd species in species_symbol_array
#p = CSV.read("data\\15pmTF_TFPIAT.csv",DataFrame)
#plot(T[1:7001],Thrombin[1:7001])


text:
using DICOMClient

using Documenter
using FileIO
using ImageMagick
using Images
using JSON3
using HTTP
using ReferenceTests
using Test

DICOMCLIENT_AUTHENTICATED_TESTS = get(ENV, "DICOMCLIENT_AUTHENTICATED_TESTS", "") == "true"
@info("DICOMCLIENT_AUTHENTICATED_TESTS: $(DICOMCLIENT_AUTHENTICATED_TESTS)")

include("test-utilities.jl")

@testset "DICOMCLIENT.jl" begin
    @testset "Unit tests" begin
        include("unit-tests.jl")
    end
    @testset "Doctests" begin
        doctest(DICOMClient)
    end
    @testset "Integration tests" begin
        @testset "Integration tests: unauthenticated" begin
            include("integration-tests/unauthenticated-integration-tests/basic-test.jl")
        end
        if DICOMCLIENT_AUTHENTICATED_TESTS
            @info("DICOMCLIENT_AUTHENTICATED_TESTS is true. Running the authenticated tests now...")
            @testset "Integration tests: authenticated" begin
            end
        else
            @info("DICOMCLIENT_AUTHENTICATED_TESTS is not true, so we will not run the authenticated tests.")
        end
    end
end


text:
"""
    TwilioSMS

A small package that provides a function `sendsms` for sending text messages using Twilio.
"""
module TwilioSMS

using HTTP
using JSON

export sendsms

"""
    sendsms(msg; from, to, sid, auth_token)

Send a text message using Twilio.

## Arguments
- `msg` The message to send
- `from` The phone number from which to send
- `to` The recipient's phone number
- `sid` The Account SID of your Twilio account. If not specified try to read the environment variable `TWILIO_ACCOUNT_SID`
- `auth_token` The Auth Token of your Twilio account. If not specified try to read the environment variable `TWILIO_AUTH_TOKEN`

## Examples
```julia-repl
julia> sendsms("Hello, World 🌍!", from="+1-202-555-0107", to="+1-512-555-0180")

julia> sendsms("With explicit sid and access token", from="+1-202-555-0107", to="+1-512-555-0180",
                sid="d5d7b014b05d54fa52647820a035089e", auth_token="bed7fe1af0d18dfd-100972663d58f196")
```
"""
function sendsms(msg; from, to, sid=_from_env("TWILIO_ACCOUNT_SID"), auth_token=_from_env("TWILIO_AUTH_TOKEN"))

    payload = HTTP.URIs.escapeuri(Dict(:To => to, :From => from, :Body => msg))

    url = "api.twilio.com/2010-04-01/Accounts/$sid/Messages.json"
    full_url = "https://$sid:$auth_token@$url"

    response = HTTP.post(full_url, ["Content-Type" => "application/x-www-form-urlencoded"], payload, basic_authorization=true)
    response.status == 201 || error("Something went wrong. Respone was:\n$(string(response))")

    return nothing
end

function _from_env(key)

    !haskey(ENV, key) && error("Environment variable \"$key\" could not be found, must be manually specified.")
    return ENV[key]
end

end # module


text:
using GPUifyLoops: @launch, @loop
using Oceananigans.Architectures: device
using Oceananigans.Utils: @loop_xyz

"""
    set!(model; kwargs...)

Set velocity and tracer fields of `model`. The keyword arguments
`kwargs...` take the form `name=data`, where `name` refers to one of the
fields of `model.velocities` or `model.tracers`, and the `data` may be an array,
a function with arguments `(x, y, z)`, or any data type for which a
`set!(ϕ::AbstractField, data)` function exists.

Example
=======
```julia
model = Model(grid=RegularCartesianGrid(size=(32, 32, 32), length=(1, 1, 1))

# Set u to a parabolic function of z, v to random numbers damped
# at top and bottom, and T to some silly array of half zeros,
# half random numbers.

u₀(x, y, z) = z/model.grid.Lz * (1 + z/model.grid.Lz)
v₀(x, y, z) = 1e-3 * rand() * u₀(x, y, z)

T₀ = rand(size(model.grid)...)
T₀[T₀ .< 0.5] .= 0

set!(model, u=u₀, v=v₀, T=T₀)
```
"""
function set!(model; kwargs...)
    for (fldname, value) in kwargs
        if fldname ∈ propertynames(model.velocities)
            ϕ = getproperty(model.velocities, fldname)
        elseif fldname ∈ propertynames(model.tracers)
            ϕ = getproperty(model.tracers, fldname)
        else
            throw(ArgumentError("name $fldname not found in model.velocities or model.tracers."))
        end
        set!(ϕ, value)
    end
    return nothing
end

function set!(Φ::NamedTuple; kwargs...)
    for (fldname, value) in kwargs
        ϕ = getproperty(Φ, fldname)
        set!(ϕ, value)
    end
    return nothing
end

set!(u::Field, v::Number) = @. u.data = v

set!(u::Field{X, Y, Z, A}, v::Field{X, Y, Z, A}) where {X, Y, Z, A} =
    @. u.data.parent = v.data.parent

"Set the CPU field `u` to the array `v`."
function set!(u::AbstractCPUField, v::Array)
    for k in 1:u.grid.Nz, j in 1:u.grid.Ny, i in 1:u.grid.Nx
        u[i, j, k] = v[i, j, k]
    end
    return nothing
end

# Set the GPU field `u` to the array `v`.
@hascuda function set!(u::AbstractGPUField, v::Array)
    v_field = Field(location(u), CPU(), u.grid)
    set!(v_field, v)
    set!(u, v_field)
    return nothing
end

# Set the GPU field `u` to the CuArray `v`.
@hascuda function set!(u::AbstractGPUField, v::CuArray)
    @launch device(GPU()) config=launch_config(u.grid, :xyz) _set_gpu!(u.data, v, u.grid)
    return nothing
end

function _set_gpu!(u, v, grid)
	@loop_xyz i j k grid begin
        u[i, j, k] = v[i, j, k]
    end
    return nothing
end

# Set the GPU field `u` data to the CPU field data of `v`.
@hascuda set!(u::AbstractGPUField, v::AbstractCPUField) = copyto!(u.data.parent, v.data.parent)

# Set the CPU field `u` data to the GPU field data of `v`.
@hascuda set!(u::AbstractCPUField, v::AbstractGPUField) = u.data.parent .= Array(v.data.parent)

"Set the CPU field `u` data to the function `f(x, y, z)`."
set!(u::AbstractCPUField, f::Function) = interior(u) .= f.(nodes(u)...)

# Set the GPU field `u` data to the function `f(x, y, z)`.
@hascuda function set!(u::AbstractGPUField, f::Function)
    u_cpu = Field(location(u), CPU(), u.grid)
    set!(u_cpu, f)
    set!(u, u_cpu)
    return nothing
end