## ----setup, include=FALSE-------------------------------------------------------------------------
knitr::opts_chunk$set(collapse=TRUE, cache=FALSE)

## ----set-options, echo=FALSE----------------------------------------------------------------------
options(width=1e2)

## ----eval=FALSE, include=FALSE--------------------------------------------------------------------
# # We can use this to type 'rmd' in the console, and the markdown html document is then knitted fast.
# render.rmd("rmd")

## -------------------------------------------------------------------------------------------------
## Load libraries
library(ctsmTMB)
library(ggplot2) ## plots

## Create model
model <- newModel()
model$addSystem(dx ~ theta * (t*u^2-cos(t*u) - x) * dt + sigma_x*dw)
model$addObs(y ~ x)
model$setVariance(y ~ sigma_y^2)
model$addInput(u)

## Set parameter values
## note: not strictly necessary to set lower/upper bounds
model$setParameter(
  theta   = c(initial = 2, lower = 0,    upper = 100),
  sigma_x = c(initial = 0.2, lower = 1e-5, upper = 5),
  ## fix sigma_y to 0.05 by not giving any upper/lower bounds
  sigma_y = c(initial = 5e-2)
)

## Set initial state mean and covariance
## note: diag(1) is not strictly needed
model$setInitialState(list(1, 1e-1*diag(1)))

## ----collapse=TRUE--------------------------------------------------------------------------------
## set true parameters, and create data
true.pars <- c(theta=20, sigma_x=1, sigma_y=0.05)
dt.sim <- 1e-3
t.sim <- seq(0, 1, by=dt.sim)

## seed for input creation
set.seed(20)
u.sim <- cumsum(rnorm(length(t.sim),sd=0.1))
df.sim <- data.frame(t=t.sim, y=NA, u=u.sim)

## ----collapse=TRUE--------------------------------------------------------------------------------
## Set rng seeds for C++ states and observations
cpp.seeds <- c(20,20)

## perform simulation
sim <- model$simulate(data=df.sim, 
                      pars=true.pars, 
                      k.ahead = nrow(df.sim)-1, ## default
                      n.sims = 2,
                      cpp.seeds = cpp.seeds)


## -------------------------------------------------------------------------------------------------
mat <- as.matrix(data.frame(sim$times$i0, x=sim$states$x$i0, y=sim$observations$y$i0))
head(mat)
tail(mat)

## ----collapse=TRUE--------------------------------------------------------------------------------
## Extract all observations
y.sim <- sim$observations$y$i0[,1]

## Only select every tenth, to reduce the available information
iobs <- seq(1, length(t.sim), by=10)
t.obs <- t.sim[iobs]
y.obs <- y.sim[iobs]
u.obs <- u.sim[iobs]

## Create data for re-estimation
df.obs <- data.frame(
  t = t.obs,
  u = u.obs,
  y = y.obs
)

## Try to estimate the parameters
fit <- model$estimate(df.obs)

## ----collapse=TRUE--------------------------------------------------------------------------------
fit

## ----fig.height=5, fig.width=9, out.width="100%", fig.align='center'------------------------------
sim <- model$simulate(data=df.sim, pars=true.pars, n.sims=5, cpp.seeds = cpp.seeds)
x <- sim$states$x$i0
t <- sim$times$i0
matplot(t[,"t.j"], x, type="l", lty="solid", ylim=c(-4,4), xlab="Time")

## ----fig.height=5, fig.width=9, out.width="100%", fig.align='center'------------------------------
new.pars <- true.pars
new.pars["sigma_x"] <- 4
sim <- model$simulate(data=df.sim, pars=new.pars, n.sims=5, cpp.seeds=cpp.seeds)
x <- sim$states$x$i0
t <- sim$times$i0
matplot(t[,"t.j"], x, type="l", lty="solid", ylim=c(-4,4), xlab="Time")

## ----fig.height=5, fig.width=9, out.width="100%", fig.align='center'------------------------------
new.pars["sigma_x"] <- 4
sim <- model$simulate(data=df.sim, pars=new.pars, n.sims=100, cpp.seeds=cpp.seeds)

## quantiles
p <- c(0.01, 0.05, seq(0.1,0.9,by=0.1), 0.95, 0.99)
Q <- t(apply(sim$states$x$i0, 1,function(x) quantile(x, probs=p)))

## create data for distribution plot
p.center <- p[-length(p)] + diff(p)/2
col.ids <- c(col(Q[,-1]))
row.ids <- c(row(Q[,-1]))
fan.df <- data.frame(
  t = sim$times$i0[row.ids,"t.j"],
  ymin = c(Q[,-ncol(Q)]),
  ymax = c(Q[,-1]),
  ids = col.ids,
  ## symmetric colors around median
  fill.value = abs(0.5 - p.center[col.ids])
)

## Create plot
ggplot() +
  geom_ribbon(data=fan.df, aes(x=t, ymin=ymin, ymax=ymax, fill=fill.value, group=ids)) +
  geom_line(aes(x=sim$times$i0[,"t.j"], y=Q[,"50%"]), color="black", linewidth=0.3) +
  scale_fill_gradientn(colors=c("red","yellow")) +
  coord_cartesian(ylim=c(-4,4)) +
  guides(fill="none") +
  labs(x="Time",y="") +
  theme_minimal()

## ----eval=FALSE-----------------------------------------------------------------------------------
# model$simulate(data,
#                pars = NULL,
#                use.cpp = TRUE,
#                cpp.seeds = NULL,
#                method = "ekf",
#                ode.solver = "rk4",
#                ode.timestep = diff(data$t),
#                simulation.timestep = diff(data$t),
#                k.ahead = nrow(data)-1,
#                return.k.ahead = 0:min(k.ahead, nrow(data)-1),
#                n.sims = 100,
#                ukf.hyperpars = c(1, 0, 3),
#                initial.state = self$getInitialState(),
#                estimate.initial.state = private$estimate.initial,
#                silent = FALSE,
#                ...)