text:
const CPetscOptions = Ptr{Cvoid}

"""
    AbstractOptions{PetscLib <: PetscLibType}

Abstract type of PETSc solver options.
"""
abstract type AbstractOptions{PetscLib <: PetscLibType} end

"""
    GlobalOptions{PetscLib <: PetscLibType}

The PETSc global options database.
"""
struct GlobalOptions{PetscLib} <: AbstractOptions{PetscLib} end
Base.cconvert(::Type{CPetscOptions}, obj::GlobalOptions) = C_NULL

"""
    Options{PetscLib <: PetscLibType}(kw -> arg, ...)
    Options(petsclib, kw -> arg, ...)

Create a PETSc options data structure for the `petsclib`.

For construction a set of keyword argment pairs should be given. If the option
has no value it should be set to `nothing`

# Examples
```julia-repl
julia> opt = PETScX.Options(
           petsclib,
           ksp_monitor = true,
           ksp_view = true,
           pc_type = "mg",
           pc_mg_levels = 1,
       )
#PETSc Option Table entries:
-ksp_monitor true
-ksp_view true
-pc_mg_levels 1
-pc_type mg
#End of PETSc Option Table entries

julia> opt["ksp_monitor"]
"true"

julia> opt["pc_type"]
"mg"

julia> opt["pc_type"] = "ilu"
"ilu"

julia> opt["pc_type"]
"ilu"

julia> opt["bad_key"]
ERROR: KeyError: key "bad_key" not found
```

Manual: [`PetscOptionsCreate`](https://petsc.org/release/docs/manualpages/Sys/PetscOptionsCreate.html)
"""
mutable struct Options{T} <: AbstractOptions{T}
    ptr::CPetscOptions
end
Base.cconvert(::Type{CPetscOptions}, obj::Options) = obj.ptr
Base.unsafe_convert(::Type{Ptr{CPetscOptions}}, obj::Options) =
    convert(Ptr{CPetscOptions}, pointer_from_objref(obj))

Options(petsclib; kwargs...) = Options(petsclib, kwargs...)
function Options(petsclib, ps::Pair...)
    opts = Options(petsclib)
    for (k, v) in ps
        opts[k] = v
    end
    return opts
end

@for_libpetsc function Options(::$UnionPetscLib)
    opts = Options{$PetscLib}(C_NULL)
    @assert Initialized($PetscLib)
    @chk ccall(
        (:PetscOptionsCreate, $petsc_library),
        PetscErrorCode,
        (Ptr{CPetscOptions},),
        opts,
    )
    finalizer(Finalize, opts)
    return opts
end

@for_libpetsc function Finalize(opts::Options{$PetscLib})
    Finalized($PetscLib) || @chk ccall(
        (:PetscOptionsDestroy, $petsc_library),
        PetscErrorCode,
        (Ptr{CPetscOptions},),
        opts,
    )
    return nothing
end

@for_libpetsc function Base.push!(
    ::GlobalOptions{$PetscLib},
    opts::Options{$PetscLib},
)
    @chk ccall(
        (:PetscOptionsPush, $petsc_library),
        PetscErrorCode,
        (CPetscOptions,),
        opts,
    )
    return nothing
end

@for_libpetsc function Base.pop!(::GlobalOptions{$PetscLib})
    @chk ccall((:PetscOptionsPop, $petsc_library), PetscErrorCode, ())
    return nothing
end

@for_libpetsc function Base.setindex!(
    opts::AbstractOptions{$PetscLib},
    val,
    key,
)
    @chk ccall(
        (:PetscOptionsSetValue, $petsc_library),
        PetscErrorCode,
        (CPetscOptions, Cstring, Cstring),
        opts,
        string('-', key),
        isnothing(val) ? C_NULL : string(val),
    )
end

@for_libpetsc function Base.getindex(opts::AbstractOptions{$PetscLib}, key)
    val = Vector{UInt8}(undef, 256)
    set_ref = Ref{PetscBool}()
    @chk ccall(
        (:PetscOptionsGetString, $petsc_library),
        PetscErrorCode,
        (CPetscOptions, Cstring, Cstring, Ptr{UInt8}, Csize_t, Ptr{PetscBool}),
        opts,
        C_NULL,
        string('-', key),
        val,
        sizeof(val),
        set_ref,
    )
    val = GC.@preserve val unsafe_string(pointer(val))
    set_ref[] == PETSC_TRUE || throw(KeyError(key))
    return val
end

@for_libpetsc function view(
    opts::AbstractOptions{$PetscLib},
    viewer::AbstractViewer{$PetscLib} = ViewerStdout($PetscLib, MPI.COMM_SELF),
)
    @chk ccall(
        (:PetscOptionsView, $petsc_library),
        PetscErrorCode,
        (CPetscOptions, CPetscViewer),
        opts,
        viewer,
    )
    return nothing
end

@for_libpetsc GlobalOptions(::$UnionPetscLib) = GlobalOptions{$PetscLib}()

Base.show(io::IO, opts::AbstractOptions) = _show(io, opts)

"""
    with(f, opts::Options)

Call `f()` with the [`Options`](@ref) `opts` set temporarily (in addition to any
global options).
"""
function with(f, opts::Options{T}) where {T}
    global_opts = GlobalOptions{T}()
    push!(global_opts, opts)
    try
        f()
    finally
        pop!(global_opts)
    end
end


text:
using Memento

function set_logger(game::Game)
    logger = getlogger(current_module())
    logger.handlers = Dict{Any,Memento.Handler}()   # Remove console handler
    push!(logger, DefaultHandler(string(game.botname, "_", game.id, ".log")))

    setlevel!(logger, "info")

    for i in values(logger.handlers)
        i.fmt = DefaultFormatter("{msg}")   # Rewrite logging message format
    end

    return logger
end

ext: jl

max_stars_repo_name: Cerbes/Halite-II

max_issues_repo_name: Cerbes/Halite-II

max_forks_repo_name: Cerbes/Halite-II


text:
using Random
using OrdinaryDiffEq: Tsit5
using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Plots, AdvancedHMC, MCMCChains
using JLD, StatsPlots

# ----- START: OPTIMIZATION TEST FUNCTIONS ----- #
# glucose, steps, bolus, meals
function optim_all_vars_no_exo(bin_size::Int64, batch_size::Int=4, maxiters::Int=150, lr::Float64=0.01; hypo_id::Int64=1, log_dir="")
  u0, ode_data = prepare_all_data_meals_hypo(bin_size, hypo_id=hypo_id)
  _, _, _, pl = model(ode_data, batch_size, maxiters, lr, log_dir=log_dir)
  pl
end

# glucose, steps, bolus
function optim_two_vars_no_exo(bin_size::Int64, batch_size::Int=4, maxiters::Int=150, lr::Float64=0.01)
  u0, ode_data = prepare_all_data(bin_size)
  samples, losses, best_pred = model(ode_data, batch_size, maxiters, lr)
  samples, losses, best_pred
end

# (glucose, steps) OR (glucose, bolus)
function optim_one_var_no_exo(var::String, bin_size::Int64, batch_size::Int=4, maxiters::Int=150, lr::Float64=0.01)
  u0, ode_data = prepare_data(var, bin_size)
  samples, losses, best_pred = model(ode_data, batch_size, maxiters, lr)
  samples, losses, best_predb
end

# glucose, steps, bolus, meals
function optim_all_vars_exo(bin_size::Int64, batch_size::Int=4, maxiters::Int=150, lr::Float64=0.01; hypo_id::Int64=2, log_dir="")
  @show hypo_id
  u0, ode_data = prepare_all_data_meals_hypo(bin_size, hypo_id=hypo_id)
  non_exo_data = ode_data[1:2, :]
  exo_data = ode_data[3:4, :]
  # samples, losses, best_pred = model_exo(non_exo_data, exo_data, batch_size, maxiters, lr)
  model_exo(non_exo_data, exo_data, batch_size, maxiters, lr, log_dir=log_dir)
end

# glucose, steps, bolus
function optim_two_vars_exo(bin_size::Int64, batch_size::Int=4, maxiters::Int=150, lr::Float64=0.01)
  u0, ode_data = prepare_all_data(bin_size)
  non_exo_data = ode_data[1:2, :]
  exo_data = ode_data[3, :]
  samples, losses, best_pred = model_exo(non_exo_data, exo_data, batch_size, maxiters, lr)
end

# (glucose, steps) OR (glucose, bolus)
function optim_one_var_exo(var::String, bin_size::Int64, batch_size::Int=4, maxiters::Int=150, lr::Float64=0.01)
  u0, ode_data = prepare_data(var, bin_size)
  if var == "basis_steps"
    samples, losses, best_pred = model(ode_data, batch_size, maxiters, lr)
  else
    non_exo_data = ode_data[1, :]
    exo_data = ode_data[2, :]
    samples, losses, best_pred = model_exo(non_exo_data, exo_data, batch_size, maxiters, lr)
  end
end

# ----- END: OPTIMIZATION TEST FUNCTIONS ----- #

# ----- START: NEURAL BAYES TEST FUNCTIONS ----- #

# glucose, steps, bolus, meals
function bayes_all_vars_no_exo(bin_size::Int64;
                               numwarmup::Int = 500,
                               numsamples::Int = 500,
                               initstepsize::Float64 = 0.45,
                               odesolver = Tsit5,
                               log_dir="")
  u0, ode_data = prepare_all_data_meals_hypo(bin_size)
  ode_data = ode_data[1:4, :]
  samples, losses, predict_neuralode = model_bayes(ode_data, numwarmup, numsamples, initstepsize, odesolver, log_dir=log_dir)
  tsteps = range(0.0, 1.0, length=size(ode_data)[2])

  ################### RETRODICTED PLOTS: TIME SERIES #################
  pl = plot(tsteps, ode_data[1,:], color = :red, label = "Data: CGM", xlabel = "t", title = "CGM, Steps, and Bolus")
  plot!(tsteps, ode_data[2,:], color = :blue, label = "Data: Steps")
  plot!(tsteps, ode_data[3,:], color = :green, label = "Data: Bolus")
  plot!(tsteps, ode_data[4,:], color = :purple, label = "Data: Meals")

  for k in 1:300
    resol = predict_neuralode(samples[400:end][rand(1:100)])
    plot!(tsteps,resol[1,:], alpha=0.1, color = :red, label = "")
    plot!(tsteps,resol[2,:], alpha=0.1, color = :blue, label = "")
    plot!(tsteps,resol[3,:], alpha=0.1, color = :green, label = "")
    plot!(tsteps,resol[4,:], alpha=0.1, color = :purple, label = "")
  end

  idx = findmin(losses)[2]
  prediction = predict_neuralode(samples[idx])

  plot!(tsteps, prediction[1,:], color = :black, w = 2, label = "C")
  plot!(tsteps, prediction[2,:], color = :black, w = 2, label = "S")
  plot!(tsteps, prediction[3,:], color = :black, w = 2, label = "B")
  plot!(tsteps, prediction[4,:], color = :black, w = 2, label = "M")
  savefig(pl, "bayes_full_data_no_exo_bin_$(bin_size)_i_$(rand(1:100)).png")
end
# glucose, steps, bolus
function bayes_two_vars_no_exo(bin_size::Int64)
  u0, ode_data = prepare_all_data(bin_size)
  samples, losses = model_bayes(ode_data)
  samples, losses
end

# (glucose, steps) OR (glucose, bolus)
function bayes_one_var_no_exo(var::String)
  u0, ode_data = prepare_data(var, bin_size)
  samples, losses = model_bayes(ode_data)
  samples, losses
end

# glucose, steps, bolus, meals
function bayes_all_vars_exo(bin_size::Int64;
                            numwarmup::Int = 500,
                            numsamples::Int = 500,
                            initstepsize::Float64 = 0.45,
                            odesolver = Tsit5,
                            log_dir="")
  u0, ode_data = prepare_all_data_meals_hypo(bin_size)
  non_exo_data = ode_data[1:2, :]
  exo_data = ode_data[3:4, :]
  samples, losses, predict_neuralode = model_bayes_exo(non_exo_data, exo_data, numwarmup, numsamples, initstepsize, odesolver, log_dir=log_dir)
  tsteps = range(0.0, 1.0, length=size(ode_data)[2])
  ################### RETRODICTED PLOTS: TIME SERIES #################
  pl = plot(tsteps, ode_data[1,:], color = :red, label = "Data: CGM", xlabel = "t", title = "CGM, Steps, and Bolus")
  plot!(tsteps, ode_data[2,:], color = :blue, label = "Data: Steps")
  plot!(tsteps, ode_data[3,:], color = :green, label = "Data: Bolus")
  plot!(tsteps, ode_data[4,:], color = :purple, label = "Data: Meals")

  for k in 1:300
    resol = predict_neuralode(samples[400:end][rand(1:100)])
    plot!(tsteps,resol[1,1,:], alpha=0.4, color = :red, label = "")
    plot!(tsteps,resol[2,1,:], alpha=0.4, color = :blue, label = "")
  end

  idx = findmin(losses)[2]
  prediction = predict_neuralode(samples[idx])

  plot!(tsteps, prediction[1,:], color = :black, w = 2, label = "C")
  plot!(tsteps, prediction[2,:], color = :black, w = 2, label = "S")
  savefig(pl, "bayes_full_data_exo_bin_$(bin_size)_i_$(rand(1:100)).png")
end

function bayes_all_vars_exo_with_init(bin_size::Int64;
                                      numwarmup::Int = 500,
                                      numsamples::Int = 500,
                                      initstepsize::Float64 = 0.45,
                                      odesolver = Tsit5,
                                      log_dir="",
				                              initid=81)
  u0, ode_data = prepare_all_data_meals_hypo(bin_size)
  non_exo_data = ode_data[1:2, :]
  exo_data = ode_data[3:4, :]
  samples, losses, predict_neuralode, num = model_bayes_exo_with_init(non_exo_data, exo_data, numwarmup, numsamples, initstepsize, odesolver, log_dir=log_dir, initid=initid)
  tsteps = range(0.0, Float64(size(ode_data)[2]) - 1, length=size(ode_data)[2])
  ################### RETRODICTED PLOTS: TIME SERIES #################
  pl = plot(tsteps, ode_data[1,:], color = :red, label = "Data: CGM", xlabel = "t", title = "CGM, Steps, and Bolus")
  plot!(tsteps, ode_data[2,:], color = :blue, label = "Data: Steps")
  plot!(tsteps, ode_data[3,:], color = :green, label = "Data: Bolus")
  plot!(tsteps, ode_data[4,:], color = :purple, label = "Data: Meals")

  for k in 1:300
    resol = predict_neuralode(samples[Int(round(numsamples * 0.8)):end][rand(1:(Int(round(numsamples * 0.2)) - 1))])
    if (abs(maximum(resol[1,:])) < 30)
      plot!(tsteps,resol[1,:], alpha=0.1, color = :red, label = "")
    end
    if (abs(maximum(resol[2,:])) < 30)
      plot!(tsteps,resol[2,:], alpha=0.1, color = :red, label = "")
    end
  end

  idx = findmin(losses)[2]
  prediction = predict_neuralode(samples[idx])

  plot!(tsteps, prediction[1,:], color = :black, w = 2, label = "C")
  plot!(tsteps, prediction[2,:], color = :black, w = 2, label = "S")
  savefig(pl, joinpath(log_dir, "bayes_full_data_exo_with_init_bin_$(bin_size)_initid_$(initid)_i_$(num).png"))
