text:
function pythonprint(
    emod :: ERA5Hourly,
    evar :: SingleVariable,
    ereg :: ERA5Region,
)

end

function pythonprint(
    emod :: ERA5Monthly,
    evar :: SingleVariable,
    ereg :: ERA5Region,
)

end

function pythonprint(
    emod :: ERA5MonthlyHour,
    evar :: SingleVariable,
    ereg :: ERA5Region,
)

end

function pythonprint(
    emod :: ERA5Hourly,
    evar :: PressureVariable,
    ereg :: ERA5Region,
)

end

function pythonprint(
    emod :: ERA5Monthly,
    evar :: PressureVariable,
    ereg :: ERA5Region,
)

end

function pythonprint(
    emod :: ERA5MonthlyHour,
    evar :: PressureVariable,
    ereg :: ERA5Region,
)

end


text:
using Tmatrix
using ProfileView

# inputs
rx, rz = 1e-6, 1.3e-6
n_max = 2
k1_r = 1e7; k1_i = 0.0
k2_r = 1.5e7; k2_i = 1e3

@profiler calculate_Tmatrix_for_spheroid(rx, rz, n_max, Complex(k1_r, k1_i), Complex(k2_r, k2_i))
@profiler Tmatrix.calculate_Tmatrix_for_spheroid_SeparateRealImag(rx, rz, n_max, k1_r, k1_i, k2_r, k2_i)

ext: jl


text:
# Slack `chat.postMessage` recommends keeping the character length around 4000:
# https://api.slack.com/methods/chat.postMessage#formatting
const MAX_MESSAGE_LENGTH = 3800

"""
    message_count_suffix(i, n)

If `n==1` returns `""`, otherwise `"\n\n[\$i / \$n]"`.
"""
function message_count_suffix(i, n)
    return n == 1 ? "" : "\n\n[$i / $n]"
end

"""
    combine_texts(texts; max_length=$MAX_MESSAGE_LENGTH,
                  message_count_suffix=message_count_suffix) -> Vector{String}

Combine `N` strings in `texts` into `M <= N` strings, by concatenating consecutive strings, such that
each of the resulting strings has length `<= max_length`, unless the input string was already
longer than `max_length` , in which case it is passed on.

The purpose of this is for nicely splitting long lists of attachments.

The keyword argument `message_count_suffix` defaults to [`SlackThreads.message_count_suffix`](@ref), and is used
for formatting the suffix of messages in the case that a message needs to split into multiple messages.
"""
function combine_texts(texts; max_length=MAX_MESSAGE_LENGTH,
                       message_count_suffix=message_count_suffix)
    messages = String[]
    current_message = ""
    current_length = 0
    for text in texts
        new_length = length(text)
        if current_length + new_length > max_length
            # Finish current message; start new one
            if !isempty(current_message)
                push!(messages, current_message)
            end
            current_message = text
            current_length = new_length
        else
            # Add to current message
            current_message *= text
            current_length += new_length
        end
    end
    # Add last message
    if !isempty(current_message)
        push!(messages, current_message)
    end

    # Add [i / n] tags
    n = length(messages)
    for i in eachindex(messages)
        messages[i] *= message_count_suffix(i, n)
    end
    return messages
end


text:
using NPZ
using JuMP
using Ipopt
using LinearAlgebra
using WoodburyMatrices
using .Threads
using Juniper
using Cbc
using Gurobi

X1 = npzread("Kernels/G1.npy")
X2 = npzread("Kernels/G2.npy")
X3 = npzread("Kernels/G3.npy")

lnorm = -log(ϵ)*(2*size(X1,2))

