## This code contains the input parameters figures
source("branching_functions.R")
library(ggplot2)
library(reshape2)
library(popbio)


## Survival
sCaveBat = c(0.75, 0.95)
## probability of "n" births (0, 1, 2, ... N)
## Rows should sum to one
pCaveBat = matrix(
    c(0.9, 0.1,
      0.60, 0.40),
    byrow = TRUE,
    nrow = 2)

## Years to simulate and mortality vector
nYears = 30

## output matrix and loop runing model over mu values
## projection matrix
pCaveBat
projCaveBat = matrix(c(
    pCaveBat[,2],
    sCaveBat),
    byrow = TRUE, nrow = 2, ncol = 2)

popCaveBat <- matrix(0, nrow = 2, ncol = (nYears +1))
popCaveBat[,1] <- c(10, 20)
for(year in 1:nYears){
    popCaveBat[,year + 1] <- projCaveBat %*% popCaveBat[,year]
}

popCaveBat
projCaveBat
generation.time(projCaveBat)
caveBatEigen <- round(eigen(projCaveBat)$values,4)
caveBatEigen

## Tree bat assessment
pTreeBat = matrix(
    c(16/24, 4/24, 3/24, 1/24, 0/24, 0/24,
      14/24, 1/24, 3/24, 3/24, 2/24, 1/24),
    byrow = TRUE,
    nrow = 2)

apply(pTreeBat, 1, sum)
pTreeBat
sTreeBat = c(0.65, 0.85)

## projection matrix
projTreeBat = matrix(c(
    pTreeBat %*% 0:5,
    sTreeBat),
    byrow = TRUE, nrow = 2, ncol = 2)

projTreeBat

round(eigen(projTreeBat)$values,4)
round(eigen(projCaveBat)$values,4)
treeBatEigen <- eigen(projCaveBat)$values



## Compare 1-stage vs 2-stages
## Two Stage Passerine
##sPassTwo = c(0.180, 0.57)
sPassTwo = c(0.50, 0.55)
## probability of "n" births (0, 1, 2, ... N)
## Rows should sum to one
x = 0:8
pPassTwo = matrix(
    c(round(dpois(x = x, lambda = 1.3), 4),
      round(dpois(x = x, lambda = 1.8), 4)),
    byrow = TRUE,
    nrow = 2)

pPassTwo
apply(pPassTwo, 1, sum)
pPassTwo %*% x

## Years to simulate and mortality vector
nYears = 30
mu = seq(0,1, by = 0.01)


## projection matrix
pPassTwo
projPassTwo = matrix(c(
    pPassTwo %*% 0:8,
    sPassTwo),
    byrow = TRUE, nrow = 2, ncol = 2)

projPassTwo
generation.time(projPassTwo)

passTwoEigen <- round(eigen(projPassTwo)$values,4)
passTwoEigen



## Eagle assessment
pEagle = matrix(
    c(c(1, rep(0, 3)),
      c(1, rep(0, 3)),
      c(1, rep(0, 3)),
      c(9/12, 2/12, 1/12, 0),
      c(7/12, 3/12, 1.5/12, 0.5/12)),
    byrow = TRUE,
    nrow = 5)

pEagle
sEagle = c(0.60, 0.95, 0.95, 0.95, 0.90)
projEagle = matrix(
    c(0, 0, 0, sum(pEagle[4,] * 0:3), sum(pEagle[5,] * 0:3),
      sEagle[1], 0, 0, 0, 0,
      0, sEagle[2], 0, 0, 0,
      0, 0, sEagle[3], 0, 0,
      0, 0, 0, sEagle[4], sEagle[5]),
    byrow = TRUE, nrow = 5, ncol = 5)

projEagle

round(eigen(projEagle)$value, 3)
eagleEigen <- eigen(projEagle)$value

colnames(pEagle) <- 0:(dim(pEagle)[2]- 1)
rownames(pEagle) <- paste("Life Stage", 1:dim(pEagle)[1])
pEagleGG <- melt(pEagle)
colnames(pEagleGG) <- c("LifeStage", "Births", "Probability")
pEagleGG$Species <- "Eagle"
ggplot(data = pEagleGG, aes( x = Births, y = Probability)) +
    geom_bar(stat = 'identity') + facet_grid(LifeStage~.)


colnames(pTreeBat) <- 0:(dim(pTreeBat)[2]- 1)
rownames(pTreeBat) <- paste("Life Stage", 1:dim(pTreeBat)[1])
pTreeBatGG <- melt(pTreeBat)
colnames(pTreeBatGG) <- c("LifeStage", "Births", "Probability")
pTreeBatGG$Species <- "Tree Bat"

colnames(pCaveBat) <- 0:(dim(pCaveBat)[2]- 1)
rownames(pCaveBat) <- paste("Life Stage", 1:dim(pCaveBat)[1])
pCaveBatGG <- melt(pCaveBat)
colnames(pCaveBatGG) <- c("LifeStage", "Births", "Probability")
pCaveBatGG$Species <- "Cave Bat"


pCaveBatGG <- melt(pCaveBat)
colnames(pCaveBatGG) <- c("LifeStage", "Births", "Probability")
pCaveBatGG$Species <- "Cave Bat"

rownames(pPassTwo) <- paste("Life Stage", 1:dim(pPassTwo)[1])
colnames(pPassTwo) <- 0:(dim(pPassTwo)[2]- 1)
pPassGG <- melt(pPassTwo)
colnames(pPassGG) <- c("LifeStage", "Births", "Probability")
pPassGG$Species <- "Grassland\nSongbird"

pPassGG$Survival = sPassTwo
pEagleGG$Survival = sEagle
pTreeBatGG$Survival = sTreeBat
pCaveBatGG$Survival = sCaveBat

pAlls <- rbind(pPassGG,
               pCaveBatGG,
               pTreeBatGG,
               pEagleGG)

inputParmBW <-  ggplot(data = pAlls,
                     aes( x = Births, y = Probability,
                         fill = Survival)) +
    geom_bar(stat = 'identity') + facet_grid(LifeStage~Species) +
    theme_bw() + xlab(
        "Number of females born and surviving first year\nper adult female") +
    scale_fill_gradient(low = 'grey60', high = 'black',
                        guide = guide_legend(title = "Annual\nsurvival"))
inputParmBW
inputParmCol <-  ggplot(data = pAlls,
                     aes( x = Births, y = Probability,
                         fill = Survival)) +
    geom_bar(stat = 'identity') + facet_grid(LifeStage~Species) +
    theme_bw() + xlab(
        "Number of females born and surving first year\nper adult female") +
    scale_fill_gradient(low = 'red', high = 'blue',
                        guide = guide_legend(title = "Annual\nsurvival"))
inputParmCol
windows()
print(inputParmCol)

ggsave('inputParmCol.pdf', inputParmCol, width = 6, height = 6)

ggsave('inputParmBW.pdf', inputParmBW, width = 6, height = 6)

#########
## Generation times
generation.time(projCaveBat)
generation.time(projTreeBat)
generation.time(projPassTwo)
generation.time(projEagle)