library('HiddenMarkov')
library('data.table')

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

x <- read.table("data/lf-annual-cy-af.txt") 
x <- ts(x[,2],start=x[1,1])/10^6
xdt <- data.table(date=as.vector(time(x)),x=as.vector(x))

climate_m <- as.data.table(read.csv('data/winter_climate_indicies.csv'))
climate <- as.data.table(dcast(climate_m,date~variable))

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

options(warn=-1)
pars <- get.named.parlist(x, compare_state, 'gamma', 'same.sd')
hmm <- dthmm(x, Pi_init(compare_state), delta_init(compare_state), 'gamma', pars)


sink('/dev/null')
hmm <- BaumWelch(hmm)
sink()
#options(warn=2)

# The HMM model will randomly switch which number it assigns to a state (1 or 2)
# but the decoding will always be the same shape. Force it to do:
# 1 - low state 
# 2 - high state
decoding <- data.table(state=Viterbi(hmm),date=seq(1906,by=1,len=length(x)))
if(decoding$state[1] == 1){
    decoding[state==1]$state <- 99L
    decoding[state==2]$state <- 1L
    decoding[state==99]$state <- 2L
}
decoding <- lagged_dt(decoding,1)

#decoding_lag <- lagged_dt(decoding,3)
#ydt_lag <- lagged_dt(ydt,3)

#dt <- merge(climate,ydt_lag,by='date')
#dt <- merge(dt,decoding_lag,by='date')
dt <- merge(climate,decoding,by='date')
dt <- merge(dt,xdt,by='date')
dbWriteTable(con, intable, dt, row.names = FALSE, overwrite=TRUE)

dbDisconnect(con)