text:
# TODO: SALTbootstrap
# TODO: convenience Simulation modifiers to take all bc's to outgoing.
"""
    module Lasing

For solution of the nonlinear SALT equation.
Exports `HelmholtzSALT` and `MaxwellSALT` constructor and extends methods of
`NLsolve`, which must be separately imported.
"""
module Lasing

export HelmholtzSALT
export MaxwellSALT

using ..Common
using ..Spectral
using Formatting
using LinearAlgebra
using NLsolve
using RecipesBase
using SparseArrays

import ..Common.INDEX_OFFSET

################################################################################
# SALT object definitions and constructors

# works for Helmholtz and Maxwell, field `operator` is accessed through
# properties `helmholtz` or `maxwell`
struct SALTProblem{NMODES,M,TO,DIM}
    operator::TO
    ωs::Vector{Float64}
    ψs::VectorField{DIM,M}
    x::Vector{Float64}
    n::Int
    m::Int
    res::Vector{ComplexF64}
    index::Int
    solved::Ref{Bool}
    converged::Ref{Bool}
end

# main SALTProblem constructor
function SALTProblem{NMODES,MCOMPONENTS}(sim::Simulation{DIM}) where {DIM,NMODES,MCOMPONENTS}
    M = Helmholtz(sim; m=NMODES)
    n = length(sim)
    ωs = Vector{Float64}(undef, NMODES)
    ψs = VectorField{MCOMPONENTS}(sim.x, NMODES)
    n÷2>INDEX_OFFSET || throw("INDEX_OFFSET=$(INDEX_OFFSET÷MCOMPONENTS) requires number of sites in ψs_init ($n) to be greater than $(2INDEX_OFFSET÷MCOMPONENTS)")
    x = Vector{Float64}(undef,2n*NMODES)
    res = Vector{ComplexF64}(undef,n)
    index = n÷2+INDEX_OFFSET
    return SALTProblem{NMODES,MCOMPONENTS,typeof(M),DIM}(M,ωs,ψs,x,n,NMODES,res,index,Ref(false),Ref(false))
end

HelmholtzSALT{NMODES} = SALTProblem{NMODES,1}
"""
    HelmholtzSALT(sim::Simulation, m::Integer) -> salt

SALT object with `m` fields for Helmholtz problem
"""
HelmholtzSALT(sim::Simulation,m::Integer) = HelmholtzSALT{m}(sim)

MaxwellSALT{NMODES} = SALTProblem{NMODES,3}
"""
    MaxwellSALT(sim::Simulation, m::Integer) -> salt

SALT object with `m` fields for Maxwell problem
"""
MaxwellSALT(sim::Simulation,m::Integer) = MaxwellSALT{m}(sim)

function Base.getproperty(s::SALTProblem{N,M}, sym::Symbol) where {N,M}
    if Base.sym_in(sym,(:ω,:freq,:Freq,:frequencies,:frequency,:Frequencies,:Frequency,:omega,:omegas))
        return getfield(s,:ωs)
    elseif Base.sym_in(sym,(:ψ,:field,:Field,:fields,:Fields,:psi,:psis))
        return getfield(s,:ψs)
    elseif Base.sym_in(sym,(:sol,:Sol,:solution,:Solution))
        return (getfield(s,:ωs),getfield(s,:ψs))
    elseif Base.sym_in(sym,(:hemlholtz,:maxwell))
        return getfield(s,:operator)
    elseif sym == :simulation
        return getfield(getfield(s,:operator),:simulation)
    elseif Base.sym_in(sym,(:helmholtz,:maxwell))
        return getfield(s,:operator)
    else
        return getfield(s,sym)
    end
end

function Base.propertynames(::SALTProblem,private=false)
    if private
        return fieldnames(SALTProblem)
    else
        return (:ωs,:ψs,:solution,:simulation)
    end
end

################################################################################
# NLsolve extension

fnames = (:nlsolve,:fixedpoint)

for fn ∈ fnames

    @eval NLsolve.$(fn)(ms::SALTProblem; kwargs...) = $(fn)(ms,ms.ωs,ms.ψs; kwargs...)

    @eval NLsolve.$(fn)(ms::SALTProblem, ωs_init::Real, ψs_init::VectorField; kwargs...) = $(fn)(ms,[ωs_init],ψs_init; kwargs...)

    @eval begin function NLsolve.$(fn)(ms::SALTProblem{N,M}, ωs_init::Vector, ψs_init::VectorField; kwargs...) where {N,M}
            n,m = size(ψs_init)
            N==m || throw("number of modes in ψs_init ($m) must be the same as given in Maxwell_SALT ($N)")
            n==M*length(ms.simulation) || throw("number of sites in ψs_init ($(n÷M)) must be the same as given in simulation ($(length(ms.simulation)))")
            n÷2>INDEX_OFFSET || throw("INDEX_OFFSET=$(INDEX_OFFSET) requires number of sites in ψs_init ($(n÷M)) to be greater than $(2INDEX_OFFSET÷M)")
            ωψ_to_x!(ms.x,ωs_init,ψs_init)
            x_to_ωψ!(ms)
            F0 = similar(ms.x)
            J0 = initialize_J(ms)
            df = OnceDifferentiable(ms,ms,ms,ms.x,F0,J0)
            results = $(fn)(df,ms.x; kwargs...)
            # results = $(fn)(ms,ms.x; kwargs...)
            x_to_ωψ!(ms,results.zero)
            ms.solved[] = true
            ms.converged[] = results.f_converged || results.x_converged
            return results
        end
    end
end

"""
    nlsolve(salt; kwargs...) -> results

    nlsolve(salt, init_ωs, init_ψs; kwargs...) -> results

solve nonlinear SALT equation where `salt` is a `HelmholtzSALT` or `MaxwellSALT`.
Results also stored `salt`.

`init_ψs` must be a `ScalarField` (for Helmholtz) or `ElectricField` (for Maxwell)
"""
nlsolve

"""
    fixedpoint(salt; kwargs...) -> results

    fixedpoint(salt, init_ωs, init_ψs; kwargs...) -> results

Convenience wrapper for `NLsolve`'s `fixedpoint` wrapper.
Solves nonlinear SALT equation where `scpa` is a `HelmholtzSALT` or `MaxwellSALT`.
Results also stored `scpa`.

`init_ψs` must be a `ScalarField` (for Helmholtz) or `ElectricField` (for Maxwell)
"""
fixedpoint

################################################################################
# General nonlinear objective + jacobian
(s::SALTProblem)(F::Vector,x::Vector) = s(F,nothing,x)
(s::SALTProblem)(J::AbstractMatrix,x::Vector) = s(nothing,J,x)

################################################################################
# Helmholtz nonlinear objective + jacobian

# Helmholtz Jacobians/Objective
@inline function (ms::HelmholtzSALT)(F,J,x::Vector)
    x_to_ωψ!(ms,x)
    n,m = ms.n,ms.m
    nm = n*m
    index = ms.index

    if !isnothing(J) ψ = ms.ψs end

    ψμ = ms.ψs[:,1]
    @inbounds for μ ∈ 1:m
        Aμ = ms.helmholtz(ms.ωs[μ],ms.ωs,ms.ψs)
        @inbounds for i ∈ 1:n ψμ[i] = ms.ψs[i,μ] end

        if !isnothing(F)
            mul!(ms.res,Aμ,ψμ)
            rdiv!(ms.res,ψμ[index])
            @inbounds @simd for i ∈ 1:n F[(μ-1)n+i], F[nm+(μ-1)n+i] = reim(ms.res[i]) end
        end

        if !isnothing(J)
            fill!(J,0)
            ωμϕ = ms.ωs[μ]/ψμ[index]
            ω²μϕ = ms.ωs[μ]^2/ψμ[index]
            @inbounds for ν ∈ 1:m # derivative wrt to ν, evaluated for field μ
                ρνμω²μ = ω²μϕ*ms.ψs[index,ν]
                @fastmath @inbounds @simd for j ∈ 1:n
                    rr,ir = reim(ρνμω²μ*ms.helmholtz.dFχdψr[j,μ,ν]*ψ[j,μ])
                    J[j+(μ-1)n   , j+(ν-1)n] = rr
                    J[j+(μ-1)n+nm, j+(ν-1)n] = ir

                    rr,ir = reim(ρνμω²μ*ms.helmholtz.dFχdψi[j,μ,ν]*ψ[j,μ])
                    J[j+(μ-1)n   , j+(ν-1)n+nm] = rr
                    J[j+(μ-1)n+nm, j+(ν-1)n+nm] = ir
                end

                @fastmath @inbounds @simd for i ∈ 1:n
                    J[(μ-1)n+i,(ν-1)n+index], J[nm+(μ-1)n+i,(ν-1)n+index] = reim(ω²μϕ*(ms.helmholtz.dFχdϕ[i,μ,ν]*ψμ[i]))
                    J[(μ-1)n+i,nm+(ν-1)n+index], J[nm+(μ-1)n+i,nm+(ν-1)n+index] = reim(ω²μϕ*ms.helmholtz.dFχdω[i,μ,ν]*ψμ[i])
                end
                if μ == ν
                    δω = 1e-9
                    ∂∇²ψ∂ω = reduce(+,map((f,σ) -> σ*(f(ms.ωs[μ]+δω)-f(ms.ωs[μ]-δω))*ψμ./(2δω*ψμ[index]),ms.simulation.self_energy.f,ms.helmholtz.Σs[1:end-1]))
                    @fastmath @inbounds for i ∈ 1:n
                        tr, ti = reim(2ωμϕ*(ms.helmholtz.αεpFχ.nzval[i]*ψμ[i]))
                        J[     (μ-1)n+i,nm+(ν-1)n+index] += tr + real(∂∇²ψ∂ω[i])
                        J[nm + (μ-1)n+i,nm+(ν-1)n+index] += ti + imag(∂∇²ψ∂ω[i])
                    end
                end
            end

            rows = rowvals(Aμ)
            vals = nonzeros(Aμ)
            @inbounds for col ∈ 1:n
                if col!==index
                    @inbounds for j ∈ nzrange(Aμ, col)
                        row = rows[j]
                        Ar, Ai = reim(Aμ.nzval[j])
                        J[(μ-1)n+row   , (μ-1)n+col  ] += Ar
                        J[(μ-1)n+nm+row, (μ-1)n+col  ] += Ai
                        J[(μ-1)n+row   , (μ-1)n+nm+col] -= Ai
                        J[(μ-1)n+nm+row, (μ-1)n+nm+col] += Ar
                    end
                end
            end
        end
    end
    return nothing
end

# Helmholtz initialize Jacobian
function initialize_J(ms::HelmholtzSALT)
    n,m = ms.n,ms.m
    nm = n*m
    index = ms.index

    max_count = (16m^2+6m+27m)n
    count = 1
    rows = Vector{Float64}(undef,max_count)
    cols = Vector{Float64}(undef,max_count)
    vals = zeros(Float64,max_count)
    for μ ∈ 1:m
        Aμ = ms.helmholtz(ms.ωs[μ],ms.ωs,ms.ψs)
        for ν ∈ 1:m # derivative wrt to ν, evaluated for field μ
            for j ∈ 1:n
                k = j

                rows[count] = k+(μ-1)n
                cols[count] = j+(ν-1)n          ; count +=1
                rows[count] = k+(μ-1)n+nm
                cols[count] = j+(ν-1)n          ; count +=1

                rows[count] = k+(μ-1)n
                cols[count] = j+(ν-1)n+nm       ; count +=1
                rows[count] = k+(μ-1)n+nm
                cols[count] = j+(ν-1)n+nm       ; count +=1
            end

            for i ∈ 1:n
                rows[count] = (μ-1)n+i
                cols[count] = (ν-1)n+index      ; count +=1

                rows[count] = nm+(μ-1)n+i
                cols[count] = (ν-1)n+index      ; count +=1

                rows[count] = (μ-1)n+i
                cols[count] = nm+(ν-1)n+index   ; count +=1

                rows[count] = nm+(μ-1)n+i
                cols[count] = nm+(ν-1)n+index   ; count +=1
            end
            if μ == ν
                for i ∈ 1:n
                    rows[count] = (μ-1)n+i
                    cols[count] = nm+(ν-1)n+index ; count +=1

                    rows[count] = nm+(μ-1)n+i
                    cols[count] = nm+(ν-1)n+index ; count +=1
                end
            end
        end

        rowsA = rowvals(Aμ)
        for col ∈ 1:n
            if col!==index
                for j ∈ nzrange(Aμ, col)
                    row = rowsA[j]

                    rows[count] = (μ-1)n+row
                    cols[count] = (μ-1)n+col    ; count +=1

                    rows[count] = (μ-1)n+nm+row
                    cols[count] = (μ-1)n+col    ; count +=1

                    rows[count] = (μ-1)n+row
                    cols[count] = (μ-1)n+nm+col ; count +=1

                    rows[count] = (μ-1)n+nm+row
                    cols[count] = (μ-1)n+nm+col ; count +=1
                end
            end
        end
    end
    c = count - 1
    return sparse(rows[1:c],cols[1:c],vals[1:c],2nm,2nm)
end

################################################################################
# Maxwell nonlinear objective + jacobian

