Merge branch 'dev' into UnitTests

.gitignore

@@ -25,3 +25,4 @@ tests/testthat/old dat files/*
 docs
 doc
 Meta
+Rpath.Rproj

R/ecopath.R

@@ -121,7 +121,7 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
 
   #Special case where B and EE are known then need to solve for BA
   #living[BEE == 1, b := b - (Biomass * PB * EE)]
-  #living[BEE  == 1, diag.a := 0] #Need to wrk on this solution
+  #living[BEE  == 1, diag.a := 0] #Need to work on this solution
 
   A       <- matrix(0, nliving, nliving)
   diag(A) <- living[, diag.a]
@@ -323,20 +323,24 @@ return(path.model)
 rpath.stanzas <- function(Rpath.params){
   #Need to define variables to eliminate check() note about no visible binding
   StGroupNum <- First <- StanzaNum <- VBGF_d <- VBGF_Ksp <- Last <- GroupNum <- NULL
-  WageS <- age <- QageS <- Survive <- Z <- survive_L <- bs.num <- qs.num <- Leading <- NULL
+  WageS <- age <- QageS <- Survive <- Z <- survive_L <- bs.num <- qs.num <- Leading <- Oldest <- NULL
   Group <- Biomass <- R <- NageS <- bs.denom <- bs <- qs.denom <- qs <- Cons <- NULL
   QB <- BAB <- Ex <- NULL
-  
+
   #Determine the total number of groups with multistanzas
   Nsplit     <- Rpath.params$stanza$NStanzaGroups
   groupfile  <- Rpath.params$stanza$stgroups
   stanzafile <- Rpath.params$stanza$stindiv
 
+  # Add Oldest column so age groups will be properly calculated
+  stanzafile[Last == max(Last), Oldest := T, by = StGroupNum]
+  stanzafile[is.na(Oldest), Oldest := F]
+
   #Need to add vector of stanza number
   lastmonth <- rep(NA,Nsplit)
   for(isp in 1:Nsplit){
     #Put the stanzas in order for each split species
-    stnum <- order(stanzafile[StGroupNum == isp, First])
+    stnum <- order(stanzafile[StGroupNum == isp, First])    
     stanzafile[StGroupNum == isp, StanzaNum := stnum]
 
     #Calculate the last month for the final ("leading") stanza
@@ -350,7 +354,7 @@ rpath.stanzas <- function(Rpath.params){
     #(maybe data table has a better way...)
     stmax <- max(stanzafile[StGroupNum == isp, StanzaNum])
     st <- stanzafile[StGroupNum == isp & StanzaNum==stmax,]
-    
+
     gp <- groupfile[isp,]
     #Max age class in months should be one less than a multiple of 12
     #(trying 5999 - probably overkill but for safety)
@@ -370,16 +374,17 @@ rpath.stanzas <- function(Rpath.params){
   groupfile[, last := lastmonth]
 
   for(isp in 1:Nsplit){
-    nstanzas <- groupfile[StGroupNum == isp, nstanzas]
-    stanzafile[StGroupNum == isp & Leading, Last := lastmonth[isp]]
+    nstanzas <- groupfile[
+      StGroupNum == isp, nstanzas]
+    stanzafile[StGroupNum == isp & Oldest, Last := lastmonth[isp]]
 
     #Grab ecopath group codes
     group.codes <- stanzafile[StGroupNum == isp, GroupNum]
 
     #Grab index for first and last months for stanzas
     first   <- stanzafile[StGroupNum == isp, First]
     second  <- stanzafile[StGroupNum == isp, Last]
-        
+
     #Calculate weight and consumption at age    
     StGroup <- data.table(age = stanzafile[StGroupNum == isp & StanzaNum == 1, First]:
                             lastmonth[isp])

src/ecosim.cpp

@@ -156,8 +156,7 @@ List rk4_run (List params, List instate, List forcing, List fishing, List stanza
         // RSK added forced biomass logic
         NumericVector bforce = force_bybio((y-1) * STEPS_PER_YEAR + m, _);
         cur_Biomass = ifelse(bforce>B_BaseRef * EPSILON, bforce, cur_Biomass);
-
-
+
         // KYA 8/9/17 one of the NA or NaN flags is reading back as a negative integer (-2^32)
         // Not sure why.  This sets any negative biomass (assuming this means NaN) to NA_REAL
         cur_Biomass = ifelse((cur_Biomass<0),NA_REAL,cur_Biomass);
@@ -621,7 +620,20 @@ int sp, links, prey, pred, gr, egr, dest, isp, ist, ieco;
    ActiveRespLoss = FoodGain  * ActiveRespFrac  * force_byactresp(dd,_);  												 
    MzeroLoss      = MzeroMort * state_Biomass;
 
-   NumericVector NetProd = FoodGain - UnAssimLoss - ActiveRespLoss - MzeroLoss - FoodLoss;
+
+   // Add mortality forcing
+   for (int i=1; i<=NUM_DEAD+NUM_LIVING; i++){
+     FoodLoss[i]  *= force_bymort(dd, i);
+     MzeroLoss[i] *= force_bymort(dd, i);
+   }
+
+   // Add migration forcing
+   MigrateLoss = clone(state_Biomass);
+   for (int i=1; i<=NUM_DEAD+NUM_LIVING; i++){
+     MigrateLoss[i]  *= force_bymigrate(dd, i);
+   }
+
+   NumericVector NetProd = FoodGain - UnAssimLoss - ActiveRespLoss - MzeroLoss - FoodLoss - MigrateLoss;
 
 // FISHING FUNCTIONS (multiple options depending on fishing method)
 
@@ -696,10 +708,9 @@ int sp, links, prey, pred, gr, egr, dest, isp, ist, ieco;
     NumericVector FORCE_F = (NumericVector)FORCED_FRATE(y,_);
     //  Special "CLEAN" fisheries assuming q=1, so specified input is Frate
         for (sp=1; sp<=NUM_LIVING+NUM_DEAD; sp++){
-          // RSK testing...
-          // caught = FORCED_CATCH(y, sp) + FORCE_F[sp] * state_Biomass[sp];
-          caught =                          FORCE_F[sp] * state_Biomass[sp];
-          // KYA Aug 2011 removed terminal effort option to allow negative fishing pressure 
+
+             caught = FORCED_CATCH(y, sp) + FORCE_F[sp] * state_Biomass[sp];
+             // KYA Aug 2011 removed terminal effort option to allow negative fishing pressure 
                 // if (caught <= -EPSILON) {caught = TerminalF[sp] * state_Biomass[sp];}
              // KYA 10/6/17 Added productivity to Biomass limit for F>1 species (salmon inspired)
                 if (caught >= state_Biomass[sp] + NetProd[sp]){caught = (1.0 - EPSILON) * (state_Biomass[sp] + NetProd[sp]);}
@@ -753,6 +764,7 @@ int sp, links, prey, pred, gr, egr, dest, isp, ist, ieco;
    for (int i=1; i<=NUM_DEAD+NUM_LIVING; i++){
      MigrateLoss[i]  *= force_bymigrate(dd, i);
    }
+
 
 // Sum up derivitive parts (vector sums)
 // Override for group 0 (considered "the sun", never changing)