---
title: "Applying HMM to tracks for SSFs"
author: "Rachel Mawer"
date: "2023-06-23"
output: 
  html_document:
    toc: true
    toc_float: true
    theme: cerulean
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F,message=F)

#main packages #delete/add as needed
library(dplyr)
library(ggplot2)
library(knitr)
library(ggpubr)
library(plotly)
library(pracma) #for standard error
library(reactable)
library(momentuHMM)
library(amt)
library(sf)
```

## About

This file is to apply HMMs to the tracks used in previous SSF analysis. HMMs will be made for all fish of a species, then used to assign behavioural states to tracks.

In this file:

* calculating straightness indexes for steps 
* fitting HMMs to the data, and find the best model
* assigning states to fish tracks
* saving resulting files for SSF modelling

```{r read in data, calculate straightness indexes + nsd}
# will read in all data, and form 2 dataframes - barbel and grayling
files <- list.files("data/SSF dataframes",full.names = T)


#loop, read in, calculate SI, save resulting file and then add to a barbel or grayling df
#for the latter, can do a "if df exists, else" thing


#for testing
l <-1
d <-1

#windows for SI, 1min, 5min and 10min
windows <- list(60/20,300/20,600/20)


for(l in 1:length(files)){
  file_to_read <- files[l]
  dat_a <- read.csv(file_to_read) %>% filter(case_==TRUE)
  
  dat <- dat_a
  
  dat$time <- as.POSIXct(dat$t1_)


  for(i in 1:length(windows)){
     window <- windows[[i]]
  
  # Get n to find SI from middle-point of window. 
  n <- floor(window/2)
  
  if(i==1){
    dat$SI <- NA               # Initiate SI-value in df
    dat$nsd <- NA
    dat$displ <- NA
     X <- split(dat, list(dat$fish_id,dat$track,dat$burst_),drop=T)
  }else{
    dat_with_SI$SI <- NA              # Initiate SI-value in df IF adding multiple    to the tracks
    dat_with_SI$nsd <- NA
    dat_with_SI$displ <- NA
    X <- split(dat_with_SI, list(dat_with_SI$fish_id,dat_with_SI$track,dat_with_SI$burst_),drop=T)
  }
  

  for (d in 1:length(X)) { #
    track <- data.frame(X[d])                  # Call individual track data
    if(i==1){
      colnames(track) <- colnames(dat)    # Rename columns for continuity
    }else{
      colnames(track) <- colnames(dat_with_SI)
    }

    #doing it this way instead of via SL so things are in the right order
    track$step <- sqrt((track$x1_-lag(track$x1_))^2+(track$y1_-lag(track$y1_))^2)
  
    #dont need below if repeating
    if(i==1){
    track <- track %>% st_as_sf(coords = c('x1_', 'y1_'), crs = 'EPSG:32632', remove=F) # Convert to sf-object # using start coords
    }
    
    print(paste('Track', d, '/', length(X)))
    pb <- txtProgressBar(0, nrow(track), style=3)
    # Loop for even number of detections in window
    if ((window %% 2) == 0) {
      for (c in n+1:(nrow(track)-n)) {
        if(nrow(track)>=(2*n+1)){ #if track is smaller than n+8 e.g. max mid point, do nothing, otherwise yeah
          setTxtProgressBar(pb, c)
          TDist <- sum(track$step[(c-(n)):(c+n-1)],na.rm=T)                       
          EDist <- as.numeric(st_distance(track$geometry[c-(n)], track$geometry[c+n-1])) # Calculate effective distance
          track$SI[c] <- EDist / TDist # Calculate SI
          track$nsd[c] <- EDist^2
          track$displ[c] <-EDist
        }else{
          track$SI <- NA
        }
      }
    }
    # Loop for uneven number of detections in window
    if ((window %% 2) == 1) {
      for (c in n+1:(nrow(track)-n)){
        if(nrow(track)>=(2*n+1)){ #if track is smaller than n+8 e.g. max mid point, do nothing, otherwise yeah
          setTxtProgressBar(pb, c)
          TDist <- sum(track$step[(c-n):(c+(n))],na.rm=T)
          EDist <- st_distance(track$geometry[c-n], track$geometry[(c+(n))])
          track$SI[c] <- EDist / TDist
          track$nsd[c] <- EDist^2
          track$displ[c] <-EDist
        } else{
          track$SI <- NA
          track$nsd <- NA
          track$displ <- NA
        }
      }
    }
    close(pb)
    
    # Loop to combine individual tracks into 1 df
    if (d == 1) {
      dat_with_SI <- track
    }
    else {
      dat_with_SI <- rbind(dat_with_SI, track)
    }
  }
  
  if(i==1){
    dat_with_SI <- rename(dat_with_SI,SI_1min=SI,displ_1min=displ,nsd_1min = nsd)
  }
  if(i==2){
    dat_with_SI <- rename(dat_with_SI,SI_5min=SI,displ_5min=displ,nsd_5min = nsd)
  }
  if(i==3){
    dat_with_SI <- rename(dat_with_SI,SI_10min=SI,displ_10min=displ,nsd_10min = nsd)

  }
  
}

  #save csv to df for later individual state assignment
  fish_id <- unique(dat$fish_id)
  #make sep dat_with_SI without geometry
  dat_to_save <- dat_with_SI %>% select(-geometry)
  
  write.csv(dat_to_save,paste0("data/SSF dataframes with SI etc/",fish_id,".csv"),row.names = F)
  
  species <- unique(dat$species)

  if(species == "Barbel"){
    if(exists("barb_all")==FALSE){
    barb_all <- dat_to_save
    } else{
      barb_all <- rbind(barb_all,dat_to_save)
    }
  }
  
  if(species == "Grayling"){
  if(exists("gray_all")==FALSE){
  gray_all <- dat_to_save
  } else{
    gray_all <- rbind(gray_all,dat_to_save)
  }
}

  
}

