diff --git a/doc/codes-vis-readme.md b/doc/codes-vis-readme.md
index b18d0ebc..07b723bf 100644
--- a/doc/codes-vis-readme.md
+++ b/doc/codes-vis-readme.md
@@ -6,6 +6,8 @@ on the ROSS webpage.
 
 There are currently 4 types of instrumentation: GVT-based, real time sampling, virtual time sampling, and event tracing.
 See the ROSS documentation for more info on the specific options or use `--help` with your model.
+The binary layout of the model-stats output files (and a reference Python parser) is documented in
+[model-stats-binary-format.md](model-stats-binary-format.md).
 To collect data about the simulation engine, no changes are needed to model code for any of the instrumentation modes.
 Some additions to the model code is needed in order to turn on any model-level data collection.
 See the "Model-level data sampling" section on [ROSS Instrumentation blog post](http://ross-org.github.io/instrumentation/instrumentation.html).
@@ -109,7 +111,7 @@ st_model_types svr_model_types[] = {
 
 static void svr_register_model_types()
 {
-    st_model_type_register("ns-lp", &svr_model_types[0]);
+    st_model_type_register("nw-lp", &svr_model_types[0]);
 }
 
 int main(int argc, char **argv)
@@ -135,12 +137,18 @@ modes are collecting model-level data as well.
 ### CODES LPs that currently have event type collection implemented:
 If you're using any of the following CODES models, you don't have to add anything, unless you want to change the data that's being collected.
 - nw-lp (model-net-mpi-replay.c)
+- ping pong tutorial server LP (doc/example/tutorial-synthetic-ping-pong.c) — also collects
+  model-level stats in all sampling modes; a good minimal example to copy from
 - dfly server LP (model-net-synthetic.c)
 - custom dfly server LP (model-net-synthetic-custom-dfly.c)
 - fat tree server LP (model-net-synthetic-fattree.c)
 - slimfly server LP (model-net-synthetic-slimfly.c)
 - original dragonfly router and terminal LPs (dragonfly.c)
 - dragonfly custom router and terminal LPs (dragonfly-custom.cxx)
+- dragonfly dally router and terminal LPs (dragonfly-dally.cxx) — richest model-stats
+  payload (packet rates, routing-path mix, stalls, occupancy snapshots); layout in
+  [model-stats-binary-format.md](model-stats-binary-format.md)
+- dragonfly plus router and terminal LPs (dragonfly-plus.cxx)
 - slimfly router and terminal LPs (slimfly.c)
 - fat tree switch and terminal LPs (fat-tree.c)
 - model-net-base-lp (model-net-lp.c)
diff --git a/doc/example/tutorial-ping-pong-surrogate.conf.in b/doc/example/tutorial-ping-pong-surrogate.conf.in
index 8040d921..424a1fed 100644
--- a/doc/example/tutorial-ping-pong-surrogate.conf.in
+++ b/doc/example/tutorial-ping-pong-surrogate.conf.in
@@ -41,7 +41,7 @@ PARAMS
 # bandwidth in GiB/s for compute node-router channels
    cn_bandwidth="2.0";
 # ROSS message size
-   message_size="440";
+   message_size="456";
 # number of compute nodes connected to router, dictated by dragonfly config
 # file
    num_cns_per_router="2";
diff --git a/doc/example/tutorial-ping-pong-surrogate.yaml.in b/doc/example/tutorial-ping-pong-surrogate.yaml.in
index 1f621040..9bbca239 100644
--- a/doc/example/tutorial-ping-pong-surrogate.yaml.in
+++ b/doc/example/tutorial-ping-pong-surrogate.yaml.in
@@ -50,7 +50,7 @@ topology:
 # scheduler options (alternative: round-robin)
     modelnet_scheduler: fcfs
 # ROSS message size
-    message_size: 440
+    message_size: 456
 # folder path to store packet latency from terminal to terminal, if no value is given it won't save anything
     save_packet_latency_path: "${PACKET_LATENCY_TRACE_PATH}"
 # folder path to store router-local timing rows for router queueing-delay training
diff --git a/doc/example/tutorial-ping-pong.conf.in b/doc/example/tutorial-ping-pong.conf.in
index 042ee417..3cf071cf 100644
--- a/doc/example/tutorial-ping-pong.conf.in
+++ b/doc/example/tutorial-ping-pong.conf.in
@@ -1,3 +1,7 @@
+# To also collect model-level stats through the ROSS instrumentation layer, add e.g.
+#   --model-stats=1 --num-gvt=8 --stats-path=model-stats-out
+# to the run command (see doc/model-stats-binary-format.md for all modes and the
+# output format; scripts/parse-ross-model-stats.py parses the result).
 LPGROUPS
 {
    MODELNET_GRP
@@ -38,7 +42,7 @@ PARAMS
 # bandwidth in GiB/s for compute node-router channels
    cn_bandwidth="2.0";
 # ROSS message size
-   message_size="440";
+   message_size="456";
 # number of compute nodes connected to router, dictated by dragonfly config
 # file
    num_cns_per_router="2";
diff --git a/doc/example/tutorial-ping-pong.yaml.in b/doc/example/tutorial-ping-pong.yaml.in
index ad0dddf9..819f76f5 100644
--- a/doc/example/tutorial-ping-pong.yaml.in
+++ b/doc/example/tutorial-ping-pong.yaml.in
@@ -9,6 +9,11 @@ schema_version: 1
 # as in the .conf. The compiler derives repetitions = num_groups * num_routers and
 # lays out the [workload, terminal, router] LPs, so nw-lp / modelnet_dragonfly_dally*
 # are emitted for you.
+#
+# To also collect model-level stats through the ROSS instrumentation layer, add e.g.
+#   --model-stats=1 --num-gvt=8 --stats-path=model-stats-out
+# to the run command (see doc/model-stats-binary-format.md for all modes and the
+# output format; scripts/parse-ross-model-stats.py parses the result).
 
 components:
   compute_host:
@@ -49,7 +54,7 @@ topology:
 # scheduler options (alternative: round-robin)
     modelnet_scheduler: fcfs
 # ROSS message size
-    message_size: 440
+    message_size: 456
 # router buffer occupancy snapshots
     router_buffer_snapshots: [ ${BUFFER_SNAPSHOTS} ]
 # folder path to store packet latency from terminal to terminal, if no value is given it won't save anything
diff --git a/doc/example/tutorial-synthetic-ping-pong.c b/doc/example/tutorial-synthetic-ping-pong.c
index ee85d8c6..b24dbb58 100644
--- a/doc/example/tutorial-synthetic-ping-pong.c
+++ b/doc/example/tutorial-synthetic-ping-pong.c
@@ -4,6 +4,9 @@
  * See LICENSE notice in top-level directory
  */
 
+#include <stdint.h>
+#include <string.h>
+
 #include "codes/model-net.h"
 #include "codes/codes_mapping.h"
 #include "codes/surrogate/init.h" // just needed for stats on surrogate-mode
@@ -39,6 +42,30 @@ struct svr_msg {
     model_net_event_return event_rc; //helper to encode data relating to CODES rng usage
     // Used for rollback
     tw_stime previous_ts;
+    // Virtual time at which the PING carrying this message was created. It is echoed back
+    // unchanged in the PONG so the original sender can compute the round trip time (RTT).
+    // NOTE: PING/PONG messages are built as stack structs below, so this field must be
+    // explicitly assigned everywhere a message is created or it carries garbage.
+    tw_stime ping_send_ts;
+    // RTT computed when a PONG is received; stored in the message so the reverse event
+    // handler can subtract the exact same value on rollback (see handle_pong_rev_event)
+    tw_stime saved_rtt;
+};
+
+/* Sample structure filled in by svr_sample_fn() for the virtual-time (analysis LP) sampling
+ * mode of the ROSS instrumentation layer (--model-stats=3 or 4). All fields are 8 bytes wide
+ * so the struct has no padding holes -- the analysis LPs write these structs to
+ * ross-stats-analysis-lps.bin verbatim, so a well-defined layout keeps the file parseable. */
+struct svr_sample {
+    tw_lpid svr_id;     /* which server this sample belongs to */
+    int64_t pings_sent; /* all following counters are deltas since the previous sample */
+    int64_t pings_recvd;
+    int64_t pongs_sent;
+    int64_t pongs_recvd;
+    int64_t bytes_sent;
+    double rtt_sum;    /* sum (ns) of PING->PONG round trips completed in this interval */
+    int64_t rtt_count; /* mean RTT for the interval = rtt_sum / rtt_count */
+    tw_stime end_time; /* virtual time at which this sample was taken */
 };
 
 struct svr_state {
@@ -50,6 +77,29 @@ struct svr_state {
     tw_stime start_ts;        /* time that this LP started sending requests */
     tw_stime end_ts;          /* time that this LP ended sending requests */
     int payload_sum;          /* the running sum of all payloads received */
+
+    /* Shadow accumulators for the ROSS instrumentation layer (model-level stats).
+     *
+     * These record per-sampling-interval deltas: forward event handlers increment them,
+     * reverse event handlers decrement them (so they stay correct across rollbacks in
+     * optimistic mode), and the collection callbacks zero them after copying them out.
+     * Zeroing inside the GVT-mode callback is safe because GVT-time collection only ever
+     * runs on committed simulation time -- no event before the GVT can be rolled back.
+     *
+     * There are TWO independent accumulator sets, following the convention used by the
+     * dragonfly models: the scalars below feed the GVT/real-time sampling callback
+     * (svr_model_stat_collect, --model-stats=1/2) while vt_sample feeds the virtual-time
+     * analysis-LP callback (svr_sample_fn, --model-stats=3). Both paths zero their
+     * accumulators at each sample, so sharing one set would make the two modes steal
+     * deltas from each other when running with --model-stats=4 (all modes at once). */
+    uint32_t pings_sent_sample;
+    uint32_t pings_recvd_sample;
+    uint32_t pongs_sent_sample;
+    uint32_t pongs_recvd_sample;
+    int64_t bytes_sent_sample; /* payload bytes handed to model_net_event() */
+    double rtt_sum_sample;     /* sum (ns) of PING->PONG round trips completed */
+    uint32_t rtt_count_sample;
+    struct svr_sample vt_sample; /* accumulator twin for the virtual-time sampling mode */
 };
 
 /* declaration of functions */
@@ -66,6 +116,125 @@ tw_lptype svr_lp = {
     (commit_f)NULL,   (final_f)svr_finalize, (map_f)codes_mapping, sizeof(svr_state),
 };
 
+/* ---- ROSS instrumentation layer: model-level stats + event tracing -----------------------
+ *
+ * The callbacks below let the ROSS instrumentation layer collect model-defined data
+ * alongside its simulation-engine data. Nothing here runs unless instrumentation is
+ * requested on the command line, e.g.:
+ *
+ *   --model-stats=1              collect model data at GVT (every --num-gvt GVTs, default 10)
+ *   --model-stats=2              collect at real-time intervals (--rt-interval, ms)
+ *   --model-stats=3              collect at virtual-time intervals via analysis LPs
+ *                                (--vt-interval / --vt-samp-end). NOTE: pass an explicit
+ *                                --vt-samp-end for this model! It defaults to g_tw_ts_end,
+ *                                which main() sets to 24h of virtual time, so the analysis
+ *                                LPs would keep scheduling sample events long after the
+ *                                ping pong traffic ends.
+ *   --model-stats=4              all of the above
+ *   --event-trace=1              trace events (svr_event_collect below)
+ *   --stats-path=<dir>           output directory (default: stats-output)
+ *
+ * GVT/RT-mode records go to <stats-path>/ross-stats-model.bin, virtual-time samples to
+ * <stats-path>/ross-stats-analysis-lps.bin. See doc/model-stats-binary-format.md for the
+ * exact record layouts and a reference parser (scripts/parse-ross-model-stats.py). */
+
+/* Event-type collection for ROSS event tracing (--event-trace) */
+static void svr_event_collect(svr_msg* m, tw_lp* lp, char* buffer, int* collect_flag) {
+    (void)lp;
+    (void)collect_flag;
+    int type = (int)m->svr_event_type;
+    memcpy(buffer, &type, sizeof(type));
+}
+
+/* GVT- and real-time-mode collection callback (--model-stats=1 or 2). ROSS hands us a
+ * buffer sized to the mstat_sz we declare in svr_model_types below; we memcpy our shadow
+ * accumulators into it field by field (byte-serial, so the output has no padding holes)
+ * and zero each accumulator to start the next interval. */
+static void svr_model_stat_collect(svr_state* s, tw_lp* lp, char* buffer) {
+    (void)lp;
+    int index = 0;
+
+    memcpy(&buffer[index], &s->svr_id, sizeof(s->svr_id));
+    index += sizeof(s->svr_id);
+
+    memcpy(&buffer[index], &s->pings_sent_sample, sizeof(s->pings_sent_sample));
+    index += sizeof(s->pings_sent_sample);
+    s->pings_sent_sample = 0;
+
+    memcpy(&buffer[index], &s->pings_recvd_sample, sizeof(s->pings_recvd_sample));
+    index += sizeof(s->pings_recvd_sample);
+    s->pings_recvd_sample = 0;
+
+    memcpy(&buffer[index], &s->pongs_sent_sample, sizeof(s->pongs_sent_sample));
+    index += sizeof(s->pongs_sent_sample);
+    s->pongs_sent_sample = 0;
+
+    memcpy(&buffer[index], &s->pongs_recvd_sample, sizeof(s->pongs_recvd_sample));
+    index += sizeof(s->pongs_recvd_sample);
+    s->pongs_recvd_sample = 0;
+
+    memcpy(&buffer[index], &s->bytes_sent_sample, sizeof(s->bytes_sent_sample));
+    index += sizeof(s->bytes_sent_sample);
+    s->bytes_sent_sample = 0;
+
+    memcpy(&buffer[index], &s->rtt_sum_sample, sizeof(s->rtt_sum_sample));
+    index += sizeof(s->rtt_sum_sample);
+    s->rtt_sum_sample = 0;
+
+    memcpy(&buffer[index], &s->rtt_count_sample, sizeof(s->rtt_count_sample));
+    index += sizeof(s->rtt_count_sample);
+    s->rtt_count_sample = 0;
+}
+
+/* Size of the payload svr_model_stat_collect() writes: keep in sync with the fields above */
+#define SVR_MSTAT_SZ                                                                               \
+    (sizeof(tw_lpid) + 4 * sizeof(uint32_t) + sizeof(int64_t) + sizeof(double) + sizeof(uint32_t))
+
+/* Virtual-time-mode sampling callback (--model-stats=3): the analysis LP for this KP calls
+ * this at each --vt-interval with a sample buffer of sample_struct_sz bytes. We copy the
+ * vt_sample accumulator out and zero it. Unlike GVT collection this CAN be rolled back,
+ * which is why there is a reverse callback restoring the accumulator from the sample. */
+static void svr_sample_fn(svr_state* s, tw_bf* b, tw_lp* lp, struct svr_sample* sample) {
+    (void)b;
+    *sample = s->vt_sample;
+    sample->svr_id = s->svr_id;
+    sample->end_time = tw_now(lp);
+
+    memset(&s->vt_sample, 0, sizeof(s->vt_sample));
+}
+
+static void svr_sample_rc_fn(svr_state* s, tw_bf* b, tw_lp* lp, struct svr_sample* sample) {
+    (void)b;
+    (void)lp;
+    /* exact reversal: the accumulator was zeroed after the sample was taken, so adding the
+     * sampled values back restores it (any post-sample activity is already re-included by
+     * the reverse handlers of the rolled-back events themselves) */
+    s->vt_sample.pings_sent += sample->pings_sent;
+    s->vt_sample.pings_recvd += sample->pings_recvd;
+    s->vt_sample.pongs_sent += sample->pongs_sent;
+    s->vt_sample.pongs_recvd += sample->pongs_recvd;
+    s->vt_sample.bytes_sent += sample->bytes_sent;
+    s->vt_sample.rtt_sum += sample->rtt_sum;
+    s->vt_sample.rtt_count += sample->rtt_count;
+}
+
+st_model_types svr_model_types[] = {{(ev_trace_f)svr_event_collect, sizeof(int),
+                                     (model_stat_f)svr_model_stat_collect, SVR_MSTAT_SZ,
+                                     (sample_event_f)svr_sample_fn,
+                                     (sample_revent_f)svr_sample_rc_fn, sizeof(struct svr_sample)},
+                                    {NULL, 0, NULL, 0, NULL, NULL, 0}};
+
+static const st_model_types* svr_get_model_stat_types(void) {
+    return (&svr_model_types[0]);
+}
+
+/* Attach our instrumentation callbacks to the "nw-lp" LP type; codes_mapping_setup() then
+ * hands them to ROSS (st_model_settype) for every LP registered under that name */
+static void svr_register_model_types(void) {
+    st_model_type_register("nw-lp", svr_get_model_stat_types());
+}
+/* ---- end of ROSS instrumentation additions ---------------------------------------------- */
+
 const tw_optdef app_opt[] = {
     TWOPT_GROUP("Model net synthetic traffic "),
     TWOPT_UINT("num_messages", num_msgs, "Number of PING messages to be generated per terminal "),
@@ -98,6 +267,16 @@ static void svr_init(svr_state* s, tw_lp* lp) {
                                                  0); /* turns the LP Global ID into the server ID */
     s->payload_sum = 0;
 
+    /* zero the instrumentation shadow accumulators (both sampling paths) */
+    s->pings_sent_sample = 0;
+    s->pings_recvd_sample = 0;
+    s->pongs_sent_sample = 0;
+    s->pongs_recvd_sample = 0;
+    s->bytes_sent_sample = 0;
+    s->rtt_sum_sample = 0.0;
+    s->rtt_count_sample = 0;
+    memset(&s->vt_sample, 0, sizeof(s->vt_sample));
+
     //Now we create and send a self KICKOFF message - this is a PDES coordination event and thus doesn't need to be injected into the connected network
     //so we won't use model_net_event(), that's reserved for stuff we want to send across the network
 
@@ -143,6 +322,7 @@ static void handle_kickoff_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp* lp)
     ping_msg.svr_event_type = PING; //set it to type PING
     ping_msg.payload_value =
         tw_rand_integer(lp->rng, 1, 10); //encode a random payload value to it from [1,10]
+    ping_msg.ping_send_ts = tw_now(lp);  //stamp the send time so the PONG can carry it back
 
     codes_mapping_get_lp_info(
         lp->gid, group_name, &group_index, lp_type_name, &lp_type_index, NULL, &rep_id,
@@ -150,6 +330,12 @@ static void handle_kickoff_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp* lp)
     global_dest =
         codes_mapping_get_lpid_from_relative(local_dest, group_name, lp_type_name, NULL, 0);
     s->ping_msg_sent_count++;
+    /* instrumentation deltas: every increment here has a matching decrement in the reverse
+     * handler below, keeping the accumulators correct across optimistic-mode rollbacks */
+    s->pings_sent_sample++;
+    s->bytes_sent_sample += PAYLOAD_SZ;
+    s->vt_sample.pings_sent++;
+    s->vt_sample.bytes_sent += PAYLOAD_SZ;
     m->event_rc = model_net_event(net_id, "test", global_dest, PAYLOAD_SZ, 0.0, sizeof(svr_msg),
                                   (const void*)&ping_msg, 0, NULL, lp);
 }
@@ -159,6 +345,10 @@ static void handle_kickoff_rev_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp*
     model_net_event_rc2(lp,
                         &m->event_rc); //undo any model_net_event calls encoded into this message
     s->ping_msg_sent_count--; //undo the increment of the ping_msg_sent_count in the server state
+    s->pings_sent_sample--;   //undo the instrumentation deltas
+    s->bytes_sent_sample -= PAYLOAD_SZ;
+    s->vt_sample.pings_sent--;
+    s->vt_sample.bytes_sent -= PAYLOAD_SZ;
     tw_rand_reverse_unif(lp->rng); //reverse the rng call for creating a payload value;
     tw_rand_reverse_unif(lp->rng); //reverse the rng call for getting a local_dest
 }
@@ -166,6 +356,8 @@ static void handle_kickoff_rev_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp*
 static void handle_ping_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp* lp) {
     (void)b;
     s->ping_msg_recvd_count++; //increment the counter for ping messages received
+    s->pings_recvd_sample++;   //instrumentation deltas (reversed in handle_ping_rev_event)
+    s->vt_sample.pings_recvd++;
 
     int original_sender = m->sender_id; //this is the server we need to send a PONG message back to
     s->payload_sum += m->payload_value; //increment our running sum of payload values received
@@ -173,6 +365,7 @@ static void handle_ping_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp* lp) {
     svr_msg pong_msg;
     pong_msg.sender_id = s->svr_id;
     pong_msg.svr_event_type = PONG;
+    pong_msg.ping_send_ts = m->ping_send_ts; //echo the PING's send time back to the sender
     // only ping messages contain a payload value - not every value in a message struct must be utilized by all messages!
 
     codes_mapping_get_lp_info(
@@ -181,6 +374,10 @@ static void handle_ping_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp* lp) {
     tw_lpid global_dest =
         codes_mapping_get_lpid_from_relative(original_sender, group_name, lp_type_name, NULL, 0);
     s->pong_msg_sent_count++;
+    s->pongs_sent_sample++;
+    s->bytes_sent_sample += PAYLOAD_SZ;
+    s->vt_sample.pongs_sent++;
+    s->vt_sample.bytes_sent += PAYLOAD_SZ;
     m->event_rc = model_net_event(net_id, "test", global_dest, PAYLOAD_SZ, 0.0, sizeof(svr_msg),
                                   (const void*)&pong_msg, 0, NULL, lp);
 }
@@ -192,10 +389,29 @@ static void handle_ping_rev_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp* lp)
     s->pong_msg_sent_count--;
     s->payload_sum -= m->payload_value; //undo the increment of the payload sum
     s->ping_msg_recvd_count--; //undo the increment of the counter for ping messages received
+    s->pings_recvd_sample--;   //undo the instrumentation deltas
+    s->pongs_sent_sample--;
+    s->bytes_sent_sample -= PAYLOAD_SZ;
+    s->vt_sample.pings_recvd--;
+    s->vt_sample.pongs_sent--;
+    s->vt_sample.bytes_sent -= PAYLOAD_SZ;
 }
 
 static void handle_pong_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp* lp) {
     s->pong_msg_recvd_count++; //increment the counter for ping messages received
+    s->pongs_recvd_sample++;
+    s->vt_sample.pongs_recvd++;
+
+    /* this PONG closes the loop on a PING we sent earlier: the PING's send time made the
+     * round trip with the message, so the RTT is simply now minus then. The computed value
+     * is stored in the message so the reverse handler can subtract exactly what was added.
+     * This happens before the early return below -- every PONG carries an RTT, even the
+     * last one. */
+    m->saved_rtt = tw_now(lp) - m->ping_send_ts;
+    s->rtt_sum_sample += m->saved_rtt;
+    s->rtt_count_sample++;
+    s->vt_sample.rtt_sum += m->saved_rtt;
+    s->vt_sample.rtt_count++;
 
     if (s->ping_msg_sent_count >=
         num_msgs) //if we've sent enough ping messages, then we stop and don't send any more
@@ -213,6 +429,7 @@ static void handle_pong_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp* lp) {
     ping_msg.sender_id = s->svr_id; //encode our server ID into the new ping message
     ping_msg.svr_event_type = PING; //set it to type PING
     ping_msg.payload_value = tw_rand_integer(lp->rng, 1, 10); //encode a random payload value to it
+    ping_msg.ping_send_ts = tw_now(lp); //stamp the send time so the PONG can carry it back
 
     codes_mapping_get_lp_info(
         lp->gid, group_name, &group_index, lp_type_name, &lp_type_index, NULL, &rep_id,
@@ -220,6 +437,10 @@ static void handle_pong_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp* lp) {
     tw_lpid global_dest =
         codes_mapping_get_lpid_from_relative(send_to, group_name, lp_type_name, NULL, 0);
     s->ping_msg_sent_count++;
+    s->pings_sent_sample++;
+    s->bytes_sent_sample += PAYLOAD_SZ;
+    s->vt_sample.pings_sent++;
+    s->vt_sample.bytes_sent += PAYLOAD_SZ;
     m->event_rc = model_net_event(net_id, "test", global_dest, PAYLOAD_SZ, 0.0, sizeof(svr_msg),
                                   (const void*)&ping_msg, 0, NULL, lp);
 }
@@ -229,12 +450,25 @@ static void handle_pong_rev_event(svr_state* s, tw_bf* b, svr_msg* m, tw_lp* lp)
         model_net_event_rc2(
             lp, &m->event_rc); //undo any model_net_event calls encoded into this message
         s->ping_msg_sent_count--;
+        s->pings_sent_sample--;
+        s->bytes_sent_sample -= PAYLOAD_SZ;
+        s->vt_sample.pings_sent--;
+        s->vt_sample.bytes_sent -= PAYLOAD_SZ;
         tw_rand_reverse_unif(lp->rng); //undo the rng for the new payload value
         tw_rand_reverse_unif(lp->rng); //undo the rng for the new server to send a ping to
         b->c1 = 0;
     }
 
     s->pong_msg_recvd_count--; //undo the increment of the counter for ping messages received
