Added SSB and recruit monthly output (first cut)

Author: kerim aydin

R/ecosim.R

@@ -165,6 +165,12 @@ rsim.run <- function(Rsim.scenario, method = 'RK4', years = 1:100) {
   colnames(rout$annual_QB) <- sps
   colnames(rout$annual_Qlink)<-1:(length(rout$annual_Qlink[1,]))
 
+  colnames(rout$out_SSB)  <- scene$stanzas$Oldest
+  colnames(rout$out_eggs) <- scene$stanzas$Oldest
+  colnames(rout$out_Ninf) <- scene$stanzas$Oldest
+  colnames(rout$out_Winf) <- scene$stanzas$Oldest
+  colnames(rout$out_Nrec) <- scene$stanzas$Groups       
+  colnames(rout$out_Wrec) <- scene$stanzas$Groups   
   # drop the last row (should be always 0; negative index is entry to drop)
   #lastone <- length(rout$annual_Catch[,1])
   #rout$annual_Catch <-    rout$annual_Catch[-lastone,]
@@ -772,6 +778,7 @@ rsim.stanzas <- function(Rpath.params, state, params){
     rstan$Nsplit      <- juvfile$NStanzaGroups
     #Need leading zeros (+1 col/row) to make indexing in C++ easier
     rstan$Nstanzas    <- c(0, juvfile$stgroups$nstanzas)
+    rstan$Totstanzas  <- sum(rstan$Nstanzas)
     rstan$EcopathCode <- matrix(NA, rstan$Nsplit + 1, max(rstan$Nstanzas) + 1)
     rstan$Age1        <- matrix(NA, rstan$Nsplit + 1, max(rstan$Nstanzas) + 1)
     rstan$Age2        <- matrix(NA, rstan$Nsplit + 1, max(rstan$Nstanzas) + 1)
@@ -866,6 +873,23 @@ rsim.stanzas <- function(Rpath.params, state, params){
     rstan$RscaleSplit    <- c(0, rep(1, rstan$Nsplit))
     #sPred1<-sPred+1; rstan$stanzaPred     <- sPred1-1
     sPred1<-sPred+1; rstan$baseStanzaPred <- sPred1-1
+
+    # KYA 12/5/25 - added max.age to sim stanza object for indexing ssb during run
+    rstan$Oldest <- c("Outside", rep(NA, rstan$Nsplit))
+    rstan$Groups <- c("Outside", rep(NA, rstan$Totstanzas))
+    rstan$MaxAge     <- c(0, rep(0,  rstan$Nsplit))
+    sind <- 1
+    for(isp in 1:rstan$Nsplit){
+      for(ist in 1:rstan$Nstanzas[isp + 1]){
+        sind <- sind + 1
+        ieco <- 1 + rstan$EcopathCode[isp+1, ist+1] 
+        rstan$MaxAge[isp+1] <- max(rstan$MaxAge[isp+1], rstan$Age2[isp+1,ist+1])
+        if (rstan$MaxAge[isp+1] == rstan$Age2[isp+1,ist+1]){
+          rstan$Oldest[isp+1] <- params$spname[ieco]
+        }
+        rstan$Groups[sind] <- params$spname[ieco]
+      }
+    }
 
     # SplitSetPred(rstan, state)
   }
@@ -875,6 +899,7 @@ rsim.stanzas <- function(Rpath.params, state, params){
     rstan$Nsplit      <- juvfile$NStanzaGroups
     #Need leading zeros (+1 col/row) to make indexing in C++ easier
     rstan$Nstanzas       <- c(0, 0)
+    rstan$Totstanzas     <- 0
     rstan$EcopathCode    <- matrix(rep(0, 4), 2, 2)
     rstan$Age1           <- matrix(rep(0, 4), 2, 2)
     rstan$Age2           <- matrix(rep(0, 4), 2, 2)
@@ -902,7 +927,9 @@ rsim.stanzas <- function(Rpath.params, state, params){
     rstan$baseSpawnBio   <- c(0, 0)
     rstan$RscaleSplit    <- c(0, 0)
     rstan$baseStanzaPred <- c(0, 0)
-    
+    rstan$Oldest <- c("Outside", "None")
+    rstan$Groups <- c("Outside", "None")
+    rstan.MaxAge         <- c(0,0)
   }
 
 

src/ecosim.cpp

@@ -224,7 +224,7 @@ List Adams_run (List params, List instate, List forcing, List fishing, List stan
                  int StartYear, int EndYear, List InitDeriv){
 
 int y, m, dd; 
-
+//std::cout << " x1 ";
 // Get some basic needed numbers from the params List
    const int NUM_BIO = as<int>(params["NUM_LIVING"]) + as<int>(params["NUM_DEAD"]);
    const int NumPredPreyLinks          = as<int>(params["NumPredPreyLinks"]);
@@ -254,15 +254,27 @@ int y, m, dd;
 
 // Parameters from stanzas
    const int Nsplit         = as<int>(stanzas["Nsplit"]);
-
+   const int Totstanzas     = as<int>(stanzas["Totstanzas"]);
+   const NumericVector Nstanzas = as<NumericVector>(stanzas["Nstanzas"]);
+   const NumericVector MaxAge   = as<NumericVector>(stanzas["MaxAge"]);
+   NumericMatrix Age1           = as<NumericMatrix>(stanzas["Age1"]);
+   //std::cout << " x1a ";
+   // stanza outputs
+      NumericMatrix out_SSB(EndYear * 12, Nsplit + 1);
+      NumericMatrix out_eggs(EndYear * 12, Nsplit + 1);
+      NumericMatrix out_Winf(EndYear * 12, Nsplit + 1);
+      NumericMatrix out_Ninf(EndYear * 12, Nsplit + 1);
+      NumericMatrix out_Wrec(EndYear * 12, Totstanzas + 1);
+      NumericMatrix out_Nrec(EndYear * 12, Totstanzas + 1);      
+      //std::cout << " x1b ";   
 // Parameter need to track catch by Gear
    const NumericVector FishFrom = as<NumericVector>(params["FishFrom"]);
 
 // Monthly output matrices                     
    NumericMatrix out_Biomass( EndYear * 12, NUM_BIO + 1);           
    NumericMatrix out_Catch( EndYear * 12, NUM_BIO + 1);          
-   NumericMatrix out_SSB(EndYear * 12, NUM_BIO + 1);        
-   NumericMatrix out_rec(EndYear * 12, NUM_BIO + 1);
+   //NumericMatrix out_SSB(EndYear * 12, NUM_BIO + 1);        
+   //NumericMatrix out_rec(EndYear * 12, NUM_BIO + 1);
    NumericMatrix out_Gear_Catch(EndYear*12, NumFishingLinks+1);
 // Annual output matrices
    NumericMatrix annual_Catch(EndYear, NUM_BIO+1);
@@ -274,7 +286,12 @@ int y, m, dd;
 // Use Clone to make sure state/stanzas are copies of instate/instanzas, not pointers   
    List state = clone(instate);
    //List stanzas = clone(instanzas);
-   
+   //std::cout << " x1c "; 
+   NumericMatrix NageS          = as<NumericMatrix>(state["NageS"]);
+   NumericMatrix WageS          = as<NumericMatrix>(state["WageS"]);
+   NumericVector SpawnBio       = as<NumericVector>(state["SpawnBio"]);
+   NumericVector EggsStanza     = as<NumericVector>(state["EggsStanza"]);
+   //std::cout << " x1d "; 
    // Update sums of split groups to total biomass for derivative calcs
    if(Nsplit > 0){
      SplitSetPred(stanzas, state); 
@@ -284,7 +301,7 @@ int y, m, dd;
       List dyt = InitDeriv;
 
    dd = StartYear * STEPS_PER_YEAR;
-
+   //std::cout << "x2";
 // MAIN LOOP STARTS HERE
 // ASSUMES STEPS_PER_MONTH will always be 1.0, took out divisions     
    for (y = StartYear; y <= EndYear; y++){
@@ -392,8 +409,8 @@ int y, m, dd;
 
      // Write to output matricies     				          									                    
         out_Biomass( dd, _) = old_Biomass;
-        out_SSB(dd, _) = old_Biomass;
-        out_rec(dd, _) = old_Biomass;
+        //out_SSB(dd, _) = old_Biomass;
+        //out_rec(dd, _) = old_Biomass;
         out_Catch( dd, _) = new_Catch;
         out_Gear_Catch(dd, _) = new_Gear_Catch;
         annual_Catch(y-1, _) = annual_Catch(y-1, _) + new_Catch;
@@ -402,12 +419,22 @@ int y, m, dd;
           annual_QB(y-1, _)    = FoodGain/old_Biomass;
           annual_Qlink(y-1, _) = Qlink;
         }
-        
-     //NOJUV    for (i = 1; i <= juv_N; i++){
-     //NOJUV    out_SSB(dd, JuvNum[i]) = 0.0;
-     //NOJUV 		out_SSB(dd, AduNum[i]) = SpawnBio[i];
-     //NOJUV 		out_rec(dd, AduNum[i]) = NageS(firstMoAdu[i], i) * WageS(firstMoAdu[i], i);
-
+        //std::cout << "x3";      
+     // Write stanza outputs
+        int sind, isp, ist;
+        sind = 0;
+        for (isp=1; isp<=Nsplit; isp++){
+          out_Winf(dd,isp) = WageS(MaxAge[isp], isp);
+          out_Ninf(dd,isp) = NageS(MaxAge[isp], isp); 
+          out_SSB(dd,isp)  = SpawnBio[isp];
+          out_eggs(dd,isp) = EggsStanza[isp];  
+          for (ist=1; ist<=Nstanzas[isp]; ist++){
+            sind++;
+            out_Nrec(dd,sind) = NageS(Age1(isp,ist), isp);
+            out_Wrec(dd,sind) = WageS(Age1(isp,ist), isp);
+          }
+        }
+        //std::cout << "x4";
      }  // End of main months loop
 
    }// End of years loop
@@ -432,7 +459,13 @@ int y, m, dd;
      _["annual_Catch"]=annual_Catch,
      _["annual_Biomass"]=annual_Biomass,
      _["annual_QB"]=annual_QB,
-     _["annual_Qlink"]=annual_Qlink,     
+     _["annual_Qlink"]=annual_Qlink,
+     _["out_SSB"]=out_SSB,
+     _["out_eggs"]=out_eggs,
+     _["out_Winf"]=out_Winf,
+     _["out_Ninf"]=out_Ninf,
+     _["out_Nrec"]=out_Nrec,
+     _["out_Wrec"]=out_Wrec,
      _["end_state"]=state,
      _["crash_year"]=CRASH_YEAR,
      _["dyt"]=dyt);