write.csv(barb_all,"data/large df with SI/barbel.csv",row.names=F)
write.csv(gray_all,"data/large df with SI/grayling.csv",row.names=F)

```

## Fit HMM to each species


### Grayling

```{r fit hmm to grayling with SI 10,warning=F}

### prepare data 
gray_dat_to_prep <- read.csv("data/large df with SI/grayling.csv") %>% select(-step)


gray_dat_to_prep <- gray_dat_to_prep %>% mutate(ID = paste(fish_id, approach, sep = '-')) %>% dplyr::arrange(ID, t1_) #needs the arrange in order to prep

gr_data <- prepData(gray_dat_to_prep, coordNames = c('x1_', 'y1_'), type = 'UTM') %>% # Prepping data for HMM
  mutate(fish_id = as.factor(fish_id)) 


dist <- list(SI_1min  = 'beta',
             SI_5min  = 'beta',
             SI_10min = 'beta',
             disp_10min = 'gamma',
             disp_5min = 'gamma')

niter <- 50
gr_list_1min  <- vector('list', length = niter)
gr_list_5min  <- vector('list', length = niter)
gr_list_10min <- vector('list', length = niter)


for (n in 1:niter) {
  SI1_mean <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999)) # randomly picking mean for si for the 2 states
  SI1_sd <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999)) #randomly picking sd for si for the 2 states
  
  Par0 <- list(SI_1min = c(SI1_mean, SI1_sd)) #list parameters
  
  gr_list_1min[[n]] <- fitHMM(gr_data, nbStates = 2,
                              dist[c(1)], Par0)
}

for (n in 1:niter) {
  SI5_mean <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999))
  SI5_sd <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999))
  
  Par0 <- list(SI_5min = c(SI5_mean, SI5_sd))
  
  gr_list_5min[[n]] <- fitHMM(gr_data, nbStates = 2,
                              dist[c(2)], Par0)
}

for (n in 1:niter) {
  SI10_mean <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999))
  SI10_sd <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999))
  
  Par0 <- list(SI_10min = c(SI10_mean, SI10_sd))
  
  gr_list_10min[[n]] <- fitHMM(gr_data, nbStates = 2,
                               dist[c(3)], Par0)
}


gr_SI1LL  <- unlist(lapply(gr_list_1min, function(m) m$mod$minimum))
gr_SI5LL  <- unlist(lapply(gr_list_5min, function(m) m$mod$minimum))
gr_SI10LL <- unlist(lapply(gr_list_10min, function(m) m$mod$minimum))

gr_SI1 <- gr_list_1min[[which.min(gr_SI1LL)]]
gr_SI5 <- gr_list_5min[[which.min(gr_SI5LL)]]
gr_SI10 <- gr_list_10min[[which.min(gr_SI10LL)]]

saveRDS(gr_SI10,"HMMs/gr_SI10.rds")
saveRDS(gr_SI1,"HMMs/gr_SI1.rds")
saveRDS(gr_SI5,"HMMs/gr_SI5.rds")

```

```{r compare AIC of grayling models}
gr_SI10 <-readRDS("HMMs/gr_SI10.rds")
gr_SI1 <- readRDS("HMMs/gr_SI1.rds")
gr_SI5 <- readRDS("HMMs/gr_SI5.rds")
AIC(gr_SI1,gr_SI5,gr_SI10)