function d_oed(X1::Matrix{T}, X2::Matrix{T}, X3::Matrix{T}, B, ϵ, λ, lnorm; maxiter=1000) where T <: Number
    X = [X3 zeros(T, size(X3));
         zeros(T, size(X3)) X3;
         X1 X2]
    #Xn = maximum(sqrt.(sum(X.^2, dims=2)))
    Xn = sqrt.(sum(X.^2, dims=2)) #normalize all rows to be 1?
    X ./= Xn

    n = size(X3, 1)
    ndas = size(X1, 1)

    function Ml(μ, X, ϵ)
        M = ϵ*I + X'*Diagonal(μ)*X
        Mc = cholesky(Hermitian(M))
        l = μ.^2 .*diag(X*(Mc\X'))
        return (M, Mc, l)
    end

    Ml(μ) = Ml(μ, X, ϵ)

    function f(μ, λ, lnorm)
        M, Mc, l = Ml(μ)
        return -logdet(M) + λ*maximum(l) - lnorm
    end

    f(μ) = f(μ, λ, lnorm)

    grad_f(y) = f(vcat([y[i] for i in 1:length(y)], [y[i] for i in 1:length(y)], ones(T, ndas)))

    jump_wrapper_f(x...) = f(vcat([x[i] for i in 1:length(x)], [x[i] for i in 1:length(x)], ones(T, ndas)))

    function gi(i, Mc, μ, u, lmind, X, λ, T)
        ld_term = @views -X[i,:]'*(Mc\X[i,:])
        tmp_A = SymWoodbury(Mc, X[i,:], -μ[i])
        tmp_b = @views tmp_A \ X[i,:]
        tmp_1 = (u'*tmp_b)^2
        tmp_d = @views X[i,:]'*tmp_b
        tmp_den = (1+μ[i]*tmp_d)^2
        lev_term_1 = -μ[lmind]^2*((1+μ[i]*tmp_d)*tmp_1 - μ[i]*tmp_1*tmp_d) / tmp_den
        lev_term_2 = (i==lmind) ? 2*μ[lmind]*(u'*(Mc\u)) : zero(T)
        return ld_term + λ*(lev_term_1 + lev_term_2)
    end

    gi(i, Mc, μ, u, lmind) = gi(i, Mc, μ, u, lmind, X, λ, T)

    function g!(gvec, μ, X, n)
        M, Mc, l = Ml(μ)
        lm, lmind = findmax(l)
        lmind=lmind[1]
        u = @views X[lmind, :]
        @threads for i = 1:n
            g1 = gi(i, Mc, μ, u, lmind) # in respect to x component
            g2 = gi(i+n, Mc, μ, u, lmind) # in respect to y component
            gvec[i] = g1+g2
        end
    end

    g!(gvec, μ) = g!(gvec, μ, X, n)

    jump_wrapper_g!(gvec, x...) = g!(gvec, vcat([x[i] for i in 1:length(x)], [x[i] for i in 1:length(x)], ones(T, ndas)))

    # optimizer = Juniper.Optimizer
    # nl_solver= optimizer_with_attributes(Ipopt.Optimizer, "print_level"=>1)
    # mip_solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1)
    # model = Model(optimizer_with_attributes(optimizer, "nl_solver"=>nl_solver, "mip_solver"=>mip_solver, "log_levels"=>[:Table, :Info], "time_limit"=>10*3600, "feasibility_pump_time_limit"=>300,"strong_branching_time_limit"=>300, "traverse_strategy"=>:BFS,"branch_strategy"=>:Reliability,
    # "registered_functions" => [
    #     Juniper.register(:f, n, jump_wrapper_f, jump_wrapper_g!)
    # ]))
    model = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter"=>maxiter, "tol"=>1e-2))
    JuMP.register(model, :f, n, jump_wrapper_f, jump_wrapper_g!)

    @variable(model, 0.0<=x[1:n]<= 1.0)
    #@variable(model, x[1:n], binary=true)

    @NLconstraint(model, +(x...) <= B)

    @NLobjective(model, Min, f(x...))

    optimize!(model)

    return (model, x)
end

ϵ = 1e-3
lnorm = -log(ϵ)*(2*size(X1,2))
model, x = d_oed(X1, X2, X3, 32, ϵ, 1/ϵ, lnorm, maxiter=1000)
println("Model Objective = $(objective_value(model))")
npzwrite("res_mixed.npy", value.(x))
# ;

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
export AbstractForwardSensFspMatrixSparse, ForwardSensFspMatrixSparse, matvec!

abstract type AbstractForwardSensFspMatrixSparse <: ForwardSensFspMatrix end

"""
The time-varying matrix that defines the right-hand side of the truncated forward sensitivity Chemical Master Equation.
This type works with `StateSpaceSparse` subtypes and it uses `SparseMatrixCSC` type to internally represent component matrices.
"""
mutable struct ForwardSensFspMatrixSparse{NS,NR,IntT<:Integer,RealT<:AbstractFloat} <: AbstractForwardSensFspMatrixSparse
    fspmatrix::FspMatrixSparse{NS,NR,IntT,RealT}

    propensity_gradients::Vector{<:PropensityGradient}

    timeinvariant_matdiffs::Vector{SparseMatrixCSC}

    separabletv_matdiff_sparsity_pattern::SparseMatrixCSC{Bool}
    separabletv_nzmatdiffs::Vector{SparseMatrixCSC}
    separabletv_mat_map::Vector{IntT}

    jointtv_matdiff_sparsity_pattern::SparseMatrixCSC{Bool}
    jointtv_nzmatdiffs::Vector{SparseMatrixCSC}
end
get_propensity_gradients(SA::ForwardSensFspMatrixSparse) = SA.propensity_gradients
get_timeinvariant_matdiffs(SA::ForwardSensFspMatrixSparse) = SA.timeinvariant_matdiffs
get_separabletv_sparsity_pattern(SA::ForwardSensFspMatrixSparse) = SA.separabletv_matdiff_sparsity_pattern
get_separabletv_nzmatdiffs(SA::ForwardSensFspMatrixSparse) = SA.separabletv_nzmatdiffs
get_separabletv_mat_map(SA::ForwardSensFspMatrixSparse) = SA.separabletv_mat_map
get_jointtv_sparsity_pattern(SA::ForwardSensFspMatrixSparse) = SA.jointtv_matdiff_sparsity_pattern
get_jointtv_nzmatdiffs(SA::ForwardSensFspMatrixSparse) = SA.jointtv_nzmatdiffs

function ForwardSensFspMatrixSparse{RealT}(model::CmeModelWithSensitivity{IntT}, statespace::StateSpaceSparse{NS,NR,IntT,SizeT}) where {NS,NR,IntT<:Integer,SizeT<:Integer,RealT<:AbstractFloat}
    parameter_count = get_parameter_count(model)
    grad_sparsity_patterns = get_gradient_sparsity_patterns(model)
    propensity_gradients = get_propensity_gradients(model)
    parameters = get_parameters(model)

    fspmatrix = FspMatrixSparse{RealT}(statespace, get_propensities(model), parameters = get_parameters(model))
    timeinvariant_ids = get_timeinvariant_propensity_ids(fspmatrix)
    jointtv_ids = get_jointtv_propensity_ids(fspmatrix)
    separabletv_ids = get_separabletv_propensity_ids(fspmatrix)

    n = get_rowcount(fspmatrix)
    timeinvariant_matdiffs = []
    tinvariant_sparsity = grad_sparsity_patterns[timeinvariant_ids, :]
    rowvals = SparseArrays.getrowval(tinvariant_sparsity)
    for ip in 1:parameter_count
        rids = Vector{SizeT}()
        cids = Vector{SizeT}()
        vals = Vector{RealT}()
        for i in nzrange(tinvariant_sparsity, ip)
            r = timeinvariant_ids[rowvals[i]]
            Is, Js, Vs = _generate_sparsematrix_entries(statespace, get_single_pardiff(propensity_gradients[r], ip), r, parameters, valtype = RealT)
            append!(rids, Is)
            append!(cids, Js)
            append!(vals, Vs)
        end
        push!(timeinvariant_matdiffs, sparse(rids, cids, vals, n, n))
    end

    separabletv_sparsity = grad_sparsity_patterns[separabletv_ids, :]
    separabletv_nzmatdiffs = []
    separabletv_mat_map = []
    rowvals = SparseArrays.getrowval(separabletv_sparsity)
    for ip in 1:parameter_count
        for i in nzrange(separabletv_sparsity, ip)
            r = separabletv_ids[rowvals[i]]
            Is, Js, Vs = _generate_sparsematrix_entries(statespace, get_single_statefactor_pardiff(propensity_gradients[r], ip), r, parameters, valtype = RealT)
            push!(separabletv_nzmatdiffs, sparse(Is, Js, Vs, n, n))
            push!(separabletv_mat_map, findfirst(separabletv_ids .== r))
        end
    end

    jointtv_sparsity = grad_sparsity_patterns[jointtv_ids, :]
    jointtv_nzmatdiffs = []
    rowvals = SparseArrays.getrowval(jointtv_sparsity)
    for ip in 1:parameter_count
        for i in nzrange(jointtv_sparsity, ip)
            r = jointtv_ids[rowvals[i]]
            Is, Js, Vs = _generate_sparsematrix_entries(statespace, nothing, r, parameters, valtype = RealT)
            push!(jointtv_nzmatdiffs, sparse(Is, Js, Vs, n, n))
        end
    end

    return ForwardSensFspMatrixSparse{NS,NR,IntT,RealT}(
        fspmatrix,
        propensity_gradients,
        timeinvariant_matdiffs,
        separabletv_sparsity,
        separabletv_nzmatdiffs,
        separabletv_mat_map,
        jointtv_sparsity,
        jointtv_nzmatdiffs
    )
end

function matvec!(out::AbstractVector{RealT}, t::AbstractFloat, SA::ForwardSensFspMatrixSparse, vs::AbstractVector{RealT}) where {RealT<:AbstractFloat}
    A = SA.fspmatrix
    θ = get_parameters(A)
    parameter_count = length(θ)
    propensities = get_propensities(A)
    propensity_gradients = get_propensity_gradients(SA)

    jointtv_ids = get_jointtv_propensity_ids(A)
    separabletv_ids = get_separabletv_propensity_ids(A)

    ti_matdiffs = get_timeinvariant_matdiffs(SA)

    septv_pattern = get_separabletv_sparsity_pattern(SA)
    septv_map = get_separabletv_mat_map(SA)
    septv_matdiffs = get_separabletv_nzmatdiffs(SA)

    jointtv_pattern = get_jointtv_sparsity_pattern(SA)
    jointtv_matdiffs = get_jointtv_nzmatdiffs(SA)

    n = get_rowcount(A)
    matvec!(view(out, 1:n), t, A, view(vs, 1:n))
    for ip in 1:parameter_count
        outview = view(out, ip*n+1:(ip+1)*n)
        vsview = view(vs, ip*n+1:(ip+1)*n)
        matvec!(outview, t, A, vsview)

        mul!(outview, ti_matdiffs[ip], view(vs, 1:n), 1, 1)
        for i in nzrange(septv_pattern, ip)
            r = separabletv_ids[SparseArrays.getrowval(septv_pattern)[i]]

            mul!(outview, septv_matdiffs[i], view(vs, 1:n), propensities[r].tfactor(t, θ), 1.0)

            j = septv_map[i]
            tfactordiff = get_single_tfactor_pardiff(propensity_gradients[r], ip)
            mul!(outview, A.separabletv_factormatrices[j], view(vs, 1:n), tfactordiff(t, θ), 1.0)
        end

        for i in nzrange(jointtv_pattern, ip)
            r = jointtv_ids[SparseArrays.getrowval(jointtv_pattern)[i]]
            pdiff = get_single_pardiff(propensity_gradients[r], ip)
            _update_sparsematrix!(jointtv_matdiffs[i], A.states, pdiff, t, θ)
            mul!(outview, jointtv_matdiffs[i], view(vs, 1:n), 1, 1)
        end
    end
    return nothing
end


text:
using RubiSymbolics: syntax_repl_dict, function_repl_dict, @rule, Prewalk, PassThrough, identifier, FreeQ
import Metatheory.Syntax: makepattern
#--- global variables
rubi_rules_path = joinpath(@__DIR__, "../Rubi/IntegrationRules")
rules_root_path = joinpath(@__DIR__, "rules")
rule_re = r"(?<=\n)(Antiderivative\(.*,\h*~x\)\h*:=\h*.*?)(?=\n)"
#--- Rules
# r_cond = @rule (~a <-- ~b) --> (~a where ~b)    # keep as <-- in *.jl file ; 'where' has some issues at printing
#--- utility Functions
stem(path) = splitext(path)[1]
function rule_subs(m)
    try
        expr = Meta.parse(m)
        @assert expr.head == :(:=) && length(expr.args) == 2
        # make rhs using lhs variables
        lhs_pvars = Symbol[]
        makepattern(expr.args[1], lhs_pvars, [])    # fills lhs_pvars
        rhs_str = replace(string(expr.args[2]),
            [identifier(sym) => "(~"*string(sym)*")" for sym in lhs_pvars]...)
        rhs = Meta.parse(rhs_str)
        expr = Expr(:call, :(=>), expr.args[1], rhs)
        # expr = Prewalk(PassThrough(r_cond))(expr)
        return "@apply_utils "*string(expr)
    catch e
        printstyled("Parse error:", e, "\n", color=:yellow)
        println(m)
        return "# Missing due to parse error: "*m
    end
end

# rule_subs(m) = m

for (root, dirs, files) in walkdir(rubi_rules_path)
    rel_path = relpath(root, rubi_rules_path)
    mkpath(joinpath(rules_root_path, rel_path))
    root_include_file = joinpath(rules_root_path, rel_path, ".jl")
    touch(root_include_file)
    open(root_include_file, "w") do f
        for dir in dirs
            write(f, "include(\"$(joinpath(dir, ".jl"))\")\n")
        end
        for file in files
            write(f, "include(\"$(stem(file)*".jl")\")\n")
        end
    end
    for file in files
        rules_file = joinpath(rules_root_path, rel_path, stem(file)*".jl")
        touch(rules_file)
        open(rules_file, "w") do f
            orig_file_content = open(joinpath(root, file)) do rf
                read(rf, String)
            end
            out_file_content = replace(orig_file_content,
                "\r" => "",  # remove windows line end
                r"(\w[a-zA-Z0-9]*)_\."  => s"(~\1')",          # optional symbols
                r"(\w[a-zA-Z0-9]*)_(?![a-zA-Z])" => s"(~\1)",
                "x_Symbol" => "~x",
                syntax_repl_dict...,
                function_repl_dict...
            )
            out_file_content = replace(out_file_content,  # TODO: fix before
                "=#*" => "=#",
                "*#=" => "#=",)
            out_file_content = replace(out_file_content, rule_re => rule_subs)
            out_file_content = replace(out_file_content,
                r"(\n(?!\n|$))"   => "\n"*" "^4,    # indent
                r"^(?!\n)"      => " "^4,
                )
            th_number = replace(split(file)[1], "." => "_")
            write(f, "int_rules_$th_number = @theory begin\n")
            write(f, out_file_content)
            write(f, "end\n")
        end
    end
end


text:
@testset "1.1.1 Linear" begin
    include("1.1.1.2 (a+b x)^m (c+d x)^n.jl")
    include("1.1.1.3 (a+b x)^m (c+d x)^n (e+f x)^p.jl")
    include("1.1.1.4 (a+b x)^m (c+d x)^n (e+f x)^p (g+h x)^q.jl")
    include("1.1.1.5 P(x) (a+b x)^m (c+d x)^n.jl")
    include("1.1.1.6 P(x) (a+b x)^m (c+d x)^n (e+f x)^p.jl")
    include("1.1.1.7 P(x) (a+b x)^m (c+d x)^n (e+f x)^p (g+h x)^q.jl")
end

max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/.jl


text:
__precompile__()

module VCmicrobiome

using ConjugatePriors
using Rmath
using StatsBase
using Distributions

export

      vctest,
      vctestnullsim,
      microvctest

include("vctest.jl")
include("vctestnullsim.jl")
include("microvctest.jl")
# include("MultiKernelInput.jl")
# include("MMicrobiomeVCTest.jl")
end # module


text:
using Base.Threads
using DataStructures
import Base: print, show

"""
Gaussian likelihood Graph Fused Elastic Net Model
"""
mutable struct GaussianGFEN
    # trail info
    ptr::Vector{Int}
    brks::Vector{Int}
    lambdasl1::Vector{Float64}
    lambdasl2::Vector{Float64}
    num_nodes::Int
    num_slacks::Int
    num_trails::Int

    # prior info
    prior_mean::Float64
    prior_weight::Float64

    # model parameters and computed quantities
    beta::Vector{Float64}
    slack::Vector{Float64}
    dual::Vector{Float64}

    # convergence and steps
    trained::Bool
    steps::Int
    maxsteps::Int
    converged::Bool
    check_convergence_every::Int
    abstol::Float64
    reltol::Float64
    save_norms::Bool
    save_loss::Bool
    prim_norms::Vector{Float64}
    dual_norms::Vector{Float64}
    loss::Vector{Float64}

    # admm associated parameters
    admm_penalty::Float64
    admm_init_penalty::Float64
    admm_residual_balancing::Bool
    admm_residual_balancing_gap::Float64
    admm_adaptive_inflation::Bool
    admm_balance_every::Int
    admm_fixed_inflation_factor::Float64
    admm_penalty_max_inflation::Float64
    admm_min_penalty::Float64

    # Acceleration parameters
    accelerate::Bool
    accelerate_min_improvement::Float64
    accelerate_use_restarts::Bool

    # main constructor
    function GaussianGFEN(
        ptr::Vector{Int},
        brks::Vector{Int},
        lambdasl1::Vector{Float64},
        lambdasl2::Vector{Float64};
        prior_mean::Float64 = 0.0,
        prior_weight::Float64 = 1e-6,
        maxsteps::Int = Int(1e5),
        check_convergence_every::Int = 1,
        abstol::Float64 = 1e-9,
        reltol::Float64 = 1e-3,
        save_norms::Bool = true,
        save_loss::Bool = true,
        admm_init_penalty::Float64 = 1.0,
        admm_min_penalty::Float64 = 0.01,
        admm_fixed_inflation_factor::Float64 = 2.,
        admm_balance_every::Int = 10,
        admm_residual_balancing_gap::Float64 = 10.,
        admm_penalty_max_inflation::Float64 = 100.,
        admm_residual_balancing::Bool = true,
        admm_adaptive_inflation::Bool = true,
        accelerate::Bool = false,
        accelerate_min_improvement::Float64 = 0.999,
        accelerate_use_restarts::Bool = true
    )
        # Model Validation
        @assert maxsteps > 0
        @assert abstol >= 0.
        @assert reltol > 0.
        @assert admm_init_penalty > 0.
        @assert admm_fixed_inflation_factor > 0.
        @assert admm_balance_every > 0
        @assert admm_residual_balancing_gap > 0.
        @assert check_convergence_every > 0
        @assert admm_penalty_max_inflation > 0.
        @assert prior_weight > 0.

        # new object
        x = new()

        # assign trails
        x.ptr = ptr
        x.brks = brks
        x.lambdasl1 = lambdasl1
        x.lambdasl2 = lambdasl2
        x.num_nodes = maximum(ptr)
        x.num_trails = length(brks) - 1
        x.num_slacks = length(ptr)
        x.prior_mean = prior_mean
        x.prior_weight = prior_weight

        # set values for variables to be trained
        x.beta = zeros(x.num_nodes)
        x.slack = zeros(x.num_slacks)
        x.dual = zeros(x.num_slacks)

        # convergence and steps
        x.trained = false
        x.steps = 0
        x.maxsteps = maxsteps
        x.converged  = false
        x.check_convergence_every = check_convergence_every
        x.abstol = abstol
        x.reltol = reltol
        x.save_norms = save_norms
        x.save_loss = save_loss
        x.prim_norms = Float64[]
        x.dual_norms = Float64[]
        x.loss = Float64[]

        # admm associated parameters
        x.admm_penalty = admm_init_penalty
        x.admm_init_penalty = admm_init_penalty
        x.admm_residual_balancing = admm_residual_balancing
        x.admm_residual_balancing_gap = admm_residual_balancing_gap
        x.admm_adaptive_inflation = admm_adaptive_inflation
        x.admm_balance_every = admm_balance_every
        x.admm_fixed_inflation_factor = admm_fixed_inflation_factor
        x.admm_penalty_max_inflation = admm_penalty_max_inflation
        x.admm_min_penalty = admm_min_penalty

        # Acceleration parameters
        x.accelerate = accelerate
        x.accelerate_min_improvement = accelerate_min_improvement
        x.accelerate_use_restarts = accelerate_use_restarts

        # end constructor
        return x
    end
end
Base.show(io::IO, x::GaussianGFEN) = print(io, "GaussianGFEN")
Base.print(io::IO, x::GaussianGFEN) = print(io, "GaussianGFEN")

function init_vars(
    model::GaussianGFEN, # model
    y::AbstractVector{<:Union{Float64, Missing}}
)
    ptr = model.ptr
    brks = model.brks
    α = model.admm_init_penalty
    edgewts = model.lambdasl1 / α
    edgewts2 = model.lambdasl2 / α

    # hot start beta
    β = zeros(model.num_nodes)
    # k = findfirst(x -> !ismissing(x), y)
    # curr::Float64 = y[k]
    # @simd for i in k:model.num_nodes
    #     if !ismissing(y[i])
    #         curr = y[i]
    #     end
    #     if !ismissing(curr)
    #         β[i] = curr
    #     end
    # end
    # @simd for i in reverse(1:model.num_nodes)
    #     if !ismissing(y[i])
    #         curr = y[i]
    #     end
    #     if !ismissing(curr)
    #         β[i] = curr
    #     end
    # end

    # Initiate other ADMM variables
    z = zeros(model.num_slacks) # slacks
    z2 = zeros(model.num_slacks) # slacks
    u = zeros(model.num_slacks) # scaled dual
    u2 = zeros(model.num_slacks) # scaled dual
    Δz = zeros(model.num_slacks) # used for momentum and convergence
    Δu = zeros(model.num_slacks) # used for momentum
    Δz2 = zeros(model.num_slacks) # used for momentum
    Δu2 = zeros(model.num_slacks) # used for momentum

    # initialize beta

    # trail visits to each node
    num_visits = zeros(Int, model.num_nodes)
    @simd for j in 1:model.num_slacks
        z[j] = β[ptr[j]] # best guess for z
        z2[j] = β[ptr[j]] # best guess for z
        num_visits[ptr[j]] += 1
    end

    εabs_p = sqrt(2.0 * model.num_slacks) * model.abstol # max dual error for the absolute convergence
    εabs_d = sqrt(model.num_nodes) * model.abstol # max prim error for the absolute convergence

    β, z, u, Δz, Δu, z2, u2, Δz2, Δu2, α, edgewts, edgewts2, num_visits, εabs_p, εabs_d
end

# sorted(i::Int, j::Int) = (min(i, j), max(i, j))
# function cnt_reps(ptr::Vector{Int}, brks::Vector{Int})
#     edges = [sorted(ptr[i], ptr[i+1]) for i in 1:length(ptr) if (i + 1) ∉ brks]
#     ctr = counter(edges)
#     edge_cnt = [ctr[sorted(ptr[i], ptr[i+1])] for i in 1:length(ptr) - 1]
#     return edge_cntmodel.
# end

function update_primal!(
    model::GaussianGFEN,
    β::Vector{Float64},
    z::Vector{Float64},
    u::Vector{Float64},
    z2::Vector{Float64},
    u2::Vector{Float64},
    num_visits::Vector{Int},
    y::AbstractVector{<:Union{Float64, Missing}},
    α::Float64
)
    ptr, brks = model.ptr, model.brks

    # admm pseudovalue for beta
    r = zeros(model.num_nodes)
    @simd for j = 1:model.num_slacks # can't do simd here :(
        r[ptr[j]] += z[j] - u[j] + z2[j] - u2[j]
    end

    @simd for i = 1:model.num_nodes
        if !ismissing(y[i])
            H = 1.0 + 2.0 * α * num_visits[i] + model.prior_weight
            β[i] = (y[i] + α * r[i] + model.prior_mean) / H # this is very important for stability!
        else
            H = 2.0 * α * num_visits[i] + model.prior_weight
            β[i] = (α * r[i] + model.prior_mean) / H # this is very important for stability!
        end
    end
end

@inline function update_slack_and_dual!(
    model::GaussianGFEN,
    z::Vector{Float64},
    u::Vector{Float64},
    Δz::Vector{Float64},
    Δu::Vector{Float64},
    z2::Vector{Float64},
    u2::Vector{Float64},
    Δz2::Vector{Float64},
    Δu2::Vector{Float64},
    β::Vector{Float64},
    edgewts::Vector{Float64},
    edgewts2::Vector{Float64},
    parallel::Bool
)
    ptr = model.ptr
    brks = model.brks

    @simd for j in 1:model.num_slacks
        Δu[j] = -u[j] # we'll add u(k+1) in last loop
        Δz[j] = -z[j] # we'll add z(k+1) in last loop
        Δu2[j] = -u2[j] # we'll add u2(k+1) in last loop
        Δz2[j] = -z2[j] # we'll add z2(k+1) in last loop
        u[j] += β[ptr[j]] # u(k+1) += beta(k+1) - z(k+1), z is subtracted in last loop
        u2[j] += β[ptr[j]] # u2(k+1) += beta(k+1) - z2(k+1), z2 is subtracted in last loop
    end

    zbuff = u # at this point u is u(k) + β(k) because of previous loop
    zbuff2 = u2 # at this point u is u(k) + β(k) because of previous loop
    if parallel
        @threads for t = 1:model.num_trails
            from = brks[t]
            to = brks[t + 1] - 1
            filter1D!(z, zbuff, edgewts, from, to)
            filter1Dl2!(z2, zbuff2, edgewts2, from, to)
        end
    else
        for t = 1:model.num_trails
            from = brks[t]
            to = brks[t + 1] - 1
            filter1D!(z, zbuff, edgewts, from, to)
            filter1Dl2!(z2, zbuff2, edgewts2, from, to)
        end
    end

    @simd for j in 1:model.num_slacks
        u[j] -= z[j] # now u(k+1) += beta(k+1) - z(k+1)
        Δz[j] += z[j] # now Δz(k+1) = z(k+1) - z(k)
        Δu[j] += u[j] # now Δu(k+1) = u(k+1) - u(k)
        u2[j] -= z2[j] # now u2(k+1) += beta(k+1) - z2(k+1)
        Δz2[j] += z2[j] # now Δz2(k+1) = z2(k+1) - z2(k)
        Δu2[j] += u2[j] # now Δu2(k+1) = u2(k+1) - u2(k)
    end
end

@inline function compute_negll(
    model::GaussianGFEN,
    β::Vector{Float64},
    y::AbstractVector{<:Union{Float64, Missing}}
)
    ptr = model.ptr
    brks = model.brks
    ll = 0.0
    @simd for i = 1:model.num_nodes
        if !ismissing(y[i])
            ll += (y[i] - β[i])^2
        end
    end
    return 0.5ll
end

@inline function compute_loss(
    model::GaussianGFEN,
    β::Vector{Float64},
    y::AbstractVector{<:Union{Float64, Missing}}
)
    ptr = model.ptr
    brks = model.brks
    lambdasl1 = model.lambdasl1
    lambdasl2 = model.lambdasl2

    ll = compute_negll(model, β, y)
    tv = 0.0
    @simd for j in 1:model.num_slacks
        if !(j + 1 in brks)
            δ = β[ptr[j + 1]] - β[ptr[j]]
            tv += lambdasl1[j] * abs(δ) + 0.5 * lambdasl2[j] * δ^2
        end
    end
    return ll + tv
end

@inline function residual_norms!(
    model::GaussianGFEN,
    β::Vector{Float64},
    z::Vector{Float64},
    u::Vector{Float64},
    Δz::Vector{Float64},
    z2::Vector{Float64},
    u2::Vector{Float64},
    Δz2::Vector{Float64},
    α::Float64,
    εabs_p::Float64,
    εabs_d::Float64
)
    ptr = model.ptr

    prim_norm = 0.0
    dual_norm = 0.0
    prim_size = 0.0
    dual_size = 0.0
    @simd for i = 1:model.num_slacks
        prim_norm += (β[ptr[i]] - z[i])^2 + (β[ptr[i]] - z2[i])^2
        dual_norm += Δz[i]^2 + Δz2[i]^2
        prim_size += β[ptr[i]]^2
        dual_size += u[i]^2 + u2[i]^2
    end

    # convergence in terms of relative norms
    prim_size = √prim_size
    dual_size = α * √dual_size
    prim_norm =  √prim_norm
    dual_norm = α * √dual_norm

    # convergence in terms of relative norms
    if prim_size > 0. && dual_size > 0.
        ε_p = εabs_p / prim_size + model.reltol
        ε_d = εabs_d / dual_size + model.reltol
        converged  = (prim_norm < ε_p) && (dual_norm < ε_d)
    else
        converged = (prim_size == 0) # this is a trivial extreme case
    end

    # store norms if necessary
    if model.save_norms
        push!(model.prim_norms, prim_norm)
        push!(model.dual_norms, dual_norm)
    else
        model.prim_norm = [prim_norm]
        model.dual_norm = [dual_norm]
    end

    converged, prim_norm, dual_norm, prim_size, dual_size
end

@inline function inflate_penalty!(
    model::GaussianGFEN,
    u::Vector{Float64},
    u2::Vector{Float64},
    currentα::Float64,
    edgewts::Vector{Float64},
    edgewts2::Vector{Float64},
    prim_norm::Float64,
    dual_norm::Float64,
    prim_size::Float64,
    dual_size::Float64,
)
    maxgap = model.admm_residual_balancing_gap

    if model.admm_adaptive_inflation
        ρ = sqrt(prim_norm / dual_norm)
        κ = max(ρ, 1.0 / ρ)
    else
        κ = model.admm_fixed_inflation_factor
    end
    if currentα * κ < model.admm_min_penalty
        κ = model.admm_min_penalty / currentα
    end

    α = currentα
    if (prim_norm / prim_size) > maxgap * (dual_norm / dual_size)
        α *= κ
        for j in eachindex(u)
            u[j] /= κ
            u2[j] /= κ
            edgewts[j] /= κ
            edgewts2[j] /= κ
        end
    elseif (dual_norm / dual_size) > maxgap * (prim_norm / prim_size)
        α /= κ
        for j in eachindex(u)
            u[j] *= κ
            u2[j] *= κ
            edgewts[j] *= κ
            edgewts2[j] *= κ
        end
    end
    return α
end

"""
ADMM algorithm for graph fused lasso
"""
function fit!(
        model::GaussianGFEN,
        y::AbstractVector{<:Union{Float64, Missing}};
        steps::Int = typemax(Int),
        walltime::Float64 = Inf,
        parallel::Bool = false)

    # init algorithm variables, see init_vars code for symbol descriptions
    β, z, u, Δz, Δu, z2, u2, Δz2, Δu2, α, edgewts, edgewts2, num_visits, εabs_p, εabs_d = init_vars(model, y)

    step = 0
    converged = model.converged
    maxsteps = min(model.maxsteps - model.steps, steps)

    @fastmath @inbounds begin # macros for speed
    traintime = 0.0
    while !converged && step < maxsteps && traintime < walltime
        # add to train time this iteration time
        traintime += @elapsed begin

        # update variables
        step += 1
        update_primal!(model, β, z, u, z2, u2, num_visits, y, α)
        update_slack_and_dual!(model, z, u, Δz, Δu, z2, u2, Δz2, Δu2, β, edgewts, edgewts2, parallel)

        # evaluate convergence and balance residuals
        if step % model.check_convergence_every == 0
            if model.save_loss
                loss = compute_loss(model, β, y)
                push!(model.loss, loss)
            end
            converged, prim_norm, dual_norm, prim_size, dual_size = residual_norms!(model, β, z, u, Δz, z2, u2, Δz2, α, εabs_p, εabs_d)
            converged && break

            if model.admm_residual_balancing && step % model.admm_balance_every == 0
                α = inflate_penalty!(model, u, u2, α, edgewts, edgewts2, prim_norm, dual_norm, prim_size, dual_size)
            end
        end

        end # add iteration time to train_time
    end # closes while loop
    end # closes  macros for speed

    model.trained = true
    model.beta = β
    model.slack = z
    model.dual = u
    model.admm_penalty = α
    model.steps += step
    model.converged  = converged

    converged
end

function predict(model::GaussianGFEN)
    @assert model.trained "model hasn't been trained, use fit! first"
    return model.beta
end


text:
#=
    This is a script to test the hole simulation chain with all standard detector types.
    The script also produces some output plots but only very basic ones, since it
    is only to test the core functionality of the package. (So not detector specific plots)
=#

outputdir = joinpath(ENV["HOME"], "tmp/test_solidstatedetectors.jl/")
mkpath(outputdir)
@info "Test output dir: $outputdir"

@info "Loading packages"
using Plots; pyplot()
using SolidStateDetectors

T = Float32
@info "Testing now for Float32:"

plot() # creates a plot so that the plots during the following loop pop up.

key = :InvertedCoax

for key in  [:InvertedCoax, :BEGe, :Coax, :CGD]
# for key in keys(SSD_examples)
    @info "Now test detector type: $key"

    det = SolidStateDetector{T}(SSD_examples[key])
    S = SSD.get_coordinate_system(det)

    simulation = Simulation(det);

    SSD.apply_initial_state!(simulation)
    p = if key == :CGD
        plot(simulation.electric_potential, y = 0.002)
    else
        plot(simulation.electric_potential)
    end
    savefig(joinpath(outputdir, "$(key)_0_init_setup"))

    for nrefs in [0, 1, 2]
        SSD.calculate_electric_potential!(simulation, max_refinements = nrefs)
        p = if key == :CGD
            plot(simulation.electric_potential, y = 0.002)
        else
            plot(simulation.electric_potential)
        end
        savefig(joinpath(outputdir, "$(key)_1_Electric_Potential_$(nrefs)_refinements"))
    end
    SSD.calculate_electric_potential!(simulation, max_refinements = 3)
    p = if key == :CGD
        plot(simulation.electric_potential, y = 0.002)
    else
        plot(simulation.electric_potential)
    end
    savefig(joinpath(outputdir, "$(key)_1_Electric_Potential_$(3)_refinements"))

    n_contacts = length(simulation.detector.contacts)
    for contact in simulation.detector.contacts
        SSD.calculate_weighting_potential!(simulation, contact.id, max_refinements = key == :Coax ? 1 : 2)
    end
    wp_plots = if key != :CGD
        [ plot(simulation.weighting_potentials[contact.id]) for contact in simulation.detector.contacts ]
    else
        [ plot(simulation.weighting_potentials[contact.id], y = 0.002) for contact in simulation.detector.contacts ]
    end
    plot( wp_plots..., size = (1000, 1000))
    savefig(joinpath(outputdir, "$(key)_2_Weighting_Potentials"))

    SSD.calculate_electric_field!(simulation)

    plot( simulation.electric_field.grid[1], simulation.electric_field.grid[3], SSD.get_electric_field_strength(simulation.electric_field)[:, div(length(simulation.electric_field.grid[2].ticks), 2), :]',
          st=:heatmap, title = "Electric Field Streng [V / m]", xlabel = "x / m", ylabel = "x / m", aspect_ratio = 1, size = (900, 900))
    savefig(joinpath(outputdir, "$(key)_3_Electric_Field_strength"))

    if S == :cylindrical
        plot_electric_field(simulation, φ=deg2rad(30), spacing = 8)
        savefig(joinpath(outputdir, "$(key)_3_1_Electric_Field_Lines"))
    end

    SSD.set_charge_drift_model!(simulation, ADLChargeDriftModel())

    SSD.apply_charge_drift_model!(simulation)

    pos = if key == :InvertedCoax
        CylindricalPoint{T}[ CylindricalPoint{T}( 0.02, deg2rad(10), 0.025 ) ]
    elseif key == :CGD
        CartesianPoint{T}[ CartesianPoint{T}( 0.006, 0.005, 0.005  ) ] # this point should be inside all test detectors
    elseif key == :BEGe
        CylindricalPoint{T}[ CylindricalPoint{T}( 0.016, deg2rad(10), 0.015  ) ] # this point should be inside all test detectors
    elseif key == :Coax
        CylindricalPoint{T}[ CylindricalPoint{T}( 0.016, deg2rad(10), 0.005  ) ] # this point should be inside all test detectors
    end
    energy_depos = T[1460]
    @assert in(pos[1], simulation.detector) "Test point $(pos[1]) not inside the detector $(key)."

    event = SSD.Event(pos, energy_depos)
    SSD.simulate!(event, simulation)

    plot(simulation.detector)
    plot!(event.drift_paths)
    savefig(joinpath(outputdir, "$(key)_4_charge_drift"))

    # signals[:, 2] *= -1
    plot(event.signals, size = (1200, 600), lw = 1.5)
    savefig(joinpath(outputdir, "$(key)_5_induced_signals"))

end

@info "Finished testing."
@info "Test output saved in: $outputdir"


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

using Markdown
using InteractiveUtils

# ╔═╡ f88b15de-cec6-4bc8-9b68-2a407b5aeded
begin
    using Pkg: Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(; name="Manifolds", version="0.7"),
        Pkg.PackageSpec(; name="Manopt", version="0.3"),
        Pkg.PackageSpec(; name="FiniteDiff", version="2.8"),
        Pkg.PackageSpec(; name="ReverseDiff", version="1.9"),
        # Pkg.PackageSpec(; name="PlutoUI"),
    ])
end

# ╔═╡ 856f336c-e232-4f1f-b1ac-759b4558acd1
using Manifolds, Manopt, Random, LinearAlgebra

# ╔═╡ b0769dfa-28cf-440e-9ba2-1ef488f171a9
using FiniteDiff, ReverseDiff

# ╔═╡ 0213d26a-18ac-11ec-03fd-ada5992bcea8
md"""
# Using (Euclidean) AD in Manopt.jl
"""

# ╔═╡ f3bc91ee-5871-4cba-ac89-190deb71ad0f
md"""
Since [Manifolds.jl](https://juliamanifolds.github.io/Manifolds.jl/latest/) 0.7 the support of automatic differentiation support has been extended.

This tutorial explains how to use Euclidean tools to derive a gradient for a real-valued function ``F\colon \mathcal M → ℝ``. We will consider two methods: an intrinsic variant and a variant employing the embedding. These gradients can then be used within any gradient based optimisation algorithm in [Manopt.jl](https://manoptjl.org).

While by default we use [FiniteDifferences.jl](https://juliadiff.org/FiniteDifferences.jl/latest/), you can also use [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl), [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/), [ReverseDiff.jl](https://juliadiff.org/ReverseDiff.jl/), or  [Zygote.jl](https://fluxml.ai/Zygote.jl/).
"""

# ╔═╡ d9be6c2f-65fd-4685-9005-da22bf985e28
md"""
In this Notebook we will take a look at a few possibilities to approximate or derive the gradient of a function ``f:\mathcal M \to ℝ`` on a Riemannian manifold, without computing it yourself. There is mainly two different philosophies:

1. Working _instrinsically_, i.e. stay on the manifold and in the tangent spaces. Here, we will consider approximating the gradient by forward differences.

2. Working in an embedding – there we can use all tools from functions on Euclidean spaces – finite differences or automatic differenciation – and then compute the corresponding Riemannian gradient from there.

Let's first load all packages we need.
"""

# ╔═╡ 18d7459f-eed6-489b-a096-ac77ccd781af
md"""
## 1. (Intrinsic) Forward Differences

A first idea is to generalise (multivariate) finite differences to Riemannian manifolds. Let ``X_1,\ldots,X_d ∈ T_p\mathcal M`` denote an orthonormal basis of the tangent space ``T_p\mathcal M`` at the point ``p∈\mathcal M`` on the Riemannian manifold.

We can generalise the notion of a directional derivative, i.e. for the “direction” ``Y∈T_p\mathcal M`` let ``c\colon [-ε,ε]``, ``ε>0``, be a curve with ``c(0) = p``, ``\dot c(0) = Y`` and we obtain

```math
	Df(p)[Y] = \frac{\mathrm{d}}{\mathrm{d}t} f(c(t)) = \lim_{h \to 0} \frac{1}{h}(f(\exp_p(hY))-f(p))
```

We can approximate ``Df(p)[X]`` by a finite difference scheme for an ``h>0`` as

```math
DF(p)[Y] ≈ G_h(Y) := \frac{1}{h}(f(\exp_p(hY))-f(p))
```

Furthermore the gradient ``\operatorname{grad}f`` is the Riesz representer of the differential, ie.

```math
	Df(p)[Y] = g_p(\operatorname{grad}f(p), Y),\qquad \text{ for all } Y ∈ T_p\mathcal M
```

and since it is a tangent vector, we can write it in terms of a basis as

```math
	\operatorname{grad}f(p) = \sum_{i=1}^{d} g_p(\operatorname{grad}f(p),X_i)X_i
	= \sum_{i=1}^{d} Df(p)[X_i]X_i
```

and perform the approximation from above to obtain
```math
	\operatorname{grad}f(p) ≈ \sum_{i=1}^{d} G_h(X_i)X_i
```
for some suitable step size ``h``.This comes at the cost of ``d+1`` function evaluations and ``d`` exponential maps.
"""

# ╔═╡ a3df142e-94df-48d2-be08-d1f1f3854c76
md"""
This is the first variant we can use. An advantage is, that it is _intrinsic_ in the sense that it does not require any embedding of the manifold.
"""

# ╔═╡ 9a030ac6-1f44-4fa6-8bc9-1c0278e97fe2
md""" ### An Example: The Rayleigh Quotient

The Rayleigh quotient is concerned with finding Eigenvalues (and Eigenvectors) of a symmetric matrix $A\in ℝ^{(n+1)×(n+1)}$. The optimisation problem reads

```math
F\colon ℝ^{n+1} \to ℝ,\quad F(\mathbf x) = \frac{\mathbf x^\mathrm{T}A\mathbf x}{\mathbf x^\mathrm{T}\mathbf x}
```

Minimizing this function yields the smallest eigenvalue ``\lambda_1`` as a value and the corresponding minimizer ``\mathbf x^*`` is a corresponding eigenvector.

Since the length of an eigenvector is irrelevant, there is an ambiguity in the cost function. It can be better phrased on the sphere ``𝕊^n`` of unit vectors in ``\mathbb R^{n+1}``, i.e.

```math
\operatorname*{arg\,min}_{p \in 𝕊^n} f(p) = \operatorname*{arg\,min}_{p \in 𝕊^n} p^\mathrm{T}Ap
```

We can compute the Riemannian gradient exactly as

```math
\operatorname{grad} f(p) = 2(Ap - pp^\mathrm{T}Ap)
```

so we can compare it to the approximation by finite differences.
"""

# ╔═╡ 19747159-d383-4547-9315-0ed2494904a6
begin
    Random.seed!(42)
    n = 200
    A = randn(n + 1, n + 1)
    A = Symmetric(A)
    M = Sphere(n)
    nothing
end

# ╔═╡ 41c204dd-6e4e-4a70-8f06-209a469e0680
f1(p) = p' * A'p

# ╔═╡ 2e33de5e-ffaa-422a-91d9-61f588ed1211
gradf1(p) = 2 * (A * p - p * p' * A * p)

# ╔═╡ bbd9a010-1981-45b3-bf7d-c04bcd2c2128
md"""Manifolds provides a finite difference scheme in Tangent spaces, that you can introduce to use an existing framework (if the wrapper is implemented) form Euclidean space. Here we use `FiniteDiff.jl`."""

# ╔═╡ 08456b40-74ec-4319-93e7-130b5cf70ac3
r_backend = Manifolds.TangentDiffBackend(Manifolds.FiniteDiffBackend())

# ╔═╡ 12327b62-7e79-4381-b6a7-f85b08a8251b
gradf1_FD(p) = Manifolds.gradient(M, f1, p, r_backend)

# ╔═╡ 07f9a630-e53d-45ea-b109-3d4de190723d
begin
    p = zeros(n + 1)
    p[1] = 1.0
    X1 = gradf1(p)
    X2 = gradf1_FD(p)
    norm(M, p, X1 - X2)
end

# ╔═╡ 8e5f677d-dafa-49b9-b678-3f129be31dcf
md"We obtain quite a good approximation of the gradient."

# ╔═╡ 77769eab-54dd-41dc-8125-0382e5ef0bf1
md"""
## 2. Conversion of an Euclidean Gradient in the Embedding to a Riemannian Gradient of an (not necessarily isometrically) embedded Manifold

Let ``\tilde f\colon\mathbb R^m \to \mathbb R`` be a function in the embedding of an ``n``-dimensional manifold ``\mathcal M \subset \mathbb R^m`` and ``f\colon \mathcal M \to \mathbb R`` denote the restriction of ``\tilde f`` to the manifold ``\mathcal M``.

Since we can use the push forward of the embedding to also embed the tangent space ``T_p\mathcal M``, ``p\in \mathcal M``, we can similarly obtain the differential ``Df(p)\colon T_p\mathcal M \to \mathbb R`` by restricting the differential ``D\tilde f(p)`` to the tangent space.

If both ``T_p\mathcal M`` and ``T_p\mathcal R^m`` have the same inner product, or in other words the manifold is isometrically embedded in ``R^m`` (like for example the sphere ``\mathbb S^n\subset\mathbb R^{m+1}`` then this restriction of the differential directly translates to a projection of the gradient, i.e.

```math
\operatorname{grad}f(p) = \operatorname{Proj}_{T_p\mathcal M}(\operatorname{grad} \tilde f(p))
```

More generally we might have to take a change of the metric into account, i.e.

```math
\langle  \operatorname{Proj}_{T_p\mathcal M}(\operatorname{grad} \tilde f(p)), X \rangle
= Df(p)[X] = g_p(\operatorname{grad}f(p), X)
```

or in words: we have to change the Riesz representer of the (restricted/projected) differential of ``f`` (``\tilde f``) to the one with respect to the Riemannian metric. This is done using [`change_representer`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/metric.html#Manifolds.change_representer-Tuple{AbstractManifold,%20AbstractMetric,%20Any,%20Any}).
"""

# ╔═╡ 57cda07f-e432-46af-b771-5e5a3067feac
md"""
### A continued Example
We continue with the Rayleigh Quotient from before, now just starting with the defintion of the Euclidean case in the embedding, the function ``F``.
"""

# ╔═╡ c3f3aeba-2849-4715-94e2-0c44613a2ce9
F(x) = x' * A * x / (x' * x);

# ╔═╡ 786fce04-53ef-448d-9657-31208b35fb7e
md"The cost function is the same by restriction"

# ╔═╡ c1341fef-adec-4574-a642-a1a8a9c1fee5
f2(M, p) = F(p);

# ╔═╡ 0818a62f-1bef-44f7-a33f-1ab0054e853c
md"The gradient is now computed combining our gradient scheme with ReverseDiff."

# ╔═╡ 89cd6b4b-f9ef-47ac-afd3-cf9aacf43256
function grad_f2_AD(M, p)
    return Manifolds.gradient(
        M, F, p, RiemannianProjectionBackend(Manifolds.ReverseDiffBackend())
    )
end

# ╔═╡ 7c5a8a17-6f63-4587-a94a-6936bdd3cec6
X3 = grad_f2_AD(M, p)

# ╔═╡ b3e7f57f-d87a-47c5-b8ad-48b6d205fa73
norm(M, p, X1 - X3)

# ╔═╡ 893db402-283f-4e3e-8bf7-c6f22e485efb
md"""
### An Example for a nonisometrically embedded Manifold

on the manifold ``\mathcal P(3)`` of symmetric positive definite matrices.
"""

# ╔═╡ 8494a0d6-dbf2-4eb0-a555-f00e446fbe38
md"""
The following function computes (half) the distance squared (with respect to the linear affine metric) on the manifold ``\mathcal P(3)`` to the identity, i.e. $I_3$. denoting the unit matrix we consider the function

```math
	G(q) = \frac{1}{2}d^2_{\mathcal P(3)}(q,I_3) = \lVert \operatorname{Log}(q) \rVert_F^2,
```
where $\operatorname{Log}$ denotes the matrix logarithm and ``\lVert \cdot \rVert_F`` is the Frobenius norm.
This can be computed for symmetric positive definite matrices by summing the squares of the ``\log``arithms of the eigenvalues of ``q`` and divide by two:
"""

# ╔═╡ c93eb2da-89df-4751-b086-62be604d41e6
G(q) = sum(log.(eigvals(Symmetric(q))) .^ 2) / 2

# ╔═╡ e2bf6f55-7235-4d75-8bee-a325434e32ad
md"""
We can also interpret this as a function on the space of matrices and apply the Euclidean finite differences machinery; in this way we can easily derive the Euclidean gradient. But when computing the Riemannian gradient, we have to change the representer (see again [`change_representer`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/metric.html#Manifolds.change_representer-Tuple{AbstractManifold,%20AbstractMetric,%20Any,%20Any})) after projecting onto the tangent space ``T_p\mathcal P(n)`` at ``p``.

Let's first define a point and the manifold ``N=\mathcal P(3)``.
"""

# ╔═╡ 153378ca-703d-4a84-bc63-22347399a160
rotM(α) = [1.0 0.0 0.0; 0.0 cos(α) sin(α); 0.0 -sin(α) cos(α)]

# ╔═╡ 699f0177-2c5b-434b-9eca-b6fc573e497f
q = rotM(π / 6) * [1.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 3.0] * transpose(rotM(π / 6))

# ╔═╡ 9b52d05c-4cba-4922-85a4-1a2a3c74823e
N = SymmetricPositiveDefinite(3)

# ╔═╡ cb3fe7aa-1262-48f2-9ebc-3e959c72a33e
is_point(N, q)

# ╔═╡ 13712c64-48fd-4f2a-9ee4-1949e51d316f
md"""We could first just compute the gradient using `FiniteDiff.jl`, but this yields the Euclidean gradient:"""

# ╔═╡ 64beb3dd-9507-4792-be02-ae1405704690
FiniteDiff.finite_difference_gradient(G, q)

# ╔═╡ 2be4f9e8-0331-44ac-839f-7bb71d9edef9
md"""Instead, we use the [`RiemannianProjectedBackend`](https://juliamanifolds.github.io/Manifolds.jl/latest/features/differentiation.html#Manifolds.RiemannianProjectionBackend) of `Manifolds.jl`, which in this case internally uses `FiniteDiff.jl` to compute a Euclidean gradient but then uses the conversion explained above to derive the Riemannian gradient.

We define this here again as a function `grad_G_FD` that could be used in the `Manopt.jl` framework within a gradient based optimisation.
"""

# ╔═╡ 6f1d748f-27ce-496b-8561-f16972da50cc
function grad_G_FD(N, q)
    return Manifolds.gradient(
        N, G, q, RiemannianProjectionBackend(Manifolds.FiniteDiffBackend())
    )
end

# ╔═╡ 7dd656ea-08de-4172-8a92-87ad2228ce69
G1 = grad_G_FD(N, q)

# ╔═╡ 219573d2-283f-456c-a5c3-fadd734fc157
md"""
Now, we can agaon compare this to the (known) solution of the gradient, namely the gradient of (a half) the distance suqared, i.e. ``G(q) = \frac{1}{2}d^2_{\mathcal P(3)}(q,I_3)`` is given by ``\operatorname{grad} G(q) = -\operatorname{log}_q I_3``, where ``\operatorname{log}`` is the [logarithmic map](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/symmetricpositivedefinite.html#Base.log-Tuple{SymmetricPositiveDefinite,%20Vararg{Any,%20N}%20where%20N}) on the manifold.
"""

# ╔═╡ e28a2752-877c-4ab4-a253-8d26fa9a73c2
G2 = -log(N, q, Matrix{Float64}(I, 3, 3))

# ╔═╡ 25c65878-1be6-4fec-b65e-9c1741320a41
md"""Both terms agree up to ``1.2×10^{-10}``:"""

# ╔═╡ 9a66d4f3-508d-4285-9a93-df1323575202
norm(G1 - G2)

# ╔═╡ c07fb3d0-d12f-44d7-bcab-7a0d39e6af8d
isapprox(M, q, G1, G2; atol=2 * 1e-10)

# ╔═╡ f47c70b6-ca05-498f-9e10-58c3839ca427
md""" In this case we can not use `ReverseDiff.jl`, since it can not handle the `eigvals!` function that is called internally."""

# ╔═╡ 32d8d025-3993-4d31-9eea-3463e0af1c12
md"""
## Summary

This tutorial illustrates how to use tools from Euclidean spaces, finite differences or automatic differentiation, to compute gradients on Riemannian manifolds. The scheme allows to use _any_ differentiation framework within the embedding to derive a Riemannian gradient.
"""

# ╔═╡ Cell order:
# ╟─0213d26a-18ac-11ec-03fd-ada5992bcea8
# ╟─f3bc91ee-5871-4cba-ac89-190deb71ad0f
# ╟─d9be6c2f-65fd-4685-9005-da22bf985e28
# ╟─f88b15de-cec6-4bc8-9b68-2a407b5aeded
# ╠═856f336c-e232-4f1f-b1ac-759b4558acd1
# ╠═b0769dfa-28cf-440e-9ba2-1ef488f171a9
# ╟─18d7459f-eed6-489b-a096-ac77ccd781af
# ╟─a3df142e-94df-48d2-be08-d1f1f3854c76
# ╟─9a030ac6-1f44-4fa6-8bc9-1c0278e97fe2
# ╠═19747159-d383-4547-9315-0ed2494904a6
# ╠═41c204dd-6e4e-4a70-8f06-209a469e0680
# ╠═2e33de5e-ffaa-422a-91d9-61f588ed1211
# ╟─bbd9a010-1981-45b3-bf7d-c04bcd2c2128
# ╠═08456b40-74ec-4319-93e7-130b5cf70ac3
# ╠═12327b62-7e79-4381-b6a7-f85b08a8251b
# ╠═07f9a630-e53d-45ea-b109-3d4de190723d
# ╟─8e5f677d-dafa-49b9-b678-3f129be31dcf
# ╟─77769eab-54dd-41dc-8125-0382e5ef0bf1
# ╟─57cda07f-e432-46af-b771-5e5a3067feac
# ╠═c3f3aeba-2849-4715-94e2-0c44613a2ce9
# ╟─786fce04-53ef-448d-9657-31208b35fb7e
# ╠═c1341fef-adec-4574-a642-a1a8a9c1fee5
# ╟─0818a62f-1bef-44f7-a33f-1ab0054e853c
# ╠═89cd6b4b-f9ef-47ac-afd3-cf9aacf43256
# ╠═7c5a8a17-6f63-4587-a94a-6936bdd3cec6
# ╠═b3e7f57f-d87a-47c5-b8ad-48b6d205fa73
# ╟─893db402-283f-4e3e-8bf7-c6f22e485efb
# ╟─8494a0d6-dbf2-4eb0-a555-f00e446fbe38
# ╠═c93eb2da-89df-4751-b086-62be604d41e6
# ╟─e2bf6f55-7235-4d75-8bee-a325434e32ad
# ╠═153378ca-703d-4a84-bc63-22347399a160
# ╠═699f0177-2c5b-434b-9eca-b6fc573e497f
# ╠═9b52d05c-4cba-4922-85a4-1a2a3c74823e
# ╠═cb3fe7aa-1262-48f2-9ebc-3e959c72a33e
# ╟─13712c64-48fd-4f2a-9ee4-1949e51d316f
# ╠═64beb3dd-9507-4792-be02-ae1405704690
# ╟─2be4f9e8-0331-44ac-839f-7bb71d9edef9
# ╠═6f1d748f-27ce-496b-8561-f16972da50cc
# ╠═7dd656ea-08de-4172-8a92-87ad2228ce69
# ╟─219573d2-283f-456c-a5c3-fadd734fc157
# ╠═e28a2752-877c-4ab4-a253-8d26fa9a73c2
# ╟─25c65878-1be6-4fec-b65e-9c1741320a41
# ╠═9a66d4f3-508d-4285-9a93-df1323575202
# ╠═c07fb3d0-d12f-44d7-bcab-7a0d39e6af8d
# ╟─f47c70b6-ca05-498f-9e10-58c3839ca427
# ╟─32d8d025-3993-4d31-9eea-3463e0af1c12


text:
function Primes(N::Int64)

	isprime = ones(Bool,N)

	for i = 2:int64(sqrt(N))
		if isprime[i]
			for j in (i*i):i:N
				isprime[j] = false
			end
		end
	end
	return filter(x -> isprime[x], 2:N)
end

function Euler7(N::Int64)
	primes = Primes(int(N*log(N)+N*log(log(N))))
	println(primes[N])
end

@time Euler7(10001)

max_stars_repo_name: gribeill/ProjectEuler

max_issues_repo_name: gribeill/ProjectEuler

max_forks_repo_name: gribeill/ProjectEuler


text:
using Interact, Mux
using Plots

x = y = 0:0.1:30

freqs = OrderedDict(zip(["pi/4", "π/2", "3π/4", "π"], [π/4, π/2, 3π/4, π]))

mp = @manipulate for freq1 in freqs, freq2 in slider(0.01:0.1:4π; label="freq2")
    y = @. sin(freq1*x) * sin(freq2*x)
    plot(x, y)
end

WebIO.webio_serve(page("/", req -> mp), 8000) # serve on a random port

max_stars_repo_name: raphbacher/julia-intro

max_issues_repo_name: raphbacher/julia-intro

max_forks_repo_name: raphbacher/julia-intro


text:
using Seqeval
using Test

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


text:
# Gaussian One-D engines #
using Dates

struct NoneProbabilities <: GaussianProbabilities end
struct NaiveProbabilities <: GaussianProbabilities end
struct DigitalProbabilities <: GaussianProbabilities end

struct Gaussian1DSwaptionEngine{G <: Gaussian1DModel, Y <: YieldTermStructure, P <: GaussianProbabilities} <: PricingEngine
  model::G
  integrationPoints::Int
  stddevs::Float64
  extrapolatePayoff::Bool
  flatPayoffExtrapolation::Bool
  discountCurve::Y
  probabilities::P

  Gaussian1DSwaptionEngine{G, Y, P}(model::G,
                          integrationPoints::Int,
                          stddevs::Float64,
                          extrapolatePayoff::Bool,
                          flatPayoffExtrapolation::Bool,
                          discountCurve::Y,
                          probabilities::P) where {G, Y, P} =
    new{G, Y, P}(model, integrationPoints, stddevs, extrapolatePayoff, flatPayoffExtrapolation, discountCurve, probabilities)
end

Gaussian1DSwaptionEngine(model::G,
                        integrationPoints::Int = 64,
                        stddevs::Float64 = 7.0,
                        extrapolatePayoff::Bool = true,
                        flatPayoffExtrapolation::Bool = false,
                        discountCurve::Y = NullYieldTermStructure(),
                        probabilities::P = NoneProbabilities()) where {G <: Gaussian1DModel, Y <: YieldTermStructure, P <: GaussianProbabilities} =
                          Gaussian1DSwaptionEngine{G, Y, P}(model, integrationPoints, stddevs, extrapolatePayoff, flatPayoffExtrapolation, discountCurve, probabilities)
# methods #
function _calculate!(pe::Gaussian1DSwaptionEngine, swaption::Swaption)
  isa(swaption.delivery, SettlementPhysical) || error("cash-settled swaptions not yet implemented")
  settlement = reference_date(pe.model.ts)

  if swaption.exercise.dates[end] <= settlement
    swaption.results.value = 0.0 # swaption is expired
    return swaption
  end
  idx = length(swaption.exercise.dates)
  # minIdxAlive = searchsortedlast(swaption.exercise.dates, settlement) + 1
  # minIdxAlive = minIdxAlive == 1 ? 2 : minIdxAlive
  minIdxAlive = upper_bound(swaption.exercise.dates, settlement)

  swap = swaption.swap
  optType = isa(swap.swapT, Payer) ? Call() : Put()
  fixedSchedule = swap.fixedSchedule
  floatingSchedule = swap.floatSchedule

  npv0 = zeros(2 * pe.integrationPoints + 1)
  npv1 = zeros(2 * pe.integrationPoints + 1)
  z = y_grid(pe.model, pe.stddevs, pe.integrationPoints)
  p = zeros(length(z))

  # for probability computation
  npvp0 = Vector{Vector{Float64}}()
  npvp1 = Vector{Vector{Float64}}()
  if !isa(pe.probabilities, NoneProbabilities)
    @inbounds @simd for i = 1:(idx - (minIdxAlive - 1) + 1)
      npvTmp0 = zeros(2 * pe.integrationPoints + 1)
      npvTmp1 = zeros(2 * integrationPoints + 1)
      push!(npvp0, npvTmp0)
      push!(npvp1, npvTmp1)
    end
  end
  # end probability computation

  expiry0 = expiry1 = Date(0)
  expiry0Time = expiry1Time = -1.0

  while true
    if idx == minIdxAlive - 1
      expiry0 = settlement
    else
      expiry0 = swaption.exercise.dates[idx]
    end

    expiry0Time = max(time_from_reference(pe.model.ts, expiry0), 0.0)
    j1 = upper_bound(fixedSchedule.dates, expiry0 - Dates.Day(1))
    k1 = upper_bound(floatingSchedule.dates, expiry0 - Dates.Day(1))

    # for multi-threading / multi-processing, we need to trigger computation here
    # if expiry0 > settlement
    #   for l = k1:length(get_floating_coupons(swaption.swap))
    #     forward_rate(pe.model, swaption.swap.legs[2].coupons[l].fixingDate, expiry0, 0.0, swaption.swap.iborIndex)
    #     zerobond(pe.model, get_floating_pay_dates(swaption.swap)[l], expiry0, 0.0, pe.discountCurve)
    #   end
    #
    #   for l = j1:length(get_fixed_coupons(swaption.swap))
    #     zerobond(pe.model, get_fixed_pay_dates(swaption.swap)[l], expiry0, 0.0, pe.discountCurve)
    #   end
    #
    #   numeraire(pe.model, expiry0Time, 0.0, pe.discountCurve)
    # end

    @inbounds @simd for k = 1:(expiry0 > settlement ? length(npv0) : 1)
      price = 0.0

      if expiry1Time != -1.0
        yg = y_grid(pe.model, pe.stddevs, pe.integrationPoints, expiry1Time, expiry0Time, expiry0 > settlement ? z[k] : 0.0)
        payoff0 = CubicInterpolation(Spline(), Lagrange(), Lagrange(), z, npv1)
        for i in eachindex(yg)
          p[i] = payoff0(yg[i])
        end

        payoff1 = CubicInterpolation(Spline(), Lagrange(), Lagrange(), z, p)
        for i = 1:length(z) - 1
          price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[i], payoff1.b[i], payoff1.a[i], p[i], z[i], z[i], z[i + 1])
        end

        if pe.extrapolatePayoff
          if pe.flatPayoffExtrapolation
            price += gaussian_shifted_polynomial_integral(pe.model, 0.0, 0.0, 0.0, 0.0, p[length(z) - 1], z[end - 1], z[end], 100.0)
            price += gaussian_shifted_polynomial_integral(pe.model, 0.0, 0.0, 0.0, 0.0, p[1], z[1], -100.0, z[1])
          else
            if isa(optType, Call)
              price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[length(z) - 1], payoff1.b[length(z) - 1], payoff1.a[length(z) - 1], p[length(z) - 1], z[end - 1], z[end], 100.0)
            else
              price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[1], payoff1.b[1], payoff1.a[1], p[1], z[1], -100.0, z[1])
            end
          end
        end
      end
      # println("k: $k, price: $price")
      npv0[k] = price
      # for probability computation
      if !isa(pe.probabilities, NoneProbabilities)
        for m = 1:length(npvp0)
          price = 0.0
          if expiry1Time != -1.0
            yg = y_grid(pe.model, pe.stddevs, pe.integrationPoints, expiry1Time, expiry0Time, expiry0 > settlement ? z[k] : 0.0)
            payoff0 = CubicInterpolation(Spline(), Lagrange(), Lagrange(), z, npvp1[m])
            for i in eachindex(yg)
              p[i] = payoff0(yg[i])
            end

            payoff1 = CubicInterpolation(Spline(), Lagrange(), Lagrange(), z, p)
            for i = 1:length(z) - 1
              price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[i], payoff1.b[i], payoff1.a[i], p[i], z[i], z[i], z[i + 1])
            end

            if pe.extrapolatePayoff
              if pe.flatPayoffExtrapolation
                price += gaussian_shifted_polynomial_integral(pe.model, 0.0, 0.0, 0.0, 0.0, p[length(z) - 1], z[end - 1], z[end], 100.0)
                price += gaussian_shifted_polynomial_integral(pe.model, 0.0, 0.0, 0.0, 0.0, p[1], z[1], -100.0, z[1])
              else
                if isa(optType, Call)
                  price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[length(z) - 1], payoff1.b[length(z) - 1], payoff1.a[length(z) - 1], p[length(z) - 1], z[end - 1], z[end], 100.0)
                else
                  price += gaussian_shifted_polynomial_integral(pe.model, 0.0, payoff1.c[1], payoff1.b[1], payoff1.a[1], p[1], z[1], -100.0, z[1])
                end
              end
            end
          end
          npvp0[m][k] = price
        end
      end
      # end probability computation
      if expiry0 > settlement
        floatingLegNPV = 0.0
        for l = k1:length(get_floating_coupons(swaption.swap))
          floatingLegNPV += swaption.swap.nominal * get_floating_accrual_times(swaption.swap)[l] *
            (get_floating_spreads(swaption.swap)[l] + forward_rate(pe.model, swaption.swap.legs[2].coupons[l].fixingDate, expiry0, z[k], swaption.swap.iborIndex)) *
            zerobond(pe.model, get_floating_pay_dates(swaption.swap)[l], expiry0, z[k], pe.discountCurve)
        end
        # println("floatingLeg: ", floatingLegNPV)
        fixedLegNPV = 0.0
        for l = j1:length(get_fixed_coupons(swaption.swap))
          fixedLegNPV += get_fixed_coupons(swaption.swap)[l] * zerobond(pe.model, get_fixed_pay_dates(swaption.swap)[l], expiry0, z[k], pe.discountCurve)
        end
        # println("fixedLeg: ", fixedLegNPV)
        exerciseValue = (isa(optType, Call) ? 1.0 : -1.0) * (floatingLegNPV - fixedLegNPV) / numeraire(pe.model, expiry0Time, z[k], pe.discountCurve)
        # probability computation
        if !isa(pe.probabilities, NoneProbabilities)
          if idx == length(swaption.exercise.dates) # if true we are at the latest date so we init the no call prob
            npvp0[end][k] = isa(pe.probabilities, NaiveProbabilities) ? 1.0 : 1.0 / (zerobond(pe.model, expiry0Time, 0.0, 0.0, pe.discountCurve) * numeraire(pe.model, expiry0, z[k], pe.discountCurve))
          end

          if exerciseValue >= npv0[k]
            npvp0[idx - (minIdxAlive - 1)][k] = isa(pe.probabilities, NaiveProbabilities) ? 1.0 : (zerobond(pe.model, expiry0Time, 0.0, 0.0, pe.discountCurve) * numeraire(pe.model, expiry0Time, z[k], pe.discountCurve))
            for ii = idx - (minIdxAlive - 1) + 1:length(npvp0)
              npvp0[ii][k] = 0.0
            end
          end
        end
        # end probability computation
        npv0[k] = max(npv0[k], exerciseValue)
      end
    end
    npv1, npv0 = npv0, npv1

    # for probability computation
    if !isa(pe.probabilities, NoneProbabilities)
      for i in eachindex(npvp0)
        npvp1[i], npvp0[i] = npvp0[i], npvp1[i]
      end
    end
    # end probability computation

    expiry1 = expiry0
    expiry1Time = expiry0Time
    idx -= 1
    if idx < minIdxAlive - 1
      break
    end
  end
  # println("NPV1: ", npv1)

  swaption.results.value = npv1[1] * numeraire(pe.model, 0.0, 0.0, pe.discountCurve)

  # for probability computation
  if !isa(pe.probabilities, NoneProbabilities)
    prob = zeros(length(npvp0))
    for i in eachindex(npvp0)
      prob[i] = npvp[i][1] * (isa(pe.probabilities, NaiveProbabilities) ? 1.0 : numeraire(pe.model, 0.0, 0.0, pe.discountCurve))
    end
    swaption.results.additionalResults["probabilities"] = prob
  end

  return swaption