# Maxwell
@inline function (ms::MaxwellSALT)(F,J,x::Vector)
    x_to_ωψ!(ms,x)
    n,m = ms.n,ms.m
    nm = n*m
    index = ms.index

    if !isnothing(J) ψx, ψy, ψz = ms.ψs.x, ms.ψs.y, ms.ψs.z end

    ψμ = ms.ψs[:,1]
    @inbounds for μ ∈ 1:m
        Aμ = ms.maxwell(ms.ωs[μ],ms.ωs,ms.ψs)
        @inbounds for i ∈ 1:n ψμ[i] = ms.ψs[i,μ] end

        if !isnothing(F)
            mul!(ms.res,Aμ,ψμ)
            rdiv!(ms.res,ψμ[index])
            @inbounds @simd for i ∈ 1:n F[(μ-1)n+i], F[nm+(μ-1)n+i] = reim(ms.res[i]) end
        end

        if !isnothing(J)
            fill!(J,0)
            ωμϕ = ms.ωs[μ]/ψμ[index]
            ω²μϕ = ms.ωs[μ]^2/ψμ[index]
            @inbounds for ν ∈ 1:m # derivative wrt to ν, evaluated for field μ
                ρνμω²μ = ω²μϕ*ms.ψs[index,ν]
                @fastmath @inbounds @simd for j ∈ 1:n
                    k = mod1(j,n÷3)

                    rr,ir = reim(-ρνμω²μ*ms.maxwell.dFχdψr[j,μ,ν]*ψx[k,μ])
                    J[k+0n÷3+(μ-1)n   , j+(ν-1)n] = rr
                    J[k+0n÷3+(μ-1)n+nm, j+(ν-1)n] = ir

                    rr,ir = reim(-ρνμω²μ*ms.maxwell.dFχdψr[j,μ,ν]*ψy[k,μ])
                    J[k+1n÷3+(μ-1)n   , j+(ν-1)n] = rr
                    J[k+1n÷3+(μ-1)n+nm, j+(ν-1)n] = ir

                    rr,ir = reim(-ρνμω²μ*ms.maxwell.dFχdψr[j,μ,ν]*ψz[k,μ])
                    J[k+2n÷3+(μ-1)n   , j+(ν-1)n] = rr
                    J[k+2n÷3+(μ-1)n+nm, j+(ν-1)n] = ir

                    rr,ir = reim(-ρνμω²μ*ms.maxwell.dFχdψi[j,μ,ν]*ψx[k,μ])
                    J[k+0n÷3+(μ-1)n   , j+(ν-1)n+nm] = rr
                    J[k+0n÷3+(μ-1)n+nm, j+(ν-1)n+nm] = ir

                    rr,ir = reim(-ρνμω²μ*ms.maxwell.dFχdψi[j,μ,ν]*ψy[k,μ])
                    J[k+1n÷3+(μ-1)n   , j+(ν-1)n+nm] = rr
                    J[k+1n÷3+(μ-1)n+nm, j+(ν-1)n+nm] = ir

                    rr,ir = reim(-ρνμω²μ*ms.maxwell.dFχdψi[j,μ,ν]*ψz[k,μ])
                    J[k+2n÷3+(μ-1)n   , j+(ν-1)n+nm] = rr
                    J[k+2n÷3+(μ-1)n+nm, j+(ν-1)n+nm] = ir
                end

                @fastmath @inbounds @simd for i ∈ 1:n
                    J[(μ-1)n+i,(ν-1)n+index], J[nm+(μ-1)n+i,(ν-1)n+index] = reim(-ω²μϕ*(ms.maxwell.dFχdϕ[i,μ,ν]*ψμ[i]))
                    J[(μ-1)n+i,nm+(ν-1)n+index], J[nm+(μ-1)n+i,nm+(ν-1)n+index] = reim(-ω²μϕ*ms.maxwell.dFχdω[i,μ,ν]*ψμ[i])
                end
                if μ == ν
                    @fastmath @inbounds for i ∈ 1:n
                        tr, ti = reim(-2ωμϕ*(ms.maxwell.αεpFχ.nzval[i]*ψμ[i]))
                        J[     (μ-1)n+i,nm+(ν-1)n+index] += tr
                        J[nm + (μ-1)n+i,nm+(ν-1)n+index] += ti
                    end
                end
            end

            rows = rowvals(Aμ)
            vals = nonzeros(Aμ)
            @inbounds for col ∈ 1:n
                if col!==index
                    @inbounds for j ∈ nzrange(Aμ, col)
                        row = rows[j]
                        Ar, Ai = reim(Aμ.nzval[j])
                        J[(μ-1)n+row   , (μ-1)n+col  ] += Ar
                        J[(μ-1)n+nm+row, (μ-1)n+col  ] += Ai
                        J[(μ-1)n+row   , (μ-1)n+nm+col] -= Ai
                        J[(μ-1)n+nm+row, (μ-1)n+nm+col] += Ar
                    end
                end
            end
        end
    end
    return nothing
end

# Maxwell initialize Jacobian
function initialize_J(ms::MaxwellSALT)
    n,m = ms.n,ms.m
    nm = n*m
    index = ms.index

    max_count = (16m^2+6m+27m)n
    count = 1
    rows = Vector{Float64}(undef,max_count)
    cols = Vector{Float64}(undef,max_count)
    vals = zeros(Float64,max_count)
    for μ ∈ 1:m
        Aμ = ms.maxwell(ms.ωs[μ],ms.ωs,ms.ψs)
        for ν ∈ 1:m # derivative wrt to ν, evaluated for field μ
            for j ∈ 1:n
                k = mod1(j,n÷3)

                rows[count] = k+0n÷3+(μ-1)n
                cols[count] = j+(ν-1)n          ; count +=1
                rows[count] = k+0n÷3+(μ-1)n+nm
                cols[count] = j+(ν-1)n          ; count +=1

                rows[count] = k+1n÷3+(μ-1)n
                cols[count] = j+(ν-1)n          ; count +=1
                rows[count] = k+1n÷3+(μ-1)n+nm
                cols[count] = j+(ν-1)n          ; count +=1

                rows[count] = k+2n÷3+(μ-1)n
                cols[count] = j+(ν-1)n          ; count +=1
                rows[count] = k+2n÷3+(μ-1)n+nm
                cols[count] = j+(ν-1)n          ; count +=1

                rows[count] = k+0n÷3+(μ-1)n
                cols[count] = j+(ν-1)n+nm       ; count +=1
                rows[count] = k+0n÷3+(μ-1)n+nm
                cols[count] = j+(ν-1)n+nm       ; count +=1

                rows[count] = k+1n÷3+(μ-1)n
                cols[count] = j+(ν-1)n+nm       ; count +=1
                rows[count] = k+1n÷3+(μ-1)n+nm
                cols[count] = j+(ν-1)n+nm       ; count +=1

                rows[count] = k+2n÷3+(μ-1)n
                cols[count] = j+(ν-1)n+nm       ; count +=1
                rows[count] = k+2n÷3+(μ-1)n+nm
                cols[count] = j+(ν-1)n+nm       ; count +=1
            end

            for i ∈ 1:n
                rows[count] = (μ-1)n+i
                cols[count] = (ν-1)n+index      ; count +=1

                rows[count] = nm+(μ-1)n+i
                cols[count] = (ν-1)n+index      ; count +=1

                rows[count] = (μ-1)n+i
                cols[count] = nm+(ν-1)n+index   ; count +=1

                rows[count] = nm+(μ-1)n+i
                cols[count] = nm+(ν-1)n+index   ; count +=1
            end
            if μ == ν
                for i ∈ 1:n
                    rows[count] = (μ-1)n+i
                    cols[count] = nm+(ν-1)n+index ; count +=1

                    rows[count] = nm+(μ-1)n+i
                    cols[count] = nm+(ν-1)n+index ; count +=1
                end
            end
        end

        rowsA = rowvals(Aμ)
        for col ∈ 1:n
            if col!==index
                for j ∈ nzrange(Aμ, col)
                    row = rowsA[j]

                    rows[count] = (μ-1)n+row
                    cols[count] = (μ-1)n+col    ; count +=1

                    rows[count] = (μ-1)n+nm+row
                    cols[count] = (μ-1)n+col    ; count +=1

                    rows[count] = (μ-1)n+row
                    cols[count] = (μ-1)n+nm+col ; count +=1

                    rows[count] = (μ-1)n+nm+row
                    cols[count] = (μ-1)n+nm+col ; count +=1
                end
            end
        end
    end
    c = count - 1
    return sparse(rows[1:c],cols[1:c],vals[1:c],2nm,2nm)
end

################################################################################
# dictionary between real data vector and fields/frequencies

function ωψ_to_x(ωs::Vector,ψs::VectorField)
    x = Vector{Float64}(undef,2prod(size(ψs)))
    ωψ_to_x!(x,ωs,ψs)
    return x
end

ωψ_to_x!(ms::SALTProblem) = ωψ_to_x!(ms.x,ms.ωs,ms.ψs)
@inline function ωψ_to_x!(x::Vector,ωs::Vector,ψs::VectorField)
    n,m = size(ψs)
    nm = n*m
    index = n÷2+INDEX_OFFSET
    @inbounds for μ ∈ 1:m
        @inbounds for i ∈ 1:n
            if i == index
                x[(μ-1)n+index] = abs(ψs[index,μ])
                x[nm+(μ-1)n+index] = ωs[μ]
            else
                iszero(ψs[index,μ]) ? throw("field vanishes at $index, perhaps due to symmetry") : nothing
                x[(μ-1)n+i], x[nm+(μ-1)n+i] = reim(ψs[i,μ]/ψs[index,μ])
            end
        end
    end
    return nothing
end

x_to_ωψ!(ms::SALTProblem) = x_to_ωψ!(ms.ωs,ms.ψs,ms.x)
x_to_ωψ!(ms::SALTProblem,x) = x_to_ωψ!(ms.ωs,ms.ψs,x)
@inline function x_to_ωψ!(ωs::Vector,ψs::VectorField,x::Vector)
    n,m = size(ψs)
    nm = n*m
    index = n÷2+INDEX_OFFSET
    @inbounds for μ ∈ 1:m
        @inbounds for i ∈ 1:n ψs[i,μ] = complex(x[(μ-1)n+i],x[nm+(μ-1)n+i])*x[(μ-1)n+index] end
        ψs[index,μ] = complex(x[(μ-1)n+index],0)
        ωs[μ] = x[nm+(μ-1)n+index]
    end
    return nothing
end

################################################################################
# the salt bootstrap *UNDER CONSTRUCTION*

struct ThresholdFinder{TSIM,TLS,TKW}
    sim::TSIM
    tls::TLS
    nev::Int
    maxit::Int
    kwargs::TKW
end

function (tf::ThresholdFinder)(D₀::Real)
    tf.tls.D₀ = D₀
    nep = maxwell_nep(tf.sim)
    ω, _ = maxwell_eigen(nep, tf.tls.ωa; nev=tf.nev, maxit=tf.maxit, tf.kwargs...)
    return maximum(imag(ω))
end

function SALTbootstrap(
    sim::Simulation{N},
    tls::TwoLevelSystem;
    nev::Integer=5,
    maxit::Integer=50,
    δ::Real=.1,
    D0::Real=0,
    xatol::Real=1e-6,
    xrtol::Real=0,
    verbose::Bool=false,
    warn::Bool=true,
    iterations::Integer=150,
    kwargs...
    ) where N

    ωs = Vector{Float64}[]
    Ds = Float64[]
    ψs = ElectricField{N}[]

    ωᵗʰ, ψᵗʰ = findthreshold!(tls,sim,D0,δ,warn,verbose,xatol,xrtol,nev,maxit; kwargs...)
    push!(ωs,[ωᵗʰ])
    push!(ψs,0ψᵗʰ)
    push!(Ds,tls.D₀)

    if verbose
        print(stdout,"solving for ")
        printstyled(stdout,"D₀=",fmt("3.3f",tls.D₀),color=PRINTED_COLOR_NUMBER)
        println(stdout)
    end
    tls.D₀ += δ
    salt = maxwell_salt(sim,1)
    nlsolve(salt, ωs[1], δ*ψᵗʰ; iterations=iterations)

    nep = maxwell_nep(sim,salt.ωs,salt.ψs)
    ω, ψ = maxwell_eigen(nep,tls.ωa; nev=nev, maxit=maxit)
    imω,ind = findmax(imag(ω))

    if imω < xatol
        push!(ωs,salt.ωs)
        push!(ψs,salt.ψs)
        push!(Ds,tls.D₀)
        nmodes = 1
    else
        nmodes = 2
    end
    salt = maxwell_salt(sim,nmodes)
    salt.ωs[:] = ω[1]
    salt.ψs[:]
    nlsolve(salt; iterations=iterations)

    nep = maxwell_nep(sim,salt.ωs,salt.ψs)
    ω, ψ = maxwell_eigen(nep,tls.ωa; nev=nev, maxit=maxit)
    imω = maximum(imag(ω))

    return Ds,ωs,ψs
end

function findthreshold!(
            tls::TwoLevelSystem,
            sim::Simulation,
            D0::Real,
            δ::Real,
            warn::Bool,
            verbose::Bool,
            xatol::Real,
            xrtol::Real,
            nev::Integer,
            maxit::Integer;
            kwargs...
            )

    verbose ? println(stdout,"Computing threshold") : nothing
    tls.D₀ = D0
    nep = maxwell_nep(sim)
    ω, _ = maxwell_eigen(nep,tls.ωa; nev=nev, maxit=maxit)
    Δ = abs(-(extrema(real(ω))...))
    if warn && Δ<tls.γp
        printstyled(stdout,"\twarning",color=PRINTED_COLOR_WARN,bold=true)
        printstyled(stdout," frequency spread ", fmt("3.2f",Δ), " < ", fmt("3.2f",2*tls.γp)," = 2γ⟂,",color=PRINTED_COLOR_WARN)
        printstyled(stdout," consider increasing ",color=PRINTED_COLOR_WARN)
        printstyled(stdout,"nev\n",color=PRINTED_COLOR_VARIABLE)
    end
    imω1 = maximum(imag(ω)); imω2 = imω1
    iteration = 1
    verbose ? println(stdout,"\trough estimation") : nothing
    thresh = ThresholdFinder(sim,tls,nev,maxit,kwargs)
    while imω2 < 0
        D0 += δ
        imω1 = imω2
        if verbose
            print(stdout,"\t\titeration ",iteration)
            printstyled(stdout," D₀=",fmt("3.3f",tls.D₀),color=PRINTED_COLOR_NUMBER)
            println(stdout)
            iteration += 1
        end
        imω2 = thresh(D0)
    end
    verbose ? println(stdout,"\trefining threshold") : nothing
    tls.D₀ = find_zero(thresh,(tls.D₀-δ,tls.D₀); xatol=xatol,xrtol=xrtol)
    nep = maxwell_nep(sim)
    ω, ψ = maxwell_eigen(nep,tls.ωa; nev=nev, maxit=maxit, kwargs...)
    _, ind = findmax(imag(ω))
    return real(ω[ind]), ψ(ind)
