########################################
#Weithman et al. 
# Seasonal movements of post-breeding royal terns (Thalasseus maximus) in Virginia 
#########################################
# Copyright (c) [2024] [Daniel H. Catlin]
# 
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#   
#   The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
# 
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
install.packages("momentuHMM")
library(momentuHMM)

delta <- 500
knownStates<-rep(NA,nrow(HMM_multi$miData[[1]]))
center1Ind<-as.numeric(HMM_multi$miData[[1]]$center1.dist) < delta  &  as.numeric(HMM_multi$miData[[1]]$step) < 500 
center2Ind<-as.numeric(HMM_multi$miData[[1]]$center1.dist < delta) #these are used to code the known loafing states. Could consider using other definitions (e.g., airspeed)

knownStates[which(center1Ind==1)]<- 1


for(j in 1:100){
HMM_multi$miData[[j]]$Time_Away[1] <- 0 
HMM_multi$miData[[j]]$Days_Away[1] <- 0 
for (i in 2:length(HMM_multi$miData[[j]]$Time)) {
  if(center2Ind[i] == 0 & HMM_multi$miData[[j]]$ID[i] == HMM_multi$miData[[j]]$ID[i-1] ){
    HMM_multi$miData[[j]]$Time_Away[i] <-  HMM_multi$miData[[j]]$Time_Away[i-1] + 1
  } else {
    HMM_multi$miData[[j]]$Time_Away[i] <- 0
  }
  HMM_multi$miData[[j]]$Days_Away[i] <-  HMM_multi$miData[[j]]$Time_Away[i]/24
}
HMM_multi$miData[[j]]$TempZ <- (HMM_multi$miData[[j]]$Temp - mean(HMM_multi$miData[[j]]$Temp) )/ sd(HMM_multi$miData[[j]]$Temp)
HMM_multi$miData[[j]]$airspeedkm <- HMM_multi$miData[[j]]$airspeed *1.85 #converts from kts to km/hr
}
###################################################################
stateNames <- c("Resting", "Exploratory")
###################################################################
##################################################################

nbStates<-2


dist = list(step = "gamma", angle = "wrpcauchy")

# initial parameters for a null model
Par0_null <- list(step=c(100,10000,
                       100,10000),
                
                angle=c(00.01,0.01))

          ind_matrix <- matrix(c( NA,NA
                                  )  , nrow= 1 )
# fit basic model
null <- fitHMM(HMM_multi$miData[[1]], nbStates = nbStates, 
              dist = dist, Par0 = Par0_null, 
              fixPar = list(beta =  ind_matrix),
              estAngleMean = list(angle=F), 
              circularAngleMean=list(angle=T),
              knownStates = knownStates,
              stateNames = stateNames)


mi.null <- MIfitHMM(HMM_multi$miData, nbStates = nbStates, 
               dist = dist, Par0 = Par0_null, 
               fixPar = list(beta =  ind_matrix),
               estAngleMean = list(angle=F), 
               circularAngleMean=list(angle=T),
               knownStates = knownStates,
               stateNames = stateNames)
mi.null.pool <- MIpool(mi.null$HMMfits)