slwu89

# from question: https://discourse.julialang.org/t/speeding-up-jump-model-creation-with-sets-that-depend-on-other-indexes/107333

using Catlab, JuMP
using BenchmarkTools

@present NetSch(FreeSchema) begin
    (V,E)::Ob
    src::Hom(E,V)
    tgt::Hom(E,V)
    TimeType::AttrType

slwu89

using StochasticAD, Distributions, StaticArrays, Plots
using Zygote, ForwardDiff

# map rates to probabilities
function rate_to_proportion(r, t)
    1-exp(-r*t)
end

# dynamic part of SIR model
function sir_dyn_mod(x, u0, p, δt)

slwu89

rpwexp <- function(n, failRates) {
  end_time <- cumsum(failRates$duration)
  n_rate <- nrow(failRates)
  t <- 0
  tau <- Tk <- rep(0, n)
  Pk <- rexp(n = n)
  notfired <- rep(TRUE, n)
  for (i in 1:n_rate) {
    delta_tk <- (Pk[notfired] - Tk[notfired]) / failRates$rate[i]
    tau[notfired] <- t + delta_tk

slwu89

using MeasureTheory
using Plots
using Random

function counting(x)
    Dirac([x[1]+1, x[2]-1])
end

mc = Chain(x -> counting(x), Dirac([0,0]))
r = rand(mc)

slwu89

# 1. ----------
# "guess what an appropriate container eltype would be for storing results of `f`..."
Base.promote_op

f() = true
Base.promote_op(f)

# 2. ----------
# take apart whatever comes in
dump(f)

slwu89

rm(list=ls());gc()

ad <- function(x,diff = c(1,1)) {
  ## create class "ad" object. diff[1] is length of grad
  ## diff[2] is element of grad to set to 1.
  grad <- rep(0,diff[1])
  if (diff[2]>0 && diff[2]<=diff[1]) grad[diff[2]] <- 1
  attr(x,"grad") <- grad
  class(x) <- "ad"
  x

slwu89

using Distributions
using ForwardDiff
using Optim

using LineSearches # https://github.com/JuliaNLSolvers/Optim.jl/issues/713

# logdet
using LinearAlgebra

# get the data

slwu89

using Random, Distributions
using LinearAlgebra
using Statistics
using Plots
using JuMP, Ipopt
using RCall

# initial fixup and optional (not here) scaling https://github.com/cran/mgcv/blob/a534170c82ce06ccd8d76b1a7d472c50a2d7bbd2/R/smooth.r#L1602
# scaling here: https://github.com/cran/mgcv/blob/a534170c82ce06ccd8d76b1a7d472c50a2d7bbd2/R/smooth.r#L3757
function penalty_cr(D,B)

slwu89

using JuMP
using HiGHS
using LinearAlgebra
using RCall
using Plots

# example from pcls help, see https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/pcls.html

R"""
library(mgcv)

slwu89

#include <vector>
#include <list>
#include <iostream>
#include <algorithm>

// if we needed containers that took > 2 template args (set?) we could use variadic templates.
template<typename T, template<typename, typename> class C>
struct MyVariable {
    std::vector<C<T, std::allocator<T>>> values;