From 163dc48022cc9dae1347a397de04cf8b32b4ef17 Mon Sep 17 00:00:00 2001
From: milesmedina <44986194+milesmedina@users.noreply.github.com>
Date: Fri, 6 Sep 2024 19:32:26 -0400
Subject: [PATCH] phase space embeddings for pyro & chl

---
 R/embed_chl_otb.R           |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_chl_otb_ce.R        |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_chl_otb_cw.R        |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_chl_otb_ne.R        |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_chl_otb_nw.R        |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_chl_otb_se.R        |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_chl_otb_sw.R        |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_pyro_otb-ce.R       |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_pyro_otb-cw.R       |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_pyro_otb-ne.R       |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_pyro_otb-nw.R       |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_pyro_otb-se.R       |  94 ++++++++++++++++++++++++++++++++++++
 R/embed_pyro_otb-sw.R       |  94 ++++++++++++++++++++++++++++++++++++
 data/embed_chl_otb.RData    | Bin 0 -> 5717 bytes
 data/embed_chl_otb_ce.RData | Bin 0 -> 5722 bytes
 data/embed_chl_otb_cw.RData | Bin 0 -> 5695 bytes
 data/embed_chl_otb_ne.RData | Bin 0 -> 5723 bytes
 data/embed_chl_otb_nw.RData | Bin 0 -> 5722 bytes
 data/embed_chl_otb_se.RData | Bin 0 -> 5716 bytes
 data/embed_chl_otb_sw.RData | Bin 0 -> 5729 bytes
 20 files changed, 1222 insertions(+)
 create mode 100644 R/embed_chl_otb.R
 create mode 100644 R/embed_chl_otb_ce.R
 create mode 100644 R/embed_chl_otb_cw.R
 create mode 100644 R/embed_chl_otb_ne.R
 create mode 100644 R/embed_chl_otb_nw.R
 create mode 100644 R/embed_chl_otb_se.R
 create mode 100644 R/embed_chl_otb_sw.R
 create mode 100644 R/embed_pyro_otb-ce.R
 create mode 100644 R/embed_pyro_otb-cw.R
 create mode 100644 R/embed_pyro_otb-ne.R
 create mode 100644 R/embed_pyro_otb-nw.R
 create mode 100644 R/embed_pyro_otb-se.R
 create mode 100644 R/embed_pyro_otb-sw.R
 create mode 100644 data/embed_chl_otb.RData
 create mode 100644 data/embed_chl_otb_ce.RData
 create mode 100644 data/embed_chl_otb_cw.RData
 create mode 100644 data/embed_chl_otb_ne.RData
 create mode 100644 data/embed_chl_otb_nw.RData
 create mode 100644 data/embed_chl_otb_se.RData
 create mode 100644 data/embed_chl_otb_sw.RData