end

################################################################################
# Pretty Printing

import ..Common.PRINTED_COLOR_LIGHT
import ..Common.PRINTED_COLOR_WARN
import ..Common.PRINTED_COLOR_VARIABLE
import ..Common.PRINTED_COLOR_GOOD
import ..Common.PRINTED_COLOR_BAD
import ..Common.PRINTED_COLOR_INSTRUCTION

function Base.show(io::IO,ms::SALTProblem)
    if !get(io, :SCPA, false)
        print(io,ms.m, " mode")
        ms.m>1 ? print(io,"s") : nothing
        if typeof(ms)<:MaxwellSALT
            printstyled(io," MaxwellSALT",color = PRINTED_COLOR_LIGHT)
        else
            printstyled(io," HelmholtzSALT",color = PRINTED_COLOR_LIGHT)
        end
    end
    if !ms.solved[]
        printstyled(io," (",color=PRINTED_COLOR_INSTRUCTION)
        printstyled(io,"unsolved",color=PRINTED_COLOR_WARN)
        printstyled(io,", pass to ",color=PRINTED_COLOR_INSTRUCTION)
        printstyled(io,"NLsolve.nlsolve",color=PRINTED_COLOR_VARIABLE)
        printstyled(io,")",color=PRINTED_COLOR_INSTRUCTION)
    elseif ms.converged[]
        printstyled(io," (",color=PRINTED_COLOR_INSTRUCTION)
        printstyled(io,"solved",color=PRINTED_COLOR_GOOD)
        printstyled(io,", solution in fields ",color=PRINTED_COLOR_INSTRUCTION)
        printstyled(io,"ωs",color=:cyan)
        printstyled(io,", ",color=PRINTED_COLOR_INSTRUCTION)
        printstyled(io,"ψs",color=:cyan)
        printstyled(io,")",color=PRINTED_COLOR_INSTRUCTION)
    else
        printstyled(io," (",color=PRINTED_COLOR_INSTRUCTION)
        printstyled(io,"solve attempted",color=PRINTED_COLOR_BAD)
        printstyled(io,", ",color=PRINTED_COLOR_INSTRUCTION)
        printstyled(io,"no convergence",color=PRINTED_COLOR_BAD)
        printstyled(io,")",color=PRINTED_COLOR_INSTRUCTION)
    end
end

################################################################################
# Plotting

@recipe f(ms::SALTProblem; by=abs2) = ms, by
@recipe f(by::Function,ms::SALTProblem) = ms, by
@recipe f(ms::SALTProblem, by::Function) = ms.ψs, by

@recipe f(sim::Simulation, ms::SALTProblem; by=abs2) = sim, ms, by
@recipe f(ms::SALTProblem, sim::Simulation; by=abs2) = sim, ms, by

@recipe f(ms::SALTProblem, sim::Simulation, by::Function) = sim, ms, by
@recipe f(sim::Simulation, by::Function, ms::SALTProblem) = sim, ms, by
@recipe f(ms::SALTProblem, by::Function, sim::Simulation) = sim, ms, by
@recipe f(by::Function, sim::Simulation, ms::SALTProblem) = sim, ms, by
@recipe f(by::Function, ms::SALTProblem, sim::Simulation) = sim, ms, by

@recipe f(sim::Simulation, ms::SALTProblem, by::Function) = sim, ms.ψs, by

end # module


text:
using Traceur, Test

warns_for(ws, x) = any(w -> occursin(x, w.message), ws)
warns_for(ws, x, xs...) = warns_for(ws, x) && warns_for(ws, xs...)

y = 1
const cy = 2
naive_relu(x) = x < 0 ? 0 : x

randone() = rand() < 0.5 ? Int(1) : rand() < 0.5 ? Float64(1) : BigInt(1)

function naive_sum(xs)
  s = 0
  for x in xs
    s += x
  end
  return s
end

f(x) = x+y

function f2(x)
  foo = y
  sin(x)+y
end

g(x) = x+cy

function test(warnings)
  ws = warnings(() -> naive_relu(1))
  @test isempty(ws)

  ws = warnings(() -> naive_relu(1.0))
  @test warns_for(ws, "returns")

  ws = warnings(() -> randone())
  @test warns_for(ws, "returns")

  ws = warnings(() -> naive_sum([1]))
  @test isempty(ws)

  ws = warnings(() -> naive_sum([1.0]))
  @test warns_for(ws, "assigned", "returns")

  ws = warnings(() -> f(1))
  @test warns_for(ws, "global", "dispatch", "returns")

  ws = warnings(() -> f2(1))
  @test warns_for(ws, "global", "dispatch", "returns")

  ws = warnings(() -> g(1))
  @test isempty(ws)
end

@testset "Traceur" begin
  @testset "Dynamic" begin
    test(Traceur.warnings)
  end
  # @testset "Static" begin
  #   test(Traceur.warnings_static)
  # end
  # @test_nowarn @trace naive_sum(1.0)
  # @test_nowarn @trace_static naive_sum(1.0)
end

ext: jl


text:
#!/usr/bin/env julia

while true
	if length(ARGS) == 0
		println("y")
	else
		println(ARGS[1])
	end
end

max_stars_repo_name: NotWearingPants/yes

max_issues_repo_name: NotWearingPants/yes

max_forks_repo_name: NotWearingPants/yes


text:
#=
     The basic idea is to find the smallest N such that
       string(rounded(lowerbound(x), N)) == string(rounded(upperbound(x), N))
       where rounded rounds fractional values to N fractional places

    Some intervals are too wide to do that.
    Some are not too wide but straddle an integer in a way that does not
      round to match values with the desired fractional resolution.

=#
# digitsRequired(bitsOfPrecision) = ceil(Int, bitsOfPrecision*0.3010299956639812)

function smartarbstring(x::ArbFloat{P}) where P
     digts = digitsRequired(P)
     if isexact(x)
        if isinteger(x)
            return string(x, digts)
        else
            s = rstrip(string(x, digts),'0')
            if s[end]=='.'
               s = string(s, "0")
            end
            return s
        end
     elseif radius(x) > abs(midpoint(x))
        return "0"
     end

     return smartarbstring_inexact(x)
end

function smartarbstring_inexact(x::ArbFloat{P}) where P
     digts = digitsRequired(P)
     lb, ub = bounds(x)
     lbs = string_exact(lb, digts)
     ubs = string_exact(ub, digts)
     if lbs[end]==ubs[end] && lbs==ubs
         return ubs
     end
     for i in (digts-2):-2:4
         lbs = string_exact(lb, i)
         ubs = string_exact(ub, i)
         if lbs[end]==ubs[end] && lbs==ubs # tests rounding to every other digit position
            us = string_exact(ub, i+1)
            ls = string_exact(lb, i+1)
            if us[end] == ls[end] && us==ls # tests rounding to every digit position
               ubs = lbs = us
            end
            break
         end
     end
     if lbs != ubs
        ubs = stringcompact(x)
     end
     s = rstrip(ubs,'0')
     if s[end]=='.'
         s = string(s, "0")
     end
     return s
end

function smartvalue(x::ArbFloat{P}) where {P}
    s = smartarbstring(x)
    ArbFloat{P}(s)
end

function smartstring(x::ArbFloat{P}) where {P}
    s = smartarbstring(x)
    a = ArbFloat{P}(s)
    if notexact(x)
       s = string(s,upperbound(x) < a ? '₋' : (lowerbound(x) > a ? '₊' : 'ₒ'))
    end
    return s
end

function smarterstring(x::T) where {T <: ArbFloat}
    absx   = abs(x)
    sa_str = smartarbstring(absx)  # smart arb string
    sa_val = (T)(absx)             # smart arb value
    if notexact(absx)
        md,rd = midpoint(absx), radius(absx)
        lo,hi = bounds(absx)
        if     sa_val <= lo
            if lo-sa_val >= ufp2(rd)
                sa_str = string(sa_str,"⁺")
            else
                sa_str = string(sa_str,"₊")
            end
        elseif sa_val > hi
            if sa_val-hi >= ufp2(rd)
                sa_str = string(sa_str,"⁻")
            else
                sa_str = string(sa_str,"₋")
            end
        else
            if sa_val-hi >= ufp2(rd)
                sa_str = string(sa_str,"ᵒ")
            else
                sa_str = string(sa_str,"ₒ")
            end
        end
    end
    return sa_str
end

#=

function stringTrimmed{P}(x::ArbFloat{P}, ndigitsremoved::Int)
   n = max(1, digitsRequired(P) - max(0, ndigitsremoved))
   cstr = ccall(@libarb(arb_get_str), Ref{UInt8}, (Ref{ArbFloat{P}}, Int, UInt), x, n, UInt(2))
   s = unsafe_string(cstr)
   # ccall(@libflint(flint_free), Void, (Ref{UInt8},), cstr)
   s
end

=#