end


text:
using Bukdu # render HTML Router
using .Bukdu.System: MissingController
import .Bukdu.System: not_found

function not_found(c::MissingController)
    c.conn.request.response.status = 404 # 404 Not Found
    render(HTML, "custom 404 not found message")
end

if PROGRAM_FILE == basename(@__FILE__)

Bukdu.start(8080)

Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if


text:
"""
    ⊗(A, B)

Compute the kronecker product, ``A\\otimes B``.

# Examples
```julia-repl
julia> [1 0; 0 1]⊗[1 0; 0 -1]
4×4 Array{Int64,2}:
 1   0  0   0
 0  -1  0   0
 0   0  1   0
 0   0  0  -1
```
"""
⊗(A, B) = kron(A, B)

"""
    ⊗(A, B, C...)

Compute the kronecker product, ``A\\otimes B\\otimes C \\cdots``.
"""
function ⊗(A, B, C...)
    A = kron(A, B)
    for Ci in C
        A = kron(A, Ci)
    end
    return A
end


text:
# ------------------------------------
# Secure the array for explicit and implicit
# ------------------------------------
function common_allocation(cellxmax, cellymax, nval)
    # define at cell center
    Qbase  = zeros(cellxmax, cellymax, nval)           # primitive variables
    #Qbase_ave = zeros(cellxmax, cellymax, nval)        # average primitive variables
    volume = zeros(cellxmax, cellymax)                 # volume

    cellcenter = zeros(cellxmax, cellymax, 2)          # cell center

    wally = zeros(cellxmax, cellymax)                  # dis from wall
    yplus = ones(cellxmax, cellymax) * 1e5             # yplus

    Qcon     = zeros(cellxmax, cellymax, nval)         # conserved variables
    Qcon_hat = zeros(cellxmax, cellymax, nval)         # conserved variables in the general coordinate system

    mu     = zeros(cellxmax, cellymax)                 # viscosity
    lambda = zeros(cellxmax, cellymax)                 # thermal Conductivity

    RHS = zeros(cellxmax, cellymax, nval)              # right hand side

    # define at cell boundaries
    dx = zeros(cellxmax+1, cellymax)                   # distance from cell center
    dy = zeros(cellxmax, cellymax+1)                   # distance from cell center

    E_adv_hat = zeros(cellxmax+1,   cellymax, nval)    # flux of advection in the x-direction
    F_adv_hat = zeros(  cellxmax, cellymax+1, nval)    # flux of advection in the y-direction

    E_vis_hat = zeros(cellxmax+1,   cellymax, nval)    # flux of viscosity in the x-direction
    F_vis_hat = zeros(  cellxmax, cellymax+1, nval)    # flux of viscosity in the y-direction

    return Qbase, volume, cellcenter, wally, yplus, dx, dy, Qcon, Qcon_hat, mu, lambda,
            E_adv_hat, F_adv_hat, E_vis_hat, F_vis_hat, RHS
