Add A+ cohort

A_PLUS_APPROACHES_COMPARISON.md

@@ -0,0 +1,82 @@
+# A+ (Plus Group) Implementation: Two Approaches Compared
+
+## Overview
+
+Both the `aplus` and `aplus3` branches implement the same fisheries "plus group" (A+)
+accumulator — the oldest age class where survivors and newly-graduated individuals pool
+together rather than continuing to age into new bins — but they take structurally different
+approaches to scheduling and synchronizing A+ within the existing MPI cohort-parallel driver
+(`main_cohort.cpp`). Both have been verified numerically equivalent: with A+ enabled they
+produce identical checksums to each other, and with `-no-aplus` both exactly reproduce the
+pre-A+ baseline, satisfying the conservation check `no-aplus checksum == (A+-on normal
+checksum) + (A+-on A+ checksum)` to the last printed digit, on the real `skj_Fat.xml` config
+(99/100 age groups, 480 months, 10 ranks).
+
+## `aplus`: graph-based scheduling
+
+A+ is folded directly into the same `SeapodymCohortDependencyAnalyzer` task graph that already
+governs every ordinary cohort. Each calendar step of A+ becomes its own one-step task, appended
+after the normal task IDs, with two explicit dependency edges: on A+'s own previous step, and on
+the specific cohort step that "graduates" into it that month. Newborn cohorts get a third kind
+of edge added automatically — a dependency on the matching A+ step, so a spawning cohort can
+never be dispatched before A+'s density for that month is available. All of this is dispatched
+by the same `TaskStepManager`/`TaskStepWorker` pair used for everything else; A+ tasks are just
+another branch inside the existing `taskFunction`. No dedicated MPI rank, no extra communicator
+split, no custom message protocol — correctness follows directly from the dependency graph's
+existing "don't dispatch a task until every dependency is in the `completed` set" logic, a
+mechanism that was already in use and already exercised long before A+ existed.
+
+## `aplus3`: dedicated worker with ping-pong handoff
+
+A+ is removed from the dependency graph entirely and instead runs as its own independent loop
+on a dedicated MPI rank (`runAPlusWorker`). Cohorts that reach the oldest normal age class
+("feeders") hand their graduating density directly to that rank via a blocking send, embedding
+the payload in the message itself; A+ acknowledges once it has actually processed that calendar
+step, and only after that ACK does the feeder tell the manager it's done — which is what
+prevents anything downstream from running ahead of A+. Because feeders can complete out of
+order, A+ carries a small reorder buffer that holds early arrivals and drains them strictly by
+calendar step. This is a real custom protocol with its own state machine, and finding it correct
+took real debugging effort — we found and fixed a genuine out-of-order feeding bug here that the
+graph-based design was never exposed to, because it doesn't have a hand-rolled synchronization
+protocol to get wrong in the first place.
+
+## Pros and cons
+
+| | `aplus` (graph) | `aplus3` (dedicated worker) |
+|---|---|---|
+| Minimum ranks | 2 | 3 |
+| Extra machinery | None — reuses existing graph/dispatch | Dedicated worker loop, message tags, reorder buffer, ACK gating |
+| Ordering guarantee | Declarative (graph edges) | Operational (blocking ACK + buffer) |
+| A+ compute isolation | Shares dispatch queue with all cohorts | Own rank, never waits on farm backlog |
+| Rank utilization | Every rank does real cohort work | A+'s rank is mostly idle (480 single-step advances vs. 50+ cohorts/rank elsewhere) |
+| Debugging surface | Small — one well-tested mechanism | Larger — a protocol we had to design and validate ourselves |
+| Verified correctness | Yes | Yes |
+
+`aplus`'s case is really about simplicity and efficiency: it adds no new synchronization
+primitive, needs no extra rank, and every rank stays busy with real cohort work. Its downside is
+that A+'s own advancement rides in the same queue as everything else, so if the farm is ever
+fully saturated, A+ could in principle wait behind unrelated backlog — though the graph
+guarantees it's never *wrong*, just potentially not the very next thing scheduled.
+
+`aplus3`'s case is about isolation: A+ gets guaranteed, uncontended throughput on its own rank,
+decoupled from farm scheduling. The cost is a dedicated rank that, based on the timing captured
+during testing, is doing meaningfully less work than any other rank (a single one-step-per-month
+accumulator vs. dozens of full multi-step cohort trajectories), plus a bespoke protocol that's
+more code to maintain and was the source of the one genuinely subtle bug found in this project
+(the feeder-ordering gap).
+
+## Recommendation
+
+Converge on `aplus`. For the workload as it exists today — a single accumulator doing one step
+of ordinary adult dynamics per calendar month — there's no evidence A+ is compute-bound enough
+to need its own rank, and the graph-based design gets the same correctness guarantee with less
+code, fewer ranks, and a synchronization mechanism that's already been proven out across the
+rest of the scheduler rather than a new one built and debugged specifically for this feature.
+
+`aplus3` isn't wrong, and it's a legitimate design if A+'s per-step cost ever grows to actually
+compete for farm time — a much heavier per-step computation, or multiple accumulator bins each
+needing real throughput — but that's not the situation today, and paying for a dedicated rank
+plus a bespoke protocol ahead of that need isn't buying anything. Worth keeping `aplus3` around
+as a documented alternative rather than deleting it, since the two are now proven equivalent and
+it remains a reasonable fallback if requirements change, but it shouldn't be the branch new work
+builds on by default.