text:
# wind correlation matrix
windCorrMat = [[1.0000 0.0900 -0.1671 0.0236 -0.2321 -0.0796 0.0386 0.2338 0.0254 0.0316 0.4308 0.0595 -0.0562 -0.1167 -0.0145 -0.0480 -0.0158 0.0574 -0.0170 -0.1551 0.0253 -0.1666 -0.3569 0.1767 0.0129 ];[0.0900 1.0000 -0.1955 -0.1956 -0.1528 -0.0547 -0.1986 -0.2806 0.1595 0.0184 0.1455 -0.2267 0.0615 0.1421 -0.2527 0.1115 -0.0125 -0.1935 -0.1131 0.1806 -0.0653 0.0120 -0.0134 -0.1175 0.1454 ];[-0.1671 -0.1955 1.0000 0.0738 -0.0027 0.0260 -0.0849 0.0715 -0.1539 -0.2518 -0.1267 -0.0457 -0.1182 -0.1653 0.1285 -0.0695 0.1521 0.2525 0.1654 0.0629 -0.0783 0.0978 0.2266 0.1822 -0.0016 ];[0.0236 -0.1956 0.0738 1.0000 0.0601 -0.0933 0.1587 0.1864 0.0099 -0.0522 0.0627 -0.1610 0.0019 -0.0916 0.1308 -0.1083 -0.1169 0.0702 -0.1382 0.0948 0.1405 -0.0604 -0.2338 -0.2503 -0.1195 ];[-0.2321 -0.1528 -0.0027 0.0601 1.0000 0.2109 -0.0879 0.2101 -0.0902 -0.0949 -0.1608 -0.1721 -0.0025 0.2146 0.0225 -0.0867 0.0488 0.0493 0.2465 0.0798 0.0607 0.2289 0.2221 -0.1012 0.2730 ];[-0.0796 -0.0547 0.0260 -0.0933 0.2109 1.0000 0.0893 0.0659 -0.3595 -0.0634 -0.1994 -0.2182 -0.2812 -0.1187 0.0728 -0.1836 -0.0392 0.1901 0.0868 0.0389 -0.0898 0.1094 0.1126 0.1334 0.2584 ];[0.0386 -0.1986 -0.0849 0.1587 -0.0879 0.0893 1.0000 0.2063 -0.2206 -0.0775 0.0506 0.1327 -0.1301 -0.1140 0.0935 -0.0286 -0.1008 0.1660 0.1411 -0.2069 0.1801 0.1161 -0.0300 -0.0097 -0.0660 ];[0.2338 -0.2806 0.0715 0.1864 0.2101 0.0659 0.2063 1.0000 -0.2555 -0.1014 -0.0875 -0.2149 -0.0423 -0.0284 -0.1237 -0.1518 -0.2133 0.1075 0.0694 -0.1516 -0.1077 0.2447 -0.1045 0.0495 0.1214 ];[0.0254 0.1595 -0.1539 0.0099 -0.0902 -0.3595 -0.2206 -0.2555 1.0000 0.0604 0.0715 0.0484 0.2200 0.0518 0.0973 0.1897 0.0071 -0.4850 -0.2098 0.3179 0.1750 -0.2224 -0.0445 -0.3726 0.0734 ];[0.0316 0.0184 -0.2518 -0.0522 -0.0949 -0.0634 -0.0775 -0.1014 0.0604 1.0000 0.0426 0.2193 0.1119 -0.0416 0.0903 -0.0397 0.0113 0.2000 -0.2587 -0.0341 0.1276 -0.2528 -0.0365 0.1113 0.0417 ];[0.4308 0.1455 -0.1267 0.0627 -0.1608 -0.1994 0.0506 -0.0875 0.0715 0.0426 1.0000 0.0757 -0.0275 0.1014 0.0560 0.1139 -0.0387 -0.1229 0.0691 -0.0877 0.0423 0.0202 -0.2556 0.2543 0.1074 ];[0.0595 -0.2267 -0.0457 -0.1610 -0.1721 -0.2182 0.1327 -0.2149 0.0484 0.2193 0.0757 1.0000 0.1498 -0.0024 -0.0643 0.0365 0.1767 0.0570 -0.1035 -0.0288 0.0045 -0.2481 -0.3069 0.1217 -0.1585 ];[-0.0562 0.0615 -0.1182 0.0019 -0.0025 -0.2812 -0.1301 -0.0423 0.2200 0.1119 -0.0275 0.1498 1.0000 -0.3105 -0.0724 0.2968 -0.0996 -0.0006 -0.1236 0.2789 -0.1327 -0.1194 -0.3212 -0.3220 -0.2096 ];[-0.1167 0.1421 -0.1653 -0.0916 0.2146 -0.1187 -0.1140 -0.0284 0.0518 -0.0416 0.1014 -0.0024 -0.3105 1.0000 -0.0318 -0.0834 0.2158 -0.2320 -0.0176 0.0139 -0.0673 0.1386 0.3092 0.2395 0.1904 ];[-0.0145 -0.2527 0.1285 0.1308 0.0225 0.0728 0.0935 -0.1237 0.0973 0.0903 0.0560 -0.0643 -0.0724 -0.0318 1.0000 -0.0014 0.0032 0.0600 0.0404 -0.0496 0.1949 -0.0531 0.0489 0.0482 -0.1634 ];[-0.0480 0.1115 -0.0695 -0.1083 -0.0867 -0.1836 -0.0286 -0.1518 0.1897 -0.0397 0.1139 0.0365 0.2968 -0.0834 -0.0014 1.0000 -0.2575 -0.2776 0.0310 0.0761 0.1641 0.0143 -0.1511 0.0566 -0.2025 ];[-0.0158 -0.0125 0.1521 -0.1169 0.0488 -0.0392 -0.1008 -0.2133 0.0071 0.0113 -0.0387 0.1767 -0.0996 0.2158 0.0032 -0.2575 1.0000 0.2479 -0.1291 -0.0763 -0.1233 0.0714 0.1789 0.1735 0.1685 ];[0.0574 -0.1935 0.2525 0.0702 0.0493 0.1901 0.1660 0.1075 -0.4850 0.2000 -0.1229 0.0570 -0.0006 -0.2320 0.0600 -0.2776 0.2479 1.0000 0.1322 -0.0408 -0.1321 0.0500 0.1147 0.0838 -0.0394 ];[-0.0170 -0.1131 0.1654 -0.1382 0.2465 0.0868 0.1411 0.0694 -0.2098 -0.2587 0.0691 -0.1035 -0.1236 -0.0176 0.0404 0.0310 -0.1291 0.1322 1.0000 -0.1243 -0.0681 0.1800 0.1639 0.0806 0.0964 ];[-0.1551 0.1806 0.0629 0.0948 0.0798 0.0389 -0.2069 -0.1516 0.3179 -0.0341 -0.0877 -0.0288 0.2789 0.0139 -0.0496 0.0761 -0.0763 -0.0408 -0.1243 1.0000 -0.1579 -0.0846 -0.0250 -0.1005 -0.0721 ];[0.0253 -0.0653 -0.0783 0.1405 0.0607 -0.0898 0.1801 -0.1077 0.1750 0.1276 0.0423 0.0045 -0.1327 -0.0673 0.1949 0.1641 -0.1233 -0.1321 -0.0681 -0.1579 1.0000 -0.0836 -0.1160 -0.0274 -0.0943 ];[-0.1666 0.0120 0.0978 -0.0604 0.2289 0.1094 0.1161 0.2447 -0.2224 -0.2528 0.0202 -0.2481 -0.1194 0.1386 -0.0531 0.0143 0.0714 0.0500 0.1800 -0.0846 -0.0836 1.0000 0.2955 0.1056 -0.0666 ];[-0.3569 -0.0134 0.2266 -0.2338 0.2221 0.1126 -0.0300 -0.1045 -0.0445 -0.0365 -0.2556 -0.3069 -0.3212 0.3092 0.0489 -0.1511 0.1789 0.1147 0.1639 -0.0250 -0.1160 0.2955 1.0000 0.1159 0.0133 ];[0.1767 -0.1175 0.1822 -0.2503 -0.1012 0.1334 -0.0097 0.0495 -0.3726 0.1113 0.2543 0.1217 -0.3220 0.2395 0.0482 0.0566 0.1735 0.0838 0.0806 -0.1005 -0.0274 0.1056 0.1159 1.0000 0.0422 ];[0.0129 0.1454 -0.0016 -0.1195 0.2730 0.2584 -0.0660 0.1214 0.0734 0.0417 0.1074 -0.1585 -0.2096 0.1904 -0.1634 -0.2025 0.1685 -0.0394 0.0964 -0.0721 -0.0943 -0.0666 0.0133 0.0422 1.0000 ]]

max_stars_repo_name: rohitkannan/DCOPF-reserve-saturation

max_issues_repo_name: rohitkannan/DCOPF-reserve-saturation

max_forks_repo_name: rohitkannan/DCOPF-reserve-saturation


text:
@testset "Test deprecated accessors for BatteryEMS" begin
    #TODO: use PSCB test system after PR#12 is merge and delete battery device
    # c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat_ems");
    # batt = get_component(BatteryEMS, c_sys5_bat, "Bat2");
    batt = BatteryEMS(;
        name = "Bat2",
        prime_mover = PrimeMovers.BA,
        available = true,
        bus = Bus(nothing),
        initial_energy = 5.0,
        state_of_charge_limits = (min = 0.10, max = 7.0),
        rating = 7.0,
        active_power = 2.0,
        input_active_power_limits = (min = 0.0, max = 2.0),
        output_active_power_limits = (min = 0.0, max = 2.0),
        efficiency = (in = 0.80, out = 0.90),
        reactive_power = 0.0,
        reactive_power_limits = (min = -2.0, max = 2.0),
        base_power = 100.0,
        storage_target = 0.2,
        operation_cost = StorageManagementCost(
            variable = VariableCost(0.0),
            fixed = 0.0,
            start_up = 0.0,
            shut_down = 0.0,
            energy_shortage_cost = 50.0,
            energy_surplus_cost = 40.0,
        ),
    )
    dummy_value = 99.9
    @test @test_deprecated(set_energy_value!(batt, dummy_value)) == dummy_value
    @test @test_deprecated(set_penalty_cost!(batt, dummy_value)) == dummy_value

    @test @test_deprecated(get_energy_value(batt)) ==
          get_energy_surplus_cost(get_operation_cost(batt))
    @test @test_deprecated(get_penalty_cost(batt)) ==
          get_energy_shortage_cost(get_operation_cost(batt))
end


text:
#
# Kalman Filter Types and Methods
#

abstract type AbstractKalmanFilter{T,S} <: SequentialEstimator{T,S} end

const AbstractDiscreteKalmanFilter{T} =
    AbstractKalmanFilter{T,UncertainDiscreteState{T}} where T
const AbstractContinuousKalmanFilter{T} =
    AbstractKalmanFilter{T,UncertainContinuousState{T}} where T

"""
    KalmanFilter(sys::LinearSystem, obs::LinearObserver,
                 estimate::AbstractState[, consider_states::Vector{UInt}])

Linear Kalman Filter type. The internal filter model takes on a specific form
depending on the type of the initial estimate. For a `DiscreteState`,

\$x_k = A x_{k-1} + w_k\$
\$y_k = H x_k + v_k\$

where \$w_k \\sim N(0, Q)\$ and \$v_k \\sim N(0, R)\$. For a `ContinuousState`,

\$\\dot{x}(t_k) = A x(t_k) + w(t_k)\$
\$y(t_k) = H x(t_k) + v(t_k)\$

Construction of a KalmanFilter requires an initial estimate. Internally, the
initial estimate is an uncertain state type. The constructor automatically
converts absolute initial estimates to uncertain states with zero covariance.

The `consider_states` input contains a list of indices of state elements to be
considered in the filtering process and not updated.
"""
immutable KalmanFilter{T,S<:AbstractUncertainState{T}} <: AbstractKalmanFilter{T,S}
    sys::LinearSystem{T}
    obs::LinearObserver{T}
    estimate::S
    consider_states::Vector{UInt16}

    function KalmanFilter(sys::LinearSystem{T}, obs::LinearObserver{T},
                          estimate::S, consider_states::Vector
                          ) where {T, S<:AbstractUncertainState{T}}
        assert_compatibility(sys, estimate)
        assert_compatibility(obs, estimate)
        if !allunique(consider_states)
            throw(ArgumentError("Consider state indices must be unique"))
        end
        for idx = 1:length(consider_states)
            if consider_states[idx] > length(estimate.x)
                throw(DimensionMismatch(
                    "Consider indices extend beyond length of initial state."))
            end
        end
        new{T,S}(sys, obs, estimate, consider_states)
    end
end
KalmanFilter(sys, obs, estimate) = KalmanFilter(sys, obs, estimate, [])
KalmanFilter(sys, obs, estimate::AbstractAbsoluteState, consider_states) =
    KalmanFilter(sys, obs, make_uncertain(estimate), consider_states)

const DiscreteKalmanFilter{T} =
    KalmanFilter{T,UncertainDiscreteState{T}} where T
const ContinuousKalmanFilter{T} =
    KalmanFilter{T,UncertainContinuousState{T}} where T

"""
    kalman_predict(kf::KalmanFilter, t)
"""
function kalman_predict(kf::KalmanFilter, t)
    xk = predict(kf.estimate, kf.sys, t)
    yk = predict(xk, kf.obs)
    Pxy = xk.P*kf.obs.H'
    return xk, yk, Pxy
end

"""
    kalman_update!(kf, yk, zk)
"""
function kalman_update!(kf::AbstractKalmanFilter, yk::AbstractUncertainState,
                        zk::AbstractState, Pxy::Matrix)
    Kk = Pxy*inv(yk.P)
    if isempty(kf.consider_states)
        kf.estimate.x += Kk*(zk.x - yk.x)
        kf.estimate.P -= Kk*yk.P*Kk'
    else
        prev_estimate = deepcopy(kf.estimate)
        kf.estimate.x += Kk*(zk.x - yk.x)
        kf.estimate.P -= Kk*yk.P*Kk'
        reset_consider_states!(prev_estimate, kf.estimate, kf.consider_states)
    end
end

"""
    reset_consider_states!(prev_state, new_state, consider_states)
"""
function reset_consider_states!(prev_state::AbstractUncertainState,
                                new_state::AbstractUncertainState,
                                consider_states::Vector)
   new_state.x[consider_states] .= prev_state.x[consider_states]
   for i = 1:length(new_state.x), j = 1:length(new_state.x)
       if (i in consider_states) && (j in consider_states)
           new_state.P[i,j] .= prev_state.P[i,j]
       end
   end
end

"""
    process!(kf::AbstractKalmanFilter, zk::AbstractAbsoluteState[, archive::EstimatorHistory])

Process a measurement with an arbitrary Kalman Filter. Predicts the filter's
internal estimate to the time of the input measurement and performs a
correction step.

The user may optionally provide an `EstimatorHistory` archive variable to store
the incremental state data produced during the process step.
"""
function process!(::AbstractKalmanFilter) end
function process!(kf::AbstractDiscreteKalmanFilter, zk::DiscreteState)
    xk, yk, Pxy = kalman_predict(kf, zk.t)
    kf.estimate .= xk
    kalman_update!(kf, yk, zk, Pxy)
    return nothing
end
function process!(kf::AbstractDiscreteKalmanFilter, zk::DiscreteState,
                  archive::EstimatorHistory)
    if length(archive.states) == 0
      push!(archive.states, deepcopy(kf.estimate))
    end

    xk, yk, Pxy = kalman_predict(kf, zk.t)
    kf.estimate .= xk
    kalman_update!(kf, yk, zk, Pxy)

    push!(archive.states, deepcopy(kf.estimate))
    push!(archive.residuals, UncertainDiscreteState(zk.x - yk.x, yk.P, zk.t))
    return nothing
end
function process!(kf::AbstractContinuousKalmanFilter, zk::ContinuousState)
    xk, yk, Pxy = kalman_predict(kf, zk.t)
    kf.estimate .= xk
    kalman_update!(kf, yk, zk, Pxy)
    return nothing
end
function process!(kf::AbstractContinuousKalmanFilter, zk::ContinuousState,
                  archive::EstimatorHistory)
    if length(archive.states) == 0
      push!(archive.states, deepcopy(kf.estimate))
    end

    xk, yk, Pxy = kalman_predict(kf, zk.t)
    kf.estimate .= xk
    kalman_update!(kf, yk, zk, Pxy)

    push!(archive.states, deepcopy(kf.estimate))
    push!(archive.residuals, UncertainContinuousState(zk.x - yk.x, yk.P, zk.t))
    return nothing
end

"""
    kalman_filter!(state, system, observer, measurement[, archive])

Kalman Filter convenience function. Performs a single Kalman Filter iteration,
updating the provided state in-place. This function is not as efficient as
constructing creating a static `KalmanFilter` and using the `process!` method,
but allows for increased flexibility.

Use `kalman_filter!` for not-in-place state updates.
"""
function kalman_filter!(state, system, observer, measurement)
    kf = KalmanFilter(system, observer, state)
    process!(kf, measurement)
end
function kalman_filter!(state, system, observer, measurement, archive)
    kf = KalmanFilter(system, observer, state)
    process!(kf, measurement, archive)
end

"""
kalman_filter(state, system, observer, measurement[, archive])

Kalman Filter convenience function. Performs a single Kalman Filter iteration.
This function is not as efficient as constructing creating a static
`KalmanFilter` and using the `process!` method, but allows for increased
flexibility.

Use `kalman_filter!` for in-place state updates.
"""
function kalman_filter(state, system, observer, measurement)
    out_state = deepcopy(state)
    kalman_filter!(out_state, system, observer, measurement)
    return out_state