end

function allocation_implicit(cellxmax, cellymax, nval)
    # define at cell center
    Qbasen     = zeros(cellxmax, cellymax, nval)        # primitive variables for inner iteration
    Qconn      = zeros(cellxmax, cellymax, nval)        # conserved variables for inner iteration
    Qconn_hat  = zeros(cellxmax, cellymax, nval)        # conserved variables in general coordinate system for inner iteration
    Qbasem     = zeros(cellxmax, cellymax, nval)        # primitive variables for inner iteration

    dtau         = zeros(cellxmax, cellymax)            # computational time steps

    A_adv_hat_p = zeros(cellxmax, cellymax, nval, nval) # Jacobian matrix A+ for one-wave approximation
    A_adv_hat_m = zeros(cellxmax, cellymax, nval, nval) # Jacobian matrix A- for one-wave approximation
    B_adv_hat_p = zeros(cellxmax, cellymax, nval, nval) # Jacobian matrix B+ for one-wave approximation
    B_adv_hat_m = zeros(cellxmax, cellymax, nval, nval) # Jacobian matrix B- for one-wave approximation
    A_beta_shig = zeros(cellxmax, cellymax)             # A sigma for one-wave approximation
    B_beta_shig = zeros(cellxmax, cellymax)             # B sigma for one-wave approximation

    jalphaP = zeros(cellxmax, cellymax)                 # Jacobian matrix for viscosity
    jbetaP  = zeros(cellxmax, cellymax)                 # Jacobian matrix for viscosity

    delta_Q      = zeros(cellxmax, cellymax, nval)      # delta Q for lusgs
    delta_Q_temp = zeros(cellxmax, cellymax, nval)      # temporary delta Q for lusgs

    D  = zeros(cellxmax, cellymax)                      # Diagonal of LHS
    Lx = zeros(cellxmax, cellymax, nval, nval)          # lower of LHS
    Ly = zeros(cellxmax, cellymax, nval, nval)          # lower of LHS
    Ux = zeros(cellxmax, cellymax, nval, nval)          # upper of LHS
    Uy = zeros(cellxmax, cellymax, nval, nval)          # upper of LHS

    LdQ = zeros(cellxmax, cellymax, nval)               # lower of LHS
    UdQ = zeros(cellxmax, cellymax, nval)               # upper of LHS

    RHS_temp = zeros(cellxmax, cellymax, nval)          # temporary RHS

    # define at cell boundaries
    lambda_facex = zeros(cellxmax+1, cellymax)          # lambda for computational time steps
    lambda_facey = zeros(cellxmax, cellymax+1)          # lambda for computational time steps

    # misc
    norm2 = zeros(nval)                                 # Residuals by norm-2
    I = zeros(nval, nval)                               # identity matrix
    for l in 1:nval
        I[l,l] = 1.0
    end

    return Qbasen, Qconn, Qconn_hat, Qbasem, dtau, lambda_facex, lambda_facey,
            A_adv_hat_m, A_adv_hat_p, B_adv_hat_m, B_adv_hat_p, A_beta_shig, B_beta_shig,
            jalphaP, jbetaP, delta_Q, delta_Q_temp, D, Lx, Ly, Ux, Uy, LdQ, UdQ, RHS_temp,
            norm2, I