example-configs/density_skj_fat/plot_workers_details.py

@@ -61,7 +61,7 @@ def parse_logs(pattern):
                         task_id = int(m.group(1))
                         task_ids.append(task_id)
                         worker_ids.append(worker_id)
-                        phases.append('coh init')
+                        phases.append('init')
                         t_starts.append(t_start)
                         t_ends.append(ts)
                         steps.append(-1)

example-configs/density_skj_fat/rungenoa.sl

@@ -20,4 +20,4 @@ module load gimkl/2020a
 rm log_*.txt
 
 # Run the MPI job
-srun ../../bin/seapodym_cohort -s skj_fat.xml >& results/n${SLURM_NTASKS}.txt
+srun ../../bin/seapodym_cohort_aplus -s skj_fat.xml >& results/n${SLURM_NTASKS}.txt

example-configs/density_skj_fat/skj_fat.xml

@@ -244,7 +244,7 @@
 </life_stage>
 
 <nb_cohort_life_stage>
-        <skj>1 2 47</skj>
+        <skj>1 2 97</skj>
 </nb_cohort_life_stage>
 
 <age_recruit skj="3"/>
@@ -255,23 +255,39 @@
 		  30 30 30 30 30 30 30 30 30 30 30 30 
 		  30 30 30 30 30 30 30 30 30 30 30 30 
 		  30 30 30 30 30 30 30 30 30 30 30 30 
-		  30 330</skj>
+		  30 30 30 30 30 30 30 30 30 30 30 30 
+		  30 30 30 30 30 30 30 30 30 30 30 30 		  
+		  30 30 30 30 30 30 30 30 30 30 30 30 
+		  30 30 30 30 30 30 30 30 30 30 30 30 
+		  30 30 30 2970
+	</skj>
 </sp_unit_cohort>
 
 <length>
 	<skj>3     4.5   7.52  11.65 16.91 21.83 26.43 30.72 34.73 38.49 41.99 45.27 
 		 48.33 51.19 53.86 56.36 58.70 60.88 62.92 64.83 66.61 68.27 69.83 71.28 
 		 72.64 73.91 75.10 76.21 77.25 78.22 79.12 79.97 80.76 81.50 82.19 82.83 
 		 83.50 83.50 83.50 83.50 83.50 83.50 83.50 83.50 83.50 83.50 83.50 83.50
