---
output: html_document
runtime: shiny
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
	echo = FALSE,
	message = FALSE,
	warning = FALSE
)
library(shiny)
library(RGeostats)
library(ranger)
library(e1071)

rm(list=ls())
constant.define("asp",1)
default.nd = 100
default.seed = 123456
default.rot = 45
default.diffmax = 5
default.custom1 = "cos(db$x1/20) - sin(db$x2/20)"
default.custom2 = "cos(db$x1/20 - db$x2/20)"
msg.reset = "Click 'Reset' for reseting the input database"
msg.calculate = "Click 'Calculate' for updating the results"
msg.up.to.date = "Results are 'up-to-date'"
```

<style type="text/css">
.main-container {
  max-width: 100%;
  margin-left: auto;
  margin-right: auto;
}
</style>

```{r Initialization phase}

st.init.db <- function(nd=default.nd,
                       theta=default.rot,
                       custom1=default.custom1,
                       custom2=default.custom2,
                       nbtuba=1000)
{
  set.seed(default.seed)
  rot = util.ang2mat(2, theta)
  theta = theta * pi / 180
  db = db.create(nx = c(nd, nd))
  
  #Ingredients
  mod = model.create("Cubic", range = c(30,30))
  db = simtub(NA, db, mod, nbsimu = 1, nbtuba = nbtuba)
  db = db.rename(db, names = db$natt, newnames = "y1")
  mod = model.create("Spherical", range = c(30,60), aniso.rotmat = rot)
  db = simtub(NA, db, mod, nbsimu = 1, nbtuba = nbtuba)
  db = db.rename(db, names = db$natt, newnames = "y2")
  
  # Anisotropic linear variable
  mod = model.create("Linear", range = c(5,50), aniso.rotmat = rot)
  db = simtub(NA, db, mod, nbsimu = 1, nbtuba = nbtuba)
  db = db.rename(db, names = db$natt, newnames = "V.Aniso.Linear")
  
  # Anisotropic multi-structured variable 1
  w = apply(rbind(db$y1, db$y2), 2, function(a) 3*a[1] + 2*a[2])
  db = db.add(db, V.Aniso.Multi.1 = w)
  
  # Anisotropic multi-structured variable 2
  w = apply(rbind(db$y1, db$y2), 2, function(a) -a[1] + 5*a[2])
  db = db.add(db, V.Aniso.Multi.2 = w)
  
  # Circle arcs
  w = apply(rbind(db$x1, db$x2), 2, function(a) sqrt(sum(a^2))*cos(sqrt(sum(a^2))*20/nd)) / nd
  db = db.add(db, V.Circles = w)
  
  # Along theta trend
  if (tan(theta-pi/2) > 10)
    w = db$x1 / nd
  else if (tan(theta-pi/2) < -10)
    w = (nd - db$x1) / nd
  else
    w = apply(rbind(db$x1, db$x2), 2, function(a) a[2] - tan(theta-pi/2)*a[1]) / nd
  db = db.add(db, V.Trend.Along = w)
  
  # Perpendicular to theta trend
  if (tan(theta) > 10)
    w = db$x1 / nd
  else if (tan(theta) < -10)
    w = (nd - db$x1) / nd
  else
    w = apply(rbind(db$x1, db$x2), 2, function(a) a[2] - tan(theta)*a[1]) / nd
  db = db.add(db, V.Trend.Across = w)
  
  # Customized Variable 1
  w = eval(parse(text=custom1))
  db = db.add(db, V.Custom.1 = w)
  
  # Customized Variable 2
  w = eval(parse(text=custom2))
  db = db.add(db, V.Custom.2 = w)

  # Cleaning ingredients
  db = db.delete(db,"y*")
  
  # Return the created db
  db
}