end


text:
module Ex5_11
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.11 - Solution of an AFDIP
println("Example 5.11")

# define s as an improper transfer function
s = rtf('s');
# define Gu(s), Gf(s) and Gw(s)
Gu = [(s+1)/(s+2); (s+2)/(s+3)];     # enter Gu(s)
Gw = [1/(s+2); 0];                   # enter Gw(s)
Gf = [(s+1)/(s+2) 0; 0 1];           # enter Gf(s)
p = 2; mu = 1; mw = 1; mf = 2;       # enter dimensions
S = eye(mf) .> 0

# Procedure AFDI

# Step 1): choose the left nullspace as Q1 = [I -Gu] and
# form Rf1 = Q1*[Gf;0] and Rw1 = Q1*[Gw;0]
Q1 = dss([eye(p) -Gu]); Rf1 = dss(Gf); Rw1 = dss(Gw);
atol = 1.e-7;

# Step 2): determine Q[i] and corresponding Rf[i] and Rw[i]

# initialize overall filter Q and fault detection system Rf
nb = size(S,1);      # number of necessary filters
Qt = similar(Vector{typeof(Q1)},nb)
Rfwt = similar(Vector{typeof(Q1)},nb)
scale = [1, -1]  # set appropriate scaling to match example

for i = 1:nb
    # perform Procedure AFD or EFD to compute Q[i]
    indd = Vector(1:mf)[S[i,:] .== false]
    Qi1 = glnull(Rf1[:,indd])[1]
    # initialize Q[i], Rf[i] and Rw[i]
    Qi = Qi1*Q1; Rfi = Qi1*Rf1; Rwi = Qi1*Rw1;

    if norm(evalfr(Rwi,rand())) > 0.0001
       # compute the quasi-co-outer-co-inner factorization
       Rwi, Rwo, = goifac(Rwi;atol)
       # update
       Qi = Rwo\Qi; Rfi = Rwo\Rfi
    end
    # update the solution if [Q[i], Rf[i], Rw[i]] is improper or unstable
    Qi_Rfi_Rwi, M = glcf([Qi Rfi Rwi], sdeg = -3, smarg = -3, mindeg = true)

    # adjust denominator M to unit infinity norm to match example
    Mnorm = ghinfnorm(M)[1]
    Qt[i] = (scale[i]/Mnorm)*Qi_Rfi_Rwi[:,1:p+mu]
    Rfwt[i] = (scale[i]/Mnorm)*Qi_Rfi_Rwi[:,p+mu+1:end]

    println("Q[$i] = $(dss2rm(Qt[i], atol = 1.e-7))")
    println("Rf[$i] = $(dss2rm(Rfwt[i][:,1:mf], atol = 1.e-7))")
    println("Rw[$i] = $(dss2rm(Rfwt[i][:,mf+1:end], atol = 1.e-7))")
