library(timeSeries) # to be able to use this function : colProds (the product of column elements
library(plyr)
library(reshape)
library(lubridate)

library(zoo)


rm(list=ls(all=T))

HART <- read.csv("surv.table.CARIBOU.20170220_pj.csv")  #CARIBOU SURVIVAL DATA FILE
# assign starting dates
MONTHS <- data.frame()

HART$s_dd <- as.Date(HART$start,  origin="1899-12-30")

HART$end2 <- HART$end
HART$end2[is.na(HART$end2)] <- 41429

# assign deaths from disparate data bases.
HART$e_dd <- as.Date(HART$end2,  origin="1899-12-30")
HART$dead  <- 0
HART$dead [HART$DEATH..0.OR.1. == "Mortality"] <- 1
HART$dead [HART$DEATH..0.OR.1. == "dead"] <- 1
HART$dead [HART$DEATH..0.OR.1. == "Dead"] <- 1
HART$dead [HART$DEATH..0.OR.1. == "Dropped/dead not recovered"] <- 1 
HART$dead [HART$DEATH..0.OR.1. == "MORT"] <- 1 
HART$dead [HART$DEATH..0.OR.1. == "Lost"] <- 1 
HART$dead [HART$DEATH..0.OR.1. == 1] <- 1 

# CORRECT MINOR DATA BASE ERRORS
HART2 <-  subset(HART, !(CARIBOU_ID %in% c(66) & SUBPOPULATION == 'ColNorth')) # dbase error
HART2 <-  subset(HART2, !(CARIBOU_ID %in% c(67) & SUBPOPULATION == 'ColNorth')) # dbase error
HART2 <-  subset(HART2, !(CARIBOU_ID %in% c(67) & SUBPOPULATION == 'ColSouth')) # dbase error
HART2 <- subset(HART2 , select = -c(Comment.about.major.monitoring.gaps..if.animals.should.be.cencored.within.this.interval ) ) # remove this column from dataframe

####AGE CLASS
HART3 <-  subset(HART2, AGE != "C") # remove calves from survival analyses

### SELECT APPROPRIATE SUBPOPULATION
HART5 <- subset (HART3, SUBPOPULATION== "ColNorth" )# |"WellsGrey" SUBPOPULATION== "ColSouth" | SUBPOPULATION== "PurcellsSouth" | SUBPOPULATION== "Nakusp" | SUBPOPULATION== "Kennedy Siding" | SUBPOPULATION== "Moberly" | SUBPOPULATION== "Quintette") #Nakusp PurcellsSouth Kennedy Siding


## NOTE FOR MOOSE SURVIVAL, SUPOPULATIONS ARE "Columbia" and "Headwaters"
## NOTE FOR MOOSE SURVIVAL, CHANGE CARIBOU_ID TO MOOSE_ID THROUGHOUT

# select dates to bound the survival analysis. SEE PAPER FOR APPROPRIATE YEAR BOUNDARIES

search_start <- '2004-04-01' #CN Before TP (time period) 1
search_end <- '2006-03-31'



#add days_intersect column - NUMBER OF DAYS IN RISK PERIOD
HART5$days_intersect <- 0

#add died_in_range column NUMBER OF DEATHS IN RISK PERIOD

HART5$died_in_range <- 0


#<march_2015>

#march 2015, add 4 columns as per:
#  the columns "days_intersect" and "died_in_range" -  
#  Subdivided into months 5 to 10, and months 11 to 4.  
#  So the days_intersect and died_in_range columns should stay the same, 
#  but the subdivided data (days_intersect_5_10 and days_intersect_11_4) should sum exactly to  days_intersect (similar logic for "died_in_range" 
#  

HART5$days_intersect_5_10 <- 0
HART5$days_intersect_11_4 <- 0

HART5$died_in_range_5_10 <- 0
HART5$died_in_range_11_4 <- 0


#LOOP THROUGH EACH CARIBOU TO SEE IF IT DIED IN RISK PERIOD

