model_edo.Rmd

---
title: "Équations différentielles avec R"
author: "State of the R"
date: "24-28/08/2020"
output: html_document
---

[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/stateofther/finistR2020/website?filepath=finistR2020_website%2Fedo.Rmd)

# Differential Equations

## CRAN Task View

Le CRAN propose une [Task View](https://cran.r-project.org/web/views/DifferentialEquations.html) sur les équations différentielles.

Une équation différentielle est une équation entre une ou plusieurs fonctions inconnues et leurs dérivées. Elle décrit comment une fonction varie par rapport à une ou plusieurs variables (souvent le temps et/ou l'espace) et par rapport à ses dérivées.

Il y a différentes facon de classifier les équations différentielles :

- les équations peuvent être stochastique (la quantité inconnue est aléatoire) ou déterministe (la quantité inconnue est déterministe).

- les équations peuvent porter sur des fonctions à une seule variable (équation différentielle ordinaire) ou à plusieurs variables (équation aux dérivées partielles).

- les équations peuvent inclure des fonctions dont la dérivée à un certain pas de temps dépend de la dérivée à un pas de temps précédent (équation différentielle à retard ou differential equations delay). Elles peuvent aussi inclure des relations algébriques entre les variables (équation différentielle algébrique).


Il existe plusieurs packages R permettant de résoudre ces équations et d'ajuster ces modèles à de la donnée. Ici, seuls les package deSolve et diffeqr sont utilisés pour résoudre des ED.

## Résolution numérique d'EDO et d'EDP

Le package [`deSolve`](https://cran.r-project.org/web/packages/deSolve/) permet de résoudre numériquement des équations simples.

### Installation et chargement

```{r, eval = FALSE}
install.packages("deSolve")
```

```{r, message = FALSE}
library(tidyverse)
library(GGally)
library(deSolve)
library(fda)
theme_set(theme_bw())
```

### EDO simple: la désintégration atomique

On cherche à résoudre $y' = ay$ avec condition initiale $y(0) = y_0$. On commence par coder l'équation différentielle:
- `t` représente le temps courant
- `Y` représente l'état courant du système
- `parameters` stocke les paramètres du modèle (ici $a$)

```{r}
model <- function(t, Y, parameters) {
  with(as.list(parameters), {
    dy = -a * Y
    list(dy)
  })
}
```

On renseigne ensuite la jacobienne $\frac{\partial y'}{\partial y}$

```{r}
jac <- function(t, Y, parameters) {
  with(as.list(parameters), {
    PD[1, 1] <- a
    return(PD)
  })
}
```

On peut ensuite résoudre l'EDO pour $a = 1$ et $y_0 = 1$ sur l'intervalle $[0, 1]$ des pas de temps de longeur $0.01$ comme suit:
```{r}
params <- c(a = 1)
y0     <- c(1)
times  <- seq(0, 1, by = 0.01)
PD     <- matrix(0, nrow = 1, ncol = 1)
out_atome <- ode(y0, times, model, parms = params, jacfun = jac)
```

Le résultat est une matrice:
- une colonne pour le temps (reprend les valeurs de `times`)
- une colonne par dimension dans le système d'équations différentielles

On peut vérifier que la solution numérique (en bleu) est confondue avec la solution analytique (en rouge).

```{r}
plot_data <-
  data.frame(out_atome) %>%
  rename(numeric = X1) %>%
  mutate(analytic = exp(-time)) %>%
  pivot_longer(cols = -time,
               names_to = "type",
               values_to = "y")
ggplot(plot_data, aes(x = time, y = y, color = type)) +
  geom_line() +
  ylim(0, 1) +
  theme(legend.position = c(0.95, 0.95),
        legend.justification = c(1, 1),
        legend.background = element_rect(fill = NA))
```

`ode` utilise la méthode de Runge-Kutta pour calculer $y(t)$ et renvoie uniquement les valeurs $y(0), y(0.01), \dots, y(1)$ mais fait néanmoins appel à des points intérmédiaires lors du calcul. On peut s'en convaincre en comparant les valeurs finales ($y(1)$) obtenues avec `times = seq(0, 1, by = 0.01)` et `times = seq(0, 1, by = 1)`.

On peut diagnostiquer la réussite de l'intégration numérique via la commande `diagnostics`

```{r}
diagnostics(out_atome)
```

### EDO classique: le modèle SIR

Source <https://kinglab.eeb.lsa.umich.edu/480/nls/de.html>
Nous considérons le modèle différentiel classique en épidémiologie.

```{r}
# Define differential System
closed.sir.model <- function (t, state, parameters) {
  ## first extract the state variables
  S <- state[1]
  I <- state[2]
  R <- state[3]
  ## now extract the parameters
  beta <- parameters[1]
  gamma <- parameters[2]
  N <- S + I + R
  ## now code the model equations
  dSdt <- -beta * S * I/N
  dIdt <- beta * S * I/N - gamma * I
  dRdt <- gamma * I
  ## combine results into a single vector
  dxdt <- c(dSdt,dIdt,dRdt)
  ## return result as a list!
  list(dxdt)
}
```

Define parameters and times of evaluation

```{r}
parms <- c(beta = 400,gamma =  365/13)
# times stamps
times <- seq(from = 0,to = 60/365,by = 1/365/4)
# initial conditions
xstart <- c(S = 999,I = 1,R = 0)
```

Solving with ODE
```{r}
out.SIR <-
  ode(func = closed.sir.model,
      y = xstart,
      times = times,
      parms = parms,
      method = 'lsodar') %>%
  as.data.frame()
```
And plot
```{r}
out.SIR %>%
  gather(variable,value,-time) %>%
  ggplot(aes(x = time,y = value,color = variable)) +
  geom_line(size = 2) +
  labs(x = 'time (yr)',y = 'number of individuals')
```

### EDO classique: Modèle de Lorenz

Il s'agit d'une modélisation idéalisée de l'atmosphère. X, Y et Z représentent respectivement les variations verticale et horizontale de la température et le flux de convection.

$$
\left\{
\begin{align}
X' &= aX + YZ \\
Y' &= b(Y-Z) \\
Z' &= -XY + cY -Z \\
\end{align}
\right.
$$

```{r}
times <- seq(0, 100, by = 0.01)

a <- -8/3
b <- -10
c <- 28

lorenz <- function(t, y, parms) {
  with(as.list(y), {
    dX <- a * X + Y * Z
    dY <- b * (Y - Z)
    dZ <- -X * Y + c * Y - Z
    list(c(dX, dY, dZ))
  })
}

out_lorenz <- ode(y = c(X = 1, Y = 1, Z = 1), times = times, func = lorenz, parms = NULL)

df_lorenz <-
  out_lorenz %>%
  as_tibble() %>%
  mutate_all(as.numeric)

df_lorenz

ggpairs(df_lorenz,
        lower = list(continuous = wrap("points", size = 0.1)))
```

### EDO classique: Réaction de Belousov-Zhabotinskii

Il s'agit d'une récation chimique exhibant des motifs périodiques.
<center>
<iframe width="560" height="315" src="https://www.youtube.com/embed/PpyKSRo8Iec" frameborder="0" allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
</center>

$$
\begin{align}
A+B \xrightarrow{k_A} 2A \\
B+C \xrightarrow{k_B} 2B \\
C+A \xrightarrow{k_C} 2C \\
\end{align}
$$

$$
\left\{
\begin{align}
\partial_t n_A = n_A (k_A n_B - k_C n_C) \\
\partial_t n_B = n_B (k_B n_C - k_A n_A) \\
\partial_t n_C = n_C (k_C n_A - k_B n_B) \\
\end{align}
\right.
$$

```{r}
times <- seq(0, 30, by = 0.01)

k_A <- 1
k_B <- 1
k_C <- 2

bz <- function(t, y, params) {
  with(as.list(y), {
    dA <- A * (k_A * B - k_C * C)
    dB <- B * (k_B * C - k_A * A)
    dC <- C * (k_C * A - k_B * B)
    list(c(dA, dB, dC))
  })
}

out_bz <- ode(y = c(A = 0.6, B = 0.6, C = 0.3), times = times, func = bz, parms = NULL)

out_bz %>%
  as_tibble() %>%
  mutate_all(as.numeric) %>%
  pivot_longer(-time, names_to = "variable", values_to = "value") %>%
  ggplot() +
  aes(x = time, y = value, color = variable) +
  geom_line()
```

### EDP classique

La résolution d'EDP est plus complexe et consiste à transformer les EDPs en EDOs en utilisant la méthode des différences finies.

On illustre la démarche sur un modèle de diffusion de pestes (des aphides) sur une rangée de plantes positionnées entre $x = 0$ et $x = 60$.

$$
\frac{\partial N}{\partial t} = - \frac{\partial F}{\partial x} + g \times N
$$

où le flux diffusif est donné par $F = -D \frac{\partial N}{\partial x}$ et les contraintes stipulent que la densité d'aphides tombe à $0$ aux deux bouts de la rangée de plantes: $\forall t \geq 0, N(x=0, t) = N(x=60, t) = 0$.

Au temps initial, les aphides ne sont présentes qu'au milieu de la rangée de plantes:
$$
N(x, t = 0) = \begin{cases} 1 & \text{if} \quad x = 30 \\ 0 & \text{else} \end{cases}
$$

La méthodes des différences finies consiste à diviser le domaine en boîtes et à discrétiser l'équation comme suit
$$
\frac{dN_i}{dt} = - \frac{F_{i, i+1} - F_{i-1, i}}{\Delta x_i} + g\times N_i
$$
où $N_i$ représente la densité d'aphide au milieu de la boîte tandis que les flux sont définies aux interfaces entre boîtes:
$$
F_{i-1, i} = -D_{i, i-1} \times \frac{N_i - N_{i-1}}{\Delta x_{i-1, i}}
$$
pour se ramener à un système d'EDOs.

On commence par définir les équations du modèle

```{r}
Aphid <- function(t, APHIDS, parameters) {
  with(as.list(parameters), {
    ## taille des boîtes
    deltax  <- c(0.5, rep(1, numboxes - 1), 0.5)
    ## valeurs du flux aux interfaces
    Flux    <- -D * diff(c(0, APHIDS, 0)) / deltax
    ## valeurs des dérivées
    dAPHIDS <- -diff(Flux) / delx + APHIDS * r
    # the return value
    list(dAPHIDS)
  })
}
```

puis les paramètres du modèle et de la grille de discretisation:

```{r}
params <- list(
    D        <- 0.3,    # m2/day  diffusion rate
    r        <- 0.01,   # /day    net growth rate
    delx     <- 1,      # m       thickness of boxes
    numboxes <- 60,
    # distance of boxes on plant, m, 1 m intervals
    Distance <- seq(from = 0.5, by = delx, length.out = numboxes)
)
```

et enfin les conditions intiales
```{r}
# Initial conditions:  # ind/m2
APHIDS        <- rep(0, times = numboxes)
APHIDS[30:31] <- 1
# initialise state variables
state         <- c(APHIDS = APHIDS)
```

on peut alors laisser la densité évoluer sur 200 jours avec un rendu par jour:
```{r}
out_aphide <- ode.1D(state, times = 0:200, Aphid, parms = params,
              nspec = 1, names = "Aphid")
```

Le résultat est comme auparavant une matrice avec la première colonne pour le temps et une colonne par EDO dans la discrétisation.

On peut utiliser ce résultat pour voir la diffusion des aphides dans les plantes au cours du temps:
```{r}
out_aphide %>%
  data.frame %>%
  pivot_longer(cols = -time,
               names_to = "location",
               names_pattern = "APHIDS(.*)",
               values_to = "density") %>%
  mutate(location = as.integer(location)) %>%
  ggplot(aes(x = time, y = location, fill = density)) +
  geom_tile() +
  scale_fill_viridis_c()
```

### Equation différentielle à retard

Cet exemple est tiré du document : https://cran.r-project.org/web/packages/deSolve/vignettes/deSolve.pdf .

Le modèle logistique est un modèle de dynamique de population simple dans lequel une population donnée évolue jusqu'à atteindre un plateau (la capacité biotique du milieu noté M). La dynamique de la population peut être représentée par l'ED suivante :  $x'(t)=r.x(t).[\frac{1-x(t)}{M}]$. Le terme densité-dépendant $\frac{1-x(t)}{M}$ peut avoir un effet sur la population avec un temps de retard. L'ED peut alors être réécrite : $x'(t)=r.x(t).[\frac{1-x(t - \tau)}{M}]$. Sous cette forme l'ED est une équation différentielle à retard (EDR).

Dans le cas d'une population de lemming, le lag est fixé à 9 mois (0.74 année), soit le temps de croissance d'un lemming jusqu'au stade adulte. Le paramètre r est fixé sur la base de données expérimentales et la capacité biotique du milieu (M) est fixée arbitrairement à 19.

Les lignes de codes qui suivent permettent de simuler et de visualiser l'évolution d'une population de lemming ($1^{ère}$ figure). La deuxième figure représente l'évolution de $x(t)$ en fonction de $x(t-\tau)$.

```{r}
# DDE function
derivs <- function(t, y, parms) {
  if (t < 0) {
    lag <- 19
  } else {
    lag <- lagvalue(t - 0.74)
  }
  dy <- r * y * (1 - lag / m)
  list(dy, dy = dy)
}

# parameters
r <- 3.5; m <- 19

# initial values and times
yinit <- c(y = 19.001)
times <- seq(-0.74, 40, by = 0.01)

yout <- dede(y = yinit, times = times, func = derivs,
             parms = NULL, atol = 1e-10)

df_out <-
  yout %>%
  as_tibble() %>%
  mutate_all(as.numeric) %>%
  mutate(yend = lead(y), dy.yend = lead(dy.y))

ggplot(df_out) +
  aes(x = time, y = y) +
  geom_line()

ggplot(df_out) +
  aes(x = y, y = dy.y, xend = yend, yend = dy.yend) +
  geom_segment() +
  labs(x = "y", y = "y'")
```


## diffeqr

Il s'agit d'une interface R pour utiliser le package Julia [DifferentialEquations.jl](https://diffeq.sciml.ai/dev/). Le package est disponible sur le [CRAN](https://CRAN.R-project.org/package=diffeqr) et sur [GitHub](https://github.com/SciML).

### Installation

Pour l'utiliser il faut installer Julia, en le téléchargeant [ici](https://julialang.org/downloads/) et le package DifferentialEquations.jl. Attention le package R JuliaCall qui permet d'appeler des fonctions julia depuis R ne fonctionne pas avec la version 4.0 de R au moment du test donc utilisez la version 3.6.

## Installation des packages Julia

Je conseille de tester le package Julia dans un premier temps. L'instruction
qui suit permet d'exécuter du code julia dans un fichier Rmarkdown:

```{r}
julia <- JuliaCall::julia_markdown_setup()
```
- Installation de DifferentialEquations.jl et Plots.jl

```{julia}
using Pkg
Pkg.add(["Plots", "DifferentialEquations"])
```

- Résolution du système de Lorenz
```{julia}
using DifferentialEquations, Plots

function lorenz(du,u,p,t)
     du[1] = p[1]*(u[2]-u[1])
     du[2] = u[1]*(p[2]-u[3]) - u[2]
     du[3] = u[1]*u[2] - p[3]*u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(lorenz,u0,tspan,p)
sol = solve(prob, saveat=0.01)
plot(sol,vars=(1,2,3))
```

### Modèle de Lorenz

Le package diffeqr permet d'utiliser le package julia depuis R:
```{r}
de <- diffeqr::diffeq_setup()
```

Définition de la fonction dérivée:

```{r}
f <- function(u,p,t) {
  du1 = p[1]*(u[2]-u[1])
  du2 = u[1]*(p[2]-u[3]) - u[2]
  du3 = u[1]*u[2] - p[3]*u[3]
  return(c(du1,du2,du3))
}
```

Les paramètres sont contenus dans un vecteur `p`.

```{r}
u0 <- c(1.0,0.0,0.0)
tspan <- c(0.0,100.0)
p <- c(10.0,28.0,8/3)
prob <- de$ODEProblem(f, u0, tspan, p)
sol <- de$solve(prob, saveat=0.01)
```

Les méthodes numériques pour la résolution sont très nombreuses. Voir la page [ODE solvers](https://diffeq.sciml.ai/dev/solvers/ode_solve/).

La solution `sol$u` est une liste de vecteurs, et `sol$u[i]` est le vecteur u[i] au temps `sol$t[i]`. On peut le transformer en matrice R avec `sapply`:

```{r}
mat <- sapply(sol$u,identity)
```

Chaque ligne est une série temporelle, on peut transformer la solution en data.frame

```{r}
udf <- as.data.frame(t(mat))
matplot(sol$t,udf,"l",col=1:3)
```

```{r}
plotly::plot_ly(udf, x = ~V1, y = ~V2, z = ~V3, type = 'scatter3d', mode = 'lines')
```

### Ameliorer les performances

- Compilation de la fonction dérivée: Chaque calcul de la dérivée est un appel d'une fonction R. `diffeqr` propose une fonction qui permet de compiler cette fonction. La fonction R sera traduite en Julia et compilée dans le langage LLVM qui est un dialecte du langage C. Ce langage et sa compilation sont le "moteur" de Julia qui lui permet d'avoir des performances comparables aux langages compilés comme le C++ ou le Fortran.

```{r}
prob2 <- diffeqr::jitoptimize_ode(de,prob)
sol <- de$solve(prob2, de$Tsit5() )
```
On fixe ici la méthode de résolution avec Tsit5 (Tsitouras 5/4 Runge-Kutta).

Pour améliorer un peu plus les performances, on peut utiliser la version Julia de la fonction dérivée. Transformons toutes les variables R en variables Julia:

```{r}
julf <- JuliaCall::julia_eval("
function julf(du,u,p,t)
  du[1] = 10.0*(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
  du[3] = u[1]*u[2] - (8/3)*u[3]
end")
JuliaCall::julia_assign("u0", u0)
JuliaCall::julia_assign("p", p)
JuliaCall::julia_assign("tspan", tspan)
prob3 <- JuliaCall::julia_eval("ODEProblem(julf, u0, tspan, p)")
sol <- de$solve(prob3,de$Tsit5())
```

```{r}
system.time({ for (i in 1:100){ de$solve(prob, de$Tsit5()) }})
```

```{r}
system.time({ for (i in 1:100){ de$solve(prob2, de$Tsit5()) }})
```

```{r}
system.time({ for (i in 1:100){ de$solve(prob3, de$Tsit5()) }})
```

### Stochastic Differential Equation (SDE) Examples

Les performances de `diffeqr` seront plus intéressantes dans le cas des EDS. En effet, pour résoudre numériquement ce type d'équation, nous avons besoin de calculer la solution d'une EDO plusieurs milliers de fois. Le gain qui peut paraître minime pour une EDO devient mille fois plus intéressant pour une EDS.

On utilise deux fonctions `f` et `g`, où `du = f(u,t)dt + g(u,t)dW_t`

## Résolution d'une SDE

```{r}
f <- function(u,p,t) {
  du1 = p[1]*(u[2]-u[1])
  du2 = u[1]*(p[2]-u[3]) - u[2]
  du3 = u[1]*u[2] - p[3]*u[3]
  return(c(du1,du2,du3))
}
g <- function(u,p,t) {
  return(c(0.3*u[1],0.3*u[2],0.3*u[3]))
}
u0 <- c(1.0,0.0,0.0)
tspan <- c(0.0,1.0)
p <- c(10.0,28.0,8/3)

tspan <- c(0.0,100.0)
prob <- de$SDEProblem(f,g,u0,tspan,p)
fastprob <- diffeqr::jitoptimize_sde(de,prob)
sol <- de$solve(fastprob,saveat=0.005)
udf <- as.data.frame(t(sapply(sol$u,identity)))
plotly::plot_ly(udf, x = ~V1, y = ~V2, z = ~V3, type = 'scatter3d', mode = 'lines')
```

### Version avec `desolve` 

```{r}
library(deSolve)
Lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dX <-  a * X + Y * Z
    dY <-  b * (Y - Z)
    dZ <- -X * Y + c * Y - Z
    list(c(dX, dY, dZ))
  })
}

parameters <- c(a = -8/3, b = -10, c = 28)
state      <- c(X = 1, Y = 1, Z = 1)
times      <- seq(0, 100, by = 0.01)
out <- ode(y = state, times = times, func = Lorenz, parms = parameters)

lorenz_solve <- function (i){
  state      <- c(X = runif(1), Y = runif(1), Z = runif(1))
  parameters <- c(a = -8/3 * runif(1), b = -10 * runif(1), c = 28 * runif(1))
  out <- ode(y = state, times = times, func = Lorenz, parms = parameters)
}
```
```{r}
system.time({ lapply(1:1000,lorenz_solve) })
```

```{r}
prob_func <- function (prob,i,rep){
  de$remake(prob,u0=runif(3)*u0,p=runif(3)*p)
}
ensembleprob = de$EnsembleProblem(fastprob, prob_func = prob_func, safetycopy=FALSE)
```

```{r}
system.time({ de$solve(ensembleprob,de$Tsit5(),de$EnsembleSerial(),trajectories=1000,saveat=0.01) })
```

Le package `diffeqr` fonctionne également sur GPU, vous pouvez voir un exemple dans ce [billet de blog](http://www.stochasticlifestyle.com/gpu-accelerated-ode-solving-in-r-with-julia-the-language-of-libraries/).

## Estimation des paramètres

On chercher mainenant à estimer les paramètres structurels d'un système diférentiel à partir d'observations bruitées. On l'applique sur le SIR.

```{r}
sigma <-  20
noisy <- out.SIR
noisy$S  <- pmax(out.SIR$S  + rnorm(length(out.SIR$time),0,sigma),0)
noisy$I  <- pmax(out.SIR$I  + rnorm(length(out.SIR$time),0,sigma),0)
noisy$R  <- pmax(out.SIR$R  + rnorm(length(out.SIR$time),0,sigma),0)


noisy %>%
  gather(variable,value,-time) %>%
  ggplot(aes(x = time,y = value,color = variable)) +
  geom_point(size = 2) +
  theme_classic() +
  labs(x = 'time (yr)',y = 'number of individuals')

observ <- noisy[,2:4]
```

Utiliser une méthode des moindres carrés habituelles est rendue compliquée par le fait que la fonction de régression $f$ est solution d'une quat diff donc n'a pas de solution explicite. De plus, chaque évaluation de la fonction (par un schéma numérique) est couteux d'un point de vue computationnel.

Une autre approche est d'exprimer $f$ dans une base (spline par exemple), cette fonction devra d'une part s'ajuster aux donnnées et d'autre part être solution du système (qui est introduit comme une pénalité).
Cette méthode qui date de 2007, est implémentée dans le package `pcode`
La méthode est décrite dans la [vignette](https://cran.r-project.org/web/packages/pCODE/vignettes/pcode-vignette.html) du pacakge.

On définit dabord une base de décomposition pour les fonctions (package fda)
```{r}
#" basis list
knots <- seq(0,max(times),length.out=21)
#order of basis functions
norder <- 4
#number of basis funtions
nbasis <- length(knots) + norder - 2
#creating Bspline basis
basis_dim1  <- create.bspline.basis(c(0,max(times)),nbasis,norder,breaks = knots)
basis = list(basis_dim1,basis_dim1,basis_dim1)
```

Puis on optimise: 

```{r, eval=FALSE}
pcode.result <- pcode(data = observ, time = times, ode.model = closed.sir.model,
                      par.initial = c(100,1), par.names = c('beta','gamma'),state.names = c('S','I','R'),
                      basis.list = basis, lambda = 1e2)
```

## Jeton de reproductilité

```{r session-information}
sessionInfo()
```