st.init <- function(env,input)
{
  env$message = msg.calculate
  env$db0 = st.init.db(nd=input$nd,
                       theta=input$rotation,
                       custom1=input$custom1,
                       custom2=input$custom2)
  env
}
```

```{r Some functions}
st.eval.debug <- function(env,input)
{
  if (! miss(input$ligne))
    eval(parse(text = input$ligne))
}

st.update.msg <- function(env,reset=TRUE)
{
  if (reset)
    env$message = msg.reset
  else if (env$message != msg.reset)
    env$message = msg.calculate
  env
}

st.print <- function(env)
{
  if (! miss(env$message))
    cat(env$message)
}

st.check <- function(env,input)
{
  env$message = NA
  if (length(input$auxvars) > 2)
  {
    env$message <- "Error: Maximum 2 covariables authorized"
  }
  if (input$mainvar %in% input$auxvars)
  {
    env$message <- "Error: Main variable cannot be a covariable"
 } 
  env
}

st.calculate.train <- function(input,db)
{
  ns = db$nsamples
  if (input$regular == 2)
  {
    by = round(1/input$ratio, 0)
    tidx = seq(0, ns, by)
  }
  else if (input$regular == 3)
  {
    by = round(sqrt(1/input$ratio), 0)
    sx = seq(0, input$nd-1, by)
    tidx = c()
    for (s1 in sx)
      for (s2 in sx)
        tidx = c(tidx, s1*input$nd + s2)
  }
  else # input$regular = 1
  {
    set.seed(input$seed)
    nt = ns * input$ratio
    tidx = sample(ns, nt)
  }
  train = rep(FALSE, ns)
  train[tidx] = TRUE
  train
}

local.krige <- function(dbin,dbout,auxnames,radix,verbose=FALSE)
{
  nbaux = length(auxnames)
  uc = c("1")
  
  if (nbaux > 0) 
  {
    dbin = db.locate(dbin,auxnames,"f")
    dbout = db.locate(dbout,auxnames,"f")
    if (nbaux == 1) uc = c("1","f1")
    if (nbaux == 2) uc = c("1","f1","f2")
  }
  
  vario = vario.calc(dbin, dirvect = c(0,45,90,135), lag = 2, tolang = 5, nlag = 20, uc=uc)
  model = model.auto(vario, struct = c(1,2,3,5,12), draw = FALSE, tolsigma = 10)
  if (nbaux == 0)
    krig.message = paste("KO for", db.getname(dbin, "z", NA), "using best model:", model$summary, "...\n")
  else
    krig.message = paste("KED for", db.getname(dbin, "z", NA), "with", paste(db.getname(dbin,"f",NA),collapse=' '), 
                    "using best model:", model$summary, "...\n")
  if (verbose) cat(krig.message)

  neigh = neigh.create(ndim = 2, type = 2, nmaxi = 20, radius = 100)
  dbout = kriging(dbin, dbout, model, neigh, flag.std = FALSE, uc=uc)
  
  list(dbout=dbout, vario=vario, model=model)
}

st.calculate <- function(env,input)
{
  if (env$message == msg.reset)
    st.init(env,input)
  
  aux.ranks = as.integer(input$auxvars)
  main.rank = as.integer(input$mainvar)
  env = st.check(env,input)
  
  if (! miss(env$db0) && miss(env$message))
  {
    db         = env$db0
    db         = db.locerase(db,"z")
    vars       = db.ident(db,"V*")
    db         = db.locate(db,vars[main.rank],"z")
    auxnames   = vars[aux.ranks]
    main.var   = vars[main.rank]
    
    if (input$lognormal)
      db[,main.var] = exp(db[,main.var] / 5)
    
    # Input dataset extraction
    dbi = db.add(db, st.calculate.train(input, db), locnames = "sel")

    # Random forest
    db = ml.rf(dbi,db,auxnames=auxnames,radix="rf",
               num.trees=c(100,500,1000,2000),
               mtry=seq(1,length(auxnames)+2),  # X1, X2, V, AUXs
               replace=input$rf_replace,
               sample.fraction=ifelse(input$rf_replace,1,input$rf_ratio),
               verbose=TRUE)
    db = db.rename(db, names = "*rf", newnames = "rf")

    # SVM
    db = ml.svm(dbi,db,auxnames=auxnames,radix="svm",
                kernel=c("radial", "linear", "sigmoid", "polynomial"),
                cost=c(1,10,100),
                verbose=TRUE)
    db = db.rename(db, names = "*svm", newnames = "svm")
    
    # Ordinary co-kriging
    kres = local.krige(dbi,db,auxnames=auxnames,radix="krig",verbose=TRUE)
    db = db.rename(kres$dbout, names = "Kriging*", newnames = "krig")
    
    env$message = msg.up.to.date
    env$dbres = db
    env$model = kres$model
    env$vario = kres$vario
  }
  env
}