end

# glucose, steps, bolus
function bayes_two_vars_exo(bin_size::Int64)
  u0, ode_data = prepare_all_data(bin_size)
  non_exo_data = ode_data[1:2, :]
  exo_data = ode_data[3, :]
  samples, losses = model_bayes_exo(non_exo_data, exo_data)
end

# (glucose, steps) OR (glucose, bolus)
function bayes_one_var_exo(var::String)
  u0, ode_data = prepare_data(var, bin_size)
  if var == "basis_steps"
    samples, losses = model_bayes_exo(ode_data)
  else
    non_exo_data = ode_data[1, :]
    exo_data = ode_data[2, :]
    samples, losses = model_bayes_exo(non_exo_data, exo_data)
  end
end

# ----- END: NEURAL BAYES TEST FUNCTIONS ----- #

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using Base: @deprecate, depwarn
import Dates: Date, DateTime

# BEGIN Intervals 1.X.Y deprecations

export Inclusivity, inclusivity
include("inclusivity.jl")

@deprecate Date(interval::Interval{Date}) convert(Date, interval)
@deprecate DateTime(interval::Interval{DateTime}) convert(DateTime, interval)
@deprecate Date(interval::AnchoredInterval{P, Date} where P) convert(Date, interval)
@deprecate DateTime(interval::AnchoredInterval{P, DateTime} where P) convert(DateTime, interval)

function Endpoint{T,D}(ep::T, included::Bool) where {T,D}
    B = bound_type(included)
    depwarn("`Endpoint{T,D}(ep, $included)` is deprecated, use `Endpoint{T,D,$(repr(B))}(ep)` instead.", :Endpoint)
    return Endpoint{T,D,B}(ep)
end

function LeftEndpoint(ep, included::Bool)
    B = bound_type(included)
    depwarn("`LeftEndpoint(ep, $included)` is deprecated, use `LeftEndpoint{$(repr(B))}(ep)` instead.", :LeftEndpoint)
    return LeftEndpoint{B}(ep)
end

function RightEndpoint(ep, included::Bool)
    B = bound_type(included)
    depwarn("`RightEndpoint(ep, $included)` is deprecated, use `RightEndpoint{$(repr(B))}(ep)` instead.", :RightEndpoint)
    return RightEndpoint{B}(ep)
end

# intervals.jl
function Interval{T,L,R}(f::T, l::T, inc::Inclusivity) where {T,L,R}
    left_inc = bound_type(first(inc))
    right_inc = bound_type(last(inc))
    if L !== left_inc || R !== right_inc
        throw(ArgumentError("Specified inclusivity ($(repr(left_inc)), $(repr(right_inc))) doesn't match bound types ($(repr(L)), $(repr(R)))"))
    end
    depwarn("`Interval{T,$(repr(L)),$(repr(R))}(f, l, $(repr(inc)))` is deprecated, use `Interval{T,$(repr(L)),$(repr(R))}(f, l)` instead.", :Interval)
    return Interval{T,L,R}(f, l)
end

function Interval{T}(f, l, inc::Inclusivity) where T
    L = bound_type(first(inc))
    R = bound_type(last(inc))
    depwarn("`Interval{T}(f, l, $(repr(inc)))` is deprecated, use `Interval{T,$(repr(L)),$(repr(R))}(f, l)` instead.", :Interval)
    return Interval{T,L,R}(f, l)
end

function Interval{T}(f, l, x::Bool, y::Bool) where T
    L = bound_type(x)
    R = bound_type(y)
    depwarn("`Interval{T}(f, l, $x, $y)` is deprecated, use `Interval{T,$(repr(L)),$(repr(R))}(f, l)` instead.", :Interval)
    return Interval{T,L,R}(f, l)
end

function Interval(f, l, inc::Inclusivity)
    L = bound_type(first(inc))
    R = bound_type(last(inc))
    depwarn("`Interval(f, l, $(repr(inc)))` is deprecated, use `Interval{T,$(repr(L)),$(repr(R))}(f, l)` instead.", :Interval)
    return Interval{L,R}(f, l)
end

function Interval(f, l, x::Bool, y::Bool)
    L = bound_type(x)
    R = bound_type(y)
    depwarn("`Interval(f, l, $x, $y)` is deprecated, use `Interval{T,$(repr(L)),$(repr(R))}(f, l)` instead.", :Interval)
    return Interval{L,R}(f, l)
end

function inclusivity(interval::AbstractInterval{T,L,R}) where {T,L,R}
    depwarn("`inclusivity(interval)` is deprecated and has no direct replacement. See `bounds_types(interval)` for similar functionality.", :inclusivity)
    return Inclusivity(L === Closed, R === Closed; ignore_depwarn=true)
end

# anchoredintervals.jl
function AnchoredInterval{P,T,L,R}(anchor::T, inc::Inclusivity) where {P,T,L,R}
    left_inc = bound_type(first(inc))
    right_inc = bound_type(last(inc))
    if L !== left_inc || R !== right_inc
        throw(ArgumentError("Specified inclusivity ($(repr(left_inc)), $(repr(right_inc))) doesn't match bound types ($(repr(L)), $(repr(R)))"))
    end
    depwarn("`AnchoredInterval{P,T,$(repr(L)),$(repr(R))}(anchor, $(repr(inc)))` is deprecated, use `AnchoredInterval{P,T,$(repr(L)),$(repr(R))}(anchor)` instead.", :AnchoredInterval)
    return AnchoredInterval{P,T,L,R}(anchor)
end

function AnchoredInterval{P,T}(anchor, inc::Inclusivity) where {P,T}
    L = bound_type(first(inc))
    R = bound_type(last(inc))
    depwarn("`AnchoredInterval{P,T}(anchor, $(repr(inc)))` is deprecated, use `AnchoredInterval{P,T,$(repr(L)),$(repr(R))}(anchor)` instead.", :AnchoredInterval)
    return AnchoredInterval{P,T,L,R}(anchor)
end

function AnchoredInterval{P,T}(anchor, x::Bool, y::Bool) where {P,T}
    L = bound_type(x)
    R = bound_type(y)
    depwarn("`AnchoredInterval{P,T}(anchor, $x, $y)` is deprecated, use `AnchoredInterval{P,T,$(repr(L)),$(repr(R))}(anchor)` instead.", :AnchoredInterval)
    return AnchoredInterval{P,T,L,R}(anchor)
end

function AnchoredInterval{P}(anchor, inc::Inclusivity) where P
    L = bound_type(first(inc))
    R = bound_type(last(inc))
    depwarn("`AnchoredInterval{P}(anchor, $(repr(inc)))` is deprecated, use `AnchoredInterval{P,$(repr(L)),$(repr(R))}(anchor)` instead.", :AnchoredInterval)
    return AnchoredInterval{P,L,R}(anchor)
end

function AnchoredInterval{P}(anchor, x::Bool, y::Bool) where P
    L = bound_type(x)
    R = bound_type(y)
    depwarn("`AnchoredInterval{P}(anchor, $x, $y)` is deprecated, use `AnchoredInterval{P,$(repr(L)),$(repr(R))}(anchor)` instead.", :AnchoredInterval)
    return AnchoredInterval{P,L,R}(anchor)
end

function HourEnding(anchor, x::Bool, y::Bool)
    L = bound_type(x)
    R = bound_type(y)
    depwarn("`HourEnding(anchor, $x, $y)` is deprecated, use `HourEnding{$(repr(L)),$(repr(R))}(anchor)` instead.", :HourEnding)
    return HourEnding{L,R}(anchor)
end

function HourEnding(anchor, inc::Inclusivity)
    L = bound_type(first(inc))
    R = bound_type(last(inc))
    depwarn("`HourEnding(anchor, $(repr(inc)))` is deprecated, use `HourEnding{$(repr(L)),$(repr(R))}(anchor)` instead.", :HourEnding)
    return HourEnding{L,R}(anchor)
end

function HourBeginning(anchor, x::Bool, y::Bool)
    L = bound_type(x)
    R = bound_type(y)
    depwarn("`HourBeginning(anchor, $x, $y)` is deprecated, use `HourBeginning{$(repr(L)),$(repr(R))}(anchor)` instead.", :HourBeginning)
    return HourBeginning{L,R}(anchor)
end

function HourBeginning(anchor, inc::Inclusivity)
    L = bound_type(first(inc))
    R = bound_type(last(inc))
    depwarn("`HourBeginning(anchor, $(repr(inc)))` is deprecated, use `HourBeginning{$(repr(L)),$(repr(R))}(anchor)` instead.", :HourBeginning)
    return HourBeginning{L,R}(anchor)
end

function HE(anchor, x::Bool, y::Bool)
    L = bound_type(x)
    R = bound_type(y)
    if !x && y
        depwarn("`HE(anchor, $x, $y)` is deprecated, use `HE(anchor)` instead.", :HE)
    else
        depwarn("`HE(anchor, $x, $y)` is deprecated, use `HourEnding{$(repr(L)),$(repr(R))}(ceil(anchor, Hour))` instead.", :HE)
    end
    return HourEnding{L,R}(ceil(anchor, Hour))
end

function HE(anchor, inc::Inclusivity)
    L = bound_type(first(inc))
    R = bound_type(last(inc))
    if !first(inc) && last(inc)
        depwarn("`HE(anchor, $(repr(inc)))` is deprecated, use `HE(anchor)` instead.", :HE)
    else
        depwarn("`HE(anchor, $(repr(inc)))` is deprecated, use `HourEnding{$(repr(L)),$(repr(R))}(ceil(anchor, Hour))` instead.", :HE)
    end
    return HourEnding{L,R}(ceil(anchor, Hour))
end

function HB(anchor, x::Bool, y::Bool)
    L = bound_type(x)
    R = bound_type(y)
    if x && !y
        depwarn("`HB(anchor, $x, $y)` is deprecated, use `HB(anchor)` instead.", :HB)
    else
        depwarn("`HB(anchor, $x, $y)` is deprecated, use `HourBeginning{$(repr(L)),$(repr(R))}(floor(anchor, Hour))` instead.", :HB)
    end
    return HourBeginning{L,R}(floor(anchor, Hour))
end

function HB(anchor, inc::Inclusivity)
    L = bound_type(first(inc))
    R = bound_type(last(inc))
    if first(inc) && !last(inc)
        depwarn("`HB(anchor, $(repr(inc)))` is deprecated, use `HB(anchor)` instead.", :HB)
    else
        depwarn("`HB(anchor, $(repr(inc)))` is deprecated, use `HourBeginning{$(repr(L)),$(repr(R))}(floor(anchor, Hour))` instead.", :HB)
    end
    return HourBeginning{L,R}(floor(anchor, Hour))
end

# END Intervals 1.X.Y deprecations


text:
using RecipesBase
using StatPlots

@userplot MeanPlot
@userplot AutocorPlot
@userplot HistogramPlot
@userplot DensityPlot
@userplot MixedDensityPlot
@userplot TracePlot

@recipe function f(p::MeanPlot)
  c, i = p.args
  seriestype := :line
  xaxis := "Iteration"
  yaxis := "Mean"
  MCMCChain.cummean(c.value[:, i, :])
end

@recipe function f(p::AutocorPlot; maxlag = nothing)
  c, i = p.args
  lags = 0:(maxlag === nothing ? round(Int, 10 * log10(length(c.range))) : maxlag)
  ac = MCMCChain.autocor(c, lags=collect(lags))
  xaxis := "Lag"
  yaxis := "Autocorrelation"
  ac.value[i, :, :]
end

@recipe function f(p::HistogramPlot)
  c, i = p.args
  seriestype := :histogram
  xaxis := "Value"
  yaxis := "Value"
  c.value[:, i, :]
end

@recipe function f(p::DensityPlot)
  c, i = p.args
  seriestype := :density
  xaxis := "Value"
  yaxis := "Density"
  c.value[:, i, :]
end

@recipe function f(p::TracePlot)
  c, i = p.args
  seriestype := :line
  xaxis := "Iteration"
  yaxis := "Density"
  c.value[:, i, :]
end

@recipe function f(p::MixedDensityPlot)
  c, i = p.args
  discrete = MCMCChain.indiscretesupport(c, barbounds)
  seriestype := discrete[i] ? :histogram : :density
  xaxis := "Value"
  yaxis := "Density"
  c.value[:, i, :]
end

@recipe function f(c::MCMCChain.AbstractChains,
                   ptypes::AbstractArray;
                   barbounds = (0, Inf))
  nrows, nvars, nchains = size(c.value)
  ntypes = length(ptypes)
  layout := (nvars, ntypes)
  indices = reshape(1:nvars*ntypes, ntypes, nvars)'

  legend --> false

  for (j, ptype) in enumerate(ptypes)
    for i in 1:nvars
      @series begin
          subplot := indices[i, j]
          title := c.names[i]
          labels := map(k -> "Chain $k", 1:nchains)
          ptype([c, i])
      end
    end
  end
end


text:
# ============================================================================ #
macro libpath(prefix, dir...)
    @static if Sys.islinux()
        ext = ".so"
    elseif Sys.iswindows()
        ext = ".dll"
    elseif Sys.isapple()
        # nothing to do here
        ext = ""
    else
        ext = ".so"
    end

    return quote
        # path to lib is now relative to the containing file
        # (i.e. the file that is invoking this macro)
        joinpath($(dir...), $prefix * $ext)
    end
