inst/spectraluncertainty-shiny/app.R

In EarthSystemDiagnostics/spectraluncertainty:

# Source external files and modules
library(spectraluncertainty)
library(orbitalforcing)
library(tidyverse)
#library(plotly)
library(shinyWidgets)

# Functions for Shiny app
#' Plot proxy record
#'
#' @param data
#'
#' @export
#'
#' @examples
ggproxyPlot <- function(data) {

  ggplot(data, aes(x = age, y = value, colour = resolution)) +
    geom_line(data = filter(data, resolution == "Smoothed")) +
    geom_ribbon(data = filter(data, resolution == "Smoothed"),
                aes(ymin = value - Total.error, ymax = value + Total.error,
                    fill = "+-1 SE"), colour = "White", alpha = 0.5) +
    geom_point(data = filter(data, resolution == "Raw data"), alpha = 0.5) +
    geom_rug(data = filter(data, resolution == "Raw data"),
             sides = "b", colour = "Grey") +
    scale_fill_manual("Error", values = c("+-1 SE" = "#3B9AB2")) +
    scale_colour_manual("Data", values = c("Raw data" = "Black", "Smoothed" = "Blue"),
                        guide = guide_legend(
                          override.aes = list(linetype = c(NA,1),
                                              shape = c(16, NA)))) +
    theme_bw() +
    theme(legend.position = "right")
}


#' Plot world map with crosshairs.
#'
#' @param lat Latitude
#' @param lon Longitude
#'
#' @return
#' @import maps
#' @export
#'
#' @examples
PlotWorld <- function(lat, lon){
  world <- map_data("world")
  worldmap <- ggplot(world, aes(x = long, y = lat, group = group)) +
    geom_polygon(fill = "Grey") +
    coord_quickmap() +
    scale_x_continuous(limits = c(-180, 180), breaks = seq(-180, 180, 60))+
    scale_y_continuous(limits = c(-90, 90), breaks = c(-45, 0, 45)) +
    geom_vline(xintercept = lon, linetype = 2) +
    geom_hline(yintercept = lat, linetype = 2) +
    labs(x = "Longitude", y = "Latitude") +
    theme_bw()
  worldmap
}


depth.opts <- c(breitkreuz.depth.tbl$depth.upr[1:7])


