Skip to content

Commit 9f98598

Browse files
author
kerim aydin
committed
Adding by-species derivative diagnostics to rsim.run
1 parent a3ead6a commit 9f98598

5 files changed

Lines changed: 75 additions & 24 deletions

File tree

‎R/RcppExports.R‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
11
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
22
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
33

4-
rk4_run <- function(params, instate, forcing, fishing, stanzas, StartYear, EndYear) {
5-
.Call('_Rpath_rk4_run', PACKAGE = 'Rpath', params, instate, forcing, fishing, stanzas, StartYear, EndYear)
4+
rk4_run <- function(params, instate, forcing, fishing, stanzas, StartYear, EndYear, spnum, spstanza) {
5+
.Call('_Rpath_rk4_run', PACKAGE = 'Rpath', params, instate, forcing, fishing, stanzas, StartYear, EndYear, spnum, spstanza)
66
}
77

8-
Adams_run <- function(params, instate, forcing, fishing, stanzas, StartYear, EndYear, InitDeriv) {
9-
.Call('_Rpath_Adams_run', PACKAGE = 'Rpath', params, instate, forcing, fishing, stanzas, StartYear, EndYear, InitDeriv)
8+
Adams_run <- function(params, instate, forcing, fishing, stanzas, StartYear, EndYear, InitDeriv, spnum, spstanza) {
9+
.Call('_Rpath_Adams_run', PACKAGE = 'Rpath', params, instate, forcing, fishing, stanzas, StartYear, EndYear, InitDeriv, spnum, spstanza)
1010
}
1111

1212
deriv_vector <- function(params, state, forcing, fishing, stanzas, inyear, m, tt) {

‎R/ecosim.R‎

Lines changed: 21 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -118,9 +118,23 @@ rsim.scenario <- function(Rpath, Rpath.params, years = 1:100){
118118
#'
119119
#'@export
120120
#'
121-
rsim.run <- function(Rsim.scenario, method = 'RK4', years = 1:100) {
121+
rsim.run <- function(Rsim.scenario, method = 'RK4', years = 1:100, spname = NULL) {
122+
123+
scene <- copy(Rsim.scenario)
124+
125+
# Look up C code species number from species name
126+
if (length(spname)>1){stop("Only one name allowed for species name.")}
127+
if(is.null(spname)){
128+
sp.num <- 0
129+
} else {
130+
sp.num <- as.numeric(scene$params$spnum[spname])
131+
if(is.na(sp.num)){sp.num <- 0; warning("Group ",spname," for diagnostics output.")}
132+
}
133+
134+
# Check if the species code is associated with a stanza
135+
sp.stanza <- which(scene$stanzas$EcopathCode %in% sp.num) %% (scene$stanzas$Nsplit+1) - 1
136+
if(length(sp.stanza)!=1){sp.stanza <- 0}
122137

123-
scene <- copy(Rsim.scenario)
124138

125139
# Perform argument checks: Check method name and figure out starting and ending years for run
126140
if (method != 'RK4' && method != 'AB') {
@@ -143,7 +157,7 @@ rsim.run <- function(Rsim.scenario, method = 'RK4', years = 1:100) {
143157
if (method == 'RK4') {
144158
rout <- rk4_run(scene$params, scene$start_state,
145159
scene$forcing, scene$fishing,
146-
scene$stanzas, syear, eyear)
160+
scene$stanzas, syear, eyear, sp.num, sp.stanza)
147161
} else if (method == 'AB') {
148162
#Run initial derivative
149163
derv <- deriv_vector(scene$params, scene$start_state,
@@ -154,7 +168,7 @@ rsim.run <- function(Rsim.scenario, method = 'RK4', years = 1:100) {
154168
#Run Adams-Bashforth Algorithm
155169
rout <- Adams_run(scene$params, scene$start_state,
156170
scene$forcing, scene$fishing,
157-
scene$stanzas, syear, eyear, derv)
171+
scene$stanzas, syear, eyear, derv, sp.num, sp.stanza)
158172
}
159173
# Nicely Name output vectors
160174
sps <- scene$params$spname[1:(1+scene$params$NUM_BIO)]
@@ -927,9 +941,9 @@ rsim.stanzas <- function(Rpath.params, state, params){
927941
rstan$baseSpawnBio <- c(0, 0)
928942
rstan$RscaleSplit <- c(0, 0)
929943
rstan$baseStanzaPred <- c(0, 0)
930-
rstan$Oldest <- c("Outside", "None")
931-
rstan$Groups <- c("Outside", "None")
932-
rstan.MaxAge <- c(0,0)
944+
rstan$Oldest <- c("Outside")
945+
rstan$Groups <- c("Outside")
946+
rstan$MaxAge <- c(0,0)
933947
}
934948

935949

‎R/ecosim_multistep.R‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ rsim.step <- function(Rsim.scenario, Rsim.output, method = 'AB',year.end){
2929
next.run <- Adams_run(Rsim.scenario$params, full.run$end_state,
3030
Rsim.scenario$forcing, Rsim.scenario$fishing,
3131
Rsim.scenario$stanzas, step.start, step.end,
32-
full.run$dyt)
32+
full.run$dyt, 0, 0)
3333
}
3434
#Merge runs
3535
last <- nrow(full.run$out_Biomass)

‎src/RcppExports.cpp‎

Lines changed: 12 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -11,8 +11,8 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
1111
#endif
1212

1313
// rk4_run
14-
List rk4_run(List params, List instate, List forcing, List fishing, List stanzas, int StartYear, int EndYear);
15-
RcppExport SEXP _Rpath_rk4_run(SEXP paramsSEXP, SEXP instateSEXP, SEXP forcingSEXP, SEXP fishingSEXP, SEXP stanzasSEXP, SEXP StartYearSEXP, SEXP EndYearSEXP) {
14+
List rk4_run(List params, List instate, List forcing, List fishing, List stanzas, int StartYear, int EndYear, int spnum, int spstanza);
15+
RcppExport SEXP _Rpath_rk4_run(SEXP paramsSEXP, SEXP instateSEXP, SEXP forcingSEXP, SEXP fishingSEXP, SEXP stanzasSEXP, SEXP StartYearSEXP, SEXP EndYearSEXP, SEXP spnumSEXP, SEXP spstanzaSEXP) {
1616
BEGIN_RCPP
1717
Rcpp::RObject rcpp_result_gen;
1818
Rcpp::RNGScope rcpp_rngScope_gen;
@@ -23,13 +23,15 @@ BEGIN_RCPP
2323
Rcpp::traits::input_parameter< List >::type stanzas(stanzasSEXP);
2424
Rcpp::traits::input_parameter< int >::type StartYear(StartYearSEXP);
2525
Rcpp::traits::input_parameter< int >::type EndYear(EndYearSEXP);
26-
rcpp_result_gen = Rcpp::wrap(rk4_run(params, instate, forcing, fishing, stanzas, StartYear, EndYear));
26+
Rcpp::traits::input_parameter< int >::type spnum(spnumSEXP);
27+
Rcpp::traits::input_parameter< int >::type spstanza(spstanzaSEXP);
28+
rcpp_result_gen = Rcpp::wrap(rk4_run(params, instate, forcing, fishing, stanzas, StartYear, EndYear, spnum, spstanza));
2729
return rcpp_result_gen;
2830
END_RCPP
2931
}
3032
// Adams_run
31-
List Adams_run(List params, List instate, List forcing, List fishing, List stanzas, int StartYear, int EndYear, List InitDeriv);
32-
RcppExport SEXP _Rpath_Adams_run(SEXP paramsSEXP, SEXP instateSEXP, SEXP forcingSEXP, SEXP fishingSEXP, SEXP stanzasSEXP, SEXP StartYearSEXP, SEXP EndYearSEXP, SEXP InitDerivSEXP) {
33+
List Adams_run(List params, List instate, List forcing, List fishing, List stanzas, int StartYear, int EndYear, List InitDeriv, int spnum, int spstanza);
34+
RcppExport SEXP _Rpath_Adams_run(SEXP paramsSEXP, SEXP instateSEXP, SEXP forcingSEXP, SEXP fishingSEXP, SEXP stanzasSEXP, SEXP StartYearSEXP, SEXP EndYearSEXP, SEXP InitDerivSEXP, SEXP spnumSEXP, SEXP spstanzaSEXP) {
3335
BEGIN_RCPP
3436
Rcpp::RObject rcpp_result_gen;
3537
Rcpp::RNGScope rcpp_rngScope_gen;
@@ -41,7 +43,9 @@ BEGIN_RCPP
4143
Rcpp::traits::input_parameter< int >::type StartYear(StartYearSEXP);
4244
Rcpp::traits::input_parameter< int >::type EndYear(EndYearSEXP);
4345
Rcpp::traits::input_parameter< List >::type InitDeriv(InitDerivSEXP);
44-
rcpp_result_gen = Rcpp::wrap(Adams_run(params, instate, forcing, fishing, stanzas, StartYear, EndYear, InitDeriv));
46+
Rcpp::traits::input_parameter< int >::type spnum(spnumSEXP);
47+
Rcpp::traits::input_parameter< int >::type spstanza(spstanzaSEXP);
48+
rcpp_result_gen = Rcpp::wrap(Adams_run(params, instate, forcing, fishing, stanzas, StartYear, EndYear, InitDeriv, spnum, spstanza));
4549
return rcpp_result_gen;
4650
END_RCPP
4751
}
@@ -93,8 +97,8 @@ END_RCPP
9397
}
9498

9599
static const R_CallMethodDef CallEntries[] = {
96-
{"_Rpath_rk4_run", (DL_FUNC) &_Rpath_rk4_run, 7},
97-
{"_Rpath_Adams_run", (DL_FUNC) &_Rpath_Adams_run, 8},
100+
{"_Rpath_rk4_run", (DL_FUNC) &_Rpath_rk4_run, 9},
101+
{"_Rpath_Adams_run", (DL_FUNC) &_Rpath_Adams_run, 10},
98102
{"_Rpath_deriv_vector", (DL_FUNC) &_Rpath_deriv_vector, 8},
99103
{"_Rpath_SplitSetPred", (DL_FUNC) &_Rpath_SplitSetPred, 2},
100104
{"_Rpath_SplitUpdate", (DL_FUNC) &_Rpath_SplitUpdate, 6},

‎src/ecosim.cpp‎

Lines changed: 37 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,10 @@
66
// Currently does not contain aged-structured species.
77
// [[Rcpp::export]]
88
List rk4_run (List params, List instate, List forcing, List fishing, List stanzas,
9-
int StartYear, int EndYear){
9+
int StartYear, int EndYear, int spnum, int spstanza){
1010

1111
int y, m, dd, t;
12-
// Input rates are in units of years or years^-1. Integration is wri tten so
12+
// Input rates are in units of years or years^-1. Integration is written so
1313
// that integration timesteps always line up with months, for data reasons.
1414
// STEPS_PER_YEAR should be 12 (for months), and STEPS_PER_MONTH sets the
1515
// rk4 integration timestep. So effective integration timestep with respect
@@ -221,7 +221,7 @@ List outdat = List::create(
221221
// Currently does not contain aged-structured species.
222222
// [[Rcpp::export]]
223223
List Adams_run (List params, List instate, List forcing, List fishing, List stanzas,
224-
int StartYear, int EndYear, List InitDeriv){
224+
int StartYear, int EndYear, List InitDeriv, int spnum, int spstanza){
225225

226226
int y, m, dd;
227227
//std::cout << " x1 ";
@@ -265,7 +265,11 @@ int y, m, dd;
265265
NumericMatrix out_Winf(EndYear * 12, Nsplit + 1);
266266
NumericMatrix out_Ninf(EndYear * 12, Nsplit + 1);
267267
NumericMatrix out_Wrec(EndYear * 12, Totstanzas + 1);
268-
NumericMatrix out_Nrec(EndYear * 12, Totstanzas + 1);
268+
NumericMatrix out_Nrec(EndYear * 12, Totstanzas + 1);
269+
270+
// Species diagnostic outputs
271+
NumericMatrix out_species_deriv(EndYear * 12, 10);
272+
269273
//std::cout << " x1b ";
270274
// Parameter need to track catch by Gear
271275
const NumericVector FishFrom = as<NumericVector>(params["FishFrom"]);
@@ -324,6 +328,19 @@ int y, m, dd;
324328
NumericVector biomeq = as<NumericVector>(dyt["biomeq"]);
325329
NumericVector FishingLoss = as<NumericVector>(dyt["FishingLoss"]);
326330
NumericVector Qlink = as<NumericVector>(dyt["Qlink"]);
331+
332+
// Added for per-species diagnostics
333+
NumericVector FoodLoss = as<NumericVector>(dyt["FoodLoss"]);
334+
//"FoodGain",
335+
NumericVector DetritalGain = as<NumericVector>(dyt["DetritalGain"]);
336+
NumericVector FishingGain = as<NumericVector>(dyt["FishingGain"]);
337+
NumericVector UnAssimLoss = as<NumericVector>(dyt["UnAssimLoss"]);
338+
NumericVector ActiveRespLoss = as<NumericVector>(dyt["ActiveRespLoss"]);
339+
NumericVector MzeroLoss = as<NumericVector>(dyt["MzeroLoss"]);
340+
//"FishingLoss",
341+
NumericVector DetritalLoss = as<NumericVector>(dyt["DetritalLoss"]);
342+
NumericVector MigrateLoss = as<NumericVector>(dyt["MigrateLoss"]);
343+
327344

328345
// Now Update the new State Biomass using Adams-Basforth
329346
NumericVector new_Biomass =
@@ -434,6 +451,20 @@ int y, m, dd;
434451
out_Wrec(dd,sind) = WageS(Age1(isp,ist), isp);
435452
}
436453
}
454+
455+
// Write diagnostic species outputs
456+
//FoodGain, DetritalGain, FishingGain, FoodLoss, UnAssimLoss, ActiveRespLoss,
457+
//MzeroLoss, FishingLoss, DetritalLoss, MigrateLoss
458+
out_species_deriv(dd, 0) = FoodGain[spnum];
459+
out_species_deriv(dd, 1) = DetritalGain[spnum];
460+
out_species_deriv(dd, 2) = FishingGain[spnum];
461+
out_species_deriv(dd, 3) = FoodLoss[spnum];
462+
out_species_deriv(dd, 4) = UnAssimLoss[spnum];
463+
out_species_deriv(dd, 5) = ActiveRespLoss[spnum];
464+
out_species_deriv(dd, 6) = MzeroLoss[spnum];
465+
out_species_deriv(dd, 7) = FishingLoss[spnum];
466+
out_species_deriv(dd, 8) = DetritalLoss[spnum];
467+
out_species_deriv(dd, 9) = MigrateLoss[spnum];
437468
//std::cout << "x4";
438469
} // End of main months loop
439470

@@ -460,6 +491,7 @@ int y, m, dd;
460491
_["annual_Biomass"]=annual_Biomass,
461492
_["annual_QB"]=annual_QB,
462493
_["annual_Qlink"]=annual_Qlink,
494+
_["out_species_deriv"]=out_species_deriv,
463495
_["out_SSB"]=out_SSB,
464496
_["out_eggs"]=out_eggs,
465497
_["out_Winf"]=out_Winf,
@@ -835,6 +867,7 @@ int sp, links, prey, pred, gr, egr, dest, isp, ist, ieco;
835867
_["MzeroLoss"]=MzeroLoss,
836868
_["FishingLoss"]=FishingLoss,
837869
_["DetritalLoss"]=DetritalLoss,
870+
_["MigrateLoss"]=MigrateLoss,
838871
_["FishingThru"]=FishingThru,
839872
//_["PredSuite"]=PredSuite,
840873
//_["HandleSuite"]=HandleSuite,

0 commit comments

Comments
 (0)