end
# ============================================================================ #
macro calllib(ifile, idx, chan_type, jtyp)#ctyp, jtyp)

    ctyp = Symbol("c", jtyp)
    # jtyp = Symbol("SMR" * type_name * "Channel")

    fread = "read_" * chan_type * "_channel"
    ffree = "free_" * chan_type * "_channel"

    return quote
        if splitext($(esc(ifile)))[2] != ".smr"
            error("Input file is not an smr file")
        end

        ptr = ccall(($fread, LIBSMR), Ptr{$ctyp}, (Cstring, Cint),
            $(esc(ifile)), Cint($(esc(idx))))

        if ptr != C_NULL
            out = $jtyp(unsafe_load(ptr))
            ccall(($ffree, LIBSMR), Cvoid, (Ptr{$ctyp},), ptr)
        else
            error("call to " * string($fread) * " failed")
        end
        out
    end
end
# ============================================================================ #

max_stars_repo_name: scottiealexander/libsmr

max_issues_repo_name: scottiealexander/libsmr

max_forks_repo_name: scottiealexander/libsmr


text:
include("../src/utils.jl")

@testset "utils" begin
  @test filter_handler(nothing) == nothing
  @test filter_handler(Dict("doi" => "10.1038/nature12373")) == "doi:10.1038/nature12373,"

end

max_stars_repo_name: naustica/Pitaya

max_issues_repo_name: naustica/Pitaya

max_forks_repo_name: naustica/Pitaya


text:
import Base: length, size, eltype, getindex, setindex!, similar, copy, hash,
    isequal, ==, convert, promote_rule, isapprox, show, @propagate_inbounds
import Compat.LinearAlgebra: checksquare, diag, complex, rank, normalize!,
    normalize, ishermitian, issymmetric, isdiag, triu, tril
import Compat: norm, opnorm
if VERSION < v"1.0.0-"
    import Base: full
end
# Special QuObject methods
data(A::QuObject) = A.data
dims(A::QuObject) = A.dims
isnormalized(A::QuVector) = norm(A) ≈ 1
isnormalized(A::QuMatrix) = trace(A) ≈ 1
isdensityop(A::QuMatrix) = isapproxhermitian(A) && isnormalized(A)
function checkdensityop(A::QuMatrix)
    herm = isapproxhermitian(A); norm = isnormalized(A)
    if !herm || !norm
        a = !herm ? "is not Hermitian" : ""
        b = !herm && !norm ? " and " : ""
        c = !norm ? "does not have unity trace (trace: $(real(trace(A))))" : ""
        throw(ArgumentError(string("not a density operator, matrix ",a,b,c)))
    end
end
dense(x::Ket) = Ket(Vector(x.data),x.dims)
dense(x::Bra) = Bra(Vector(x.data),x.dims)
dense(A::Operator) = Operator(Matrix(A.data),A.dims)
dimsmatch(A::QuObject,B::QuObject) =
    A.dims==B.dims || throw(DimensionMismatch("subspace dimensions do not match"))
function dimsmatch(As::AbstractVecOrTuple{QuObject},Bs::AbstractVecOrTuple{QuObject})
    for A ∈ As, B ∈ Bs; dimsmatch(A,B); end
end
dimsmatch(As::AbstractVecOrTuple{QuObject},B::QuObject) = dimsmatch(As,(B,))
dimsmatch(A::QuObject,Bs::AbstractVecOrTuple{QuObject}) = dimsmatch((A,),Bs)

# Tensored indexing methods
@propagate_inbounds getindex(A::QuVector,t::Tuple) =
    getindex(A.data,tensored_sub2ind(dims(A),t))
@propagate_inbounds getindex(A::QuMatrix,ti::Tuple,tj::Tuple) =
    getindex(A.data,tensored_sub2ind(dims(A),ti),tensored_sub2ind(dims(A),tj))
@propagate_inbounds setindex!(A::QuVector,v,t::Tuple) =
    setindex!(A.data,v,tensored_sub2ind(dims(A),t))
@propagate_inbounds setindex!(A::QuMatrix,v,ti::Tuple,tj::Tuple) =
    setindex!(A.data,v,tensored_sub2ind(dims(A),ti),tensored_sub2ind(dims(A),tj))

# Translate basic Base array methods to QuObjects
@propagate_inbounds getindex(A::QuObject,idx...) = getindex(A.data,idx...)
@propagate_inbounds setindex!(A::QuObject,v,idx...) = setindex!(A.data,v,idx...)
length(A::QuObject) = length(A.data)
size(A::QuObject,d) = size(A.data,d)
size(A::QuObject) = size(A.data)
checksquare(A::QuMatrix) = checksquare(A.data)
eltype(A::QuObject) = eltype(A.data)
diag(A::QuMatrix,k::Int=0) = diag(A.data,k)
triu(A::QuMatrix) = Operator(triu(A.data),A.dims)
tril(A::QuMatrix) = Operator(tril(A.data),A.dims)
full(A::QuObject) = Array(A.data)
complex(x::Ket) = Ket(complex(x.data),x.dims)
complex(x::Bra) = Bra(complex(x.data),x.dims)
complex(A::Operator) = Operator(complex(A.data),A.dims)
norm(x::QuVector,p::Integer=2) = norm(x.data,p)
norm(A::QuMatrix,p::Integer=2) = norm(A.data,p)
opnorm(A::QuMatrix,p::Integer=2) = opnorm(A.data,p)
trace(A::QuMatrix) = trace(A.data)
rank(A::QuMatrix) = rank(A.data)
normalize!(x::QuVector) = (normalize!(x.data,2);x)
normalize!(A::QuMatrix) = (rmul!(A.data,1/trace(A.data));A)
normalize(x::QuObject) = normalize!(copy(x))
scale!(A::QuObject,b::Number) = (rmul!(A.data,b);A)
scale(A::QuObject,b::Number) = scale!(copy(A))
ishermitian(A::Operator) = ishermitian(A.data)
isapproxhermitian(A::Operator) = isapproxhermitian(A.data)
issymmetric(A::QuMatrix) = issymmetric(A.data)
isunitary(A::QuMatrix) = isunitary(A.data)
isapproxunitary(A::QuMatrix) = isapproxunitary(A.data)
isdiag(A::QuMatrix) = isdiag(A.data)
similar(A::T) where {T<:QuObject}= T(similar(A.data),A.dims)
copy(A::T) where {T<:QuObject} = T(copy(A.data),A.dims)
hash(A::T,h::UInt) where {T<:QuObject} = hash(hash(A.data,hash(A.dims,hash(T))),h)
isequal(A::T,B::T) where {T<:QuObject} = isequal(A.dims,B.dims)&&isequal(A.data,B.data)
==(A::Ket,B::Ket) = isequal(A.dims,B.dims)&&(A.data==B.data)
==(A::Bra,B::Bra) = isequal(A.dims,B.dims)&&(A.data==B.data)
==(A::Operator,B::Operator) = isequal(A.dims,B.dims)&&(A.data==B.data)
isapprox(A::QuObject,B::QuObject;kwargs...) = isequal(A.dims,B.dims)&&isapprox(A.data,B.data;kwargs...)

# Conversion and promotion rules
convert(::Type{Ket{T,D}},x::Ket{S,D}) where {T,S,D} = Ket(convert(T,x.data),x.dims)
convert(::Type{Bra{T,D}},x::Bra{S,D}) where {T,S,D} = Bra(convert(T,x.data),x.dims)
convert(::Type{Operator{T,D}},A::Operator{S,D}) where {T,S,D} = Operator(convert(T,A.data),A.dims)
promote_rule(::Type{Ket{T,D}},::Type{Ket{S,D}}) where {T,S,D} = Ket{promote_type(T,S),D}
promote_rule(::Type{Bra{T,D}},::Type{Bra{S,D}}) where {T,S,D} = Bra{promote_type(T,S),D}
promote_rule(::Type{Operator{T,D}},::Type{Operator{S,D}}) where {T,S,D} = Operator{promote_type(T,S),D}

# Misc. conversions
import Compat.LinearAlgebra: Hermitian, Symmetric
Hermitian(A::Operator) = error("Use the `hermitian` function to make an Operator Hermitian")
hermitian(A::Operator) = Operator(Hermitian(A.data),A.dims)
Symmetric(A::Operator) = error("Use the `symmetric` function to make an Operator Symmetric")
symmetric(A::Operator) = Operator(Symmetric(A.data),A.dims)

# Show methods
function show(io::IO, A::T) where T<:QuMatrix
    n = size(A,1)
    print(io, "$n×$n $T with dimensions ", join(A.dims,'⊗'))
end
function show(io::IO, ::MIME"text/plain", A::T) where T<:QuMatrix
    show(io, A)
    println(io); print_array(IOContext(io,:compact=>true), A.data); println(io)
end
show(io::IO, ψ::QuVector) = print(io, braket(ψ))
function show(io::IO, ::MIME"text/plain", ψ::T) where T<:QuVector
    n = size(ψ,1)
    println(io, "$n-d $T with dimensions ", join(ψ.dims,'⊗'))
    println(io, braket(ψ))
end

# Bra-ket printing
function braket(ψ::QuVector, N::Int = 5)
    # N is the max number of kets/bras to show
    idx = findall(x->abs2(x)>2.5E-5,ψ.data) # keep amplitudes larger than 0.005
    isempty(idx) && return "0"
    perm = sortperm(ψ.data[idx], by=abs2, rev=true)
    idx = length(perm)>N ? idx[perm[1:N]] : idx[perm]
    coeffs = map(Array(ψ.data[idx])) do x
        @sprintf("%.2f∠%d°", abs(x), rad2deg(angle(x)))
    end
    labels = join.(tensored_ind2sub.((dims(ψ),), idx.-1), ",")
    return prettybraket(ψ,coeffs,labels) * (length(perm)>N ? " +…" : "")
end

prettybraket(::Ket,coeffs,labels) = join(string.(coeffs,"|",labels,"⟩")," + ")
prettybraket(::Bra,coeffs,labels) = join(string.(coeffs,"⟨",labels,"|")," + ")


text:
"""
Runway fileformat

Example:
#;RUNWAY;2;446;20181206;20190102;241483;EAR_P
06/12/2018;BIAR;01;0000;2359;Y;5;1;1;T
06/12/2018;BIAR;19;0000;2359;N;5;1;1;T
06/12/2018;BIBA;03;0000;2359;N;10;10;5;T
06/12/2018;BIBA;12;0000;2359;Y;10;10;5;T
06/12/2018;BIBA;21;0000;2359;N;10;10;5;T
06/12/2018;BIBA;30;0000;2359;N;10;10;5;T
06/12/2018;BIBD;04;0000;2359;Y;10;10;5;T
06/12/2018;BIBD;22;0000;2359;N;10;10;5;T
06/12/2018;BIBF;18;0000;2359;N;10;10;5;T
06/12/2018;BIBF;36;0000;2359;Y;10;10;5;T
06/12/2018;BIBL;03;0000;2359;Y;10;10;5;T
06/12/2018;BIBL;21;0000;2359;N;10;10;5;T
"""
module Runway

export read

using ..util
using CSV
using Dates
using DataFrames

const ddmmYYYY = DateFormat("dd/mm/YYYY")
const hhmm = DateFormat("HHMM")

const fileformat = Dict(1=>String, 2=>String, 3=>String, 4=>String,
5=>String, 6=>String, 7=>Int64, 8=>Int64, 9=>Int64, 10=>String)

const header_format = ["DATEACTIVE", "AIRPORT", "RWY", "TIMEOPEN",
"TIMECLOSED", "YN", "NUM1", "NUM2", "NUM3", "T"]

function read(file)
    df = CSV.read(file, delim=';', types=fileformat, header=header_format,
    copycols=true, datarow=2)
    reformat!(df)
    return df
end

function reformat!(df)
    df[:,:TEMP] = format_date.(df[:,:DATEACTIVE], ddmmYYYY)
    select!(df, Not(:DATEACTIVE))
    df[:,:DATEACTIVE] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_time.(df[:,:TIMEOPEN], hhmm)
    select!(df, Not(:TIMEOPEN))
    df[:,:TIMEOPEN] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_time.(df[:,:TIMECLOSED], hhmm)
    select!(df, Not(:TIMECLOSED))
    df[:,:TIMECLOSED] =  df[:,:TEMP]
    select!(df, Not(:TEMP))
end

end # module


text:
using GeoData, Test, ArchGDAL
using GeoData: reproject, convertmode

