#!/usr/bin/env Rscript

#suppressPackageStartupMessages(library("optparse"))
#suppressPackageStartupMessages(library("rjson"))

options_list <- list(
    make_option(c("-s", "--save_log"), action="store_true", default=FALSE, help="Save R logs"),
    make_option(c("-m", "--output_r_log"), action="store", help="Output dataset for R logs"),
    make_option(c("-o", "--output"), action="store", help="Output dataset")
)

parser <- OptionParser(usage="%prog [options] file", options_list)
args <- parse_args(parser, positional_arguments=TRUE)
opt <- args$options


daylength = function(L){
    # from Forsythe 1995
    p = 0.8333
    dl <- NULL
    for (i in 1:365) {
        theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186)))
        phi <- asin(0.39795 * cos(theta))
        dl[i] <- 24 - 24/pi * acos((sin(p * pi/180) + sin(L * pi/180) * sin(phi))/(cos(L * pi/180) * cos(phi)))
    }
    # return a vector of daylength in 365 days
    dl
}


# source("daylength.R")
hourtemp = function(L,date){
    # L = 37.5 specify this in main program
    # base development threshold for BMSB
    threshold <- 12.7
    # threshold2 <- threshold/24 degree hour accumulation
    #expdata <- tempdata[1:365,11:13] # Use daily max, min, mean
    # daily minimum
    dnp <- expdata[date,2]
    # daily maximum
    dxp <- expdata[date,3]
    dmean <- 0.5 * (dnp + dxp)
    #if (dmean>0) {
        #dnp <- dnp - k1 * dmean  
        #dxp <- dxp + k2 * dmean 
    #} else {
        #dnp <- dnp + k1 * dmean  
        #dxp <- dxp - k2 * dmean
    #}
    dd <- 0  # initialize degree day accumulation

    if (dxp<threshold) {
        dd <- 0
    } else {
        # extract daylength data for entire year
        dlprofile <- daylength(L)
        # initialize hourly temperature
        T <- NULL
        #initialize degree hour vector
        dh <- NULL
        # calculate daylength in given date
        # date <- 200
        y <- dlprofile[date]
        # night length
        z <- 24 - y
        # lag coefficient
        a <- 1.86
        # night coefficient
        b <- 2.20
        #import raw data set
        #tempdata <- read.csv("tempdata.csv")
        # Should be outside function otherwise its redundant
        # sunrise time
        risetime <- 12 - y/2
        # sunset time
        settime <- 12 + y/2
        ts <- (dxp - dnp) * sin(pi * (settime-5)/(y + 2 * a)) + dnp
        for (i in 1:24) {
            if (i > risetime && i < settime) {
                # number of hours after Tmin until sunset
                m <- i - 5
                T[i] = (dxp - dnp) * sin(pi * m/(y + 2 * a)) + dnp
                if (T[i]<8.4) {
                    dh[i] <- 0
                } else {
                    dh[i] <- T[i] - 8.4
                }
            } else if (i>settime) {
                n <- i - settime
                T[i] = dnp + (ts - dnp) * exp(-b * n/z)
                if (T[i]<8.4) {
                    dh[i] <- 0
                } else {
                    dh[i] <- T[i] - 8.4
                }
            } else {
                n <- i + 24 - settime
                T[i] = dnp + (ts - dnp) * exp(-b * n / z)
                if (T[i]<8.4) {
                    dh[i] <- 0
                } else {
                    dh[i] <- T[i] - 8.4
                }
            }
        }
        dd <- sum(dh) / 24
    }
    return = c(dmean, dd)
    return
}


mortality.egg = function(temperature) {
    if (temperature < 12.7) {
        mort.prob = 1
    } else {
        # 100% mortality if <12.7
        mort.prob = 0.8 - temperature / 40
        if (mort.prob<0) {
            mort.prob = 0.01
        }
    }
    return = mort.prob
    return
}