```


```{r Display functions}
st.display <- function(env,input,row,col)
{
  db0       = env$db0
  dbres     = env$dbres
  aux.ranks = as.integer(input$auxvars)
  main.rank = as.integer(input$mainvar)
  nbaux     = length(aux.ranks)

  if (! miss(db0))
  {
    vars  = db.ident(db0,"V*")
    main.var = vars[main.rank]
    if (input$lognormal)
    {
      main.val = exp(db0[,main.var] / 5)
    }
    else
    {
      main.val = db0[,main.var]
    }
    zlim = range(main.val)
  }
  
  if (! miss(db0))
  {
    if (row == 1)
    {
      if (input$training)
      {
        sel.train = st.calculate.train(input,db0)
        dbs = db.add(db0, s = sel.train, locnames = "sel")
      }
      
      if (col == 1)
      {
        nm = db.getname(db0,NA,main.var)
        zz = paste(round(zlim,2), collapse=',')
        title = paste0(nm, " [", zz , "]")
        if (input$lognormal)
        {
          dbl = db.add(db0, main.val)
          plot(dbl, title = title, zlim = zlim)
        }
        else
        {
          plot(db0, title = title, name = main.var, zlim = zlim)
        }
        if (input$training)
          plot(dbs, name="s", flag.aspoint=TRUE, add=TRUE)
      }
      if (col == 2 && nbaux > 0)
      {
        plot(db0,name=vars[aux.ranks[1]])
      }
      if (col == 3 && nbaux > 1)
      {
        plot(db0,name=vars[aux.ranks[2]])
      }
    }
  }
  if (! miss(dbres) && env$message == msg.up.to.date)
  {
    if (row == 2)
    {
      if (col == 1)
        plot(dbres,name="rf",   title="Random Forest (RF)",               zlim = zlim)
      if (col == 2)
        plot(dbres,name="svm",  title="Support Vector Machine (SVM)",     zlim = zlim)
      if (col == 3)
      {
        if (nbaux > 0)
          plot(dbres,name="krig", title="Kriging with Ext. Drifts (KED)", zlim = zlim)
        else
          plot(dbres,name="krig", title="Ordinary Kriging (KO)",          zlim = zlim)
      }
        
    }
    else if (row == 3 && input$match)
    {
      if (col == 1)
      {
        d = abs(dbres$rf-main.val)^2
        db = db.add(dbres, d)
        plot(db, title = paste("Diff RF (m=",round(mean(db$d),4),")",sep=""), zlim = c(0, input$diffmax), col = heat.colors(50))
      }
      if (col == 2)
      {
        d = abs(dbres$svm-main.val)^2
        db = db.add(dbres, d)
        plot(db, title = paste("Diff SVM (m=",round(mean(db$d),4),")",sep=""), zlim = c(0, input$diffmax), col = heat.colors(50))
      }
      if (col == 3)
      {
        d = abs(dbres$krig-main.val)^2
        db = db.add(dbres, d)
        if (nbaux > 0)
          plot(db, title = paste("Diff KED (m=",round(mean(db$d),4),")",sep=""), zlim = c(0, input$diffmax), col = heat.colors(50))
        else
          plot(db, title = paste("Diff KO (m=",round(mean(db$d),4),")",sep=""), zlim = c(0, input$diffmax), col = heat.colors(50))
      }
    }
  }
}
```


```{r}
ui <- shinyUI(
  navbarPage("", id="navbar",
             tabPanel("RF vs SVM vs KED",
                      sidebarLayout(
                        sidebarPanel(
                          navbarPage("",id="navbar2",
                                     tabPanel("Reset",
                                              numericInput("nd","Grid Size (NX or NY)",default.nd),
                                              numericInput("rotation","Anisotropic Rotation (°)",default.rot),
                                              textInput("custom1", label = "Customized Input Variable#1",
                                                        value = default.custom1,
                                                        placeholder="Enter a valid equation"),
                                              textInput("custom2", label = "Customized Input Variable#2",
                                                        value = default.custom2,
                                                        placeholder="Enter a valid equation"),
                                              # Toggle "display data"
                                              actionButton("reset","Reset")
                                     ),
                                     tabPanel("Calculate",
                                              sliderInput("ratio","Sampling Ratio",
                                                          min=0.001,max=0.1,value=0.025),
                                              numericInput("seed","Sampling Seed",default.seed),
                                              radioButtons("regular", label = "Sampling mode",
                                                           choices = list("Irregular"  = 1,
                                                                          "1D Regular" = 2,
                                                                          "2D Regular" = 3), 
                                                           selected = 1),
                                              selectInput("mainvar", label = "Main Variable", 
                                                          choices = list("Aniso.Linear"  = 1, 
                                                                         "Aniso.Multi.1" = 2, 
                                                                         "Aniso.Multi.2" = 3,
                                                                         "Circles"       = 4,
                                                                         "Trend.Along"   = 5,
                                                                         "Trend.Across"  = 6,
                                                                         "Custom.Var.1"  = 7,
                                                                         "Custom.Var.2"  = 8),
                                                          selected = 1),
                                              checkboxInput("lognormal", label = "Lognormal Main Variable",
                                                            value = FALSE),
                                              checkboxGroupInput("auxvars", label = "Co-variables", 
                                                                 choices = list("Aniso.Linear"  = 1, 
                                                                                "Aniso.Multi.1" = 2, 
                                                                                "Aniso.Multi.2" = 3,
                                                                                "Circles"       = 4,
                                                                                "Trend.Along"   = 5,
                                                                                "Trend.Across"  = 6,
                                                                                "Custom.Var.1"  = 7,
                                                                                "Custom.Var.2"  = 8),
                                                                 selected = 0),
                                              actionButton("calculate","Calculate")
                                     ),
                                     tabPanel("Display",
                                              checkboxInput("training", label = "Show Training Data",
                                                            value = TRUE),
                                              checkboxInput("match", label = "Show Matching Criteria",
                                                            value = TRUE),
                                              numericInput("diffmax","Matching Criteria Maximum Value (white)",
                                                           default.diffmax)
                                     ),
                                     tabPanel("Random Forest",
                                              checkboxInput("rf_replace", label = "Replace",
                                                            value = TRUE),
                                              conditionalPanel(
                                                condition = "input.rf_replace == false",
                                                sliderInput("rf_ratio","Sampling Ratio",
                                                            min=0.1,max=0.95,value=0.63)
                                              ),
                                              actionButton("calculate2","Calculate")
                                     )
                          )
                        ),
                        mainPanel(
                          fluidPage(
                            fluidRow(
                              column(4, plotOutput(outputId="graph1.1", height="300px"), style='padding:0px; margin-top:-20px;'),
                              column(4, plotOutput(outputId="graph1.2", height="300px"), style='padding:0px; margin-top:-20px;'),
                              column(4, plotOutput(outputId="graph1.3", height="300px"), style='padding:0px; margin-top:-20px;')),
                            fluidRow(
                              column(4, plotOutput(outputId="graph2.1", height="300px"), style='padding:0px; margin-top:-40px;'),
                              column(4, plotOutput(outputId="graph2.2", height="300px"), style='padding:0px; margin-top:-40px;'),
                              column(4, plotOutput(outputId="graph2.3", height="300px"), style='padding:0px; margin-top:-40px;')),
                           fluidRow(
                              column(4, plotOutput(outputId="graph3.1", height="300px"), style='padding:0px; margin-top:-40px;'),
                              column(4, plotOutput(outputId="graph3.2", height="300px"), style='padding:0px; margin-top:-40px;'),
                              column(4, plotOutput(outputId="graph3.3", height="300px"), style='padding:0px; margin-top:-40px;')),
                            verbatimTextOutput("mesarea")
                          )
                        )
                      )
             ),
             tabPanel("Debug",
                      sidebarLayout(
                        sidebarPanel(
                          textInput(inputId = "ligne", 
                                    label = "Enter the Code Line to be executed", 
                                    value = "env$db0") 
                        ),
                        mainPanel(
                          verbatimTextOutput("debug.text"),
                          plotOutput(outputId="debug.plot")
                        )
                      )
             ),
             tabPanel("Quit",
                      icon = icon("circle-o-notch"),
                      value = "quit"
             )
  )
)