@testset "reproject" begin
    cea = ProjString("+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
    wktcea = convert(WellKnownText, cea)
    projcea = convert(ProjString, cea)
    wkt4326 = convert(WellKnownText, EPSG(4326))
    proj4326 = convert(ProjString, EPSG(4326))

    @test reproject(proj4326, cea, Y(), -30.0) ≈ -3.658789324855012e6
    @test reproject(wktcea, EPSG(4326), Y(), -3.658789324855012e6) ≈ -30.0
    @test reproject(projcea, wkt4326, Y(), -3.658789324855012e6) .≈ -30.0
    @test reproject(cea, proj4326, Y(), [-3.658789324855012e6]) ≈ [-30.0]

    @test reproject(proj4326, cea, X(), 180.0) ≈ 1.7367530445161372e7
    @test reproject(cea, EPSG(4326), X(), 1.7367530445161372e7) ≈ 180.0
    @test reproject(cea, wkt4326, X(), 1.7367530445161372e7) .≈ 180.0
    @test reproject(cea, proj4326, X(), [1.7367530445161372e7]) ≈ [180.0]
end

@testset "convertmode" begin
    projcea = ProjString("+proj=cea")
    proj4326 = convert(ProjString, EPSG(4326))

    lonstart, lonend = 0.5, 179.5
    cealonstart, cealonend = reproject(proj4326, projcea, X(), [lonstart, lonend])
    cealonrange = LinRange(cealonstart, cealonend, 180)
    lon = X(cealonrange; mode=Projected(Ordered(), Regular(step(cealonrange)),
              Intervals(Center()), projcea, proj4326))
    convertedlon = convertmode(Mapped, lon)
    @test all(isapprox.(bounds(convertedlon), (0.0, 180.0); atol=1e-10))
    @test val(convertedlon) ≈ 0.5:179.5

    projectedlon = convertmode(Projected, convertedlon)
    @test val(projectedlon) ≈ val(lon)
    @test all(isapprox.(bounds(projectedlon), bounds(lon); atol=1e-10))

    # Y in not linear
    latbounds = -90.0, 90.0
    ceabounds = reproject(proj4326, projcea, Y(), latbounds)
    # Vals are the bounding range without the end bound
    cealatrange = LinRange(ceabounds..., 181)[1:180]
    lat = Y(cealatrange; mode=Projected(Ordered(), Regular(step(cealatrange)),
              Intervals(Start()), projcea, proj4326))
    convertedlat = convertmode(Mapped, lat)
    @test first(convertedlat) ≈ -90.0 atol=1e-5
    @test all(isapprox.(bounds(convertedlat), (-90.0, 90.0); atol=1e-5))
    @test val(lat) ≈ reproject(proj4326, projcea, Y(), val(convertedlat))

    projectedlat = convertmode(Projected, convertedlat)
    @test val(projectedlat) ≈ val(lat)
    @test all(bounds(projectedlat) .≈ bounds(lat))

    A = DimArray(zeros(length(lon), length(lat)), (lon, lat))
    Aconv = convertmode(Mapped, A)

    @test index(Aconv) == (index(convertedlon), index(convertedlat))
    @test val.(span(Aconv)) == val.(span.((convertedlon, convertedlat)))
end


text:
module PanelGeometry

import Base: +, -, zero

using StaticArrays
using LinearAlgebra
using Rotations
using CoordinateTransformations

using ..AeroMDAO: Point2D, affine_2D, rotation, inverse_rotation, structtolist, sine_spacing, cosine_spacing, partition

## Panel setup
#==========================================================================================#

abstract type AbstractPanel end

include("2d_panels.jl")
include("3d_panels.jl")

# Methods on panels
p1(p :: AbstractPanel) = p.p1
p2(p :: AbstractPanel) = p.p2
p3(p :: AbstractPanel) = p.p3
p4(p :: AbstractPanel) = p.p4

transform_panel(panel_1 :: AbstractPanel, panel_2 :: AbstractPanel) = transform_panel(panel_1, collocation_point(panel_2))

panel_dist(panel_1 :: AbstractPanel, panel_2 :: AbstractPanel) = norm(collocation_point(panel_2) - collocation_point(panel_1))

# split_panels(panels :: Vector{<: Panel}) = collect.(span(panel -> panel_location(panel) == "upper", panels))

end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
# define new aggregator algorithms here and specify their properties

# Direct type methods

"""
Gillespie, Daniel T. (1976). A General Method for Numerically Simulating the
Stochastic Time Evolution of Coupled Chemical Reactions. Journal of
Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
"""
struct Direct <: AbstractAggregatorAlgorithm end

"""
Gillespie, Daniel T. (1976). A General Method for Numerically Simulating the
Stochastic Time Evolution of Coupled Chemical Reactions. Journal of
Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
"""
struct DirectFW <: AbstractAggregatorAlgorithm end

"""
- A. Slepoy, A.P. Thompson and S.J. Plimpton, A constant-time kinetic Monte
  Carlo algorithm for simulation of large biochemical reaction networks, Journal
  of Chemical Physics, 128 (20), 205101 (2008). doi:10.1063/1.2919546

- S. Mauch and M. Stalzer, Efficient formulations for exact stochastic
  simulation of chemical systems, ACM Transactions on Computational Biology and
  Bioinformatics, 8 (1), 27-35 (2010). doi:10.1109/TCBB.2009.47
"""
struct DirectCR <: AbstractAggregatorAlgorithm end

"""
J. M. McCollum, G. D. Peterson, C. D. Cox, M. L. Simpson and N. F. Samatova, The
  sorting direct method for stochastic simulation of biochemical systems with
  varying reaction execution behavior, Computational Biology and Chemistry,
  30 (1), 39049 (2006). doi:10.1016/j.compbiolchem.2005.10.007
"""
struct SortingDirect <: AbstractAggregatorAlgorithm end

"""
- V. H. Thanh, C. Priami and R. Zunino, Efficient rejection-based simulation of
  biochemical reactions with stochastic noise and delays, Journal of Chemical
  Physics, 141 (13), 134116 (2014). doi:10.1063/1.4896985

- V. H. Thanh, R. Zunino and C. Priami, On the rejection-based algorithm for
  simulation and analysis of large-scale reaction networks, Journal of Chemical
  Physics, 142 (24), 244106 (2015). doi:10.1063/1.4922923
"""
struct RSSA <: AbstractAggregatorAlgorithm end

"""
V. H. Thanh, R. Zunino, and C. Priami, Efficient Constant-Time Complexity
Algorithm for Stochastic Simulation of Large Reaction Networks, IEEE/ACM
Transactions on Computational Biology and Bioinformatics, Vol. 14, No. 3,
657-667 (2017).
"""
struct RSSACR <: AbstractAggregatorAlgorithm end

struct RDirect <: AbstractAggregatorAlgorithm end

# NRM-based methods

"""
Gillespie, Daniel T. (1976). A General Method for Numerically Simulating the
Stochastic Time Evolution of Coupled Chemical Reactions. Journal of
Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
"""
struct FRM <: AbstractAggregatorAlgorithm end

"""
Gillespie, Daniel T. (1976). A General Method for Numerically Simulating the
Stochastic Time Evolution of Coupled Chemical Reactions. Journal of
Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
"""
struct FRMFW <: AbstractAggregatorAlgorithm end

"""
M. A. Gibson and J. Bruck, Efficient exact stochastic simulation of chemical
systems with many species and many channels, Journal of Physical Chemistry A,
104 (9), 1876-1889 (2000). doi:10.1021/jp993732q
"""
struct NRM <: AbstractAggregatorAlgorithm end

"""
Elf, Johan and Ehrenberg, M, Spontaneous separation of bi-stable biochemical
systems into spatial domains of opposite phases,Systems Biology, 2004 vol. 1(2)
pp. 230-236. doi:10.1049/sb:20045021
"""
struct NSM <: AbstractAggregatorAlgorithm end

"""
Constant-complexity stochastic simulation algorithm with optimal binning,
Kevin R. Sanft and Hans G. Othmer, Journal of Chemical Physics 143, 074108 (2015);
doi: 10.1063/1.4928635
"""
struct DirectCRDirect <: AbstractAggregatorAlgorithm end

const JUMP_AGGREGATORS = (Direct(),DirectFW(),DirectCR(),SortingDirect(),RSSA(),FRM(),FRMFW(),NRM(),RSSACR(),RDirect())

# For JumpProblem construction without an aggregator
struct NullAggregator <: AbstractAggregatorAlgorithm end

# true if aggregator requires a jump dependency graph
needs_depgraph(aggregator::AbstractAggregatorAlgorithm) = false
needs_depgraph(aggregator::DirectCR) = true
needs_depgraph(aggregator::SortingDirect) = true
needs_depgraph(aggregator::NRM) = true
needs_depgraph(aggregator::RDirect) = true

# true if aggregator requires a map from solution variable to dependent jumps.
# It is implicitly assumed these aggregators also require the reverse map, from
# jumps to variables they update.
needs_vartojumps_map(aggregator::AbstractAggregatorAlgorithm) = false
needs_vartojumps_map(aggregator::RSSA) = true
needs_vartojumps_map(aggregator::RSSACR) = true

is_spatial(aggregator::AbstractAggregatorAlgorithm) = false
is_spatial(aggregator::NSM) = true
is_spatial(aggregator::DirectCRDirect) = true


text:
## MODEL TRAITS

# model trait names:
const MODEL_TRAITS = [:input_scitype, :output_scitype, :target_scitype,
                      :is_pure_julia, :package_name, :package_license,
                      :load_path, :package_uuid,
                      :package_url, :is_wrapper, :supports_weights,
                      :supports_online, :docstring,
                      :name, :is_supervised, :prediction_type,
                      :implemented_methods, :hyperparameters,
                      :hyperparameter_types]
const SUPERVISED_TRAITS = filter(MODEL_TRAITS) do trait
    !(trait in [:output_scitype,])
end
const UNSUPERVISED_TRAITS = filter(MODEL_TRAITS) do trait
    !(trait in [:target_scitype, :prediction_type, :supports_weights])
end

# fallback trait declarations:
input_scitype(::Type) = Unknown
output_scitype(::Type) = Unknown
target_scitype(::Type) = Unknown  # used for measures too
is_pure_julia(::Type) = missing
package_name(::Type) = "unknown"
package_license(::Type) = "unknown"
load_path(::Type) = "unknown"
package_uuid(::Type) = "unknown"
package_url(::Type) = "unknown"
is_wrapper(::Type) = false
supports_online(::Type) = false
supports_weights(::Type) = false  # used for measures too
docstring(M::Type) = string(M)
docstring(M::Type{<:MLJType}) = name(M)
docstring(M::Type{<:Model}) =
    "$(name(M)) from $(package_name(M)).jl.\n"*
    "[Documentation]($(package_url(M)))."

# "derived" traits:
name(M::Type) = string(M)
name(M::Type{<:MLJType}) = split(string(coretype(M)), '.')[end] |> String
is_supervised(::Type{<:Type}) = false
is_supervised(::Type{<:Supervised}) = true
prediction_type(::Type) = :unknown # used for measures too
prediction_type(::Type{<:Deterministic}) = :deterministic
prediction_type(::Type{<:Probabilistic}) = :probabilistic
prediction_type(::Type{<:Interval}) = :interval
implemented_methods(M::Type{<:MLJType}) = map(f->f.name, methodswith(M))
hyperparameters(M::Type) = collect(fieldnames(M))
_fieldtypes(M) = [fieldtype(M, fld) for fld in fieldnames(M)]
hyperparameter_types(M::Type) = string.(_fieldtypes(M))
# function hyperparmeter_defaults(M::Type)
#     try
#         model = M()
#         return [Meta.parse(string(getproperty(model, fld)))
#                 for fld in fieldnames(M)]
#     catch
#         return []
#     end
# end

# following 5 lines commented out because they dissallow precompilation:
# for trait in MODEL_TRAITS
#     eval(quote
#         $trait(object) = $trait(typeof(object))
#     end)
# end

input_scitype(object) = input_scitype(typeof(object))
output_scitype(object) = output_scitype(typeof(object))
target_scitype(object) = target_scitype(typeof(object))
is_pure_julia(object) = is_pure_julia(typeof(object))
package_name(object) = package_name(typeof(object))
package_license(object) = package_license(typeof(object))
load_path(object) = load_path(typeof(object))
package_uuid(object) = package_uuid(typeof(object))
package_url(object) = package_url(typeof(object))
is_wrapper(object) = is_wrapper(typeof(object))
supports_weights(object) = supports_weights(typeof(object))
supports_online(object) = supports_online(typeof(object))
docstring(object) = docstring(typeof(object))
name(object) = name(typeof(object))
is_supervised(object) = is_supervised(typeof(object))
prediction_type(object) = prediction_type(typeof(object))
implemented_methods(object) = implemented_methods(typeof(object))
hyperparameters(object) = hyperparameters(typeof(object))
hyperparameter_types(object) = hyperparameter_types(typeof(object))


text:
"""
 Centermanifold computation and simplest normal form of Hopf bifurcation.
 This code is available to compute centermanifold and normal form of only semisimple Jacobian case (i.e. diagonalizable Jacobian).
 For nonsemisimple case, algorithm should be modified.

 SymEngine package is used.

 For theoretical background of Centermanifold computation see:
    http://www.orcca.on.ca/TechReports/TechReports/2000/TR-00-15.pdf
 for simplest normal form of Hopf bifurcation see:
    https://iopscience.iop.org/article/10.1088/0951-7715/16/1/317
"""
module HopfNormalForm

using SymEngine, OffsetArrays, StaticArrays

export calculate_normal_form, verify, Systems, Flutter_opt

################### Basic functions for symbolic computation in Julia ###################
function coeff2(f::Basic, e::Array{Basic,1}, N)::Basic #coefficient of multivariate Polynomial for univariate case just use diff
    # f is function of plinomials of variable e and N is index vector of power series. for example,
    # if f=e(1)^j*e(2)^k -> I=[j,k]
    v = f
    for i in eachindex(N)
        if N[i] == 0
            v *= 1
        else
            v = diff(v, e[i], N[i]) / factorial(N[i])
        end
    end
    for j in eachindex(N)
        v = subs(v, e[j] => 0)
    end
    return v
end

function recursive_apply(f, x::Basic)::Basic
    cls = SymEngine.get_symengine_class(x)
    args = [recursive_apply(f, arg) for arg in get_args(x)]
    if cls == :Add
        return sum(args)
    elseif cls == :Mul
        return prod(args)
    elseif cls == :Pow
        return args[1]^args[2]
    elseif cls in [:RealDouble, :RealMPFR, :ComplexDouble, :ComplexMPFR]
        return f(x)
    else
        return x
    end
end

function chop_expr(x::Basic, zero_tol::Real)::Basic #chopping very small coefficients to simplify the equaion expression
    return recursive_apply(x -> (abs(N(x)) <= zero_tol) ? zero(typeof(x)) : x, x)
end

function recursive_real(x::Basic)::Basic # Real part of complex symbols. All symbol variables are assumed to be real.
    # To define a complex symbol variable c, express it as c=a+b*I.
    return recursive_apply(real, x)
end

function recursive_imag(x::Basic)::Basic # Imaginary part of complex symbols. All symbol variables are assumed to be real.
    # To define a complex symbol variable c, express it as c=a+b*I.
    return recursive_apply(imag, x)
end

function verify(
    f::Array{Basic},
    T::Array{Basic},
    v::Array{Basic},
    norder::Integer,
    cm::Integer,
    m3::Integer,
    n::Integer,
)::Array{Basic}
    # verify the result of Centermanifold computation upto norder
    # f: function input, T: computed centermanifold, v: state-space variables, norder: oreder of computation,
    # cm: dimension of centermanifold, m3: dimension of stable-unstable manifold
    G2 = Array{Basic}(undef, m3)
    dhf = Array{Basic}(undef, m3)
    dh = Array{Basic}(undef, m3)
    F = Array{Basic}(undef, cm)
    veri = Array{Basic}(undef, m3)
    veri_o = Array{Basic}(undef, m3, norder - 1)
    @vars ϵ
    # Compute F
    for i = 1:cm
        F[i] = f[i] - Eign[i] * v[i]
        for m = 1:n
            if m > cm
                F[i] = subs(F[i], v[m] => T[m-cm])
                F[i] = expand(F[i])
                F[i] = chop_expr(F[i], 1e-10)
            end
        end
    end
    # Compute dh*f
    for i = 1:m3
        dhf[i] = 0
        for l = 1:cm
            dhf[i] += diff(T[i], v[l]) * F[l]
        end
    end
    # Compute g-dh*f
    for k = 1:m3
        G[k] = f[cm+k] - Eign[cm+k] * v[cm+k]
        for m = 1:n
            if m > cm
                G[k] = subs(G[k], v[m] => T[m-cm])
            end
        end
        G2[k] = G[k] - dhf[k]
    end
    # compute dh(x)J_cx-J_s*h(x)
    for i = 1:m3
        dh[i] = 0
        for l = 1:cm
            dh[i] += diff(T[i], v[l]) * Eign[l] * v[l]
        end
        dh[i] = dh[i] - Eign[i+cm] * T[i]
    end
    # veri[i]= g-dh*f - {dh(x)J_cx-J_s*h(x)}
    # We have to check order j=2,3,…,norder of veri[i] is zero
    for i = 1:m3
        veri[i] = dh[i] - G2[i]
        veri[i] = expand(veri[i])
        veri[i] = chop_expr(veri[i], 1e-10)
    end
    # returning value veri_o[i,j] is j-th order of veri[i]
    for j = 2:norder
        for i = 1:m3
            for m = 1:cm
                veri[i] = subs(veri[i], v[m] => ϵ * v[m])
            end
            veri_o[i, j-1] = subs(diff(veri[i], ϵ, j), ϵ => 0) / factorial(j)
        end
    end
    return veri_o
end

#==
Note that complex conjugate pair of eigenvalue corresponding to center subspace is transfromed to +-i.
This is equivalent to rescailing time.

################### Define the index for recursive computation ###################
To compute centermanifold of higher dimension n>3, you should modify this part.
index of coefficient arrays are [k,l,m,n] for centermanifold of dimension 3
k is for k-th part of function g=[g_1,g_2,…,g_k,…,g_m3], l,m,n is for coefficient of polynomial of v_1^l* v_2^m* v_3^n
of power series of g_k
==#

function calculate_normal_form(
    f::Vector{Basic},
    v::Vector{Basic},
    norder::Integer,
    cm::Integer,
    m1::Integer,
    m2::Integer,
    m3::Integer,
    n::Integer,
    ω_h::AbstractFloat,
    zero_tol::AbstractFloat
)
    G = Array{Basic}(undef, m3);
    T = Array{Basic}(undef, m3);
    d = Array{Basic}(undef, m3);
    F = Array{Basic}(undef, cm);

    a = Array{Basic}(undef, cm)
    RD = Array{Basic}(undef, cm, norder); # Reduced dynamics on the center manifold by order
    RD_sum = Array{Basic}(undef, cm); # Reduced dynamics on the center manifold summation of all orders

    @vars ϵ, ν
    # ϵ is symbolic variable which is used in the algorithm which does not have any physical meaning

    if cm == 2 # set multidimensional array for saving coefficints of the polynomial
    # for cm>=3 you should modify this part
        B = OffsetArray{Basic}(undef, 1:m3, 0:norder, 0:norder)
        H = OffsetArray{Basic}(undef, 1:m3, 0:norder, 0:norder)
    elseif cm == 3
        B = OffsetArray{Basic}(undef, 1:m3, 0:norder, 0:norder, 0:norder)
        H = OffsetArray{Basic}(undef, 1:m3, 0:norder, 0:norder, 0:norder)
        A = OffsetArray{Basic}(undef, 1:2, 0:norder, 0:norder, 0:norder)
    end

    Eign = Array{Basic}(undef, n); # Array to save eigenvalues

    ########## Collect Eigenvalues ############
    # Eigenvalues are collected from the input function
    for i = 1:n
        Eign[i] = diff(f[i], v[i])
        for k = 1:n
            Eign[i] = subs(Eign[i], v[k] => 0)
        end
    end

    ########## Initializing the centermanifold and function d (i.e. Dh(x)f(x,h(x))) ######
    for k = 1:m3
        T[k] = 0
        d[k] = 0
    end

    ############ Compute g(x,h(x)) ##############
    for j = 2:norder
        for k = 1:m3
            G[k] = f[cm+k] - Eign[cm+k] * v[cm+k]
            for m = 1:n
                if m > cm
                    G[k] = subs(G[k], v[m] => T[m-cm])
                end
            end
            G[k] = G[k] - d[k]
            for m = 1:cm
                G[k] = subs(G[k], v[m] => ϵ * v[m])
            end
            G[k] = subs(diff(G[k], ϵ, j), ϵ => 0) / factorial(j) # Collect the j-th order of the G[k]
            G[k] = expand(G[k])
            G[k] = chop_expr(G[k], zero_tol)

            if cm == 2
                ind_l = ((k, l, j - l) for l = 0:j)
            elseif cm == 3
                ind_l = ((k, l, m, j - l - m) for l = 0:j for m = 0:j-l)
            end
            for arg in ind_l
                xs = one(Basic)
                δ = zero(Basic)
                B[CartesianIndex(arg)] = coeff2(G[arg[1]], v, arg[2:end])
                for i = 1:cm
                    B[CartesianIndex(arg)] =
                        subs(B[CartesianIndex(arg)], v[i] => 0)
                end
                for i = 1:cm
                    δ += Eign[i] * arg[1+i]
                    xs *= v[i]^arg[1+i]
                end
                δ = δ - Eign[cm+arg[1]]
                H[CartesianIndex(arg)] = B[CartesianIndex(arg)] / δ
                H[CartesianIndex(arg)] = expand(H[CartesianIndex(arg)])
                H[CartesianIndex(arg)] =
                    chop_expr(H[CartesianIndex(arg)], zero_tol)
                T[arg[1]] += H[CartesianIndex(arg)] * xs
                T[arg[1]] = expand(T[arg[1]])
            end
        end
    ############ Compute Dh(x)f(x,h(x)) ##############
        for i = 1:cm
            F[i] = f[i] - Eign[i] * v[i]
            for m = 1:n
                if m > cm
                    F[i] = subs(F[i], v[m] => T[m-cm])
                    F[i] = expand(F[i])
                    F[i] = chop_expr(F[i], zero_tol)
                end
            end
            d[i] = 0
        end
        for i = 1:m3, l = 1:cm
            d[i] += diff(T[i], v[l]) * F[l]
            d[i] = expand(d[i])
            d[i] = chop_expr(d[i], zero_tol)
        end
    end
    ########## Computing reduced dynamics on the centermanifold ##########
    for k = 1:cm, j = 1:norder
        if j == 1
            RD[k, j] = Eign[k] * v[k]
        else
            RD[k, j] = f[k] - Eign[k] * v[k]
            for m = 1:n
                if m > cm
                    RD[k, j] = subs(RD[k, j], v[m] => T[m-cm])
                end
            end
            for m = 1:cm
                RD[k, j] = subs(RD[k, j], v[m] => ϵ * v[m])
            end
            RD[k, j] = subs(diff(RD[k, j], ϵ, j), ϵ => 0) / factorial(j) # Collect the j-th order of Reduced dynmics
            RD[k, j] = expand(RD[k, j])
            RD[k, j] = chop_expr(RD[k, j], zero_tol)
        end
    end
    for n = 1:cm
        RD_sum[n] = sum(RD[n, :])
    end
    ########## Simplest Normal form computation of Hopf bifurcation ##########

    A[1, 1, 0, 1] = recursive_real(coeff2(RD_sum[2], v, [1, 1, 0]));
    A[1, 2, 1, 0] = recursive_real(coeff2(RD_sum[2], v, [0, 2, 1]));
    A[1, 2, 0, 0] = recursive_real(coeff2(RD_sum[2], v, [0, 2, 0]));
    A[2, 1, 1, 0] = recursive_imag(coeff2(RD_sum[2], v, [0, 1, 1]));
    A[2, 2, 0, 0] = recursive_imag(coeff2(RD_sum[2], v, [0, 2, 0]));
    A[1, 1, 1, 0] = recursive_real(coeff2(RD_sum[2], v, [0, 1, 1]));

    S = A[1, 2, 1, 0] - A[1, 2, 0, 0] * A[2, 1, 1, 0] - A[2, 2, 0, 0] * A[1, 1, 1, 0];
    S = expand(S);

    return (CM=T,RD=RD_sum,SLC=S,LCO_R2=-A[1, 1, 0, 1] / S* ν)
end

# Include example systems
include("Systems.jl")
include("Flutter_opt.jl")

end # module


text:
# ---
# title: 1561. Maximum Number of Coins You Can Get
# id: problem1561
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Sort
# link: <https://leetcode.com/problems/maximum-number-of-coins-you-can-get/description/>
# hidden: true
# ---
#
# There are 3n piles of coins of varying size, you and your friends will take
# piles of coins as follows:
#
#   * In each step, you will choose **any  **3 piles of coins (not necessarily consecutive).
#   * Of your choice, Alice will pick the pile with the maximum number of coins.
#   * You will pick the next pile with maximum number of coins.
#   * Your friend Bob will pick the last pile.
#   * Repeat until there are no more piles of coins.
#
# Given an array of integers `piles` where `piles[i]` is the number of coins in
# the `ith` pile.
#
# Return the maximum number of coins which you can have.
#
#
#
# **Example 1:**
#
#
#
#     Input: piles = [2,4,1,2,7,8]
#     Output: 9
#     Explanation: Choose the triplet (2, 7, 8), Alice Pick the pile with 8 coins, you the pile with **7** coins and Bob the last one.
#     Choose the triplet (1, 2, 4), Alice Pick the pile with 4 coins, you the pile with **2** coins and Bob the last one.
#     The maximum number of coins which you can have are: 7 + 2 = 9.
#     On the other hand if we choose this arrangement (1, **2** , 8), (2, **4** , 7) you only get 2 + 4 = 6 coins which is not optimal.
#
#
# **Example 2:**
#
#
#
#     Input: piles = [2,4,5]
#     Output: 4
#
#
# **Example 3:**
#
#
#
#     Input: piles = [9,8,7,6,5,1,2,3,4]
#     Output: 18
#
#
#
#
# **Constraints:**
#
#   * `3 <= piles.length <= 10^5`
#   * `piles.length % 3 == 0`
#   * `1 <= piles[i] <= 10^4`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Documenter
using MixedModelsMakie

makedocs(;
    sitename = "MixedModelsMakie",
    doctest = true,
    checkdocs = :exports,
    pages = [
        "index.md",
        "api.md",
    ],
)

deploydocs(;repo = "github.com/palday/MixedModelsMakie.jl.git", devbranch = "main", push_preview = true)


text:
# Coefficients are taken from https://users.ece.cmu.edu/~koopman/lfsr/index.html
#
# Index by the number of bits needed. Using fewer than 4 bits is apparently not supported,
# which is fine because usually one would want longer sequences anyways.
#
# Shorter sequences are simply generated using a longer LFSR and iterating multiple times.
const FEEDBACK = (
    0x0,
    0x0,
    0x0,
    0x9,
    0x12,
    0x21,
    0x41,
    0x8E,
    0x108,
    0x204,
    0x402,
    0x829,
    0x100D,
    0x2015,
    0x4001,
    0x8016,
    0x10004,
    0x20013,
    0x40013,
    0x80004,
    0x100002,
    0x200001,
    0x400010,
    0x80000D,
    0x1000004,
    0x2000023,
    0x4000013,
    0x8000004,
    0x10000002,
    0x20000029,
    0x40000004,
    0x80000057,
    0x100000029,
    0x200000073,
    0x400000002,
    0x80000003B,
    0x100000001F,
    0x2000000031,
    0x4000000008,
    0x800000001C,
    0x10000000004,
    0x2000000001F,
    0x4000000002C,
    0x80000000032,
    0x10000000000D,
    0x200000000097,
    0x400000000010,
    0x80000000005B,
    0x1000000000038,
    0x200000000000E,
    0x4000000000025,
    0x8000000000004,
    0x10000000000023,
    0x2000000000003E,
    0x40000000000023,
    0x8000000000004A,
    0x100000000000016,
    0x200000000000031,
    0x40000000000003D,
    0x800000000000001,
    0x1000000000000013,
    0x2000000000000034,
    0x4000000000000001,
    0x800000000000000D,
)


text:
ourpath(rows, cols, withna) = joinpath(@__DIR__, "..", "data", string("rows_", rows, "_cols_", cols, "_na_", withna))

ns = [100, 10_000, 1_000_000];
cs = [20, 200];
samples = 5;
missing_share = 0.5;

max_stars_repo_name: davidanthoff/csv-comparison

max_issues_repo_name: davidanthoff/csv-comparison

max_forks_repo_name: davidanthoff/csv-comparison


text:
# streamplots

#makie
using CairoMakie
let
    nonStablePoincare(x, y) = Point2f(x * (x^2 + y^2 - 1) - y * (x^2 + y^2 + 1), y * (x^2 + y^2 - 1) + x * (x^2 + y^2 + 1))
    stableVanDerPaul(x, y) = Point2f(y, (1 - x^2) * y - x)
    semiStable(x, y) = Point2f(-y + x * (-1 + x^2 + y^2)^2, x + y * (-1 + x^2 + y^2)^2)
    titles = ["non-stable", "stable", "semi-stable"]
    functions = [nonStablePoincare, stableVanDerPaul, semiStable]
    cmaps = [:plasma, :viridis, :gnuplot2]

    fig = Figure(resolution = (1000, 400), fontsize = 22)
    axs = [Axis(fig[1, i], xlabel = "x", ylabel = "y", title = titles[i],
        aspect = 1, backgroundcolor = :black) for i in 1:3]
    [streamplot!(axs[i], functions[i], -4 .. 4, -4 .. 4, colormap = Reverse(cmaps[i]),
        gridsize = (32, 32), arrow_size = 10) for i in 1:3]
    [hideydecorations!(axs[2], grid = false, ticks = false) for i in 2:3]
    [hideydecorations!(axs[3], grid = false, ticks = false) for i in 3:3]
    [limits!(axs[i], -4, 4, -4, 4) for i in 1:3]
    display(fig)
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
"""
  Description

    Methods implemented by components that can supply the values of Universal Constants.
  Units

    avogadroConstant => 1/mol
    boltzmannConstant => J/K
    idealGasStateReferencePressure => Pa
    molarGasConstant => J/mol/K
    speedOfLightInVacuum => m/s
    standardAccelerationOfGravity => m/s2
  Methods

    getuniversalconstant()
    getuniversalconstantlist()
"""
module ICapeThermoUniversalConstants
    import ..ECapeExceptions: ECapeInvalidArgument, ECapeUnknown
    export getuniversalconstant,getuniversalconstantlist
    universalConstants=Dict{String,Float64}("avogadroConstant"=>6.0221419947E23,"boltzmannConstant"=>1.380650324E-23,"idealGasStateReferencePressure"=>101325,"molarGasConstant"=>8.31447215,"speedOfLightInVacuum"=>2.997924581E8,"standardAccelerationOfGravity"=>9.80665)

    """
      Retrieves the value of a Universal Constant.

        [in] constantId::String
        [out] constantValue::Float64
    """
    function getuniversalconstant(
      constantId::String)
      try
        return constantValue=universalConstants[constantId]
      catch err
        if(isa(err,KeyError))
          throw(ECapeInvalidArgument())
        else
          throw(ECapeUnknown())
        end
      end
    end

    """
      Returns the identifiers of the supported Universal Constants.

        [out] constantId::Vector{String}
    """
    function getuniversalconstantlist()
      [key::String for key in keys(universalConstants)]
    end

end


text:
#=
architecture.jl
- Julia version 1.0
- Author anderscs
- Date 2018-11-26
=#
using LinearAlgebra, Distributed, Printf, NPZ

function group_argmax(array::Array, n::Int)
    queue = zeros(n)
    max_group = zeros(n)
    for i in 1:length(array)
        a = array[i]
        qmin = minimum(queue)
        if a > qmin
            amin = argmin(queue)
            queue[amin] = a
            max_group[amin] = i
        end
    end
    return max_group
end

function greedy(constraint::Function, indexes::Array{Int}, m_l::Int, parallel::Bool=false, avalible_cores::Int=4)
    """
    Greedy selection of nodes
    """

    selected = Array{Int}([])
    plot = false
    choices = copy(indexes)

    for i in 1:m_l
        @printf("i = %d\n", i)

        calc(node::Int) = constraint(union(selected, node))

        @time values = pmap(calc, choices)

        greedy_choice = choices[argmax(values)]

        if plot
            values = sort(values)
            oplt.plot(values)
            oplt.show()
            # current_best = np.max(values)
        end

        selected = union(selected, [greedy_choice])
        choices = setdiff(choices, [greedy_choice])
        #@printf("selected = %s; choice = %s;\n", selected, greedy_choice)
    end

    return selected
end

function compute_index_set(layer, m_l, weights)
    data = load_data(model_name = "mnist_vgg2",
    prefix = "",
    group_name = "inter_layer_covariance")

    cov_list = data["cov"]
    cov = cov_list[layer + 1, :, :]
    # println("cov for $layer", cov[1:10,1:10])

    shape = size(cov)
    indexes = Array{Int}(1:shape[1])

    theta = 0.5
    W = weights
    #R_z = W' * inv(W * W') * W

    function C(j::Array{Int})
        f = setdiff(indexes, j)
        n = length(f)
        sig_inv = inv(cov[j, j])

        Id = Matrix{Float64}(I, n, n)
        R_z = Matrix{Float64}(I, n, n) # Projection matrix

        ch = theta * Id + (1 - theta) * R_z

        difference = tr(ch * cov[f, j] * sig_inv * cov[j, f])
        normalizer = tr(ch * cov[f, f])
        return difference / normalizer
    end

    j = greedy(C, indexes, m_l)

    return j
end

const results_dir = "results"

function load_data(group_name::String, model_name::String = "mnist_vgg2", prefix::String = "")

    data = Dict()
    dir = joinpath(results_dir, model_name, prefix, group_name)
    for (root, dirs, files) in walkdir(dir)
        for file in files
            parts = split(file, ".")
            println("loading $file in $root")
            if parts[2] == "npy"
                data[parts[1]] = npzread(joinpath(root, file))
            end
        end
    end

    return data
end

function save_data(data::Dict, group_name::String, model_name = "mnist_vgg2", prefix = "")
    dir = joinpath(results_dir, model_name, prefix, group_name)
    for (name, value) in data
        npzwrite(joinpath(root, name), value)
    end

    return data
end

function test()
    cov = Matrix{Float64}(I, 500, 500)
    m_l = 200
    weights = Matrix{Float64}(I, 500, 500)
    compute_index_set(cov, m_l, weights)
end

max_stars_repo_name: Anderssorby/odin

max_issues_repo_name: Anderssorby/odin

max_forks_repo_name: Anderssorby/odin


text:
const YP_y0(::Type{Float64}) = (
    -1.84950800436986690637E16, 4.42733268572569800351E16, -3.46628303384729719441E15,
    8.75906394395366999549E13, -9.82136065717911466409E11, 5.43526477051876500413E9,
    -1.46639295903971606143E7, 1.55924367855235737965E4
)
const YQ_y0(::Type{Float64}) = (
    2.50596256172653059228E17, 3.17157752842975028269E15, 2.02979612750105546709E13,
    8.64002487103935000337E10, 2.68919633393814121987E8, 6.26107330137134956842E5,
    1.04128353664259848412E3, 1.00000000000000000000E0
)
const PP_y0(::Type{Float64}) = (
    9.99999999999999997821E-1, 5.30324038235394892183E0, 8.74716500199817011941E0,
    5.44725003058768775090E0, 1.23953371646414299388E0, 8.28352392107440799803E-2,
    7.96936729297347051624E-4
)
const PQ_y0(::Type{Float64}) = (
    1.00000000000000000218E0, 5.30605288235394617618E0, 8.76190883237069594232E0,
    5.47097740330417105182E0, 1.25352743901058953537E0, 8.56288474354474431428E-2,
    9.24408810558863637013E-4
)
const QP_y0(::Type{Float64}) = (
    -6.05014350600728481186E0, -5.14105326766599330220E1, -1.47077505154951170175E2,
    -1.77681167980488050595E2, -9.32060152123768231369E1, -1.95539544257735972385E1,
    -1.28252718670509318512E0, -1.13663838898469149931E-2
)
const QQ_y0(::Type{Float64}) = (
    2.42005740240291393179E2, 2.06209331660327847417E3, 5.93072701187316984827E3,
    7.24046774195652478189E3, 3.88240183605401609683E3, 8.56430025976980587198E2,
    6.43178256118178023184E1, 1.00000000000000000000E0
)
const PP_j1(::Type{Float64}) = (
    1.00000000000000000254E0, 5.21451598682361504063E0,
    8.42404590141772420927E0, 5.11207951146807644818E0,
    1.12719608129684925192E0, 7.31397056940917570436E-2,
    7.62125616208173112003E-4
)
const PQ_j1(::Type{Float64}) = (
    9.99999999999999997461E-1, 5.20982848682361821619E0,
    8.39985554327604159757E0, 5.07386386128601488557E0,
    1.10514232634061696926E0, 6.88455908754495404082E-2,
    5.71323128072548699714E-4
)
const QP_j1(::Type{Float64}) = (
    2.52070205858023719784E1, 2.11688757100572135698E2,
    5.97489612400613639965E2, 7.10856304998926107277E2,
    3.66779609360150777800E2, 7.58238284132545283818E1,
    4.98213872951233449420E0, 5.10862594750176621635E-2
)
const QQ_j1(::Type{Float64}) = (
    3.36093607810698293419E2, 2.82619278517639096600E3,
    7.99704160447350683650E3, 9.56231892404756170795E3,
    4.98641058337653607651E3, 1.05644886038262816351E3,
    7.42373277035675149943E1, 1.00000000000000000000E0
)
const YP_y1(::Type{Float64}) = (
    -7.78877196265950026825E17, 2.02439475713594898196E17,
    -8.12770255501325109621E15, 1.14509511541823727583E14,
    -6.47355876379160291031E11, 1.26320474790178026440E9
)
const YQ_y1(::Type{Float64}) = (
    3.97270608116560655612E18, 6.87141087355300489866E16,
    6.20557727146953693363E14, 3.88231277496238566008E12,
    1.87601316108706159478E10, 7.34811944459721705660E7,
    2.35564092943068577943E5, 5.94301592346128195359E2,
    1.00000000000000000000E0
)
const besseli0_small_coefs(::Type{Float32}) = (
    1.0000000416172086f0, 0.2499995567758189f0, 0.027778554423089f0,
    0.0017355879591086688f0, 6.96204574858338f-5, 1.8959194799355527f-6,
    4.298760429567388f-8, 3.8884953805853935f-10, 1.483799986601676f-11
)
const besseli0_small_coefs(::Type{Float64}) = (
    0.9999999999999998, 0.2500000000000052, 0.027777777777755364,
    0.001736111111149161, 6.94444444107536e-5, 1.9290123635366806e-6,
    3.9367592765038015e-8, 6.151201574092085e-10, 7.593827956729909e-12,
    7.596677643342155e-14, 6.255282299620455e-16, 4.470993793303175e-18,
    2.1859737023077178e-20, 2.0941557335286373e-22
)
const besseli0_med_coefs(::Type{Float32}) = (
    0.3989423094522082f0, 0.049857687814784155f0, 0.028606362029333827f0,
    0.018926572151825607f0, 0.11423960556773063f0
)
const besseli0_med_coefs(::Type{Float64}) = (
    0.3989422804014326, 0.04986778505064754, 0.028050628512954097,
    0.02921968830978531, 0.04466889626137549, 0.10220642174207666,
    -0.9937439085650689, 91.25330271974727, -4901.408890977662,
    199209.2752981982, -6.181516298413396e6, 1.4830278710991925e8,
    -2.7695254643719645e9, 4.0351394830842026e10, -4.5768930327229974e11,
    4.0134844243070063e12, -2.6862476523182016e13, 1.3437999451218112e14,
    -4.856333741437621e14, 1.1962791200680235e15, -1.796269414464399e15,
    1.239942074380968e15
)
const besseli0_large_coefs(::Type{Float32}) =  (
   0.398942288181429f0, 0.04986085645388595f0, 0.028961207830173957f0
)
const besseli0_large_coefs(::Type{Float64}) = (
    0.3989422804014327, 0.04986778505003549, 0.028050629664438713,
    0.029218603406797997, 0.045199090067282434
)
const besseli1_small_coefs(::Type{Float32}) = (
    0.08333336088599108f0, 0.0069442679189687774f0, 0.00034740666242346436f0,
    1.1502079282185252f-5, 2.8884784673820464f-7, 3.686308746029452f-9,
    1.2271327751884605f-10
)
const besseli1_med_coefs(::Type{Float32}) = (
    3.98942115977513013f-1, -1.49581264836620262f-1, -4.76475741878486795f-2,
    -2.65157315524784407f-2, -1.47148600683672014f-1
)
const besseli1_small_coefs(::Type{Float64}) = (
    0.08333333333333334, 0.006944444444444374, 0.00034722222222248526,
    1.1574074073690356e-5, 2.7557319253050506e-7, 4.920949730519126e-9,
    6.834656365321179e-11, 7.593985414952446e-13, 6.904652315442046e-15,
    5.2213850252454655e-17, 3.405120412140281e-19, 1.6398527256182257e-21,
    1.3161876924566675e-23
)
const besseli1_med_coefs(::Type{Float64}) = (
    0.39894228040143276, -0.149603355151029, -0.04675104787903509,
    -0.04090746353279043, -0.05744911840910781,-0.12283724006390319,
    1.0023632527650936, -94.90954045770921, 5084.06899084327,
    -206253.5613716743, 6.387439538535799e6, -1.529244018923123e8,
    2.849523551208316e9, -4.141884344471782e10, 4.6860149658304974e11,
    -4.097852944580042e12, 2.7345051110005453e13, -1.3634833112030364e14,
    4.909983186948892e14, -1.2048200837913132e15, 1.8014682382937435e15,
    -1.2377987428989558e15
)
const besseli1_large_coefs(::Type{Float64}) = (
    0.39894228040143265, -0.14960335515036183, -0.0467510491854453,
    -0.04090618769936314, -0.05808395099511361
)
const besseli1_large_coefs(::Type{Float32}) = (
    0.3989422695655856f0, -0.1495936968314777f0, -0.04802187089304704f0
)
const JP_j0(::Type{Float32}) = (
    -1.729150680240724f-1, 1.332913422519003f-2, -3.969646342510940f-4,
    6.388945720783375f-6, -6.068350350393235f-8
)
const MO_j0(::Type{Float32}) = (
    7.978845717621440f-1, -3.355424622293709f-6, -4.969382655296620f-2,
    -3.560281861530129f-3, 1.197549369473540f-1, -2.145007480346739f-1,
    1.864949361379502f-1, -6.838999669318810f-2
)
const PH_j0(::Type{Float32}) = (
    -1.249992184872738f-1,  6.490598792654666f-2, -1.939906941791308f-1,
    1.001973420681837f0, -4.974978466280903f0, 1.756221482109099f1,
    -3.630592630518434f1, 3.242077816988247f1,
)
const JP32 = (
    -3.405537384615824f-2, 1.937383947804541f-3, -4.541343896997497f-5,
    6.009061827883699f-7, -4.878788132172128f-9
)
const YP32 = (
    4.202369946500099f-2, -2.641785726447862f-3, 6.719543806674249f-5,
    -9.496460629917016f-7, 8.061978323326852f-9
)
const MO132 = (
    7.978845453073848f-1, 4.976029650847191f-6, 1.493389585089498f-1,
    5.435364690523026f-3, -2.102302420403875f-1, 3.138238455499697f-1,
    -2.284801500053359f-1,6.913942741265801f-2,
)
const PH132 = (
    3.749989509080821f-1, -1.637986776941202f-1, 3.503787691653334f-1,
    -1.544842782180211f0, 7.222973196770240f0, -2.485774108720340f1,
    5.073465654089319f1, -4.497014141919556f1,
)
const YP_y1(::Type{Float32}) = (
    1.707584643733568f-1, -1.584289289821316f-2, 5.344486707214273f-4,
    -9.413212653797057f-6, 9.454583683980369f-8
)
const MO_y1(::Type{Float32}) = (
    7.978845717621440f-1, -3.355424622293709f-6, -4.969382655296620f-2,
    -3.560281861530129f-3, 1.197549369473540f-1, -2.145007480346739f-1,
    1.864949361379502f-1, -6.838999669318810f-2
)
const PH_y1(::Type{Float32}) = (
    -1.249992184872738f-1,  6.490598792654666f-2, -1.939906941791308f-1,
    1.001973420681837f0, -4.974978466280903f0, 1.756221482109099f1,
    -3.630592630518434f1, 3.242077816988247f1,
)
const P1_k0(::Type{Float32}) = (
    -1.372508979104259711f-1, 2.622545986273687617f-1, 5.047103728247919836f-3
)
const Q1_k0(::Type{Float32}) = (
    1.000000000000000000f0, -8.928694018000029415f-2, 2.985980684180969241f-3
)
const P2_k0(::Type{Float32}) = (
    1.159315158f-1, 2.789828686f-1, 2.524902861f-2,
    8.457241514f-4, 1.530051997f-5
)
const P3_k0(::Type{Float32}) = (
    2.533141220f-1, 5.221502603f-1,
    6.380180669f-2, -5.934976547f-2
)
const Q3_k0(::Type{Float32}) = (
    1.000000000f0, 2.679722431f0,
    1.561635813f0, 1.573660661f-1
)
const P1_k0(::Type{Float64}) = (
    -1.372509002685546267e-1, 2.574916117833312855e-1,
    1.395474602146869316e-2, 5.445476986653926759e-4,
    7.125159422136622118e-6
)
const Q1_k0(::Type{Float64}) = (
    1.000000000000000000e+00, -5.458333438017788530e-02,
    1.291052816975251298e-03, -1.367653946978586591e-05
)
const P2_k0(::Type{Float64}) = (
    1.159315156584124484e-01, 2.789828789146031732e-01,
    2.524892993216121934e-02, 8.460350907213637784e-04,
    1.491471924309617534e-05, 1.627106892422088488e-07,
    1.208266102392756055e-09, 6.611686391749704310e-12
)
const P3_k0(::Type{Float64}) = (
    2.533141373155002416e-1, 3.628342133984595192e0,
    1.868441889406606057e1, 4.306243981063412784e1,
    4.424116209627428189e1, 1.562095339356220468e1,
    -1.810138978229410898e0, -1.414237994269995877e0,
    -9.369168119754924625e-2
)
const Q3_k0(::Type{Float64}) = (
    1.000000000000000000e0, 1.494194694879908328e1,
    8.265296455388554217e1, 2.162779506621866970e2,
    2.845145155184222157e2, 1.851714491916334995e2,
    5.486540717439723515e1, 6.118075837628957015e0,
    1.586261269326235053e-1
)
const Y_k0 = 1.137250900268554688

const Y_k1(::Type{Float32}) = 8.695471287f-2
const Y_k1(::Type{Float64}) = 8.69547128677368164e-2

const Y2_k1(::Type{Float32}) = 1.450342178f0
const Y2_k1(::Type{Float64}) = 1.45034217834472656

const P1_k1(::Type{Float32}) = (
    -3.621379531f-3, 7.131781976f-03, -1.535278300f-5
)
const P1_k1(::Type{Float64}) = (
    -3.62137953440350228e-3, 7.11842087490330300e-3,
    1.00302560256614306e-5, 1.77231085381040811e-6
)
const Q1_k1(::Type{Float32}) = (
    1.000000000f0, -5.173102701f-2, 9.203530671f-4
)
const Q1_k1(::Type{Float64}) = (
    1.00000000000000000e0, -4.80414794429043831e-2,
    9.85972641934416525e-4, -8.91196859397070326e-6
)
const P2_k1(::Type{Float32}) = (
    -3.079657469f-1, -8.537108913f-2,
    -4.640275408f-3, -1.156442414f-4
)
const P2_k1(::Type{Float64}) = (
    -3.07965757829206184e-1, -7.80929703673074907e-02,
    -2.70619343754051620e-3, -2.49549522229072008e-5
)
const Q2_k1(::Type{Float32}) = (
    one(Float32), zero(Float32)
)
const Q2_k1(::Type{Float64}) = (
    1.00000000000000000e0, -2.36316836412163098e-2,
    2.64524577525962719e-4, -1.49749618004162787e-6
)

const P3_k1(::Type{Float32}) = (
    -1.970280088f-1, 2.188747807f-2,
    7.270394756f-1, 2.490678196f-1
)
const P3_k1(::Type{Float64}) = (
    -1.97028041029226295e-1, -2.32408961548087617e0,
    -7.98269784507699938e0, -2.39968410774221632e0,
    3.28314043780858713e1, 5.67713761158496058e1,
    3.30907788466509823e1, 6.62582288933739787e0,
    3.08851840645286691e-1
)
const Q3_k1(::Type{Float32}) = (
    1.000000000f0, 2.274292882f0,
    9.904984851f-1, 4.585534549f-2
)
const Q3_k1(::Type{Float64}) = (
    1.00000000000000000e0, 1.41811409298826118e1,
    7.35979466317556420e1, 1.77821793937080859e2,
    2.11014501598705982e2, 1.19425262951064454e2,
    2.88448064302447607e1, 2.27912927104139732e0,
    2.50358186953478678e-2
)


text:
@doc """
toFunctionMXelements(as::Vector{Vector{Matrix{T}}}) where T <: Real

Return an array of discretized matrix-time function in the for of an array of matrices of discretised functions (as elements)

## Arguments
- `as::Vector{Vector{Matrix{T}}}`: Vector of discretised matrix function
""" toFunctionMXelements
function toFunctionMXelements(as::AbstractVector{<:AbstractVector{<:AbstractArray{T}}},::ItoIsometryMethod{K}) where {T <: Real,K}
    idxs=CartesianIndices(as[1][1])
    return [[SVector{K}([
        aₜ[idx]
        for aₜ in aᵢ]) for idx in idxs] for aᵢ in as]
end
function toFunctionMXelements(as::AbstractVector{<:AbstractArray{T}},::ItoIsometryMethod{K}) where {T <: Real,K}
    idxs=CartesianIndices(as[1])
    return [SVector{K}([
        aₜ[idx]
        for aₜ in as]) for idx in idxs]
end

###############################################################################
############################### CovMX vs CovVec ###############################
###############################################################################

@doc """
MxToCovVec(mx::AbstractArray,SDmatrixSize::Int64)

Create covariance vector from covariance matrix
- `mx::AbstractArray`: covariance matrix to convert
- `SDmatrixSize::Int64`: size of the large mapping matrix
""" MxToCovVec
function MxToCovVec(mx::AbstractArray, SDmatrixSize::Int64)
#take raw excitation matrix and turn it into the semi discretized 2nd Moment vector
    mxdim = size(mx, 1);
    ri = reverse(0:SDmatrixSize)
    vdim = sum(ri)
    result = spzeros(vdim)
    idxs = idxTuples(mxdim)
    rstrt = 0;
    rstp = mxdim;
    strt = 0;
    stp = mxdim;
    for i in 1:mxdim
        result[rstrt + 1:rstp] = (x -> mx[x...]).(idxs[strt + 1:stp])
        strt = stp;
        stp += mxdim - i;
        rstrt = rstrt + ri[i]
        rstp = rstp + ri[i + 1]
    end
    return result
end

@doc """
    VecToCovMx(CovV::AbstractArray,dims::Int64)

Create covariance matrix from covariance vector
- `CovV::AbstractArray`: covariance vector to convert
- `dims::Int64`: size of the state space
""" VecToCovMx
function VecToCovMx(CovV::AbstractArray, SDmatrixSize::Int64)
  #Covariance vector to Covariance matrix
          dim = SDmatrixSize;  from = 1;  to = SDmatrixSize;  i = 0;
    result = diagm(i => vec(CovV[from:to]));
    while (dim > 0)
                    from += dim; dim -= 1; to = from + dim - 1; i += 1;
        result[:,:] .+= diagm(i => vec(CovV[from:to]));
    end
    return Symmetric(result)
end

function idxTuples(nN::Int64)
    idx = [[(i, i - j) for i in 1 + j:nN] for j in 0:nN - 1]
    idxrng = 1:length(idx)
    result = Array{Tuple{Int64,Int64},1}(undef, sum(idxrng))
    iIDX = prepend!(accumulate(+, reverse(idxrng)), 0)
    for i in 1:length(idx)
        result[iIDX[i] + 1:iIDX[i + 1]] .= idx[i]
    end
    return result
end

###############################################################################
############################ M2 Mapping Utilities #############################
###############################################################################
struct CovVecIdx
    sectionStarts::Vector{Int64} #Every diagonal's start-1 in the covariance vector
    function CovVecIdx(siz::Int64)
        resvec = Vector{Int64}(undef, siz + 1);
        resvec[1] = 0;
        for (i, n) in enumerate(siz:-1:1)
            resvec[i + 1] = resvec[i] + n;
        end
        new(resvec)
    end
    CovVecIdx(MX::AbstractArray) = CovVecIdx(size(MX, 1))
end

function (cvIdx::CovVecIdx)(i::Int64, j::Int64)
    cvIdx.sectionStarts[abs(i - j) + 1] + min(i, j)
end

# Squared process generator #######################################################
#
function calculate_noise_mxelems(rst::AbstractResult{d}) where d
    ns = zeros(Int64,rst.LDDEP.w);
    for stsubMXsᵢ in rst.stsubMXs
        ns[stsubMXsᵢ[1].nID] += d^2 * length(stsubMXsᵢ[1].MXfun)
    end
    ns
end
function calculate_noise_velems(rst::AbstractResult{d}) where d
    ns = zeros(Int64,rst.LDDEP.w);
    for stsubVᵢ in rst.stsubVs
        ns[stsubVᵢ[1].nID] += d
    end
    ns
end

function M2_Mapping_from_Sparse(SM::SparseMatrixCSC{TSMX,Int64},rst::AbstractResult) where TSMX<: Union{<:AbstractVector{T},<:T} where T <:Real
    idx = CovVecIdx(SM)
    (_Is, _Js, _Vs) = findnz(SM)
    elemnum = sum(eachindex(_Vs))
    Is = Vector{Int64}(undef, elemnum)
    Js = Vector{Int64}(undef, elemnum)
    Vs = Vector{T}(undef, elemnum)
    k = 1;
    for i in eachindex(_Is)
        for j in i:length(_Is)
            Is[k] = idx(_Is[i], _Is[j])
            Js[k] = idx(_Js[i], _Js[j])
            if Is[k] <= idx.sectionStarts[2] && Js[k] > idx.sectionStarts[2]
                Vs[k] = 2 * rst.itoisometrymethod(_Vs[i],_Vs[j])
            else
                Vs[k] = rst.itoisometrymethod(_Vs[i],_Vs[j])
            end
            k += 1
        end
    end
    sparse(Is, Js, Vs, idx.sectionStarts[end], idx.sectionStarts[end])
end

function M2_Mapping_from_Sparse(SM::SparseMatrixCSC{TSMX,Int64},idxs::AbstractVector{<:Integer},rst::AbstractResult) where TSMX<: Union{<:AbstractVector{T},<:T} where T <:Real
    M2_Mapping_from_Sparse(SM[idxs,idxs],rst)
end

function M1toM2_Mapping_Generator_from_Sparse(FM::SparseMatrixCSC{TSMX,Int64}, FV::SparseVector{TSV,Int64},rst::AbstractResult) where TSMX<: Union{<:AbstractVector{T},<:T} where TSV<: Union{<:AbstractVector{T},<:T} where T <:Real # Condition: FM and FV have the same dimension (FM ∈ Rⁿˣⁿ, FV ∈ Rⁿ)
    idx = CovVecIdx(FV)
    (_MIs, _MJs, _MVs) = findnz(FM)
    (_FIs, _FVs) = findnz(FV)
    elemnum = length(_MVs) * length(_FVs)
    Is = Vector{Int64}(undef, elemnum)
    Js = Vector{Int64}(undef, elemnum)
    Vs = Vector{T}(undef, elemnum)
    k = 1
    for i in eachindex(_MIs)
        for j in eachindex(_FIs)
            Is[k] = idx(_MIs[i], _FIs[j])
            Js[k] = _MJs[i]
            if _MIs[i] == _FIs[j]
                Vs[k] = 2 * rst.itoisometrymethod(_MVs[i],_FVs[j])
                # Vs[k] = 2 * _MVs[i] * _FVs[j]
            else
                Vs[k] = rst.itoisometrymethod(_MVs[i],_FVs[j])
                # Vs[k] = _MVs[i] * _FVs[j]
            end
            k += 1
        end
    end
    return sparse(Is, Js, Vs, idx.sectionStarts[end], idx.sectionStarts[2])
end

function M1toM2_Mapping_Generator_from_Sparse(FM::SparseMatrixCSC{TSMX,Int64}, FV::SparseVector{TSV,Int64},idxs::AbstractVector{<:Integer},rst::AbstractResult)  where TSMX<: Union{<:AbstractVector{T},<:T} where TSV<: Union{<:AbstractVector{T},<:T} where T <:Real # Condition: FM and FV have the same dimension (FM ∈ Rⁿˣⁿ, FV ∈ Rⁿ)
    M1toM2_Mapping_Generator_from_Sparse(FM[idxs,idxs], FV[idxs], rst)
end

function M2_Additive_from_Sparse(FV::SparseVector{TSMX,Int64},rst::AbstractResult) where TSMX<: Union{<:AbstractVector{T},<:T} where T <:Real
    idx = CovVecIdx(FV)
    (_Is, _Vs) = findnz(FV)
    elemnum = sum(1:length(_Vs))
    Is = Vector{Int64}(undef, elemnum)
    Vs = Vector{T}(undef, elemnum)
    k = 1
    for i in eachindex(_Is)
        for j in i:length(_Is)
            Is[k] = idx(_Is[i], _Is[j])
            Vs[k] = rst.itoisometrymethod(_Vs[i], _Vs[j])
            k += 1
        end
    end
    sparsevec(Is, Vs, idx.sectionStarts[end])
end
function M2_Additive_from_Sparse(FV::SparseVector{TSMX,Int64},idxs::AbstractVector{<:Integer},rst::AbstractResult) where TSMX<: Union{<:AbstractVector{T},<:T} where T <:Real
    M2_Additive_from_Sparse(FV[idxs],rst)
end

# Reduction of the 2nd moment mapping vectors on a time interval
function reduce_additive(M2s::Vector{TM2}, M12s::Vector{TM12}, M1s::Vector{TM1}, V2s::Vector{TV2}, V1s::Vector{TV1}, V1st::TV1st) where {TM2<:AbstractMatrix, TM12<:AbstractMatrix, TM1<:AbstractMatrix, TV2<:AbstractVector, TV1<:AbstractVector, TV1st<:AbstractVector}
    mV1 = V1st;
    mV2 = M12s[1]*mV1 .+ V2s[1];

    # mV2 = V2s[1];
    for i in 2:length(V2s)
        mV1 .= M1s[i] * mV1 .+ V1s[i-1]
        mV2 .= M2s[i] * mV2 .+ M12s[i]*mV1 .+ V2s[i]

        # mV2 .= M2s[i] * mV2 .+ M12s[i]*mV1 .+ V2s[i]
        # mV1 .= M1s[i-1] * mV1 .+ V1s[i-1]
    end
    mV2
end


text:
# Helper functions

function get_all_puzzles(d::Array{Int,2})
    all_puzzles = d[:,Int(RL_Fit.c_puzzle)]
    sort(unique(all_puzzles))
end
function get_all_episodes_in_puzzle(d::Array{Int,2}, i_puzzle::Int)
    idx = d[:,Int(RL_Fit.c_puzzle)].== i_puzzle
    all_epi = d[idx,Int(RL_Fit.c_episode)]
    sort(unique(all_epi))
end
function get_data_by_puzzle_episode(d::Array{Int,2}, i_puzzle::Int, i_episode::Int)
    idx = (d[:,Int(RL_Fit.c_puzzle)].== i_puzzle) .& (d[:,Int(RL_Fit.c_episode)].==i_episode)
    d[idx,:]
end
function get_data_stats_and_dimensions(data)
    data_stats = get_data_stats(data)
    nr_states = data_stats.max_state_id
    nr_actions = data_stats.max_action_id
    (nr_samples,) = size(data)
    data_stats, nr_states, nr_actions, nr_samples
end

converttoInt(variable) = typeof(variable) == Float64 ? Int(round(variable)) : variable


text:
using Test
using Distributions
using Random

using EnsembleKalmanProcesses.DataStorage

@testset "DataStorage" begin
    seed = 2021
    Random.seed!(seed)

    parameter_samples = rand(MvNormal(2, 0.1), 10) #10 samples of 4D params
    data_samples = rand(MvNormal(12, 2), 10) #10 samples of 12D data
    data_samples_short = data_samples[:, 1:(end - 1)]

    new_parameter_samples = rand(MvNormal(2, 0.1), 10) #10 samples of 4D params
    new_data_samples = rand(MvNormal(12, 2), 10) #10 samples of 12D data
    new_data_samples_short = new_data_samples[:, 1:(end - 1)]

    #test DataContainer
    parameter_samples_T = permutedims(parameter_samples, (2, 1))
    idata = DataContainer(parameter_samples)
    idata_T = DataContainer(parameter_samples_T, data_are_columns = false)

    @test get_data(idata) == parameter_samples
    @test get_data(idata_T) == parameter_samples
    @test size(idata) == size(parameter_samples)

    #test Paired Container
    idata = DataContainer(parameter_samples)
    odata = DataContainer(data_samples)
    odata_short = DataContainer(data_samples_short)

    @test_throws DimensionMismatch PairedDataContainer(parameter_samples, data_samples_short)
    iopairs = PairedDataContainer(parameter_samples, data_samples)

    @test_throws DimensionMismatch PairedDataContainer(idata, odata_short)
    iopairs2 = PairedDataContainer(idata, odata)

    @test get_data(iopairs) == get_data(iopairs2)
    @test get_inputs(iopairs) == parameter_samples
    @test get_outputs(iopairs) == data_samples

end


text:
using HSA
using HSA.Execution
using FactCheck

using HSA.Execution: prepare_args, cleanup_args, get_or_init_defaults, clear_defaults,
                     build_kernel, allocate_args, kernel_cache, build_dispatch

include("kernels.jl")

facts("The execution framework") do
    context("manages the configuration details for kernel invocations") do
        # does not hold a runtime reference initially
        @fact_throws HSA.status_string(4107)

        cfg = get_or_init_defaults()

        @fact HSA.status_string(4107) --> not(nothing) # does not throw
        @fact cfg.queue --> not(nothing)
        @fact cfg.agent --> not(nothing)

        # disposes runtime reference
        clear_defaults()

        @fact_throws HSA.status_string(4107)
    end

    context("can automatically convert kernel arguments") do
        rt = Runtime()

        context("Arrays to pointers, registering the Memory with the runtime") do
            args = Array[ Array(Int, 10), Array(Float64, 5) ]

            kargs = prepare_args(args)

            for ka in kargs
                @fact ka --> x -> isa(x, Ptr)
            end

            cleanup_args(args)
        end

        finalize(rt)
    end

    context("can build a kernel and cache the executable") do
        cfg = get_or_init_defaults()
        @fact length(kernel_cache) --> 0

        kernel_info = build_kernel(cfg.agent, vadd, [Ptr{Float64},Ptr{Float64},Ptr{Float64}])

        @fact kernel_info --> not(nothing)
        @fact length(kernel_cache) --> 1

        kernel_info2 = build_kernel(cfg.agent, vadd, [Ptr{Float64}, Ptr{Float64}, Ptr{Float64}])

        @fact kernel_info2 --> kernel_info
        @fact length(kernel_cache) --> 1

        @pending HSA.executable_get_symbol(
            kernel_info.executable,
            "&vadd") --> not(nothing)

    end

    context("can allocate memory for arguments") do
        cfg = get_or_init_defaults()
        kernel_info = build_kernel(cfg.agent, vadd, [Ptr{Float64},Ptr{Float64},Ptr{Float64}])

        args = Any[Ptr{Float64}(0), Ptr{Float64}(0), Ptr{Float64}(0)]

        karg_memory = allocate_args(cfg.agent, kernel_info, args)

        @fact karg_memory.ptr --> not(0)
    end

    context("can build the dispatch packet") do
        cfg = get_or_init_defaults()
        kernel_info = build_kernel(cfg.agent, vadd, [Ptr{Float64},Ptr{Float64},Ptr{Float64}])
        args = Any[Ptr{Float64}(1), Ptr{Float64}(2), Ptr{Float64}(3)]
        karg_memory = allocate_args(cfg.agent, kernel_info, args)

        signal = Signal(value = 1)

        packet = build_dispatch((100,), kernel_info, karg_memory, signal)

        @fact packet --> x -> isa(x, HSA.KernelDispatchPacket)

        finalize(karg_memory)
        finalize(signal)
    end

    context("can execute a vector copy") do
        const N = 1000

        a = Array(Int,N); rand!(a)
        b = Array(Int,N); rand!(b)
        expected = copy(a)

        @hsa (N) vcopy(a,b)

        @fact b --> expected
    end

    context("can execute a float vector copy") do
        const N = 1000

        a = Array(Float64,N); rand!(a)
        b = Array(Float64,N); rand!(b)
        expected = copy(a)

        @hsa (N) vcopy(a,b)

        @fact b --> expected
    end

    context("can execute a vector add on pointers") do
        const n = 100
        args = Array[ Array(Int, n), Array(Int, n), Array(Int, n) ]
        kargs = [ pointer(a) for a in args ]

        for a in args
            rand!(a)
        end

        expected = args[1] + args[2]

        @hsa (n) vadd(kargs[1], kargs[2], kargs[3])

        @fact args[3] --> expected
    end

    context("can execute a kernel calling another function, as long as it's \
    inlined") do
        const n = 100
        args = Array[ Array(Int, n), Array(Int, n), Array(Int, n) ]

        for a in args
            rand!(a)
        end

        expected = args[1] + args[2]

        @hsa (n) vadd_fun(args[1], args[2], args[3])

        @fact args[3] --> expected
    end

    context("can execute a simple matrix multiplication") do
        const arows = 3
        const acols = 4

        a = Array(Float64, arows, acols); rand!(a)
        b = Array(Float64, acols, arows); rand!(b)
        c = Array(Float64, arows, arows); rand!(c)
        c_expected = a * b

        @hsa (arows) mmul(a,b,c,acols)

        @fact c --> c_expected
    end

    context("can execute a 2d matrix multiplication") do
        const arows = 3
        const acols = 4

        a = Array(Float64, arows, acols); rand!(a)
        b = Array(Float64, acols, arows); rand!(b)
        c = Array(Float64, arows, arows); rand!(c)
        c_expected = a * b

        @hsa (arows, arows) mmul2d(a,b,c,acols)

        @fact c --> c_expected
    end
end


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

using Test, Libdl, LibSSH2_jll

@testset "LibSSH2_jll" begin
    # We use a `startswith()` here because when built from source, this returns "1.9.0_DEV"
    vn = startswith(unsafe_string(ccall((:libssh2_version, libssh2), Cstring, (Cint,), 0)), "1.9.0")
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
struct SPSSDataFrame
    vars::VariableDictionary
    data::Vector{Any}
    producer::String
    creation_date::String
    creation_time::String
    label::String
end

function read_sav(io::IO)
    comp, bias, casesz, ncases, producer, cdate, ctime, flabel = readheader(io)
    bias = Int(bias)
    rectyp = read(io, Int32)
    @assert rectyp == 2 "variable dictionary records should immediately follow header"
    rectyp, vdict = readvariabledictionary(io, max(casesz, 5))
    casesz = length(vdict.nms)
    ibuf = Vector{UInt8}(casesz)
    fbuf = Vector{Union{Float64,Missing}}(casesz)
    if rectyp == 3   # read value labels, if any
        rectyp = readvaluelabels!(io, vdict, comp == 1, bias)
    end
    data = [missings(T, Int64(ncases)) for T in vdict.typs]
    if rectyp == 6
        rectyp = readdocumentrecord(io)
    end
    if rectyp == 7
        rectyp, infodict = readextensionrecords(io, vdict)
    end
    if rectyp == 999
        skip(io, 4)
        if comp == 1
            typs = vdict.typs
            @inbounds for k in 1:ncases
                readbytecoderec!(io, ibuf, fbuf, bias)
                for j in 1:length(typs)
                    datj = data[j]
                    ibj = ibuf[j]
                    fbj = fbuf[j]
                    T = typs[j]
                    if T == String
                        if !ismissing(fbj)
                            datj[k] = strip(String(reinterpret(UInt8, [fbj])))
                        end
                    elseif T == Float64
                        datj[k] = fbj
                    elseif T == UInt8
                        datj[k] = ibj == 0xff ? missing : ibj
                    else
                        throw(ArgumentError("Unexpected type $T encountered in typs"))
                    end
                end
            end
        end
    end
    if !isempty(vdict.valuelabels)
        for (k,v) in vdict.valuelabels
            ind = vdict.invnms[k]
            data[ind] = droplevels!(recode(categorical(data[ind], true, ordered = true), v...))
        end
    end
    for (i,v) in enumerate(data)
        if !any(ismissing, v)
            data[i] = disallowmissing(v)
        end
    end
    SPSSDataFrame(vdict, data, producer, cdate, ctime, flabel), infodict
end


text:
using LinearAlgebra
using Test

using SetProg, SetProg.Sets
using Polyhedra
using MultivariatePolynomials
import DynamicPolynomials

using JuMP
const MOIT = MOI.Test

@testset "Switched System" begin
    config = MOIT.Config()
    @testset "Ellipsoid" begin
        @testset "Feasible" begin
            optimize!(mock) = MOIU.mock_optimize!(
                mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, zeros(3)))
            Tests.feasible_switched_system_ell_test(bridged_mock(optimize!), config)
        end
        @testset "Infeasible" begin
            optimize!(mock) = MOIU.mock_optimize!(
                mock, MOI.INFEASIBLE, MOI.NO_SOLUTION, MOI.INFEASIBILITY_CERTIFICATE,
                (MOI.VectorOfVariables, SetProg.SumOfSquares.PositiveSemidefinite2x2ConeTriangle) => [[0.00282108, 0.0, 0.00282108]],
                (MOI.VectorAffineFunction{Float64}, SetProg.SumOfSquares.PositiveSemidefinite2x2ConeTriangle) => [[0.00117187, 0.0, 2.82044], [2.82044, 0.0, 0.00117187]]
            )
            Tests.infeasible_switched_system_ell_test(bridged_mock(optimize!), config)
        end
    end
    @testset "Quadratic" begin
        @testset "Feasible" begin
            optimize!(mock) = MOIU.mock_optimize!(
                mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, zeros(3)))
            Tests.feasible_switched_system_quad_test(bridged_mock(optimize!), config)
        end
        # Blocked by a bug in MockOptimizer which does not implement support for MomentMatrixAttribute correctly