end
function kalman_filter(state, system, observer, measurement, archive)
    out_state = deepcopy(state)
    kalman_filter!(out_state, system, observer, measurement, archive)
    return out_state
end


text:
using Vinyl: @reset, shift

# ---------------------------------------------------------------------------- #

# Simulate Functions
k = @reset 1 + shift(k -> k)

# This is equivalent to (x -> 1 + x)
# The reset block stops running at `shift`, and continues when we provide a
# value to `k`.

k(2) # 3
k(4) # 5

# ---------------------------------------------------------------------------- #

# Turning callbacks into straight-line code.
# Imagine this callback retreives a web page asynchronously.
load(cb, page) = (sleep(1);cb(page))

# Here's a JavaScript-style function that loads two pages. Awful!
function f(cb)
  load("page1") do page1
    load("page2") do page2
      cb((page1,page2))
    end
  end
end

f(print) # ("page1", "page2")

# Here's a version with continuations.
function f()
  @reset begin
    page1 = shift(k -> load(k, "page1"))
    page2 = shift(k -> load(k, "page2"))
    return (page1, page2)
  end
end

f() # ("page1", "page2")

# ---------------------------------------------------------------------------- #

# What if Julia didn't have exceptions?

throwit(x) = (_ = shift(_ -> x))

function wrong()
  println("You'll see this")
  _ = throwit(ErrorException("An error occured!"))
  println("You won't see this")
end

x = @reset wrong()

# ---------------------------------------------------------------------------- #

# Hijacking Julia's control flow; we can explore both branches of an `if`
# statement.

quantum_predicate() = shift(k -> (k(true), k(false)))

function foo(x)
  quantum_predicate() && (x = -x)
  2x
end

@reset foo(5) # (-10, 10)

# Cartesian product of two bits:
@reset (quantum_predicate(), quantum_predicate())


text:
# ------------------------------------------------------------------------------------------
# # Timeseries Prediction
#
# ** EVERYTHING IN THIS TUTORIAL IS VERY NEW! CONSIDER IT ON AN ALPHA PHASE!**
#
# **APIs are likely to change!**
#
# **Examples shown in this tutorial are also in the repo
# `TimeseriesPrediction.jl/examples`**
#
# Topics:
#
# * Nature of prediction models of **DynamicalSystems.jl**
# * Local model prediction
# * Multi-variate local model Prediction
# * Spatio-temporal Timeseries Prediction
# * Docstrings
#
# ## Nature of prediction models
# Suppose you have a scalar or multi-variate timeseries and you want to predict its future
# behaviour.
#
# You can either take your *neural-network/machine-learning hammer* and lots of computing
# power **or** you can use methods from nonlinear dynamics and chaos.
#
# **DynamicalSystems.jl** follows the second approach. This road is not only surprisingly
# powerful, but also much, **much** simpler.
#
# ---
#
# # Local Model Prediction
#
# Local model prediction does something very simple: it makes a prediction of a state, by
# finding the future of similar (*neighboring*) states! Then it uses the predicted state as
# a new state from which other predictions can be made!
#
# Yeap, that simple.
#
# Let's see how well this method fares in a simple system, the Roessler system (3D &
# chaotic):
#
# $$
# \begin{aligned}
# \dot{x} &= -y-z \\
# \dot{y} &= x+ay \\
# \dot{z} &= b + z(x-c)
# \end{aligned}
# $$
# ------------------------------------------------------------------------------------------

using DynamicalSystems

# This initial condition gives a good prediction:
u0_good = [0.065081, 0.917503, 0.300242]

ross = Systems.roessler(u0_good)

# ------------------------------------------------------------------------------------------
# Let's get a "measurement" from the roessler system
# ------------------------------------------------------------------------------------------

dt = 0.1 # sampling rate
tf = 1000.0 # final time
tr = trajectory(ross, tf; dt = dt)

# This is the measurement
s = tr[50:end, 2] # we skip the first points, they are transient
# This is the accompanying time vector:
timevec = collect(0:dt:tf)[50:end];

# ------------------------------------------------------------------------------------------
# How does this timeseries look?
# ------------------------------------------------------------------------------------------

using PyPlot; figure(figsize = (8,4))
plot(timevec, s, lw = 1.0);

# ------------------------------------------------------------------------------------------
# Please note: these are chaotic oscillations, the system is *not* periodic for the chosen
# (default) parameter values!
#
# Alright, so we have a recorded some timeseries of length:
# ------------------------------------------------------------------------------------------

length(s)

# ------------------------------------------------------------------------------------------
# And now we want to predict!
#
# Let's see the prediction function in action! The function to use is
# ```julia
# localmodel_tsp(s, D::Int, τ, p::Int; kwargs...)
# ```
# Here `s` is the timeseries to be predicted. `D, τ` are the values of the `Reconstruction`
# that has to be made from `s`. The last argument `p` is simply the amount of points to
# predict!
#
# The `Reconstruction` idea and functions were introduced in the tutorial "Delay Coordinates
# Embedding".
#
# This local model prediction method assumes that the system is on some kind of chaotic
# attractor. This is why it is crucial to reconstruct a signal before using the method!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's use only a first part of the timeseries as a "training set"
# ------------------------------------------------------------------------------------------

N = length(s)
N_train = 1000
s_train = s[1:N_train];

# ------------------------------------------------------------------------------------------
# and then use the rest of the timeseries to compare with the prediction
# ------------------------------------------------------------------------------------------

s_test = s[N_train+1:end];

# ------------------------------------------------------------------------------------------
# Here we define the parameters to make a prediction:
# ------------------------------------------------------------------------------------------

τ = 17
D = 3
p = 500

s_pred = localmodel_tsp(s_train, D, τ, p)
# prediction always includes last point of `s_train`

# ------------------------------------------------------------------------------------------
# Let's plot!
# ------------------------------------------------------------------------------------------

figure(figsize=(8,4))
past = 100
plot(timevec[N_train-past:N_train+1], s[N_train-past:N_train+1], color = "C1", label = "timeseries")
plot(timevec[N_train:N_train+p], s[N_train:N_train+p], color = "C3", label = "real future")
plot(timevec[N_train:N_train+p], s_pred, color = "C0", linestyle = "dashed", alpha = 0.5, label = "prediction")
legend(); xlabel("\$t\$"); ylabel("\$y\$")
println("Prediction of $(p) points from $(N_train) points. i.c.: $(get_state(ross))")

# ------------------------------------------------------------------------------------------
# Of course the prediction depends strongly on:
#
# * Choosing proper `Reconstruction` parameters
# * The initial condition
#
# How did I know that the value of `τ=17` was good?
# ------------------------------------------------------------------------------------------

estimate_delay(s, "first_zero")

# ------------------------------------------------------------------------------------------
# The function `localmodel_tsp` also accepts some keyword arguments which I did not discuss.
# These are:
#
#   * `method = AverageLocalModel(2)` : Subtype of [`AbstractLocalModel`](@ref).
#   * `ntype = FixedMassNeighborhood(2)` : Subtype of [`AbstractNeighborhood`](@ref).
#   * `stepsize = 1` : Prediction step size.
#
# We already know what does the `ntype` keyword does: it chooses a neighborhood type.
#
# The `method` keyword chooses the method of the local prediction. There are two methods,
# the `AverageLocalModel`, which we already used by default, as well as the
# `LinearLocalModel`. Their docstrings are at the end of this tutorial.
#
# Without explanations, their call signatures are:
# ```julia
# AverageLocalModel(n::Int)
# LinearLocalModel(n::Int, μ::Real)
# LinearLocalModel(n::Int, s_min::Real, s_max::Real)
# ```
# ------------------------------------------------------------------------------------------

using BenchmarkTools
@btime localmodel_tsp($s_train, $D, $τ, $p)
println("Time for predicting $(p) points from $(N_train) points.")

# ------------------------------------------------------------------------------------------
# Let's bundle all the production-prediction-plotting process into one function and play
# around!
# ------------------------------------------------------------------------------------------

function predict_roessler(N_train, p, method, u0 = rand(3); ntype = FixedMassNeighborhood(5))

    ds = Systems.roessler(u0)
    dt = 0.1
    tr = trajectory(ds, (N_train+p)÷dt; dt = dt)

    s = tr[:, 2] # actually, any of the 3 variables of the Roessler work well

    s_train = s[1:N_train]
    s_test = s[N_train+1:end]

    # parameters to predict:
    τ = 17
    D = 3

    s_pred = localmodel_tsp(s_train, D, τ, p; method = method, ntype = ntype)

    figure(figsize=(8,4))
    past = 100
    plot(timevec[N_train-past:N_train+1], s[N_train-past:N_train+1], color = "C1", label = "timeseries")
    plot(timevec[N_train:N_train+p], s[N_train:N_train+p], color = "C3", label = "real future")
    plot(timevec[N_train:N_train+p], s_pred, color = "C0", linestyle = "dashed", alpha = 0.5, label = "prediction")
    legend(); xlabel("\$t\$"); ylabel("\$y\$")
    mprint = Base.datatype_name(typeof(method))
    println("N_train = $(N_train), p = $(p), method = $(mprint), u0 = $(u0)")
    return
end

predict_roessler(3000, 500, LinearLocalModel(2, 5.0))
# Linear Local model is slower than Average local model, and in general not that
# much more powerful.

# ------------------------------------------------------------------------------------------
# # Multi-Variate Prediction
#
# On purpose I was always referring to `s` as "timeseries". There is no reason for `s` to be
# scalar though, this prediction method works just as well when predicting multiple
# timeseries. And the call signature does not change at all!
#
# The following example demonstrates the prediction of the Lorenz96 model
#
# $$
# \frac{dx_i}{dt} = (x_{i+1}-x_{i-2})x_{i-1} - x_i + F
# $$
#
# a system that displays high-dimensional chaos and is thus very difficult to predict!
# ------------------------------------------------------------------------------------------

using DynamicalSystems, PyPlot

#Generate timeseries set
ds = Systems.lorenz96(5; F=8.)
ic = get_state(ds)
Δt = 0.05
s = trajectory(ds, 2100; dt=Δt)[:,1:2]

#Set Training and Test Set
N_train = 40000
p = 200
s_train = s[1:N_train,1:2]
s_test  = s[N_train:N_train+p,1:2]

#Embedding Parameters
D = 5; # total dimension of reconstruction is D*2 ! ! !
x = s[:, 1]
τ = estimate_delay(x, "first_zero")
println("Delay time estimation: $(τ)")

#Prediction
method = LinearLocalModel(2, 2.5)
method = AverageLocalModel(2)
ntype = FixedMassNeighborhood(5)
s_pred  = localmodel_tsp(s_train, D, τ, p; method = method, ntype = ntype)

figure(figsize=(12,4))
ax = subplot(121)
plot((N_train:N_train+p)*Δt, s_test[:,1], label="signal")
plot((N_train:N_train+p)*Δt, s_pred[:,1], label="prediction")
ylabel("\$x_1(n)\$")
xlabel("\$t\$")
legend()
ax = subplot(122)
plot((N_train:N_train+p)*Δt, s_test[:,2], label="signal")
plot((N_train:N_train+p)*Δt, s_pred[:,2], label="prediction")
ylabel("\$x_2(n)\$")
xlabel("\$t\$")
legend()
tight_layout();
println("Prediction p=$p of Lorenz96 Model (5 Nodes) from $N_train points")
println("i.c.: $ic")

# ------------------------------------------------------------------------------------------
# # Spatio Temporal Timeseries prediction
#
# Spatio-temporal systems are systems that depend on both space and time, i.e. *fields*
# (like Partial Differential Equations). These systems can also be predicted using these
# local model methods!
#
# In the following sections we will see 3 examples, but there won't be any code shown for
# the last example (because its humongus).
#
# See `TimeseriesPrediction.jl/examples` repository for more examples!
#
# ## Barkley Model
# The Barkley model consists of 2 coupled fields each having 2 spatial dimensions, and is
# considered one of the simplest spatio-temporal systems
#
# $$
# \begin{align}
# \frac{\partial u }{\partial t} =& \frac{1}{\epsilon} u (1-u)\left(u-\frac{v+b}{a}\right) +
# \nabla^2 u\nonumber \\
# \frac{\partial v }{\partial t} =& u - v
# \end{align}
# $$
#
# ------------------------------------------------------------------------------------------

# This Algorithm of evolving the Barkley model is taken from
# http://www.scholarpedia.org/article/Barkley_model
function barkley(T, Nx, Ny)
    a = 0.75; b = 0.02; ε = 0.02

    u = zeros(Nx, Ny); v = zeros(Nx, Ny)
    U = Vector{Array{Float64,2}}()
    V = Vector{Array{Float64,2}}()

    #Initial state that creates spirals
    u[40:end,34] = 0.1; u[40:end,35] = 0.5
    u[40:end,36] = 5; v[40:end,34] = 1
    u[1:10,14] = 5; u[1:10,15] = 0.5
    u[1:10,16] = 0.1; v[1:10,17] = 1
    u[27:36,20] = 5; u[27:36,19] = 0.5
    u[27:36,18] = 0.1; v[27:36,17] = 1

    h = 0.75; Δt = 0.1; δ = 0.001
    Σ = zeros(Nx, Ny, 2)
    r = 1; s = 2

    function F(u, uth)
        if u < uth
            u/(1-(Δt/ε)*(1-u)*(u-uth))
        else
            (u + (Δt/ε)*u*(u-uth))/(1+(Δt/ε)*u*(u-uth))
        end
    end

    for m=1:T
        for i=1:Nx, j=1:Ny
            if u[i,j] < δ
                u[i,j] = Δt/h^2 * Σ[i,j,r]
                v[i,j] = (1 - Δt)* v[i,j]
            else
                uth = (v[i,j] + b)/a
                v[i,j] = v[i,j] + Δt*(u[i,j] - v[i,j])
                u[i,j] = F(u[i,j], uth) + Δt/h^2 *Σ[i,j,r]
                Σ[i,j,s] -= 4u[i,j]
                i > 1  && (Σ[i-1,j,s] += u[i,j])
                i < Nx && (Σ[i+1,j,s] += u[i,j])
                j > 1  && (Σ[i,j-1,s] += u[i,j])
                j < Ny && (Σ[i,j+1,s] += u[i,j])
            end
            Σ[i,j,r] = 0
        end
        r,s = s,r
        #V[:,:,m] .= v
        #U[:,:,m] .= u
        push!(U,copy(u))
        push!(V,copy(v))
    end
    return U,V