server <- shinyServer(function(input, output, session) 
{
  env = reactiveValues(db0 = st.init.db(), dbres = NA,
                       message = msg.calculate,
                       model = NA, vario = NA)

  observeEvent(input$reset,      { env = st.init(env,input) })
  observeEvent(input$calculate,  { env = st.calculate(env,input) })
  observeEvent(input$calculate2, { env = st.calculate(env,input) })
  
  observeEvent(input$nd,         { env = st.update.msg(env,TRUE) })
  observeEvent(input$rotation,   { env = st.update.msg(env,TRUE) })
  observeEvent(input$custom1,    { env = st.update.msg(env,TRUE) })
  observeEvent(input$custom2,    { env = st.update.msg(env,TRUE) })
  
  observeEvent(input$ratio,      { env = st.update.msg(env,FALSE) })
  observeEvent(input$seed,       { env = st.update.msg(env,FALSE) })
  observeEvent(input$regular,    { env = st.update.msg(env,FALSE) })
  observeEvent(input$lognormal,  { env = st.update.msg(env,FALSE) })
  observeEvent(input$mainvar,    { env = st.update.msg(env,FALSE) })
  observeEvent(input$auxvars,    { env = st.update.msg(env,FALSE) })
  
  observeEvent(input$rf_replace, { env = st.update.msg(env,FALSE) })
  observeEvent(input$rf_ratio,   { env = st.update.msg(env,FALSE) })
  
  output$graph1.1  <- renderPlot ({ st.display(env,input,1,1) })
  output$graph1.2  <- renderPlot ({ st.display(env,input,1,2) })
  output$graph1.3  <- renderPlot ({ st.display(env,input,1,3) })
  output$graph2.1  <- renderPlot ({ st.display(env,input,2,1) })
  output$graph2.2  <- renderPlot ({ st.display(env,input,2,2) })
  output$graph2.3  <- renderPlot ({ st.display(env,input,2,3) })
  output$graph3.1  <- renderPlot ({ st.display(env,input,3,1) })
  output$graph3.2  <- renderPlot ({ st.display(env,input,3,2) })
  output$graph3.3  <- renderPlot ({ st.display(env,input,3,3) })
  
  output$mesarea   <- renderPrint({ st.print(env)            })
  output$debug.text<- renderPrint({ st.eval.debug(env,input) })
  output$debug.plot<- renderPlot ({ st.eval.debug(env,input) })

  observe({ if (input$navbar == "quit") stopApp()})
  
})

shinyApp(ui=ui,server=server,options=list(height=1000))
```