library('data.table')
library('reshape2')
library('plyr')
library('HiddenMarkov')
library('ggplot2')
library('RSQLite')
library('VGAM') 
library('forecast')
library('sm')

source('lib.R')

compare_state <- 2
nsims <- 1200
nblocks <- 1
dbname <- 'hmm_sim.db'
simtable <- 'simulations'
intable <- 'input_data'

m <- dbDriver('SQLite')
con <- dbConnect(m, dbname = dbname)
if(dbExistsTable(con, simtable)) dbRemoveTable(con, simtable)

dt <- as.data.table(dbReadTable(con, intable))

np <- nrow(dt)
sim_length <- np

options(warn=-1)
pars <- get.named.parlist(dt$x, 2, 'gamma', 'same.sd')
hmm <- dthmm(dt$x, Pi_init(2), c(.4,.6), 'gamma', pars)

# have to be careful here, the Algorithm will randomly switch which number
# its assigns to the high or low state,
sink('/dev/null')
hmm <- BaumWelch(hmm)
sink()

#fitglm = glm(state-1 ~ amo + pdo, data=dt, family=binomial())
vglmfit = vglm(state ~ amo + pdo + enso + state_lag1, multinomial, data=dt)
#fitglm = glm(state-1 ~ amo+pdo+enso+pdsi, data=dt, family=binomial())

#bestmodel = fitglm
#sink('/dev/null')
#bestmodel = stepAIC(fitglm)
#sink()

#statepred <- predict(fitglm,dt,type="response")
#statepred <- statepred/max(statepred)
statep <- predict(vglmfit, type="response")
#if(compare_state == 2){
#    maxp <- max(statep[,1])
#    statep[,1] <- statep[,1]/maxp
#    statep[,2] <- 1-statep[,1]
#}
#####################################
# Verify the predictions are working
#beta = coefficients(vglmfit,1)
#X=t(data.frame(ones=1,amo=dt$amo,pdo=dt$pdo))
#exp(beta[,1] %*% X)/(1+exp(beta[,1] %*% X)+ exp(beta[,2] %*% X))
####################################

rank1 <- function(x,r)rank(c(r,cumsum(x)))[1]
arfit <- ar(dt$x)

printf('Simulating...')
pb <- txtProgressBar(1,nsims,style=3)
for(j in 1:nsims){
    setTxtProgressBar(pb,j)

    sim_states <- apply(statep,1,rank1,runif(1))
    sims_climate <- sim_hmm(hmm, sim_length, sim_states)$x

    sims_imposed <- sim_hmm(hmm, sim_length, dt$state)$x
    sims_hmm_obj <- sim_hmm(hmm, sim_length)
    sims_hmm_state <- sims_hmm_obj$y
    sims_hmm <- sims_hmm_obj$x

    sims_ar <- simulate(arfit,nrow(dt))


    dbWriteTable(con,'simulations', row.names = FALSE, append=TRUE,
        rbind(
            data.table(sim=j, i=1:sim_length, value=sims_climate, state=sim_states, kind='HMC'),
            data.table(sim=j, i=1:sim_length, value=sims_imposed, state=dt$state, kind='HMI'),
            data.table(sim=j, i=1:sim_length, value=sims_hmm, state=sims_hmm_state, kind='HMG'),
            data.table(sim=j, i=1:sim_length, value=sims_ar, state=NA, kind='AR')
        )
    )

}

close(pb)

simulations <- as.data.table(dbReadTable(con,'simulations'))

printf('Simulation PDFs...')
dens <- ddply(simulations,.(sim,kind),summarise,
    density=sm.density(value,eval.points=seq(0,30,.5),display='none')$estimate,
    .progress=progress_text())
dens$x<- seq(0,30,.5)

dbWriteTable(con,'simulation_pdfs', dens, row.names = FALSE, append=FALSE, overwrite=TRUE)