```

Taking SI_10 model

```{r assign states}
#to edit this code, from Jelger
gr_SI10 <-readRDS("HMMs/gr_SI10.rds")

gr_data$State <- as.factor(viterbi(gr_SI10))


write.csv(gr_data,"data/steps with states/grayling_all.csv")

```

```{r split gr data into individual fish}

fish_ids <- unique(gr_data$fish_id)

for(i in 1:length(fish_ids)){
  id <- fish_ids[i]
  dat <- gr_data %>% filter(fish_id==id)
  write.csv(dat,paste0("data/steps with states/individual fish/",id,".csv"))


}
```

### Barbel

```{r fit hmm to barbel with SI 10,warning=F}

### prepare data 
barb_dat_to_prep <- read.csv("data/large df with SI/barbel.csv") %>% select(-step)


barb_dat_to_prep <- barb_dat_to_prep %>% mutate(ID = paste(fish_id, approach, sep = '-')) %>% dplyr::arrange(ID, t1_) #needs the arrange in order to prep

barb_data <- prepData(barb_dat_to_prep, coordNames = c('x1_', 'y1_'), type = 'UTM') %>% # Prepping data for HMM
  mutate(fish_id = as.factor(fish_id)) 


dist <- list(SI_1min  = 'beta',
             SI_5min  = 'beta',
             SI_10min = 'beta',
             disp_10min = 'gamma',
             disp_5min = 'gamma')

niter <- 50
barb_list_1min  <- vector('list', length = niter)
barb_list_5min  <- vector('list', length = niter)
barb_list_10min <- vector('list', length = niter)


for (n in 1:niter) {
  SI1_mean <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999)) # randomly picking mean for si for the 2 states
  SI1_sd <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999)) #randomly picking sd for si for the 2 states
  
  Par0 <- list(SI_1min = c(SI1_mean, SI1_sd)) #list parameters
  
  barb_list_1min[[n]] <- fitHMM(barb_data, nbStates = 2,
                              dist[c(1)], Par0)
}

for (n in 1:niter) {
  SI5_mean <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999))
  SI5_sd <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999))
  
  Par0 <- list(SI_5min = c(SI5_mean, SI5_sd))
  
  barb_list_5min[[n]] <- fitHMM(barb_data, nbStates = 2,
                              dist[c(2)], Par0)
}

for (n in 1:niter) {
  SI10_mean <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999))
  SI10_sd <- runif(n=2, min=c(0.001, 0.5), max=c(0.5, 0.999))
  
  Par0 <- list(SI_10min = c(SI10_mean, SI10_sd))
  
  barb_list_10min[[n]] <- fitHMM(barb_data, nbStates = 2,
                               dist[c(3)], Par0)
}


barb_SI1LL  <- unlist(lapply(barb_list_1min, function(m) m$mod$minimum))
barb_SI5LL  <- unlist(lapply(barb_list_5min, function(m) m$mod$minimum))
barb_SI10LL <- unlist(lapply(barb_list_10min, function(m) m$mod$minimum))

barb_SI1 <- barb_list_1min[[which.min(barb_SI1LL)]]
barb_SI5 <- barb_list_5min[[which.min(barb_SI5LL)]]
barb_SI10 <- barb_list_10min[[which.min(barb_SI10LL)]]

saveRDS(barb_SI10,"HMMs/barb_SI10.rds")
saveRDS(barb_SI1,"HMMs/barb_SI1.rds")
saveRDS(barb_SI5,"HMMs/barb_SI5.rds")

```

```{r compare AIC of barbling models}
barb_SI10 <-readRDS("HMMs/barb_SI10.rds")
barb_SI1 <- readRDS("HMMs/barb_SI1.rds")
barb_SI5 <- readRDS("HMMs/barb_SI5.rds")
AIC(barb_SI1,barb_SI5,barb_SI10)


```

Taking the SI_10 model.

```{r assign states}
#to edit this code, from Jelger
barb_SI10 <-readRDS("HMMs/barb_SI10.rds")

barb_data$State <- as.factor(viterbi(barb_SI10))


write.csv(barb_data,"data/steps with states/barbel_all.csv")

```

```{r split barb data into individual fish}
barb_data<- read.csv("data/steps with states/barbel_all.csv")

fish_ids <- unique(barb_data$fish_id)



for(i in 1:length(fish_ids)){
  id <- fish_ids[i]
  
  
  dat <- barb_data %>% filter(fish_id==id)
  write.csv(dat,paste0("data/steps with states/individual fish/",id,".csv"))


}
```