-		 85.00 87.96 </skj>
+		 85.00 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96
+		 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96
+		 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96
+		 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96
+		 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96 87.96		 
+		 </skj>
 </length>
 
 <weight>
 	<skj>0.0003 0.00109 0.00571 0.023 0.08  0.18  0.32  0.53  0.78  1.09  1.44  1.84
 		 2.27   2.73    3.21    3.72  4.23  4.76  5.30  5.83  6.36  6.89  7.40  7.91 
 		 8.41   8.89    9.36    9.81  10.25 10.66 11.07 11.45 11.82 12.17 12.51 12.83 
 		 13.5   13.5    13.5    13.5  13.50 13.50 13.50 13.50 13.50 13.50 13.50 13.50
-		 14.0	15.56  </skj>
+		 14.0	15.56   15.56   15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56
+		 15.56  15.56   15.56   15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56
+		 15.56  15.56   15.56   15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56
+		 15.56  15.56   15.56   15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56
+		 15.56  15.56   15.56   15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56
+		 15.56  15.56   15.56   15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56 15.56		 		 
+		 </skj>
 </weight>
 
 

src/SeapodymCohort.cpp

@@ -3,6 +3,8 @@
 #include "Date.h"
 #include "sys/stat.h"
 #include <chrono>
+#include <cmath>
+#include <iostream>
 #include "DistDataCollector.h"
 #include "DataProvider.h"
 
@@ -30,15 +32,15 @@ std::vector<double> SeapodymCohort::GetCohortDensity()
 	return vec;
 }
 
-void SeapodymCohort::InitializeCohort(dvar_vector& x, DistDataCollector& dataCollector, const bool writeoutputfiles) 
+void SeapodymCohort::InitializeCohort(dvar_vector& x, DistDataCollector& dataCollector, int numTimeSteps, const bool writeoutputfiles)
 {
 
-    double t_all = MPI_Wtime();
-    double t_rs = 0.0;
-    double t_reset = MPI_Wtime();
-    //Reset model parameters: moved to worker init in main_cohort.cpp
-    //reset(x);
-    time_xreset += MPI_Wtime() - t_reset;	
+double t_all = MPI_Wtime();
+double t_rs = 0.0;
+double t_reset = MPI_Wtime();
+	//Reset model parameters:
+	reset(x);
+time_xreset += MPI_Wtime() - t_reset;	
 
 	//----------------------------------------------//
 	//	ALLOCATE AND INITIALIZE COHORT DENSITY	//
@@ -88,6 +90,47 @@ double t_c = MPI_Wtime();
 			}
 t_copy_acc += MPI_Wtime() - t_c;
 		}
+
+		if (aPlusEnabled){
+			// Plus group's density as of the end of t-1. SpawningBiomass_comp's
+			// own loop (a < param->sp_nb_cohorts[sp]) always included the
+			// oldest age class, A+ or not; nb_age_class here IS exactly that
+			// oldest (A+) index, since it excludes A+ from the "normal" range
+			// by construction. A+ publishes its density under this same
+			// chunk-id formula each step - see computeAPlusChunkId() and its
+			// use in main_cohort.cpp's taskFunction().
+			int chunk_id = nb_age_class * numTimeSteps + (tstart_cohort - 1);
+			dataCollector.getAsync(chunk_id, data.data());
+
+double t_f2 = MPI_Wtime();
+			dataCollector.flush();
+t_flush_acc += MPI_Wtime() - t_f2;
+
+			if (std::isnan(data[0])) {
+				// NaN (DistDataCollector::BAD_VALUE) here means chunk_id was
+				// never written before this cohort tried to read it - a
+				// missing/out-of-order dependency, not a numerical blow-up.
+				// Note: BAD_VALUE == NaN can never be detected with ==
+				// (NaN != NaN by definition), hence std::isnan() here.
+				std::cerr << "ERROR: cohort " << cohort_id
+					<< " (tstart_cohort=" << tstart_cohort << ") read BAD_VALUE (NaN) "
+					<< "from A+ chunk " << chunk_id << " (t=" << (tstart_cohort - 1) << ")\n";
+				MPI_Abort(MPI_COMM_WORLD, 1);
+			}
+
+double t_c2 = MPI_Wtime();
+			int index = 0;
+			for (int i = map.imin1; i <= map.imax1; i++){
+				const int jmin1 = map.jinf1[i];
+				const int jmax1 = map.jsup1[i];
+				for (int j = jmin1 ; j <= jmax1; j++){
+					mat.dvarDensity(0,nb_age_class).elem_value(i,j) = data[index];
+					index++;
+				}
+			}
+t_copy_acc += MPI_Wtime() - t_c2;
+		}
+
 		dataCollector.endEpoch(); // should be as late as possible
 
 double t_total = MPI_Wtime() - t_block;
