Sampling Contamination 3D.Rmd

---
title: "Sampling Contamination 3D"
author: "Xianbin Cheng"
date: "January 29, 2019"
output: html_document
---

# Method #

1. Load in libraries, visualization functions and check the session info.

```{r, warning = FALSE, message= FALSE}
source(file = "Sampling_libraries.R")
source(file = "Sampling_contamination.R")
source(file = "Sampling_visualization.R")
library(plotly)
```

2. Define the input paramaters and create a simulation function.

* `n_contam` = the number of contamination points  
* `x_lim` = the limits of the x-axis  
* `y_lim` = the limits of the y-axis  
* `z_lim` = the limits of the z-axis
* `x` = the horizontal coordinate of the contamination center, which follows a uniform distribution by default
* `y` = the vertical coordinate of the contamination center, which follows a uniform distribution by default
* `z` = the vertical coordinate of the contamination center, which follows a uniform distribution by default
* `cont_level` = a vector that indicates the mean contamination level (logCFU/g or logCFU/mL) and the standard deviation in a log scale, assuming contamination level follows a log normal distribution $ln(cont\_level)$~$N(\mu, \sigma^2)$. 
* `dis_level` = a vector that indicates the mode (ppb) and the lower bound (ppb), assuming contamination level follows $lb+Gamma(\alpha, \theta=\frac {mode-20}{\alpha-1})$

**Discrete Spread**

* `n_affected` = the number of affected kernels near the contamination spot, which follows a Poisson distribution ($Pois(\lambda = 5)$)   
* `covar_mat` = covariance matrix of `x` and `y`, which defines the spread of contamination. Assume the spread follows a 2D normal distribution with var(X) = 0.25, var(Y) = 0.25 and cov(X, Y) = 0  


```{r}
## The input parameters
n_contam = rpois(n = 1, lambda = 3)
x_lim = c(0, 10)
y_lim = c(0, 10)
z_lim = c(0, 10)
lims = list(xlim = x_lim, ylim = y_lim, zlim = z_lim)
cont_level = c(7, 1)
dis_level = c("mode"= 40000, "lb" = 20)
spread = "discrete"

### Discrete
n_affected = rpois(n = 1, lambda = 5)
covar_mat = make_covar_mat(spread = spread, varx = 0.25, vary = 0.25, varz = 0.25, covxy = 0, covxz = 0, covyz = 0)

### Continuous
spread_radius = 1
LOC = 10^(-3)
```

3. Generate contamination.

```{r}
sim_contam_new
```

```{r}
# Continuous
contam_xy1 = sim_contam_new(n_contam = n_contam, lims = list(xlim = c(0,10), ylim = c(0,10)), spread = "continuous", covar = covar_mat, n_affected = n_affected, spread_radius = spread_radius, cont_level = cont_level, dis_level = dis_level)

# Discrete
contam_xy2 = sim_contam_new(n_contam = n_contam, lims = lims, spread = "discrete", covar = covar_mat, n_affected = n_affected, spread_radius = spread_radius, cont_level = cont_level, dis_level = dis_level)
```

# Result

1. Simulated data.

```{r}
kable_styling(kable(contam_xy1, format = "html"), full_width = TRUE)
kable_styling(kable(contam_xy2, format = "html"), full_width = TRUE)
```

2. Visualization.

```{r}
# Continuous spread
contam_draw(data = contam_xy1, spread = "continuous", xlim = lims$xlim, ylim = lims$ylim)
```

```{r, size = "50%"}
# Discrete spread
# 2D projection
contam_draw(data = contam_xy2, spread = "discrete", xlim = lims$xlim, ylim = lims$ylim)

# 3D view
scatter3D(x = contam_xy2$X, y = contam_xy2$Y, z = contam_xy2$Z, type = "p", colvar = as.numeric(contam_xy2$dis_level), 
          xlim = lims$xlim, ylim = lims$ylim, zlim = lims$zlim, phi = 20, theta = 40 )
```