for (i in 1:length(HART5$CARIBOU_ID)) 
{
	if (HART5[i,]$s_dd <= search_end && HART5[i,]$e_dd >= search_start)
	{
		#we have intersection...now calculate
		max_start <- as.Date(HART5[i,]$s_dd)
		if (search_start > max_start)
		{
			max_start <- as.Date(search_start)
		}

		min_end <- as.Date(HART5[i,]$e_dd)
		if (search_end < min_end)
		{
			min_end <- as.Date(search_end)
		}

		days <- (min_end - max_start) + 1
		HART5[i, ]$days_intersect <- days
    
		#<march_2015>
    if (days >= 0)  #the seq function will throw an error for these strange cases.
    {
  		tmp_intersect_5_10 <- 0
  		tmp_intersect_11_4 <- 0
  		my_day_seq <- as.yearmon(seq(max_start, min_end, "day"))
  
      for (j in 1:length(my_day_seq))
  		{
        #print (format(my_day_seq[j], format="%m"))
  		  if (format(my_day_seq[j], format="%m") %in% c("05", "06", "07", "08", "09", "10"))
  		  {
  		    tmp_intersect_5_10 <- tmp_intersect_5_10 + 1
  		  }
  		  else if (format(my_day_seq[j], format="%m") %in% c("11", "12", "01", "02", "03", "04"))
  		  {
  		    tmp_intersect_11_4 <- tmp_intersect_11_4 + 1
  		  }
  		}
      HART5[i,]$days_intersect_5_10 <- tmp_intersect_5_10
  		HART5[i,]$days_intersect_11_4 <- tmp_intersect_11_4
    }
		#</march_2015>


		#check if died in range
		#if there is a 1 for death and the e_dd is contained in the search span
		if (HART5[i,]$dead == 1)
		{
			if (HART5[i,]$e_dd >= search_start && HART5[i,]$e_dd <= search_end)
			{
				HART5[i, ]$died_in_range <- 1
        
        #<march_2015>
        #we know it's dead, thus it must have died in one of two ranges.
				if (format(HART5[i,]$e_dd, format="%m") %in% c("05", "06", "07", "08", "09", "10"))
				{
				  HART5[i,]$died_in_range_5_10 <- 1
				}
				else if (format(HART5[i,]$e_dd, format="%m") %in% c("11", "12", "01", "02", "03", "04"))
				{
          #strictly speaking, this could be just an else, but may as well check.
				  HART5[i,]$died_in_range_11_4 <- 1
				}
				#</march_2015>
			}
		}
	}
}

HART5

HART6 <- HART5#  

days_surv_rate <- 1- sum(HART6$died_in_range)/sum(HART6$days_intersect) 
sum(HART6$days_intersect)
sum(HART6$died_in_range)
annual <- days_surv_rate ^ 364.25
annual
nobou <- length(HART6$days_intersect[HART6$days_intersect>=1]) # no. bou collared
nobou 

# split by the 2 RISK periods


days_surv_rateS <- 1- sum(HART6$died_in_range_5_10)/sum(HART6$days_intersect_5_10)  # S means summer
sum(HART6$days_intersect_5_10)
sum(HART6$died_in_range_5_10)
S <- days_surv_rateS ^ 184   # SUMMER SURVIVAL RATE; 184 IS NUMBER OF DAYS IN SUMMER RISK PERIOD
S
nobou <- length(HART6$days_intersect[HART6$days_intersect>=1]) # no. bou collared
nobou 

# W means winter

days_surv_rateW <- 1- sum(HART6$died_in_range_11_4)/sum(HART6$days_intersect_11_4)  # W means winter
sum(HART6$days_intersect_11_4)
sum(HART6$died_in_range_11_4)
W <- days_surv_rateW ^ 181.25 # WINTER SURVIVAL RATE ; 181.25 IS NUMBER OF DAYS IN WINTER  RISK PERIOD
W

annual_2per <- W*S          # WINTER SURVIVAL RATE X SUMMER SURVIVAL RATE
annual_2per                # final annual suvival value for 2 risk periods
annual



###########################################################
###########################################################
###########################