@@ -119,26 +162,70 @@ double t_sp = MPI_Wtime();
 time_spawning += MPI_Wtime() - t_sp;
 	}
 
+	FinishInitialize(writeoutputfiles);
+
+time_init_cohort_restart += t_rs;
+time_init_cohort_spawning += MPI_Wtime() - t_all - t_rs;
+}
+
+void SeapodymCohort::InitializeAPlus(dvar_vector& x, const std::vector<double>& mergedDensity, bool seedFromFile)
+{
+	//Reset model parameters:
+	reset(x);
+
+	//----------------------------------------------//
+	//	ALLOCATE AND INITIALIZE COHORT DENSITY	//
+	//----------------------------------------------//
+	dvarCohortDensity.allocate(map.imin1, map.imax1, map.jinf1, map.jsup1);
+	if (seedFromFile){
+		// t=0: no upstream task dependency; the plus group's opening balance
+		// is just the initial condition for its age bin (age_start ==
+		// nb_age_class), read from file exactly like any other initial cohort.
+		RestoreDistributions(mat.nb_age_built);
+		dvarCohortDensity = mat.init_density_species(0,age_start);
+	} else {
+		// t>0: mergedDensity already holds (previous A+ pool) + (individuals
+		// that just graduated into the top age class this step), flattened
+		// in the same map.imin1..imax1/jinf1..jsup1 order GetCohortDensity()
+		// uses. This is the density stepForward() will apply this step's
+		// mortality/movement/feeding-habitat dynamics to - the AgePlus merge
+		// happens here, before dynamics, not instead of them.
+		int index = 0;
+		for (int i = map.imin1; i <= map.imax1; i++){
+			const int jmin1 = map.jinf1[i];
+			const int jmax1 = map.jsup1[i];
+			for (int j = jmin1 ; j <= jmax1; j++){
+				dvarCohortDensity.elem_value(i,j) = mergedDensity[index];
+				index++;
+			}
+		}
+	}
+
+	FinishInitialize(false);
+}
+
+void SeapodymCohort::FinishInitialize(bool writeoutputfiles)
+{
 	if (writeoutputfiles){
 		if (!param->gcalc())
 			ConsoleOutput(0,0);
 	}
 
 	//----------------------------------------------//
 	// PRECOMPUTE VARIABLES USED IN COHORT MODELING	//
-	//----------------------------------------------//	
+	//----------------------------------------------//
 	//precompute thermal habitat parameters
 	for (int sp=0; sp < nb_species; sp++){
 		func.Vars_at_age_precomp(*param,sp);
 		func.mortality_range_age_comp(*param,mat,sp);
-	
+
 		//precompute seasonal switch function
 		if (param->seasonal_migrations[sp]){
 			func.Seasonal_switch_year_precomp(*param,mat,map,
 						value(param->dvarsSpawning_season_peak[sp]),
 						value(param->dvarsSpawning_season_start[sp]),sp);
 		}
-	}	
+	}
 
 	getDate(jday, t_count);// In order to get qtr
 	if (param->type_oxy==1)
@@ -149,9 +236,6 @@ time_spawning += MPI_Wtime() - t_sp;
 	age = age_start;
 	past_month = month;
 	past_qtr = qtr;
-
-time_init_cohort_restart += t_rs;
-time_init_cohort_spawning += MPI_Wtime() - t_all - t_rs;
 }
 
 void SeapodymCohort::stepForward(bool writeoutputfiles)