mortality.nymph = function(temperature) {
    if (temperature<12.7) {
        mort.prob = 0.03
    } else {
        # at low temperature
        mort.prob = -temperature * 0.0008 + 0.03
    }
    return = mort.prob
    return
}


mortality.adult = function(temperature) {
    if (temperature < 12.7) {
        mort.prob = 0.002
    } else {
        mort.prob = -temperature * 0.0005 + 0.02
    }
    return = mort.prob
    return
}


# model initialization
# TODO: add tool params for the following options.
# start with 1000 individuals
n <- 1000
# Generation, Stage, DD, T, Diapause
vec.ini <- c(0,3,0,0,0)
# overwintering, previttelogenic, DD = 0, T = 0, no-diapause
vec.mat <- rep(vec.ini,n)
# complete matrix for the population
vec.mat <- t(matrix(vec.mat, nrow = 5))
# latitude for Asheville NC
L <- 35.58
# complete photoperiod profile in a year, requires daylength function
ph.p <- daylength(L)

# load temperature data@location/year
load("asheville2014.Rdat")

# time series of population size
tot.pop <- NULL

# gen.0 pop size
gen0.pop <- rep(0, 365)
gen1.pop <- rep(0, 365)
gen2.pop <- rep(0, 365)

# aggregate
S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, 365)
g0.adult <- g1.adult <- g2.adult <- rep(0, 365)

# birth death adults
N.newborn <- N.death <- N.adult <- rep(0, 365)

# degree-day accumulation
dd.day <- rep(0, 365)

# start tick
ptm  <-  proc.time()

