---
title: "Optimizer"
author: "Maximilian Pichler"
date: "15 6 2020"
output: 
  html_document: 
    keep_md: yes
    toc: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
set.seed(42)
library(sjSDM)
library(dplyr)

communities = lapply(c(25, 50, 100), function(sp) lapply(1:10, function(i) simulate_SDM(env = 3L, sites = 100L, species = sp)))

optims = c("adamax", "rmsprop", "accsgd", "adabound", "sgd", "diffgrad")
scheduler = c(TRUE, FALSE)
lr = c(0.01, 0.005, 0.002, 0.001)
epochs = c(100, 200)
bs = c(5, 10, 30)
sp = c(25, 50, 100)
iter = 1:10


test = data.frame(expand.grid(optims, scheduler, lr, epochs, bs, sp, iter))
colnames(test) = c("optim", "scheduler","lr","epochs", "bs", "sp", "iter")
test$acc = NA
test$grads = NA
test$rmse = NA
test = test[order(test$sp),]
```


```{r}
n = 12
cuts = cut(1:nrow(test), breaks = n)
levels(cuts) = 1:n

tests = lapply(1:n, function(i) test[cuts==paste0(i),])

library(snow)
cl = snow::makeCluster(n)
nodes = unlist(snow::clusterEvalQ(cl, paste(Sys.info()[['nodename']], Sys.getpid(), sep='-')))
control = snow::clusterEvalQ(cl, {library(sjSDM)})
snow::clusterExport(cl, list("communities", "test", "tests", "nodes"))

results = parLapply(cl, 1:n, function(n) {
  tmp = tests[[n]]
  for(i in 1:nrow(tmp)) {
    if(tmp[i,]$sp == 25) j = 1
    else if(tmp[i,]$sp == 50) j = 2
    else j = 3
    com = communities[[j]][[tmp[i,]$iter]]
    
    opt = 
    switch (as.character(tmp[i,]$optim),
      adamax = Adamax(),
      rmsprop = RMSprop(),
      accsgd = AccSGD(),
      adabound = AdaBound(),
      sgd = SGD(),
      diffgrad = DiffGrad()
    )
    
    if(n %in% 1:4) device = 0
    else if(n %in% 5:8) device = 1
    else if(n %in% 9:12) device = 2
    
    
    m = sjSDM(com$response, com$env_weights, 
              learning_rate = tmp[i,]$lr, 
              step_size = tmp[i,]$bs, 
              iter = tmp[i,]$epochs, 
              control = sjSDMControl(opt, tmp[i,]$scheduler),
              device = device)
    
    tmp[i,]$acc = com$corr_acc(getCov(m))
    tmp[i,]$grads = m$model$params[[2]][[1]]$grad$pow(2.0)$mean()$sqrt()$item()
    tmp[i,]$rmse = sqrt(mean( ( t(coef(m)[[1]])[-1,] - com$species_weights )^2 ))
    rm(m)
  }
  return(tmp)
  
})

```

```{r}
res = do.call(rbind, results)
summary(lm(acc~as.factor(sp) + lr + epochs + optim, data = res))

res2 = 
  res %>% 
    group_by(optim, scheduler, lr, epochs, bs, sp) %>% 
    summarise(acc = mean(acc), grads = mean(grads), rmse = mean(rmse))

```


## Generally
```{r}
par(mfrow=c(1,2))
boxplot(res2$acc~res2$optim, ylim = c(0.5, 1.0))
boxplot(res2$rmse~res2$optim, ylim = c(0.0, 0.6))
```
RMSprop and AdaBound are the most stable, Adamax the worst


## Learning rate
```{r}
par(mfrow=c(1,3))
for(i in c(0.01, 0.005, 0.001)) boxplot(res2$acc[res2$lr==i]~res2$optim[res2$lr==i], ylim = c(0.5, 1.0))
```
AdaBound is influenced the least by the learning rate!


## Epochs
```{r}
par(mfrow=c(1,2))
for(i in c(100,200)) boxplot(res2$acc[res2$epochs==i]~res2$optim[res2$epochs==i], ylim = c(0.5, 1.0))
```

## Batch size
```{r}
par(mfrow=c(1,3))
for(i in bs) boxplot(res2$acc[res2$bs==i]~res2$optim[res2$bs==i], ylim = c(0.5, 1.0))
```
RMSprop and AdaBound are insensitive against the batch size


## With or without scheduler:
```{r}
par(mfrow=c(1,2))
boxplot(res2$acc[ !res2$scheduler]~res2$optim[!res2$scheduler], ylim = c(0.5, 1.0))
boxplot(res2$acc[ res2$scheduler]~res2$optim[res2$scheduler], ylim = c(0.5, 1.0))
```


## RMSprop vs AdaBound
```{r}
summary(lm(acc~., data=res[res$optim=="rmsprop",-c(1,9,10)]))
summary(lm(acc~., data=res[res$optim=="adabound",-c(1,9,10)]))
summary(lm(acc~., data=res[res$optim=="adamax",-c(1,9,10)]))
```