end

Nx = 50
Ny = 50
Tskip = 100
Ttrain = 500
p = 20
T = Tskip + Ttrain + p

U,V = barkley(T, Nx, Ny);

Vtrain = V[Tskip + 1:Tskip + Ttrain]
Vtest  = V[Tskip + Ttrain :  T]

D = 2
τ = 1
B = 2
k = 1
c = 200

Vpred = localmodel_stts(Vtrain, D, τ, p, B, k; boundary=c)
err = [abs.(Vtest[i]-Vpred[i]) for i=1:p+1];
println("Maximum error: $(maximum(maximum.(err)))")

# ------------------------------------------------------------------------------------------
# Plotting the real evolution, prediction, and error side by side (plotting takes a *lot* of
# time) with `Tskip = 200, Ttrain = 1000, p = 200` produces:
#
# ![Barkley prediction](barkley_stts_prediction.gif)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Kuramoto-Sivashinsky
#
# This system consists of only a single field and a single spatial dimension
#
# $$
# y_t =  −y y_x − y_{xx} − y_{xxxx}
# $$
# where the subscripts denote partial derivatives.
#
# We will not show code for this example, but only the results. The code is located at
# `TimeseriesPrediction.jl/examples/KSprediction.jl` and is quite large for a Jupyter
# notebook.
#
# ### Prediction
#
# A typical prediction with parameters `D = 2, τ = 1, Β = 5, k = 1` and system parameters `L
# = 150` (see file for more) looks like this:
#
# ![Prediction of KS system](KS150_tr40000_D2_tau1_B5_k1.png)
#
# The vertical axis, which is *time*, is measured within units of the maximal Lyapunov
# exponent of the system Λ, which is around 0.1.
#
# ### Mean Squared Error of prediction
#
# We now present a measure of the error of the prediction, by averaging the error values
# across all timesteps and across all spatial values, and then normalizing properly
#
# ![Error of prediction](KS_NRMSE_L6_Q64_D1_tau1_B5_k1_nn4_nw3PndWDWSt.png)
#
# The above curves are also *averaged* over 10 different initial conditions and subsequent
# predictions!
#
# The parameters used for the prediction
#
# You can compare it with e.g. figure 5 of this [Physical Review Letters
# article](https://doi.org/10.1103/PhysRevLett.120.024102).
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Docstrings
# ------------------------------------------------------------------------------------------

?localmodel_tsp

?AbstractLocalModel

?localmodel_stts

?STReconstruction

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
""" Day19 Advent of Code 2021 """
function day19()
    part = [0, 0]

    scannerblocks = split(read("AoCdata/AoC_2021_day19.txt", String), "\n\n")
    scanners = [Vector{Int}[] for _ in 1:length(scannerblocks)]
    for (i, block) in enumerate(scannerblocks), line in strip.(split(block, "\n"))
        !isempty(line) && line[end] != '-' && push!(scanners[i], parse.(Int, split(line, ",")))
    end
    converted, checked, distances = Set(1), Set{Int}(), Dict{Int, Tuple{Int, Int, Int}}()
    sort!(scanners[1])

    # simulate all configurations of scanner j and then look for a set of 12 between the two
    # which differ by a constant xyz vector such as [0, 0, 0].
    while length(converted) < length(scanners)
        for i in converted
            i in checked && continue
            push!(checked, i)
            mat = sort(scanners[i])
            for j in 1:length(scanners)
                j ∈ converted && continue
                print("Compare $i to $j...\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b")
                for xsign in [1, -1], ysign in [1, -1], zsign in [1, -1], ax in ([1,2,3], [2,1,3], [3,2,1], [1,3,2], [2,3,1], [3,1,2])
                    newmat = deepcopy(scanners[j])
                    for i in eachindex(newmat)
                        newmat[i] .*= [xsign, ysign, zsign]
                    end
                    for (i, t) in enumerate(newmat)
                        newmat[i] .= [t[ax[1]], t[ax[2]], t[ax[3]]]
                    end
                    for b in mat, t in newmat
                        diff = t .- b
                        minindex = min(length(mat), length(newmat))
                        samediffs = count(newmat[n] - mat[m] == diff for n in 1:minindex, m in 1:minindex)
                        if samediffs > 4
                            push!(converted, j)
                            scanners[j] .= [p .-= diff for p in sort!(newmat)]
                            distances[j] = (diff[1], diff[2], diff[3])
                            break
                        end
                    end
                end
            end
        end
    end

    part[1] = length(unique(reduce(vcat, scanners)))
    part[2] = maximum(sum(abs.(d1 .- d2)) for d1 in values(distances), d2 in values(distances))

    return part
end

part = day19()

println("\nPart 1: ", part[1]) # 425
println("Part 2: ", part[2]) # 13354


text:
using Pkg

function get_previous_julia_version(current::VersionNumber=VERSION)::VersionNumber
    search_dir = "$(homedir())/.julia/environments"
    versions = sort(VersionNumber.(readdir(search_dir)))
    if versions[end] == current
        pop!(versions)
    end
    return versions[end]
end

function get_installed_packages(version::VersionNumber=VERSION)::Vector{String}
    pkgdir = "$(homedir())/.julia/environments/v$(version.major).$(version.minor)/"
    project_file = readlines("$pkgdir/Project.toml")
    pkgs = String[]
    pkg_list_active = false
    for (i, line) in enumerate(project_file)
        if line == "[deps]"
            pkg_list_active = true
            continue
        elseif pkg_list_active && line == ""
            break
        end
        if pkg_list_active
            pkgs = [pkgs; split(line, ' ')[1]]
        end
    end
    return pkgs
end

function transfer_pkgs()
    currently_installed = keys(Pkg.installed())
    previously_installed = get_installed_packages(get_previous_julia_version())
    for pkg in previously_installed
        print("\nLooking for package $pkg...")
        if pkg in currently_installed
            print("already installed.\n")
        else
            print("attempting installation.\n")
            try
                Pkg.add(pkg);
            catch
                continue
            end
        end
    end
    return nothing
end

transfer_pkgs()

max_stars_repo_name: dysonance/Preferences

max_issues_repo_name: dysonance/Preferences

max_forks_repo_name: dysonance/Preferences


text:
using Test, Documenter, DocumenterGHActionContainer
DocMeta.setdocmeta!(DocumenterGHActionContainer, :DocTestSetup, :(using DocumenterGHActionContainer), recursive = true)

@testset "DocumenterGHActionContainer" begin
    doctest(DocumenterGHActionContainer)
end


text:
using JuLIP
using JuLIP.ASE
using JuLIP.Potentials

using JuLIPMaterials.Elasticity: elastic_moduli

println("=========================")
println(" Testing Elasticity      ")
println("=========================")

println("Testing `elastic_moduli`")
at = bulk("Cu")
calc = lennardjones(r0 = rnn("Cu"))
C = elastic_moduli(calc, at)

C2d = C[1:2,1:2,1:2,1:2]
println( reshape(C2d, 4, 4) )


text:
using Distributions

using Profile
using BenchmarkTools

X = 40 * rand(30, 1) .+ 20;
Y = (2 * X .+ randn(30, 1)).^0.5
Y = reshape(Y, length(Y))

k = SquaredExponential(Uniform());
g = BoxCox(Uniform(1, 3))
f = Identity()
mdl = Model(k, g, f, X, Y, 30, 30);

x = [20.0]

show(equalinterval(mdl, 0.95, x))
display(plotdensity(mdl, 0.95, x, 20))

@benchmark equalinterval(mdl, 0.95, x)

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
@testset "283.move-zeroes.jl" begin
    nums = [0, 1, 0, 3, 12]
    move_zeros!(nums)
    @test nums == [1, 3, 12, 0, 0]

    nums = [0, 0, 1, 2]
    move_zeros!(nums)
    @test nums == [1, 2, 0, 0]

    nums = [0]
    move_zeros!(nums)
    @test nums == [0]
end


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   ReciprocalPath{P}
#   --> P is the reciprocal point type (<: AbstractReciprocalPoint{D})
#       --> D is the dimension of embedding space
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#       - TODO interface testing
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
abstract type AbstractReciprocalPath{P <: AbstractReciprocalPoint{D} where {D}} <: AbstractArray{P,1} end

# export the type
export AbstractReciprocalPath

# Array functionality for path type

# Indexing
function Base.IndexStyle(::Type{<:AbstractReciprocalPath})
    return IndexLinear()
end

# size
function Base.size(path::P) where {RP, P<:AbstractReciprocalPath{RP}}
    return size(points(path))
end

# getindex
function Base.getindex(path::P, i::Int64) where {RP, P<:AbstractReciprocalPath{RP}}
    return getindex(points(path), i)
end

# setindex
function Base.setindex!(path::P, v, i::Int64) where {RP, P<:AbstractReciprocalPath{RP}}
    return setindex!(points(path), v, i)
end

################################################################################
#
#	INTERFACING / ACCESSING RECIPROCAL PATHS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# default constructor interface
# used for creation of new reciprocal paths
function newReciprocalPath(
            ::Type{RPA},
            points :: Vector{<:RPO}
        ) :: RPA where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO}}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'newReciprocalPath' for concrete reciprocal path type " *
            string(RPA) * " with reciprocal point type " * string(RPO) *
            " in " * string(D) * " spatial dimensions" )
end

# export the newUnitcell interface
export newReciprocalPath

# access the list of points
function points(
            p :: RPA
        ) :: Vector{RPO} where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO}}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'points' for concrete reciprocal path type " * string(RPA) )
end
# set the list of points
function points!(
            p :: RPA,
            points :: Vector{<:RPO}
        ) where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO}}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'points!' for concrete reciprocal path type " * string(RPA) )
end

# export access to points
export points, points!

function saveReciprocalPath(
        rp :: RPA,
        fn :: AbstractString,
        group :: AbstractString = "reciprocal_path"
        ;
        append :: Bool = false
    ) where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO}}

end

function loadReciprocalPath(
        ::Type{RPA},
        fn :: AbstractString,
        group :: AbstractString = "reciprocal_path"
    ) where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO}}

end

# export access to points
export saveReciprocalPath, loadReciprocalPath

# some more beauty interface
# builds on interface defined above but can also be overwritten

# number of reciprocal points
function numPoints(
            p :: RPA
        ) :: Int64 where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO}}

    # return length of the point array that the reciprocal path type implements
    return length(points(p))
end

# export the function
export numPoints

# access a specific reciprocal point
function point(
            path  :: RPA,
            index :: Int64
        ) :: RPO where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO}}

    # return the respective reciprocal point
    return points(path)[index]
end

# export the function
export point

# splatted constructor interface
# used for creation of new reciprocal paths based on number of points
# (does not have to be overwritten)
function newReciprocalPath(
            ::Type{RPA},
            points :: RPO ...
        ) :: RPA where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO}}

    # return the usual constructor
    return newReciprocalPath(RPA, RPO[p for p in points])
end
# set the list of points with splatted input
function points!(
            p :: RPA,
            points :: RPO ...
        ) where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO}}

    # return the fitting function
    return points!(p, RPO[p for p in points])
end

################################################################################
#
#	OVERWRITING BASE.SHOW FOR RECIPROCAL PATHS
#
################################################################################
function show(io::IO, p::RPA) where {D, RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO}}
    print(io, "Reciprocal Path of type ", string(RPA), " including ", numPoints(p), " points.\n")
    print(io, "1) ", point(p,1))
    for i in 2:numPoints(p)
        print("\n    |\n")
        print(io, i, ") ", point(p,i))
    end
end

################################################################################
#
#	INTERFACING FOR BUILDING RECIPROCAL PATHS
#
################################################################################

# CONSTRUCT A PATH BY PASSED POINTS (as Symbols)
function getReciprocalPath(
            ::Type{RPA},
            unitcell :: U,
            points :: Symbol ...
        ) :: RPA where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO},S,B,U<:AbstractUnitcell{S,B}}

    # pass to the specific function
    return newReciprocalPath(RPA, [getReciprocalPoint(RPO, unitcell, p) for p in points])
end
function getReciprocalPath(
            ::Type{RPA},
            unitcell :: U,
            points :: Tuple{Symbol,Int64} ...
        ) :: RPA where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO},S,B,U<:AbstractUnitcell{S,B}}

    # pass to the specific function
    return newReciprocalPath(RPA, [getReciprocalPoint(RPO, unitcell, p[1],p[2]) for p in points])
end
function getReciprocalPath(
            ::Type{RPA},
            unitcell :: U,
            points :: Union{Symbol, Tuple{Symbol,Int64}} ...
        ) :: RPA where {D,RPO<:AbstractReciprocalPoint{D}, RPA<:AbstractReciprocalPath{RPO},S,B,U<:AbstractUnitcell{S,B}}

    # pass to the specific function
    return newReciprocalPath(RPA, [getReciprocalPoint(RPO, unitcell, p) for p in points])
end