+    s->pongs_recvd_sample--;
+    s->vt_sample.pongs_recvd--;
+
+    /* undo the RTT recorded by the forward handler (unconditional -- it runs before the
+     * early-return branch above) using the exact value saved in the message */
+    s->rtt_sum_sample -= m->saved_rtt;
+    s->rtt_count_sample--;
+    s->vt_sample.rtt_sum -= m->saved_rtt;
+    s->vt_sample.rtt_count--;
 }
 
 static void svr_finalize(svr_state* s, tw_lp* lp) {
@@ -325,6 +559,12 @@ int main(int argc, char** argv) {
     model_net_register();
     svr_add_lp_type();
 
+    /* attach the ROSS instrumentation callbacks when any instrumentation mode was requested
+     * on the command line (--model-stats / --event-trace); must happen before
+     * codes_mapping_setup(), which is where the callbacks get handed to ROSS per LP */
+    if (g_st_ev_trace || g_st_model_stats || g_st_use_analysis_lps)
+        svr_register_model_types();
+
     codes_mapping_setup();
 
     net_ids = model_net_configure(&num_nets);
diff --git a/doc/model-stats-binary-format.md b/doc/model-stats-binary-format.md
new file mode 100644
index 00000000..2c06030c
--- /dev/null
+++ b/doc/model-stats-binary-format.md
@@ -0,0 +1,237 @@
+# Binary format of the ROSS model-level stats output
+
+This document specifies the on-disk layout of the model-level data that CODES models emit
+through the ROSS instrumentation layer, so external tools (e.g. NetMaestro) can parse it
+without reading the C source. It covers the payloads written by the ping pong tutorial's
+server LP (`doc/example/tutorial-synthetic-ping-pong.c`) and the dragonfly-dally terminal
+and router LPs (`src/networks/model-net/dragonfly-dally.cxx`).
+
+For how a model *registers* these callbacks, see [codes-vis-readme.md](codes-vis-readme.md).
+A reference parser lives at `scripts/parse-ross-model-stats.py` (stdlib-only Python).
+
+## Producing the files
+
+Model-level stats are enabled purely from the ROSS command line — no config-file changes:
+
+| flag | meaning |
+|---|---|
+| `--model-stats=1` | sample at GVT, every `--num-gvt`-th GVT (default 10) |
+| `--model-stats=2` | sample at wall-clock intervals of `--rt-interval` ms |
+| `--model-stats=3` | sample at virtual-time intervals of `--vt-interval` (until `--vt-samp-end`) via analysis LPs. `--vt-samp-end` defaults to the simulation end time — set it explicitly for models with a far-future `g_tw_ts_end` (the ping pong tutorial uses 24 h of virtual time), or the analysis LPs keep scheduling empty sample events long after the traffic ends |
+| `--model-stats=4` | all of the above |
+| `--event-trace=N` | event tracing (separate `evtrace` file; not covered here, see the ROSS docs) |
+| `--stats-path=<dir>` | output directory, default `stats-output` |
+| `--stats-prefix=<name>` | filename prefix, default `ross-stats` |
+| `--buffer-size=<bytes>` | per-collection-type circular buffer size |
+
+Example (ping pong tutorial, from the build directory):
+
+```
+mpirun -np 3 doc/example/tutorial-synthetic-ping-pong --sync=3 --num_messages=100 \
+    --model-stats=1 --num-gvt=8 --stats-path=model-stats-out -- doc/example/tutorial-ping-pong.conf
+```
+
+Outputs:
+
+- GVT and real-time samples → `<stats-path>/<prefix>-model.bin`
+- virtual-time samples → `<stats-path>/<prefix>-analysis-lps.bin`
+
+If `<stats-path>` already exists, ROSS appends `-<pid>-<timestamp>` to the directory name
+instead of reusing it — scripts should run in a fresh directory or glob for the suffix.
+
+## General properties (both files)
+
+- Everything is written in the **native byte order and type sizes of the machine that ran
+  the simulation** — in practice little-endian LP64 (`long` = 8 bytes). There is no
+  endianness marker in the file.
+- Records are packed back to back with **no alignment padding or framing between records**;
+  read until EOF. All record and payload layouts below have no internal padding holes.
+- In parallel runs each PE flushes its own circular buffer at GVT into one shared file via
+  MPI-IO, so records are **not globally sorted by time**; sort on the timestamp field after
+  parsing.
+- If a PE's buffer fills up, further records are **silently dropped** (a warning with the
+  missed byte count is printed at the end of the run); increase `--buffer-size` if that
+  happens.
+- Delta-style fields cover the interval since the previous sample of the same mode. The
+  stretch between the last sample and the end of the run is **not emitted**, so sums over
+  samples slightly undercount end-of-run totals.
+- Under optimistic sync (`--sync=3`), real-time-mode samples are taken mid-window, so a
+  later rollback can make individual RT deltas **slightly negative** (sums across samples
+  stay correct). GVT and VT samples only cover committed/rolled-forward time and don't
+  have this artifact; prefer GVT mode for clean per-interval data.
+- With `--model-stats=4`, the GVT and RT modes share one accumulator (ROSS invokes the
+  same collection callback for both), so each delta lands in whichever sample fires first
+  — GVT records and RT records each hold a partial view and only their **combined** sums
+  are complete. The virtual-time path keeps its own accumulator, so VT records are
+  complete on their own. Don't sum GVT+VT (that double-counts); pick one mode's records
+  per analysis.
+
+## `ross-stats-model.bin` — GVT and real-time samples
+
+Each record is: `sample_metadata` (24 B) + `model_metadata` (24 B) + `model_sz` payload
+bytes. The structs come from ROSS (`core/instrumentation/st-instrumentation.h`).
+
+`sample_metadata` — Python struct format `<iidd`:
+
+| off | type | field | notes |
+|---|---|---|---|
+| 0 | i32 | flag | always 3 (`MODEL_TYPE`) — validate this |
+| 4 | i32 | sample_sz | always 24 (= sizeof model_metadata) |
+| 8 | f64 | ts | virtual time of the sampled LP — use this for time axes |
+| 16 | f64 | real_time | wall-clock seconds since the run started |
+
+`model_metadata` — Python struct format `<IIIfiI`:
+
+| off | type | field | notes |
+|---|---|---|---|
+| 0 | u32 | peid | PE (MPI rank) |
+| 4 | u32 | kpid | KP |
+| 8 | u32 | lpid | **global LP id** (truncated to 32 bits) |
+| 12 | f32 | gvt | GVT at collection; low precision — prefer `ts` |
+| 16 | i32 | stats_type | 1 = GVT sample, 2 = real-time sample |
+| 20 | u32 | model_sz | payload size in bytes that follows |
+
+Dispatch the payload decoder on `model_sz` (all three payload sizes are distinct for any
+one run configuration) and cross-check with the id in the first payload field.
+
+### Payload: ping pong server LP (`nw-lp`) — 44 bytes, `<Q4IqdI`
+
+All counters are deltas over the sampling interval.
+
+| off | type | field |
+|---|---|---|
+| 0 | u64 | svr_id (server-relative id) |
+| 8 | u32 | pings_sent |
+| 12 | u32 | pings_recvd |
+| 16 | u32 | pongs_sent |
+| 20 | u32 | pongs_recvd |
+| 24 | i64 | bytes_sent (payload bytes handed to model-net) |
+| 32 | f64 | rtt_sum (ns; sum of PING→PONG round trips completed) |
+| 40 | u32 | rtt_count (mean RTT for the interval = rtt_sum / rtt_count) |
+
+### Payload: dragonfly-dally terminal — `72 + 16R + 8RQ` bytes
+
+R = `num_rails`, Q = `num_qos_levels` from the network config (tutorial config: R=1, Q=1
+→ **96 bytes**). Fields are deltas unless marked *snapshot* (instantaneous value at the
+sampling instant).
+
+| off | type | field |
+|---|---|---|
+| 0 | u64 | terminal_id (terminal-relative id) |
+| 8 | i64 | fin_chunks (chunks delivered to this terminal) |
+| 16 | i64 | data_size (payload bytes delivered) |
+| 24 | i64 | fin_hops (total hops over delivered chunks) |
+| 32 | f64 | fin_chunks_time (ns; total end-to-end latency of delivered chunks — mean chunk latency = fin_chunks_time / fin_chunks) |
+| 40 | f64 ×R | busy_time (ns the injection link spent busy, per rail) |
+| 40+8R | i64 | packets_gen (packets injected) |
+| 48+8R | i64 | packets_fin (packets fully delivered) |
+| 56+8R | i64 | min_fin (delivered chunks that took a minimal route) |
+| 64+8R | i64 | nonmin_fin (delivered chunks that took a non-minimal route) |
+| 72+8R | u64 ×R | stalled_chunks (chunks that could not inject due to full VC, per rail) |
+| 72+16R | i32 ×R·Q | vc_occupancy — *snapshot*, bytes occupying the injection VCs, rail-major (`[rail][qos]`) |
+| 72+16R+4RQ | i32 ×R·Q | terminal_length — *snapshot*, bytes queued at the NIC, rail-major |
+
+Note: while a terminal is frozen by the network surrogate (surrogate configs only), the
+pointer-backed fields (busy_time, stalled_chunks, vc_occupancy, terminal_length) read as
+zeros.
+
+### Payload: dragonfly-dally router — `8 + 32·radix` bytes
+
+radix = ports per router (tutorial config: 7 → **232 bytes**). Port order is fixed by the
+dragonfly connection manager: intra-group ports, then global ports, then terminal ports
+(see `codes/network-manager/dragonfly-network-manager.h`); mapping port index → link
+requires the topology/connection files used by the run.
+
+| off | type | field |
+|---|---|---|
+| 0 | u64 | router_id (router-relative id) |
+| 8 | (f64,i64) ×radix | interleaved per port: busy_time (ns), link_traffic (bytes forwarded) |
+| 8+16·radix | u64 ×radix | stalled_chunks (chunks diverted to the overflow queue, per port) |
+| 8+24·radix | i32 ×radix | vc_occupancy_sum — *snapshot*, bytes in all VCs of the port |
+| 8+28·radix | i32 ×radix | queued_count — *snapshot*, bytes in the port's overflow queue |
+
+## `ross-stats-analysis-lps.bin` — virtual-time samples
+
+Each record is: `lp_metadata` (48 B) + `sample_sz` payload bytes. Only LPs with a non-zero
+VT sample size appear (in a ping pong run: the server, terminal and router LPs).
+
+`lp_metadata` (`core/instrumentation/ross-lps/analysis-lp.h`) — Python struct format
+`<QQQddii` (note the id order differs from `model_metadata`):
+
+| off | type | field |
+|---|---|---|
+| 0 | u64 | lpid (global LP id, full width) |
+| 8 | u64 | kpid |
+| 16 | u64 | peid |
+| 24 | f64 | ts (virtual time of the sampling point) |
+| 32 | f64 | real_time |
+| 40 | i32 | sample_sz (payload bytes that follow) |
+| 44 | i32 | flag (always 3 = model) |
+
+The payloads are the C sample structs written verbatim, so they **contain pointer members
+whose on-disk bytes are meaningless** — parsers must skip those 8-byte slots (marked below).
+The variable-length arrays the pointers refer to are appended directly after the struct,
+in the listed order. Dispatch on `sample_sz` (72 / 152 / 296 at the tutorial config).
+
+### Payload: ping pong server `struct svr_sample` — 72 bytes, `<Qqqqqqdqd`
+
+svr_id, pings_sent, pings_recvd, pongs_sent, pongs_recvd, bytes_sent (all i64 deltas
+except u64 svr_id), rtt_sum (f64, ns), rtt_count (i64), end_time (f64 virtual time of the
+sample). No pointers, no appended arrays.
+
+### Payload: dally terminal `struct dfly_cn_ross_sample` — `128 + 16R + 8RQ` bytes
+
+| off | type | field |
+|---|---|---|
+| 0 | u64 | terminal_id |
+| 8 | i64 | fin_chunks |
+| 16 | i64 | data_size |
+| 24 | f64 | fin_hops (double here, unlike the GVT/RT payload) |
+| 32 | f64 | fin_chunks_time (ns) |
+| 40 | i64 | packets_gen |
+| 48 | i64 | packets_fin |
+| 56 | i64 | min_fin |
+| 64 | i64 | nonmin_fin |
+| 72 | f64 | end_time (virtual time of the sample) |
+| 80 | i64 | fwd_events (forward events processed) |
+| 88 | i64 | rev_events (events rolled back) |
+| 96 | 4×8 B | pointer slots — skip |
+| 128 | f64 ×R | busy_time (ns per rail) |
+| 128+8R | u64 ×R | stalled_chunks per rail |
+| 128+16R | i32 ×R·Q | vc_occupancy *snapshot*, rail-major |
+| 128+16R+4RQ | i32 ×R·Q | terminal_length *snapshot*, rail-major |
+
+Tutorial config: **152 bytes**.
+
+### Payload: dally router `struct dfly_router_ross_sample` — `72 + 32·radix` bytes
+
+| off | type | field |
+|---|---|---|
+| 0 | u64 | router_id |
+| 8 | f64 | end_time |
+| 16 | i64 | fwd_events |
+| 24 | i64 | rev_events |
+| 32 | 5×8 B | pointer slots — skip |
+| 72 | f64 ×radix | busy_time (ns per port) |
+| 72+8·radix | i64 ×radix | link_traffic (bytes per port) |
+| 72+16·radix | u64 ×radix | stalled_chunks per port |
+| 72+24·radix | i32 ×radix | vc_occupancy_sum *snapshot* |
+| 72+28·radix | i32 ×radix | queued_count *snapshot* |
+
+Tutorial config: **296 bytes**.
+
+## Mapping LP ids to model entities
+
+The first payload field of every record is the **model-relative id** (svr_id /
+terminal_id / router_id) — prefer it over the global `lpid`. For reference, the tutorial
+config lays LPs out in 36 repetitions of [2 `nw-lp`, 2 terminals, 1 router], so global
+`lpid = 5*rep + offset` with offsets 0–1 = servers, 2–3 = terminals, 4 = router
+(model-net's internal LP ordering within a repetition follows the config file's group
+listing order).
+
+## Format stability
+
+These payloads changed (fields appended, VT structs reworked) when the vis-oriented stats
+were added to dragonfly-dally; any parser of the pre-2026-07 dally payloads needs
+updating. Treat the layouts in this file as the source of truth and update it in the same
+PR as any payload change.
diff --git a/scripts/parse-ross-model-stats.py b/scripts/parse-ross-model-stats.py
new file mode 100755
index 00000000..02ab90bf
--- /dev/null
+++ b/scripts/parse-ross-model-stats.py
@@ -0,0 +1,367 @@
+#!/usr/bin/env python3
+"""Reference parser for the ROSS model-level stats binaries produced by CODES runs.
+
+Parses the two files written when a simulation runs with --model-stats:
+
+  <stats-path>/ross-stats-model.bin         GVT (stats_type=1) and real-time (2) samples
+  <stats-path>/ross-stats-analysis-lps.bin  virtual-time samples from the analysis LPs
+
+and decodes the payloads of the ping pong tutorial server LP and the dragonfly-dally
+terminal/router LPs. The full byte-level format is documented in
+doc/model-stats-binary-format.md -- keep the two in sync.
+
+Payloads are dispatched on their size, which depends on the network configuration; pass
+--num-rails/--num-qos/--radix for non-default configs (defaults match the ping pong
+tutorial config: 1 rail, 1 QoS level, radix 7). Unknown payload sizes are reported (and
+dumped as hex with --dump-unknown) but do not fail the parse.
+
+Output: one CSV per LP type to stdout (or to <prefix>-<type>.csv with --csv-prefix).
+Exits non-zero on framing errors or truncated input so tests can assert parseability.
+
+stdlib only -- no numpy/pandas required.
+"""
+
+import argparse
+import csv
+import struct
+import sys
+
+MODEL_TYPE_FLAG = 3  # lp_metadata/sample_metadata flag value for model data
+
+SAMPLE_METADATA = struct.Struct("<iidd")  # flag, sample_sz, ts, real_time
+MODEL_METADATA = struct.Struct("<IIIfiI")  # peid, kpid, lpid, gvt, stats_type, model_sz
+LP_METADATA = struct.Struct("<QQQddii")  # lpid, kpid, peid, ts, real_time, sample_sz, flag
+
+STATS_TYPE_NAMES = {1: "gvt", 2: "rt", 3: "vt"}
+
+
+class FormatError(Exception):
+    pass
+
+
+def _decode_svr_model(payload):
+    f = struct.Struct("<Q4IqdI")
+    vals = f.unpack_from(payload)
+    return dict(
+        zip(
+            (
+                "svr_id",
+                "pings_sent",
+                "pings_recvd",
+                "pongs_sent",
+                "pongs_recvd",
+                "bytes_sent",
+                "rtt_sum_ns",
+                "rtt_count",
+            ),
+            vals,
+        )
+    )
+
+
+def _decode_terminal_model(payload, rails, qos):
+    off = 0
+    row = {}
+    (row["terminal_id"], row["fin_chunks"], row["data_size"], row["fin_hops"]) = (
+        struct.unpack_from("<Q3q", payload, off)
+    )
+    off += 32
+    (row["fin_chunks_time_ns"],) = struct.unpack_from("<d", payload, off)
+    off += 8
+    for i in range(rails):
+        (row[f"busy_time_ns_r{i}"],) = struct.unpack_from("<d", payload, off)
+        off += 8
+    (row["packets_gen"], row["packets_fin"], row["min_fin"], row["nonmin_fin"]) = (
+        struct.unpack_from("<4q", payload, off)
+    )
+    off += 32
+    for i in range(rails):
+        (row[f"stalled_chunks_r{i}"],) = struct.unpack_from("<Q", payload, off)
+        off += 8
+    for i in range(rails):
+        for j in range(qos):
+            (row[f"vc_occupancy_r{i}q{j}"],) = struct.unpack_from("<i", payload, off)
+            off += 4
+    for i in range(rails):
+        for j in range(qos):
+            (row[f"terminal_length_r{i}q{j}"],) = struct.unpack_from("<i", payload, off)
+            off += 4
+    return row
+
+
+def _decode_router_model(payload, radix):
+    off = 0
+    row = {}
+    (row["router_id"],) = struct.unpack_from("<Q", payload, off)
+    off += 8
+    for i in range(radix):
+        row[f"busy_time_ns_p{i}"], row[f"link_traffic_p{i}"] = struct.unpack_from(
+            "<dq", payload, off
+        )
+        off += 16
+    for i in range(radix):
+        (row[f"stalled_chunks_p{i}"],) = struct.unpack_from("<Q", payload, off)
+        off += 8
+    for i in range(radix):
+        (row[f"vc_occupancy_sum_p{i}"],) = struct.unpack_from("<i", payload, off)
+        off += 4
+    for i in range(radix):
+        (row[f"queued_count_p{i}"],) = struct.unpack_from("<i", payload, off)
+        off += 4
+    return row
+
+
+def _decode_svr_vt(payload):
+    vals = struct.Struct("<Q5qdqd").unpack_from(payload)
+    return dict(
+        zip(
+            (
+                "svr_id",
+                "pings_sent",
+                "pings_recvd",
+                "pongs_sent",
+                "pongs_recvd",
+                "bytes_sent",
+                "rtt_sum_ns",
+                "rtt_count",
+                "end_time",
+            ),
+            vals,
+        )
+    )
+
+
+def _decode_terminal_vt(payload, rails, qos):
+    row = {}
+    (
+        row["terminal_id"],
+        row["fin_chunks"],
+        row["data_size"],
+        row["fin_hops"],
+        row["fin_chunks_time_ns"],
+        row["packets_gen"],
+        row["packets_fin"],
+        row["min_fin"],
+        row["nonmin_fin"],
+        row["end_time"],
+        row["fwd_events"],
+        row["rev_events"],
+    ) = struct.unpack_from("<Qqqddqqqqdqq", payload, 0)
+    off = 96 + 4 * 8  # scalars + 4 pointer slots (opaque on disk)
+    for i in range(rails):
+        (row[f"busy_time_ns_r{i}"],) = struct.unpack_from("<d", payload, off)
+        off += 8
+    for i in range(rails):
+        (row[f"stalled_chunks_r{i}"],) = struct.unpack_from("<Q", payload, off)
+        off += 8
+    for i in range(rails):
+        for j in range(qos):
+            (row[f"vc_occupancy_r{i}q{j}"],) = struct.unpack_from("<i", payload, off)
+            off += 4
+    for i in range(rails):
+        for j in range(qos):
+            (row[f"terminal_length_r{i}q{j}"],) = struct.unpack_from("<i", payload, off)
+            off += 4
+    return row
+
+
+def _decode_router_vt(payload, radix):
+    row = {}
+    (row["router_id"], row["end_time"], row["fwd_events"], row["rev_events"]) = (
+        struct.unpack_from("<Qdqq", payload, 0)
+    )
+    off = 32 + 5 * 8  # scalars + 5 pointer slots (opaque on disk)
+    for i in range(radix):
+        (row[f"busy_time_ns_p{i}"],) = struct.unpack_from("<d", payload, off)
+        off += 8
+    for i in range(radix):
+        (row[f"link_traffic_p{i}"],) = struct.unpack_from("<q", payload, off)
+        off += 8
+    for i in range(radix):
+        (row[f"stalled_chunks_p{i}"],) = struct.unpack_from("<Q", payload, off)
+        off += 8
+    for i in range(radix):
+        (row[f"vc_occupancy_sum_p{i}"],) = struct.unpack_from("<i", payload, off)
+        off += 4
+    for i in range(radix):
+        (row[f"queued_count_p{i}"],) = struct.unpack_from("<i", payload, off)
+        off += 4
+    return row
+
+
+def build_dispatch(rails, qos, radix):
+    """Map payload size -> (lp type name, decoder) for both files."""
+    model = {
+        44: ("svr", _decode_svr_model),
+        72 + 16 * rails + 8 * rails * qos: (
+            "terminal",
+            lambda p: _decode_terminal_model(p, rails, qos),
+        ),
+        8 + 32 * radix: ("router", lambda p: _decode_router_model(p, radix)),
+    }
+    vt = {
+        72: ("svr", _decode_svr_vt),
+        128 + 16 * rails + 8 * rails * qos: (
+            "terminal",
+            lambda p: _decode_terminal_vt(p, rails, qos),
+        ),
+        72 + 32 * radix: ("router", lambda p: _decode_router_vt(p, radix)),
+    }
+    if len(model) != 3 or len(vt) != 3:
+        raise FormatError(
+            "payload sizes collide for this rails/qos/radix combination; "
+            "cannot dispatch on size alone"
+        )
+    return model, vt
+
+
+def parse_model_file(path, dispatch, rows, unknown):
+    """Parse ross-stats-model.bin (GVT + RT records)."""
+    with open(path, "rb") as f:
+        data = f.read()
+    off = 0
+    n = 0
+    while off < len(data):
+        if len(data) - off < SAMPLE_METADATA.size + MODEL_METADATA.size:
+            raise FormatError(f"{path}: truncated record header at byte {off}")
+        flag, sample_sz, ts, real_time = SAMPLE_METADATA.unpack_from(data, off)
+        if flag != MODEL_TYPE_FLAG or sample_sz != MODEL_METADATA.size:
+            raise FormatError(
+                f"{path}: bad sample_metadata at byte {off} "
+                f"(flag={flag}, sample_sz={sample_sz}) -- lost framing"
+            )
+        off += SAMPLE_METADATA.size
+        peid, kpid, lpid, gvt, stats_type, model_sz = MODEL_METADATA.unpack_from(data, off)
+        off += MODEL_METADATA.size
+        if len(data) - off < model_sz:
+            raise FormatError(f"{path}: truncated payload at byte {off}")
+        payload = data[off : off + model_sz]
+        off += model_sz
+        n += 1
+        meta = {
+            "stats_type": STATS_TYPE_NAMES.get(stats_type, str(stats_type)),
+            "ts": ts,
+            "real_time": real_time,
+            "gvt": gvt,
+            "peid": peid,
+            "kpid": kpid,
+            "lpid": lpid,
+        }
+        entry = dispatch.get(model_sz)
+        if entry is None:
+            unknown.append((path, model_sz, payload))
+            continue
+        lp_type, decoder = entry
+        row = dict(meta)
+        row.update(decoder(payload))
+        rows.setdefault(lp_type, []).append(row)
+    return n
+
+
+def parse_vt_file(path, dispatch, rows, unknown):
+    """Parse ross-stats-analysis-lps.bin (virtual-time records)."""
+    with open(path, "rb") as f:
+        data = f.read()
+    off = 0
+    n = 0
+    while off < len(data):
+        if len(data) - off < LP_METADATA.size:
+            raise FormatError(f"{path}: truncated record header at byte {off}")
+        lpid, kpid, peid, ts, real_time, sample_sz, flag = LP_METADATA.unpack_from(data, off)
+        if flag != MODEL_TYPE_FLAG:
+            raise FormatError(
+                f"{path}: bad lp_metadata at byte {off} (flag={flag}) -- lost framing"
+            )
+        off += LP_METADATA.size
+        if len(data) - off < sample_sz:
+            raise FormatError(f"{path}: truncated payload at byte {off}")
+        payload = data[off : off + sample_sz]
+        off += sample_sz
+        n += 1
+        meta = {
+            "stats_type": "vt",
+            "ts": ts,
+            "real_time": real_time,
+            "gvt": "",
+            "peid": peid,
+            "kpid": kpid,
+            "lpid": lpid,
+        }
+        entry = dispatch.get(sample_sz)
+        if entry is None:
+            unknown.append((path, sample_sz, payload))
+            continue
+        lp_type, decoder = entry
+        row = dict(meta)
+        row.update(decoder(payload))
+        rows.setdefault(lp_type, []).append(row)
+    return n
+
+
+def write_csv(rows, csv_prefix):
+    for lp_type, entries in sorted(rows.items()):
+        entries.sort(key=lambda r: (r["ts"], r["lpid"]))
+        # union of keys across entries, preserving first-seen order (vt rows have extras)
+        fields = []
+        for r in entries:
+            for k in r:
+                if k not in fields:
+                    fields.append(k)
+        if csv_prefix:
+            out = open(f"{csv_prefix}-{lp_type}.csv", "w", newline="")
+        else:
+            out = sys.stdout
+            out.write(f"# --- {lp_type} ({len(entries)} records) ---\n")
+        writer = csv.DictWriter(out, fieldnames=fields, restval="")
+        writer.writeheader()
+        writer.writerows(entries)
+        if csv_prefix:
+            out.close()
+
+
+def main():
+    ap = argparse.ArgumentParser(
+        description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
+    )
+    ap.add_argument("files", nargs="+", help="ross-stats-model.bin and/or ross-stats-analysis-lps.bin")
+    ap.add_argument("--num-rails", type=int, default=1, help="dragonfly-dally num_rails (default 1)")
+    ap.add_argument("--num-qos", type=int, default=1, help="dragonfly-dally num_qos_levels (default 1)")
+    ap.add_argument("--radix", type=int, default=7, help="dragonfly-dally router radix (default 7)")
+    ap.add_argument("--csv-prefix", help="write <prefix>-<lptype>.csv files instead of stdout")
+    ap.add_argument("--dump-unknown", action="store_true", help="hex-dump undecodable payloads")
+    args = ap.parse_args()
+
+    model_dispatch, vt_dispatch = build_dispatch(args.num_rails, args.num_qos, args.radix)
+
+    rows = {}
+    unknown = []
+    total = 0
+    try:
+        for path in args.files:
+            if "analysis-lps" in path:
+                total += parse_vt_file(path, vt_dispatch, rows, unknown)
+            else:
+                total += parse_model_file(path, model_dispatch, rows, unknown)
+    except (FormatError, OSError) as e:
+        print(f"error: {e}", file=sys.stderr)
+        return 1
+
+    write_csv(rows, args.csv_prefix)
+
+    counts = {t: len(e) for t, e in sorted(rows.items())}
+    print(f"# parsed {total} records: {counts}, unknown payloads: {len(unknown)}", file=sys.stderr)
+    if unknown:
+        sizes = sorted({sz for (_, sz, _) in unknown})
+        print(
+            f"# undecoded payload sizes {sizes} -- check --num-rails/--num-qos/--radix "
+            "against the run's network config (see doc/model-stats-binary-format.md)",
+            file=sys.stderr,
+        )
+        if args.dump_unknown:
+            for path, sz, payload in unknown:
+                print(f"# {path} sz={sz}: {payload.hex()}", file=sys.stderr)
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/src/networks/model-net/dragonfly-dally.cxx b/src/networks/model-net/dragonfly-dally.cxx
index 26e12202..ac2212a6 100644
--- a/src/networks/model-net/dragonfly-dally.cxx
+++ b/src/networks/model-net/dragonfly-dally.cxx
@@ -470,6 +470,49 @@ struct dfly_cn_sample {
     long rev_events;
 };
 
+/* Sample structs for the ROSS instrumentation virtual-time (analysis LP) sampling mode
+ * (--model-stats=3/4). These are deliberately separate from dfly_cn_sample and
+ * dfly_router_sample above: those also define the on-disk format of the CODES-native
+ * sampling files (*-cn-sampling*.bin), so extending them would silently change that
+ * format too. Doubling as (a) the accumulator held in LP state (ross_sample /
+ * ross_rsample; array members heap-allocated, snapshot members unused/NULL) and (b) the
+ * record written to ross-stats-analysis-lps.bin (array members point directly after the
+ * struct in the analysis-LP buffer). All scalar members are 8 bytes so there are no
+ * padding holes; pointer members sit at the end of the struct and are meaningless on
+ * disk (parsers skip them). See doc/model-stats-binary-format.md. */
+struct dfly_cn_ross_sample {
+    tw_lpid terminal_id;
+    long fin_chunks_sample;
+    long data_size_sample;
+    double fin_hops_sample;
+    tw_stime fin_chunks_time;
+    int64_t packets_gen_sample;
+    int64_t packets_fin_sample;
+    int64_t min_fin_sample;    /* chunks that finished having taken a minimal path */
+    int64_t nonmin_fin_sample; /* chunks that finished having taken a non-minimal path */
+    tw_stime end_time;
+    long fwd_events;
+    long rev_events;
+    /* appended after the struct in the analysis-LP buffer, in this order: */
+    tw_stime* busy_time_sample;      /* [num_rails] */
+    uint64_t* stalled_chunks_sample; /* [num_rails] */
+    int* vc_occupancy;               /* [num_rails * num_qos_levels] snapshot */
+    int* terminal_length;            /* [num_rails * num_qos_levels] snapshot */
+};
+
+struct dfly_router_ross_sample {
+    tw_lpid router_id;
+    tw_stime end_time;
+    long fwd_events;
+    long rev_events;
+    /* appended after the struct in the analysis-LP buffer, in this order: */
+    tw_stime* busy_time;             /* [radix] */
+    int64_t* link_traffic_sample;    /* [radix] */
+    uint64_t* stalled_chunks_sample; /* [radix] */
+    int* vc_occupancy_sum;           /* [radix] snapshot, summed over VCs */
+    int* queued_count;               /* [radix] snapshot, bytes in overflow queues */
+};
+
 struct dfly_qhash_entry {
     struct dfly_hash_key key;
     char* remote_event_data;
@@ -744,7 +787,12 @@ struct terminal_state {
     long fin_hops_ross_sample;
     tw_stime fin_chunks_time_ross_sample;
     tw_stime* busy_time_ross_sample;
-    struct dfly_cn_sample ross_sample;
+    int64_t packets_gen_ross_sample;
+    int64_t packets_fin_ross_sample;
+    int64_t min_fin_ross_sample;
+    int64_t nonmin_fin_ross_sample;
+    uint64_t* stalled_chunks_ross_sample; //[rail_id]
+    struct dfly_cn_ross_sample ross_sample;
 
     // Variables to recover latency of packets sent to other terminals
     // Sent packets (to be populated at by commit handler of packet sender)
@@ -834,7 +882,8 @@ struct router_state {
     /* following used for ROSS model-level stats collection */
     tw_stime* busy_time_ross_sample;
     int64_t* link_traffic_ross_sample;
-    struct dfly_router_sample ross_rsample;
+    uint64_t* stalled_chunks_ross_sample; //[port]
+    struct dfly_router_ross_sample ross_rsample;
     tw_stime last_time;
 };
 
@@ -845,20 +894,21 @@ void custom_dally_dragonfly_event_collect(terminal_dally_message* m, tw_lp* lp,
 void custom_dally_dragonfly_model_stat_collect(terminal_state* s, tw_lp* lp, char* buffer);
 void custom_dally_dfly_router_model_stat_collect(router_state* s, tw_lp* lp, char* buffer);
 static void ross_dally_dragonfly_rsample_fn(router_state* s, tw_bf* bf, tw_lp* lp,
-                                            struct dfly_router_sample* sample);
+                                            struct dfly_router_ross_sample* sample);
 static void ross_dally_dragonfly_rsample_rc_fn(router_state* s, tw_bf* bf, tw_lp* lp,
-                                               struct dfly_router_sample* sample);
+                                               struct dfly_router_ross_sample* sample);
 static void ross_dally_dragonfly_sample_fn(terminal_state* s, tw_bf* bf, tw_lp* lp,
-                                           struct dfly_cn_sample* sample);
+                                           struct dfly_cn_ross_sample* sample);
 static void ross_dally_dragonfly_sample_rc_fn(terminal_state* s, tw_bf* bf, tw_lp* lp,
-                                              struct dfly_cn_sample* sample);
+                                              struct dfly_cn_ross_sample* sample);
 
 st_model_types custom_dally_dragonfly_model_types[] = {
     {(ev_trace_f)custom_dally_dragonfly_event_collect, sizeof(int),
      (model_stat_f)custom_dally_dragonfly_model_stat_collect,
-     sizeof(tw_lpid) + sizeof(long) * 2 + sizeof(double) + sizeof(tw_stime) * 2,
+     0, //updated in terminal_dally_init() since it's based on num_rails/num_qos_levels
      (sample_event_f)ross_dally_dragonfly_sample_fn,
-     (sample_revent_f)ross_dally_dragonfly_sample_rc_fn, sizeof(struct dfly_cn_sample)},
+     (sample_revent_f)ross_dally_dragonfly_sample_rc_fn,
+     0}, //updated in terminal_dally_init() since it's based on num_rails/num_qos_levels
     {(ev_trace_f)custom_dally_dragonfly_event_collect, sizeof(int),
      (model_stat_f)custom_dally_dfly_router_model_stat_collect,
      0, //updated in router_dally_init() since it's based on the radix
@@ -881,10 +931,20 @@ void custom_dally_dragonfly_event_collect(terminal_dally_message* m, tw_lp* lp,
 void custom_dally_dragonfly_model_stat_collect(terminal_state* s, tw_lp* lp, char* buffer) {
     (void)lp;
 
+    const dragonfly_param* p = s->params;
+    int const num_qos_levels = p->num_qos_levels;
     int index = 0;
     tw_lpid id = 0;
     long tmp = 0;
     tw_stime tmp2 = 0;
+    int64_t tmp3 = 0;
+    uint64_t tmp4 = 0;
+    int tmp5 = 0;
+
+    /* NOTE on NULL guards: while a terminal is frozen in surrogate mode its state is
+     * memset and only select members are carried over (see
+     * dragonfly_dally_terminal_highdef_to_surrogate), so pointer members here may be NULL
+     * mid-run; emit zeros for those fields rather than crashing GVT/RT collection. */
 
     id = s->terminal_id;
     memcpy(&buffer[index], &id, sizeof(id));
@@ -910,11 +970,56 @@ void custom_dally_dragonfly_model_stat_collect(terminal_state* s, tw_lp* lp, cha
     index += sizeof(tmp2);
     s->fin_chunks_time_ross_sample = 0;
 
-    for (int i = 0; i < s->params->num_rails; i++) {
-        tmp2 = s->busy_time_ross_sample[i];
+    for (int i = 0; i < p->num_rails; i++) {
+        tmp2 = s->busy_time_ross_sample ? s->busy_time_ross_sample[i] : 0;
         memcpy(&buffer[index], &tmp2, sizeof(tmp2));
         index += sizeof(tmp2);
-        s->busy_time_ross_sample[i] = 0;
+        if (s->busy_time_ross_sample)
+            s->busy_time_ross_sample[i] = 0;
+    }
+
+    tmp3 = s->packets_gen_ross_sample;
+    memcpy(&buffer[index], &tmp3, sizeof(tmp3));
+    index += sizeof(tmp3);
+    s->packets_gen_ross_sample = 0;
+
+    tmp3 = s->packets_fin_ross_sample;
+    memcpy(&buffer[index], &tmp3, sizeof(tmp3));
+    index += sizeof(tmp3);
+    s->packets_fin_ross_sample = 0;
+
+    tmp3 = s->min_fin_ross_sample;
+    memcpy(&buffer[index], &tmp3, sizeof(tmp3));
+    index += sizeof(tmp3);
+    s->min_fin_ross_sample = 0;
+
+    tmp3 = s->nonmin_fin_ross_sample;
+    memcpy(&buffer[index], &tmp3, sizeof(tmp3));
+    index += sizeof(tmp3);
+    s->nonmin_fin_ross_sample = 0;
+
+    for (int i = 0; i < p->num_rails; i++) {
+        tmp4 = s->stalled_chunks_ross_sample ? s->stalled_chunks_ross_sample[i] : 0;
+        memcpy(&buffer[index], &tmp4, sizeof(tmp4));
+        index += sizeof(tmp4);
+        if (s->stalled_chunks_ross_sample)
+            s->stalled_chunks_ross_sample[i] = 0;
+    }
+
+    /* instantaneous occupancy/backlog snapshots (not deltas, so nothing is zeroed) */
+    for (int i = 0; i < p->num_rails; i++) {
+        for (int j = 0; j < num_qos_levels; j++) {
+            tmp5 = s->vc_occupancy ? s->vc_occupancy[i][j] : 0;
+            memcpy(&buffer[index], &tmp5, sizeof(tmp5));
+            index += sizeof(tmp5);
+        }
+    }
+    for (int i = 0; i < p->num_rails; i++) {
+        for (int j = 0; j < num_qos_levels; j++) {
+            tmp5 = s->terminal_length ? s->terminal_length[i][j] : 0;
+            memcpy(&buffer[index], &tmp5, sizeof(tmp5));
+            index += sizeof(tmp5);
+        }
     }
 
     return;
@@ -929,6 +1034,8 @@ void custom_dally_dfly_router_model_stat_collect(router_state* s, tw_lp* lp, cha
     tw_lpid id = 0;
     tw_stime tmp = 0;
     int64_t tmp2 = 0;
+    uint64_t tmp3 = 0;
+    int tmp4 = 0;
 
     id = s->router_id;
     memcpy(&buffer[index], &id, sizeof(id));
@@ -945,6 +1052,27 @@ void custom_dally_dfly_router_model_stat_collect(router_state* s, tw_lp* lp, cha
         index += sizeof(tmp2);
         s->link_traffic_ross_sample[i] = 0;
     }
+
+    for (i = 0; i < p->radix; i++) {
+        tmp3 = s->stalled_chunks_ross_sample[i];
+        memcpy(&buffer[index], &tmp3, sizeof(tmp3));
+        index += sizeof(tmp3);
+        s->stalled_chunks_ross_sample[i] = 0;
+    }
+
+    /* instantaneous per-port snapshots (not deltas, so nothing is zeroed) */
+    for (i = 0; i < p->radix; i++) {
+        tmp4 = 0;
+        for (int j = 0; j < p->num_vcs; j++)
+            tmp4 += s->vc_occupancy[i][j];
+        memcpy(&buffer[index], &tmp4, sizeof(tmp4));
+        index += sizeof(tmp4);
+    }
+    for (i = 0; i < p->radix; i++) {
+        tmp4 = s->queued_count[i];
+        memcpy(&buffer[index], &tmp4, sizeof(tmp4));
+        index += sizeof(tmp4);
+    }
     return;
 }
 
@@ -966,7 +1094,7 @@ static void custom_dally_router_register_model_types(st_model_types* base_type)
 /*** END of ROSS event tracing additions */
 
 static void ross_dally_dragonfly_rsample_fn(router_state* s, tw_bf* bf, tw_lp* lp,
-                                            struct dfly_router_sample* sample) {
+                                            struct dfly_router_ross_sample* sample) {
     (void)lp;
     (void)bf;
 
@@ -977,12 +1105,22 @@ static void ross_dally_dragonfly_rsample_fn(router_state* s, tw_bf* bf, tw_lp* l
     sample->end_time = tw_now(lp);
     sample->fwd_events = s->ross_rsample.fwd_events;
     sample->rev_events = s->ross_rsample.rev_events;
-    sample->busy_time = (tw_stime*)((&sample->rev_events) + 1);
-    sample->link_traffic_sample = (int64_t*)((&sample->busy_time[0]) + p->radix);
+    /* variable-length arrays live directly after the struct in the analysis-LP buffer */
+    sample->busy_time = (tw_stime*)(sample + 1);
+    sample->link_traffic_sample = (int64_t*)(sample->busy_time + p->radix);
+    sample->stalled_chunks_sample = (uint64_t*)(sample->link_traffic_sample + p->radix);
+    sample->vc_occupancy_sum = (int*)(sample->stalled_chunks_sample + p->radix);
+    sample->queued_count = sample->vc_occupancy_sum + p->radix;
 
     for (; i < p->radix; i++) {
         sample->busy_time[i] = s->ross_rsample.busy_time[i];
         sample->link_traffic_sample[i] = s->ross_rsample.link_traffic_sample[i];
+        sample->stalled_chunks_sample[i] = s->ross_rsample.stalled_chunks_sample[i];
+        /* instantaneous per-port snapshots read straight from live state */
+        sample->vc_occupancy_sum[i] = 0;
+        for (int j = 0; j < p->num_vcs; j++)
+            sample->vc_occupancy_sum[i] += s->vc_occupancy[i][j];
+        sample->queued_count[i] = s->queued_count[i];
     }
 
     /* clear up the current router stats */
@@ -992,20 +1130,24 @@ static void ross_dally_dragonfly_rsample_fn(router_state* s, tw_bf* bf, tw_lp* l
     for (i = 0; i < p->radix; i++) {
         s->ross_rsample.busy_time[i] = 0;
         s->ross_rsample.link_traffic_sample[i] = 0;
+        s->ross_rsample.stalled_chunks_sample[i] = 0;
     }
 }
 
 static void ross_dally_dragonfly_rsample_rc_fn(router_state* s, tw_bf* bf, tw_lp* lp,
-                                               struct dfly_router_sample* sample) {
+                                               struct dfly_router_ross_sample* sample) {
     (void)lp;
     (void)bf;
 
     const dragonfly_param* p = s->params;
     int i = 0;
 
+    /* the snapshot arrays (vc_occupancy_sum, queued_count) touched no LP state, so there
+     * is nothing to restore for them */
     for (; i < p->radix; i++) {
         s->ross_rsample.busy_time[i] = sample->busy_time[i];
         s->ross_rsample.link_traffic_sample[i] = sample->link_traffic_sample[i];
+        s->ross_rsample.stalled_chunks_sample[i] = sample->stalled_chunks_sample[i];
     }
 
     s->ross_rsample.fwd_events = sample->fwd_events;
@@ -1013,43 +1155,88 @@ static void ross_dally_dragonfly_rsample_rc_fn(router_state* s, tw_bf* bf, tw_lp
 }
 
 static void ross_dally_dragonfly_sample_fn(terminal_state* s, tw_bf* bf, tw_lp* lp,
-                                           struct dfly_cn_sample* sample) {
+                                           struct dfly_cn_ross_sample* sample) {
     (void)lp;
     (void)bf;
 
+    const dragonfly_param* p = s->params;
+    int const num_qos_levels = p->num_qos_levels;
+
+    /* NULL guards: a terminal frozen in surrogate mode has most pointer members unset
+     * (see dragonfly_dally_terminal_highdef_to_surrogate); emit zeros for those fields */
+
     sample->terminal_id = s->terminal_id;
     sample->fin_chunks_sample = s->ross_sample.fin_chunks_sample;
     sample->data_size_sample = s->ross_sample.data_size_sample;
     sample->fin_hops_sample = s->ross_sample.fin_hops_sample;
     sample->fin_chunks_time = s->ross_sample.fin_chunks_time;
-    sample->busy_time_sample = (tw_stime*)((&sample->busy_time_sample[0] + s->params->num_rails));
-    for (int i = 0; i < s->params->num_rails; i++)
-        sample->busy_time_sample[i] = s->ross_sample.busy_time_sample[i];
+    sample->packets_gen_sample = s->ross_sample.packets_gen_sample;
+    sample->packets_fin_sample = s->ross_sample.packets_fin_sample;
+    sample->min_fin_sample = s->ross_sample.min_fin_sample;
+    sample->nonmin_fin_sample = s->ross_sample.nonmin_fin_sample;
     sample->end_time = tw_now(lp);
     sample->fwd_events = s->ross_sample.fwd_events;
     sample->rev_events = s->ross_sample.rev_events;
 
+    /* variable-length arrays live directly after the struct in the analysis-LP buffer */
+    sample->busy_time_sample = (tw_stime*)(sample + 1);
+    sample->stalled_chunks_sample = (uint64_t*)(sample->busy_time_sample + p->num_rails);
+    sample->vc_occupancy = (int*)(sample->stalled_chunks_sample + p->num_rails);
+    sample->terminal_length = sample->vc_occupancy + p->num_rails * num_qos_levels;
+
+    for (int i = 0; i < p->num_rails; i++) {
+        sample->busy_time_sample[i] =
+            s->ross_sample.busy_time_sample ? s->ross_sample.busy_time_sample[i] : 0;
+        sample->stalled_chunks_sample[i] =
+            s->ross_sample.stalled_chunks_sample ? s->ross_sample.stalled_chunks_sample[i] : 0;
+        /* instantaneous snapshots read straight from live state */
+        for (int j = 0; j < num_qos_levels; j++) {
+            sample->vc_occupancy[i * num_qos_levels + j] =
+                s->vc_occupancy ? s->vc_occupancy[i][j] : 0;
+            sample->terminal_length[i * num_qos_levels + j] =
+                s->terminal_length ? s->terminal_length[i][j] : 0;
+        }
+    }
+
     s->ross_sample.fin_chunks_sample = 0;
     s->ross_sample.data_size_sample = 0;
     s->ross_sample.fin_hops_sample = 0;
+    s->ross_sample.packets_gen_sample = 0;
+    s->ross_sample.packets_fin_sample = 0;
+    s->ross_sample.min_fin_sample = 0;
+    s->ross_sample.nonmin_fin_sample = 0;
     s->ross_sample.fwd_events = 0;
     s->ross_sample.rev_events = 0;
     s->ross_sample.fin_chunks_time = 0;
-    for (int i = 0; i < s->params->num_rails; i++)
-        s->ross_sample.busy_time_sample[i] = 0;
+    for (int i = 0; i < p->num_rails; i++) {
+        if (s->ross_sample.busy_time_sample)
+            s->ross_sample.busy_time_sample[i] = 0;
+        if (s->ross_sample.stalled_chunks_sample)
+            s->ross_sample.stalled_chunks_sample[i] = 0;
+    }
 }
 
 static void ross_dally_dragonfly_sample_rc_fn(terminal_state* s, tw_bf* bf, tw_lp* lp,
-                                              struct dfly_cn_sample* sample) {
+                                              struct dfly_cn_ross_sample* sample) {
     (void)lp;
     (void)bf;
 
-    for (int i = 0; i < s->params->num_rails; i++)
-        s->ross_sample.busy_time_sample[i] = sample->busy_time_sample[i];
+    /* the snapshot arrays (vc_occupancy, terminal_length) touched no LP state, so there
+     * is nothing to restore for them */
+    for (int i = 0; i < s->params->num_rails; i++) {
+        if (s->ross_sample.busy_time_sample)
+            s->ross_sample.busy_time_sample[i] = sample->busy_time_sample[i];
+        if (s->ross_sample.stalled_chunks_sample)
+            s->ross_sample.stalled_chunks_sample[i] = sample->stalled_chunks_sample[i];
+    }
     s->ross_sample.fin_chunks_time = sample->fin_chunks_time;
     s->ross_sample.fin_hops_sample = sample->fin_hops_sample;
     s->ross_sample.data_size_sample = sample->data_size_sample;
     s->ross_sample.fin_chunks_sample = sample->fin_chunks_sample;
+    s->ross_sample.packets_gen_sample = sample->packets_gen_sample;
+    s->ross_sample.packets_fin_sample = sample->packets_fin_sample;
+    s->ross_sample.min_fin_sample = sample->min_fin_sample;
+    s->ross_sample.nonmin_fin_sample = sample->nonmin_fin_sample;
     s->ross_sample.fwd_events = sample->fwd_events;
     s->ross_sample.rev_events = sample->rev_events;
 }
@@ -4038,6 +4225,23 @@ static void dragonfly_dally_terminal_highdef_to_surrogate(terminal_state* s, tw_
     s->data_size_sample = frozen_state->data_size_sample;
     s->ross_sample.data_size_sample = frozen_state->ross_sample.data_size_sample;
     s->data_size_ross_sample = frozen_state->data_size_ross_sample;
+    s->packets_gen_ross_sample = frozen_state->packets_gen_ross_sample;
+    s->packets_fin_ross_sample = frozen_state->packets_fin_ross_sample;
+    s->min_fin_ross_sample = frozen_state->min_fin_ross_sample;
+    s->nonmin_fin_ross_sample = frozen_state->nonmin_fin_ross_sample;
+    s->ross_sample.packets_gen_sample = frozen_state->ross_sample.packets_gen_sample;
+    s->ross_sample.packets_fin_sample = frozen_state->ross_sample.packets_fin_sample;
+    s->ross_sample.min_fin_sample = frozen_state->ross_sample.min_fin_sample;
+    s->ross_sample.nonmin_fin_sample = frozen_state->ross_sample.nonmin_fin_sample;
+    /* keep the instrumentation delta buffers alive while frozen so ROSS model-stats
+     * collection (--model-stats) can run without dereferencing NULL; the hidden network
+     * state (vc_occupancy, terminal_length, ...) is intentionally left unset and
+     * collection emits zeros for it (see the NULL guards in
+     * custom_dally_dragonfly_model_stat_collect / ross_dally_dragonfly_sample_fn) */
+    s->busy_time_ross_sample = frozen_state->busy_time_ross_sample;
+    s->stalled_chunks_ross_sample = frozen_state->stalled_chunks_ross_sample;
+    s->ross_sample.busy_time_sample = frozen_state->ross_sample.busy_time_sample;
+    s->ross_sample.stalled_chunks_sample = frozen_state->ross_sample.stalled_chunks_sample;
     s->total_msg_size = frozen_state->total_msg_size;
     s->finished_msgs = frozen_state->finished_msgs;
     s->rank_tbl_pop = frozen_state->rank_tbl_pop;
@@ -4083,6 +4287,14 @@ static void dragonfly_dally_terminal_surrogate_to_highdef(terminal_state* s, tw_
     frozen_state->data_size_sample = s->data_size_sample;
     frozen_state->ross_sample.data_size_sample = s->ross_sample.data_size_sample;
     frozen_state->data_size_ross_sample = s->data_size_ross_sample;
+    frozen_state->packets_gen_ross_sample = s->packets_gen_ross_sample;
+    frozen_state->packets_fin_ross_sample = s->packets_fin_ross_sample;
+    frozen_state->min_fin_ross_sample = s->min_fin_ross_sample;
+    frozen_state->nonmin_fin_ross_sample = s->nonmin_fin_ross_sample;
+    frozen_state->ross_sample.packets_gen_sample = s->ross_sample.packets_gen_sample;
+    frozen_state->ross_sample.packets_fin_sample = s->ross_sample.packets_fin_sample;
+    frozen_state->ross_sample.min_fin_sample = s->ross_sample.min_fin_sample;
+    frozen_state->ross_sample.nonmin_fin_sample = s->ross_sample.nonmin_fin_sample;
     frozen_state->total_msg_size = s->total_msg_size;
     frozen_state->finished_msgs = s->finished_msgs;
     frozen_state->rank_tbl_pop = s->rank_tbl_pop;
@@ -4622,8 +4834,31 @@ static void terminal_dally_init(terminal_state* s, tw_lp* lp) {
 
 
     s->ross_sample.busy_time_sample = (tw_stime*)calloc(s->params->num_rails, sizeof(tw_stime));
+    s->ross_sample.stalled_chunks_sample =
+        (uint64_t*)calloc(s->params->num_rails, sizeof(uint64_t));
     s->busy_time_ross_sample = (tw_stime*)calloc(s->params->num_rails, sizeof(tw_stime));
+    s->stalled_chunks_ross_sample = (uint64_t*)calloc(s->params->num_rails, sizeof(uint64_t));
     s->busy_time_sample = (tw_stime*)calloc(s->params->num_rails, sizeof(tw_stime));
+    s->packets_gen_ross_sample = 0;
+    s->packets_fin_ross_sample = 0;
+    s->min_fin_ross_sample = 0;
+    s->nonmin_fin_ross_sample = 0;
+
+    /* The model-stats payload and virtual-time sample sizes depend on num_rails and
+     * num_qos_levels, so they are set here at runtime rather than in the static
+     * custom_dally_dragonfly_model_types initializer (all terminals share one params set
+     * per annotation, so every terminal writes the same values). Keep these formulas in
+     * sync with custom_dally_dragonfly_model_stat_collect() and
+     * ross_dally_dragonfly_sample_fn(). */
+    if (g_st_model_stats)
+        lp->model_types->mstat_sz = sizeof(tw_lpid) + 3 * sizeof(long) + sizeof(tw_stime) +
+                                    4 * sizeof(int64_t) +
+                                    p->num_rails * (sizeof(tw_stime) + sizeof(uint64_t)) +
+                                    2 * p->num_rails * num_qos_levels * sizeof(int);
+    if (g_st_use_analysis_lps && g_st_model_stats)
+        lp->model_types->sample_struct_sz = sizeof(struct dfly_cn_ross_sample) +
+                                            p->num_rails * (sizeof(tw_stime) + sizeof(uint64_t)) +
+                                            2 * p->num_rails * num_qos_levels * sizeof(int);
 
     /*if(s->terminal_id == 0)
         {
@@ -4768,17 +5003,22 @@ static void router_dally_init(router_state* r, tw_lp* lp) {
     r->busy_time = (tw_stime*)calloc(p->radix, sizeof(tw_stime));
     r->busy_time_sample = (tw_stime*)calloc(p->radix, sizeof(tw_stime));
 
-    /* set up for ROSS stats sampling */
+    /* set up for ROSS stats sampling; keep the size formulas in sync with
+     * custom_dally_dfly_router_model_stat_collect() and ross_dally_dragonfly_rsample_fn() */
     r->link_traffic_ross_sample = (int64_t*)calloc(p->radix, sizeof(int64_t));
     r->busy_time_ross_sample = (tw_stime*)calloc(p->radix, sizeof(tw_stime));
+    r->stalled_chunks_ross_sample = (uint64_t*)calloc(p->radix, sizeof(uint64_t));
     if (g_st_model_stats)
         lp->model_types->mstat_sz =
-            sizeof(tw_lpid) + (sizeof(int64_t) + sizeof(tw_stime)) * p->radix;
+            sizeof(tw_lpid) +
+            (sizeof(int64_t) + sizeof(tw_stime) + sizeof(uint64_t) + 2 * sizeof(int)) * p->radix;
     if (g_st_use_analysis_lps && g_st_model_stats)
         lp->model_types->sample_struct_sz =
-            sizeof(struct dfly_router_sample) + (sizeof(tw_stime) + sizeof(int64_t)) * p->radix;
+            sizeof(struct dfly_router_ross_sample) +
+            (sizeof(tw_stime) + sizeof(int64_t) + sizeof(uint64_t) + 2 * sizeof(int)) * p->radix;
     r->ross_rsample.busy_time = (tw_stime*)calloc(p->radix, sizeof(tw_stime));
     r->ross_rsample.link_traffic_sample = (int64_t*)calloc(p->radix, sizeof(int64_t));
+    r->ross_rsample.stalled_chunks_sample = (uint64_t*)calloc(p->radix, sizeof(uint64_t));
 
     rc_stack_create(&r->st);
     rc_stack_create(&r->cc_st);
@@ -4956,6 +5196,8 @@ static void packet_generate_predicted_rc(terminal_state* s, tw_bf* bf, terminal_
     s->packet_counter--;
     s->total_gen_size -= msg->packet_size;
     s->packet_gen--;
+    s->packets_gen_ross_sample--;
+    s->ross_sample.packets_gen_sample--;
     packet_gen--;
 }
 
@@ -4964,6 +5206,8 @@ static void packet_generate_predicted(terminal_state* s, tw_bf* bf, terminal_dal
                                       tw_lp* lp) {
     packet_gen++;
     s->packet_gen++;
+    s->packets_gen_ross_sample++;
+    s->ross_sample.packets_gen_sample++;
     s->total_gen_size += msg->packet_size;
 
     assert(lp->gid != msg->dest_terminal_lpid);
@@ -5108,6 +5352,8 @@ static void packet_generate_rc(terminal_state* s, tw_bf* bf, terminal_dally_mess
 
     s->total_gen_size -= msg->packet_size;
     s->packet_gen--;
+    s->packets_gen_ross_sample--;
+    s->ross_sample.packets_gen_sample--;
     packet_gen--;
     s->packet_counter--;
 
@@ -5172,6 +5418,8 @@ static void packet_generate(terminal_state* s, tw_bf* bf, terminal_dally_message
     packet_gen++;
 
     s->packet_gen++;
+    s->packets_gen_ross_sample++;
+    s->ross_sample.packets_gen_sample++;
     s->total_gen_size += msg->packet_size;
     msg->saved_processing_time = tw_now(lp);
 
@@ -5540,6 +5788,8 @@ static void packet_send_rc(terminal_state* s, tw_bf* bf, terminal_dally_message*
     if (bf->c1) {
         s->in_send_loop[msg->rail_id] = msg->saved_send_loop;
         s->stalled_chunks[msg->rail_id]--;
+        s->stalled_chunks_ross_sample[msg->rail_id]--;
+        s->ross_sample.stalled_chunks_sample[msg->rail_id]--;
         if (bf->c3)
             s->last_buf_full[msg->rail_id] = msg->saved_busy_time;
 
@@ -5605,6 +5855,8 @@ static void packet_send(terminal_state* s, tw_bf* bf, terminal_dally_message* ms
         bf->c1 = 1;
         s->in_send_loop[msg->rail_id] = 0;
         s->stalled_chunks[msg->rail_id]++;
+        s->stalled_chunks_ross_sample[msg->rail_id]++;
+        s->ross_sample.stalled_chunks_sample[msg->rail_id]++;
         if (!s->last_buf_full[msg->rail_id]) {
             bf->c3 = 1;
             msg->saved_busy_time = s->last_buf_full[msg->rail_id];
@@ -6069,10 +6321,16 @@ static void packet_arrive_rc(terminal_state* s, tw_bf* bf, terminal_dally_messag
     if (g_congestion_control_enabled)
         cc_terminal_send_ack_rc(s->local_congestion_controller);
 
-    if (msg->path_type == MINIMAL)
+    if (msg->path_type == MINIMAL) {
         minimal_count--;
-    if (msg->path_type == NON_MINIMAL)
+        s->min_fin_ross_sample--;
+        s->ross_sample.min_fin_sample--;
+    }
+    if (msg->path_type == NON_MINIMAL) {
         nonmin_count--;
+        s->nonmin_fin_ross_sample--;
+        s->ross_sample.nonmin_fin_sample--;
+    }
 
     N_finished_chunks--;
     s->finished_chunks--;
@@ -6112,6 +6370,8 @@ static void packet_arrive_rc(terminal_state* s, tw_bf* bf, terminal_dally_messag
 
     if (bf->c1) {
         s->packet_fin--;
+        s->packets_fin_ross_sample--;
+        s->ross_sample.packets_fin_sample--;
         packet_fin--;
         stat->recv_count--;
         stat->recv_bytes -= msg->packet_size;
@@ -6249,11 +6509,17 @@ static void packet_arrive(terminal_state* s, tw_bf* bf, terminal_dally_message*
 
     uint64_t const num_chunks = num_chunks_for(msg->packet_size, s->params->chunk_size);
 
-    if (msg->path_type == MINIMAL)
+    if (msg->path_type == MINIMAL) {
         minimal_count++;
+        s->min_fin_ross_sample++;
+        s->ross_sample.min_fin_sample++;
+    }
 
-    if (msg->path_type == NON_MINIMAL)
+    if (msg->path_type == NON_MINIMAL) {
         nonmin_count++;
+        s->nonmin_fin_ross_sample++;
+        s->ross_sample.nonmin_fin_sample++;
+    }
 
     if (msg->path_type != MINIMAL && msg->path_type != NON_MINIMAL)
         printf("\n Wrong message path type %d ", msg->path_type);
@@ -6286,6 +6552,8 @@ static void packet_arrive(terminal_state* s, tw_bf* bf, terminal_dally_message*
     if (msg->chunk_id == num_chunks - 1) {
         bf->c1 = 1;
         s->packet_fin++;
+        s->packets_fin_ross_sample++;
+        s->ross_sample.packets_fin_sample++;
         packet_fin++;
 
         stat->recv_count++;
@@ -7082,6 +7350,8 @@ static void router_packet_receive_rc(router_state* s, tw_bf* bf, terminal_dally_
     }
     if (bf->c4) {
         s->stalled_chunks[output_port]--;
+        s->stalled_chunks_ross_sample[output_port]--;
+        s->ross_rsample.stalled_chunks_sample[output_port]--;
         if (bf->c22) {
             s->last_buf_full[output_port] = msg->saved_busy_time;
         }
@@ -7288,6 +7558,8 @@ static void router_packet_receive(router_state* s, tw_bf* bf, terminal_dally_mes
     } else {
         bf->c4 = 1;
         s->stalled_chunks[output_port]++;
+        s->stalled_chunks_ross_sample[output_port]++;
+        s->ross_rsample.stalled_chunks_sample[output_port]++;
         cur_chunk->msg.saved_vc = msg->vc_index;
         cur_chunk->msg.saved_channel = msg->output_chan;
         assert(output_chan < s->params->num_vcs && output_port < s->params->radix);
@@ -8232,6 +8504,7 @@ static void save_terminal_state(terminal_state* into, terminal_state const* from
     // from->sample_stat
     // from->ross_sample
     // from->busy_time_ross_sample
+    // from->stalled_chunks_ross_sample
     // from->busy_time_sample
 
     memcpy(into, from, sizeof(terminal_state));
@@ -8410,6 +8683,10 @@ static bool check_terminal_state(terminal_state* before, terminal_state* after)
     is_same &= (before->data_size_ross_sample == after->data_size_ross_sample);
     is_same &= (before->fin_hops_ross_sample == after->fin_hops_ross_sample);
     is_same &= (before->fin_chunks_time_ross_sample == after->fin_chunks_time_ross_sample);
+    is_same &= (before->packets_gen_ross_sample == after->packets_gen_ross_sample);
+    is_same &= (before->packets_fin_ross_sample == after->packets_fin_ross_sample);
+    is_same &= (before->min_fin_ross_sample == after->min_fin_ross_sample);
+    is_same &= (before->nonmin_fin_ross_sample == after->nonmin_fin_ross_sample);
     is_same &= (before->last_in_queue_time == after->last_in_queue_time);
 
     // Compare string buffers
@@ -8718,7 +8995,7 @@ static void print_terminal_state(FILE* out, char const* prefix, terminal_state*
         state
             ->busy_time_ross_sample); // ingnoring as this part of the code is never used. Originally part of instrumentation
     fprintf(
-        out, "%s  |              ross_sample = <dfly_cn_sample object>\n",
+        out, "%s  |              ross_sample = <dfly_cn_ross_sample object>\n",
         prefix); // ingnoring as this part of the code is never used. Originally part of instrumentation
 
     // modified outside of process and reverse computation (at commit and at surrogate change)
@@ -8762,6 +9039,7 @@ static void save_router_state(router_state* into, router_state const* from) {
     // from->busy_time_sample
     // from->link_traffic_sample
     // from->link_traffic_ross_sample
+    // from->stalled_chunks_ross_sample
 
     memcpy(into, from, sizeof(router_state));
 
@@ -9217,7 +9495,7 @@ static void print_router_state(FILE* out, char const* prefix, router_state* stat
         fprintf(out, "%s  |  ]\n", prefix);
     }
 
-    fprintf(out, "%s  |             ross_rsample = <dfly_router_sample object>\n", prefix);
+    fprintf(out, "%s  |             ross_rsample = <dfly_router_ross_sample object>\n", prefix);
     fprintf(out, "%s  |                last_time = %g\n", prefix, state->last_time);
 }
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 9c586328..e7eb2b2d 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -357,6 +357,11 @@ set(test-shell-files
     example-ping-pong-surrogate-1.sh
     example-ping-pong-surrogate-2.sh
     example-ping-pong-surrogate-3.sh
+    # ROSS instrumentation model-stats smoke tests: custom assertions on the
+    # binary stats output (file presence + reference-parser round trip), so
+    # they stay as scripts.
+    example-ping-pong-model-stats.sh
+    example-ping-pong-model-stats-modes.sh
     )
 
 configure_file(conf/union-milc-jacobi-workload/dfdally-72-par.conf.in conf/union-milc-jacobi-workload/dfdally-72-par.conf.in @ONLY)
diff --git a/tests/example-ping-pong-model-stats-modes.sh b/tests/example-ping-pong-model-stats-modes.sh
new file mode 100755
index 00000000..8564f328
--- /dev/null
+++ b/tests/example-ping-pong-model-stats-modes.sh
@@ -0,0 +1,74 @@
+#!/bin/bash
+# Test: ROSS instrumentation model-level stats, all sampling modes at once
+# (--model-stats=4 = GVT + real-time + virtual-time).
+# Checks the virtual-time (analysis LP) path produces a parseable
+# ross-stats-analysis-lps.bin alongside ross-stats-model.bin. Real-time records
+# are wall-clock dependent (a fast run may emit few or none), so their presence
+# is not asserted -- only that the run and the files stay well-formed.
+# Format spec: doc/model-stats-binary-format.md
+
+if [[ -z $bindir ]] ; then
+    echo bindir variable not set
+    exit 1
+fi
+
+export PACKET_SIZE=4096
+export CHUNK_SIZE=64
+export PACKET_LATENCY_TRACE_PATH=packet-latency-highdef/
+cat "$bindir/doc/example"/tutorial-ping-pong.template.conf.in | envsubst > tutorial-ping-pong.conf
+
+# NOTE: --vt-samp-end is required here: it defaults to g_tw_ts_end, which the
+# tutorial sets to 24h of virtual time, so without it the analysis LPs would
+# keep self-scheduling sample events long after the ping pong traffic ends.
+# 2e6 ns with 1e5 ns intervals covers the active phase with ~20 sample points.
+mpirun -np 3 "$bindir/doc/example/tutorial-synthetic-ping-pong" --sync=3 \
+    --num_messages=10 --payload_sz=8192 \
+    --model-stats=4 --num-gvt=100 --rt-interval=1 \
+    --vt-interval=1e5 --vt-samp-end=2e6 \
+    --stats-path=model-stats-out \
+    -- tutorial-ping-pong.conf \
+    > model-output.txt 2> model-output-error.txt
+err=$?
+[[ $err -ne 0 ]] && exit $err
+
+grep 'Net Events Processed' model-output.txt
+err=$?
+[[ $err -ne 0 ]] && exit $err
+
+# Both output files exist and are non-empty
+for f in ross-stats-model.bin ross-stats-analysis-lps.bin ; do
+    if [[ ! -s model-stats-out/$f ]] ; then
+        >&2 echo "model-stats-out/$f missing or empty"
+        exit 1
+    fi
+done
+
+if command -v python3 > /dev/null ; then
+    scriptdir="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+    python3 "$scriptdir/../scripts/parse-ross-model-stats.py" \
+        model-stats-out/ross-stats-model.bin \
+        model-stats-out/ross-stats-analysis-lps.bin \
+        --csv-prefix model-stats 2> parser-summary.txt
+    err=$?
+    if [[ $err -ne 0 ]] ; then
+        >&2 cat parser-summary.txt
+        >&2 echo "reference parser failed"
+        exit $err
+    fi
+    cat parser-summary.txt
+
+    grep 'unknown payloads: 0' parser-summary.txt
+    err=$?
+    if [[ $err -ne 0 ]] ; then
+        >&2 echo "parser hit unknown payload sizes (payload/format drift?)"
+        exit 1
+    fi
+
+    # the virtual-time path produced records (first CSV column is stats_type)
+    if ! cut -d, -f1 model-stats-svr.csv | grep -q '^vt$' ; then
+        >&2 echo "no virtual-time svr records decoded"
+        exit 1
+    fi
+else
+    echo "python3 not found; skipping reference-parser check"
+fi
diff --git a/tests/example-ping-pong-model-stats.sh b/tests/example-ping-pong-model-stats.sh
new file mode 100755
index 00000000..1be0d4af
--- /dev/null
+++ b/tests/example-ping-pong-model-stats.sh
@@ -0,0 +1,75 @@
+#!/bin/bash
+# Test: ROSS instrumentation model-level stats, GVT mode (--model-stats=1).
+# Runs the ping pong tutorial with GVT-based model-stats collection and checks
+# that ross-stats-model.bin is produced, non-empty, and parseable by the
+# reference parser (scripts/parse-ross-model-stats.py), with records for all
+# three LP payloads (svr, dally terminal, dally router).
+# Format spec: doc/model-stats-binary-format.md
+
+if [[ -z $bindir ]] ; then
+    echo bindir variable not set
+    exit 1
+fi
+
+export PACKET_SIZE=4096
+export CHUNK_SIZE=64
+export PACKET_LATENCY_TRACE_PATH=packet-latency-highdef/
+cat "$bindir/doc/example"/tutorial-ping-pong.template.conf.in | envsubst > tutorial-ping-pong.conf
+
+# GVT-based sampling needs a synchronization protocol that computes intermediate
+# GVTs, so run conservative/optimistic (sync=3), not sequential. This run computes
+# a few thousand GVTs; sampling every 100th keeps the output around a megabyte
+# (every GVT would be ~50x that) while still yielding dozens of samples per LP.
+mpirun -np 3 "$bindir/doc/example/tutorial-synthetic-ping-pong" --sync=3 \
+    --num_messages=10 --payload_sz=8192 \
+    --model-stats=1 --num-gvt=100 --stats-path=model-stats-out \
+    -- tutorial-ping-pong.conf \
+    > model-output.txt 2> model-output-error.txt
+err=$?
+[[ $err -ne 0 ]] && exit $err
+
+# Simulation completed
+grep 'Net Events Processed' model-output.txt
+err=$?
+[[ $err -ne 0 ]] && exit $err
+
+# Model stats file exists and is non-empty
+if [[ ! -s model-stats-out/ross-stats-model.bin ]] ; then
+    >&2 echo "model-stats-out/ross-stats-model.bin missing or empty"
+    exit 1
+fi
+
+# Parse it with the reference parser (soft-skip if python3 is unavailable;
+# python3 is not a hard test dependency of this repo)
+if command -v python3 > /dev/null ; then
+    scriptdir="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+    python3 "$scriptdir/../scripts/parse-ross-model-stats.py" \
+        model-stats-out/ross-stats-model.bin \
+        --csv-prefix model-stats 2> parser-summary.txt
+    err=$?
+    if [[ $err -ne 0 ]] ; then
+        >&2 cat parser-summary.txt
+        >&2 echo "reference parser failed on ross-stats-model.bin"
+        exit $err
+    fi
+    cat parser-summary.txt
+
+    # every payload size decoded (config dims match the parser defaults)
+    grep 'unknown payloads: 0' parser-summary.txt
+    err=$?
+    if [[ $err -ne 0 ]] ; then
+        >&2 echo "parser hit unknown payload sizes (payload/format drift?)"
+        exit 1
+    fi
+
+    # at least one record of each LP type (csv has a header line)
+    for lptype in svr terminal router ; do
+        nrec=$(( $(wc -l < model-stats-${lptype}.csv) - 1 ))
+        if [[ $nrec -lt 1 ]] ; then
+            >&2 echo "no ${lptype} records in ross-stats-model.bin"
+            exit 1
+        fi
+    done
+else
+    echo "python3 not found; skipping reference-parser check"
+fi