for (n.sim in 1:1000) {
    # loop through 1000 simulations
    for (day in 1:365) {
        # loop through 365 day/yr
        photoperiod <- ph.p[day] 
        # photoperiod in the day
        temp.profile <- hourtemp(L,day) 
        # temperature profile
        mean.temp <- temp.profile[1] 
        # mean temp
        dd.temp <- temp.profile[2] 
        # degree-day
        dd.day[day] <- dd.temp
        death.vec <- NULL
        # trash bin for death
        birth.vec <- NULL
        # record new born
        for (i in 1:n) { 
            # loop through all individual
            vec.ind <- vec.mat[i,]
            # find individual record
            # first of all, still alive?
            if (vec.ind[2] == 0) {
                # egg
                death.prob = mortality.egg(mean.temp)
            } else if (vec.ind[2] == 1 | vec.ind[2] == 2) {
                # nymph
                death.prob = mortality.nymph(mean.temp)
            }  else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) {
                # for adult
                death.prob = mortality.adult(mean.temp) 
            }
            u.d <- runif(1)
            if (u.d<death.prob) {
                death.vec <- c(death.vec,i)
            } else {
                # aggregrate index of dead bug
                # event 1 end of diapause
                if (vec.ind[1] == 0 && vec.ind[2] == 3) {
                    # overwintering adult (previttelogenic)
                    if (photoperiod>13.5 && vec.ind[3] > 68 && day < 180) {
                        # add 68C to become fully reproductively matured
                        # transfer to vittelogenic
                        vec.ind <- c(0,4,0,0,0)
                        vec.mat[i,] <- vec.ind
                    } else {
                         # add to DD
                        vec.ind[3] <- vec.ind[3] + dd.temp
                        vec.ind[4] <- vec.ind[4] + 1 # add 1 day in current stage
                        vec.mat[i,] <- vec.ind
                    }
                }
                if (vec.ind[1]!=0 && vec.ind[2] == 3) {
                    # NOT overwintering adult (previttelogenic)
                    current.gen <- vec.ind[1]
                    if (vec.ind[3]>68) {
                        # add 68C to become fully reproductively matured
                        # transfer to vittelogenic
                        vec.ind <- c(current.gen,4,0,0,0)
                       vec.mat[i,] <- vec.ind
                    } else {
                        # add to DD
                        vec.ind[3] <- vec.ind[3] + dd.temp
                        # add 1 day in current stage
                        vec.ind[4] <- vec.ind[4] + 1
                        vec.mat[i,] <- vec.ind
                    }
                }
                # event 2 oviposition -- where population dynamics comes from
                # vittelogenic stage, overwintering generation
                if (vec.ind[2] == 4 && vec.ind[1] == 0 && mean.temp>10) {
                    if (vec.ind[4] == 0) {
                        # just turned in vittelogenic stage
                        n.birth = round(runif(1,10,20))
                    } else {
                        p.birth = 1/4/75 
                        # prob of birth
                        u1 <- runif(1)
                        if (u1<p.birth) {
                            n.birth = n.birth
                        }
                    }
                    # add to DD
                    vec.ind[3] <- vec.ind[3] + dd.temp
                    # add 1 day in current stage
                    vec.ind[4] <- vec.ind[4] + 1
                    vec.mat[i,] <- vec.ind
                    if (n.birth>0) {
                        # add new birth -- might be in different generations
                        # generation  + 1
                        new.gen <- vec.ind[1] + 1
                        # egg profile
                        new.ind <- c(new.gen,0,0,0,0)
                        new.vec <- rep(new.ind,n.birth)
                        # update batch of egg profile
                        new.vec <- t(matrix(new.vec,nrow = 5))
                        # group with total eggs laid in that day
                        birth.vec <- rbind(birth.vec,new.vec)
                    }
                }
                # event 2 oviposition -- for gen 1.
                if (vec.ind[2] == 4 && vec.ind[1] == 1 && mean.temp>12.5 && day<222) {
                    # vittelogenic stage, 1st generation
                    if (vec.ind[4] == 0) {
                        # just turned in vittelogenic stage
                        n.birth = round(runif(1,10,20))
                    } else {
                        p.birth = 1/4/75 
                        u1 <- runif(1)
                        if (u1<p.birth) {
                            n.birth = n.birth
                        }
                    }
                    # add to DD
                    vec.ind[3] <- vec.ind[3] + dd.temp
                    # add 1 day in current stage
                    vec.ind[4] <- vec.ind[4] + 1
                    vec.mat[i,] <- vec.ind
                    if (n.birth>0) {
                        # add new birth -- might be in different generations
                        # generation + 1
                        new.gen <- vec.ind[1] + 1
                        # egg profile
                        new.ind <- c(new.gen,0,0,0,0)
                        new.vec <- rep(new.ind,n.birth)
                        # update batch of egg profile
                        new.vec <- t(matrix(new.vec,nrow = 5))
                        # group with total eggs laid in that day
                        birth.vec <- rbind(birth.vec,new.vec)
                    }
                }
                # event 3 development (with diapause determination)
                # event 3.1 egg development to young nymph (vec.ind[2] = 0 -> egg)
                # egg stage
                if (vec.ind[2] == 0) {
                    # add to DD
                    vec.ind[3] <- vec.ind[3] + dd.temp
                    # from egg to young nymph
                    if (vec.ind[3] >= 53.30 && -0.9843 * dd.temp + 33.438>0) {
                        current.gen <- vec.ind[1]
                        # transfer to young nym stage
                        vec.ind <- c(current.gen,1,0,0,0)
                    } else {
                        # add 1 day in current stage
                        vec.ind[4] <- vec.ind[4] + 1
                    }
                    vec.mat[i,] <- vec.ind
                }
                # event 3.2 young nymph to old nymph (vec.ind[2] = 1 -> young nymph: determines diapause)
                # young nymph stage
                if (vec.ind[2] == 1) {
                    # add to DD
                    vec.ind[3] <- vec.ind[3] + dd.temp
                    # from young to old nymph
                    if (vec.ind[3] >= 537/2 && -0.45 * dd.temp + 18>0) {
                        current.gen <- vec.ind[1]
                        # transfer to old nym stage
                        vec.ind <- c(current.gen,2,0,0,0)
                        # prepare for diapausing
                        if (photoperiod<13.5 && day > 180) {
                            vec.ind[5] <- 1
                        }
                    } else {
                        # add 1 day in current stage
                        vec.ind[4] <- vec.ind[4] + 1
                    }
                    vec.mat[i,] <- vec.ind
                }
                # event 3.3 old nymph to adult: previttelogenic or diapausing?
                # old nymph stage
                if (vec.ind[2] == 2) {
                    # add to DD
                    vec.ind[3] <- vec.ind[3] + dd.temp
                    # from old to adult
                    if (vec.ind[3] >= 537/2 && -0.50 * dd.temp + 22>0) {
                        current.gen <- vec.ind[1]
                        # non-diapausing adult -- previttelogenic
                        if (vec.ind[5] == 0) {
                            vec.ind <- c(current.gen,3,0,0,0)
                        # diapausing 
                        } else {
                            vec.ind <- c(current.gen,5,0,0,1)
                        }
                    } else {
                        # add 1 day in current stage
                        vec.ind[4] <- vec.ind[4] + 1
                    }
                    vec.mat[i,] <- vec.ind
                }
                # event 4 growing of diapausing adult (unimportant, but still necessary)## 
                if (vec.ind[2] == 5) {
                    vec.ind[3] <- vec.ind[3] + dd.temp
                    vec.ind[4] <- vec.ind[4] + 1
                    vec.mat[i,] <- vec.ind
                }
            } # else if it is still alive
        } # end of the individual bug loop
        # find how many died
        n.death <- length(death.vec)
        if (n.death>0) {
            vec.mat <- vec.mat[-death.vec, ]
        }
        # remove record of dead
        # find how many new born  
        n.newborn <- length(birth.vec[,1])
        vec.mat <- rbind(vec.mat,birth.vec)
        # update population size for the next day
        n <- n-n.death + n.newborn 

        # aggregate results by day
        tot.pop <- c(tot.pop,n)
        # egg
        s0 <- sum(vec.mat[,2] == 0)
        # young nymph
        s1 <- sum(vec.mat[,2] == 1)
        # old nymph
        s2 <- sum(vec.mat[,2] == 2)
        # previtellogenic
        s3 <- sum(vec.mat[,2] == 3)
        # vitellogenic
        s4 <- sum(vec.mat[,2] == 4)
        # diapausing
        s5 <- sum(vec.mat[,2] == 5)
        # overwintering adult
        gen0 <- sum(vec.mat[,1] == 0)
        # first generation
        gen1 <- sum(vec.mat[,1] == 1)
        # second generation
        gen2 <- sum(vec.mat[,1] == 2)
        # sum of all adults
        n.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5)
        # gen.0 pop size
        gen0.pop[day] <- gen0
        gen1.pop[day] <- gen1
        gen2.pop[day] <- gen2
        S0[day] <- s0
        S1[day] <- s1
        S2[day] <- s2
        S3[day] <- s3
        S4[day] <- s4
        S5[day] <- s5
        g0.adult[day] <- sum(vec.mat[,1] == 0)
        g1.adult[day] <- sum((vec.mat[,1] == 1 & vec.mat[,2] == 3) | (vec.mat[,1] == 1 & vec.mat[,2] == 4) | (vec.mat[,1] == 1 & vec.mat[,2] == 5))
        g2.adult[day] <- sum((vec.mat[,1] == 2 & vec.mat[,2] == 3) | (vec.mat[,1] == 2 & vec.mat[,2] == 4) | (vec.mat[,1] == 2 & vec.mat[,2] == 5))
        N.newborn[day] <- n.newborn
        N.death[day] <- n.death
        N.adult[day] <- n.adult
    }   
#print(n.sim)
}

proc.time() - ptm
dd.cum <- cumsum(dd.day)
save(dd.day, dd.cum, S0, S1, S2, S3, S4, S5, N.newborn, N.death, N.adult, tot.pop, gen0.pop, gen1.pop, gen2.pop, g0.adult, g1.adult, g2.adult, file=opt$output)