---
title: "briskaR Demonstrations"
author: "Jean-françois Rey and Samuel Soubeyrand"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{briskaR Demonstration}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r, echo=F}
library("briskaR")
```

## Landscape

> __Commands to load and display a real landscape__
  
```{r,echo=TRUE,fig.width=7,fig.height=7,eval=FALSE}
## load the maize_65 data set incorporated in briskaR
## in which polygons are either maize or non-maize fields
data(maize_65)

## then, randomly assign the GMO and non-GMO marks to maize fields
## non-maize fields being considered as receptors
initial_types=maize_65@thelandscape@data
nb_maize=sum(initial_types[,1])
nb_non_maize=sum(initial_types[,2])
types=as.data.frame(cbind(0,NA,initial_types[,2]))
names(types)=c("sources","neutral","receptors")
GM_fields=sample((1:nrow(types))[initial_types[,1]==1],
0.4*(nb_maize+nb_non_maize),replace=FALSE)
types[GM_fields,1]=1
types[,2]=1*(types[,1]==0 & types[,3]==0)

## generate land1 with the new marks for polygons
land1=maize_65
land1@thelandscape@data=types

## transform receptors at the border of the study domain
## into neutral cells
border_size=200
for(i in 1:land1@n){
  if(land1@thelandscape@data[i,3]==1){
    x0=land1@thelandscape@polygons[[i]]@Polygons[[1]]@coords
    if(sum(x0[,1]<land1@xmin+border_size | x0[,1]>land1@xmax-border_size |
      x0[,2]<land1@ymin+border_size | x0[,2]>land1@ymax-border_size)>0){
      land1@thelandscape@data[i,3]=0
      land1@thelandscape@data[i,2]=1
    }
  }
}
## display the landscape
plot(land1)
```
```{r,fig.width=8,fig.height=8,echo=FALSE}
data(maize_65)
plot(maize.landscape)
```

>__Commands to generate and display a simulated landscape:__
```{r,fig.width=7,fig.height=7}
## generate the Voronoi tesselation
land=simulateInitialPartition(n=200, prop=0.4, range=10,
xmin=0, xmax=5000, ymin=0, ymax=5000)
## generate receptors at the borders of Voronoi cells
land2=simulateThickMargins(land, border_size=200, prob=0.5,
mean_thickness=5, v_thickness=50)
## display the landscape
plot(land2)
```

## Individuals

  > __Commands to generate and display exposed individuals for the real landscape `land1`:__
  
```{r,eval=FALSE}
## tuning parameters and variables
nb_ind=100 ; time_min=1 ; time_max=61
birth_dates=sample(time_min:time_max, size=nb_ind, replace=TRUE)
life_expectancies=rep(20, nb_ind)
toxic_gap=rep(15, nb_ind)
## generate exposed individuals
ind=simulateIndividuals(land1, n=nb_ind, mintime=time_min,
                        maxtime=time_max, dob=birth_dates,
            life_duration=life_expectancies, toxic_threshold=toxic_gap)
## plot individuals in the landscape
plot(land1,ind)
```
```{r,fig.width=8,fig.height=8,echo=FALSE}
plot(maize.landscape,maize.individuals)
```

## Toxic Intensity

  >__Commands to compute the spatio--temporal quantities of GM pollen present on host plant leaves across the real landscape `land1` and the study period:__


```{r,fig.width=7,fig.height=7,eval=FALSE}
## generate emission functions for every sources
list_sources=getSPSources(land1)
emitted_pollen=create.pollen.sources(length(list_sources),time_max)
emitted_pollen=data.frame(t=emitted_pollen,
row.names=row.names(list_sources))
## load precipitation data and restrict them to June and July 2013
data(Precipitation)
precip=Precipitation$RR[Precipitation$AN==2013 & (Precipitation$MOIS==6 | Precipitation$MOIS==7)]
plot(1:61, precip, pch=19, xlab="Day", ylab="Rainfall intensity (in mm)")
```
```{r,fig.width=7,fig.height=7,echo=FALSE}
data(Precipitation)
precip=Precipitation$RR[Precipitation$AN==2013 & (Precipitation$MOIS==6 | Precipitation$MOIS==7)]
plot(1:61, precip, pch=19, xlab="Day", ylab="Rainfall intensity (in mm)")
```

```{r,eval=FALSE}
## compute concentration of pollen across space and time
tox=toxicIntensity(objectL=land1, toxic_emission=emitted_pollen,
                mintime=time_min, maxtime=time_max, size_raster=2^10,
            alpha=list(minalpha=0.1, maxalpha=0.95, covariate_threshold=30,
                                          simulate=FALSE, covariate=precip))
## display concentration of pollen at day 61
plot(x=land1, time=61, objectT=tox)
```
```{r,fig.width=7,fig.height=7,echo=FALSE,eval=TRUE}
tox<-array(0,c(61,1024,1024))
tox[61,,]<-as.matrix(maize.toxicintensity61)
attr(tox,"class") <- "ToxicIntensityRaster"
plot(x=maize.landscape, time=61, objectT=tox)
```

## EcoToxicological

  >__Commands to compute and display the internal concentration of Bt toxin and its impact for every host:__

```{r,eval=FALSE}
ind_impact=ecoToxic(land1, ind, tox, time_min, time_max,kin=4.19e-7, kout=0.1, deltat=0.1)
```

  >Displaying the cumulative impact of the Bt toxin at time 61 can be done as follows:
```{r,fig.width=7,fig.height=7,eval=FALSE}
#plot(x=land1, y=ind_impact, time=30, objectT=tox)
plot(x=land1, y=ind_impact, time=time_max, objectT=tox)
```
```{r,fig.width=7,fig.height=7,echo=TRUE}
plot(maize.landscape,maize.individuals,time=61,objectT=tox)
```

  >The histogram of maximum internal concentration is obtained as follows:
```{r,fig.width=7,fig.height=7,eval=FALSE}
temp=apply(ind_impact@intern_toxic,1,max)
hist(temp, main="Real landscape", xlab="Maximum internal concentration",
breaks=seq(0,max(temp)+5,5))
abline(v=15,col=3)
```
```{r,fig.width=8,fig.height=8,echo=FALSE}
temp=apply(maize.individuals@intern_toxic,1,max)
hist(temp, main="Real landscape", xlab="Maximum internal concentration",
breaks=seq(0,max(temp)+5,5))
abline(v=15,col=3)
```

  >Displaying the time series of environmental and internal concentration of toxins can be done as follows:
```{r,fig.width=7,fig.height=7,eval=FALSE}
ID=21
plot(x=land1, y=ind_impact, objectT=tox, numind=ID)
```
![time series](toxic21.png)

