remove .R~

.gitignore

 *.o
 *.so
 *.R#
+*.R~
\ No newline at end of file

Rscript.R~

-
-########################################################################
-########################################################################
-#INITIALIZATION AND DATA LOAD
-########################################################################
-########################################################################
-
-library(TMB)
-library(dplyr)
-library(ggplot2)
-library(pracma)
-
-#clear workspace / set path
-rm(list=ls())
-this_path = rstudioapi::getActiveDocumentContext()$path
-setwd(dirname(this_path))
-
-# Compile and load neg. log-likelihood function
-#Kalman method
-compile("TMB_kalman.cpp")
-dyn.load(dynlib("TMB_kalman"))
-#Mixed-effects method
-compile("TMB_standard.cpp")
-dyn.load(dynlib("TMB_standard"))
-
-## Data Load
-df = readRDS("data.rds")
-cp = readRDS("cp.rds")
-
-# Setup data
-Y  = df[,8:12]
-Q  = df[,3:7]
-#
-houses = 1:5
-n      = length(houses)
-#
-NAobs       = apply(Y,2,FUN=function(x) as.numeric(!is.na(x)))
-NAints      = apply(NAobs,1,sum,na.rm=T)
-Y2          = Y
-Y2[is.na(Y2)] = 0
-obsVar      = 1^2
-Vmat        = 0*Q+obsVar
-
-# Upper and lower parameter bounds
-parUB      = c(500,1500,2,2)
-parLB      = c(1,1,1e-5,1e-5)
-
-# TMB data
-tmbdata = list(
-  t      = df$t,
-  Y      = as.matrix(Y2),
-  Q      = as.matrix(Q),
-  Tg     = df$Tg,
-  parUB  = parUB,
-  parLB  = parLB,
-  X0     = rep(50,n+1),
-  P0     = 15*diag(n+1),
-  Vmat   = as.matrix(Vmat),
-  NAobs  = NAobs,
-  NAints = NAints
-)
-
-# TMB parameters
-tmbpars = list(
-  Cpars = rep(logit((20-parLB[1])/parUB[1]),n),
-  Rpars = rep(logit((200-parLB[2])/parUB[2]),n),
-  Spars = c( rep(logit((1e-2-parLB[3])/parUB[3]),n), 
-             logit((1e-2-parLB[4])/parUB[4]) )
-)
-
-# Create objective function from c++ template
-nll    = MakeADFun(tmbdata,
-                   tmbpars,
-                   DLL="TMB_kalman",
-                   silent=F)
-
-# Optimize to estimate parameters
-opt  = nlminb(nll$par,nll$fn,nll$gr,control=list(iter.max=500,eval.max=500))
-pars = opt$par
-
-# Obtain smoothed state estimates with using TMB
-iobs    = which(!is.na(Y),TRUE)
-tmbdata2 = list(
-  t          = df$t,
-  iobs_row_Y = iobs[,1]-1,
-  iobs_col_Y = iobs[,2]-1,
-  iobs_row_X = iobs[,1]-1,
-  iobs_col_X = iobs[,2]-1,
-  Y          = as.matrix(Y),
-  Q          = as.matrix(Q),
-  Tg         = df$Tg,
-  sY         = as.matrix(sqrt(Vmat)),
-  parUB      = parUB,
-  parLB      = parLB,
-  simulate   = 0,
-  X0         = rep(60,n+1),
-  V0         = diag(rep(100),n+1)
-)
-tmbpars2 = list(
-  Cpars=pars[1:n],
-  Rpars=pars[(n+1):(2*n)],
-  Spars=pars[(2*n+1):(3*n+1)],
-  X    = matrix(50,nrow=nrow(Q),ncol=ncol(Q)+1)
-)
-nll2    = MakeADFun(tmbdata2,
-                    tmbpars2,
-                    random=c("X"),
-                    DLL="TMB_standard",
-                    silent=F)
-# Optimize random effects
-nll2$fn()
-sdr2   = sdreport(nll2, ignore.parm.uncertainty=T)
-pred   = summary(sdr2,"random")
-states   = c()
-statesSd = c()
-m        = nrow(tmbpars2$X)
-for(i in 1:(n+1)){
-  ran       = ((i-1) * m + 1):(i*m)
-  states    = cbind(states,pred[ran,1])
-  statesSd  = cbind(statesSd,pred[ran,2])
-}
-t = df$t
-
-p = ggplot() +
-  geom_polygon(aes(x=c(t,rev(t)),y=c(states[,6]+2*statesSd[,6],rev(states[,6]-2*statesSd[,6]))), fill="blue", alpha = 0.4) +
-  geom_line(aes(x=t,y=states[,6],color="Street"),lwd=1)
-for(i in 1:n){
-  p = p +
-    geom_point(aes_(x=t,y=Y[,i],color="House",shape="Observations"),size=0.25) +
-    geom_line(aes_(x=t,y=states[,i],color="House"),size=0.25)
-}
-p = p +
-  labs(
-    title = "Smoothed State Estimates",
-    x = "Time",
-    y = "Degrees Celcius",
-    color = "",
-    shape = ""
-  ) +
-  theme(
-    plot.title = element_text(hjust=0.5,size=20), #plot title
-    plot.background = element_rect(color="white",fill="white"), #around inner window
-    plot.margin = margin(r=10),
-    panel.background = element_rect(fill = "white", colour = "white"), #inner window
-    axis.line.y = element_line(color="black"),
-    axis.line.x = element_line(color="black"),
-    legend.key = element_rect(fill="white"),
-    legend.title.align = 0.5,
-    legend.text = element_text(size=10),
-    legend.box = "horizontal",
-    legend.position = "top",
-    legend.key.width = unit(10,"pt"),
-    legend.key.height = unit(10,"pt"),
-    legend.key.size = unit(3,"line"),
-    legend.margin = margin(t=0,r=0,b=0,l=0,unit="pt"),
-    legend.background = element_rect(fill=alpha('white', 0)),
-    axis.title.x = element_text(size=10, margin = margin(t=0, r=0, b=0, l=0),color="black"),
-    axis.title.y = element_text(size=10, margin = margin(t=0, r=5, b=0, l=5),color="black"),
-    axis.text = element_text(size=10,color="black"), #the tick labels
-    axis.ticks = element_line(size=0.5,color="black") #the actual tick lines
-  )
-print(p)
-