end
Q = FDIFilter(gminreal.(Qt; atol), p, mu)
R = FDIFilterIF(gminreal.(Rfwt; atol), 0, 0, mf, mw)
println("Q = ")
display(Q)
println("R = ")
display(R)
println("gap = $(fdif2ngap(R,S)[1])")

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
export FOTemp, gettemperature, gettemperatures, logtemperature, logtemperatures

include("FOTempLowLevel.jl")


text:
module IpyJL
using JSON
"""
## Cell(::Any, ::String, ::Any, ::Dict, ::Integer)
The cell type is just a Julian equivalent to the JSON data that is read in
for Jupyter cells.
### Properties
- outputs::Any - Output of the cell
- ctype::String - Cell type (code/md)
- cont::Any - The content of the cell
- meta::Dict - Meta information for a given cell.
- n::Integer - The execution position of the cell.
### Constructors
- Cell(::Dict) Constructs cells from a dictionary of cell-data.
"""
mutable struct Cell
        outputs::Any
        ctype::String
        cont::Any
        meta::Dict
        n::Integer
        function Cell(nb_cdict::Dict)
            n = 0
            outputs = ""
            try
                n = nb_cdict["execution_count"]
            catch
                n = 0
            end
            if isnothing(n)
                n = 0
            end
            try
                outputs = nb_cdict["outputs"]
            catch
                outputs = ""
            end
            new(outputs, nb_cdict["cell_type"], nb_cdict["source"],
             nb_cdict["metadata"], n)
    end
end
"""
## read_ipynb(f::String) -> Vector{Cell}
Reads an IPython notebook into a vector of cells.
### example
read_ipynb("helloworld.ipynb")
"""
function read_ipynb(f::String)
    file = open(f)
    j = JSON.parse(file)
    [Cell(cell) for cell in j["cells"]]