#        @testset "Infeasible" begin
#            function optimize!(mock)
#                MOI.set(mock, MOI.TerminationStatus(), MOI.INFEASIBLE)
#                MOI.set(mock, MOI.PrimalStatus(), MOI.NO_SOLUTION)
#                MOI.set(mock, MOI.DualStatus(), MOI.INFEASIBILITY_CERTIFICATE)
#                DynamicPolynomials.@polyvar x[1:2]
#                for (F, S) in MOI.get(mock, MOI.ListOfConstraints())
#                    if F == MOI.VectorAffineFunction{Float64}
#                        cis = MOI.get(mock, MOI.ListOfConstraintIndices{F, S}())
#                        function _moment_matrix(q)
#                            SetProg.SumOfSquares.build_moment_matrix(q, monovec(x))
#                        end
#                        MOI.set(mock, SetProg.SumOfSquares.MomentMatrixAttribute(),
#                                cis[1], _moment_matrix([0.00336272, 0.0, 5.64559]))
#                        MOI.set(mock, SetProg.SumOfSquares.MomentMatrixAttribute(),
#                                cis[2], _moment_matrix([5.64559, 0.0, 0.00336272]))
#                    end
#                end
#            end
#            Tests.infeasible_switched_system_quad_test(bridged_mock(optimize!), config)
#        end
    end
    @testset "Quartic" begin
        @testset "Feasible" begin
			α = 11.814054544955727
            optimize!(mock) = MOIU.mock_optimize!(
                mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [α, 0.0, 0.0, -α, -0.0, α]))
            Tests.feasible_switched_system_quartic_test(bridged_mock(optimize!), config)
        end
        @testset "Infeasible" begin
            # Blocked for the same reason than quad
        end
    end
end


text:
using BSON
using LinearAlgebra
using SparseArrays

include("glance.jl")

function biplot(dataset, edges, inclsecond=true)
    data = BSON.load("$dataset-xy.bson")
    B = data["B"]
    xy = data["xy"]'
    labels = data["labels"]
    nnodes = length(labels)
    src, dst = findnz(B)[1:2]
    dst .+= nnodes

    # include hyperedges as labeled
    n = length(labels)
    xy1 = xy[1:nnodes, :]
    xy2 = xy[(nnodes + 1):end, :]

    outname = "viz-biplot-$dataset"
    if edges
        outname = "$(outname)-edges"
    end
    generate_nice_biplot(xy1, xy2,
                         src, dst,
                         "$outname.png",
                         mymarkersize1=3, mymarkersize2=1,
                         mymarkeralpha1=0.5, mymarkeralpha2=0.05,
                         incledges=edges,
                         inclsecond=inclsecond,
                         labels=labels,
                         )
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: PhilChodrow/HypergraphModularity

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
#Julia Parallel Processing

addproc(1)
srand(111)
r = remotecall(rand, 2, 3, 4)
s = @spawnat 2 1 .+ fetch(r)
fetch(s)