function getReciprocalPath(
            unitcell :: U,
            points :: Symbol ...
        ) :: ReciprocalPath{ReciprocalPoint{D}} where {LS,D,S<:AbstractSite{LS,D},B,U<:AbstractUnitcell{S,B}}

    # pass to the specific function
    return newReciprocalPath(ReciprocalPath{ReciprocalPoint{D}}, [getReciprocalPoint(ReciprocalPoint{D}, unitcell, p) for p in points])
end
function getReciprocalPath(
            unitcell :: U,
            points :: Tuple{Symbol,Int64} ...
        ) :: ReciprocalPath{ReciprocalPoint{D}} where {LS,D,S<:AbstractSite{LS,D},B,U<:AbstractUnitcell{S,B}}

    # pass to the specific function
    return newReciprocalPath(ReciprocalPath{ReciprocalPoint{D}}, [getReciprocalPoint(ReciprocalPoint{D}, unitcell, p[1],p[2]) for p in points])
end
function getReciprocalPath(
            unitcell :: U,
            points :: Union{Symbol, Tuple{Symbol,Int64}} ...
        ) :: ReciprocalPath{ReciprocalPoint{D}} where {LS,D,S<:AbstractSite{LS,D},B,U<:AbstractUnitcell{S,B}}

    # pass to the specific function
    return newReciprocalPath(ReciprocalPath{ReciprocalPoint{D}}, [getReciprocalPoint(ReciprocalPoint{D}, unitcell, p) for p in points])
end

# export the reciprocal path
export getReciprocalPath


text:
#*************************************************************************#
#
# DifferenceEquation Type
#
#*************************************************************************#
"""
    DifferenceEquation(mdl::ParametricModel, yterms, uterms)

Type representing a difference equation. `mdl` is a static function.
It relates the output at a given instant with its previous values
and with the input. The `yterms` and `uterms` are vector of vectors
containing integers, each one specifying a dependence with previous
terms of `y` (or with input terms):

    y[k] = mdl(y_1[k-yterms[1][1]], y_1[k-yterms[1][2]],..., y_1[k-yterms[1][end]],
               y_2[k-yterms[2][1]], y_2[k-yterms[2][2]],..., y_2[k-yterms[2][end]],
               ...
               y_end[k-yterms[end][1]], y_end[k-yterms[end][2]],..., y_1[k-yterms[end][end]],
               u_1[k-uterms[1][1]], u_1[k-uterms[1][2]],..., u_1[k-uterms[1][end]],
               u_2[k-uterms[2][1]], u_2[k-uterms[2][2]],..., u_2[k-uterms[2][end]],
               ...
               u_end[k-uterms[end][1]], u_end[k-uterms[end][2]],..., u_1[k-yterms[end][end]], Θ)

The following values are internally stored:

| Atributes | Type                  | Brief Description                                     |
|:--------- |:--------------------- |:----------------------------------------------------- |
| mdl       | `StaticModel`         | Static model.                                         |
| yterms    | `Vector{Vector{Int}}` | Contains dependencies on previous terms of `y`.       |
| uterms    | `Vector{Vector{Int}}` | Contains dependencies on previous terms of `u`.       |
| ny        | `Int`                 | Number of inputs `ny = length(yterms)`.               |
| nu        | `Int`                 | Number of outputs `nu = length(uterms)`.              |
| n         | `Int`                 | System order. Maximum value in `yterms`.              |
| iodelay   | `Int`                 | Input-output delay. Minimum value in `uterms`.        |
| m         | `Int`                 | Distance between maximum and minimum `uterms` values. |

The previous  mentioned difference equation can be compactedly
represented using mathematical form:

   y[k] = F(y[k-1],..., y[k-n], u[k-iodelay],..., u[k-iodelay-m+1], Θ)
"""
struct DifferenceEquation{M<:StaticModel}
    mdl::M
    yterms::Vector{Vector{Int}}
    uterms::Vector{Vector{Int}}
    ny::Int
    nu::Int
    n::Int
    m::Int
    iodelay::Int

    # Internal buffers
    x_buf::AbstractVector{Float64}
    z_buf::AbstractVector{Float64}
    dΘ_buf::AbstractMatrix{Float64}
    dx_buf::AbstractMatrix{Float64}
end