end
"""
## ipynbjl(ipynb_path::String, output_path::String)
Reads notebook at **ipynb_path** and then outputs as .jl Julia file to
**output_path**.
### example
ipynbjl("helloworld.ipynb", "helloworld.jl")
"""
function ipynbjl(ipynb_path::String, output_path::String)
    cells = read_ipynb(ipynb_path)
    output = cells_to_string(cells)
    open(output_path, "w") do file
           write(file, output)
    end;
end
"""
## cells_to_string(::Vector{Any}) -> ::String
Converts an array of Cell types into text.
"""
function cells_to_string(cells::Any)
    f = ""
    for cell in cells
        line = ""
        header = string("# ", cell.n, "\n")
        if cell.ctype != "code"
            line = line * string("cellmd", string(cell.n), " = \"\"\"")
            line = line * string(sep(cell.cont)) * "\"\"\"\n"
        else
            line = ""
            line = line * string(sep(cell.cont)) * "\n"
        end
        f = f * header * line
    end
    return(f)
end

"""
## sep(::Any) -> ::String
Separates and parses lines of individual cell content via an array of strings.
Returns string of concetenated text.
"""
function sep(content::Any)
    total = string()
    if length(content) == 0
        return("")
    end
    for line in content
        total = total * string(line)
    end
    total
end
export ipynbjl
end # module


text:
# # [COEXIST Multi-Generational COVID Model](@id coexist_example)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/covid/coexist/coexist.ipynb)

using AlgebraicPetri

using LabelledArrays
using OrdinaryDiffEq
using Plots
using JSON

using Catlab
using Catlab.CategoricalAlgebra
using Catlab.Graphics
using Catlab.Programs
using Catlab.Theories
using Catlab.WiringDiagrams
using Catlab.Graphics.Graphviz: run_graphviz

display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>"1"));
save_fig(g, fname::AbstractString, format::AbstractString) = begin
    open(string(fname, ".", format), "w") do io
        run_graphviz(io, g, format=format)
    end
end

# Define helper functions for defining the two types of
# reactions in an epidemiology Model. Either a state
# spontaneously changes, or one state causes another to change

ob(x::Symbol,n::Int) = codom(Open([x], LabelledReactionNet{Number,Int}(x=>n), [x])).ob;
function spontaneous_petri(transition::Symbol, rate::Number,
                           s::Symbol, s₀::Int,
                           t::Symbol, t₀::Int)
    Open(LabelledReactionNet{Number,Int}(unique((s=>s₀,t=>t₀)), (transition,rate)=>(s=>t)))
end;
function exposure_petri(transition::Symbol, rate::Number,
                        s::Symbol, s₀::Int,
                        e::Symbol, e₀::Int,
                        t::Symbol, t₀::Int)
    Open(LabelledReactionNet{Number,Int}(unique((s=>s₀,e=>e₀,t=>t₀)), (transition,rate)=>((s,e)=>(t,e))))
end;

# Set arrays of initial conditions and rates to use in functor

pop = [8044056, 7642473, 8558707, 9295024,8604251,9173465,7286777,5830635,3450616] .- (4*1000);
N = sum(pop) + length(pop)*4*1000;
social_mixing_rate =
  [[5.10316562022642,1.28725377551533,1.30332531065247,2.31497083312315,1.1221598200343,0.606327539457772,0.453266757158743,0.177712174722219,0.0111726265254263],
   [1.15949254996891,8.00118824220649,1.24977685411394,1.51298690806342,1.88877951844257,0.835804485358679,0.431371281973645,0.343104864504218,0.0324429672946592],
   [1.19314902456243,1.2701954426234,3.55182053724384,1.81286158254244,1.80561825747571,1.29108026766182,0.708613434860661,0.248559044477893,0.0215323291988856],
   [1.83125260045684,1.32872195974583,1.56648238384012,2.75491288061819,1.94613663227464,1.2348814962672,0.863177586322153,0.244623623638873,0.0394364256673532],
   [0.910395333788561,1.7011898591446,1.60014517035071,1.99593275526656,2.90894801031624,1.37683234043657,0.859519958701156,0.488960115017174,0.110509077357166],
   [0.56560186656657,0.865574490657954,1.31557291022074,1.45621698394508,1.58310342861768,1.92835669973181,0.963568493650797,0.463041280007004,0.183483677017087],
   [0.544954016221808,0.575775829452094,0.930622416907882,1.31190809759635,1.27375718214796,1.24189546255302,1.32825334016313,0.66235513907445,0.0946971569608397],
   [0.319717318035767,0.68528632728864,0.488468642570909,0.556345582530282,1.08429412751444,0.893028152305907,0.991137484161889,1.17651345255182,0.12964732712923],
   [0.201086389216809,0.648252461859761,0.423327560644352,0.897268061280577,2.4516024037254,3.54014694719397,1.41761515077768,1.29700599099082,1.0189817510854]];

fatality_rate = [0.00856164, 0.03768844, 0.02321319, 0.04282494, 0.07512237, 0.12550367, 0.167096  , 0.37953452, 0.45757006];

# Define an `oapply` function that connects the building block Petri nets to
# the operations we will use in the model.

F(ex, n) = oapply(ex, Dict(
    :exposure=>exposure_petri(Symbol(:exp_, n), 1*social_mixing_rate[n][n]/pop[n], Symbol(:S,n), pop[n], Symbol(:I,n), 1000, Symbol(:E,n), 1000),
    :exposure_e=>exposure_petri(Symbol(:exp_e, n), .01*social_mixing_rate[n][n]/pop[n], Symbol(:S,n), pop[n], Symbol(:E,n),1000, Symbol(:E,n),1000),
    :exposure_i2=>exposure_petri(Symbol(:exp_i2, n), 6*social_mixing_rate[n][n]/pop[n], Symbol(:S,n), pop[n], Symbol(:I2,n), 1000, Symbol(:E,n),1000),
    :exposure_a=>exposure_petri(Symbol(:exp_a, n), 5*social_mixing_rate[n][n]/pop[n], Symbol(:S,n), pop[n], Symbol(:A,n),1000, Symbol(:E,n),1000),
    :progression=>spontaneous_petri(Symbol(:prog_, n), .25, Symbol(:I,n), 1000, Symbol(:I2,n), 1000),
    :asymptomatic_infection=>spontaneous_petri(Symbol(:asymp_, n), .86/.14*.2, Symbol(:E,n), 1000, Symbol(:A,n), 1000),
    :illness=>spontaneous_petri(Symbol(:ill_, n), .2, Symbol(:E,n), 1000, Symbol(:I,n), 1000),
    :asymptomatic_recovery=>spontaneous_petri(Symbol(:arec_, n), 1/15, Symbol(:A,n), 1000, Symbol(:R,n), 0),
    :recovery=>spontaneous_petri(Symbol(:rec_, n), 1/6, Symbol(:I2,n), 1000, Symbol(:R,n), 0),
    :recover_late=>spontaneous_petri(Symbol(:rec2_, n), 1/15, Symbol(:R,n), 0, Symbol(:R2,n), 0),
    :death=>spontaneous_petri(Symbol(:death2_, n), (1/15)*(fatality_rate[n]/(1-fatality_rate[n])), Symbol(:I2,n), 1000, Symbol(:D,n), 0)));

# Define the COEXIST model using the `@relation` macro

coexist = @relation (s, e, i, i2, a, r, r2, d) begin
    exposure(s, i, e)
    exposure_i2(s, i2, e)
    exposure_a(s, a, e)
    exposure_e(s, e, e)
    asymptomatic_infection(e, a)
    asymptomatic_recovery(a, r)
    illness(e, i)
    progression(i, i2)
    death(i2, d)
    recovery(i2, r)
    recover_late(r, r2)
end;
display_uwd(coexist)

# Define an `oapply` function that can be used to create a model of
# cross exposure between two sets of populations

F_cx(ex, x,y) = oapply(ex, Dict(
    :exposure=>exposure_petri(Symbol(:exp_, x,y), 1*social_mixing_rate[x][y]/(pop[x]+pop[y]), Symbol(:S,x), pop[x], Symbol(:I,y), 1000, Symbol(:E,x), 1000),
    :exposure_e=>exposure_petri(Symbol(:exp_e, x,y), .01*social_mixing_rate[x][y]/(pop[x]+pop[y]), Symbol(:S,x), pop[x], Symbol(:E,y),1000, Symbol(:E,x),1000),
    :exposure_a=>exposure_petri(Symbol(:exp_a, x,y), 5*social_mixing_rate[x][y]/(pop[x]+pop[y]), Symbol(:S,x), pop[x], Symbol(:A,y),1000, Symbol(:E,x),1000),
    :exposure_i2=>exposure_petri(Symbol(:exp_i2, x,y), 6*social_mixing_rate[x][y]/(pop[x]+pop[y]), Symbol(:S,x), pop[x], Symbol(:I2,y), 1000, Symbol(:E,x),1000),
    :exposure′=>exposure_petri(Symbol(:exp_, y,x), 1*social_mixing_rate[y][x]/(pop[x]+pop[y]), Symbol(:S,y), pop[y], Symbol(:I,x), 1000, Symbol(:E,y), 1000),
    :exposure_e′=>exposure_petri(Symbol(:exp_e, y,x), .01*social_mixing_rate[y][x]/(pop[x]+pop[y]), Symbol(:S,y), pop[y], Symbol(:E,x),1000, Symbol(:E,y),1000),
    :exposure_a′=>exposure_petri(Symbol(:exp_a, y,x), 5*social_mixing_rate[y][x]/(pop[x]+pop[y]), Symbol(:S,y), pop[y], Symbol(:A,x),1000, Symbol(:E,y),1000),
    :exposure_i2′=>exposure_petri(Symbol(:exp_i2, y,x), 6*social_mixing_rate[y][x]/(pop[x]+pop[y]), Symbol(:S,y), pop[y], Symbol(:I2,x), 1000, Symbol(:E,y),1000)
  ),
  Dict(
    :s=>ob(Symbol(:S, x), pop[x]),
    :e=>ob(Symbol(:E, x), 1000),
    :a=>ob(Symbol(:A, x), 1000),
    :i=>ob(Symbol(:I, x), 1000),
    :i2=>ob(Symbol(:I2, x), 1000),
    :r=>ob(Symbol(:R, x), 0),
    :r2=>ob(Symbol(:R2, x), 0),
    :d=>ob(Symbol(:D, x), 0),
    :s′=>ob(Symbol(:S, y), pop[y]),
    :e′=>ob(Symbol(:E, y), 1000),
    :a′=>ob(Symbol(:A, y), 1000),
    :i′=>ob(Symbol(:I, y), 1000),
    :i2′=>ob(Symbol(:I2, y), 1000),
    :r′=>ob(Symbol(:R, y), 0),
    :r2′=>ob(Symbol(:R2, y), 0),
    :d′=>ob(Symbol(:D, y), 0)
  ));

# Use this new presentation to define a model
# of cross exposure between two populations

crossexposure = @relation (s, e, i, i2, a, r, r2, d, s′, e′, i′, i2′, a′, r′, r2′, d′) begin
    exposure(s, i′, e)
    exposure_i2(s, i2′, e)
    exposure_a(s, a′, e)
    exposure_e(s, e′, e)
    exposure′(s′, i, e′)
    exposure_i2′(s′, i2, e′)
    exposure_a′(s′, a, e′)
    exposure_e′(s′, e, e′)
end;
display_uwd(crossexposure)

# To combine these two models, we need to create a final relational model and
# use the `bundle_legs` function in our `oapply` that enables us to model 3
# population wires instead of each individual state as a wire. Each of these
# populations has their own COEXIST model, and interact through cross exposure

bundled_cross(x,y) = bundle_legs(F_cx(crossexposure, x, y), [tuple([1:8;]...), tuple([9:16;]...)])
bundled_coex(x) = bundle_legs(F(coexist, x), [tuple([1:8;]...)])
F_tcx(ex) = oapply(ex, Dict(
    :crossexp12=>bundled_cross(3,4),
    :crossexp13=>bundled_cross(3,5),
    :crossexp23=>bundled_cross(4,5),
    :coex1=>bundled_coex(3),
    :coex2=>bundled_coex(4),
    :coex3=>bundled_coex(5)));

threeNCoexist = @relation (pop1, pop2, pop3) begin
    crossexp12(pop1, pop2)
    crossexp13(pop1, pop3)
    crossexp23(pop2, pop3)
    coex1(pop1)
    coex2(pop2)
    coex3(pop3)
end;
display_uwd(threeNCoexist)
#-
threeNCoexist_algpetri = apex(F_tcx(threeNCoexist))
Graph(threeNCoexist_algpetri);
save_fig(Graph(threeNCoexist_algpetri), "3ncoexist_petri", "svg"); # hide

# ![3-generation COEXIST model petri net](3ncoexist_petri.svg)

# We can JSON to convert this Petri net into an
# easily shareable format

JSON.print(tables(threeNCoexist_algpetri), 2)

# We can now easily generate a solver for DifferentialEquations.jl
# because we encoded the intitial parameters and rates throughout
# the construction of the model, the final result knows its
# concentrations and rates.

tspan = (0.0,100.0);
prob = ODEProblem(vectorfield(threeNCoexist_algpetri),concentrations(threeNCoexist_algpetri),tspan,rates(threeNCoexist_algpetri));
sol = solve(prob,Tsit5());
plot(sol, xlabel="Time", ylabel="Number of people")

# If we want to model other intervention methods,
# we can simply adjust the rates of exposure to
# represent stay at home orders and mask wearing.
# Because of how we have defined our rates, we can
# simply update the social mixing rates, and
# resolve the model.

for i in 1:length(social_mixing_rate)
  for j in 1:length(social_mixing_rate[1])
    social_mixing_rate[i][j] = social_mixing_rate[i][j] / (i != j ? 10 : 5);
  end
end
threeNCoexist_algpetri = apex(F_tcx(threeNCoexist));

prob = ODEProblem(vectorfield(threeNCoexist_algpetri),concentrations(threeNCoexist_algpetri),tspan,rates(threeNCoexist_algpetri));
sol = solve(prob,Tsit5());
plot(sol, xlabel="Time", ylabel="Number of people")


text:
############################################################################
### based on -2ll
############################################################################

AIC(sem_fit) = minus2ll(sem_fit) + 2n_par(sem_fit)

############################################################################
### based on χ² - 2df
############################################################################

# AICχ²(sem_fit::SemFit{Mi, S, Mo, O} where {Mi, S, Mo <: AbstractSemSingle, O}) = χ²(sem_fit) - 2df(sem_fit)

############################################################################
### based on χ² + 2q
############################################################################

# AICq(sem_fit::SemFit{Mi, S, Mo, O} where {Mi, S, Mo <: AbstractSemSingle, O}) = χ²(sem_fit) + 2npar(sem_fit)


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))

using Printf, Dates, Query, DataFrames, CSVFiles, VegaLite
include("common.jl")

tests_to_run = [
    :textparse,
    :csvfiles,
    :textparse06,
    :csv_single,
    :csv_parallel,
    :csv06,
    :csvreader,
    # :tablereader, # Removed for now because the package is broken on Julia 1.3
    :pandas,
    :rfreads,
    :rfreadp,
    :rreadr,
    :dataframes,
    :pythonpandas,
    :pythonarrows,
    :pythonarrowp
]

runid = "master"

jl06bin = if Sys.iswindows()
    `$(joinpath(homedir(), "AppData", "Local", "Julia-0.6.4", "bin", "julia.exe"))`
elseif Sys.isapple()
    `/Applications/Julia-0.6.app/Contents/Resources/julia/bin/julia`
else
    `julia-0.6`
end

jlbin = joinpath(Sys.BINDIR, Base.julia_exename())

rbin = haskey(ENV, "R_HOME") ? joinpath(ENV["R_HOME"], "bin", Sys.iswindows() ? "RScript.exe" : Sys.isapple() ? "RScript" : "Rscript") : Sys.iswindows() ? "RScript.exe" : Sys.isapple() ? "RScript" : "Rscript"

pythonbin = "python"

platform = Sys.iswindows() ? "Windows" : Sys.isapple() ? "MacOS" : Sys.islinux() ? "Linux" : "Unknown"

if isfile(joinpath(@__DIR__), "local_config.jl")
    include("local_config.jl")
end

experiment_date = now()

output_folder_name = joinpath(@__DIR__, "..", "output", replace("run_$experiment_date", ":" => "_"))
mkpath(output_folder_name)

function run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, runtime, packagename, filename_for_label, script_filename)
    error_log_path = joinpath(output_folder_name, "errs.txt")

    open(error_log_path, "a") do f
        println(f)
        println(f, "NOW DOING $script_filename with runid=$runid, rows=$rows, cols=$cols, withna=$withna, filename=$filename, warmup_filename=$warmup_filename, samples=$samples, runtime=$runtime, packagename=$packagename, filename_for_label=$filename_for_label")
        println(f)
    end

    try
        for i in 1:samples
            if Sys.iswindows()
                path_to_esl = joinpath(@__DIR__, "..", "deps", "EmptyStandbyList.exe")
                run(`$path_to_esl`)
            elseif Sys.isapple()
                run(`sudo purge`)
            elseif Sys.islinux()
                run(`sudo sh -c 'echo 3 > /proc/sys/vm/drop_caches'`)
            end
            t1 = 0.0
            t2 = 0.0
            t3 = 0.0
            bytes1 = 0.0
            bytes2 = 0.0
            bytes3 = 0.0
            if runtime==:julia_1_single
                script_path = joinpath(@__DIR__, "packagescripts", script_filename)
                proj_env_path = joinpath(@__DIR__, "..")
                timings_as_string = readlines(pipeline(`$jlbin --startup-file=no --project=$proj_env_path $script_path $warmup_filename $filename`, stderr=error_log_path, append=true))
                t1 = parse(Float64, timings_as_string[1])
                t2 = parse(Float64, timings_as_string[2])
                t3 = parse(Float64, timings_as_string[3])
                bytes1 = parse(Float64, timings_as_string[4])
                bytes2 = parse(Float64, timings_as_string[5])
                bytes3 = parse(Float64, timings_as_string[6])
            elseif runtime==:julia_1_parallel
                script_path = joinpath(@__DIR__, "packagescripts", script_filename)
                proj_env_path = joinpath(@__DIR__, "..")
                timings_as_string = withenv("JULIA_NUM_THREADS"=>Sys.CPU_THREADS) do
                     return readlines(pipeline(`$jlbin --startup-file=no --project=$proj_env_path $script_path $warmup_filename $filename`, stderr=error_log_path, append=true))
                end
                t1 = parse(Float64, timings_as_string[1])
                t2 = parse(Float64, timings_as_string[2])
                t3 = parse(Float64, timings_as_string[3])
                bytes1 = parse(Float64, timings_as_string[4])
                bytes2 = parse(Float64, timings_as_string[5])
                bytes3 = parse(Float64, timings_as_string[6])
            elseif runtime==:julia_0_6
                script_path = joinpath(@__DIR__, "packagescripts", script_filename)
                timings_as_string = readlines(pipeline(`$jl06bin --startup-file=no $script_path $warmup_filename $filename`, stderr=error_log_path, append=true))
                t1 = parse(Float64, timings_as_string[1])
                t2 = parse(Float64, timings_as_string[2])
                t3 = parse(Float64, timings_as_string[3])
                bytes1 = parse(Float64, timings_as_string[4])
                bytes2 = parse(Float64, timings_as_string[5])
                bytes3 = parse(Float64, timings_as_string[6])
            elseif runtime==:r_project
                script_path = joinpath(@__DIR__, "packagescripts", script_filename)
                timings_as_string = readlines(pipeline(`$rbin $script_path $warmup_filename $filename`, stderr=error_log_path, append=true))
                t1 = parse(Float64, timings_as_string[1])
                t2 = parse(Float64, timings_as_string[2])
                t3 = parse(Float64, timings_as_string[3])
                bytes1 = parse(Float64, timings_as_string[4])
                bytes2 = parse(Float64, timings_as_string[5])
                bytes3 = parse(Float64, timings_as_string[6])
            elseif runtime==:python
                script_path = joinpath(@__DIR__, "packagescripts", script_filename)
                timings_as_string = readlines(pipeline(`$pythonbin $script_path $warmup_filename $filename`, stderr=error_log_path, append=true))
                t1 = parse(Float64, timings_as_string[1])
                t2 = parse(Float64, timings_as_string[2])
                t3 = parse(Float64, timings_as_string[3])
                bytes1 = parse(Float64, timings_as_string[4])
                bytes2 = parse(Float64, timings_as_string[5])
                bytes3 = parse(Float64, timings_as_string[6])
            else
                error("Incorrect julia version specified.")
            end
            push!(df, (runid=runid, attempt = :warmup, file=filename_for_label, rows=rows, cols=cols, withna=withna, package=packagename, sample=i, timing=t1, bytes=bytes1))
            push!(df, (runid=runid, attempt = :first, file=filename_for_label, rows=rows, cols=cols, withna=withna, package=packagename, sample=i, timing=t2, bytes=bytes2))
            push!(df, (runid=runid, attempt = :second, file=filename_for_label, rows=rows, cols=cols, withna=withna, package=packagename, sample=i, timing=t3, bytes=bytes3))
        end
    catch err
        @info err
    end

end

function read_specific_file(df, runid, rows, cols, withna, filename, samples)
    filename_for_label = filename
    warmup_filename = joinpath(ourpath(rows, cols, withna), "warmup_" * filename)
    filename = joinpath(ourpath(rows, cols, withna), filename)

    if :textparse in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :julia_1_single, "Julia TextParse.jl", filename_for_label, "textparse_script.jl")
    end

    if :csvfiles in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :julia_1_single, "Julia CSVFiles.jl", filename_for_label, "csvfiles_script.jl")
    end

    if :textparse06 in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :julia_0_6, "Julia v0.6 TextParse.jl", filename_for_label, "textparse_julia06_script.jl")
    end

    if :csv06 in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :julia_0_6, "Julia v0.6 CSV.jl", filename_for_label, "csv_julia06_script.jl")
    end

    if :csv_single in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :julia_1_single, "Julia CSV.jl", filename_for_label, "csv_s_script.jl")
    end

    if :csv_parallel in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :julia_1_parallel, "Julia CSV.jl parallel", filename_for_label, "csv_p_script.jl")
    end

    if :dataframes in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :julia_1_single, "Julia DataFrames.jl", filename_for_label, "dataframes_script.jl")
    end

    if :tablereader in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :julia_1_single, "Julia TableReader.jl", filename_for_label, "tablereader_script.jl")
    end

    if :csvreader in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :julia_1_single, "Julia CSVReader.jl", filename_for_label, "csvreader_script.jl")
    end

    if :pandas in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :julia_1_single, "Julia Pandas.jl", filename_for_label, "pandas_script.jl")
    end

    if :rfreads in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :r_project, "R fread", filename_for_label, "rfreads_script.R")
    end

    if :rfreadp in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :r_project, "R fread parallel", filename_for_label, "rfreadp_script.R")
    end

    if :rreadr in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :r_project, "R readr", filename_for_label, "rreadr_script.R")
    end

    if :pythonpandas in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :python, "Python Pandas", filename_for_label, "python_pandas_script.py")
    end

    if :pythonarrowp in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :python, "Python Arrow parallel", filename_for_label, "python_arrowp_script.py")
    end

    if :pythonarrows in tests_to_run
        run_script(df, runid, rows, cols, withna, filename, warmup_filename, samples, :python, "Python Arrow", filename_for_label, "python_arrows_script.py")
    end

    nothing
end

df = DataFrame(runid=String[], attempt=Symbol[], file=String[], rows=Int[], cols=Int[], withna=Bool[], package=String[], sample=Int[], timing=Float64[], bytes=Float64[])

for n in ns, c in cs, withna in [true, false]
    for filename in filter(i->endswith(i, ".csv") && !startswith(i, "warmup"), readdir(ourpath(n, c, withna)))
        @info "$filename, n=$n, c=$c, withna=$withna"
        read_specific_file(df, runid, n, c, withna, filename, samples)
    end
end

df[!, :platform] .= platform
df[!, :experiment_date] .= experiment_date

df |> save(joinpath(output_folder_name, "csvreader.csv"))

df_versions = DataFrame(package=String[], version=String[])
pkg_ctx = Pkg.Types.Context()
for p in ["CSV", "CSVFiles", "DataFrames", "TextParse", "CSVReader", "Pandas", "TableReader"]
    if haskey(pkg_ctx.env.project.deps, p)
        push!(df_versions, ("Julia " * p * ".jl", string(pkg_ctx.env.manifest[pkg_ctx.env.project.deps[p]].version)))
    end
end
push!(df_versions, ("Julia v0.6 CSV.jl", "0.2.5"))
push!(df_versions, ("Julia v0.6 TextParse.jl", "0.5.0"))

try
    vers = readlines(pipeline(`$rbin $(joinpath("packagescripts", "rfread_version.R"))`, stderr=error_log_path, append=true))[1]
    push!(df_versions, ("R fread", vers))
catch err
end

try
    vers = readlines(pipeline(`$rbin $(joinpath("packagescripts", "rreadr_version.R"))`, stderr=error_log_path, append=true))[1]
    push!(df_versions, ("R readr", vers))
catch err
end

try
    vers = readlines(pipeline(`$pythonbin $(joinpath("packagescripts", "python_arrow_version.py"))`, stderr=error_log_path, append=true))[1]
    push!(df_versions, ("Python Arrow", vers))
catch err
end

try
    vers = readlines(pipeline(`$pythonbin $(joinpath("packagescripts", "python_pandas_version.py"))`, stderr=error_log_path, append=true))[1]
    push!(df_versions, ("Python Pandas", vers))
catch err
end

df_versions |> save(joinpath(output_folder_name, "csvreaderpackageversions.csv"))

for c in cs
    p1 = df |>
    @filter(_.attempt==:first) |>
    @filter(!_.withna && _.cols == c) |>
    @vlplot(facet={row={field=:file, typ=:nominal, title=nothing}, column={field=:rows, typ=:ordinal}}, resolve={scale={x=:independent}}, title="Without missing data ($c columns)", background=:white) +
        (
            @vlplot(y={"package:n", title=nothing}, resolve={scale={x=:shared}}) +
            @vlplot(:bar,  x={"min(timing):q", title="seconds"}, color={"package:n", scale={scheme=:category20}}) +
            @vlplot(:rule, x="min(timing):q", x2="max(timing):q") +
            @vlplot(:tick, x="timing:q")
        )
    p1 |> save(joinpath(output_folder_name, "cols_$(c)_withoutna.pdf"))
    p1 |> save(joinpath(output_folder_name, "cols_$(c)_withoutna.png"))

    p2 = df |>
    @filter(_.attempt==:first) |>
    @filter(_.withna && _.cols == c) |>
    @vlplot(facet={row={field=:file, typ=:nominal, title=nothing}, column={field=:rows, typ=:ordinal}}, resolve={scale={x=:independent}}, title="With missing data ($c columns)", background=:white) +
        (
            @vlplot(y={"package:n", title=nothing}, resolve={scale={x=:shared}}) +
            @vlplot(:bar,  x={"min(timing):q", title="seconds"}, color={"package:n", scale={scheme=:category20}}) +
            @vlplot(:rule, x="min(timing):q", x2="max(timing):q") +
            @vlplot(:tick, x="timing:q")
        )
    p2 |> save(joinpath(output_folder_name, "cols_$(c)_withna.pdf"))
    p2 |> save(joinpath(output_folder_name, "cols_$(c)_withna.png"))
end

max_stars_repo_name: davidanthoff/csv-comparison

max_issues_repo_name: davidanthoff/csv-comparison

max_forks_repo_name: davidanthoff/csv-comparison