diff --git a/R/embed_chl_otb.R b/R/embed_chl_otb.R
new file mode 100644
index 0000000..3d3f2d8
--- /dev/null
+++ b/R/embed_chl_otb.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_chl_otb.RData" )
+  dat <- ssa_chl_otb$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_chl_otb <- list( x = x,
+                         t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_chl_otb, file = "../data/embed_chl_otb.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_chl_otb_ce.R b/R/embed_chl_otb_ce.R
new file mode 100644
index 0000000..3340a9e
--- /dev/null
+++ b/R/embed_chl_otb_ce.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_chl_otb_ce.RData" )
+  dat <- ssa_chl_otb_ce$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_chl_otb_ce <- list( x = x,
+                            t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_chl_otb_ce, file = "../data/embed_chl_otb_ce.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_chl_otb_cw.R b/R/embed_chl_otb_cw.R
new file mode 100644
index 0000000..746e975
--- /dev/null
+++ b/R/embed_chl_otb_cw.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_chl_otb_cw.RData" )
+  dat <- ssa_chl_otb_cw$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_chl_otb_cw <- list( x = x,
+                            t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_chl_otb_cw, file = "../data/embed_chl_otb_cw.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_chl_otb_ne.R b/R/embed_chl_otb_ne.R
new file mode 100644
index 0000000..e7c530a
--- /dev/null
+++ b/R/embed_chl_otb_ne.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_chl_otb_ne.RData" )
+  dat <- ssa_chl_otb_ne$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_chl_otb_ne <- list( x = x,
+                            t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_chl_otb_ne, file = "../data/embed_chl_otb_ne.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_chl_otb_nw.R b/R/embed_chl_otb_nw.R
new file mode 100644
index 0000000..1cbbf09
--- /dev/null
+++ b/R/embed_chl_otb_nw.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_chl_otb_nw.RData" )
+  dat <- ssa_chl_otb_nw$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_chl_otb_nw <- list( x = x,
+                            t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_chl_otb_nw, file = "../data/embed_chl_otb_nw.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_chl_otb_se.R b/R/embed_chl_otb_se.R
new file mode 100644
index 0000000..deebd44
--- /dev/null
+++ b/R/embed_chl_otb_se.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_chl_otb_se.RData" )
+  dat <- ssa_chl_otb_se$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_chl_otb_se <- list( x = x,
+                            t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_chl_otb_se, file = "../data/embed_chl_otb_se.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_chl_otb_sw.R b/R/embed_chl_otb_sw.R
new file mode 100644
index 0000000..006d5a9
--- /dev/null
+++ b/R/embed_chl_otb_sw.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_chl_otb_sw.RData" )
+  dat <- ssa_chl_otb_sw$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_chl_otb_sw <- list( x = x,
+                            t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_chl_otb_sw, file = "../data/embed_chl_otb_sw.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_pyro_otb-ce.R b/R/embed_pyro_otb-ce.R
new file mode 100644
index 0000000..7d4e911
--- /dev/null
+++ b/R/embed_pyro_otb-ce.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_pyro_otb_ce.RData" )
+  dat <- ssa_pyro_otb_ce$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_pyro_otb_ce <- list( x = x,
+                             t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_pyro_otb_ce, file = "../data/embed_pyro_otb_ce.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_pyro_otb-cw.R b/R/embed_pyro_otb-cw.R
new file mode 100644
index 0000000..498b0b1
--- /dev/null
+++ b/R/embed_pyro_otb-cw.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_pyro_otb_cw.RData" )
+  dat <- ssa_pyro_otb_cw$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_pyro_otb_cw <- list( x = x,
+                             t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_pyro_otb_cw, file = "../data/embed_pyro_otb_cw.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_pyro_otb-ne.R b/R/embed_pyro_otb-ne.R
new file mode 100644
index 0000000..1b3d3c5
--- /dev/null
+++ b/R/embed_pyro_otb-ne.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_pyro_otb_ne.RData" )
+  dat <- ssa_pyro_otb_ne$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_pyro_otb_ne <- list( x = x,
+                             t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_pyro_otb_ne, file = "../data/embed_pyro_otb_ne.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_pyro_otb-nw.R b/R/embed_pyro_otb-nw.R
new file mode 100644
index 0000000..e1d84f9
--- /dev/null
+++ b/R/embed_pyro_otb-nw.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_pyro_otb_nw.RData" )
+  dat <- ssa_pyro_otb_nw$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_pyro_otb_nw <- list( x = x,
+                             t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_pyro_otb_nw, file = "../data/embed_pyro_otb_nw.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_pyro_otb-se.R b/R/embed_pyro_otb-se.R
new file mode 100644
index 0000000..c7e7c0f
--- /dev/null
+++ b/R/embed_pyro_otb-se.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_pyro_otb_se.RData" )
+  dat <- ssa_pyro_otb_se$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_pyro_otb_se <- list( x = x,
+                             t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_pyro_otb_se, file = "../data/embed_pyro_otb_se.RData" )
+  
\ No newline at end of file
diff --git a/R/embed_pyro_otb-sw.R b/R/embed_pyro_otb-sw.R
new file mode 100644
index 0000000..1cdce21
--- /dev/null
+++ b/R/embed_pyro_otb-sw.R
@@ -0,0 +1,94 @@
+# Phase space reconstruction
+# Miles Medina, ECCO Scientific, 2024
+#
+
+rm(list=ls(all=TRUE)) 
+
+# Load libraries
+  if(!require(tseriesChaos)) { install.packages('tseriesChaos') }; library(tseriesChaos)
+  if(!require(zoo)) { install.packages('zoo') }; library(zoo)
+  if(!require(plotly)) { install.packages('plotly') }; library(plotly)
+  if(!require(rEDM)) { install.packages('rEDM') }; library(rEDM)
+
+# Load data and select column
+  load( "../data/ssa_pyro_otb_sw.RData" )
+  dat <- ssa_pyro_otb_sw$dat
+  colnames( dat )
+  var <- 'signal'
+  x <- dat[,var] |> na.omit()
+  x <- scale( x )
+  t <- dat$month
+
+# Plot the time series
+  par(mfrow=c(2,1))
+  plot( x, type = 'l', xlab = 'time', main = "Time series x" )
+  abline( v = axTicks(1), col = rgb(0,0,0,0.2) )
+  abline( h = axTicks(2), col = rgb(0,0,0,0.2) )
+
+# Find embedding parameters
+  par(mfrow=c(3,1))
+  # Set embedding delay (d) using AMI function
+  ami <- mutual( x, lag.max = 20 )  # average mutual information function
+  local.d.min <- rollapply( ami, 3, function (x) which.min(x)==2 )  # local minima
+  d <- as.numeric( which(local.d.min==TRUE)[1] )  # first local min
+  points( x = d, y = ami[d+1], cex = 2, col = rgb(1,0,0,0.7) )
+  d
+
+  # Set Theiler window (tw) using ACF
+  ac <- acf( x, lag.max = 20 )  # autocorrelation function
+  ac$acf_abs <- ac$acf |> abs()  # absolute values of acf
+  local.tw.min <- rollapply(ac$acf_abs, 3, function (x) which.min(x)==2 )  # abs local minima
+  tw <- as.numeric( which(local.tw.min==TRUE)[1] )  # first abs local min
+  points( x = tw, y = ac$acf[tw+1], cex = 2, col = rgb(1,0,0,0.7) )
+  tw
+
+  # Set embedding dimension (m) using false nearest neighbors test
+  fnn <- false.nearest( series=x, m=6, d=d, t=tw, eps=sd(x), rt=10 )
+  threshold <- 0.15
+  plot( fnn[1,], type = 'b', main = "FNN", pch = 16, cex = 2,
+        xlab = 'dim', ylab = 'proportion of false neighbors' )
+    abline( h = threshold, lty = 2, col = rgb(0,0,0,0.6) )
+  m <- as.numeric( which( fnn[1,] <= threshold )[1] )
+  if( m==2 ){ m <- 3 }
+  points( x = m, y = fnn[1,m], cex = 3, col = rgb(1,0,0,0.7) )
+  
+  
+# Time-delay embedding
+  Mx <- embedd( x, m = m, d = d ) |> as.data.frame()
+
+# Plot phase space reconstruction
+  # Rename columns
+  for(i in 1:m){
+    if(i==1){ names(Mx)[i]<-'x(t)'
+    } else {
+      names(Mx)[i] <- paste0('x(t+',d*(i-1),')')
+    }
+  }  # // end i
+  # Plotly
+  p1 <- plot_ly( Mx, x = ~Mx[,1], y = ~Mx[,2], z = ~Mx[,3],
+                 type = 'scatter3d', mode = 'lines',
+                 opacity = 0.75, line = list(width = 6, reverscale = FALSE) ) |> 
+    layout( title = 'Reconstructed attractor',
+            scene = list( xaxis = list(title=names(Mx)[1]),
+                          yaxis = list(title=names(Mx)[2]),
+                          zaxis = list(title=names(Mx)[3])
+            ) )
+  p1
+  
+  
+# Test for nonlinear stationarity with space-time separation plots
+  par(mfrow=c(1,1))
+  stp <- stplot( series = x, m = m, d = d, mdt = length(x) )
+  
+
+# Export embedding outputs
+  embed_pyro_otb_sw <- list( x = x,
+                             t = t,
+                     Mx = Mx,
+                     d = d,
+                     tw = tw,
+                     m = m,
+                     plot = p1
+                     )
+  save( embed_pyro_otb_sw, file = "../data/embed_pyro_otb_sw.RData" )
+  
\ No newline at end of file
diff --git a/data/embed_chl_otb.RData b/data/embed_chl_otb.RData
new file mode 100644
index 0000000000000000000000000000000000000000..0a0354e37100cbb72879b6f3aad55b490d5450f3
GIT binary patch
literal 5717
zcmZveXD}R$w}<uKS6gkxiV$Jd=rwwc7Q!N75k&8uD67}JT6AKSDADT@H4(k{wTRB@
zUEKHH`{B-i=ALio%x~t*nR%X>1BxTU`kx{=G$%DNoS~P+CBxw;;?6FmA^tE>R#w98
ztR*BmVki=i1<%4IqYxGs?hrwh13#E%S%^-MiN4Z5)pIx+ZZ7K3)1yt9J5EMF;ds6*
z;V}!ejrwuE>TJ1^h!R6d?xP&Lu4L~}c1~i#9YOH7H>T8h?=>>B+=&dTZ&TJrPqWq}
zQ`v{Lrb*xe;1@G>03e$Z2`+x&x_miDG`rTp^gF`F2E8TmcQacc%KR<XzDClqB!M`u
zzM;}cRNvcDbi0gY0H!KI18O26TXSS3&Cm4xd`LNvD&JDv@%&TsZBfyeg-D-BVLd0y
zrTT4NpCS&7(qr{6X4m(_(7CFI>ubF+ebO<<{I${nLO?qms=JYaVT>hEHTSQ)5jwQw
zsKu(~_YORNtoxq~<Rv60n8EHQGQl+jOk$LXN|VwsO~~VNU1l<vEka0C_$`RZDn3uT
zu0bxi1A^q7cf+Exzg%zqeh+^&@cSe;gG)rzR0`cg=bbY)aMVZFUy~%o6j{4BI^>o%
zdGD45U-{)UM~QX{);3PGkPFku(KRJOT&YxyP8o+x@Uj^tw-UOT3oA7Hi)d;jQue}S
z0CnjJs|j%`&IaUe6y2g2`h(hp-I#rSLeax2$G+!N`F9;`{Jr?OLd!?n!AaiqA>IW`
z56GRPmZa+>RBM<+UB|w!Bh3V$J#?r%Ro#N+Dz$;O-K)$m&*z@vb7X~`yM<o>CMZYW
zN2<UtmJ;c6tQ6D}p#?pru9$aCh8q~02AMRYur13Zw|K=p2(}NR;}%Z_sP31DQ91E5
z%;TNro3i-}c{}J-nn!Y!xq`L~o9Az@m{))FP?9&g-u-iltR^S0Kie%_d?%Xw4;yYH
zkWSuFa-RO78$KI&GkO2h<`5E}YV99%tiE#ocuQ+m+j^=V_(9{XtD^lQfK{anK0+hE
zfcp-jN(=`JtbgYh0Cf?NXk5H^%r{aq(3KRN+bA+7$xdLIv?2($z5zoyfudad%)C&E
zvs8nPeeIT6?S=P!ZGaZ<tM3O(s)SZHRyyZ@d#ta|b-&!D$AF5KJh~cu+NP-g@%czg
zWMQxL6>$_VDSz||_yu9{hkkglncqnyVB^c2+aZN6Wo5?333y34oHSj0ZJ>bT>M-1*
zx6vz~pjuC-dnlmt?%V2i_!(r=$#HZc(mOXz1@@Go;1uItBfWSOoz0~qm4n<7a=5f`
zt9vobGSk)8rkOUydVW?*jK8$fhCcN%g%WTNzME<4z^%E-{h{$g6ta(5?d^xjpF+(p
z#BYz6Q+K^<HlFA6hh<7fdK{@-oW2Q;&@kAzeERH4Il*W>^RT9=AgWK`SI|}O&X#hB
zORI+n($4?l^WD1Ax4pwzr3mc-y09sdyF<Rd{{1%E>zAVgHhDU66&R7VQynzG39;JO
z>M?u!C-_wKIZ#@>S~CXUf4#<3V1JPYqGP0jGQw>y+RDwZ{_@5hYX%v_ApzTx+VQtu
zK0F*G4%{O-*(a2;a#iq2(N&yy3H5Gx=WTTE4sG>1F$zF*wrHUlkJnb+0~7X<qbMgp
z+ujE^9lvd^iuo+V?|dA@GuuyRp<PUod%WJ1YFTnD|8#;R?+~mRC*=!+dP*A^EwRnB
z{4KvZpaF;ashr~2sK?boG1xdB6t)OHst2Ps%zo6<)&#l(5i^JSLGv!)q0y60n{0=;
zy)nxNWMcDB95)Ho=(7K0+;vu(Yr8br3{tMZhFqvu&SEDc_pum<*8_R|uLL;8_owP7
z&?HZ&WFe}w<u0a$<H<Yy8vC~`-nqUH#<|<?W>%Ej?@UADr*h9KKS6O`-^Au*>jrU_
z{cvo*&h2nkqcb$O8(oUfpAB@~se&V0KHitN3I{>Z*@C<!aRA--GY^Q!EbBMu><roY
zeZ%i8b9t+I-Y1zPa~;bI82JfcU;sRPJHd>9XVd|d*#4WttF6_No_{$pFGA+}JA4~B
zy%X>?DX88=rYrAq2lPG`nLVk0><8oK-^**x6aI-;>7N5l7U5KrIl}K<$9(G(Ef;Io
ztg3Ros71S(LqPq?2^kxdT(oFdhWMiw@vkz^mUoL|)(+$0KTumBltIE?ef6&+8=Hb2
zy$|05YQac9zCp3RUAAgJj~_k#0UB;2R)6CAkieFM$&xAj3%DD4loXP6uoO%BplMxk
z2kDo<{pOTv?ox<jHZ=yN9Ff<qr2R!4s|88Bk3kf-U;?OJl1wiJUeFl@v1p2P{s|pJ
zPTPHZP}Ov$jwR)BVT5BwV%P1=Eut~G5|4|ko*qz9f8DihQc^dkIRy)lk;KO2?}gth
zK88i{)pz=`%MMX{Nwgp5X{Sg<@lEph&ec^KRbLJAw||bTLpFMeQ(uCS>3q&mkzXVX
zV>>~1n{^>*A&dp@SoDEwB&V>=p`%gDkW9_MG(RwdOx9-|;j?%ABOEv6t#cPE&;AhQ
zfKQctQo6P&r=CQ?)+u?H=gKXS=wYMIQ~b);Hp-{*5n@N{+VoQVImlH|hO@J$x~yUa
z>P&7pYjT?vD*P?(_hkgmH&f3!j9|s!470I&4$Q5^^KGi;1-(zl<NJ!SsOqq?N%!jZ
zI3q#-YuMsn9~;)!K?+^YV7ue-#$~XD82oThIX@HBD%6HOcktK$F!n13y~p23p7t!h
z|EsVRV)Nm1)NdA!`hUcv-IRUI?3K`V*O%io%3}so6XrL@Q+b>-g{6Ui|6FZRTd<sg
zp^q}6J9<QfA=|X!ByVgFzF^uIwqFM4WImv|psAgwOK-zdYuyhgU-7GhT`!Ij?=F*M
z68DXtBL*O5Vf%5o&cU^t;{0Y@x3}@d0LiPJ?uMd!ff6Q({)2BT_=fqu8nk(tiAQ5B
z7BA(D<0nC}%DqsPb!J>ZEM0=A$<4eA9!8M}iIVosse8eON_zZL-dQhK?Q2c^=)dL#
zonE^9;d0F!IJ;EpJL8h3_r49K4tMaF9-j`Jceinwwz9%pRbOeGU(Zkenh2mM)X48%
zF4eVNvqqreqpNg#o<>goN`Kb0b_ynZS)IY-x?n!nn|IdLM+eqQbZ%RVz~sNRn{e#)
zf(}(R(2~^a#V&=Bv#G+_|8ci|bydkcQJlw?xS_MG>J`KDmVsA-dr8>owr%lM%exLl
z@#$O$+I`;%{&~s;Jz0NQYJX$q`fKu7Vj*oC?fxe_hO&ulfX?1B;y44gCfz(_-XF~S
z`UhdqEh{DjN7R?WS)EaWRmQva8eJa+di_vWpK93~r1$sozD;_GB`5E*h!96;pK{jA
z%1(DW0@yD@ZW4re+O3`91Xw#6J3p25+Dd=t3HwtL{ZfLDW%7{0b(h1CmBxkcd8T+#
zMMuY)q2_ikzl7V?x4B=V;C{K}78-iuHT8OUd~yU2n&UoUTG;orKgqgD<!cBQ6gRd%
z%3l_qcvKf85l8@!CD@fi2ScIc*IeU1;nHN+z~qHc6mDl!`nH5f@G3n_uG^1&5J#5o
zoKhhc+mh(ipu0H~$*`%{Jr*{M-@yozNw5?sz998xXT=7QxU;Nb!q|w~@XnRl=WvkZ
zn^&PA!Uvk4%%K@!A4^W+SU)^i`6TOvQW8qA9Ru0JvGOTz^^^Cpb+NL@uGPm46jt#%
zSYUtz%R0P;v;hb?)&qGm04pCq6%UGi>z<sCiv(<9cMJX$>2?XXz(3_luxuh&_%Z;Y
z#0jF^Nt3_A2jckB{xk}8!Fw=EhGWT*3jvN5$Bo#{K1k0Hz$}?`BmGh(366OR?!w4$
zrRg`By2Zkg%$wBRW1ofzShmz1SQ9J-3om%h<&ikoRLMlg_62w0w75tr7|WV~Zgn8x
z^;m%=mLK^sO)|Mi@IRXS#_eH$pYJ%cqwoSHf3}9Ih9dDcCA+%_T0zHv0!y4M!WJ4B
zI05lPT76$(!+x|N3G0S_O_m)$^tXZ~fxFBaF?$Hszk3(yy3Z79@NO~(rdll;4wN=_
zSAxIVxFnCu?g}_htzcb+`i@)dreWn%T(fi+K1{yChI?xz5wI=_2hq(};Rb?B2oJyH
z5&6=r{d+X_i7sw1@K`I^j|#&d!HZoJ=9?grr!ejtg`Z&ONzg&aDpw3z8x3v$6!CZH
zIKtJX8zC>f(?urNjUPIO3ne&Vz(TUPX1W0PpRC#vp|G#nqOWk*RI%lVuO+%`<nLL!
z&#)mloiyi?wc&0l_8N!`d*zKRb8HO$HE$?o$%Y*Tejw)>(dH5r<(Ajh$8KD93xPFs
zs45PU$aRf4>m(G{^5-(iDRXo)Wsd>YBXB(z=bmG2H`J2wo}v$e-@$ZDkUZkq@xs>y
z=2dj^BJ<f>2<|HASRs@WH)rb!K}QTP)!CnEt)?05k&YV@2f89PKx3?djC=|r??Cm%
zC7Aw>I_@F5Ke*3I&gwhp`+Sl~NOO{jDO>L;Z91Y>I`R+C3*GmAB)<(z@=h5xc)GSf
zTKkgych~bDfg#HTM8{eodvydzM`)DWU;s!bEH6lvH!?F89<^@|CY+3MH6zP0d83H@
zG)9_Uo+k2hK%XmNOYp6|4!a&(U3&VK(Bd{e9R87b3shy+fN%_i(Mu8{4~cF?0Dr`U
zaf?;iwKIEEqrZ~4y(lK%CTH#O5a#B7hC{Ql_sn4e2mesVI$Ps4P)}95O+kglb=l36
zrlC|gmoGKojLxtrnWE6znwo`&np;UfrnI&SJ#y<mNDoHG8)gg|8kBymvI)zdK-N9A
z=75zzt4EDw)j2q=cRabx&R#x9il_-p^F?#Q!{h26UYCxDOlp{Ooq3;Ps9yXNsTz9|
z1t0QtQ7=*wD(B4k{CebqR%rTi%<2QKhxs_AJb^$K6Yfvz<m{7h!{O_Rl7;)_{{87m
zQ`#@-voxMjc1!fu3x3og&ffG<3%fQy=;AoDcU8rdf!39~<+Dou4uJ2UsFfJwNo*zE
zij@_WxrK>fFj^k#L>vdZ2Aw#Awix3O)06DtBw_D=AM6(A1-=<#UfPn6nd1d7-V75(
zzh}$b9eLZ`G2HW2G1cDI9?z0pU*M&!VVVGL5+sr658)7<(5^JW#*i}UwDtfubs%wZ
zCXW@kplD3YT8j}8T%YJ%>g46pxwYoo=~$2?G=3*sVR2-uwGy65D0o)zjt!=w;pFg5
zb>x-0o*Fk`OgaZbO06E<v=J*j$}Y+ec&fo^3?Ukp(xKuZEtw_cf_TA)GPgt+&aA(t
zEOu;8s%6O|{K^?X0nb0J^w+ohdNYh2tuL5Mx5wG}oP=P^rGMnbzWZPL{*Sj1y}X|P
zZ)}`*8i?Rn6z}_68&1W_O$y#-I}b#0aiHs$qi?rql?X~-PUvTuFmi)$6oZrJ$5}d@
zAbaVLV=NQjzT|q<Kk6<YR6UwmY^cq$Oum2UIgpK78LYiy^ak}FS<eTNwI|&KV62H|
zA;bSc{=aJcEs~DBAg%uoS$a-Cck$}pRNhAN%oYCyVb)H39m?V;`eLvt=kxb}e_g*9
zJGoey{~n^q+CkJS)M3K9I51b!`$?{t!l`o<@&^9$<bt(=2~$<dG|;&)Am2&bAZsZ}
zjf-h6tj=#8kX!cu)AQJ^>0F(k+c6@_=z0F~W)HagtkQsELz@KY*vmEJRoBcst_^t_
zvlXP^RmL_-r`E3zWF5B@>xKf@;lqWl{`AlUW($-&#L44$yji;gm-Y19^hNG{8?K^_
z0DD7QaAY|ywn|-8I8GHYlhjUuk6?tni@5#cTnEK0(?l)ta9J9Wie{liVSxb&c<}l3
zSH^io#t`Zz+f@ft(}Q(yc!^c<VGXh1hmg*n^wSZa&~}U~sl)+qaQ2?ZDwFj5Ob&01
zQ*Y8Go%k#8&Y8+LD)KVqmU)zB?@d^feEkoM1jKp*np5O)65IGa&d~GxaU<@{Jiz=Q
zhTm)*;!(7pAM6(EO}3jzx&OM{UyZ0wUxLSY$9nRZ2ss=Y#K>yxZMfa}E29%w{<_Rk
zT6WS+@d=(+I$WdJE|sP*$+N?RiD%m0I8h-cv*<EYfn>Uz8<v8D&*ND@k#>qt5*r|}
z(yKK0So|{EDfJ|A>$wt%V0nIX9(H}rtxHi-VBXF^fJhVwJejy=Vr~<rV$VFIB<Ey%
zOse*al-yUQk|Jc<T+-Bvo0vW_bT*}lMHNIo3$9@eD%H-DI9%x_u3u%<wG6c-meLQX
z_>i=_eRTZIJg(eNqbgw(uao{NG|??WK=o5z`>4Qx#hb@ls`JHzvT}Zf*CCE^{!gyp
zcAB|%ytVs8ZH-~8#GcIKQk5YT$smTI$w+}<VrIK16Ky+luL@cdN9~^##lYb@S}Rm>
z0<DHcT>3C6^cD%v%upy(GF5<W)ubt@@)M3omT$Xr3!I;B5Y2jM)#oK)DXHxCWD!C&
zq={iqjXMKLexd}0(v)KR-&s52%{7Xq9OqBf6)V!}5a|DI-m;MH<@VO`fadm?=SEQI
z0AzAX{l&Gv>{g@XvnNy=_FQ#Y&^l-@rF~unv9~Myf@F?|iIQq5OjBLGe1~S92aoK|
zWxLQ`iZ=+sipX>}c2(hBprfPLHDm3Wi0zyXZso{e(~e<3u)^w8O71BCF`$rZ=0-3>
zWlwBr{K>1Fh7q<(CvRrO_Hl>5f+r*7lQX%}1o;H#t2<l6pZGb$Z&IxK?<9u{0q+1_
z#7b}?<ZJF!0f<xmyOhi)&*)Eo5^39H`PFclZYtSM^?&|+7x+}vAHD?A*AP+}k$)Qa
zT_;N`A!Bey|9;s<z!>6`QfJ3k9HQhkttkl99%Afe&k`}MHflZN<iHnD#$iIydVje6
zE!!U+8JHCha-_17j1X&J)tTZoq~?EPFhZsO1oVO)pL|xkkGs8qWb}ZS?8Yu*cD}<k
zC4zxVjFG!-V0fkckP^sycn7~##HxDb<Wv@*UFhJRCEwo4-gw>*<xh^`$SiG6<+bWB
zk^amH?0fNfQ+RVmL)+{6)9aDLg2if@;SIRJThr72b@__fQ*R5Z+;SB%+Ei)>0xb=^
z@18aZ-l=}`_iFmK^9sxFJE#x6Q|M;+v>F(WQgu9sQ5Xu%Fs<hqyw~e~W0~@gaG4@w
zmrfnvkjDnOJ;Fe%xVZQ+e6^M|jdFHtZt+<g*9yM-E9Nlhs*eBTv`#$&#ab=t0x?$j
z(iQ4(G{7v`u;Mo}bv&br0CI(l28$k$f^Ga&#vZ_^VAjHzo}273h7b$O-6*ec#poLn
z{pxpfJmifN3yEgC9|s%toM5nMwQ4SN0L=AI8~?WT%Y8$#m5~uFR_56J-2*SxonR`l
zWoCe+gxP3nnt{4<nW5T_i=USxfpHNbfpTo0BGKGR@54kpj}*AL8)yAPSBUHeAg>YH
zAU}eSN{dy-*&;1#Y_e8bq|8k($t;l!q!A{PERPcouS?bbtg$V)!LFraAGP9qpXy&Z
zTMpZCv57JW*yxg>SibGgC<l3&Ny=;e(x#mtqaBTRer{t+^_RIVLYc250SoKD0ELbz
A2><{9

literal 0
HcmV?d00001

diff --git a/data/embed_chl_otb_ce.RData b/data/embed_chl_otb_ce.RData
new file mode 100644
index 0000000000000000000000000000000000000000..124d9dd59c62ef7c8ec060f4e6e5f3e534416d4e
GIT binary patch
literal 5722
zcmZvdXE+=TvxaSStA-^y3DJV+ZCRq%=ye4_c!_AyR|}$q=)H^HyH!{3$!Z(Duiks-
zyyyG(UFXcdnP>jYT+cnj_yZ5)KY?>#if^bhK-qy$3a~y2gmdWue5dfjeB9Qwoy;`r
zZ4<wvIlE%aDfJulD^SKDhYDkdT$FZ0g>k`;h)6w|z7YATyga<;!c@W_iq(zHEZVpi
z<sH-0RlY~JzBhj)pGH%CHy>S*Z*d~M$`A)8jaV{<q)C)?-Q6m8N`}~fj`t)|nBFuS
zhsYW}7dGBe`W_myX0BB%Xjn;rPZ^8V<u<WaQ9tp4_BV*CeuSFaS)<p(uJXqHD_!If
z)J$~L8i&gMP4{t1Nw^`o`<B|`g^eg;Abqd(!y_sE&gTYKa!1evdP9*BB${_e+5QW{
zbtR8YEt{rac(XndgzPi8f_*QufxEYitB34;m)x0J-xj{#Ua>9`Yc_Cpi*8Vz_B{;d
za``!+UW#6y?|FFm_|bLzeI0Zu`+!+75uK4x8U&v5%sWt;FG|OTB2Ht-8Pye|d@seW
zbx+5G#s?e~5FNy<Y<&%Bnl3NQMw%94VdVlHk{(*y4e!{XqIX&r%{J6L8b%9l2d9sQ
zbr3$rOtQ7M4}t4jUJdeM8i(4U{J38czXS6Lf4CLfRok!UlPqcmPVSU&Id%j&ewlL?
zZ_Oi~{xAnctPc&qr7T8zBV`b41FrimO(fTJ(`sP2ubl3y%KeVLeUz(jb-$}uKYqD_
z>a@nX1G^`^o_<M_IxOKOQ^vG$b8r^O$>8;;&EItjL~y+|h3xlQ)EyT+w)>%1T$ebr
zjC7%AJ$}mG9^K5)>aWIW^8;en?(~jJ*z@j&I@c)C_t9UWwSKBz*;0GXU;eh~?3qC1
zZlWEpD>IsC_559L9Bh)G>&8RVkB_Bawn$33YIVB>J33HR@f^+)PrDP-`rI^0T%Qc`
z#zjjQtQuSN-r&3Kz=w}!e?{|@caeaX#F=&ku<SGbjw%70<nL;eF9?a*G%pT+s-6me
z0Fj25NqFi#u&kip?@yed>2$Aj1}D?&g2uS>7dH8su+2&t|5_&~V)ML<2Kt2?tmcsz
zM9*0uDC`wnE?FOD$d(^YhCZw=u){ot*@zeyYdQ@I1<{Ll=og$6=eoWi1QqJgWj1ta
zP`=@Vj!KqqEyKnU&s4@w%GK_!pvH#HQ|G9;9>)<BGyF2>L}bDvsy8vZ<47F_!&pw=
zO|0YQXEdNwsqtVO15h=lWnF*1l}}id<zsI^+MaJG@8JtUtbp)?g^2e_ZCoo0veUlZ
zUfbh)lwYLCsdD!naO3rho7U5hBI{@Uk0r~)=U#Tqpvdea=zSIYK&u4*SSLT=Bb`uC
zby;apcVb5@hdt)Y6TN6tqwrQKS&#l~UWR&4EN}Li9Vw>=TU~+U^He?{PTt6FeMOa*
z_}ZG*Ue!&5TbGc4P2tz4ev)aYBk&EJaP<x_U54OQ@7*Y~VMuNK$il;|rbs5ccHBT8
zm~|D$n`5#rVN}jR6POj|XvVp#K=1}^^CweSi+Kfl2Uq+yAg0@@{(DvnL0kZ-7`_*>
zr{(FC{Bp3^{vx^e`ljmIosh>t?aUJ<MzQ>K^LLgto!;rfVq<!_9@(83$JJx1kF^SI
z!c2<N2Nip6UOE<KRVg%ZnH;71UI^(H9432TQxUazCLJR;%6t3nb{XA8`?^p7Vo)x`
zbFDP~0qMc~bETD}=mYicOxTE*7IVhzq|;Dz=z|C&B82yjmj_lcMF#o;73sDy<UL2D
z$9Wapj5${NUmxrv4+`lO>m)bKM&(gw<YF3Suf}Uc$}4(siMp|^2O25dZwQm+`XMxw
z*Bm^p-H@a?2@qtab&kPpHlPe;vroT_*Lj5txtK|NAhp>qYtO!Iq~Z+U>C|Ft&YL58
zD(@l!6b~P0UaYGB(A4irnwBk_f*)_C<6Cn)btXD~P8F213*~Z{Q-j?Pk{A%`+EX5)
zHpy4&-<(u58@gT`4zDZvb&GinqyLy{Yp!aQeBPTy75m{&|9xn;0u;60j+w~Nj@5;Z
zpam`uWJtu-$}2yMf71g-H)=fivGdSKqHhVJ`JsgiG;MY+EFF8Kl>UERz?Etwt>BI^
zA_GSR^YTm`U^(%Y$L*hpffCuhHmU1JOkMZ<_ge}_QN@j5pLxfVzTyax$A|vKI?YBU
zc!_C`95lvA@L=vc4J>-z|0FTy@cIic|C&MTk(095dVrpnNcAUV`7ZhavX=J)|ETDa
z_Qzz<9ISojFa1CF&-3#9r~|UjT~*YitazFH!j0R*Aw)u?BsbW?48^d$vAergw|>v@
zNR3v6OJ|z&V??`6$o1h*e0{fAiIX6dr6cZPNiRx@&kWahF4$fioVlZWPg{9i?Jc8r
z>Yf3fpZ;FuEN@)|vvwDr*E+CzLVpnv`K8h@<4aicd;?Bg?O|Bnee+{0>gR(EKTr1V
z#~MLF(qz;$A(H+pp8c)j{O@dE%;lZECC}ky%w_SE)@t)9jW*Ziq_$~jYu-fs=}+f@
zHqD0+I7-hMs5P0$tHGUqjP3ng&|Z2ePr3MfVRWv~ruGCscacPA-exV1nz7>PSk(l^
zi0dmdFb^_VQYVE^ZLaO2TnC<0E5mOO9!!(aQa#q$$MZg|SZlHpo9#F2<(b93O=qXs
zay?oXKEOGlik84bh-XZzWtZ01l5M&2SGApp3X^?Eol>Cc`fGJhkKc~){Z`B?C)A7b
zlb;A@pLITHO4{eyU%vDmf3qRB+P+j2Rm#Alo6Wbn{N^it*KO?QhT7hLxP*u#_k66Q
zWS_vgjyeca?Zl8_nyee2Q&!RNtu%7D7}4_VM&jd*mj}$TyCsKJMfE_GemKfbdwuvQ
zS<*CCm=Yq|&US%LkS<R({@T7EM$0&~?qFx!LI8U0d%hDR4bh$TGgyC=Oxj*TZx~TE
z?46-bTk?fmN!BumD=iXv^bB<xojc41;^+L}%{yzaNuag0YntV<I6T_I^;v%H5*?2<
zM{BS1xTVG(&f4HQ72^)Dao#g^pl<v{kj;vd{dvC0e(tCwaD1KO@o~J-RI4whR-fJ@
zIqw1TW?A#{()++C-2?m}89V|)#b14$4LcIvcg^B&U|6R0$@|=k%wCJ=Vo26tv0CY_
zZ_Z{kMFm{ZgS(d-N*Fm=lHK=bbL0}CKMn+{ymRjcfXUXM&J}K=t`u=_d)FsJK)CB*
zo7Ip08@qnVbL*!YTm8#^j?_Hrv9Tt<9!@e7@`~yhp*<p3hqLw58>+Agwpw)E-3@yc
zlUn?}6?G{^)=0e}<qrnzirY_hz0|0zZ7FtF142%*jJY%Axfk{BvXNIcYdvr$8fMrs
z+ek34Q5fsLU$W@InZr9D$p7<W5e0rhzGrb*LWxocWp&XrY)N;ds6b~&vNBqj)QIie
zdA_~9-ZJJ);lTBmuqr(XdH!!Rg&I*AV-}f_#2o);w;D0HfG3cZ@CtSg1Y=fiedm&?
zb;8rzY%Da8j>MBx34djr34#FwZYwkSW8G5>LOXv*SLJQ|LmZXujx|f~BsEAD@c!9@
zLyUIky0jnOwN1_u&W*sz4pSAz2g2H<ex%{VEj3$=B*qnv*z@C(9NZutOkey#w}2yz
z8}1cF86-w3t|!N!I1OYz;xDoha|Z{`J2)PI8$bwvk!_Di_hV*ZxG=7i<YeP|(I30S
zytf7S;8@|(Ah(HbnPU6`aKEtt3<5jo(62iYSY8-Tg#6^k={d`oH<A6w2Ec4;Kho=)
z6xO+O4g%W+uoBB;0O+7TH)V<kK=Ao-kKc`P$q0{Md-=pH+TMdl<Ify1q{z<2n0hdx
z0nIqa6*+8y_xQ(xG0g!zU@<e<bt-thUoQ9HG{)7laE$C{x6uELU5QP%ITRupz=-oe
zI9L{7FD;dK@<YGgZ=Cu(E9WkN4*Oa#<_LoUaNEi>jrk}xcm#M4ki<j_b&g2?9rt`a
z8b5j|n6*QEixK1aObRm}&xv*WC@>Ni3_$yJT4Y&$zyLFzw=!i<cqSzVB)JP_ZNM^;
zGG*S~S#`mGp8{a)Gcy13iwUDykqMB!vkrMie2@ad7O%JY%Nf?haQr)lq2BTDbbvpu
zC;tkEx&QNH*MO~nHXQdR<~G7(&g!yXC!fFd*H!!jo(CLX*8949bS?spfX{Ead_;qz
z0kD@dGEde>02J%`T_?z#5qRy=Y0+Ul5VWX|_a5t>wtAlA8W~XZEQ|OWyBb00DYfDu
zjR1NR^#(b|TFMk9TdJ;$Jz=zwuT%k+LK~4MIIy=Xc=?c-uC_m^_}Ad-bP`X|6}!%4
z3}P_z)~n7yUi4?&8_>)_$Z^oW3vD$6K!Bc<D+quea9<U}ALIH}OjP|MtDy+b?D#cn
z3R`ju#uCu!edisf#u+L3k`-$lLE@B}yC4Sl0G{5IN!UyEDU*y{Km^#$9P+V-v1TUP
za9{Dw<txfNvL6p&eSW>St-iezXD^b9s=*-}9B<xRcoB@~7-n^!CLmtCzmc?|?Xmpd
zjN-o=9__#zAQEh&q2N`P^~hTD46X5YJ>QkoBHzn%p$VL;T=E(_e&>&~5h0w+{vGeq
z^Oc2#3A`w>hg3YCg`dU<2?-@lNb!n-(MoThY!_QtO}Adv35y;d$c+`HvLjMi<w|jf
ztO*#I6U*}XA~#2_F4yi-C=$p4;YQSKc8gfp4vLm7Dw@<<ir@_9SHz-+qB!`CxdNFS
zT0hB3i$62dcQgIGfHkP)#|>JybPI=y(Ck&Dj|K~9n=k-PQ<H>}vHd>X4aY|Q)QnLn
z5H{ElO-+d`oXlN2Ra<X7Ik`Ng-(G)3!+v@4?fwJlpZau-IlGZQ`&`0`%@NS5UCG$F
zIrWpeQNA2=a!5RX=}?wT(yf#+!9k}0))vMf#)9&A02sT3R<<~fTDC=y)~v<TUH!Lw
zjiVh;ASPHsICVrr`o4MADE*U7LP(<S%Sm;a77mqEcuSU@Mt+bKnaK<J;GseD!UVkp
zxpD&QLd3U77SO1(^h3F@t+UD&f%kX>)8Gr6HGJ3w++x<D9=}XCKMygJ|E6&iFYFUT
z3df&tYsz}3<`k<mBL-#)WnDEo)l%h@K(+*hV1l4T&JY-LkP!!K&D?9t1iaixtv7mn
zWSZm2nL^P1hGp9vlihjWNUb0zta$Zv=GWf53PQunaV5TZXH3=%x}!neYR%lE7aY6^
z32*SM_dQqN^*sxZ^z^fGR!GbL)R!Wsk}x;@tZpKFD#Co1y9-@#pp8A511@>1Q2BSO
zAfXKOG6<(I33B-vt2No?>U94@*5XK=kC!_y@^Is6$;a!PByHFKF<Txu<o_#If;Tq)
zD1f&}?|Od9*5`dp5qmI4KWwV>^18FoWE;e!$TddCapkX>V{=ddW%CQS-gYiumR8>Q
zBdLG4%o!Qn`S+og?hVew4QdaoiT9{KdGPORd|w=@mY{ZL>p8TV$sAqwVUw-SviNWF
z8)ltSZE%#^r2H2z{}OU<>NRYJ5>RVqh<YL?k;<+LnWK@;ussfVXzqXxZ$fHB+${(0
z@!d#C3#~mH#KDe^(9jvkFJh!mrdb6sa+X#?*KNtfRA|(v?suR9$F^_B*O#xC9Ca2&
zT{HdRH@hKE-;y(IGBUNTy`{bh`}3InQtBU~Nf?C1e><A1S0>8b=f)AGEnu0fpwbPc
zS=E2ptXr9b=}nDu{1u=wNJn5p(WA2UUF4^Ef-njA9Lcs6LF}cq+ppo*1JL@1pXt2C
zIiv&=@TvBD(6it8X%gG!U&Vq4ZAEvUzARzb3ky0tQ{83xHmkj@!Tutld$fkYf#X(?
zY?G>cb2NP&vi!)m%Is1ydpRAFQjdtO5faKZzfMh?e`~?3B5fjHA}60h?Q(qhG45ba
zt{zO7B7yaRRQ9Rqv9K<<#bCuNXB>SMmt4Eeqsa#kD2QqX>Jgqlod6AHbWHb4cof%o
zT}ywoOUpwd2o;smHJjZv!}48NOLVC%Pml3(yMu#L=&dz$ucDrFDwODP9-IouN@>6O
z&V*}S47L;JAV<bTCk>-M2DxvfB!qF;-C=0wtIHQ?{?N^zNfg)Sq$pv>>HbP$`Q_eh
zVOv;rvvp2|Rf{iZshqE__4fTJE+Vq|h9QblLG|2Hv|lqE#eqzBiJJWM9TAq$z7``A
zF})kzqNn(t)cTASC~dDkK&Q+|Cj7Fd!bXi6;ya8sM|2MT+5wJ#Hn-V*2Rh(*SRl^2
zII^-U)IL)3_~ysf6k53hu#m*_zHI4X-^+V&W8L=nN@D7UpC^&qx)ZgHjTFW~Hj-8e
zZ9Axal^8pa@fdrAGDEo7RZY6@3qEIv=G|<V0RzF9nUmj`J+?J5ReMPH-WEU$S#_;|
ztKh^8J8Nn6!7$)i$HV}oDy5RHN|-vY?K<9O_V{=CV*Q~YR4u8OU@$)rXMm+J@e;)O
z(w{A^K#e=`9i5^8+yYvub{)rOQpU1s_4a~bdNJO{ZVgR2HO<{tl|F-V*w<`QmR&4V
zW3*k6ZW&h;+;g&_<BRFG58IUQdTZ3f!$)L)_>)E{mSz&s!)KaCxb`GNr#0PNNMEjQ
zsYK2l`i-tX-=?ynVbxK(ZFO_4VXSj*`&AlY=^C-#>f;~eF7l~eJD%n~MU^Eg{rCLH
zi2CwA-|eR+w(~0EMM!^^O<d(Fn(s!=ygu(NCDBan4^?SFsS(aX7jbpXcZxi{V9hs4
zgDLX!E*)~iW90S);_nQMJ}ub~QHmwa!uVJlyqn|F+J7IAaIjKWSm*NRPUWakZjM-9
zTHjAtul!@?U>DLn?NA1(7N0ZVqn77P31cSkO3<z{=1x%7?F8Z5!hkac)c$0J@M^<i
zO*~7rcW>IRm$<9bG;&LwN|HYFFzF5ZX)7y1dgMk=fu$>>OWg^9jSh{Qp{iE*ZwrQ`
z`65Rb3g%8@YEV^>spMz2^mab7i>?HEI-GL1L^vfyyx;;Q=5X`0CH<0!*tlWMI4jZ-
zP3^%k0~2l(S;RXxk&xY1m3mEhz|e#Odm@fbgMMSRw5Vp(i#M-d)zhUUgn!%~PCH=v
zF!NHV<m~}d+2cKc+dY;W#T_}z<q=NraC5h4@k)CzQ6*8c<rtfM+rL^da@S#sV+qFC
zVF1khNCpC3IG%skNvX}2LTpO$MNVks3r4-t<}^_{dTN-2Za4g@>e@vaoG(0u4WF(0
z_%w-DiKp=ReoeC&PTn_|?0BnFNcFF{u{%wZ#ttpM3pm-jNVG^1pm@ir8><4k%bL)l
z4LX2Vbh$$v?W&llloE_5Hi`ncW{AWX>x_!t5Q>Pj6Rm09kD6CWbRPG>ca~D#4$i*T
zRMvQ@m+Qi3(Y7H}P7RG^pvy~A>LdLE554>4p`-ChAy8;$8xh#C^6JCCCsBU{yJ)!J
zIoTPgKB~i#<u*D~I6om@97ZRAw_o3fpYxHUe|Q6!>P%XgFF5%{#Q&dmhcjLNi&p-g
tl7f@v=eS=UcSok8Tg2#@3eiwK3d|6JU+F5T)WJ{B6rm<MO;H#a{{u(gT8#hz

literal 0
HcmV?d00001

diff --git a/data/embed_chl_otb_cw.RData b/data/embed_chl_otb_cw.RData
new file mode 100644
index 0000000000000000000000000000000000000000..1b3c354f727f5fa9a91a268419983f39c758461d
GIT binary patch
literal 5695
zcmZvbWmFW7w#8|rrT>Eh0!m6MFj9kbHw>XLgdox&-O?>FG|~c6(mBk~5<`e|4c#!3
zBaFkl_pSHgzPH}_cJ@9W&Tp-~*%Kb%{r8ZZyd^cIUtD6-_>9N-$;I2)$ho>PtR~8s
zM&dd6`^_f~z>|OI)INnluY4PB8Z)oHRoj>eTn9W)SRhuasV;$pM{*i4M1^u+XlD<g
z4q=vsU2DJn&$};R74!a+`CHa{DYPywDhm&ZtY4d&XVWqf9Hu^`WZj(kO8D9|*DxAj
zt^R=T;}b<fwKzD?bhloQ6qCS3*0&P3S*`b8*WYW%-)pXnRb6wq<I9DzZp^P?v!m1F
zgU&HFxVwL|pH06U&}(7NO($};e@`|gzQ^J$`b#IT*isgm-tSrvIQh;bc5TH-h{n81
z`H!N#j3vXv$Yi+|?~}V=f(aQ}_)hzdw!2P<0%qL8)%3{Fz!$P>oA2%DSZXxJ3aswj
z9$_2y2P#6me2gr>kHKux?>k2+VN(+z&jqI!{A(aPu)Ujuk=WpEUX+hRs!z-Ny@ue8
z_sgc!<<_TJ9z(uWyzVB#HPfzcKVRDO%`KTHlI!vgE;=Mgjt2)bJ9I#uSu;AAK_<;)
zt6p_%<|TD<(~Y3;#@)r{+$nEje%ICVKXGnedcUK9)>?{>%J|HWj#Y<BS!XhX*s+Z?
z6W@Cz>`z6>=5EVDz*aFlCy!fSAY6RT?~R|}aR?QhH`%n}u^3&&Bu~iiHHI*s)1`CB
z-=r<xz^-8H1oJ!n2APH^PKx5q{+8XuaRGUjd8<bS00}DObCy@|gfqGPSAHJYvS*H_
zd2RdXr9MW}F3I2TgT~2f96zF3++$fB)Pp)FKA;&Yn9Xse6OzC8!(=K0+U7L%bYL#l
z!WT5fY1kvavVfQT04U67qWFr}DQQbjdCfft^=WmKn8JL^p>mC6c2pJ^qt<H?X;4$!
zTJOt>q^|)mz)FU@a6pUhPKX9bFl41;0xPS@UNmzgB$=C*8FgMJ8C)-IqFC>S_43nD
z4K*@wFkOirSljn^_HcxEX#}<Ez!IorR5lzIpZ&SG?;v!d;Yx8ixu<b*P0E%C(WI$a
zq?`KNqDhz^EdRyO)0<ESM0X&t9zax3*oGT%jS=N#e>Q3EtCO_$QO<GtaB{~8l*@X0
zAj2wK{i(t``}r_)nQyM8$`<<KNPnV1L$?ODQ|Y3=86)n%CWvk<-h|KVNIn0K=k$%)
z4c+SUX{JIV<v#ch@^SE$uS@md&GcA!twH{AU6H{Ur-SzPQ=_2r`Q6~_cf2p`-1Y@3
zynu0NK8v(F>5$*5x8p6A+6~vo?20tL9ZBLBBsX2Iv<fZ06vGVFsiqeqE;VE9%pSKH
z`ZE=h@VGRKGi?dsf(Y_mVVt0ID;GZuB;$FV0jTr?Cty67nx@iO=c!Sz4lqO8tp(uI
z#D25J$PM$XJ<~3!X0J1!Uel(MFXiRStTn+W3wKEujy1;47N%QMcG&&b)ygLQq$PPK
z+WE=89TAp-Eo|#wEocLdV17L^!9x5hE_?GV8;NYoJ)QKcZzA06EhgD*(_@>2TGw#g
zncwaKtJOSr?7g{3@1!Mtq$0e!(`}sd7_3ipgQ?rv$of*SbjkIMf1_Rf$$fWNUqz4~
zk+;3O?2cYjp`Bzpm5PZM%mrVlv&XYB!lmc~OGJCFxPK3`{GM_5skbrS5MFJNXbkKx
z_a(8>HZpqC4HnFwxgV&vI90Yj88RzdL_jgt8`l@<eQ&1<aG45{R>DSp*KK2>UMa<G
zaXGC_G^Kje;n#MfkKHc6k6a-an_4igS&5#L%yPiM^Y*lVhFMIH^c?iM><t@+Dq7Oz
zu}pgde0Rg^H^WOjV?WJ%er}QqB>VZg_&!{<ni286nw_cSZS_9<m2@?|JS^%np)6YY
zV?(|#>?&d!c2yk_SkoDEmcr#0+H?V&u-tCCIu2JiO<!t1Y}Z`=ZolLw>2SvaY#T1t
zo=EEfR!rCP@sP9)-cSfbxo`X;*v!BJ(ooA!@q_PRyOGL8zH0@!G~N3N<hoU%WNY+s
zuJx(j{>+<E?48{#qCAekx#q-xtYsv*Q(8c&+}IaeWDPMf9BQ$-!nKj0V6wZFzhEoJ
z6VTFoV?d)T=ldq|=h<rR=ac#T+{xHY@8qqit^pzB{qy}ua?+t)zm%mg1Fjv%4pL#%
zCBAT5v=)1mxKoz1rw&zLt_O&mjhD?f>0W-4H5r7GxhbU42E8KR^J&cg4o7kY-PjFv
zuDuA6*R|A5y)8t)uBTb_myz~kg#P)+-1BXg1t#Twqm65l`3Nh*HKn^{ZVk8RC|)DP
zO`AqrEczOmfp@2oxAbk&w#f0<TbhrO`u)eG*g)I|VP3|Poa@Qj18os~s^$VK`nzNH
z8y=6urt*51Kk-QVD7<S$rpLE26n#=?>lJN@+W&b`8L1gcw8?{96Uv@<e>H3M!N*Eg
z;;{oLW8yK;Szykt@WdCi&bC)D`UvKJLXp3`?VohiLaTGcGJ{Z*qL;l$@Rn=3V+!yR
z$D+c5ZS~A*qtE5ELsW=0<JbA;Gm0+?EDNB-wZi~oUFvg~KbupUL)j@(t$%SzHdXd7
zn~h5ZDmTug6<CkH2_E5e%Rid$P_Z>gHvBzntGY8hebIKcH{PBdr6F0L|J|aU9ZHoF
zy7v@X{Hb#?Qj?`#u_92)Wnp=rWMYc@9gSH4zEiS*1vkcE3o_j(j?6T(b}6*IdVBq)
z?@#cZWIMNGq8D>!WUlXWn44*#N%s6&bR&~W4aj4O9p+Nd5Yl^|RXJ`u3F@wpkr#f)
zF{8Q;E>ct-g<&UY5IHAOlWI*BYZm=dz=?u$B8Qv2_RwK*{*dzaRf9PCiJRcf_2T`9
z6Y5TXAggF#esiyj@;3+Lj8##k?fDgAsJ-9U%r2zznSG+5=P7gKdAe_VU0?bX&Jw9H
zEbKm%^Jyj0-ge8~rz)jW!$4-@=L$=1Jupr@Y9c<9OTj$^Y8$uZQ0+dEYik#96$R|5
zXg<zTqaD2!_a33n(yYA1>}Dwb4%uJTJ^Li-H&~3FF%-|ka!#KXzkOb|XR<L~Wq7MY
z|IOpN;IM5Z+6*{$7r*bXB3MKA<sS(rGtGP}<X7XEAEZV&xAFX!YqUJxtNX>K1u-wh
zkuP#H()Ukd9r8R2)Gq7X*2No#z4g{Q!fFr#{_iXMGv<X`n|pIe1MP`hsfNtbEovwq
zr^j18u<f51<<UD4olL{*5$41ByxaFqALf^OsNJ;B0IVL*@o1Z@;5%9+Y~{-&m$p&;
z8c`%*D*gFrPCL-pD#7z{UG!uf_>mhDd6kU3r*+766pEtG++-dN=<(PV<AT|7q(&R)
zL2q+r&eitp_RGI%JWord6BgwIr)$K2uV972^DIJm#6^)mnppkj!q?nNp?OlOf4iXp
zLZxZH4n)o(5&xJE`txW#{U@L7+Vi_4{t#g(i3k9r3fUvbe1Y{3MBGw__E0T=2(F1k
zo)$4+coAlZdjg-wsMn}_JS}1z*}@-6?p)EJ5<-X~>Med3-WkR9Kk@){I{qF3oaERg
zEG8U8lT8xDjS`54Lq-;u^YI=Lo>4C_V(=k%5rU+F{3rrcC^425vjwpt=qAK^A?Tvv
zrcn*Tb41SM-8>Lq6e&Il9y_%UCCUahPk2g#V^{#=o#7V}W1=^N$?%Uvu;omUHquM{
z5Z)(b-7Lp|FwOAw2kp#$lzyxjQiv1&x+wOYG$Md9w1aVhwfB_p`T^nziico`08Vg>
z3?mH#k!BNMZNFdOjRA1+vJk>(IKkGo1`p~kT$V+I2=gBC3Q<oWO^tqtCJU=V^aDE#
z5NI$GaQDF*p@Qhb=S1B_$TI1XVC@k8|8_fBF{KbIyl!HwA%Zg+&J8(s5=Xcaw30$3
zQKSUPq;izU2VoT9R;29|82W!RQ(!6n%}l+u9pDlTU!y!`#|tC?5$REOlOGr2F%!bs
z+-cAvI(1gWINk5!Vda#CPto?@`TKchfk=p~d}ZO&{Wf(rMBr;+Wqp|(f`%yMDOxX7
zy%(V<%dPy0>MDnU0>|_(*v*Jr?Ii-N4vGdZl%TJW5GRJ_3-jf8Bh?PalVk6LiUf)E
zqDST^3*S&-mb5V2(FjSx9;zqJU{444GyEXX1JaW$85TcLG<_IY0i%wX#A6$>IJbGw
zM!q762WRX%!N-Z9E5blDSmOtw^cYb1B%TGK4?ZR$B-N&OAj_s7&6kXrC+?v|3+ABt
zP7%I%Z1~ro+(qgPqwzu1eZ&Op<UvGe!*J4q&tmL;Qxtqo_@NQ*(g;~LOzFc@N;Gr0
zE&)ygO~ZkGf$|C$CU|b!w;JZZrD7Zn$KTo>YD0A4t&?L~5ig_$)FjCU#8TtgMyMbE
zY>ZLUhDW`hR)P;6Fcl3NzQu7d-P0!A#}0-LY=Uflvo~N;OtGpp@H^G#u)2*aih3)*
zzkNr{mT>fK)Vy$a$e+CZtWk*#)3^QqA^y+Zx$m?vyML_O6XvM2<6S6mmp|jI!YbXD
z+R8l>7pfW`!5lTB9hE#w3A$ACZI;Til9w7Pl}1XPjp_$F&!?q7?CiYL&5;!YI@@V~
zzsL2>om825QO8Qtjvw#@6e8(lU?eXcUa)f;R=%6A**Q2o-95C<&~Fln)q?y}<^9U8
znY;gh{56jfJ6RPzSx&D?A%I6gntR&Xz?ddASBGhn0~k_6$hB?LV#G23uC)&~KGiBH
zxiP^BWv6)i#fJZ!0F>^sBd-$0rw)zRFrHNv;wa7bdk;4m3OSKkS(RCC<sw{N4pe&r
zqiJut`l|Sd-^v8KU93^00AS&u2e2!7o$wX%8uM=#m*v+r`hG6!N&5O)*3@HUTD~!|
z!;L|tJgieL@_a$AETQSRv<ckO%_(rb!}zL_WRiO4YuW{OkF)38&<Nlu$7l7OPL3Z&
zgI;s>b>;k=Ri!BdBmuP6zg{eGJd4fH_|<!i2B*=hKMM@f1ATiXuCu(mCX~1#`}4EI
zU;n^Ef7i7rqk#|NRn7vUqX|^Yge424KCyY)HayLGQIU4lX0+UNnzJ?Rnw=U8$vkOM
z#0gmA59KxlV`N4@)s3JHD%oliESZP<BE0*<$t+&iyz+qI)?f31G!f*Y7uPQurCAn*
zY8*fJJhQmPyM8?+SGh=99v5Vh{lTyrV!*;M9^FLF$QM!Ox%v0YySZ6j3H(QPB}U3H
zriC1x$HkWw>vYcMM(O;E<Sl7&C5*JU%|`JL%Cil7FE57+%Te>^So1GdV;4^TA%c_j
zJAYcR^DUV?m?nnqkpCUV9k;Jm|DWviuL#O%{$3+la*&qOkGi(rY>w=8?dUK%-$+vJ
zzo0VUn!HXWFW1_vf2vHUA2yCzAvvB5J~Y9>Tm6x4$P0^$`_$S$@%PaB<$yWX(+%IW
zKTt||+6Gt*s}XHAf))UTJiG9&wR;f4lk@LB{O4C#<Rt|>ef<9xqlfA5`Fd-*mYXpP
zPM;PG2FHJ!C>}3?3?Gl3{{3l<D?mptuT5JOQLaULu*ua@p2bK$?=KV9kG}D(0{v@^
zl3WR!J-I%Nw+^^jW^jBM5t5P}H}Jk~WkEvgO<Um`kCy&W2zYf=<P~%YBSTD&g?*%w
z0jP|Pq#p1`!CbfacU(uohkE@_xI0{XiXQz6V~%AMi%Row-;gprcx?0;Xg7g~k#BQx
z+wXniDzGpe3OeUc9E-7$lwm-eJoHy#Rd9R6lr1L@;dBI=>d_WgMbMID>BdrZrsM8+
zyxQZ#hKwQVg>OW1C|2W0bcD1X_ZN9iH|C4T(waQA+_?XHn?g;c_d{zAJ~3E{qav`F
zAtdiC^l;}mx-NcD$vIl~p-(9*!k77x%U{oub63QgJ{q6>wu>lHl`Pn;_f9Ei`ix6B
zc>3*BUk7S%;#ZX5<2wQPR1}DKsF1+x3C>9;{5Ny+ejsgI%mT4Mr4C0XpHZ2H%MU(<
zekBbqDx1cvSJTLO?XqL0ViUue$%o3|IHLS-Ri7ABBh)+WVoSQTXnSiNUX(s=Fsqe~
zWiV{0gwooG^LcztDT*KXU>^P**-s*gO_%N*Ogs3>R9XI^btaE;ysi#cfI5-mZa-j3
zqElAC3I^snY5fG8k%~(u=AGbd3~+I8*AQ0FNpGhV%O|6R<*y-&Vsc$#s_|4xH7B2O
z!K61g7dL>z`GO92?FRxg<{#%Ya`tl3HK9B+T+k3z39f8kPnY+(L#VlIh>F{b{demI
z@LQsbWY-!2q1RKP6p3Fl-lEKY__yMRwnCNxi^jdJq+KjmUeT9CpR*1l4>Z@ux#<~&
z@4=||gAGW2viGMugJdJ#4b=L4+Zjci0v>G!0{iA9b83mdY}Ua6>vWEvjcN(j$0O%r
zAoqahrb@IQqCz8O<Y|fPjw7w@w+-cC3WM8LevKKbKkQJ$y+o5&PD0GVWx5tEFEcJH
ztTs9WF^LgCan{AJj_>on((QeP!E$c{#7{u&1=hp|0qpwL{RJPDq;rLz@N*NpHm`7t
z`c%J9CWuWvoe3=#b%J_fL%EYGo!)GbeBInPohV{Pr8ozDJ<RnyBf?4LorBeh4cEFd
zyyAaS#@pbSLA_b`!3I#zkmG02=-J)V-8mzxlj++JZ9JFt=BpryZm53uO0Vr@90H=!
zX5$b-2H=hVTovQhtZ0l1@Oq{M6-ia%gf^wOzAG{s8l`tKR?4pXe%xL$z_gh6HmKnh
zW6&EN?KP$wGf#CDeSRaMG>>yTYbg-7)ew6DsleeZQ}N5I&f0T8f+8$_M9p9%fKrd2
z(>TCT<h77-gvryQb3q8y9;2K+*GHFXz!xv!+2HQwU9U9zi^?Qsl(zw(^v!sEtPFqM
zi%Yahct5epiNcoH9LCw~RB3d|QX+Q9`)2a^Z7tLn@aQjAF_bTCsW?mWQudqPOQ=DO
z736ckADePIfpTl-+2zV#(_#jUT4~R10b22<e?BOxv!((~ed()uss`wlUOT5Lq-yQz
zW|<q05a}n?DTR^&R2(4&Ldo{zPK~KG3VmE^<sO7&OPx49NBYGnWFN`hHh0x%B=%_c
zVLX%^sv(%BjX~hul3NTEt&k*bo*r%-ny8nWDygX8C~7wtN^UI+EtT`!VZRuv6aU5<
zf795!vvO%RA77b{ei7J03%_ADua$YTx+=g5_D-yj<DPgptJ+%aVUhh~hSdBt)bqM2
zil3bEKn@br(Yo0=8>aO`>4QP}?z*9rW4~GQvU}#%ZPB%cnA*1zrK}XXI={USay|_@
zF(@Hdb9e02ZDP~h-|#!@6+ag(cLOGB<J6uqDN_y(%B=2{v|h$z96qf<4$(@7Q~M+=
z*HdhoqJ|?8_h_F~Zga7!ES*;nD9@1|zgS=7a~qv{2=XCU2dGQ>vA5D2l{dbKy)Hkl
z)zgC9y^f7<P0U-pEXPPNDvtQALr-*|SNEeY-3P%|twEQ_WBG<xl3>+FBi`k%gX~2e
P%)MUt%$;g79^QWeP-aUK

literal 0
HcmV?d00001

diff --git a/data/embed_chl_otb_ne.RData b/data/embed_chl_otb_ne.RData
new file mode 100644
index 0000000000000000000000000000000000000000..4f8b6fd897c9a66583f1e615043e60cb3b464874
GIT binary patch
literal 5723
zcmZvbWmppcz^xe|-Q7cw?(P^VNP{9U7>%IPT@oVFprqtzP)cbSF(joy8a6;WMl(iz
zaev(B-tRmA&UwzS^PCr$K!Ek%!aK6Ye{D8M*OhMF9$}lHl|;$OA+XRd&U>|B3~WaJ
z3`vV|q7R@4I3#dO&)VudYtVE^h!^7iO@&`fKu2RFxUd@)pZ?@Go(F2pC#&PY{sPt+
zI{fGI_DbRUPQnGVa<fz5l`XnZm)le*1)^~N*p@?RH!!0})u}@r+kVYW;A5I(ua-^y
zifK>v6JDzQ2jkRPNG;JxRvNC|erd0e<j*v2ODk~L`Rkvk93t>vsjYtv22WnY?17J7
zf9zhf9-?6$=ZvlIknS@*XC$G&#D9z1jm-1%^i<EJ$ebHHedcWNSs^vMDLF`_6v#8G
z4c<)khgo|$cORl=UO<pv@~9B`!_i)Q<*p@y9fv*S(xJEf!82NUVU>?HfPJITnS&L2
zvMDdL(dFVaV?Mu8_!3-4Vth<%LU|cd@B5*#MQ6IWZ!d%{kUH)-DBnD`c;x9VI$%Oy
zo8ADySdhLK?<3yam6s{Wrh>6U%v1BM!nuwtlN)GvkY6fVYJ+%w^cuPwqS=Q%<b*NP
z3eNNmY&?I$J{`;YvTZe42lhCBW_sD$CQTO(17@r$Q7g3hO1F14><__H+Bfout!9aj
zDPObt!j?+j2-GTz^ERM<z5@rX9>LOWR=sZWV}zvu7U+id81^<z#V(YHTsw`2s_Xfg
zkonO$XrKEJi!x(&qFYU##E0Qg?WT}lIGYNMni(;s2C_20Ko&<#EoRDu1B@Z#1!37N
z(LS@vtV<u$H}lZb78Zu(>xX;%oo(*)T`m6Z7p0Jb^^1r6gR@bElEgSji{hy|CI2nA
zZCv*uN415Yl5--8$xP|=jPR_1XsaKsz~bs2=OHw@Y>Q}OtwmlIBDON=YuOg%{kyJB
z!7(4KeSB;oU_S5c<)CX;uvaanxXqD2JTP;iE5w*Nw7q|=sZ98*QG>>%@#T8MuIut2
zCrWz{v-@g78$~G++S{?4Gp*MOp|;PD2BECho)6fUJJilAz5P)+(_HsqUcDKB!7RMI
z-))7fg5!p+i72@pn|j7on#bn4$1f=pW7=}%LHTbv$&0J$oxdN-JjA<I7|Qwlb)gSV
zo2|=5WQ_0X6h7Gov{jsr|Ms8a>h|v=hez{IxOBqHzCQ^_URdiXn=nlJm9(`vy&{2-
zE4-u1N|s+7wmS?3XYVyS_+K(~whRZ{-6GnHHiDuZIwX&htDV~+?`9xzvo&hKg&IH2
z76n9o=-*Xp1dt~I^Iiy;C6l=V%zwZ&H(R_N^<hr_!=?c?%uInB4m2ko$zRHaGspxE
zU1n6A<j+>u^uoRsjpC%qj28%LWmJ(-f4UgfYWkki-QKAC$tIZe(!&2jl%F)O04cgm
z+}}2s>+zlMx14_pw>fEi-;-sV7?ovT!krc!S+9p4={!VP*mBI^^E`92^SXLq$mYRS
zezWY-v*LqhfZQU2yO_vt#GF&QP;EjxWyA+8qWTvllGAiNE}E|7KUMXF+!7gbKo;tW
zI9CThj8zrR<g<8fjRFIXOvnfJ`;-P2Myc)mDqNzCR}5R>obIPzhE2ff1+Q&Y4j3%y
zLL2XU#MCYvr6w+rUv4Y83lMb<ruvVySNR?OyD!Xa=ZBlB5p<!wT*zJwr3D)FeSw+T
zUBADdV}>Wb$;T5gaxrwTlk6jKXz!RS?sawM7uzm*)u(eXSNmaljeINIlzRy%%c!KZ
zGHyTm30iq<?{NLERq%qt%W?4iX7v<j8z$>ijpxoPg%@=Ps%ovYq?!%1X9+(O4D2BS
z@e04+ZqZ!xCK6_L{Asy+E+=f&5hJ^+LQ79+y3FDsPF{Dt+%4{~W_!VVuOs4m>#N`K
zZVA|Y$yNPaM{)han4(?GBY1XTd{6dvAz8~EMvnYRiytd=&64mZ<ummxr{TtQc1~X&
zu?&)`1!w#1%GKtExxk;dDCc!6+2-Yw7(=61Q;QevYwAWX%T~8(Aad2PDaXkp8<i6J
z++wSpL=kPSsF!0tAz{(OH|iypNo`R+`S}W3A86u;Ls;8j;|^@XyBQGMjcHQ)Pu9z`
zTNe`G?62}OlK<GdR^fL;$xhh5ep+WG!yzU{n((oi!7LKcM0NM_b>41VbsqNW-9)~I
zMaP|u^id%0PAcNbkHap8EW`Voit9gf>+#<~HiX^2JzMl-kL^gvY3tp__r~NTCC{7B
zttZcNRD@Z45kg%qGQQ0@E@QiA(Urx0H4pd#?j<`eCK?u-{?Ps33#pk95A8OI1b;JY
zLpUQPzgOKg_!MXBrWryb)vw8l3G=|5g2BkUgFEBZy7j0s=KM^9m56idZt&ftJCA|F
zd9I-gxl6jjyOiZOh}UKgEdiDTC2D3nYnF%7;c~&;iTD0ynJil>7^m-wcH-mcDVKj(
zmrw8NZ^|cT41D2S5(YvO&y6>zbebM>$d%zvK+Wst*7Dd{?^T>O<XuRF4+Ehyu*Fxr
z>!1n&&Ufw5pG&cwF;@na<p)g>=-=eVn$GA<{RSRm45lI>SE^|+=wLi}U~Fi>DHhf0
z-}!+yt1OWfz3x2L`M$XrvnNdq{&KYO&Z|-^T(Mu!&qvKC)WspR(DYhx3PwkU-?bWu
zs{~WNh;?|O4~vtAhU2TFp!E4Jx0>ThDe>0LEanE%y$hyca^2>V<FHQX4zxGGt81${
zn){(&R_47qzsIc^Sk~^_v!E8BuqCQB%%+h7VX-k~f0kc2nRMA!R)gAxjue_Hlucn@
zz`w(MCrXDay?>D!zp{v^_0tWc3I0c%hd9U+FWTl7R(VJ&d5>fQQ;KyR#N%YNDh@vo
z)2LvZhoUkSkyW_ng2_T|@(o<!shdUtyEnSN+d+ZE0x4aU8|ThfK6<8;=TXj`a<@vx
zmy2So941jl?%;p~L%6;tQyhM2BU2RdcISfs<G4SVVo&+3ZPPJNp?YHbBX^vvi2a&w
z)~<aX;x;<_g%h*azzMmzMd}T&3s`^Yq(7&!Jcv2qc5TAQ+IVSuEvuzuhpc7e3a?A!
zW--jfYOF9pA{zhXid$+hrAv2Yc|njo%h9^pm)Cs#%jY>~2^-9YxKvx*<KI_cN}xB6
zw!vTu6}CUi<wPGbUE;&E=})E%{ez!eD-Afo&Q06XJy&9tXI~4J$(LKSW$Hq|l_`f2
zAkKHHegyO*wO2iqAVy7CDq!?jn|eyXW-%(-!#^Go4dmLf1<k60z1vfAj?NNfc{I<4
zseGP}fVXa{PwR0i&)yvzMx^8xWxcnTzN7MHQffjimW?)JZ;cNVx!(Q~>AG{!K(M1~
zT9L3Z!xU7#IP!(Fy-uapf@^p=jnXvUo_2$2{6`a$LoN;{eiJBw$vAQi*t{T<V&1yI
z{!>U}j>sCkWBr27KT)wfsP?a8)t23wPGkZv-A=w_Y|++Ekq6d&NUAo#5I2V`P<pqF
z{ou_Vm@^eMgO&@Jyp3WZMAM{F^>BiOaY8wFJ3tS(Xp7W&tO*=gUgj{jWWgcUPf{dy
zDh?Kq6yuv}(xdXEjq1b$gdi4YJMm$!!?D5o`DMU#)K#gU&{3gB2}Fh)%Ch@V^hdRj
z+*=Bn3p&QXd$vojx`y3J;7FQO*iv+PuCLmTwZ>XXj-*z_!%}9x&}4VOLg1aTBP&4N
zxM$Kxa*%SxDe)`z+!=ICa?v%m9L<RiyH3x`A`!ecvXdira-3w44d9M#mx~=Th1Rr6
z6>_~RX5I}`1p?ZELlzMefHs<wI#mSW83WQq)fuaRIGAxro~e75!0|0pH}3gA6O|tA
zNM%w-QKT`uITjl+#ywR9rvUei2`OqQT;;0xH+{~O6%>xuivJuJ{d&?RA|3mIZ|F^r
zT|wys7TRmlRCNujlfqFd39>zfadMed!G#G7fg)1y<)}|YC+VZWIB2?5uZTkId+Z^>
zo&?a`B$_sx$NbwAI)kj}8cPoMB#vD(LICtg;z)vAR%OFx*&fa?>p?|5F#IECGTaQB
znNA*<$)kR=NcF7IkSTGwu{I(klC{x%iJ)f?NBN28o9;&T9w_S!hidGVo)yvQ9M%JE
zDgIGzGw}&d1Q6$rW}y&p$AM&I|J5Um4^IGy;CC=F8S=bQT>UaIa|H0vTCQjTbmAb7
z*`eLjx9Vuy|J;N10)mn3(Ts@2s>Ql(ix9zM*|rKMS?E_o2SlI%9qJK-IN?kuJ`ojB
zVcf_7zy$ui*sj}W{7sI&L<5n&ZvzrcI1g+{Al9+&vP&TL8nY*^H^V}2=Xf>=n#Y!&
zAdnZi9x;K>5~w!7-xDd~_zt@b=kN-IHttD@w&p<Y01R>NxRKYr4*+CYgn9IjvZGpf
z$YwQb5+o+byNMQg1PH+CWJXGoiy`$D)66PMoqysGO9BPWjK5Q8eye<$w-1O{l-g0Y
z_X0Mwro3L!7Ch#pyQWOI{oEh+8`BZ?d-{|ixx!O1W|NP_a6OP{&{VbnewyCRV0feC
z!rbTZUykPgEj+j>Zn2U-Idy~^eeu%D5gV@@O6E5Y+<sfI+J4{)v)<DP^usd{eA{5e
zMT5&Psl(*xP?oAqko}{PCzY4n!rJ;h<<WUmpE?=tdCHs+^K_H4_=~wn51uIHF^!wC
z(WlcAspWd>eWl(RlJQgR#EMF|)}%ykRCsZx7t~EZqNAjtV(+;s6Dv8{##7S>Wsmh<
z6E02X7-aI|6TdWTP9^J&cVJ`?rmlMCozq-@x!#0n*UwM$)N&hfC>NS3#{BG=cmt73
zO|v61sha&)&GMA`YnrdT8Mw<1%shc*Ol~Y=Cu3(Q<~Z*uI^$;PTEj=%9ZB^XcZ3$Y
zZUS9z5bvc}WQkJ~CgyT7i-tw9_$f7FSw>FaRRD&t%1bx_K>!nlC?iIRs7BYhF2B>m
z{r!aL%Gxg6tCK?`7eNJ&jft`GX;B53V4z8yo&sR6oY#KCEu>mGZM$9LhaCwe-ey#e
zQ71v%Yuev>;ZNjZ%|v4Xem>v+Zlt^KuNJ)Otk~VuQ$!uh=UMv;=)?JGt9H4ic@v~0
zusL$@ouhI`oorm7-eSX&zZCU&zf@++Zkv4-#WCb1C8oj%yd(-Dh(1mK()NwpAzc)j
zXgvDT*v@FGB~PL}5z?JmmV59mtllOhJ1|9T@|I!ui1DRJh)?+$U^F^C$ViNwh_hc>
zVy^6t-9%jAbHipV5$&&M;v1KdPU4g7Qfv<5+QqX!2uI%<(A1t(WD;<T=%?^!U;TR#
zV|?amZKE61s`lFrvcEdWnuwG>ThS_59`a!ha_5|Q|F13S4iUm>=!wty|KDu{l3ww@
z#fXfG&eh(}^DOos7Ay#e`9GzVkOz-@n8u-#!Zq%H5m80ACUN{Ds#<`6i2C7JRqKo-
zo-MHN4i73mnm)<?<5LeRUXAmQDM2+n79mElJ|sTOjP9I$?;Tr{NLtd?rNvjYwo$OO
zr?P{h|A&bG%fL;-8O^xZ=>MRiR&lFtpP?L_g89&+&Q9@$U7OC4ijb{dWj^c0;+@PD
z#{IBT1xmK|(T;JghOBP+VBYSy&^UG;Hja7~JAZh|T5XA{X=WOLEe@!{Xnd6+@{~BJ
z7HDl@+nfsc!x8+k=f(K^gh7?KA`q0@ompq|?B3CWOMI;q*&+STr1>kX7s~~f?znXQ
z<S#T>*6>-@2(7Pbh(8vDx#;64)H7Ht&{C!|92eGjgM@xi@*<Z1r(boTaI=V0JZh`D
z>w8aRijSS>NQ(btpIXv?5`=OiNL%xqzHbSld*$(F0e5>=2e+hbvcSkT13uJ1q|XOF
zAlMXvV@3})yVo9!p#blysjK-<42|6d4O~)CTQuW0D=|9{T~$-fDh*8rM(=umvEeR=
zwxl9SwgB!Lh59w*p45fVucqffwWb#(B;q(_Qvh!d?Vnq4;upMs#QdZ^)}s$7{}YM+
zledQ%gcKArYa9aXnV-!k5~i8eKFi<_(_>@eHe;ZmD2X18+lPVF`ind>bkC1-c!9ih
z>ZvgC=WOdQMt)OoDVqezcI4}m=>Dc@TEsOe)zK^?ps#xFWwozZWmH<ci$9~3*!F;E
zo;prXyp}ZNG8S`Ks}%dygz*>+=}u}pWlN&b()CrBjjtK~%U_rKZNFkHF`SYwq=mAM
zerzM9-AG%T%qjAVh2EO)XynDItrGTI39rNFhz5tdSOca_b;w8OCmqw}W>Yqk6Z%Mt
z7B*Ar3`?K_#YON6NpLyq$0p9*;bj1xg}@jpEXt;E+4cJiidn=;vCZ(&`eOX$JJQp4
zsx4BJN2(sU?Xc~6y0cA?1q<P3Rzpsx?ySgN2E~xe9yck&x8R@#!NsSUy}E6B>g7Xt
z@Ed-)zatS$6lvw7Ze<BVnWJAS(<ySp3zOKhKLqE}+O_H)1BHoGwbgjfyF|r456=h|
zahwHk<otn`nS9=~#dR5!${<9DmnZ2`>0(TYTd%?u)v0;1_>P64$HE)U!hJBXR!Cdo
z<xv$vb!kd__Q2DMW$G*Zug1MXWWHd>3zEM54s}vK69-HH+G$4q)t$;`pqs*ceYBM8
zMO8@|{ZN-)@W7HK^&0ucr3_tNB;K3`&}gt)7r!!o4|Ny#;Q^KZE&b9#IAzvD+K`-_
z-ujpJ-ZvN})XKH#ks+iD)P6YZa1x}7N*Z^0<txU>OOjRn&|ggEA3^FOL69H{eJz+K
zCi4DVH7?zbCv)UP^u1F~P#x6rpQrO?e|q}0KbKT++A`2YSKOoz<fO*;UC(*okAtsJ
z_<{I|)C;Rt_=4NpQ~X!jzi7TvmPW-f6s2nO`tgfsa2o?&J}WvMvrQ>~l*+b;+e6i!
zCW;e?ol^u#tGcZbo5bi!fqm$zO#IVyvD@pZ#dJj~qa+lF)v@<qdk7I^4B__bMn}Sv
zM^~UrC#el{un(1WxRdg;+n?GIG}Yy1pM95!O|&6s-840-0ZCYs2`v{_2I$!^`!=XF
z)vj)ishX4Ki7SfWuU`hQhHS#{v@UN3#<b~`^hs2G1=AEn?~^O9C?B=vREE@&&K+`b
zy%b&_9k`0C#B|8<w%Ur}SHLlq!qVAwN1twKsz3%mRPrAw_3mT6_Qih9JYPXnEp`f?
ze7_(E$A-!-X@#v5y><9y7(97eqf|1H(}bV-IjV5URevuaR^i}iOE~yQ1bWs%q^QK4
zXFf&>3G@BdJJND!-q$5+QZ93`qH=cRf}>BvKO5eSKF7qfZj;2M61^-n>8~-hbn9D6
zw))k)DL<n{<ejNe_0(FO*7H(6_gpO3ZzKxVwl@^7ekw6%(}I6n(z?EuZM*eON=GM|
z`z6R*4BEZrG*42ibxlr?^}L^Ae68<(s)LV*v!sXWoi9WcaY8!riZAE;der7Wleo_}
zcy8Yh7~+miD`@$+s>p!7B-jf2Uff0vG<{oV6wWNT0(@ij%Sn(|9@Z)uBw@9qo^EgM
qWPV2;bl-?ymK^Khf97ke;*W|{%bC(I8zq5%zyXi$kkut&Vf_cWNkHQO

literal 0
HcmV?d00001

diff --git a/data/embed_chl_otb_nw.RData b/data/embed_chl_otb_nw.RData
new file mode 100644
index 0000000000000000000000000000000000000000..45ee217bdbf2c42bf406c3e178b3a216a0b77986
GIT binary patch
literal 5722
zcmZuxcQ70dl(u?bB|1TbC_(g22v+ZcWTSVn`00K1-hyoOmQ^BpuvS?lqPN6GZ&|%W
zZ%fww?&fahZsy*9@6DU{-Z$U(J<b$Ty#ET~!#8A>X69@?PbChbo2wD!wyaU5AJz6)
z^l7hjw8n`N@+eX%B_t)@>bSc@+}>5@S&3}E+v+hFd?LZ1VGsNkr$J+MsUg6C>+nS>
zEU}RMygxdsxko8aDz5%r@$G3T@m=*rVb<%!B#2*6G(INnk9{U3_DC;i4;0PxVCP6I
zl}j%TolI+|V_6O7U|~^9j?q~EVcQ}2*pg9~=KH1llAf?ab;ZSZhKsKgW8JnfM{w-I
z*2~%7&Aa+Be8pLq)9Ug@YdGJv^%1~6!=TnD3V!{e7YE7-_c<Bn*zBk-BO(8A;(;H6
z=&P)qyyYyreUtVwoHmTcW#};4BE2|z9}Xh0aI9RbEU3H-^u_or8()cf<gW}|74>%>
zS+A!FVGAcF5pr}eLMiK`X023)Y=}rf(d(!^q?caAt+e6u^M^EY0`Qhz*JXyqUj$-1
zc5OF3950X+`)P*XiS%y>Z6-BJ0+k99bu$JN6H1Drb>9?XZH`m>0x7`a3`&vxxu#op
z#o*NAza6wkHZ<hk6E99MV3fQWqU#LvS*hDNMi<r5%pA*#KjEC*RbSD60-~okb-Z@&
zaWK8ld<$OH?%98}pttmaHS%Pta?<<n4CQMWC*RuMMt0zAoqzvTJgLnB^G*uImZJJJ
zAHDHHR54(=<&BZUno#F)IRyWACxt52kb$k=PSdxB0BU#T0OxNfFdiej3<jG9jt=JS
zW9smZunQkBQliG+us!MekHaALD8uVJ!^$9?2~ckkovrdJz<vq<M(=w_nuM`#o9A1n
za295ePk1(B5B5E(%T_b9Z#6=V!|fcbEnAsJGBTQgiqbSpoJrt3<xmnxqGE*IaR=WL
z;>p>juT+J1(X*3o&PZfJ?Za(eB+zI*bN~ZlD0!IdMKM{@;Rlf%Njba>DhuNenUke5
z#Pz?ly(4G&a?uNjGCo=9EsGF5Z}hmoxVXzO!gT>bghw4QEiF$z-?-l^f$c?P1MTW4
zL?}L}n0#4gu}J=9<(G`e@Hh%8@fZ2hXI_=D!g2GB(WTRC`y}FZe6s!$xk}+26Z&Yi
zwLot^ZSSD)-C=-q>F2R{*My7h<D%~NiqM40-`L7uFwl`=G+XDBPVAx-M*e)F<{?BI
z@#!*oKIgOE1>&L`b&esEMQxv>N7Y$&e!Uv?LVcAn@+|{h8l%GTF=X7An0G5o%$^2X
z-qh9^`644N5EqEXd`ZAie!#tcVM4y=zE=l^F2mR(r@rj^PX$nuQ{kAq{?R{$pRmgV
z{%kStA5!S_qQ8SCiwsxNBi*vj<pMoKko)#;@`Bu`-4DcksXpzbxVH>S^NHFd*B>sJ
z6&%MSRpNA&R~nT`W^Zle{9Ohu!S}dz#HW*WM1cC2NacX<Xfl9vM}hSEE$$E7FMXB7
zGn1f<3ikT4@^@cgm+-UcBYJiVSPvRxdvMTmdp%8p|LA;zd$kEy-NWylu|^_?H2Ftx
z-qG+$SPTjhksyc<TT1gWdJ7!bESV9RxaPvQP}w{1UhN~9h08^r&NsSUM;Jt5%Co`>
za~8+tug~$!PUe9W+w48btZcp+Zv&K5+&vWHsV%KzO<+0QZAZr%VO4?~)2TTeh&Jrn
zl6H8CcTx9g=Y?~(_r=N`XM~b*TJ6IT)L-=8(;iOQe%IVLj1$>1sjoilL6bS$yIyiU
zOo+IZbDA6Q?7^>?t@D0YCG5Uy6}fv@yI$Qkx_jY3w_Qdr^}0Z`M6$oVihksc`9EM*
zbLow)HJbXMZd0<nOYZucp2#2OEnMJsW&YkOaJK=Xt<aC0N_}*bk<Bw%%@}>a^H)~_
zZ1gibwyrCp#dx`E<(E+VgR)eC(aB0AcCEbQL?WF^?>l<?vj|r2D9&SCIpa3yEug;m
z;zw6uc@IK9kh|k(*Z?f&sK2^c)^m+LAOQS`8Z|avyW<{k-r5V6M^^=|fbi{q55a~1
zQcP8t!muQ8zOvgqYXdKxN8Wk{yNdr*ONg8=7yt0YFoTg|{nz7eB8cBb!!{kn-LW1K
zdEq$hIT@n9p+%1Isy8s#n!YXVLQSRn=l<LUseh<Vkk!<e4q%9D;1*)H7GtzLN8eJV
ze+qKj(!QNc^z^p7nTv+=lUV#&_*;60;q*5Dnbq1kndFVm)()+V8}+X9)cIj@M`$rl
zvGe!O5|QkY-xe+|{i~7B&d6?+2vB4ED2`o^83|~WybD1r{bt8{mWDH#k^s#rb?B<X
z3eqxgvdw;r*-nV#aR${Z=kagwuPzO8*WqIIMxM=k@)5N4a!L`0@87=!D+tY~-1Z5C
zRi)gb>CkjJ>mw+$tpEZv@K*3;*f^EaQq$UXdYx0X2v^&>_7)%~K*$fa+6+fL4_L1n
zJd*z%D&QhLyJ3YH(I==d<(up`Hy2TbOT=Rq3|~T0hvX7#9b^ufTg=iMKsmxIeBzn|
z8$?$zQC2(a6@N&!!<xyJf~yKze`c!g^qnjk2QII^7zh;(K3Tk4rx<^febdi@I94rc
zALRyt7lg!|e=ZS%eU883qR;Tb*)4t77#oG(toIWQy2_EMUN*`zo-TnTIgznTXRSV*
zGp?`1_S@fSS6i)heOiSwn#mx;>`NxYjWnxbyi#Au5UtKF_-wU(v;^0)^QgQIU|x+J
zaQJg1#)#PaZZYgu;qLAEt562n(B6***-VffbP?I0T7|b^XSzmN<Ju8+8ge0AAp;u|
zwd*_6x*jy5VN~u*SFdu58kiX_+1tHJH{HX6Yh3q+mu$h7Kc1vst>A-C;hXk@)jWu>
z?18nka#O_K*yN4E@-NKDab9Z!CH>ky%V1qSsnXOqPaT^mwM{(hX{5OR4w6plb^Rq7
z^|KfLWJJ7At_S${GsDl(oF+cj{dJ;W8DZ<q?FlOR#tL++x0Q=+_31j`*0c3W(Msv{
zh1%9xMqX!3{=+v{D$*@|GEP^lhL`NIkb<}j6|=4QeEiWDqD+^nsFWh&&m~1r)|e&X
zeAEp!gX7^V^ZaAFo<hcExrlLN@O60oVzNI2;+uU}=~HS<a|%2f8Efc+HWI>DKx0NU
zzOVXquqU0V_%9b^t;Ip=Cl{xM=GmoQ@e*I=z$&()MwKN?!imsCWE9DfBQNp_+OK=S
zVpQ85={Rhz*s_n?I$NFxqdex6b2ZPp<>v#76~&b+%u*IYq*uJ_`3YlvQSq9W3`C&E
zzr|sC7a}uN_Zgty$r$)Njm?A6JpK5^IK*m^wTCTbkNrsgZQ>{O`KwQ?k!=)TcYQIp
zA~NIsakU15zrs3=Geq+ed9&Q@DsU{O=UVh1*#<NP=ybH3N)mHwwBlHIo)eK{{ey)g
zZKtY2vg){-HC-{6ve-UE=}gn!HlzQ0)`xruW@N0rPjreHp)f24N48@VrAT}ooXG-Z
zinb<b&iyXxk2fFlT{IYqV;v%cP7(?L6sb3KdR2(RfkWBQdrEP(V};V$>C<M~G@CYp
zc?xgCOddc9CH8HmNiS10nBe;5)K!!s)irM>3g4Q%IqQ42zMIngLr2cuZq;=BJF2M|
zRag}6;gCPniy#>9SfbR1D465;pUiQ|h=J5!4YJXG1WP0v`%n~-IN@=istvw4?eQTO
z5Kv5rmCBUCLy?>YvZCAp%@k=ihd``RCJssx^%!6g)$?G3xHm9b5nq-#kNDbis`j~Z
z*~tTs%w<=_0S-T#slga!0!QXGF)j<d9<uWW)ndFm{2{&G6x9Zzc^%BGJoG>tGizH6
zW~07-JmnHq2n<pjLPNa(|H3wf+7RYZ$&O`vP)QcM>+RQYs3P(H5YYjY8P4Un;80P*
z8}7Adh9)S$UFcM{>K1+k(*^<SeH>^XB@10)I#=Za5#acSf}x<O9_9@YG!k&iwdSIV
zA{Xa8CjS?BG4a{sI*GKZ3xdDqnGOITp#VOZ413E(4_Ku>efIC<A0_}4Q3&r4NV2$o
zCE5KWRa%w8oz^IXt1SG7-BKkGe+hU{+Wlkj-z{qrfP_n3APM_|m=dF(>*P5D%@I$p
zjB3bK&u%B+p7EF-T0(S70j<I-CKKm8NapIoKbCLws63&J{w^`*rudcdc;ZXArfBai
z;17TfE95{Jn*j);n^V1a8p8(@1`)22H`?G`oAlcCI6aPDW+Yt&Tw6CPP+q(BX5*m<
zg9PWc@ggXW>86BKy`VwQME-pM=i<U{iY@eO6VOgPWoX;vIdK^C8tH~?V@z<-jr5vt
zpU8t_t=<YiB+z?}EYAZ#(vx{3H}uA~owGj!fK!?)BDo`9+l6ietnn_LptE>5JVy#o
z7}37wd#W^-KE>fZYJ||*A1IKxaEm*nE9uEQ^4uA{HixK6hgiWqVf8G&Jz;|jJdI+c
zH)j(f!*U7q^bZf?6>Z9G;z0Pf-^%uXe=6I%1yc^1$`<#ZJV>wdR*b_4up6!ijjX0t
z-DS0;t&9C<M*e>rB4vtCqZ6Wo#O)pr%U0+yyO@C^(Re{ywI$2peqZbecGQFFLJXej
zHW~4rQ<_(2&y0Ro<IcEnZ2nRByj;<-u1<LK%OnLgaGrKd<&}(=;q3L$m=#ISAsJSZ
zep5ywP74(!=}gYan_K(wBy|h3d$=j5)ljNwoXKn`V|@0KznUk_O;FWIo#VtCnwr|{
z<MK^sDu3eKGwGjHBMm2tU7F9bTV>PYSvrqkZsd!ME9JGXTUkk6J|Y$Q^8Y17h~kd&
zGM9se3C@+P<6l2=A}+=Ps(LZ7xj)$YaN}@O;@sS8wd}X*3oevPti)ku*89DGP8J1X
ztg6AZEW!+XoKMf#g;KKkH=4CfG^OV0hfEehDWx3OZ@n3zl(vd;)R`bIDF;%txPw-<
z^QbS=WMyj|KgX*@9b1aYfV%i*dYVYjfGzC&u47-ej$^)nmW+lRy~R&c&m_NTK6|o?
z#3?qGv7FxnI^RBdJ7a~0m`(cFNoFS8U#6HMn|gC!5=!{{@D=xhQd&kKWd@e;M+`9o
zAx=$>n$OwKdi8uNh{mEB1(>q+Bt|P@A0;IxIp|1{NgT)isP?vNG=XP=o-rE=Sd~J=
zHIwK@m<P-lIEfg-SF9Yy5SOnlku`10gFBLmS+>c=e`>2b#XlL~D)kOud>|1xE46?4
zfOy2yzS{6+{o~`Og-sX`IfES?%gOnWI4c~+njw~y#gA?va)#4u2u0Vl6K2DZx(s?T
zr||xu;4jjmaK%^L_EOJU?|fW%5Z#LHu<gs!CF|nFQQy#Di0J*~$|(KPJn8=u;eW~Q
zRAy7)|0^4A3QPt{{}xD?D+yrlPffO4JC8}fTC{u=(Z*wZ)fKHS`cHyohngVu%S^SV
zXX>9KsR`~pMfdXWuL`FNYm0na?=8TGp7$Qs7@c42zt$DME*Zox4tp!GK%RU(bVL+U
zwrAhC966FsNR0g-XZ(8<u2M1XjMDdiyz<}f$cYp*yy*;y7nrU1!+COies1=H!z$C5
zJ*XG=%0|TLtYmH<*@U@|mM10sE8b<vQQS6L_ca<-PHh-kWb>KYaDHFjOlsw|z7%C>
z3ssc7fNOY<-=mNfukbfAzG;z&!@aosXjoxGL6HeOQaM^WQqZPI3pw1Z<A~7i3#?oH
z`1VO?CxBXBU|4dW2wf&GYxo#8!r%uLQh+v)8L4d<$rdpyrm<PeL>nxG9Z#tgAy_z0
zKhO9V-5mgoYh}2;AR%!z1OVdZ_!xXoGFJJWQeU}PEf-n0`&C5>4~ZEseAB(h(S0Va
z2K<dWc59Y#$Oi<k$i7n8z25RXDS_Ny57BhY*xq`moiEr7&K2exJ-a&}7yi)b!dqHj
zZdU$?l2R@xP;B_pm3M)Yc#H#oGO_G03WsfpMOv)6m1-Pa6~dZbMU3p;^xAvws4=|9
zsYI-k2YyOXKs$l%JPrJPn5tbtw4t$<38Uy3wa2+a+Hf7FY>*}yLkRn~HU5?=q}LZC
zG5J?;6H%3@Q)7EE!g;~3n@vBp@??#0NVA8+&aJgs)=3}A10qr$)R-Czlxql@EEUL`
z2tKIcC+_zHI@I4jfk+ncetk1;s%<t_D&snKP>-WOvzv5)AV8Q{AG$;4bQd8C$&sIx
z-;r(-Uv1P!ti6Z!ze(J-K+w~4>Q>c{DJ6yExA?J%jnql4<c0Ef3o`<JI8{^YDpB<o
zXi*Kl`4`Z63y-45>hO#nd?o8IDYj^d1z8D<z2M%NUj%Ae1zt3esG=%vT|ep4%&e}I
zGdfuWwb?4YR#Oz3s=M%Gn^J&6YH=6wcR@z)(_`0l=TjN!JVoy{NjF;jDgr3e{@z)>
zKcG3g(dp3BuO15MzY&x}k3jdbXrF9_!_$P{$7-ADB<hR0V#u~iCgL@i>)w%oM%6Za
zXhxngYlesg?^#qQ-o<EDHoMk5_&vwcJMWX`)W4nT27jmiu5U&N>}y!ulhnPa2W-~@
z%J_J-<EgX{M?g%gsHQ!+#&>uY{16i9%xddY>Wk&Q2Je@7v38!z(|7_V4p^XUrKyk*
zuH+l7R8j3}-=1XHYqwpJQz<^B7oK(0oYW_p%*nF`1ho|zZ=KBHoMr>pz=$ps?n+I6
z1S<y9+LJD5hrtTk<}t0~D^XKD?o#*6p*_~fgW-y+q~|Ug+Y&;k50ng=y5>KvmimYr
z@ucIPrO+i0(6O|*7R2j(q*M~K&;+>&j*f*&%|9aeFwSIZvT~L?v*rnj)Dak#m7>Tc
zgMSv!V;0HigPr)oq8cJvHH0X;A>+7jd(}1rj!$cBAI0-;H0jihq-1KF^)isThSAT#
zcOz&k5sjABMx>YvHN7sy?nOJPiS<ufHTiGjvjjJv*F1LQQlO2wix!k^X`gZ_<yicz
z#KS8-d85-E)t^4L1Yf{qHqF$l%)sznIov~EEc=q1n+mEgP!ZR_MZ@4)+)4KND<(C`
zuy-Tg?;L4<dYcSQS?EfkX_Nb<Wa8XAbslh2Iy8Tw(Si4w)#V;OrD-$cAx;_Qmm=YF
zDi1<76DuUC+%LPDl&z|E^%Xa?Eo!_sYk%e}y2LB}Ws@^7Ax?|iLeN~z?VAMO)V*!S
z&dlp7(=H{F;E>jz0(?zPP4}8!bCVX-Gh4I2o*eL}KY~ZFXPL|!gxoh2ZK*eKJp7Cv
zeQ`cxqrM3|J->tmTKiAiV0R?9g-Tj%E^hmu=M_);-ot837$l898l3?!4f{d?9|#_(
z37Tc;FkBV?Fk<@^u=K8#V>Um)DASG(p6_=fv_>8vz3m{qpF<St`?G{0;kV1@?mx8$
zKnyn2Av(83BlBIM-lzUed$qvP(Rvne)<bqW;jCot!tR9ImO#@>V$CSIIagGaGKQzP
z)uj3PIMHgRF7Kl)q3YH)hsO&UMT{SFKFWtFN|3&t83;~EhLdKF@_g`lW~1$fU4~r$
z2pMgw`FIrQ6mJ@eK`Bu$U7}$%oc_M@&-I4W%$q1545vZ_ot$ZK6zgvW{oB*<@csiF
C>qzPV

literal 0
HcmV?d00001

diff --git a/data/embed_chl_otb_se.RData b/data/embed_chl_otb_se.RData
new file mode 100644
index 0000000000000000000000000000000000000000..779d25c6381fce29e50e97824993b48055bd8457
GIT binary patch
literal 5716
zcmZvbWmFUlw6$rZyG2yG92x{+1cl)R=^;d<LApU=m>C+Sq$DM!V?bi)8CnUEZWy{@
zM!N6&ee15f){XOTpLNduv!Bfs^8oMv56QuMQZVB*`=bGUV)0VeR<$@PMm@&RCBxCu
zQOlmCGtpK;5tk5Af&#t5KGrziudfShBe74_Y@$OCoc7`RjFs4zpEQo5OR;`_cW>>d
zd(+SVdYLETvK5zamJeJHKW>==3?Lr@XNK#!A9RJ*)QfS8_|JhtWnFhormcvq9IK^=
z9%yzM#b|#QfiqIP?$i%a;e4!{NMv<bxlz@yU3xmjltxao>UI$9=@xaB>dfb-`Z=st
zA!#W_t=}EtL-w|n-hwjo1G6}Qs-!=fK5}9Bby~uR=dNdksL=8v)b;KBPHD$tzIhPn
z{B^?_wJN7b-5*y#8EW&!=myiB&!s-lx-xcDt;prZ5O#E=fZR=SCRmcOA04qM1#fSy
zy|<6}RJ_JzJa`zdfMc|1Yu+b^b*ZKBp&OXxI1NpBWa6Rc<bkFJm~Dwt(sD!SH^tRU
zj}-5T0LRsh*V!B_-tqp0v}rklY@TDCG{8!%y1J<+H9w7lm*W7&g>kQ&qE*3`{ccQo
zCIGaoaTB~dOvT*x+UK^ajUCx=SLzbheHouXTkb4bbuq=r1wIz&IJ)^kkb(@}{B+$z
z4*|%YX(8@bn%YJlA@T+1Jv~RW|B&Z$;oj$KTEDDjQexGcTlH7vnbXJ3$}Kzg>kUQE
zrrIA83o!IU$8Ix;g=Rlo%Sl?Ym1{+7?m=pq%4emO+&rA9TVR1BA^eW>xWylPsyWsF
zG+9qoH}yUzH#3Jxw4JE)-=OU-LxMEKl~A@@jq6>n>w9iIYzh&=S0OO00q&$x!28?9
zS-K=wUP0XHQTz3_DERguz6YuDO_xc}t7l|Gss*{Rqq74P=?VybfpRj2av3-LqF5C3
zds02%#KpWIOwBT!!-?B|_YU%LEU}CEE8B%iX359bV*rO_ktchjecuF}zXafZ5OeBw
z>KS6vWlIV^r&9#2Sr_P$<+KRg?FI=9LzaxR5-FhATky!hRJ-Tobj9teYKh12=(ZoB
zBD9BTd3c${5Je7bQ*s^+<r4C=|Fk`INQ-|}fT+12eG4pXHLRPTq7eW^NIHEDQhfs|
zGJLcm9OM0``)xa+AC7y}N}|zukgVe(Ej6mY=Yv_F;Zcd)r5elHTt~C8MNe$+<R;*(
zihjR~w>Bf{(uwYFI*KL8ueS76>-5Eu`%FFKV(#)LsDC}ARI1$Cegn9x>!>$XyRM8A
zKtLmNotB-lADfS#Ln{e@V-H_QT*Z?c4`|o_mDh`Let3JR{k5TcG8ZvdQ!PVm{KhZX
zIe6k+?j9#AX`r+hQ!sK1MxG}Noz$#cesOy>S1^z_|7^6qe(mwzi&^(x`xFDq#-F%I
zj*hQByq-Zu0?*R#ZLegxt2|ZHCRgaHaxnvDiSmCY$&5LL-r#^ZzG|}?w9xlGp&S3P
z=FY`4AGHd<N9I?W7x>^dF*nVlj%&TDbw<Q`oL}kg4UT4HR`FG{1?+`q%oxXqbwpy3
z`9EF6Y5b<S!?k@yp3{CqILEL-Blj%7`gay);`>+0?>N_+3bwCGw~)ggorkh~TG6*e
zuC24sm+6+YLa5Wawh9-La#Y3aZnAk~4&!SIn&4*A4fkJ1yZaps7F$#~e_y-}p$PH~
ztjCQ}o7I;?eO%6yFK~ykp9rH6XjQkpie+4EpWQh$vtGNbllW-D^{Uf_;Hu>aW7xe^
zr<U5kd@X_4(3;S7_U`t!;0x^77^_n~<tpsr_;s}ZGZ$_EQor4!HeK!NxYR7*4_^z-
z377MmHG^<fH9(M?iZ^J5XEMKc1)rV7?5w_*VFb=~TvD9gxjCH23-2~lMd|v+Y&8d=
z=y1boSMp#B)?fT*r#VRe8;Np{9h0)l$eyitv5VRDVt6d&&vKmtq#psz&(CuHh<?2=
zUtbVqDbT<j)JQK$CY5{M?TVc?*U5T7yD{&nnUx*xf5YNkk8h1j<*odFUnvM*Y-TAf
zcnt2n>$vO<Biwutn8&`v2s=s#KRk9!AF3;2dG>?+(#~Uv*6hoaf^V$@S;v8Oa>;gn
z&$N<K@Xcx?nOy&_?VI}v6A0{gBQ&+K)1hLhWqh9Nm+q+uKT_w1HdZw9kx793CgZY~
zeOH;A!LZ)@mpdW6XLQ?t%yY%-dJmff?tgCl<B4nXhn|iH9(}vec33ekQ(VAavh}uK
zRo$|E2u1^wxgSP&_v<|L-9@`TY!@#wUh^^rpB}>ZPp_}AqVmO9R2h3eva-3&a(kws
z6dHb8llYtaPX?(~LWR6rAFq=6UY#`v^Y@P=H7*=)`DkX-^?Zs3C*Aq(2I?M!8M4de
zzkY=1eC{kG%5KPmUY@Kvt2XyyO}$TVnC-WcTNZhP+sylM3yD$a=kJGP#O@w+m>i(?
z@mf5N$4-fMl+4&!Q(*%i*{gSv4!4pOPG(=~#*csECNH%4{O)XRs5??KkzH%9;=k35
z8z~Lj5)zpH%J4rFc&R$;{_ysNC9l&@4igDPJIykr{ttAQga&L}-Zs6ZBJA^P*}VDc
z5@IF;`4E~^(jz@+AP<`wV@k_?z{q=${`pw8(e&d!RxrEu_w(SPS%gRR)^O4vtQ5q}
zFBFp)!xp&SmIY}np(sbpKfRCwrs~A*KAdoy*9I#9YYppT805V_eQ-l4g$=T2LmAEG
z!*lC3t4|JVXpKLnR64hxF+Y3uIB=QMlM?78eHnS4%68F7irk`Be;bw@x_R!sSVy)9
z^8naQb%ViK3($c^<;BsuX4mW4Q-?Sqrpa5>47aic|HrpSoApd-37(Hh=ZE`!+clM7
zkP0s7D&77Up2Y}Z08~EDHU}@QMNv^&--80Uu*ZZF$U&EZ`erYOo8!EG>ofq!&0~cI
z$s<7|Ukm4|M!WYUVsw9r1L~wRUD~*wAEox{bla^*&1RtLJ9Ym4sHAT=FgYKmhRw6g
zHf=7AWhlPqTODZhmW>70aIy@s)(8lCJ?LN-?`XFRelQT=g2PM<k%LE3-c3RL?%=mS
z7WTEzBu+eBeikouJio}keS4R_!kYHt2$sMlNiuzgH2Yz(`Dw3Zl9`E5yI+UzBS&KH
zEnD>xFZ6RxTpzEuaS`B?MYN=v+W|u|iqg*cjPL6+P-$h$MTtgkZ>d?6$wd#8Yv9&4
zV_E%PB-@_l&bmqppCwDP`?@}e(3lJqzQ8#Nch)~<H}#7~>p-gapYIfOHh8k=q@FCr
z4<5z(l8X-Y$eMCqfUq;w0xZ97uVc+Z2A;MQiyj|?dN;y;Kl)ZfQv_60aGSlYxL7KP
z{7q`p=DAZHWYaP|Bsr$ZZ2jpC1@&C+bKV3W?MU;7;HUX1j6pH#&6!{)GSEQceeve}
z^Ff7R=k;b{|EZ8A#9~anz)H>d+!l#zl00W0inb|5-}8oo;sLucyq#@77!HC1DblGS
z&t~l4aZ<2Rx*CRI3m@WZPBa;Mgdhcvi>!(K0)URkzahIghRecXcz(ojX+sWoNeq@O
zm_<Ag9*}Im8Lp2h#0L^X6=rJT{czNSbmAKTd{PQFo74zUAoL?$;`HFbCU-K!MM;Ci
zF(1%JcrWm7$ZC@D+zIblYvMN%CgW$RGH}HZE=oD}3qd$OIx4K4(Ss1%iiaY-lbR8P
ze+&T-RhC@xyhHm*<xX(y=fiD6a!J$)5G2>!=#h|RLVs3&%5-uF^$af@;Gon@Y)P7_
zb4c3AGPBm1L&i#ULx4^QStdAxqWv)=`1W`(u6-}GTBeRC(f*rbK{0K82t4<}*+Vof
zrV!7V@Sb-@91aRW;rajT%P(js1)I!nL^KlYlia*~Nz}$QbBtGn1`$@{KkEB62_yLT
z(+-`51`<^kT{^#*>zHESN1`7Q$&p{s;N9axSZ6fxrtwh>*vDA@PHloVLI}Pk=Wt1v
z+&{kgmlNieH24n*Z@@3<utGIy_-=Sm7R(cjRrYstVkM4!YxD=H+`hH(?U%=nRU?09
z@M(AXsjgh-9%I^{Je)9|^Tz=3pv+iZm}4jUm140foa!{=O#I><=?!x<&93YL2WAbP
zKzhd$4kU$0VA<dwADcaDV!^Y(U)gm~`j&U;!;z`uPS!@jsnkw^ZHIp(4nNc2qCV_E
z>x6s!f$n(2S@5~Y+vvLhM0ea6tF|I{xQ9jPwHpQ%;SmJ~lH53A00c0mt`7vcR1l?^
zUQ6Q;;5G{zvk!#r-<)tH*n^NG5(t+(BvJv}{V!`rDsrsTY~vxhG%ExlY4`|ys5mx{
zL1K!4Ch*CcFuI-gP{9L80C^g&Kykx_8Hb~YZ%~hKxG{rp1hGqA_xq^xK3~#I6!qGA
zAP5tLw@l+tgK^}v{+{Lkv`kku&X%9qaYq<zP)Ey~S^75H318zy!<Le@Cs%D#g2@6%
z<dB>BoG#`%x;MDzst(p;?L$>!A$Qxg!=jDyxQZADWJ6a@26I(WQb@9rY3cu`DE|8a
zjQ;j-W>mDlSAMLm3_>J*xZ>y1t`gRzDEM&8$}-?ApL<1`-8FFbQ-o9gyF7i4B0Y&l
z4b^at<dnl6`i8=(nHlVtCGN0pH5sDgBV(7voJ^`cc$1`ZF0|lmR-|*>;H5X;v6X&=
zDSNL;*x~|d3zft`>08#<eUh;g#dH#Lk()jddGrPG;q{T>Lzaaju}oGJ@}oQSsR&&y
zUJ4wo2v8=3y;)sU@278utCzc2{Q-wolVs^|Lc}v2&8)u!z+~5r*J_c1|D<YUl{>YW
zy-^#(=|Zc_os7@%64QR#U~BTbA2-;)&yXL_|0bhjH2mG2vs8CLN<Sa>L<%YXmzcgt
z<r((|*tZ={bb4_7WFmwqE1Qt#Rl3e+qD{8E>iaye;`{AY7ypC(=Tb6uy|uN1q@pXc
zWuss($ZCSt?>@QbVLbl&pUh-$w#)s8E+!km20RBRi8IPHPh$g#CA0nM=xgQEcZf8+
z@XKuDsdqtA6@7y9*~?)v7j={YeO$1boDW-zMjW-KT0)9FY5}!2fiHQE(Uk6W5krG-
zifNPINDFs`KmC-6v-!u6Y#bLpR%4m}vdBnCg)P1$No2x^%_f{|Bi5=co;)nHV<}qv
z12w)>jf-Y*)ipafitRp_dw%feVP`r?**gh48p}uDhr%)gbTfZ8M3gZ%kn`~>cfOGm
z&d_D!Nn`+=3)dy8)#B10Q}75KgrnQmUb7d$Gr7|1>#53Z%f~f@V)Y$P9LpHli=qwN
zCt?y5R6DCyn+;zeKW`b~7KC!%^?Fxapr(X!Fxq<mh2MVymj=v$uKXWs)MseGso>`{
znLkiu7{gQw4Z{`><=Mb6@2at@0(*}07uEpZSY3f4>yRYI&k^|gYsc(03g_>Txcj;|
zXSa7b^Iz|CT?S=u-64y8b%(k%$b_LPn<igMWc<pM_=0@oyMB+obkvl5B<3vf`TuIh
z8<fZO&;MDA4Mc&(-SSytS(QdNkZ*N!v$)8#R-jIOM*Q0599kxJujSwEdQ=L`ql9$H
z(V(g*0a3Cd-Ce)}O2Co+eWOGh<B|uV5puplIC5vv-Qc@2<&~nfJbmG!m(9GG=OWZ=
zJeB%B?7gpW_d3{V?U>OAn>HqmwPh6*1{GH@p0&hPq_hhzSH5x%OJ8p}C!WkBd7<<D
zi%?4aE1Jfz)Zfa)pI&WJLa;vbXOWL6R~th)Z_~Ju#qZ6TT37P%{GtS8cWDFXYtR$;
z1PuRlcZL=}(Bh00Yw}mI{qybQUO2Z`npmouKy!CJ^q1MTwp5Lm)hcVG3`;OO**Ku9
zxF_odWp=Ok121C2_IA%Kc6a%44V83W-ZdL3(r&?=DWox4|5K*e=!WH?b-}@ey3{gN
zItRJ=$}3nh94C<PChDKzAa>QhX$NSDd27PbS}*3SM=^LuVX@xz4R1{w{+#Y!mLGk0
z$Wpg`P1_2qe_-U$8m?Zhn<&)qq~uEkS*qG@w=I`<F}r(;?sJVoBWH^j229+&)*en}
z<VLSAlD%-?UWYOjZu`lK488dD%9b}zROd`4G~(2?KDX*e=|nu0E>cY2v<l_8OXOUX
z_~rBJc!FMfd_#-Ge#-D(%VT7k%&kS0J0g)jm?qwugTl9OP<fbY&0ZcAopN7H?T#6|
zI&Y}}0Vp>@k{;LWJLYol6!ws;PD0C<iD0O6Vw;k7U}lHqFtgEj*Ro2*!rXau9rRgo
zpP;nt+0*Uj8DjQgg(hGA?b>;YD-+NGf<bv?c+t0)_rb0^Ob@XB2<tK!*hmOFvqXyh
z;6(cowHBRCGmWh&o14-V0|F#t_%iu!ojmHTXYGfHq|VRmBWzn|!b=&ws=hqU@<i$j
zmzt7HsfX9wZ3?x!*Yc<r7=^fB`O4&2-m!^*!`T~D@)Ah}o_3}NyiaGlTz4$dp}Q2$
z5!qlbW4w{ez|2C2H@c=z91&lCs^d!s?;S2mQm-y}e|6T$aMtmh&RJ^<(*Aq6vvuTX
zqTI3slA86P;avF#XZU(tID}KQ^`F?q%*;Y^>+e@52#3MOIFWl)vdn~ZwI`wT8rJFi
zC<m64PszZy%J&=})~$8vTpZf)ugQ<kV-%mQ?}UL%T=&;nh67jPY5MF&S=Z<g#nZ`0
zMF>BT3UGAyPJE>*cGZsNQK1eYiMx-my$VZmh)Nt-F|6p}LpVo2Gwqa>O0CJ$59V`U
zL+Ay$wOzoOV(+aoJP}N$>vHy^z3VCGfzPD_Knq7nniBdIFE>;Kg>{&Vk8Qk!pa(BZ
z97|;4#wKi?Zxxt93QWm#-;;2xdEnQjCwu_#Grcf<p1i-JmJ#-H2L|#cB4D8iJ+)`d
zgo_JVqf!T9lM-b{P9%v=lXzu3mUs&E$-AKga&^^WQ@K@{Dy_eGypXi1dLXMH`$K!z
zj=`CoW#9x0K`Ez`<EzCia#aIC)vI3OvNV8^P@+D4AGrOS1}87CVJFWov9F!qztAE4
zOv|_#iTB*NCS_9vz^pwyj>0J##NaZ+;kL)7`l$l6$yTPsZw6HLe5|x&k_^?GBEjsk
zv?JZF>UCPSAZZGQeOcL0gq(k;v$~DGgG<3q`BEoLgxi3brz9IP4vZg>%5H>Nj*KzU
z(b4yL(%G-TgeN=F4P^709Ek%n3F<{od>ts~!iJ({&gx2JTE%t|>Zg<R&*8=PS<lJe
z7n7iVwt~Q6KWR4#hJZlT!CTQ~=$WpnME>D%FW0n@3v%iUS(&a@`NYDDqGXxMjV4a#
z=hZa3r95FgLRH)B72(2*vzZBRD-E8^#AGdr08^A~s(hnrfYBU?I5MMw#&PgEk3MH5
zEkmTHx!*6$fUh!K?I}~3LJ-w?5TZC+>x6*gLE-(~<!2p#QE*a(si@rBC#litfIQI?
zD(+F(c0F$q&m?9oLmQJf{F`wg%`01lZr~)0mdTDywsl7ndypBKgo6r?f6LsvA9lnP
l13KPZmd?+FSbMou$6E@pjklJiXJ}uzOS|?($;RX1{RfnIXh#45

literal 0
HcmV?d00001

diff --git a/data/embed_chl_otb_sw.RData b/data/embed_chl_otb_sw.RData
new file mode 100644
index 0000000000000000000000000000000000000000..6491a4854f7a8552f0b748223c3b8a068eb7dac3
GIT binary patch
literal 5729
zcmZ9MWmpsdv#9CrMq+7Px^qbhDX9fqQjic3Sh~AQT2e&m5EO)^mX?xkB$n=_k)@C4
ze)oCqJ!gK-Ge6#$d1nBLcxeA8aF7=Ge%dAs9T2X1ZcRSUzeO+E$iWn26{Q~+vI#i3
zL?v_MAyc_Sn!qdGMsQ|TRcQo|5tb^y`cFJ1iFh(Cig*r+G}VZb^w8twlUK({@#jw=
zDUa|%hr*MP6KTZaxZFt#!d=9D@2C9SbPaxl5c|leHlAqVPTHDIE@=5GI6_Er_BmIS
zv~_9PhZ4S2Dr~%8Q)tmLLW<N%WM!jbew4{S=T#2FJmo^^h4_wX#zr5cdmZd#gRYWV
zebnmcf>6_ER=>C!dppZYmV9)^`nW}JY38wNm17H<Bo_0+^j`_4_<w33<Kbjmku$bb
z`;|C9W+v;zqUm!U;B~GnHs(I$Ds%OnQGrUu#WT+!N^4`=1ZTxEwn5T%+dE^lhxTt>
z{)BZ1j(kzGq~Er@l~>Son!m%qOi=;Q`y!?Mvl4ZANZZ*DC)$8kAFJhG3&DHtuq{?d
zHLNJ&ZO2Ut_t>v;GmC2m__fg3Nu~m~#3VhfRK8##<S?QEmUQNuBS+{^dnTmO=1E|M
z@Z<dX(~f_P?!L?6OmgYWL?gi);eS(9V69%eZ22G@Fwv;Hl6Rh>*F8BzdCP}Lw87=@
z4rtmiG=D8+EH7%qTB~wVS62JLm_1M%MA3}o>&T}2@viea8x;V}-k%P<>dG$6tMh+)
zNPBKS@h#!)j_~u{%^cKhY>>WPs|Z){5Hq*|byg=OhAVgc#@@%#(lB$kwc18wS4rxs
z*6)qBokT|{vUn+)cWgbq>>>rBu%2?E-dNZa@=INcn`6bzcu7<Ww=%(3HfNo0z+rv8
zG+Te<$GU&}-1hVFI{*3=*~z^xe3vQ-4#GS_@fKisC4;=jK<F#jxrq6}5TzvH=MU9?
zLy?~EgnF`^9juV`LnO4M!jz!DPfWYaOQ5N<9*lx4{VT}sW9K>|RD+&w`kJw|(R*%m
zzx_nt&dn^}ciC`adTcdY@?Hwp$o2DX3(UeTbBAdYws;F1;{^&lgC%uQTOoPm^gv7|
zL=ak;7H43{=@HYAED0|3Y@s<S0$aqqUCb#)ypv)#)`y?RPl6g?L5DDKKgDm=dStD)
znkCGVX#voAm*Q;8Y!&Q-#Hq^zubJjL^^$lK)LqV2uDL^$131NrnWL(MdNr41`1?!C
zanj|EA!``<zPbsu1>G%UM(2%YM^NRSKO-LyE6URBTMq?wOfQTQF!wY16ZA)e+xUwy
zgsMQk>#zUbn0wjnc2fqL{_VkD{bYaA1|>F`0mQAOrX}&pFRzB2sRf$fZtm|RMzsex
zCI28)2_vBz2bIhD+1cxD*s~!DDR&e|;6?LMeT{gl0^g*=^sydz5~z{c?n*i=xZJSp
zIRPvP0!7SoyJ86lwr>8kJ*6-Y2J#%=Ex0#(C<T^#o5=Tb(zswu`Z_}&+t;I4zXhxv
z-Mk>dxebja{p$x)T#=!te!pP1G3FLq_)0<O$FDvwu;g%UTeCWUmFpqY)T@l3v`!>L
zhYv-z+oHalf^y*$sP>?vq?M)1ZaMk>u=qr2zU+U{m7Yd<`T$)Djt#Ke_^Pf^wmIvL
znv1y#NsxQ1|LkP-0no(!V!Lu3w?hZ2Je%&w)H+91oUT5kz5!$Wg&W(vDBwf0i7UU=
zRrB7&++CfmLDn@^%Z*k5Ip0UI%e0OA;O9x#G40;@J0sZ)Y8pL1EAbX<&*lz}4vSqX
zIzJ>GSglAq!*A=Q+5E*CG)=;<+#wkn<rHF3zab!?EWK2I^g2zF$N@5^3VSTW?Ff}T
zEvG;HyZO~ZJ|wcxb>L>)_=e=%B8u?n@35!``>9+?`oKKD{d#7wD6aW(3*+`N5mVdG
z9oowUAFI;cWYvl}|K5#qN{1qXjC{gJzwa1KJiMzNKYSg`rE1j~Zo-(GYXG6wO&Kqs
zlbr?a#h$(HKJ^C{cQ^8=fZ3BY={eEUQh$wV*KcM~B275DIVR)kMp}dkKIRO4l8q-z
z9j=a69T0omU7F&%Vg*0~K)W!IVw-6w$L}kYz2g8b*zwR6i&|ZLWfUnN|Lca3orX#}
zcC$_EBMK*(0~*qE+^GNQ$lWq${sT#^F7&>oDOzD19(TMn6x~B;sE-WXsKXMy*I$W|
zjqM1?iDM>?`-DUsGHbq+UHq$39w{{*03APONLjN{Kb&B^@hae!L%vu@>V=qfT<6mu
zsI2m*ARzAm$)FI&$NEpvffL`Re~SdVk*Mi{hporo<DQR7BumLO29P7#f5vBDoGI-s
z3iuv2@R}C4AFKp3wKDZO4iB0gB0KwT240N0rdWkv*kxRgdJzoV#QKFFJC8v1r@CG0
z-MN}i%M(4_=0iRRx?Z9Bgw0|SjES}1d&*(Y{*_vMPz0%s51Y(4C&epLo5a`>XQ95l
zt;KJcF4)}+WF0HUG?dv$PwtyK+x@!{2;<hyt$sBiea)zAR}<GKeQfcu?Q<Zl!+XzH
zTAk!cQ}Y4-BMFCbRT<@DtpU4S2(fMt#jacBjzvJ=pX0rA_Nq#Vblp1}v{q5x&e1@t
znYH{2ZBpA(Ja_$#+qrR&{pQ`pdrsAi4R?)lfEp92X~OTv*T+mNooLz8^6%z5B}fsc
zD^7td)|9Pg3t|CNZfuu*C%g>mv*O@}cbCYbI^VBTFJaq_#v(zWxB6G3f^%yd-&|wW
zcTG0yi5=u=C6P&5?b$G^R(;Xh4K*y8GQK|57a@1k_A~E;=6T90xV~b_<Oc7(n@%OJ
z5b+q#hB%O3#9t;d_1<ryxNP#X4(nFbW_9ws<G<Z7TD-hpa~w~34-Glr8|4gUb6*VV
z^gMqC6aX~drygZrG*5Yr#4-C{xQjm#);W8E+-C=WS{kzVh#ij2=FDL7>)mGhU<;>g
z>e*EzF;lPvl~Zc9JM>#E4C``LtMwmT3m_-kcU#6Ff#0A|sU0w{d0ME~(X{914Uwr}
z{1zmMz{59-l-~L5&DazgdnQ>ablRK8y#47P)T`csrt6P^y?T%@fx{&AIJ1iuSt-ZM
zxU=z-0>X79S#C9=7Yp^2l?_XlqoXd#iGRlO_uHGW1wISg7LQh~H}}{zUEc2NpwgCQ
zpWKb;LJC0DD2OEO{eo<R#W>Kl8Fe(IaO>8j!W~(M9qpr#!px_2Z?A@;05@H|u~(|(
z9S&+{ckxN%9eWgazmCkZ+B|O+=eC*Rhi;G1%W{8o?YD!CalOI=eE6+~c=*oYh%aFE
z(M0s7Oi*BDh39;n=`D0vz4)%AW-pweP1^R=l73Pl(dkVoEpr(2?rB>==q|_bQ`zu+
z{!aGk?8RB4nAS6aC3uT#!L=k#(2zDqnq$<b9S6-xc9+x;D%mVO7EaYwWrJF-tY=+?
z1Y5p)%%866HhqcNmXvM3P5g49S@b*}GshkE{KTcjPrCV;_%_ML{k~xxi=Q^KgXZ#e
zs(aW#q&7N;;WC!x(DA-N?xn*sg<XPH(n0;OyU0v5X8d!)uKCCg$}Q`zCp3TR)3&1Q
z>GO8{<j=lXi_B2T)JQY{=IsDWY1laS4d)snyn}FyylXss1Puh(Vt=M9RDKm@&mxO{
zjP4<{wy8vlW<hdW&SHY`C^sk<2Ect_8$1l7jR;cSq|KYS#pX^G2n!5n#6YQ}hKGe>
z+&ZwBV0K_{35I>fc%T^M>pF?7`Er(M%wjf)`gL4%jVsM~3S`mfnk^E<Xu&xJC?4Z7
zMrvc+OYD&<$`_xhn5DYB1O4FM!&C&|--@Mnhs8#kW1wUvV38m!)Cr3r`ZzY6$W|`d
zsvdXoL%(1*sAY?=8L%g+7=y-)?LohWi1^ocUD#)Ae`;tBcQ%jXeLe7-l`BC#`x;z{
zt!p4$p0b4Mw1p)bljLISTsI7WDTlr#-gOdQQ+CD+GVPpVyc{bEjJQWT#b<HAWG1-9
zNUcUgV0j3xaVvJ>o#U*jvV5NKqZr93`6yO6Lw(u75{oWRUqS|@-Xm11#R(!D3=ZoB
zX|xid_$Ra!=cfEbQgpRspXs`W<b8BeRKDZrZPxjVwajo-1W#ATkVl7xc71tB{}uK^
zmP<X11ex28webxjh^Yh@nct3i`WamT3x)sviRo0O%TQ?*t&I^+5Tl*(sBprMKxq^B
z4`Bn%PwS6-soa_cx-T~VLzW0xPb|hUuFpg$P!x?8Lmn5(xuA~Ofd|cT16Fq*o8ch{
zP$_R<q^Iov7Ie?JmfbDSbZQp18*Ygaz}Xk+p5~-hp6#TwS$Tr?D6qza9-7jmdxdr%
z+(kY}+D)?KeE{xq2)7_x6Y1(jd!mA_;aBpa;dZ39l5a_TUm{+2|2NP!zb-+vLLxt|
zHQ%ldwA0uylgOH4<kT8G&0c-eQAl^{55_9?X<^tkdch7WPJ29eCINUx+s8J1=LNQt
z6ZO}Zl#L07a*Ao?#>?t_$LuuM<cW82zk>TV<ZavmKblishr&``%P{wzQGF<ff}%Hp
zjJoT71OtXrh36m7GCMt*gh4U?U$W;vb=xkRz&gL@aon_nFZ=23`lX85+Nu$e4uqqg
zx`J}ZhwD8Rzbl+qe9mc#{7^yzr*yHwe$C1_=-c|KLdCBj+p8<@=%3?+B&5<5BMPRl
z_uBKfgM1L&j(qY<N3LNvdy=%m2r)YX0FO~cseaPV_Ri$e=yXvL*Oyc}T?PH$0naNa
zlO1`L>{OV}-jUGfhHZ{TWfk#%+6TuwRVeD%5$==KNPSeY_V{A@HPxALiOjORZfS-A
z-=U-oWdBDbcutQ7MUUiF7ZCw#RDtu%jo7!@MABcyCOBKu&G%0i!?M*Clse^^A6VP^
zBXf<XHJx=T`JMrlGFYMwXt=29n0Y-JxJa{~Q+v<@%by+5Z&e(AxtSoKc@M=2BYGpl
zMUiTn`iDGfIETDMr{TCX(ZK|vM|F1_Lm}1G@D~2#2-TNAIyNOF3+MAQh&Pu--!FS^
zyWtdA#hPZ^s#0T3M2@o=k*n8^7ypvqo~tvRh&yY<T5+YozEP?)|GaI20AU4i*f@iG
zNwpmr=S9sr87!*X>d`pcW2uOcj$lc~+gmz_M910YLCAT?G-tFKJ*zZw!nspp60p&E
z8V!VFyZv}HewGc<ae`SKs64&F&%u|RyGg*!`MS6JJFWJBs(EIv`yV*BXqa%q;sZOk
zF|d59f2s(n(~gp3$%1`<>x{CtX2h?gPtn>%L>~v^cX2-sy@{OhOVVgoDTPqx%ohl0
z(<}0mYg^_q8;KwtS*cZkJzoBmRMwiT)xnj%0AHbqr9n>ye^>T#m#vm$5G-w1;ym;u
z8Pv^`wDuqU{;z5OjC2F)u9K+!PdR)3xyh^R-x5haF5PFj1c~DU6S1oUE(;K46qAGx
zakmcC960Vsb}G%v=ZYExDiMm6>-0oW7`*E3qk^ZZ=Les$$wvkEpOW_}{@=JEfqJo?
zM4k-vuI#-ow#~^zu(Wjvk!7`=6L=b{)PT@`pTqy-Vbt~4@+a7(+xg$1p(cgrO)@dS
z`Qk^wDb&+*T2Tj*s=?vcwOS*XU4P?tbhWae?^Z+-)bpLF&Ny}#d38wkz;$riK7TWN
za9RNGcwux1pEd4nmRzE%$c`A>)#A2Udyts95}m(|WvI140sMyeFB=RiiO1C^cE-1>
zX*Coe#_>s#KvXSpekRh^O{}QJQqu!@7jd8`SHKN^Nv$bUL@3*(%up0D@&y_o>}M+3
z9)bn`$=~(jDQH7yFO{~m_GHU}G)^QE{mb+fuQn5&P`WrH)w`Dy3re`~YSyTLUX5$S
z4&-&@7i-EKLvHe<kII3{;YI}P13Y4-eszUHk8JM;+gTqkkN{3X@XO?;PKmyvhpQz1
zN__0z@rmUG{6Kw-JRVI>r=G#o$i_sCxt79rL3BbUc<n4P$DOwfV-yFb;k)57kq9lz
z!LM*owD(QIw;#X3uWz5O?0A=R1&+0|+0W(*_HPHQ3coU_9HHCr8o%q0(yP)-<7t#H
zON}PT#KG#@V{ulyTE)fK#<g~8+*uKY2*2itV8P9LO)eatnaIT(^OD+{&A2U(E;A13
z^gF&pFC)Q}o(r7um1M@>AY4yplC&(4&hw$urFaEMf8tnZ`qjR_UH`4U^A#ee^GRf<
zKTyno3^+-={=A9en^B>p_^fU|{e!Z>qZ7NOplQau_m>YJ7>6S{y0urmM!GN0-pFA%
zi~28DXl(Z55PcT$!3!qOqdw9!(^IW3{rqVS0%JC$%rXVEQCtQrQv{SzMtQUC4J}dQ
zK+HzRbBpG8U;R0tz`BxzRhW3t`Zd|enY{opOt=tx$xa}Zd?)T@6hlDIwG4YucJH}7
z<71+q$gxbndu6B*>A>NTk4I#`$+DDdSq*L3*J+%vXNSX{k)_d$IF%;yr(<c&@WknK
zj;sCqyztQk7VU3%cp11nY^<67RylNiTaGk(tT(i2Biq&nz&>}eh=$F)VQfMAX*K?|
zH1rvaS~y^Qq9ix5vpR8PR&}P^E(-RW&Mn(L<91&%)#s9&SZ0|<{j&H}^XvMju2&yd
zIAg4dSM>3L25)Z#QFe&|#`p7mq>Z1P7rOJraAu~b7)nSp3V$@sMcpa4hxc!YPj%((
zMIJzm_6dxH9})p4qfeHV8W4q@^Oq0y7kev9s`E$st@bc(FWzr|0dL)`Qw0Lg>4KOi
z0*lRrcw6I$D55lG0Y54tnclE35Pj9CP~=D-%*}4{PgNQIPN`+|x}vFmq<vnG-vg`R
z%tmuKuqFK&d=S}~nVvKQ)oW24=E=xCJEyZqv}z@y$t&g8c6{IHa{jEk7?LoAXE3Bf
z1m<A_%VbHa3+Q7SaTJ}6zDX&263?;eVGB0^#);qwA0l@q@hbgIYCdX+n~UavtCvgD
z6@7i!)b_EBp!N>dFBS78UbhW#c)-ZAx)rU}?T__bf(Bn7t}4u!V7y~;@_e!4Mqp$J
zjBZ6@R~8G>@O&|#7#4#!YBF(h8rh`Vl6eXP-z-!&#vI^9*&vy|wANCLkz<KMOz}&m
z(5mLD)GWzJ!0b42gL=f+T-!a*TkF+ihuF3<gbW-Pl-lK8gBg%)&Q)e~4JOAFvVYF&
z+IlrBz!hMpaIY~!u8Ak>l)YKAIz!MhXeB^$3TRJGPX46iCPt{=a(M2%N9jYw{bg9_
zgI=i<2RziNys-$uR9gefu0Cz)jB;Av-&vpBYj?6k9TZO!>BDV1pZ$$a6M+})KBl=T
zFuY^SRlIxH7`>vFt&yfu=F!XG5r@5UtR!}pbshcl;ZPfyW&2&Fl+~g)3wkA!e<7Uj
zJsbgV*&j=IaV9!v2@||4ZeCx_d9(dqTvhc0>#R+5+yvsh-WkWF7#Ax>ql695;bh(N
z+O;Y-P5%wPM5|ztAS%Y8I$g!NdCe~x%x2pW{^vkO?d#o)qpGGxRj;~8WY;Uhf_O0B
zHeR$1i71IOIMaRs*)LcPbIs9=kvFlI5|?;n*z6)k1m)i9`<7w$WiEr6x)nw2y$pe(
Qk*5f^N^hIGWHhw@0x1I|cmMzZ

literal 0
HcmV?d00001