# Define UI ----
ui <- fluidPage(
  shinyjs::useShinyjs(),
  tags$head(
    tags$style(HTML("
                    .shiny-output-error-validation {
                    color: red;
                    }
                    "))
  ),

  titlePanel("Spectral uncertainty (0.1)"),
  p(tags$div(h2(
      HTML(paste(tags$span(style="color:red", "Warning, we are still in the validation and testing phase."), sep = ""))
    )),
    br()),
  h5("Theory, assumptions and data sources."), br(),
  p("The theory behind the spectral error model and derivation of the mathematical
    expression is documented in Kunz et al., (in preparation)."),
  p("As some error terms depend on the climate signal the spectral error model
    requires some assumptions to be made about the local climate signal:"),
  p("The seasonal cycle of temperature and d18Osw is derived from Breitkreuz et al.,
    JGR (2018) and modulated on orbital time-scales assuming a linear relationship
    to the local orbital forcing."),
  p("The climate spectrum is a composite of the observed sea surface temperature
    spectrum (Laepple and Huybers, GRL 2014) for time-scales faster than 1/30 years
    and a power-law spectrum for slower time-scales."),
  p("Basing the error estimates largely on the modern climate variability is a
    reasonable assumption as the temporal changes in marine non-sea ice covered
    regions throughout the Quaternary are small compared the modern spatial
    difference (e.g. polar vs. tropical settings)."),
  p("Additionally, bioturbation is assumed to follow the Berger and Heath (JMR 1968) model."), br(),
  p(h5("References:"),
HTML("<ul>
         <li>Kunz, Dolman and Laepple.: Timescale-dependent uncertainty of
sediment proxy records. Part I: Theoretical concept. In preparation.</li>
<br/>
         <li>Berger, W. H. and Heath, G. R.: Vertical mixing in pelagic
         sediments, J. Mar. Res., 26(2), 134–143, 1968.</li>
         <li>Breitkreuz, C., Paul, A., Kurahashi‐Nakamura, T., Losch, M. and
         Schulz, M.: A Dynamical Reconstruction of the Global Monthly Mean
         Oxygen Isotopic Composition of Seawater, Journal of Geophysical Research:
         Oceans, 123(10), 7206–7219, doi:10.1029/2018JC014300, 2018.</li>
         <li>Laepple, T. and Huybers, P.: Global and regional variability in
         marine surface temperatures, Geophysical Research Letters, 41(7),
         2528–2534, doi:10.1002/2014GL059345, 2014.</li>
         </ul>")),
  hr(),
  p(
    "Please contact Dr Andrew Dolman <andrew.dolman@awi.de>, Dr Torben Kunz <tkunz@awi.de>,
    or Dr Thomas Laepple <tlaepple@awi.de>, at the Alfred-Wegener-Institute,
    Helmholtz Centre for Polar and Marine Research,
    Germany, for more information.", br(), br(),
    "This work was supported by German Federal Ministry of Education and Research
    (BMBF) as Research for Sustainability initiative",
    a("(FONA)", href = "https://www.fona.de/",  target = "_blank"),
    "through the",
    a("PalMod", href = "https://www.palmod.de/",  target = "_blank"),
    "project (FKZ: 01LP1509C).",
    br(),
    br(),
    a(
      img(src = "PalMod_Logo_RGB.png", align = "top"),
      href = "https://www.palmod.de/",
      target = "_blank"
    )
  ),
  br(),
  sidebarPanel(width = 4,
    tabsetPanel(id = "sidebar",
      tabPanel(
        "Parameters", value = "params",
        # Proxy characteristics ------
        fluidRow(
          h4("Proxy type and characteristics"),
          column(12,
                 radioButtons(
                   "proxy.type", "Proxy type",
                   choiceNames = c("T° (Mg/Ca)", "T° (Uk'37)", "d18O"),
                   choiceValues = c("T_deg_Mg_Ca", "T_deg_Uk37", "d18O")
                 ))),
        fluidRow(
          uiOutput("ui.proxy.chars")
        ),
        # Climate assumptions -----
        fluidRow(
          h4("Climate assumptions"),
          h6("Seasonal cycle"),
          checkboxInput("use.lat", "Use proxy location to determine the seasonal cycle parameters.",
                        value = TRUE),
          #helpText("Use proxy location to determine the seasonal cycle parameters."),
          fluidRow(column(width = 12,
            column(width = 6,
                   sliderInput(
                     "latitude",
                     h6("Latitude of proxy location"),
                     value = 0.5,
                     step = 1,
                     min = -78.5,
                     max = 89.5
                   )),
            column(width = 6,
                   sliderInput(
                     "longitude",
                     h6("Longitude of proxy location"),
                     value = 0.5,
                     step = 1,
                     min = -179.5,
                     max = 179.5
                   )),
            column(width = 6,
                   shinyWidgets::sliderTextInput("depth", h6("Habitat depth range [m]"),
                                                 choices = depth.opts, selected = c(0, -50),
                                                 grid = TRUE,
                                                 from_max = -550, to_min = -50,
                                                 dragRange = TRUE)))),
          fluidRow(column(width = 12,
                          uiOutput("ui.seas.cycle"))),

          fluidRow(column(width = 12,
                          uiOutput("ui.power.spec")))
          ),
        #Core Sampling ------
        # fluidRow(
        #   h4("Core sampling"),
        #   #conditionalPanel(condition = "output.fileUploaded == true", {
        #     h6("Infer temporal parameters from uploaded proxy record."),
        #     checkboxInput("infertauT",
        #                   HTML("Infer &tau;<sub>T</sub> from proxy record?"),
        #                   value = FALSE)
        #     #}
        #     #)
        #),
        # Timescale ----
        fluidRow(
          column(width = 12,
                 uiOutput("ui.infertauT")
          )
        ),
        fluidRow(
          column(width = 12,
            uiOutput("ui.timescale")
            )
        ),
        fluidRow(column(
          width = 6,
          numericInput(
            "slice.thickness",
            h6("Thickness of sediment slices [cm]"),
            value = 1,
            step = 0.1,
            min = 0,
            max = 10
          )
        )),
        # Sedimentation parameters -----
        fluidRow(
          h4("Sedimentation parameters"),
          column(
            6,
            numericInput(
              "bio.depth",
              h6("Bioturbation depth [cm]"),
              value = 10,
              step = 1,
              min = 0,
              max = 30
            )
          ),
          column(
            6,
            numericInput(
              "sed.acc.rate",
              h6("Sediment accumulation rate [cm/ka]"),
              value = 50,
              step = 1,
              min = 0,
              max = 100
            )
          )
        )
      ),
      # tabPanel("Load and save proxy data",
      #          fluidRow(h4("Upload proxy data"),
      #                   csvFileInput("datafile", "User data (.csv format)")),
      #          hr(),
      #          conditionalPanel("output.fileUploaded", {
      #            fluidRow(h4("Download proxy data with uncertainties"),
      #                     downloadButton("downloadData", "Download"))
      #          })
      # ),
      tabPanel(
        "Load and save proxy data", value = "load.save",
        # Input: Select a file ----
        fileInput("file1", "Choose CSV File",
                  multiple = FALSE,
                  accept = c("text/csv",
                             "text/comma-separated-values,text/plain",
                             ".csv")),

        checkboxInput("header", "Header", TRUE),


        selectInput("sep", "Separator",
                     choices = c("," = ",",
                                 ";" = ";",
                                 "tab" = "\t"),
                     selected = ";"),

        selectInput("quote", "Quote",
                     choices = c(None = "",
                                 "Double Quote" = '"',
                                 "Single Quote" = "'"),
                     selected = '"'),
      hr(),
       conditionalPanel("output.fileUploaded", {
          fluidRow(h4("Download proxy data with uncertainties"),
                   downloadButton("downloadData", "Download"))
        })

      ),
      # Time-slices -------
      tabPanel(
        "Define time-slices", value = "timeslice.pars",
        fluidRow(
          uiOutput("ui.timeslices")
          )
      )
      )),
  # Output panels -----
  mainPanel(
    tabsetPanel(id = "inTabset",
                tabPanel("Proxy error model",
                          h4("Seasonal cycle parameters"),
                          fluidRow(column(6, verbatimTextOutput("seas.amp.pars"))),
                         h4("Proxy record location"),
                         fluidRow(column(12, plotOutput("world"))),
                         fluidRow(column(12, htmlOutput("exp.corr"))),
                         # h4("Inferred parameters"),
                         # fluidRow(column(6, verbatimTextOutput("inf.par.table"))),
                         h4("Error table"),
                         fluidRow(column(12, tableOutput("err.table")))
                ),
                tabPanel("Error model diagnostics",
                         fluidRow(column(8,
                                         h4("Power spectra of error sources"),
                                         plotOutput("spec.plot"),
                                         hr(),
                                         h4("Timescale dependent error"),
                                         plotOutput("var.plot"),
                                         checkboxInput(
                                           "const.comps",
                                           "Plot constant errors?",
                                           value = FALSE),
                                         checkboxGroupInput(
                                           "var.comps",
                                           "Plot these variance components:",
                                           choices = error.components[
                                             error.components %in% c("Climate", "Total.error",
                                                                     "Orbital.mod.seas.bias",
                                                                     "Orbital.mod.seas.unc.",
                                                                     "Reference.climate") == FALSE],
                                           selected = error.components[
                                             error.components %in% c("Climate", "Total.error",
                                                                     "Seasonal.unc.",
                                                                     "Orbital.mod.seas.bias",
                                                                     "Orbital.mod.seas.unc.",
                                                                     "Reference.climate") == FALSE],
                                           inline = TRUE)
                                         )
                                  ),
                         hr(),
                         fluidRow(column(4,
                                         h4("Used parameters"),
                                         htmlOutput("used.pars"))),
                         value = "diag.panel"
                ),
                tabPanel("Uploaded proxy data",
                         h4("Uploaded proxy data"),
                         fluidRow(column(8, dataTableOutput("raw.data"))),
                         value = "proxy.data.panel"
                ),
                tabPanel("Smoothed proxy record + error",
                         #fluidRow(column(12, dataTableOutput("proxy.plot.data"))),
                         fluidRow(column(12, plotOutput("proxy.plot", height = 400))),
                         fluidRow(column(12, dataTableOutput("proxy.data.smoothed.err.table"))),
                         value = "pointwise.panel"
                ),
                tabPanel("Time-slices",
                         h4("Error on time-slices and error on difference between two time-slices"),
                         verbatimTextOutput("time.slice.err"),
                         value = "timeslice.panel"
                )
    ))
)



# Server logic -----
server <- function(input, output, session) {

  observe({
    if (input$sidebar == "load.save") {
      updateTabsetPanel(session, "inTabset",
                        selected = "proxy.data.panel")
    }

    if (input$sidebar == "timeslice.pars") {
      updateTabsetPanel(session, "inTabset",
                        selected = "timeslice.panel")
    }

    })


  # Create UI elements -----
  ## UI proxy type ----
  output$ui.proxy.chars <- renderUI({

    pars <- GetSpecPars(proxy.type = input$proxy.type)

    fluidPage(fluidRow(column(
      6,
      sliderInput(
        "n.months",
        h6("Months of proxy production per year"),
        value = pars$tau_p * 12,
        step = 0.1,
        min = 1,
        max = 12
      )
    ),
    column(
      width = 6,
      numericInput(
        "n.samples",
        label = h6(
          HTML("<I>N</I> - No. signal carriers (e.g. foraminifera) per sample")
        ),
        value = pars$N,
        step = 1,
        min = 1,
        max = Inf
      )
    )),
    fluidRow(column(
      6,
      sliderInput(
        "phi_c",
        h6("Phase of signal carrier production relative to peak of seasonal climate cycle"),
        value = 0,
        step = 0.1,
        min = -1,
        max = 1,
        post = " * pi"
      )
    ),
    column(
      width = 6,
      sliderInput(
        "delta_phi_c",
        label = h6(
          HTML("Uncertainty in phase of signal carrier production")
        ),
        value = 0.5,
        step = 0.1,
        min = 0,
        max = 1,
        post = " * 2pi"
      )
    )),
    fluidRow(
      h4("Noise parameters"),
      column(
        6,
        numericInput(
          "sigma.meas",
          h6("Measurement noise"),
          value = pars$sigma.meas,
          step = 0.01,
          min = 0,
          max = 2
        )
      ),
      column(
        6,
        numericInput(
          "sigma.ind",
          h6("Individual noise"),
          value = pars$sigma.ind,
          step = 0.1,
          min = 0,
          max = 5
        )
      ),
      column(
        6,
        numericInput(
          "sigma.cal",
          h6("Calibration bias"),
          value = pars$sigma.cal,
          step = 0.1,
          min = 0,
          max = 5
        )
      )
    ))
  })

  # UI timescale ----

  output$ui.infertauT <- renderUI({
       fluidRow(
        h4("Core sampling"),
          h6("Infer temporal parameters from uploaded proxy record."),
          checkboxInput("infertauT",
                        HTML("Infer &tau;<sub>T</sub> from proxy record?"),
                        value = if (is.null(input$file1)) {FALSE} else {TRUE})
      )
  })

  output$ui.timescale <- renderUI({
    if (input$infertauT == TRUE){
      inf.pars <- inferred.timescale.pars()
      fluidRow(
        column(
          width = 6,
          numericInput(
            "tau_T",
            h6(HTML("&tau;<sub>T</sub> - Length of the timeseries [kyr]")),
            value = inf.pars$tau_T,
            step = 0,
            min = inf.pars$tau_T,
            max = inf.pars$tau_T
          )
        ),
        column(
          width = 6,
          numericInput(
            "delta_t",
            h6(HTML("&Delta;<sub>t</sub> - Core sampling resolution [years]")),
            value = round(inf.pars$delta_t),
            step = 0,
            min = inf.pars$delta_t,
            max = inf.pars$delta_t
          )
        ),
        column(
          width = 6,
          sliderInput(
            "tau_r",
            label = h6(HTML("&tau;<sub>r</sub> - Interpreted resolution of proxy timeseries")),
            value = signif(1000 * inf.pars$tau_T / 10, 1),
            step = signif(1000 * inf.pars$tau_T / 100, 1),
            min = 0,
            max = signif(1000 * inf.pars$tau_T / 3, 1)
          )
        )
      )
      }else{
      fluidRow(
        column(
          width = 6,
          numericInput(
            "tau_T",
            h6(HTML("&tau;<sub>T</sub> - Length of the timeseries [kyr]")),
            value = 2e04/1e03,
            step = 1e03/1e02,
            min = 0e03/1e03,
            max = 1e06/1e03
          )
        ),
        column(
          width = 6,
          numericInput(
            "delta_t",
            h6(HTML("&Delta;<sub>t</sub> - Core sampling resolution [years]")),
            value = 100,
            step = 10,
            min = 0,
            max = 10000
          )
        ),
        column(
          width = 6,
          numericInput(
            "tau_r",
            label = h6(HTML("&tau;<sub>r</sub> - Interpreted resolution of proxy timeseries")),
            value = 100,
            step = 10,
            min = 0,
            max = 10000
          )
        )
      )
    }


  })

  ## UI seasonal cycle ----
  output$ui.seas.cycle <- renderUI({
    if (input$use.lat == TRUE){
      seas.pars <- seas.amp.pars()
      fluidRow(column(width = 12,
        column(
          width = 6,
          numericInput(
            "seas.amp",
            h6(paste0("Seasonal cycle amplitude ",
                      if(input$proxy.type != "d18O") {"[°C]"}else{"[‰ PDB]"})),
            value = seas.pars$mean.amp,
            step = 0.1,
            min = 0,
            max = 20
          )
        ),
        column(
          width = 6,
          numericInput(
            "sig.sq_a",
            h6("Proportional change in seasonal cycle amplitude over precessionary cycle."),
            value = seas.pars$sig.sq_a,
            step = 0.01,
            min = 0,
            max = 1
          )
        ),
        column(
          width = 6,
          numericInput(
            "phi_a",
            h6("phi_a"),
            value = 0.125,
            step = 0.005,
            min = 0,
            max = pi
          )
        )
        )
      )
    }else{
      pars <- GetSpecPars(proxy.type = input$proxy.type)
      fluidRow(column(width = 12,
        column(
          width = 6,
          numericInput(
            "seas.amp",
            h6(paste0("Seasonal cycle amplitude ",
                      if(input$proxy.type != "d18O") {"[°C]"}else{"[‰ PDB]"})),
            value = pars$seas.amp,
            step = 0.1,
            min = 0,
            max = 18
          )
        ),
        column(
          width = 6,
          numericInput(
            "sig.sq_a",
            h6("Proportional change in seasonal cycle amplitude over precessionary cycle."),
            value = pars$sig.sq_a,
            step = 0.01,
            min = 0,
            max = 0.4999
          )
        ),
        column(
          width = 6,
          numericInput(
            "phi_a",
            h6("phi_a"),
            value = 0.125,
            step = 0.005,
            min = 0,
            max = pi
          )
        )
      ))
    }
  })
  # UI power.spec -----
  output$ui.power.spec <- renderUI({
    pars <- GetSpecPars(proxy.type = input$proxy.type)
    if (input$use.lat == TRUE){
      fluidRow(column(12,
                      h4("Power Spectrum"),
                      column(
                        width = 6,
                        numericInput(
                          "clim.pow.slope",
                          h6("Slope of the power spectrum of the climate"),
                          value = pars$clim.spec.fun.args$b,
                          step = 0.1,
                          min = 0.5,
                          max = 1.5
                        )
                      )
      )
      )
    } else {
      fluidRow(column(12,
      h4("Power Spectrum"),

      column(
        width = 6,
        numericInput(
          "clim.pow.slope",
          h6("Slope of the power spectrum of the climate"),
          value = pars$clim.spec.fun.args$b,
          step = 0.1,
          min = 0.5,
          max = 1.5
        )
      ),
      column(
        width = 6,
        numericInput(
          "clim.pow.0.5",
          h6("Power spectral density at f = 1/2"),
          value = pars$clim.spec.fun.args$a,
          step = pars$clim.spec.fun.args$a/10,
          min = pars$clim.spec.fun.args$a / 10,
          max = pars$clim.spec.fun.args$a * 10
        )
      )
    )
    )}
  })

  # UI timeslices ------
  output$ui.timeslices <- renderUI({
    if (is.null(input$file1) == FALSE){
      dat <- proxy.data.smoothed()
      pars <- collect.pars()

      #d_t <- 1/(2 * max(spec.obj()[["proxy.error.spec"]]$nu))
      d_t <- pars$delta_t
      tau_T <- pars$tau_T
      min_t <- round(min(dat$age))
      max_t <- round(max(dat$age))

      fluidRow(column(12,
                      h4("Time-slices"),
                      column(
                        width = 12,
                        sliderInput(
                          "tau_1",
                          h6("Duration of time-slice 1"),
                          min = min_t,
                          max = max_t,
                          value = c(min_t, min_t + 3*d_t),
                          step = d_t,
                          dragRange = TRUE
                        )
                      ),
                      column(
                        width = 12,
                        sliderInput(
                          "tau_2",
                          h6("Duration of time-slice 2"),
                          min = min_t,
                          max = max_t,
                          value = c(min_t + 9*d_t, min_t + 12*d_t),
                          step = d_t,
                          dragRange = TRUE
                        )
                      )
      )
      )
    }else{
      pars <- collect.pars()

      d_t <- pars$delta_t
      tau_T <- pars$tau_T
      min_t <- 0
      max_t <- tau_T

      fluidRow(column(12,
                      h4("Time-slices"),
                      column(
                        width = 12,
                        sliderInput(
                          "tau_1",
                          h6("Duration of time-slice 1"),
                          min = min_t,
                          max = max_t,
                          value = c(min_t + tau_T/4, min_t + tau_T/4 + 10 * d_t),
                          step = d_t,
                          dragRange = TRUE
                        )
                      ),
                      column(
                        width = 12,
                        sliderInput(
                          "tau_2",
                          h6("Duration of time-slice 2"),
                          min = min_t,
                          max = max_t,
                          value = c(max_t - tau_T/4 - 10 * d_t, max_t - tau_T/4),
                          step = d_t,
                          dragRange = TRUE
                        )
                      )
      )
      )
    }

  })

  # load data ------
  proxy.data <- reactive({

    # input$file1 will be NULL initially.
    req(input$file1)

    tryCatch(
      {
        df <- read.csv(input$file1$datapath,
                       header = input$header,
                       sep = input$sep,
                       quote = input$quote)

        names(df)[1:2] <- c("age", "value")
      },
      error = function(e) {
        # return a safeError if a parsing error occurs
        stop(safeError(e))
      }
    )
    return(df)
  })

  # infer timescale parameters from data -----
  inferred.timescale.pars <- eventReactive(proxy.data(), {
    dat <- proxy.data()
    inf.pars <- list(tau_T = round(diff(range(dat$age))/1e03),
                     delta_t = round(mean(diff(dat$age))),
                     tau_r = round(mean(diff(dat$age)))
    )
    return(inf.pars)
  })

  output$inf.par.table <- renderPrint(inferred.timescale.pars())
  outputOptions(output, 'inf.par.table', suspendWhenHidden=FALSE)



  # create flag for whether a file has been uploaded
  output$fileUploaded <- reactive({
    return(!is.null(proxy.data()))
  })

  outputOptions(output, 'fileUploaded', suspendWhenHidden=FALSE)
  fileUploaded <- reactive({
    return(!is.null(proxy.data()))
  })


  proxy.data.smoothed <- reactive({

    df <- proxy.data()

    # return NULL if nothing has been uploaded
    if(is.null(df)) return(NULL)

    df <- SmoothIrregular(df, "age", "value", tau_smooth = input$tau_r)

    return(df)
  })


  proxy.data.smoothed.err <- reactive({
    df <- proxy.data.smoothed()
    # add uncertainties to proxy record
    spec.err <- pw.spec.err()
    df <- add_column(df, !!! spec.err)

    # return NULL if nothing has been uploaded
    if(is.null(proxy.data())) return(NULL)

    #df <- mutate_at(df, vars(-age), signif, digits = 3)
    return(df)
  })


  proxy.plot.dat <- eventReactive({proxy.data.smoothed.err() },{

   pd1 <- proxy.data() %>%
      mutate(resolution = "Raw data")

   pd2 <- proxy.data.smoothed.err() %>%
      rename(value = mean.value) %>%
      mutate(resolution = "Smoothed")

   data <- bind_rows(pd1, pd2) %>%
      select(age, Total.error, value, resolution)

    return(data)
  })




  # use location for seasonal cycle parameters -----
  seas.amp.pars <- reactive({

    if (input$use.lat){
      validate(
        need(paste(input$latitude, input$longitude) %in%
               unique(paste(breitkreuz.amp$latitude, breitkreuz.amp$longitude)),
             message = "Location on land")
      )

      orbital.pars <- RelativeAmplitudeModulation(latitude = input$latitude,
                                                  maxTimeKYear = 23,
                                                  minTimeKYear = 1,
                                                  bPlot = FALSE)


      amps <- AmpFromLocation(longitude = input$longitude, latitude = input$latitude,
                      depth.upr = input$depth[1], depth.lwr = input$depth[2],
                      proxy.type = input$proxy.type)

      seas.amp.pars <- list(
        mean.amp = signif(amps$mean.amp, 3), sig.sq_c = signif(amps$sig.sq_c, 3),
        sig.sq_a = round(orbital.pars$sig.sq_a, 3)
      )
    }else{
      seas.amp.pars <- list(seas.amp = input$seas.amp,
                            siq.sq_c = VarSine(input$seas.amp),
                            sig.sq_a = input$sig.sq_a)
    }

    return(seas.amp.pars)
  })

  output$seas.amp.pars <- renderPrint((seas.amp.pars()))


  # toggle input fields -----

  # observe({
  #   if (is.null(input$file1)) shinyjs::enable("infertauT")
  # })

  observeEvent(collect.pars(), {
    shinyjs::toggleState("tau_T", condition = input$infertauT == FALSE)
    shinyjs::toggleState("delta_t", condition = input$infertauT == FALSE)
  })

  observeEvent(collect.pars(), {
    shinyjs::toggleState("latitude", condition = input$use.lat == TRUE)
    shinyjs::toggleState("longitude", condition = input$use.lat == TRUE)
    shinyjs::toggleState("depth", condition = input$use.lat == TRUE)

    shinyjs::toggleState("seas.amp", condition = input$use.lat == FALSE)
    shinyjs::toggleState("phi_a", condition = input$use.lat == FALSE)
    shinyjs::toggleState("sig.sq_a", condition = input$use.lat == FALSE)

  })


  output$world <- renderPlot({
    if (input$use.lat){
      PlotWorld(lat = input$latitude, lon = input$longitude)
    }
  })


  # Collect parameters ------

  collect.pars <- reactive({

    pars <- list(
      tau_s = 1e03 * input$slice.thickness / input$sed.acc.rate,
      tau_b = 1e03 * input$bio.depth / input$sed.acc.rate,
      tau_p = input$n.months/12,
      tau_r = input$tau_r,
      tau_T = input$tau_T * 1e03,
      delta_t = input$delta_t,
      N = input$n.samples,
      n.k = 5,
      clim.spec.fun = "ModelSpectrum",
      clim.spec.fun.args = list(latitude = input$latitude,
                                beta = input$clim.pow.slope,
                                variable = input$proxy.type),
      # clim.spec.fun = "ClimPowerFunction",
      # clim.spec.fun.args = list(a = (1/2)^input$clim.pow.slope * input$clim.pow.0.5,
      #                           b = input$clim.pow.slope),
      sig.sq_a = input$sig.sq_a,
      sig.sq_c = VarSine(input$seas.amp),
      nu_a = 1/23e03,
      nu_c = 1/1,
      phi_a = input$phi_a,
      phi_c = input$phi_c * pi,
      delta_phi_c = input$delta_phi_c * 2* pi,
      sigma.meas = input$sigma.meas, sigma.ind = input$sigma.ind,
      sigma.cal = input$sigma.cal
    )
    return(pars)
  })

  output$used.pars <- renderTable({
    pars <- collect.pars()
    pars <- lapply(pars, function(p) {

      if (is.numeric(p)) {p <- signif(p, 3)}

      p

      })
    par.vec <- unlist(pars)
    data.frame(Argument = names(par.vec), value = as.vector(par.vec))
  })

  # calculate error spectrum -----
  spec.obj <- reactive( {
    pars <- collect.pars()
    spec.obj <- do.call(ProxyErrorSpectrum, pars)
    #spec.obj <- RenameSpecComponents(spec.obj)
  })

  output$spec.plot <- renderPlot({
    PlotSpecError(spec.obj())},
    width = 800, height = "auto",
    res = 150)

  # integrate error spec ------
  var.obj <- eventReactive(spec.obj(), {
    pars <- collect.pars()
    IntegrateErrorSpectra(spec.obj())
  })

  time.dep.var <- eventReactive(spec.obj(), {
    # pars <- collect.pars()
    #
    # s.res <- exp(seq(log(pars$delta_t), log(pars$tau_T/2), length.out = 100))
    # lst <- lapply(s.res, function(s) {
    #   as.list(IntegrateErrorSpectra(spec.obj(), tau_smooth = s,
    #                                 method = "sincfilter"))}
    #   )
    # do.call(rbind.data.frame, lst)
    IntegrateErrorSpectra(spec.obj())
  })

  # Get error from var -----
  spec.err <- eventReactive(var.obj(), {
    pars <- collect.pars()
    GetProxyError(var.obj(), timescale = input$tau_r, exclude = "Seasonal.unc.")
    })



  # Pointwise spectra, variances and errors ------
  pw.spec.err <- reactive({

    proxy.data <- proxy.data.smoothed()
    d_t <- proxy.data$mean.d_t

    pars.i <- collect.pars()

    spec.err.lst <- lapply(1:length(d_t), function(i) {

      if(is.finite(d_t[i])){
        pars.i$delta_t <- d_t[i]
        spec.obj <- do.call(ProxyErrorSpectrum, pars.i)
        #spec.obj <- RenameSpecComponents(spec.obj)

        var.obj <- IntegrateErrorSpectra(spec.obj, method = "sincfilter", tau_smooth = pars.i$tau_r)

        err.obj <- GetProxyError(var.obj, timescale = pars.i$tau_r, exclude = "Seasonal.unc.")

        err.obj <- err.obj %>%
          filter(component == "Total.error", smoothed.resolution == pars.i$tau_r) %>%
          select(component, exl.f.zero) %>%
          spread(component, exl.f.zero)

        return(err.obj)

      } else NA


    })

    err.df <- do.call(rbind.data.frame, spec.err.lst)

    return(err.df)
    })

  # time-slice error ------
  time.slice.err <- reactive({

    pars <- collect.pars()

    spec.obj <- spec.obj()$proxy.error.spec

    delta_t <- 1/(2 * max(spec.obj$nu))
    #delta_t <- pars$delta_t

    d_ts <- delta_t * round((mean(range(input$tau_2)) - mean(range(input$tau_1))) / delta_t)


    print(d_ts)

    err.ts <- ErrDiff2TimeSlices(
      spec.obj(),
      tau_1 = diff(range(input$tau_1)),
      tau_2 = diff(range(input$tau_2)),
      delta_ts = d_ts)

    err.ts

  })


  # outputs after calculations -----

  output$raw.data <- renderDataTable({
    proxy.data()
  }, options = list(pageLength = 30))


  output$proxy.data.smoothed.err.table <- renderDataTable({
    mutate_all(proxy.data.smoothed.err(), signif, digits = 3)
  }, options = list(pageLength = 8))


  output$proxy.plot <- renderPlot({

    data <- proxy.plot.dat()

    data <- data[order(data$value, na.last = FALSE),]
    ggproxyPlot(data)
    #plotly::ggplotly(ggproxyPlot(data))
  })

  # output$proxy.plot.data <- renderDataTable({
  #
  #   data <- proxy.plot.dat()
  #
  #   data <- data[order(data$value, na.last = FALSE),]
  #   data
  # })



  # Downloadable csv of selected dataset plus uncertainties ----
  output$downloadData <- downloadHandler(
    filename = function() {
      paste("output-data", ".csv", sep = "")
    },
    content = function(file) {
      write.csv(proxy.data.smoothed.err(), file, row.names = FALSE)
      #write.csv(proxy.plot.dat(), file, row.names = FALSE)
    }
  )


  output$err.table <- renderTable({

    spec.err()
    #mutate_if(spec.err(), is.numeric, signif, digits = 3)

  })

  output$var.plot <- renderPlot({

    # work around for graphics device issue with ° and permille
    unts <- ifelse(input$proxy.type == "d18O", "d18O", "degC")

    ylab <- switch(unts,
                   degC = expression("Error variance [degC"^2 ~"]"),
                   d18O = expression("Error variance [permille d"^18*"O"^2 ~"]"))


    PlotTSDVariance(time.dep.var(),
                    components = input$var.comps,
                    include.constant.errors = input$const.comps,
                    units = unts) +
      labs(y = ylab)
    },
    width = 800, height = "auto",
    res = 150)

  output$pointwise.error <- renderDataTable(
    mutate_all(pw.spec.err(), signif, digits = 3)
    )

  # expected correlation between replicate proxy records -----
  output$exp.corr <- renderText({
    ec <- ExpectedCorrelation(spec.obj())
    ec.out <- ec[ec$correlation.type == "proxy-proxy" & ec$smoothed.resolution == input$tau_r, 3]
    paste0("Expected correlation between replicate proxy records: r = ",
           round(ec.out, 2), "; R<sup>2</sup> = ",
           round((ec.out)^2, 2))
  })

  output$time.slice.err <- renderPrint(
    time.slice.err()
  )

}

# Run the app ----
shinyApp(ui = ui, server = server)