#
# This is a Shiny web application. You can run the application by clicking
# the 'Run App' button above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#

library(shiny)

library(rjags)
library(tidyverse)
library(rjags)
library(ggmcmc)
library(polspline)
library(propagate)



densitycm2<-seq(0,0.4,length=100)
rate<-10*densitycm2/(1+10*densitycm2*2) +rnorm(100,0,0.01)
d<-data.frame(densitycm2,rate)

#d<-read.csv("/home/aqm/course/data/buntings.csv")
#write.csv(d,"buntings.csv",row.names=FALSE)
d<-read.csv("buntings.csv")
d$densitycm2<-d$density/10000

hol_mod<-"
model{

#Likelihood
for( i in 1:n)
{
  
  y[i] ~ dgamma( ((mu[i]^2)/sig), (mu[i]/sig))
  #y[i] ~ dnorm(mu[i],inv.var)
  mu[i]<- (1-v)*a*x[i]/(1+a*x[i]*H)
  
}

#priors
a~dgamma( ((a_mu^2)/a_sig), (a_mu/a_sig))
H~dgamma( ((H_mu^2)/H_sig), (H_mu/H_sig))
v ~ dunif(vmin,1)

sig~dunif(0,10)
inv.var   ~ dgamma(0.01, 0.01)
sigma     <- 1/sqrt(inv.var)

}"

# Define UI for application that draws a histogram
ui <- fluidPage(
  tags$head(
    tags$style(
      HTML(".shiny-notification {
           height: 100px;
           width: 800px;
           position:fixed;
           top: calc(50% - 50px);;
           left: calc(50% - 400px);;
           }
           "
      )
      )
      ),
   
   # Application title
   titlePanel("Jags model for Bunting data"),
   
   # Sidebar with a slider input for number of bins 
   sidebarLayout(
      sidebarPanel(
         h3("Set the priors:"),
         sliderInput("v","Vigilance lower %:", min = 0, max = 60,value = 0,step=1),
         sliderInput("a_mu","a_mu:", min = 1, max = 50,value = 12),
         sliderInput("a_sig","a_sig:", min = 0.01, max = 1,value = 0.1,step=0.01),
         sliderInput("H_mu","H_mu:", min = 0.1, max = 3,value = 2,step=0.01),
         sliderInput("H_sig","H_sig:", min = 0.1, max = 10,value = 1,step=0.1)
         
      ),
      
      # a_mu=12,a_sig=0.1,H_mu=2,H_sig=1
      mainPanel(
        
        # Output: Tabset w/ plot, summary, and table ----
        tabsetPanel(type = "tabs",
                    tabPanel("Plot", plotOutput("Plot")),
                    tabPanel("Prior plot", plotOutput("PPlot")),
                    tabPanel("Summary", verbatimTextOutput("summary")),
                    tabPanel("Prior Summary", verbatimTextOutput("psummary")),
                    tabPanel("Diagnostics", plotOutput("diag")),
                    tabPanel("Prior diagnostics", plotOutput("pdiag"))
        )
        
      )
   )
)

# Define server logic required to draw a histogram
server <- function(input, output) {
   
  observe({
    xlabel <-expression("Density of seeds per cm"^{2})
    ylabel<-"Feeding rate seeds per second"
    
    a_mu<-input$a_mu
    a_sig<-input$a_sig
    H_mu<-input$H_mu
    H_sig<-input$H_sig
    v<-input$v/100
    
    dat<-list(n=length(d$densitycm2),
              x=d$densitycm2,
              y=d$rate,
              a_mu=a_mu,
              a_sig=a_sig,
              H_mu=H_mu,
              H_sig=H_sig,
              vmin=v)
    
    prior_dat<-list(n=1,
              x=d$densitycm2,
              y=d$rate,
              a_mu=a_mu,
              a_sig=a_sig,
              H_mu=H_mu,
              H_sig=H_sig,
              vmin=v)
    
    
    withProgress(message = 'Running MCMC sampling: Plase wait about 30 seconds:', value = 0.2, {
    mod1 <- jags.model(textConnection(hol_mod), data = dat, n.chains = 1)
    incProgress(0.1)
    prior_mod <- jags.model(textConnection(hol_mod), data = prior_dat, n.chains = 1)
    incProgress(0.1)
    update(mod1, 10000)
    incProgress(0.1)
    update(prior_mod, 10000)
    mod1_sim <- coda.samples(model = mod1, variable.names = c("a","H","v"), n.iter = 10000)
    incProgress(0.2)
    prior_mod_sim <- coda.samples(model = prior_mod, variable.names = c("a","H","v"), n.iter = 10000)
    incProgress(0.2)
    
    })
    ms <-ggs(mod1_sim) 
    ms %>% spread(Parameter, value) -> ms
    
    x<-seq(0,0.4,length=100)
    pred<-function(x=x,mt=ms){
      pred<-mt$a*x/(1+mt$a*x*mt$H)
      return(list(quantile(pred,c(0.025,0.5,0.975))))
    }
    
    pred_mod<-data.frame(x=x,do.call("rbind",sapply(x,pred)))
    names(pred_mod)<-c("x","lwr","mid","upr")
    
    g0<-ggplot(data=d,aes(x=densitycm2,y=rate)) + xlab(xlabel) + ylab(ylabel)
    g1<-g0+geom_point()+ geom_line(data=pred_mod,aes(x=x,y=lwr),col="blue") + geom_line(data=pred_mod,aes(x=x,y=upr),col="red") + geom_line(data=pred_mod,aes(x=x,y=mid),col="black")
    
    
    ####
    ms2 <-ggs(prior_mod_sim) 
    ms2 %>% spread(Parameter, value) -> ms2
    
    x<-seq(0,0.4,length=100)
    pred<-function(x=x,mt=ms2){
      pred<-mt$a*x/(1+mt$a*x*mt$H)
      return(list(quantile(pred,c(0.025,0.5,0.975))))
    }
    
    prior_mod<-data.frame(x=x,do.call("rbind",sapply(x,pred)))
    names(prior_mod)<-c("x","lwr","mid","upr")
    
   
    gp<-g0+geom_point()+ geom_line(data=prior_mod,aes(x=x,y=lwr),col="blue") + geom_line(data=pred_mod,aes(x=x,y=upr),col="red") + geom_line(data=pred_mod,aes(x=x,y=mid),col="black")
    ######
    
    
    output$Plot <- renderPlot(g1)
    output$PPlot <- renderPlot(gp)
    output$diag <- renderPlot(plot(mod1_sim))
    output$pdiag <- renderPlot(plot(prior_mod_sim))
    output$summary <- renderPrint({
    summary(mod1_sim)
    })
    output$psummary <- renderPrint({
      summary(prior_mod_sim)
    })
  })
   
}

# Run the application 
shinyApp(ui = ui, server = server)