function DifferenceEquation(mdl::ParametricModel;
                            yterms::Vector{Vector{Int}}=Vector{Vector{Int}}(0),
                            uterms::Vector{Vector{Int}}=Vector{Vector{Int}}(0))

    # Get number of inputs, outputs terms
    ny = length(yterms)
    nu = length(uterms)

    # Compute order of terms
    n = isempty(yterms) ? 0 : maximum([maximum(z) for z in yterms])
    iodelay = isempty(uterms) ? 0 : minimum([minimum(z) for z in uterms])
    m = isempty(uterms) ? 0 : maximum([maximum(z)-iodelay+1 for z in uterms])

    # Number of input to the static model
    ninputs = sum([length(z) for z in [yterms; uterms]])

    # Check coherence with static model
    if ny != 0 && ny != length_output(mdl)
        throw(ArgumentError("yterms should have $(length_output(mdl)) vectors, one
                             for each output."))
    end
    ny = length_output(mdl)

    if ninputs != length_input(mdl)
        throw(ArgumentError("The total number of inputs to the static model
                             is $ninputs when it should be $(length_input(mdl))."))
    end

    # Check yterms
    for ynterms in yterms
        for j in ynterms
            if j <= 0
                throw(ArgumentError("yterms should not contain non-positive elements."))
            end
        end
    end

    # Initialize buffers
    x_buf = Vector{Float64}(length_input(mdl))
    z_buf = z_buffer(mdl)
    dΘ_buf = dΘ_buffer(mdl)
    dx_buf = dx_buffer(mdl)

    # Call original constructor
    DifferenceEquation(mdl, yterms, uterms, ny, nu, n, m, iodelay,
                       x_buf, z_buf, dΘ_buf, dx_buf)
end

#*************************************************************************#
#
# System Simulation
#
#*************************************************************************#

"""
    simulate(diffeq::DifferenceEquation, u, Θ, y0) -> y

For a given input vector `u` and initial conditions `y0` simulate the system,
returning the output `y`.

Inputs:
------
* `diffeq` is a structure representing a difference equation.
* `u` is a matrix containing the sequence of inputs. It should have
 dimension `nu`-by-`nd`, where `nd` is the number of observations
* `Θ` is the parameter vector, containing `nparams` elements.
* `y0` is a matrix containing the initial conditions. It should have
 dimension `ny`-by-`n`.

Outputs:
-------
* `y` is the sequence of ouputs. It is a matrix with dimension `ny`-by-`nd-m+1`.

The dimension `nd` is infered from the provided vectors.
"""
function simulate{R<:Real,
                  S<:Real,
                  T<:Real}(diffeq::DifferenceEquation, u::AbstractMatrix{T},
                           Θ::AbstractVector{R}, y0::AbstractMatrix{S})

    # Get sizes
    nu, nd = size(u)
    ny = diffeq.ny
    m = diffeq.m

    # Create vector
    y = Matrix{Float64}(ny, nd-m+1)

    # Simulate
    simulate!(diffeq, u, Θ, y0, y)

    # Return the computed vector
    return y
end

"""
    simulate!(diffeq::DifferenceEquation, u, Θ, y0, y[, dΘ, dy0])

For a given input vector `u` and initial conditions `y0` simulate the system
and store the results in the given vector `y`. When the arrays `dΘ`
and `dy0` are provided they are subscribed with the partial derivatives
in relation to `Θ` and `y0`.

Inputs:
------
* `diffeq` is a structure representing a difference equation.
* `u` is a matrix containing the sequence of inputs. It should  have
 dimension `nu`-by-`nd+m-1`, where `nd` is the number of observations
* `Θ` is the parameter vector, containing `nparams` elements.
* `y0` is a matrix containing the initial conditions. It should have
 dimension `ny`-by-`n`.

Outputs:
-------
* `y` is buffer where the computed output will be written. It
should be a matrix with dimension `ny`-by-`nd`, where `nd` is
 the number of observations
* `dΘ`  is buffer where the computed derivatives `dy/dΘ`
 will be written.  It should have dimensions `ny`-by-`nd`-by-`nparams`.
* `dy0` is buffer where the computed derivatives `dy/dy0`
 will be written. It should have dimensions `ny`-by-`nd`-by-`ny`-by-`n`.

If empty arrays are provided for `y`, `dΘ` or `dy0` the respective
computations are bypassed.

The dimension `nd` is infered from the provided vectors.
"""
function simulate!(diffeq::DifferenceEquation,
                   u::AbstractMatrix{<:Real},
                   Θ::AbstractVector{<:Real},
                   y0::AbstractMatrix{<:Real},
                   y::AbstractMatrix{<:Real};
                   dΘ::AbstractArray{<:Real, 3}=Array{Float64}(0, 0, 0),
                   dy0::AbstractArray{<:Real, 4}=Array{Float64}(0, 0, 0, 0))
    # Check Arguments
    argcheck_simulate(diffeq, u, Θ, y0, y, dΘ, dy0)

    # Get Model
    mdl = diffeq.mdl

    # Get model values
    n = diffeq.n
    m = diffeq.m
    iodelay = diffeq.iodelay
    ny = diffeq.ny
    nu = diffeq.nu

    # Get data set size
    nd = size(y, 2)

    # Initialize buffers
    x = diffeq.x_buf
    z = diffeq.z_buf
    if isempty(dΘ)
        dtheta = Matrix{Float64}(0, 0)
    else
        dtheta = diffeq.dΘ_buf
    end

    if isempty(dΘ) && isempty(dy0)
        dx = Matrix{Float64}(0, 0)
    else
        dx = diffeq.dx_buf
    end

    # TODO: vectorize this computation to improve legibility and versatility of the code
    # Compute output sequence
    for k = 1:nd
        # Assemble input vector
        ind = 1
        if n != 0
            for i = 1:ny
                for yterm in diffeq.yterms[i]
                    x[ind] = k-yterm>0 ? y[i, k-yterm]:y0[i, k+n-yterm]
                    ind += 1
                end
            end
        end
        for i = 1:nu
            for uterm in diffeq.uterms[i]
                x[ind] = u[i, k+m-uterm+iodelay-1]
                ind += 1
            end
        end

        # Evaluate input function
        evaluate!(mdl, x, Θ, z, dx=dx, dΘ=dtheta)
        y[:, k] .= z

        # Compute dΘ
        if !isempty(dΘ)
            dΘ[:, k, :] = dtheta
            for i = 1:ny
                ind = 1
                for j = 1:length(diffeq.yterms)
                    for yterm in diffeq.yterms[j]
                        if k-yterm>0
                            dΘ[i, k, :] += dx[i, ind]*dΘ[j, k-yterm, :]
                        end
                        ind += 1
                    end
                end
            end
        end # End compute dΘ

        # Compute dy0
        if !isempty(dy0)
            for i = 1:ny
                for j = 1:ny
                    for p = 1:n
                        ind = 1
                        dy0[i, k, j, p] = 0
                        for l = 1:length(diffeq.yterms)
                            for yterm in diffeq.yterms[l]
                                if k-yterm>0
                                    dy0[i, k, j, p] += dx[i, ind]*dy0[l, k-yterm, j, p]
                                elseif p==k-yterm+n && j==l
                                    dy0[i, k, j, p] += dx[i, ind]
                                end
                                ind += 1
                            end
                        end
                    end
                end
            end
        end # End compute dy0

    end
end

function argcheck_simulate(diffeq::DifferenceEquation,
                           u::AbstractMatrix{<:Real},
                           Θ::AbstractVector{<:Real},
                           y0::AbstractMatrix{<:Real},
                           y::AbstractMatrix{<:Real},
                           dΘ::AbstractArray{<:Real, 3},
                           dy0::AbstractArray{<:Real, 4})

    # Get dimensions
    size_u = size(u)
    nparams = length(Θ)
    size_y0 = size(y0)
    ny, nd = size(y)
    size_dΘ = size(dΘ)
    size_dy0 = size(dy0)
    mdl = diffeq.mdl

    # Argument Check
    if ny != diffeq.ny
        throw(ArgumentError(string("The matrix `y` have $ny rows when it should ",
                                   "have $(diffeq.ny).")))
    end
    if !isempty(u) && size_u != (diffeq.nu, nd+diffeq.m-1)
        throw(ArgumentError(string("The matrix `u` have dimension $(size_u) ",
                                   "when it should have dimension ",
                                   "$((diffeq.nu, nd+diffeq.m-1)).")))
    end
    if nparams != length_params(mdl)
        throw(ArgumentError(string("Parameter vector have $nparams elements when",
                                   "it should have $(length_params(mdl)).")))
    end
    if !isempty(y0) && size_y0 != (ny, diffeq.n)
        throw(ArgumentError(string("The matrix `y0` have dimension $(size_y0) ",
                                   "when it should have dimension ",
                                   "$((ny, diffeq.n)).")))
    end
    if !isempty(dΘ) && size_dΘ != (ny, nd, nparams)
        throw(ArgumentError(string("The array `dΘ` have dimension $(size_dΘ) ",
                                   "when it should have dimension ",
                                   "$((ny, nd, nparams)).")))
    end
    if !isempty(dy0) && size_dy0 != (ny, nd, ny, diffeq.n)
        throw(ArgumentError(string("The array `dy0` have dimension $(size_dy0) ",
                                   "when it should have dimension ",
                                   "$((ny, nd, ny, diffeq.n)).")))
    end
end


text:
using GLFW
using VulkanCore
using VulkanCore.LibVulkan

@assert GLFW.VulkanSupported()

include(joinpath(@__DIR__, "vkhelper.jl"))

const WIDTH = 800
const HEIGHT = 600

enableValidationLayers = true

## init GLFW window
GLFW.WindowHint(GLFW.CLIENT_API, GLFW.NO_API)    # not to create an OpenGL context
GLFW.WindowHint(GLFW.RESIZABLE, 0)
window = GLFW.CreateWindow(WIDTH, HEIGHT, "Vulkan")

## init Vulkan
## creating instance
appInfoRef = VkApplicationInfo(
    "Application Name: Create Instance",
    v"1.0.0",
    "No Engine Name",
    v"1.0.0",
    VK_API_VERSION_1_2,
) |> Ref

extensions = GLFW.GetRequiredInstanceExtensions()
if enableValidationLayers
    push!(extensions, VK_EXT_DEBUG_UTILS_EXTENSION_NAME)  # for debugging, see message callback section
end
@assert check_extensions(extensions)

layers = ["VK_LAYER_KHRONOS_validation"]
@assert check_layers(layers)

## message callback
pfnUserCallback = @cfunction(debug_callback, VkBool32, (VkDebugUtilsMessageSeverityFlagBitsEXT, VkDebugUtilsMessageTypeFlagsEXT, Ptr{VkDebugUtilsMessengerCallbackDataEXT},Ptr{Cvoid}))
debugCreateInfoRef = VkDebugUtilsMessengerCreateInfoEXT(pfnUserCallback, C_NULL) |> Ref

## create instance
createInfoRef = VkInstanceCreateInfo(appInfoRef, debugCreateInfoRef, layers, extensions) |> Ref

instanceRef = Ref(VkInstance(C_NULL))
result = GC.@preserve appInfoRef layers extensions debugCreateInfoRef vkCreateInstance(createInfoRef, C_NULL, instanceRef)
@assert result == VK_SUCCESS "failed to create instance!"

debugMessengerRef = Ref(VkDebugUtilsMessengerEXT(C_NULL))
if CreateDebugUtilsMessengerEXT(instanceRef[], debugCreateInfoRef, C_NULL, debugMessengerRef) != VK_SUCCESS
    @error "failed to set up debug messenger!"
end

## window surface
surface = GLFW.CreateWindowSurface(instanceRef[], window, C_NULL)
surfaceRef = Ref(surface)

## physical device
physicalDeviceRef = Ref(VkPhysicalDevice(VK_NULL_HANDLE))
devices = get_devices(instanceRef[])

# device suitability checks
function is_device_suitable(device::VkPhysicalDevice, surface::VkSurfaceKHR)
    indices = find_queue_families(device, surface)
    return is_complete(indices)
end

for device in devices
    if is_device_suitable(device, surface)
        physicalDeviceRef[] = device
        break
    end
end

if physicalDeviceRef[] == VkPhysicalDevice(VK_NULL_HANDLE)
    @error "failed to find a suitable GPU!"
end

## create logical device
queuePriorityRef = Ref(Cfloat(1.0))
indices = find_queue_families(physicalDeviceRef[], surfaceRef[])
queueCreateInfos = VkDeviceUniqueQueueCreateInfos(indices, queuePriorityRef)

deviceFeaturesRef = Ref(VkPhysicalDeviceFeatures())
logicalCreateInfoRef = VkDeviceCreateInfo(length(queueCreateInfos), queueCreateInfos, deviceFeaturesRef) |> Ref

deviceRef = Ref(VkDevice())
GC.@preserve queuePriorityRef queueCreateInfoRef deviceFeaturesRef begin
    result = vkCreateDevice(physicalDeviceRef[], logicalCreateInfoRef, C_NULL, deviceRef)
end
@assert result == VK_SUCCESS "failed to logical device!"

## retrieving queue handles
graphicsQueueRef = Ref(VkQueue(C_NULL))
vkGetDeviceQueue(deviceRef[], indices.graphicsFamily, 0, graphicsQueueRef)
@assert graphicsQueueRef[] != VkQueue(C_NULL) "failed to retrieve queue handles!"

presentQueueRef = Ref(VkQueue(C_NULL))
vkGetDeviceQueue(deviceRef[], indices.presentFamily, 0, presentQueueRef)
@assert presentQueueRef[] != VkQueue(C_NULL) "failed to retrieve present queue handles!"

## main loop
while !GLFW.WindowShouldClose(window)
    GLFW.PollEvents()
end

## cleaning up
if enableValidationLayers
    DestroyDebugUtilsMessengerEXT(instanceRef[], debugMessengerRef[], C_NULL)
end

vkDestroyDevice(deviceRef[], C_NULL)
vkDestroySurfaceKHR(instanceRef[], surfaceRef[], C_NULL)
vkDestroyInstance(instanceRef[], C_NULL)
GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
"""
    ConditionalMvNormal(m)

Specialization of ConditionalDistribution for `MvNormal`s (for performance with
batches of inputs).  Does the same as ConditionalDistribution(MvNormal,m)
but for batches of inputs a `BatchMvNormal` is constructed that does
not just map over the batch but uses faster matrix multiplications.

```julia-repl
julia> m = SplitLayer(100,[100,100])
julia> p = ConditionalMvNormal(m)
julia> @time rand(p, rand(100,10000);
julia> @time rand(p, rand(100,10000);
julia> @time rand(p, rand(100,10000);
 0.047122 seconds (23 allocations: 38.148 MiB, 24.25% gc time)

julia> p = ConditionalDistribution(MvNormal, m)
julia> @time rand(p, rand(100,10000);
julia> @time rand(p, rand(100,10000);
julia> @time rand(p, rand(100,10000);
 3.626042 seconds (159.97 k allocations: 18.681 GiB, 34.92% gc time)
```

The mapping `m` must return a `Tuple` with mean and variance.
For a convenient way of doing this you can use a `SplitLayer`.

# Examples

`ConditionalMvNormal` and `SplitLayer` together support 3 different variance
configurations: fixed/unit variance, shared variance, and trained variance. The
three different configurations are explained below.

## Fixed/unit variance

Pass a function to the `SplitLayer` that returns the fixed variance with
appropriate batch dimensions
```julia-repl
julia> σ(x::Vector) = 2
julia> σ(x::Matrix) = ones(Float32,size(x,2)) .* 2
julia> m = SplitLayer(Dense(2,3), σ)
julia> p = ConditionalMvNormal(m)
julia> condition(p,rand(Float32,2)) isa DistributionsAD.TuringScalMvNormal
```
Passing a mapping with a single output array assumes unit variance.

## Shared variance

For a learned variance that is the same across the the whole batch, simply pass
a vector (or scalar) to the `SplitLayer`. The `SplitLayer` wraps vectors/scalars
into a `TrainableVector`s/`TrainableScalar`s.
```julia-repl
julia> m = SplitLayer(Dense(2,3), ones(Float32,3))
julia> p = ConditionalMvNormal(m)
julia> condition(p,rand(Float32,2)) isa DistributionsAD.TuringDiagMvNormal
```

## Trained variance

Simply pass another trainable mapping for the variance. By just supplying input
sizes to `SplitLayer` you can automatically create `Dense` layers with given
activation functions. In this example the second activation function makes sure
that the variance is always positive
```julia-repl
julia> m = SplitLayer(2,[3,1],[identity,abs])
julia> p = ConditionalMvNormal(m)
julia> condition(p,rand(Float32,2)) isa DistributionsAD.TuringScalMvNormal
```
"""
struct ConditionalMvNormal{Tm} <: AbstractConditionalDistribution
    mapping::Tm
end

function condition(p::ConditionalMvNormal, z::AbstractVector)
    (μ,σ) = mean_var(p.mapping(z))
    if length(σ) == 1
        σ = σ[1]
    end
    DistributionsAD.TuringMvNormal(μ,σ)  # for CuArrays/gradients
end

function condition(p::ConditionalMvNormal, z::AbstractMatrix)
    (μ,σ) = mean_var(p.mapping(z))
    if size(σ,1) == 1
        σ = dropdims(σ, dims=1)
    end
    BatchMvNormal(μ,σ)
end

# dispatches for different outputs from mappings
# general case
mean_var(x::Tuple) = x
# single output assumes σ=1
mean_var(x::AbstractVector) = (x, 1)
mean_var(x::AbstractMatrix) = (x, fillsimilar(x,size(x,2),1))
# fixed scalar variance
# mean_var(x::Tuple{<:AbstractVector,<:Real}) = x; is already coverged
function mean_var(x::Tuple{<:AbstractMatrix,<:Real})
    (μ,σ) = x
    (μ,fillsimilar(μ,size(μ,2),σ))
end

@functor ConditionalMvNormal


text:
abstract type Disturbances{T} end

"""
    empty disturbances
"""
struct EmptyDisturbances{T} <: Disturbances{T}
    w::Vector{T}
end

function empty_disturbances(model)
    EmptyDisturbances(zeros(model.nw))
end

function disturbances(d::EmptyDisturbances, x, t)
    return d.w
end


text:
using Test
using LinearAlgebra
using Random
using KrylovKit

const n = 10
const N = 100

const η₀   = 0.75 # seems to be necessary to get sufficient convergence for GKL iteration with Float32 precision
const cgs = ClassicalGramSchmidt()
const mgs = ModifiedGramSchmidt()
const cgs2 = ClassicalGramSchmidt2()
const mgs2 = ModifiedGramSchmidt2()
const cgsr = ClassicalGramSchmidtIR(η₀)
const mgsr = ModifiedGramSchmidtIR(η₀)

Random.seed!(76543210)

include("linalg.jl")

module PureVecs
    using Test
    using LinearAlgebra
    using Random
    using KrylovKit

    include("setcomparison.jl")

    const n = 10
    const N = 100

    const η₀   = 0.75 # seems to be necessary to get sufficient convergence for GKL iteration with Float32 precision
    const cgs = ClassicalGramSchmidt()
    const mgs = ModifiedGramSchmidt()
    const cgs2 = ClassicalGramSchmidt2()
    const mgs2 = ModifiedGramSchmidt2()
    const cgsr = ClassicalGramSchmidtIR(η₀)
    const mgsr = ModifiedGramSchmidtIR(η₀)

    wrapvec(v) = v
    unwrapvec(v) = v
    wrapop(A::AbstractMatrix) = A

    t = time()
    include("factorize.jl")
    include("linsolve.jl")
    include("eigsolve.jl")
    include("schursolve.jl")
    include("svdsolve.jl")
    include("expintegrator.jl")
    t = time() - t
    println("Julia Vector type: tests finisthed in $t seconds")
end

module MinimalVecs
    using Test
    using LinearAlgebra
    using Random
    using KrylovKit

    include("setcomparison.jl")

    const n = 10
    const N = 100

    const η₀   = 0.75 # seems to be necessary to get sufficient convergence for GKL iteration with Float32 precision
    const cgs = ClassicalGramSchmidt()
    const mgs = ModifiedGramSchmidt()
    const cgs2 = ClassicalGramSchmidt2()
    const mgs2 = ModifiedGramSchmidt2()
    const cgsr = ClassicalGramSchmidtIR(η₀)
    const mgsr = ModifiedGramSchmidtIR(η₀)

    include("minimalvec.jl")

    wrapvec(v) = MinimalVec(v)
    unwrapvec(v::MinimalVec) = getindex(v)
    wrapop(A::AbstractMatrix) =
        (v, flag=false)->wrapvec(flag ? A'*unwrapvec(v) : A*unwrapvec(v))

    t = time()
    include("factorize.jl")
    include("linsolve.jl")
    include("eigsolve.jl")
    include("schursolve.jl")
    include("svdsolve.jl")
    include("expintegrator.jl")
    t = time() - t
    println("Minimal vector type: tests finisthed in $t seconds")
end

include("recursivevec.jl")


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--opt1"
        help = "an option with an argument"
        "--opt2", "-o"
        help = "another option with an argument"
        arg_type = Int
        default = 0
        "--flag1"
        help = "an option without argument, i.e. a flag"
        action = :store_true
        "arg1"
        help = "a positional argument"
        required = true
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg, val) in parsed_args
        println("  $arg  =>  $val")
    end
end

main()


text:
# this is taken from Louis Dressel's POMDPs.jl. Thanks Louis!

# Abstract Type representing alpha vectors
abstract type Alphas end

mutable struct POMDPAlphas <: Alphas
    alpha_vectors::Matrix{Float64}
    alpha_actions::Vector{Int64}

    POMDPAlphas(av::Matrix{Float64}, aa::Vector{Int64}) = new(av, aa)

    # Constructor if no action list is given
    # Here, we 0-index actions, to match sarsop output
    function POMDPAlphas(av::Matrix{Float64})
        numActions = size(av, 1)
        alist = [0:(numActions-1)]
        return new(av, alist)
    end

    # Constructor reading policy from file
    function POMDPAlphas(filename::AbstractString)
        alpha_vectors, alpha_actions = read_pomdp(filename)
        return new(alpha_vectors, alpha_actions)
    end

    # Default constructor
    function POMDPAlphas()
        return new(zeros(0,0), zeros(Int64,0))
    end
end


text:
"""
Base module for Yao.
"""
module YaoBase

using LinearAlgebra, LuxurySparse, SparseArrays, Random

include("utils/ast_tools.jl")

include("utils/constants.jl")
include("utils/math.jl")
include("utils/interface.jl")

include("error.jl")
include("abstract_register.jl")
include("adjoint_register.jl")

include("inspect.jl")
include("instruct.jl")

# compat with older version of dependencies
include("compat.jl")

# deprecation warns
include("deprecations.jl")

end # module


text:
struct Foo{A,B,C}
    a::A
    b::B
    c::C
end

@info Foo(1,2,3)

function create_foo(a, b, c)
    Foo(a, b, c)
end

function create_foo(a, b, c, d::String)
    Foo(a[:c], b, c)
end

create_foo(Dict(:a=>100), 10, 1, "test")

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl