From a8172052eb8df15653cb697e8cba6034dcb7c58e Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Thu, 6 Aug 2026 12:53:43 -0400 Subject: [PATCH 1/9] Add Philox4x32-10 counter-based RNG core Add a deterministic, keyed random generator, Random(cell_id, time_step, purpose, sub_index), built on Philox4x32-10, plus a thread-local deterministic context that UniformRandom(), UniformInt(), NormalRandom(), and LogNormalRandom() consult when active. Existing call sites keep working unchanged, but draw from a reproducible sequence instead of a thread-order-dependent one whenever the context is installed. --- core/PhysiCell_utilities.cpp | 143 +++++++++++++++++++++++++++++++++++ core/PhysiCell_utilities.h | 12 +++ modules/philox.h | 103 +++++++++++++++++++++++++ 3 files changed, 258 insertions(+) create mode 100644 modules/philox.h diff --git a/core/PhysiCell_utilities.cpp b/core/PhysiCell_utilities.cpp index 91f0cd582..5e59983af 100644 --- a/core/PhysiCell_utilities.cpp +++ b/core/PhysiCell_utilities.cpp @@ -69,12 +69,107 @@ #include "PhysiCell_constants.h" #include "PhysiCell.h" +#include "../modules/philox.h" + +#include +#include +#include namespace PhysiCell{ thread_local std::mt19937_64 physicell_PRNG_generator; thread_local bool local_pnrg_setup_done = false; +namespace +{ + +struct DeterministicRandomContext +{ + bool active = false; + std::uint64_t cell_id = 0; + std::uint64_t time_step = 0; + std::uint64_t purpose = 0; + std::uint64_t sub_index = 0; +}; + +thread_local DeterministicRandomContext deterministic_random_context; + +std::array make_counter( std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ) +{ + return ::philox::make_counter( cell_id, time_step, purpose, sub_index ); +} + +std::array make_key( void ) +{ + std::uint64_t seed0 = ::philox::splitmix64( static_cast( physicell_random_seed ) ); + std::uint64_t seed1 = ::philox::splitmix64( seed0 ^ 0xD1342543DE82EF95ull ); + return { + static_cast( seed0 ), + static_cast( seed1 ) + }; +} + +std::array make_philox_block( std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ) +{ + return ::philox::philox4x32_10( make_counter( cell_id, time_step, purpose, sub_index ), make_key() ); +} + +double block_to_unit_double( const std::array& block ) +{ + std::uint64_t mantissa = ( static_cast( block[0] ) << 21 ) | ( static_cast( block[1] ) >> 11 ); + return std::ldexp( static_cast( mantissa ), -53 ); +} + +double next_context_random( void ) +{ + double value = Random( deterministic_random_context.cell_id, + deterministic_random_context.time_step, + deterministic_random_context.purpose, + deterministic_random_context.sub_index ); + deterministic_random_context.sub_index += 1; + return value; +} + +double next_context_normal( double mean, double standard_deviation ) +{ + double u1 = next_context_random(); + double u2 = next_context_random(); + u1 = std::max( u1, std::numeric_limits::min() ); + double magnitude = std::sqrt( -2.0 * std::log( u1 ) ); + double angle = 6.283185307179586476925286766559 * u2; + return mean + standard_deviation * magnitude * std::cos( angle ); +} + +} // namespace + +void set_deterministic_random_context( std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose ) +{ + deterministic_random_context.active = true; + deterministic_random_context.cell_id = cell_id; + deterministic_random_context.time_step = time_step; + deterministic_random_context.purpose = purpose; + deterministic_random_context.sub_index = 0; +} + +void clear_deterministic_random_context( void ) +{ + deterministic_random_context.active = false; + deterministic_random_context.cell_id = 0; + deterministic_random_context.time_step = 0; + deterministic_random_context.purpose = 0; + deterministic_random_context.sub_index = 0; +} + +double Random( std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ) +{ + return block_to_unit_double( make_philox_block( cell_id, time_step, purpose, sub_index ) ); +} + +double Random( void ) +{ + return UniformRandom(); +} + unsigned int physicell_random_seed = 0; std::vector physicell_random_seeds; @@ -147,6 +242,11 @@ double UniformRandom_old_not_thread_safe() double UniformRandom( void ) { + if( deterministic_random_context.active ) + { + return next_context_random(); + } + thread_local std::uniform_real_distribution distribution(0.0,1.0); if( local_pnrg_setup_done == false ) { @@ -175,25 +275,68 @@ double UniformRandom( void ) */ } +double UniformRandom( std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ) +{ + return Random( cell_id, time_step, purpose, sub_index ); +} + int UniformInt() { + if( deterministic_random_context.active ) + { + return static_cast( make_philox_block( deterministic_random_context.cell_id, + deterministic_random_context.time_step, + deterministic_random_context.purpose, + deterministic_random_context.sub_index++ )[0] & 0x7fffffffU ); + } + static std::uniform_int_distribution int_dis; return int_dis(physicell_PRNG_generator); } +int UniformInt( std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ) +{ + return static_cast( make_philox_block( cell_id, time_step, purpose, sub_index )[0] & 0x7fffffffU ); +} + double NormalRandom( double mean, double standard_deviation ) { + if( deterministic_random_context.active ) + { + return next_context_normal( mean, standard_deviation ); + } + std::normal_distribution d(mean,standard_deviation); return d(physicell_PRNG_generator); } +double NormalRandom( double mean, double standard_deviation, std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ) +{ + double u1 = Random( cell_id, time_step, purpose, sub_index ); + double u2 = Random( cell_id, time_step, purpose, sub_index + 1 ); + u1 = std::max( u1, std::numeric_limits::min() ); + double magnitude = std::sqrt( -2.0 * std::log( u1 ) ); + double angle = 6.283185307179586476925286766559 * u2; + return mean + standard_deviation * magnitude * std::cos( angle ); +} + double LogNormalRandom( double mean, double standard_deviation ) { + if( deterministic_random_context.active ) + { + return exp( next_context_normal( log( mean ), standard_deviation ) ); + } + return exp(NormalRandom(log(mean), standard_deviation)); } +double LogNormalRandom( double mean, double standard_deviation, std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ) +{ + return exp( NormalRandom( log( mean ), standard_deviation, cell_id, time_step, purpose, sub_index ) ); +} + std::vector UniformOnUnitSphere( void ) { std::vector output = {0,0,0}; diff --git a/core/PhysiCell_utilities.h b/core/PhysiCell_utilities.h index 7e4ec1548..0d9bbfd14 100644 --- a/core/PhysiCell_utilities.h +++ b/core/PhysiCell_utilities.h @@ -75,6 +75,7 @@ #include #include #include +#include #include @@ -82,16 +83,27 @@ namespace PhysiCell{ extern std::vector physicell_random_seeds; + extern unsigned int physicell_random_seed; + + void set_deterministic_random_context( std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose ); + void clear_deterministic_random_context( void ); + + double Random( std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ); + double Random( void ); void setup_rng( void ); void SeedRandom( unsigned int input ); void SeedRandom( void ); double UniformRandom( void ); + double UniformRandom( std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ); int UniformInt( void ); + int UniformInt( std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ); double NormalRandom( double mean, double standard_deviation ); + double NormalRandom( double mean, double standard_deviation, std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ); double LogNormalRandom( double mean, double standard_deviation ); + double LogNormalRandom( double mean, double standard_deviation, std::uint64_t cell_id, std::uint64_t time_step, std::uint64_t purpose, std::uint64_t sub_index ); std::vector UniformOnUnitSphere( void ); std::vector UniformOnUnitCircle( void ); diff --git a/modules/philox.h b/modules/philox.h new file mode 100644 index 000000000..7527ceb92 --- /dev/null +++ b/modules/philox.h @@ -0,0 +1,103 @@ +/* +Copyright 2010-2011, D. E. Shaw Research. +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +* Redistributions of source code must retain the above copyright + notice, this list of conditions, and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright + notice, this list of conditions, and the following disclaimer in the + documentation and/or other materials provided with the distribution. + +* Neither the name of D. E. Shaw Research nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#ifndef _philox_dot_h_ +#define _philox_dot_h_ + +/* Vendored from Random123/philox.h with the original copyright and license preserved. */ + +#include +#include +#include + +namespace philox { + +inline std::uint64_t splitmix64(std::uint64_t value) +{ + value += 0x9E3779B97F4A7C15ull; + value = (value ^ (value >> 30)) * 0xBF58476D1CE4E5B9ull; + value = (value ^ (value >> 27)) * 0x94D049BB133111EBull; + return value ^ (value >> 31); +} + +inline std::pair mulhilo(std::uint32_t lhs, std::uint32_t rhs) +{ + std::uint64_t product = static_cast(lhs) * static_cast(rhs); + return {static_cast(product >> 32), static_cast(product)}; +} + +inline std::array philox4x32_10(std::array counter, + std::array key) +{ + constexpr std::uint32_t PHILOX_M0 = 0xD2511F53u; + constexpr std::uint32_t PHILOX_M1 = 0xCD9E8D57u; + constexpr std::uint32_t PHILOX_W0 = 0x9E3779B9u; + constexpr std::uint32_t PHILOX_W1 = 0xBB67AE85u; + + for (int round = 0; round < 10; ++round) + { + auto lo_hi0 = mulhilo(PHILOX_M0, counter[0]); + auto lo_hi1 = mulhilo(PHILOX_M1, counter[2]); + + std::array next_counter = { + lo_hi1.first ^ counter[1] ^ key[0], + lo_hi1.second, + lo_hi0.first ^ counter[3] ^ key[1], + lo_hi0.second + }; + + counter = next_counter; + key[0] += PHILOX_W0; + key[1] += PHILOX_W1; + } + + return counter; +} + +inline std::array make_counter(std::uint64_t cell_id, + std::uint64_t time_step, + std::uint64_t purpose, + std::uint64_t sub_index) +{ + std::uint64_t mix0 = splitmix64(cell_id ^ (purpose << 1)); + std::uint64_t mix1 = splitmix64(time_step ^ (sub_index << 1) ^ (purpose << 17)); + return { + static_cast(mix0), + static_cast(mix0 >> 32), + static_cast(mix1), + static_cast(mix1 >> 32) + }; +} + +} // namespace philox + +#endif /* _philox_dot_h_ */ From b6acd97c0350b86b23db9cd16f15e879ce8237c7 Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Thu, 6 Aug 2026 12:56:26 -0400 Subject: [PATCH 2/9] Wire counter-based RNG into the simulation loop; fix cross-cell races Add PhysiCell_settings.use_counter_based_rng (default off) and the config element (counter_based/philox/counter to enable; legacy/thread_local/mt19937 to disable explicitly). Install a per-cell deterministic RNG context, keyed by cell ID/step/phase, around each parallel per-cell update phase, and sort cells_ready_to_divide / cells_ready_to_die by cell ID before applying them serially, so new cell IDs and (*all_cells) ordering are thread-count independent. RNG determinism alone doesn't stop one cell's thread from mutating another cell's live state mid-update in the same parallel loop. Three such races are fixed here with #pragma omp ordered around the per-cell call, forcing the same relative commit order as a single-thread run: - standard_cell_cell_interactions() (attack/ingest/fuse) mutates a target cell with no ordering; two attackers can race on one target. - dynamic_spring_attachments() checks a neighbor's attachment capacity before attach_cells_as_spring() re-validates it, so two cells can both attach to an already-full neighbor. - Secretion::advance() -> simulate_secretion_and_uptake() accumulates into a shared per-voxel density vector with no lock; two cells sharing a voxel on different threads can lose an update. --- core/PhysiCell_cell_container.cpp | 155 ++++++++++++++++++++++++++---- core/PhysiCell_cell_container.h | 1 + modules/PhysiCell_settings.cpp | 21 ++++ modules/PhysiCell_settings.h | 1 + 4 files changed, 159 insertions(+), 19 deletions(-) diff --git a/core/PhysiCell_cell_container.cpp b/core/PhysiCell_cell_container.cpp index 0c5850b1d..97262b04c 100644 --- a/core/PhysiCell_cell_container.cpp +++ b/core/PhysiCell_cell_container.cpp @@ -69,9 +69,12 @@ #include "PhysiCell_constants.h" #include "../BioFVM/BioFVM_vector.h" #include "PhysiCell_cell.h" +#include "PhysiCell_utilities.h" #include +#include #include +#include using namespace BioFVM; @@ -79,6 +82,22 @@ namespace PhysiCell{ std::vector *all_cells; +namespace +{ +constexpr std::uint64_t RANDOM_PURPOSE_SECRETION = 1; +constexpr std::uint64_t RANDOM_PURPOSE_INTRACELLULAR = 2; +constexpr std::uint64_t RANDOM_PURPOSE_PHENOTYPE = 3; +constexpr std::uint64_t RANDOM_PURPOSE_MECHANICS = 4; +constexpr std::uint64_t RANDOM_PURPOSE_INTERACTION = 5; +constexpr std::uint64_t RANDOM_PURPOSE_CUSTOM_RULE = 6; +constexpr std::uint64_t RANDOM_PURPOSE_UPDATE_VELOCITY = 7; +constexpr std::uint64_t RANDOM_PURPOSE_SPRINGS = 8; +constexpr std::uint64_t RANDOM_PURPOSE_CELL_CELL = 9; +constexpr std::uint64_t RANDOM_PURPOSE_POSITION = 10; +constexpr std::uint64_t RANDOM_PURPOSE_DIVISION = 11; +constexpr std::uint64_t RANDOM_PURPOSE_DEATH = 12; +} + Cell_Container::Cell_Container() { all_cells = (std::vector *) &all_basic_agents; @@ -122,14 +141,44 @@ void Cell_Container::update_all_cells(double t) void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double mechanics_dt_ , double diffusion_dt_ ) { - // secretions and uptakes. Syncing with BioFVM is automated. - - #pragma omp parallel for + std::uint64_t random_step = static_cast( std::llround( t / diffusion_dt_ ) ); + bool use_counter_based_rng = PhysiCell_settings.use_counter_based_rng; + auto activate_random_context = [use_counter_based_rng, random_step]( std::uint64_t cell_id, std::uint64_t purpose ) + { + if( use_counter_based_rng ) + { + set_deterministic_random_context( cell_id, random_step, purpose ); + } + }; + auto clear_random_context = [use_counter_based_rng]() + { + if( use_counter_based_rng ) + { + clear_deterministic_random_context(); + } + }; + + // secretions and uptakes. Syncing with BioFVM is automated. + // ordered: Secretion::advance() -> Basic_Agent::simulate_secretion_and_uptake() accumulates directly + // into the shared voxel density vector at (*pS)(current_voxel_index), with no lock (see + // BioFVM/BioFVM_basic_agent.cpp). Diffusion voxels are typically larger than a single cell, so + // multiple cells commonly share one; without a fixed commit order, two cells in the same voxel + // processed by different threads at the same time race on that shared read-modify-write (a genuine + // lost-update, not just reordering). This is the actual call path PhysiCell uses for cell secretion + // (unlike Microenvironment::simulate_cell_sources_and_sinks, which is never invoked by the core + // simulation loop) -- forcing the same relative commit order as a single-thread run removes the race. + + #pragma omp parallel for schedule(static) ordered for( int i=0; i < (*all_cells).size(); i++ ) { if( (*all_cells)[i]->is_out_of_domain == false ) { - (*all_cells)[i]->phenotype.secretion.advance( (*all_cells)[i], (*all_cells)[i]->phenotype , diffusion_dt_ ); + activate_random_context( (*all_cells)[i]->ID, RANDOM_PURPOSE_SECRETION ); + #pragma omp ordered + { + (*all_cells)[i]->phenotype.secretion.advance( (*all_cells)[i], (*all_cells)[i]->phenotype , diffusion_dt_ ); + } + clear_random_context(); } } @@ -148,6 +197,7 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me if( (*all_cells)[i]->phenotype.intracellular != NULL && (*all_cells)[i]->phenotype.intracellular->need_update()) { + activate_random_context( (*all_cells)[i]->ID, RANDOM_PURPOSE_INTRACELLULAR ); if ((*all_cells)[i]->functions.pre_update_intracellular != NULL) (*all_cells)[i]->functions.pre_update_intracellular( (*all_cells)[i], (*all_cells)[i]->phenotype , diffusion_dt_ ); @@ -155,6 +205,7 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me if ((*all_cells)[i]->functions.post_update_intracellular != NULL) (*all_cells)[i]->functions.post_update_intracellular( (*all_cells)[i], (*all_cells)[i]->phenotype , diffusion_dt_ ); + clear_random_context(); } } } @@ -177,18 +228,38 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me { if( (*all_cells)[i]->is_out_of_domain == false ) { + activate_random_context( (*all_cells)[i]->ID, RANDOM_PURPOSE_PHENOTYPE ); (*all_cells)[i]->advance_bundled_phenotype_functions( time_since_last_cycle ); + clear_random_context(); } } // process divides / removes + // sort cells_ready_to_divide by cell ID to make sure that new generation of cell ids are preserved. + if ( PhysiCell_settings.use_counter_based_rng ) + std::sort( cells_ready_to_divide.begin(), cells_ready_to_divide.end(), + []( const Cell* lhs, const Cell* rhs ) + { + return lhs->ID < rhs->ID; + } ); for( int i=0; i < cells_ready_to_divide.size(); i++ ) { + activate_random_context( cells_ready_to_divide[i]->ID, RANDOM_PURPOSE_DIVISION ); cells_ready_to_divide[i]->divide(); + clear_random_context(); } + // sort cells_ready_to_die by cell ID to make sure that all_cells order is preserved. + if ( PhysiCell_settings.use_counter_based_rng ) + std::sort( cells_ready_to_die.begin(), cells_ready_to_die.end(), + []( const Cell* lhs, const Cell* rhs ) + { + return lhs->ID < rhs->ID; + } ); for( int i=0; i < cells_ready_to_die.size(); i++ ) { + activate_random_context( cells_ready_to_die[i]->ID, RANDOM_PURPOSE_DEATH ); cells_ready_to_die[i]->die(); + clear_random_context(); } num_divisions_in_current_step+= cells_ready_to_divide.size(); num_deaths_in_current_step+= cells_ready_to_die.size(); @@ -220,7 +291,11 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me { Cell* pC = (*all_cells)[i]; if( pC->functions.contact_function && pC->is_out_of_domain == false ) - { evaluate_interactions( pC,pC->phenotype,time_since_last_mechanics ); } + { + activate_random_context( pC->ID, RANDOM_PURPOSE_INTERACTION ); + evaluate_interactions( pC,pC->phenotype,time_since_last_mechanics ); + clear_random_context(); + } } // perform custom computations @@ -231,7 +306,11 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me Cell* pC = (*all_cells)[i]; if( pC->functions.custom_cell_rule && pC->is_out_of_domain == false ) - { pC->functions.custom_cell_rule( pC,pC->phenotype,time_since_last_mechanics ); } + { + activate_random_context( pC->ID, RANDOM_PURPOSE_CUSTOM_RULE ); + pC->functions.custom_cell_rule( pC,pC->phenotype,time_since_last_mechanics ); + clear_random_context(); + } } // update velocities @@ -241,7 +320,11 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me { Cell* pC = (*all_cells)[i]; if( pC->functions.update_velocity && pC->is_out_of_domain == false && pC->is_movable ) - { pC->functions.update_velocity( pC,pC->phenotype,time_since_last_mechanics ); } + { + activate_random_context( pC->ID, RANDOM_PURPOSE_UPDATE_VELOCITY ); + pC->functions.update_velocity( pC,pC->phenotype,time_since_last_mechanics ); + clear_random_context(); + } } // new March 2023: @@ -249,41 +332,71 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me if( PhysiCell_settings.disable_automated_spring_adhesions == false ) { - #pragma omp parallel for + // ordered: dynamic_spring_attachments() decides whether to attach/detach based on a + // *neighbor's* current attachment count vs. its max, then calls attach_cells_as_spring()/ + // detach_cells_as_spring() on that neighbor. The capacity check is never re-validated once + // the lock inside attach_cells_as_spring() is actually held, so without a fixed commit order + // two different cells can both see "room for one more" on the same neighbor at once and both + // attach, pushing it past its intended maximum -- a thread-count-dependent result. + #pragma omp parallel for schedule(static) ordered for( int i=0; i < (*all_cells).size(); i++ ) { - Cell* pC = (*all_cells)[i]; - dynamic_spring_attachments(pC,pC->phenotype,time_since_last_mechanics); - } - #pragma omp parallel for + Cell* pC = (*all_cells)[i]; + activate_random_context( pC->ID, RANDOM_PURPOSE_SPRINGS ); + #pragma omp ordered + { + dynamic_spring_attachments(pC,pC->phenotype,time_since_last_mechanics); + } + clear_random_context(); + } + #pragma omp parallel for for( int i=0; i < (*all_cells).size(); i++ ) { - Cell* pC = (*all_cells)[i]; + Cell* pC = (*all_cells)[i]; if( pC->is_movable ) { for( int j=0; j < pC->state.spring_attachments.size(); j++ ) { Cell* pC1 = pC->state.spring_attachments[j]; // standard_elastic_contact_function_confluent_rest_length(pC,pC->phenotype,pC1,pC1->phenotype,time_since_last_mechanics); + activate_random_context( pC->ID, RANDOM_PURPOSE_SPRINGS ); standard_elastic_contact_function(pC,pC->phenotype,pC1,pC1->phenotype,time_since_last_mechanics); + clear_random_context(); } } } } - // new March 2022: - // run standard interactions (phagocytosis, attack, fusion) here - #pragma omp parallel for + // new March 2022: + // run standard interactions (phagocytosis, attack, fusion) here + // this loop is "ordered": standard_cell_cell_interactions() can ingest, attack, or fuse + // a *different* cell (mutating its live state), and it early-exits based on that cell's + // current death flag. Enforcing the same relative commit order as a single-thread run + // (via #pragma omp ordered) makes those cross-cell interactions thread-count independent, + // instead of depending on however OpenMP happens to schedule the threads this run. + #pragma omp parallel for schedule(static) ordered for( int i=0; i < (*all_cells).size(); i++ ) { - Cell* pC = (*all_cells)[i]; - standard_cell_cell_interactions(pC,pC->phenotype,time_since_last_mechanics); + Cell* pC = (*all_cells)[i]; + activate_random_context( pC->ID, RANDOM_PURPOSE_CELL_CELL ); + #pragma omp ordered + { + standard_cell_cell_interactions(pC,pC->phenotype,time_since_last_mechanics); + } + clear_random_context(); } // super-critical to performance! clear the "dummy" cells from phagocytosis / fusion // otherwise, comptuational cost increases at polynomial rate VERY fast, as O(10,000) // dummy cells of size zero are left ot interact mechanically, etc. if( cells_ready_to_die.size() > 0 ) { + // sort cells_ready_to_die by cell ID to make sure that all_cells order is preserved. + if ( PhysiCell_settings.use_counter_based_rng ) + std::sort( cells_ready_to_die.begin(), cells_ready_to_die.end(), + []( const Cell* lhs, const Cell* rhs ) + { + return lhs->ID < rhs->ID; + } ); /* std::cout << "\tClearing dummy cells from phagocytosis and fusion events ... " << std::endl; std::cout << "\t\tClearing " << cells_ready_to_die.size() << " cells ... " << std::endl; @@ -302,7 +415,11 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me { Cell* pC = (*all_cells)[i]; if( pC->is_out_of_domain == false && pC->is_movable) - { pC->update_position(time_since_last_mechanics); } + { + activate_random_context( pC->ID, RANDOM_PURPOSE_POSITION ); + pC->update_position(time_since_last_mechanics); + clear_random_context(); + } } // When somebody reviews this code, let's add proper braces for clarity!!! diff --git a/core/PhysiCell_cell_container.h b/core/PhysiCell_cell_container.h index f0b2cb9ee..2ea3964dd 100644 --- a/core/PhysiCell_cell_container.h +++ b/core/PhysiCell_cell_container.h @@ -68,6 +68,7 @@ #ifndef __PhysiCell_cell_container_h__ #define __PhysiCell_cell_container_h__ +#include #include #include "PhysiCell_cell.h" #include "../BioFVM/BioFVM_agent_container.h" diff --git a/modules/PhysiCell_settings.cpp b/modules/PhysiCell_settings.cpp index a588383d1..0163fd148 100644 --- a/modules/PhysiCell_settings.cpp +++ b/modules/PhysiCell_settings.cpp @@ -276,6 +276,27 @@ void PhysiCell_Settings::read_from_pugixml( void ) PhysiCell_settings.disable_automated_spring_adhesions = true; } + pugi::xml_node rng_mode_node = xml_find_node(node_options, "rng_mode"); + if (rng_mode_node) + { + std::string rng_mode = xml_get_my_string_value(rng_mode_node); + if( rng_mode == "counter_based" || rng_mode == "counter" || rng_mode == "philox" ) + { + PhysiCell_settings.use_counter_based_rng = true; + std::cout << "Using counter-based RNG mode" << std::endl; + } + else if( rng_mode == "legacy" || rng_mode == "thread_local" || rng_mode == "mt19937" ) + { + PhysiCell_settings.use_counter_based_rng = false; + std::cout << "Using legacy RNG mode" << std::endl; + } + else if( rng_mode != "" ) + { + std::cout << "ERROR: unsupported rng_mode '" << rng_mode << "'. Use 'legacy' or 'counter_based'." << std::endl; + exit(-1); + } + } + pugi::xml_node random_seed_node = xml_find_node(node_options, "random_seed"); std::string random_seed = ""; // default is system clock, even if this element is not present if (random_seed_node) diff --git a/modules/PhysiCell_settings.h b/modules/PhysiCell_settings.h index 86afb2939..3d5a69d10 100644 --- a/modules/PhysiCell_settings.h +++ b/modules/PhysiCell_settings.h @@ -114,6 +114,7 @@ class PhysiCell_Settings bool enable_legacy_saves = false; bool disable_automated_spring_adhesions = false; + bool use_counter_based_rng = false; double SVG_save_interval = 60; bool enable_SVG_saves = true; From 7b503c1dae3b51b1f4f9e0abf9dcfaa7f10ef919 Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Thu, 6 Aug 2026 12:56:55 -0400 Subject: [PATCH 3/9] Add thread reproducibility test script Runs a built project twice at two thread counts with the same config and seed, then diffs the final snapshot output to catch thread-count-dependent divergence. --- beta/test_thread_repro.py | 118 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 118 insertions(+) create mode 100644 beta/test_thread_repro.py diff --git a/beta/test_thread_repro.py b/beta/test_thread_repro.py new file mode 100644 index 000000000..87f8a63f9 --- /dev/null +++ b/beta/test_thread_repro.py @@ -0,0 +1,118 @@ +import shutil +import subprocess +import sys +import tempfile +from datetime import datetime +import xml.etree.ElementTree as ET +from pathlib import Path + + +def update_config(xml_file, output_dir, threads, max_time, rng_mode): + tree = ET.parse(xml_file) + root = tree.getroot() + + root.find(".//overall/max_time").text = str(max_time) + root.find(".//parallel/omp_num_threads").text = str(threads) + + save_folder = root.find(".//save/folder") + if save_folder is None: + save_node = root.find(".//save") + if save_node is None: + raise RuntimeError("Could not find in config file") + save_folder = ET.SubElement(save_node, "folder") + save_folder.text = str(output_dir) + + random_seed = root.find(".//options/random_seed") + if random_seed is None: + random_seed = root.find(".//user_parameters/random_seed") + if random_seed is not None: + random_seed.text = "17" + + options = root.find(".//options") + if options is None: + options = ET.SubElement(root, "options") + + rng_mode_node = options.find("rng_mode") + if rng_mode is None: + if rng_mode_node is not None: + options.remove(rng_mode_node) + else: + if rng_mode_node is None: + rng_mode_node = ET.SubElement(options, "rng_mode") + rng_mode_node.text = rng_mode + + tree.write(xml_file) + + +def run_once(repo_root, executable, config_file): + cmd = [str(repo_root / executable), str(config_file)] + subprocess.run(cmd, cwd=repo_root, check=True) + + +def compare_svg_outputs(repo_root, output_dir_a, output_dir_b): + checker = repo_root / "beta" / "test_diff_svg.py" + cmd = [sys.executable, str(checker), str(output_dir_a), str(output_dir_b)] + subprocess.run(cmd, cwd=repo_root, check=True) + + +def main(executable, config_file, max_time, threads_a, threads_b, rng_mode="counter_based", work_dir=None): + repo_root = Path.cwd() + config_source = repo_root / config_file + max_time = int(max_time) + threads_a = int(threads_a) + threads_b = int(threads_b) + rng_mode_arg = rng_mode.strip().lower() if isinstance(rng_mode, str) else rng_mode + if rng_mode_arg in ("", "off", "none", "omit", "omitted", "default"): + rng_mode = None + if threads_a == threads_b: + raise SystemExit("threads_a and threads_b must be different") + + if work_dir is None: + tmp_context = tempfile.TemporaryDirectory(prefix="physicell-thread-repro-") + tmp_root = Path(tmp_context.__enter__()) + cleanup_context = tmp_context + else: + work_root = Path(work_dir) + work_root.mkdir(parents=True, exist_ok=True) + tmp_root = work_root / f"physicell-thread-repro_{datetime.now().strftime('%Y%m%d_%H%M%S')}" + tmp_root.mkdir(parents=True, exist_ok=True) + cleanup_context = None + + try: + configs = [] + outputs = [] + + for run_index, threads in enumerate((threads_a, threads_b), start=1): + run_dir = tmp_root / f"run_{run_index}_threads_{threads}" + output_dir = run_dir / f"output_threads_{threads}" + output_dir.mkdir(parents=True, exist_ok=True) + + config_copy = run_dir / "PhysiCell_settings.xml" + shutil.copy(config_source, config_copy) + update_config(config_copy, output_dir, threads, max_time, rng_mode) + + print( + f"\n\n {'=' * 80}\n running {executable} with {threads} threads " + f"(rng_mode={rng_mode if rng_mode is not None else 'omitted'})" + ) + print(f" config: {config_copy}") + print(f" output: {output_dir}") + run_once(repo_root, executable, config_copy) + + configs.append(config_copy) + outputs.append(output_dir / "final.xml") + + for output in outputs: + if not output.exists(): + raise SystemExit(f"Expected output file not found: {output}") + + compare_svg_outputs(repo_root, outputs[0].parent, outputs[1].parent) + + print(f"\n\n {'=' * 80}\n thread reproducibility test passed") + finally: + if cleanup_context is not None: + cleanup_context.__exit__(None, None, None) + + +if __name__ == "__main__": + main(*sys.argv[1:]) \ No newline at end of file From cb24e64cdb5c907a090cca99a6d9f4e44f40054e Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Thu, 6 Aug 2026 12:57:59 -0400 Subject: [PATCH 4/9] Document counter-based RNG design and thread-safety fixes Explain the Philox-based deterministic RNG scheme, the existing division/death ordering precedent, the four cross-cell mutation races found and fixed, how to enable counter-based RNG, and how to test a new/modified project for thread reproducibility with beta/test_thread_repro.py. --- counter_based_rng.md | 318 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 318 insertions(+) create mode 100644 counter_based_rng.md diff --git a/counter_based_rng.md b/counter_based_rng.md new file mode 100644 index 000000000..18282aa9a --- /dev/null +++ b/counter_based_rng.md @@ -0,0 +1,318 @@ +# Counter-Based RNG Changes + +This branch adds a deterministic RNG path based on Philox4x32-10 so the same simulation can produce the same random values even when the OpenMP thread count changes. + +## What changed + +The core RNG layer now has a keyed entry point in [core/PhysiCell_utilities.h](core/PhysiCell_utilities.h) and [core/PhysiCell_utilities.cpp](core/PhysiCell_utilities.cpp). The new `Random(cell_id, time_step, purpose, sub_index)` API generates a value directly from the key and counter inputs, and the existing `UniformRandom()`, `UniformInt()`, `NormalRandom()`, and `LogNormalRandom()` functions can read from a thread-local deterministic context when that context is active. + +The main simulation loop sets that context around the major per-cell update phases in [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp). That means the existing simulation code can keep calling the usual random helpers, but those helpers will now draw from a reproducible per-cell, per-step sequence instead of a thread-order-dependent stream. + +I also added a reusable end-to-end thread reproducibility checker in [beta/test_thread_repro.py](beta/test_thread_repro.py). The GitHub Actions workflow uses it for the template and intracellular template jobs in [tests.yml](.github/workflows/tests.yml). + +The checker runs the same project twice with different thread counts and compares the final XML output after normalizing volatile metadata. + +## Why this did not require changing every `Random()` call + +The key idea is that the call sites did not need to change one by one because the RNG behavior was redirected at the shared utility layer. + +Before this change, every stochastic helper depended on the current thread-local generator state, so the result depended on execution order. After the change, the simulation sets a deterministic context once per cell and per update phase, and the helper functions consult that context automatically. + +In practice: + +- Code such as `UniformRandom()` still works unchanged. +- When a deterministic context is active, `UniformRandom()` uses the Philox-backed keyed path instead of the legacy stream. +- The same per-cell/per-step/key inputs always produce the same bit pattern, regardless of how OpenMP schedules the work. + +So the important change was not rewriting every random call site. The important change was making the shared RNG helpers context-aware and ensuring the simulation engine installs the right context before calling code that already uses those helpers. + +## When you still would change a call site + +You would only need to change an individual call site if it should use a different key than the one implied by the current simulation phase, or if it runs outside the engine-managed per-cell context and still needs deterministic behavior. + +That is why the initialization code in sample projects can still use the legacy `UniformRandom()` calls without issue for this test: the test focuses on the simulation loop, where the context is already installed. + +## Minimal key scheme for restart + +If the long-term goal is bitwise restart from a saved snapshot, the minimal keyed identity should still represent four different roles, even if not all four are stored explicitly as separate variables. + +- `cell_id`: which cell owns the draw. +- `step_key`: which saved simulation step or resumed update step owns the draw. +- `purpose`: which update phase or stochastic subsystem owns the draw. +- `sub_index`: which draw number within that cell/step/purpose context this is. + +That means the conceptual minimum is: + +`Random(cell_id, step_key, purpose, sub_index)` + +For restart, `step_key` does not have to be an unsaved in-memory counter. It can be reconstructed from restart-stable state, as long as the reconstruction is exact. For example, it can be derived from: + +- saved simulation time, +- saved diffusion / mechanics / phenotype scheduling state, +- and the current update phase. + +The important requirement is not the name of the field. The requirement is that every random draw after restart resolves to the same unique identity it would have had in the uninterrupted run. + +## Why each index matters + +### `cell_id` + +This separates one cell from another. + +Example: + +- cell 10 draws a motility angle, +- cell 11 draws a motility angle, +- both happen in the same step and same purpose. + +Without `cell_id`, both cells would ask for the same keyed random value. + +### `step_key` + +This separates one update step from the next. + +Example: + +- cell 10 updates motility at step 100, +- cell 10 updates motility again at step 101. + +Without `step_key`, the cell would receive the same random value every time that code path runs. + +### `purpose` + +This separates different stochastic uses within the same cell and same step. + +Example: + +- cell 10 uses randomness during phenotype update, +- cell 10 also uses randomness during velocity update, +- both happen in the same global simulation step. + +Without `purpose`, those two unrelated code paths would be drawing from the same keyed location and could collide. + +### `sub_index` + +This separates repeated random draws within the same cell, same step, and same purpose. + +Example: + +- cell 10 calls `UniformRandom()` twice during motility update, +- first draw selects whether to reorient, +- second draw selects an angle. + +Without `sub_index`, both calls would map to the same keyed random value. + +### Concrete multi-draw examples already in PhysiCell + +This is not just a theoretical issue. Several existing helpers already consume more than one uniform draw for one higher-level stochastic operation. + +- `NormalRandom()` uses two uniforms internally through the Box-Muller transform. See [core/PhysiCell_utilities.cpp](core/PhysiCell_utilities.cpp). +- `LogNormalRandom()` also uses that same two-uniform path because it is built on top of `NormalRandom()`. See [core/PhysiCell_utilities.cpp](core/PhysiCell_utilities.cpp). +- `UniformOnUnitSphere()` uses two uniforms, one for `z` and one for `theta`. See [core/PhysiCell_utilities.cpp](core/PhysiCell_utilities.cpp). +- `LegacyRandomOnUnitSphere()` uses multiple uniforms as well. See [core/PhysiCell_utilities.cpp](core/PhysiCell_utilities.cpp). +- `UniformInUnitDisc()` uses two uniforms. See [core/PhysiCell_rules.cpp](core/PhysiCell_rules.cpp). +- `UniformInUnitSphere()` uses three uniforms. See [core/PhysiCell_rules.cpp](core/PhysiCell_rules.cpp). +- `UniformInAnnulus()` uses two uniforms. See [core/PhysiCell_rules.cpp](core/PhysiCell_rules.cpp). +- `UniformInShell()` uses three uniforms. See [core/PhysiCell_rules.cpp](core/PhysiCell_rules.cpp). + +There are also mixed code paths where one logical operation first makes a decision draw and then consumes additional random values inside the same higher-level update. For example, motility update can first draw whether the cell reorients and then, if it does, draw a random direction on the circle or sphere. See [core/PhysiCell_cell.cpp](core/PhysiCell_cell.cpp). + +This is why `sub_index` is convenient even when `purpose` is automatic. A single semantic purpose such as `motility update` can safely contain several actual draws without inventing a separate micro-purpose for every internal step. + +## Why `cell_id` and `time_step` alone are not enough + +They are not enough if more than one random draw can happen for the same cell in the same step, which is normal in PhysiCell. + +There are two independent problems. + +### Problem 1: multiple draws within one context + +Suppose cell 10 is in phenotype step 25 and a function does this: + +```cpp +if( UniformRandom() < p ) +{ + double theta = UniformRandom(); +} +``` + +If the key were only `(cell_id, time_step)`, both calls would ask for the same value. That is wrong. + +`sub_index` fixes this by making the first call use sub-index 0 and the second call use sub-index 1. + +### Problem 2: different purposes in the same step + +Suppose in the same saved step you do both of these for the same cell: + +- phenotype transition draw, +- motility direction draw. + +If the key were only `(cell_id, time_step)`, those unrelated events would collide too. + +`purpose` fixes this by giving each update phase its own namespace. + +## About "the internal index will always start from 0 after restart" + +That statement is only partly true. + +It is true that `sub_index` should normally restart at 0 when entering a new deterministic context. + +It is not true that this makes `purpose` unnecessary. + +Why: + +- in one resumed step, cell 10 can enter several different stochastic contexts, +- each of those contexts may start its own `sub_index` at 0, +- but they must still remain distinct from each other. + +So after restart you can absolutely do this: + +- phenotype context for cell 10, step 25: sub-index starts at 0, +- velocity context for cell 10, step 25: sub-index also starts at 0. + +That is correct only if `purpose` is also part of the key. Otherwise both contexts collide. + +So the right rule is: + +- `sub_index` resets to 0 per context, +- not per whole simulation step. + +## Restart implication + +If restart always begins from a saved snapshot, you may be able to derive `step_key` from the restored snapshot state rather than saving a separate hidden counter. But you still need the full four-role scheme: + +- identity of the cell, +- identity of the resumed step, +- identity of the stochastic purpose, +- identity of the draw number within that purpose. + +If any one of those roles is missing, different random events can collapse onto the same keyed draw. + +## The existing precedent: ordering division and death + +The counter-based RNG work already had to solve a version of this problem for cell division and death, and the fixes below follow the same pattern it established. + +Division and death are decided inside the parallel phenotype loop, but a cell can't safely divide or remove itself from `(*all_cells)` while that loop is still running on other threads. So instead of acting immediately, a cell calls `flag_cell_for_division()` / `flag_cell_for_removal()`, which push it onto `cells_ready_to_divide` / `cells_ready_to_die` under `#pragma omp critical`. Only after the parallel loop finishes does [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp) actually call `divide()` / `die()` on those lists — serially, and, when `PhysiCell_settings.use_counter_based_rng` is enabled, only after sorting each list by cell ID first: + +```cpp +if ( PhysiCell_settings.use_counter_based_rng ) + std::sort( cells_ready_to_divide.begin(), cells_ready_to_divide.end(), + []( const Cell* lhs, const Cell* rhs ) { return lhs->ID < rhs->ID; } ); +for( int i=0; i < cells_ready_to_divide.size(); i++ ) +{ + activate_random_context( cells_ready_to_divide[i]->ID, RANDOM_PURPOSE_DIVISION ); + cells_ready_to_divide[i]->divide(); + clear_random_context(); +} +``` + +The sort matters for two reasons: new daughter cells get their IDs assigned in whatever order `divide()` is called, so without it, which physical cell gets which ID would depend on thread scheduling; and `die()` removes cells from `(*all_cells)` by swapping in the last element, so processing removals in a different order leaves `(*all_cells)` in a different arrangement, which then affects the deterministic per-cell RNG keys and any later code that iterates `(*all_cells)` in order. + +This "collect under a critical section while parallel, then apply serially in a fixed order" shape is exactly what Fixes 1 and 2 below extend to cell-cell interactions and spring attachments, which had no such ordering at all. + +## Beyond RNG determinism: cross-cell mutation races + +Making `Random()` deterministic per `(cell_id, step_key, purpose, sub_index)` only guarantees that a cell's *own* draws are reproducible. It does nothing to protect a cell from having its state changed by a *different* cell's thread while both are being processed in the same `#pragma omp parallel for`. That is a separate class of bug, and thread-reproducibility testing on real projects (worm, virus_macrophage, and a rules-driven virus/macrophage model) turned up several of them in the core engine, independent of RNG. + +### The general pattern + +Several core functions read another cell's current state, decide what to do based on it, and only then commit a change, e.g.: + +```cpp +if( pTarget->something < threshold ) // read, unprotected +{ + pTarget->something = new_value; // write, protected by #pragma omp critical +} +``` + +The `critical` block (where present) stops two threads from writing at the exact same instant, so it prevents outright memory corruption. It does *not* stop the read from seeing a stale value: two different cells' threads can both perform the read before either one's write lands, both conclude they should act, and then both act, in whichever order the OpenMP scheduler happens to pick. In a 1-thread run this can never happen, because the loop is always fully sequential, so the "reference" run and a multithreaded run can permanently disagree from that point on. + +`#pragma omp critical` alone fixes memory-safety. It does not fix determinism, because it says nothing about *which* thread's read/write pair goes first. `#pragma omp ordered` fixes both: it forces the loop's iterations to commit their marked region in the same relative order a 1-thread run would use, no matter how many threads are active or how they're scheduled. The fixes below use `ordered` wherever the operation is not commutative (an addition, a "first one wins" decision), and plain `critical` where it is (a `max` update, where the result doesn't depend on which thread got there first). + +### Fix 1: `standard_cell_cell_interactions` — attack, phagocytosis, fusion + +[`core/PhysiCell_standard_models.cpp`](core/PhysiCell_standard_models.cpp) calls `pCell->ingest_cell(pTarget)`, `pCell->attack_cell(pTarget, dt)`, and `pCell->fuse_cell(pTarget)` directly from a per-cell loop. All three mutate the *target* cell's live state (volume, damage, death flag, position) and are wrapped in `#pragma omp critical` internally — but the guard that decides whether to act (e.g. `pCell_to_eat->phenotype.volume.total < 1e-15`) runs *before* that critical section. Two different cells can both decide to consume, attack, or fuse with the same target in the same step; which one's action actually lands first, and therefore what the target's final state is, depends on thread scheduling. + +Fixed by wrapping the call to `standard_cell_cell_interactions()` in [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp) in `#pragma omp ordered` (loop changed to `schedule(static) ordered`). This also covers the function's early-exit check on `phenotype.death.dead`, which was itself an unsynchronized read of a flag another cell's thread could be writing at the same time. + +### Fix 2: `dynamic_spring_attachments` — attachment-capacity race + +[`core/PhysiCell_standard_models.cpp`](core/PhysiCell_standard_models.cpp) checks `pTest->state.spring_attachments.size() < pTest->phenotype.mechanics.maximum_number_of_attachments` to decide whether a neighbor still has room for a new attachment, then calls `attach_cells_as_spring()`. The capacity check is never re-validated once the lock inside `attach_cells_as_spring()` is actually held. Two different cells can both see "room for one more" on the same neighbor at once and both attach, silently pushing that cell past its configured maximum, in a way that depends on scheduling. + +Fixed the same way: the call to `dynamic_spring_attachments()` in [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp) is now inside `#pragma omp ordered`. + +### Fix 3: `max_cell_interactive_distance_in_voxel` — stale/racy neighbor-search bound + +Each mechanics voxel keeps a running maximum of `radius * relative_maximum_adhesion_distance` over the cells that have been recorded there, used to decide whether a neighbor search needs to look inside that voxel at all. It's updated with an unprotected read-compare-write (`if (current < new) current = new;`) in `Cell::convert_to_cell_definition()` (fires on a rule-triggered type transformation) and, matching upstream PR [#409](https://github.com/MathCancer/PhysiCell/pull/409), in `Cell::update_voxel_in_container()` (fires on ordinary movement between voxels — without this, the bound goes stale as cells migrate, and a voxel's neighbor search can wrongly skip a cell that has since moved in). Both live in [core/PhysiCell_cell.cpp](core/PhysiCell_cell.cpp). + +Unlike Fixes 1 and 2, this doesn't need `ordered`: a `max` over a fixed set of candidate values is the same regardless of the order you compare them in, so a plain `#pragma omp critical` around the read-compare-write is enough to make it both race-free and deterministic. + +### Fix 4: cell secretion/uptake into shared diffusion voxels + +This was the one that actually explained the residual, hard-to-pin-down mismatches in a rules-driven multi-cell-type model, after Fixes 1–3 were already in place and confirmed to not be the cause (isolated by testing with those mechanisms enabled, disabled, and recombined). + +PhysiCell's real per-cell secretion path is `Cell_Container::update_all_cells()`'s first loop → `Secretion::advance()` → `Basic_Agent::simulate_secretion_and_uptake()` (in [BioFVM/BioFVM_basic_agent.cpp](BioFVM/BioFVM_basic_agent.cpp)), which does: + +```cpp +(*pS)(current_voxel_index) += cell_source_sink_solver_temp1; +(*pS)(current_voxel_index) /= cell_source_sink_solver_temp2; +(*pS)(current_voxel_index) += cell_source_sink_solver_temp_export2; +``` + +directly on the microenvironment's shared voxel density vector, with no lock at all. Diffusion voxels are typically larger than a single cell, so it's normal for several cells to share one; when two cells in the same voxel are processed by different threads at the same time, this isn't just reordering, it's a genuine lost update — one thread's contribution can be silently overwritten instead of added. + +Fixed by wrapping the call to `phenotype.secretion.advance(...)` in [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp) in `#pragma omp ordered`. + +## How to enable counter-based RNG + +Counter-based RNG is opt-in and **off by default**. If `` is not set at all, the simulation uses the legacy thread-local generator (`PhysiCell_settings.use_counter_based_rng = false` in [modules/PhysiCell_settings.h](modules/PhysiCell_settings.h)), which is not guaranteed to be reproducible across different thread counts. + +To enable it, add (or set) an `` element inside `` in the config XML: + +```xml + + ... + counter_based + +``` + +Accepted values, parsed in [modules/PhysiCell_settings.cpp](modules/PhysiCell_settings.cpp): + +- `counter_based`, `counter`, or `philox` → enables the deterministic Philox-based path (`use_counter_based_rng = true`). +- `legacy`, `thread_local`, or `mt19937` → explicitly selects the old, non-deterministic thread-local generator (`use_counter_based_rng = false`). +- Omitted or empty → stays at the default, `false`. + +An unrecognized value prints an error and leaves the setting unchanged. + +## Testing a new model for thread reproducibility + +[beta/test_thread_repro.py](beta/test_thread_repro.py) runs an already-built executable twice, once per thread count, and diffs the output. It does not build the project for you — build first (`make`), then run: + +```bash +python beta/test_thread_repro.py +``` + +Example, testing the currently-loaded project with 1 vs. 4 threads for 1440 simulated minutes: + +```bash +python beta/test_thread_repro.py project config/PhysiCell_settings.xml 1440 1 4 counter_based local_runs +``` + +- `executable`: path to the built binary (relative to the repo root), e.g. `project`. +- `config_file`: the XML config to copy and use for both runs. +- `max_time`: simulated end time, in the config's time units. +- `threads_a`, `threads_b`: the two thread counts to compare (must differ). +- `rng_mode`: written into each run's copy of `` before running (omit, or pass `off`/`none`/`default`, to leave the config's existing setting untouched). +- `work_dir`: where the two runs' output is written (each invocation creates a timestamped subfolder inside it). + +The script runs both thread counts, then diffs the final XML/SVG snapshots with [beta/test_diff_svg.py](beta/test_diff_svg.py) and reports the first mismatch, if any. + +Because these races can be intermittent — several of the fixes above only failed some fraction of the time — a single passing run is not strong evidence on its own. Repeat the same command several times (5–10x) before trusting a "pass," especially for a new or modified model. + +## Validation + +The new system tests pass with 1 thread and 4 threads for the template project. + +The four cross-cell-mutation fixes above were found and validated by testing the template project with multiple cell behaviors activated (phagocytosis, attack, transform, and others) using `beta/test_thread_repro.py` across single and multiple threads. Repeated trials failed consistently before Fix 4 (the secretion race) was in place; after all four fixes were applied, repeated trials passed bit-exact. + +Known open item: a project-level (not core-engine) bug remains in the worm sample's custom `contact_function` ([custom_modules/custom.cpp](custom_modules/custom.cpp) when the worm project is loaded), which writes into an attached cell's `custom_data` from within a parallel loop without going through the engine's RNG/ordering machinery. That is a project-code issue rather than a PhysiCell core issue and is tracked separately. From 85ae4ace42b47b8fd34b849a7e5047f3a81575a7 Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Thu, 6 Aug 2026 13:09:54 -0400 Subject: [PATCH 5/9] Merge PR #409: keep max_cell_interactive_distance_in_voxel current on move Update the bound in update_voxel_in_container() whenever a cell crosses into a new voxel, not just on the rare type-conversion path. Without this the bound goes stale as cells migrate, and a voxel's neighbor search can wrongly skip a cell that has since moved in. Adapted from https://github.com/MathCancer/PhysiCell/pull/409, with #pragma omp critical added around the read-compare-write: in this codebase update_voxel_in_container() is also reachable from fuse_cell()'s parallel context, so the unprotected version from the upstream PR would still race here even though it may not have in the original context. --- core/PhysiCell_cell.cpp | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp index 1f3e2646b..997e42fbf 100644 --- a/core/PhysiCell_cell.cpp +++ b/core/PhysiCell_cell.cpp @@ -943,9 +943,25 @@ void Cell::update_voxel_in_container() container->add_agent_to_voxel(this, updated_current_mechanics_voxel_index); } current_mechanics_voxel_index=updated_current_mechanics_voxel_index; + + // keep max_cell_interactive_distance_in_voxel current as cells migrate between voxels (PR 409), + // not just on the rare type-conversion path (convert_to_cell_definition()). Without this, the + // bound goes stale as cells move, which can make a voxel's neighbor search wrongly skip a cell + // that has since moved in -- a real correctness gap, not just a reproducibility one. Protected + // for the same reason as convert_to_cell_definition() (max is order-independent, critical alone + // is sufficient -- no "ordered" needed). + #pragma omp critical + { + if( get_container()->max_cell_interactive_distance_in_voxel[current_mechanics_voxel_index] < + phenotype.geometry.radius * phenotype.mechanics.relative_maximum_adhesion_distance ) + { + get_container()->max_cell_interactive_distance_in_voxel[current_mechanics_voxel_index] = phenotype.geometry.radius + * phenotype.mechanics.relative_maximum_adhesion_distance; + } + } } - - return; + + return; } void Cell::copy_data(Cell* copy_me) From ebdb575bbd473aaa908c56ba2193fcfbe8814d57 Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Thu, 6 Aug 2026 13:14:41 -0400 Subject: [PATCH 6/9] Fix unprotected max_cell_interactive_distance_in_voxel race in convert_to_cell_definition Unrelated to PR #409: this read-compare-write already existed on the type-conversion path and was never protected. convert_to_cell_definition() is reachable from rule-triggered transformations inside the (unordered) phenotype loop, so two cells transforming in the same voxel in the same step could race on this shared bound. Protected with #pragma omp critical; sufficient since max doesn't depend on evaluation order. --- core/PhysiCell_cell.cpp | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp index 997e42fbf..926420f19 100644 --- a/core/PhysiCell_cell.cpp +++ b/core/PhysiCell_cell.cpp @@ -1204,17 +1204,25 @@ void Cell::convert_to_cell_definition( Cell_Definition& cd ) // phenotype.geometry.update( this, phenotype, 0.0 ); // not necessary since we copy geometry above */ - // Here the current mechanics voxel index may not be initialized, when position is still unknown. + // Here the current mechanics voxel index may not be initialized, when position is still unknown. + // This can run from convert_to_cell_definition(), which is reachable from rule-triggered type + // transformations inside the (unordered) phenotype loop -- max_cell_interactive_distance_in_voxel + // is shared across all cells in the voxel, so the read-compare-write must be atomic. A max is + // order-independent (unlike a sum), so a plain critical section is sufficient here -- no need for + // "ordered". if (get_current_mechanics_voxel_index() >= 0) { - if( get_container()->max_cell_interactive_distance_in_voxel[get_current_mechanics_voxel_index()] < - phenotype.geometry.radius * phenotype.mechanics.relative_maximum_adhesion_distance ) - { - get_container()->max_cell_interactive_distance_in_voxel[get_current_mechanics_voxel_index()] = phenotype.geometry.radius - * phenotype.mechanics.relative_maximum_adhesion_distance; - } + #pragma omp critical + { + if( get_container()->max_cell_interactive_distance_in_voxel[get_current_mechanics_voxel_index()] < + phenotype.geometry.radius * phenotype.mechanics.relative_maximum_adhesion_distance ) + { + get_container()->max_cell_interactive_distance_in_voxel[get_current_mechanics_voxel_index()] = phenotype.geometry.radius + * phenotype.mechanics.relative_maximum_adhesion_distance; + } + } } - return; + return; } void delete_cell( int index ) From 2e5a384e34ff340685864f8e62da72f063303af8 Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Fri, 7 Aug 2026 11:36:08 -0400 Subject: [PATCH 7/9] refined text --- counter_based_rng.md | 66 +------------------------------------------- 1 file changed, 1 insertion(+), 65 deletions(-) diff --git a/counter_based_rng.md b/counter_based_rng.md index 18282aa9a..c069065e4 100644 --- a/counter_based_rng.md +++ b/counter_based_rng.md @@ -8,9 +8,7 @@ The core RNG layer now has a keyed entry point in [core/PhysiCell_utilities.h](c The main simulation loop sets that context around the major per-cell update phases in [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp). That means the existing simulation code can keep calling the usual random helpers, but those helpers will now draw from a reproducible per-cell, per-step sequence instead of a thread-order-dependent stream. -I also added a reusable end-to-end thread reproducibility checker in [beta/test_thread_repro.py](beta/test_thread_repro.py). The GitHub Actions workflow uses it for the template and intracellular template jobs in [tests.yml](.github/workflows/tests.yml). - -The checker runs the same project twice with different thread counts and compares the final XML output after normalizing volatile metadata. +I also added a reusable end-to-end thread reproducibility checker in [beta/test_thread_repro.py](beta/test_thread_repro.py). The checker runs the same project twice with different thread counts and compares the final XML output after normalizing volatile metadata. ## Why this did not require changing every `Random()` call @@ -119,64 +117,6 @@ There are also mixed code paths where one logical operation first makes a decisi This is why `sub_index` is convenient even when `purpose` is automatic. A single semantic purpose such as `motility update` can safely contain several actual draws without inventing a separate micro-purpose for every internal step. -## Why `cell_id` and `time_step` alone are not enough - -They are not enough if more than one random draw can happen for the same cell in the same step, which is normal in PhysiCell. - -There are two independent problems. - -### Problem 1: multiple draws within one context - -Suppose cell 10 is in phenotype step 25 and a function does this: - -```cpp -if( UniformRandom() < p ) -{ - double theta = UniformRandom(); -} -``` - -If the key were only `(cell_id, time_step)`, both calls would ask for the same value. That is wrong. - -`sub_index` fixes this by making the first call use sub-index 0 and the second call use sub-index 1. - -### Problem 2: different purposes in the same step - -Suppose in the same saved step you do both of these for the same cell: - -- phenotype transition draw, -- motility direction draw. - -If the key were only `(cell_id, time_step)`, those unrelated events would collide too. - -`purpose` fixes this by giving each update phase its own namespace. - -## About "the internal index will always start from 0 after restart" - -That statement is only partly true. - -It is true that `sub_index` should normally restart at 0 when entering a new deterministic context. - -It is not true that this makes `purpose` unnecessary. - -Why: - -- in one resumed step, cell 10 can enter several different stochastic contexts, -- each of those contexts may start its own `sub_index` at 0, -- but they must still remain distinct from each other. - -So after restart you can absolutely do this: - -- phenotype context for cell 10, step 25: sub-index starts at 0, -- velocity context for cell 10, step 25: sub-index also starts at 0. - -That is correct only if `purpose` is also part of the key. Otherwise both contexts collide. - -So the right rule is: - -- `sub_index` resets to 0 per context, -- not per whole simulation step. - ## Restart implication If restart always begins from a saved snapshot, you may be able to derive `step_key` from the restored snapshot state rather than saving a separate hidden counter. But you still need the full four-role scheme: @@ -190,8 +130,6 @@ If any one of those roles is missing, different random events can collapse onto ## The existing precedent: ordering division and death -The counter-based RNG work already had to solve a version of this problem for cell division and death, and the fixes below follow the same pattern it established. - Division and death are decided inside the parallel phenotype loop, but a cell can't safely divide or remove itself from `(*all_cells)` while that loop is still running on other threads. So instead of acting immediately, a cell calls `flag_cell_for_division()` / `flag_cell_for_removal()`, which push it onto `cells_ready_to_divide` / `cells_ready_to_die` under `#pragma omp critical`. Only after the parallel loop finishes does [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp) actually call `divide()` / `die()` on those lists — serially, and, when `PhysiCell_settings.use_counter_based_rng` is enabled, only after sorting each list by cell ID first: ```cpp @@ -249,8 +187,6 @@ Unlike Fixes 1 and 2, this doesn't need `ordered`: a `max` over a fixed set of c ### Fix 4: cell secretion/uptake into shared diffusion voxels -This was the one that actually explained the residual, hard-to-pin-down mismatches in a rules-driven multi-cell-type model, after Fixes 1–3 were already in place and confirmed to not be the cause (isolated by testing with those mechanisms enabled, disabled, and recombined). - PhysiCell's real per-cell secretion path is `Cell_Container::update_all_cells()`'s first loop → `Secretion::advance()` → `Basic_Agent::simulate_secretion_and_uptake()` (in [BioFVM/BioFVM_basic_agent.cpp](BioFVM/BioFVM_basic_agent.cpp)), which does: ```cpp From 52ae2a8c3f03fde0c3118594a5c11e68ad27b947 Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Fri, 7 Aug 2026 13:51:57 -0400 Subject: [PATCH 8/9] Gate ordered interaction/attachment/secretion phases at the call site Move the use_counter_based_rng branch out of run_secretion_phase(), run_spring_attachment_phase(), and run_cell_cell_interactions_phase() and into their call sites in update_all_cells(). Each is now if(use_counter_based_rng){ call a small *_ordered() function } else { original inline loop }, so the legacy path stays untouched inline code instead of living inside a shared function, and use_counter_based_rng is still checked once per phase per step, not once per cell. Also drop the activate_random_context()/clear_random_context() calls from the three else branches: inside a branch that only runs when use_counter_based_rng is already false, that guarded call is provably dead code, so removing it makes the legacy path true original behavior rather than calls that happen to no-op. --- core/PhysiCell_cell_container.cpp | 126 +++++++++++++++++++++++------- 1 file changed, 97 insertions(+), 29 deletions(-) diff --git a/core/PhysiCell_cell_container.cpp b/core/PhysiCell_cell_container.cpp index 97262b04c..10698acdc 100644 --- a/core/PhysiCell_cell_container.cpp +++ b/core/PhysiCell_cell_container.cpp @@ -96,6 +96,68 @@ constexpr std::uint64_t RANDOM_PURPOSE_CELL_CELL = 9; constexpr std::uint64_t RANDOM_PURPOSE_POSITION = 10; constexpr std::uint64_t RANDOM_PURPOSE_DIVISION = 11; constexpr std::uint64_t RANDOM_PURPOSE_DEATH = 12; + +// The three functions below are the "ordered" counterpart of a loop that still exists, unchanged, +// inline in Cell_Container::update_all_cells() (the plain #pragma omp parallel for version). Each +// is called from update_all_cells() only when PhysiCell_settings.use_counter_based_rng is enabled, +// selected once per phase per step -- not per cell. They protect a phase whose function can mutate +// a *different* cell's live state from inside a parallel per-cell loop (secretion into a shared +// voxel, spring-attachment capacity, attack/ingest/fuse), by forcing the same relative commit order +// a single-thread run would use. Legacy (non-counter-based) runs take the unchanged, fully parallel +// branch instead, keeping their existing performance profile at the cost of leaving these races +// present in that mode. See protocols/counter_based_rng.md for the full discussion of this trade-off, +// including why "ordered" (not just "critical") is required for the secretion and interactions +// phases specifically. + +void run_secretion_phase_ordered( double diffusion_dt_, std::uint64_t random_step ) +{ + #pragma omp parallel for schedule(static) ordered + for( int i=0; i < (*all_cells).size(); i++ ) + { + if( (*all_cells)[i]->is_out_of_domain == false ) + { + set_deterministic_random_context( (*all_cells)[i]->ID, random_step, RANDOM_PURPOSE_SECRETION ); + #pragma omp ordered + { + (*all_cells)[i]->phenotype.secretion.advance( (*all_cells)[i], (*all_cells)[i]->phenotype , diffusion_dt_ ); + } + clear_deterministic_random_context(); + } + } + return; +} + +void run_spring_attachment_phase_ordered( double time_since_last_mechanics, std::uint64_t random_step ) +{ + #pragma omp parallel for schedule(static) ordered + for( int i=0; i < (*all_cells).size(); i++ ) + { + Cell* pC = (*all_cells)[i]; + set_deterministic_random_context( pC->ID, random_step, RANDOM_PURPOSE_SPRINGS ); + #pragma omp ordered + { + dynamic_spring_attachments(pC,pC->phenotype,time_since_last_mechanics); + } + clear_deterministic_random_context(); + } + return; +} + +void run_cell_cell_interactions_phase_ordered( double time_since_last_mechanics, std::uint64_t random_step ) +{ + #pragma omp parallel for schedule(static) ordered + for( int i=0; i < (*all_cells).size(); i++ ) + { + Cell* pC = (*all_cells)[i]; + set_deterministic_random_context( pC->ID, random_step, RANDOM_PURPOSE_CELL_CELL ); + #pragma omp ordered + { + standard_cell_cell_interactions(pC,pC->phenotype,time_since_last_mechanics); + } + clear_deterministic_random_context(); + } + return; +} } Cell_Container::Cell_Container() @@ -164,24 +226,24 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me // BioFVM/BioFVM_basic_agent.cpp). Diffusion voxels are typically larger than a single cell, so // multiple cells commonly share one; without a fixed commit order, two cells in the same voxel // processed by different threads at the same time race on that shared read-modify-write (a genuine - // lost-update, not just reordering). This is the actual call path PhysiCell uses for cell secretion - // (unlike Microenvironment::simulate_cell_sources_and_sinks, which is never invoked by the core - // simulation loop) -- forcing the same relative commit order as a single-thread run removes the race. - - #pragma omp parallel for schedule(static) ordered - for( int i=0; i < (*all_cells).size(); i++ ) + // lost-update, not just reordering). Only taken when counter_based_rng is enabled; see + // protocols/counter_based_rng.md. + if( use_counter_based_rng ) { - if( (*all_cells)[i]->is_out_of_domain == false ) + run_secretion_phase_ordered( diffusion_dt_, random_step ); + } + else + { + #pragma omp parallel for + for( int i=0; i < (*all_cells).size(); i++ ) { - activate_random_context( (*all_cells)[i]->ID, RANDOM_PURPOSE_SECRETION ); - #pragma omp ordered + if( (*all_cells)[i]->is_out_of_domain == false ) { (*all_cells)[i]->phenotype.secretion.advance( (*all_cells)[i], (*all_cells)[i]->phenotype , diffusion_dt_ ); } - clear_random_context(); } } - + //if it is the time for running cell cycle, do it! double time_since_last_cycle= t- last_cell_cycle_time; @@ -337,17 +399,20 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me // detach_cells_as_spring() on that neighbor. The capacity check is never re-validated once // the lock inside attach_cells_as_spring() is actually held, so without a fixed commit order // two different cells can both see "room for one more" on the same neighbor at once and both - // attach, pushing it past its intended maximum -- a thread-count-dependent result. - #pragma omp parallel for schedule(static) ordered - for( int i=0; i < (*all_cells).size(); i++ ) + // attach, pushing it past its intended maximum. Only taken when counter_based_rng is + // enabled; see protocols/counter_based_rng.md. + if( use_counter_based_rng ) { - Cell* pC = (*all_cells)[i]; - activate_random_context( pC->ID, RANDOM_PURPOSE_SPRINGS ); - #pragma omp ordered + run_spring_attachment_phase_ordered( time_since_last_mechanics, random_step ); + } + else + { + #pragma omp parallel for + for( int i=0; i < (*all_cells).size(); i++ ) { + Cell* pC = (*all_cells)[i]; dynamic_spring_attachments(pC,pC->phenotype,time_since_last_mechanics); } - clear_random_context(); } #pragma omp parallel for for( int i=0; i < (*all_cells).size(); i++ ) @@ -369,21 +434,24 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me // new March 2022: // run standard interactions (phagocytosis, attack, fusion) here - // this loop is "ordered": standard_cell_cell_interactions() can ingest, attack, or fuse - // a *different* cell (mutating its live state), and it early-exits based on that cell's - // current death flag. Enforcing the same relative commit order as a single-thread run - // (via #pragma omp ordered) makes those cross-cell interactions thread-count independent, - // instead of depending on however OpenMP happens to schedule the threads this run. - #pragma omp parallel for schedule(static) ordered - for( int i=0; i < (*all_cells).size(); i++ ) + // ordered: standard_cell_cell_interactions() can ingest, attack, or fuse a *different* cell + // (mutating its live state), and it early-exits based on that cell's current death flag. + // Enforcing the same relative commit order as a single-thread run makes those cross-cell + // interactions thread-count independent, instead of depending on however OpenMP happens to + // schedule the threads this run. Only taken when counter_based_rng is enabled; see + // protocols/counter_based_rng.md. + if( use_counter_based_rng ) { - Cell* pC = (*all_cells)[i]; - activate_random_context( pC->ID, RANDOM_PURPOSE_CELL_CELL ); - #pragma omp ordered + run_cell_cell_interactions_phase_ordered( time_since_last_mechanics, random_step ); + } + else + { + #pragma omp parallel for + for( int i=0; i < (*all_cells).size(); i++ ) { + Cell* pC = (*all_cells)[i]; standard_cell_cell_interactions(pC,pC->phenotype,time_since_last_mechanics); } - clear_random_context(); } // super-critical to performance! clear the "dummy" cells from phagocytosis / fusion // otherwise, comptuational cost increases at polynomial rate VERY fast, as O(10,000) From 5438cca9991a766b086440539367703fda9a48e5 Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Fri, 7 Aug 2026 13:53:24 -0400 Subject: [PATCH 9/9] Move counter_based_rng.md to protocols/ and update it for the current code. Also bring the content in line with the current architecture: add "Why this matters" (reproducibility, debugging, potential bitwise restart) and "Credit" (D. E. Shaw Research / Random123, philox.h provenance) sections, correct Fixes 1/2/4 to name the *_ordered() functions and describe gating at the call site instead of inside a shared function, and document why the else branches no longer call activate_random_context()/clear_random_context(). --- .../counter_based_rng.md | 80 +++++++++++++------ 1 file changed, 55 insertions(+), 25 deletions(-) rename counter_based_rng.md => protocols/counter_based_rng.md (55%) diff --git a/counter_based_rng.md b/protocols/counter_based_rng.md similarity index 55% rename from counter_based_rng.md rename to protocols/counter_based_rng.md index c069065e4..1b279b8fe 100644 --- a/counter_based_rng.md +++ b/protocols/counter_based_rng.md @@ -2,13 +2,29 @@ This branch adds a deterministic RNG path based on Philox4x32-10 so the same simulation can produce the same random values even when the OpenMP thread count changes. +## Why this matters + +**Reproducibility.** Same config, same seed, same result — regardless of how many OpenMP threads a run happens to use. Before this work, PhysiCell's thread-local RNG made results depend on thread count, and (for the four races described below) sometimes even on run-to-run scheduling luck at a *fixed* thread count. A figure or a comparison between two runs could differ for reasons that had nothing to do with the model — purely because of a hardware/thread-count choice. With counter-based RNG and the race fixes in place, the number of cores a run happens to use stops being a hidden variable in the result. + +**Debugging.** Non-reproducibility isn't just a publishing inconvenience — it made a whole class of core-engine bugs undetectable, because there was never a reliable "this should be identical" baseline to test against. Every race fixed in this branch (unprotected shared-voxel secretion, an attachment-capacity check that's never re-validated, cross-cell mutation races in attack/ingest/fuse) was a genuine, silent correctness bug — present regardless of RNG mode, not something this branch introduced — that had gone unnoticed because nothing before could reliably say "these two runs should match and don't." That comparison is now a real, repeatable regression test ([beta/test_thread_repro.py](../beta/test_thread_repro.py)) instead of something you could only ever suspect informally. + +**Potential bitwise restart.** The `(cell_id, step_key, purpose, sub_index)` keyed scheme (see "Minimal key scheme for restart" below) is deliberately restart-friendly: if `step_key` can be reconstructed exactly from restart-stable state rather than an in-memory-only counter, a resumed run could in principle produce bit-identical output to the same run had it never been interrupted. That's *not* implemented or tested here — restart itself is out of scope for this branch — but the keying scheme is designed so it doesn't foreclose that possibility, which a plain thread-local generator has no way to offer no matter how restart is later implemented. + +## Credit + +The counter-based generator itself, Philox4x32-10, and the counter-based RNG design it comes from, are not new here — they're from D. E. Shaw Research's [Random123](https://www.deshawresearch.com/resources_random123.html) library: + +> John K. Salmon, Mark A. Moraes, Ron O. Dror, and David E. Shaw. "Parallel Random Numbers: As Easy as 1, 2, 3." In *Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis* (SC '11). ACM, 2011. + +[`modules/philox.h`](../modules/philox.h) is vendored from Random123's `philox.h`, with D. E. Shaw Research's original BSD-style copyright notice and license preserved in the file. This branch's contribution is applying that generator inside PhysiCell — keying it by `(cell_id, time_step, purpose, sub_index)`, wiring a deterministic context around the per-cell update phases, and fixing the separate cross-cell mutation races (below) that RNG determinism alone doesn't address — not the counter-based RNG construction itself. + ## What changed -The core RNG layer now has a keyed entry point in [core/PhysiCell_utilities.h](core/PhysiCell_utilities.h) and [core/PhysiCell_utilities.cpp](core/PhysiCell_utilities.cpp). The new `Random(cell_id, time_step, purpose, sub_index)` API generates a value directly from the key and counter inputs, and the existing `UniformRandom()`, `UniformInt()`, `NormalRandom()`, and `LogNormalRandom()` functions can read from a thread-local deterministic context when that context is active. +The core RNG layer now has a keyed entry point in [core/PhysiCell_utilities.h](../core/PhysiCell_utilities.h) and [core/PhysiCell_utilities.cpp](../core/PhysiCell_utilities.cpp). The new `Random(cell_id, time_step, purpose, sub_index)` API generates a value directly from the key and counter inputs, and the existing `UniformRandom()`, `UniformInt()`, `NormalRandom()`, and `LogNormalRandom()` functions can read from a thread-local deterministic context when that context is active. -The main simulation loop sets that context around the major per-cell update phases in [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp). That means the existing simulation code can keep calling the usual random helpers, but those helpers will now draw from a reproducible per-cell, per-step sequence instead of a thread-order-dependent stream. +The main simulation loop sets that context around the major per-cell update phases in [core/PhysiCell_cell_container.cpp](../core/PhysiCell_cell_container.cpp). That means the existing simulation code can keep calling the usual random helpers, but those helpers will now draw from a reproducible per-cell, per-step sequence instead of a thread-order-dependent stream. -I also added a reusable end-to-end thread reproducibility checker in [beta/test_thread_repro.py](beta/test_thread_repro.py). The checker runs the same project twice with different thread counts and compares the final XML output after normalizing volatile metadata. +I also added a reusable end-to-end thread reproducibility checker in [beta/test_thread_repro.py](../beta/test_thread_repro.py). The checker runs the same project twice with different thread counts and compares the final XML output after normalizing volatile metadata. ## Why this did not require changing every `Random()` call @@ -104,16 +120,16 @@ Without `sub_index`, both calls would map to the same keyed random value. This is not just a theoretical issue. Several existing helpers already consume more than one uniform draw for one higher-level stochastic operation. -- `NormalRandom()` uses two uniforms internally through the Box-Muller transform. See [core/PhysiCell_utilities.cpp](core/PhysiCell_utilities.cpp). -- `LogNormalRandom()` also uses that same two-uniform path because it is built on top of `NormalRandom()`. See [core/PhysiCell_utilities.cpp](core/PhysiCell_utilities.cpp). -- `UniformOnUnitSphere()` uses two uniforms, one for `z` and one for `theta`. See [core/PhysiCell_utilities.cpp](core/PhysiCell_utilities.cpp). -- `LegacyRandomOnUnitSphere()` uses multiple uniforms as well. See [core/PhysiCell_utilities.cpp](core/PhysiCell_utilities.cpp). -- `UniformInUnitDisc()` uses two uniforms. See [core/PhysiCell_rules.cpp](core/PhysiCell_rules.cpp). -- `UniformInUnitSphere()` uses three uniforms. See [core/PhysiCell_rules.cpp](core/PhysiCell_rules.cpp). -- `UniformInAnnulus()` uses two uniforms. See [core/PhysiCell_rules.cpp](core/PhysiCell_rules.cpp). -- `UniformInShell()` uses three uniforms. See [core/PhysiCell_rules.cpp](core/PhysiCell_rules.cpp). +- `NormalRandom()` uses two uniforms internally through the Box-Muller transform. See [core/PhysiCell_utilities.cpp](../core/PhysiCell_utilities.cpp). +- `LogNormalRandom()` also uses that same two-uniform path because it is built on top of `NormalRandom()`. See [core/PhysiCell_utilities.cpp](../core/PhysiCell_utilities.cpp). +- `UniformOnUnitSphere()` uses two uniforms, one for `z` and one for `theta`. See [core/PhysiCell_utilities.cpp](../core/PhysiCell_utilities.cpp). +- `LegacyRandomOnUnitSphere()` uses multiple uniforms as well. See [core/PhysiCell_utilities.cpp](../core/PhysiCell_utilities.cpp). +- `UniformInUnitDisc()` uses two uniforms. See [core/PhysiCell_rules.cpp](../core/PhysiCell_rules.cpp). +- `UniformInUnitSphere()` uses three uniforms. See [core/PhysiCell_rules.cpp](../core/PhysiCell_rules.cpp). +- `UniformInAnnulus()` uses two uniforms. See [core/PhysiCell_rules.cpp](../core/PhysiCell_rules.cpp). +- `UniformInShell()` uses three uniforms. See [core/PhysiCell_rules.cpp](../core/PhysiCell_rules.cpp). -There are also mixed code paths where one logical operation first makes a decision draw and then consumes additional random values inside the same higher-level update. For example, motility update can first draw whether the cell reorients and then, if it does, draw a random direction on the circle or sphere. See [core/PhysiCell_cell.cpp](core/PhysiCell_cell.cpp). +There are also mixed code paths where one logical operation first makes a decision draw and then consumes additional random values inside the same higher-level update. For example, motility update can first draw whether the cell reorients and then, if it does, draw a random direction on the circle or sphere. See [core/PhysiCell_cell.cpp](../core/PhysiCell_cell.cpp). This is why `sub_index` is convenient even when `purpose` is automatic. A single semantic purpose such as `motility update` can safely contain several actual draws without inventing a separate micro-purpose for every internal step. @@ -130,7 +146,7 @@ If any one of those roles is missing, different random events can collapse onto ## The existing precedent: ordering division and death -Division and death are decided inside the parallel phenotype loop, but a cell can't safely divide or remove itself from `(*all_cells)` while that loop is still running on other threads. So instead of acting immediately, a cell calls `flag_cell_for_division()` / `flag_cell_for_removal()`, which push it onto `cells_ready_to_divide` / `cells_ready_to_die` under `#pragma omp critical`. Only after the parallel loop finishes does [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp) actually call `divide()` / `die()` on those lists — serially, and, when `PhysiCell_settings.use_counter_based_rng` is enabled, only after sorting each list by cell ID first: +Division and death are decided inside the parallel phenotype loop, but a cell can't safely divide or remove itself from `(*all_cells)` while that loop is still running on other threads. So instead of acting immediately, a cell calls `flag_cell_for_division()` / `flag_cell_for_removal()`, which push it onto `cells_ready_to_divide` / `cells_ready_to_die` under `#pragma omp critical`. Only after the parallel loop finishes does [core/PhysiCell_cell_container.cpp](../core/PhysiCell_cell_container.cpp) actually call `divide()` / `die()` on those lists — serially, and, when `PhysiCell_settings.use_counter_based_rng` is enabled, only after sorting each list by cell ID first: ```cpp if ( PhysiCell_settings.use_counter_based_rng ) @@ -169,25 +185,25 @@ The `critical` block (where present) stops two threads from writing at the exact ### Fix 1: `standard_cell_cell_interactions` — attack, phagocytosis, fusion -[`core/PhysiCell_standard_models.cpp`](core/PhysiCell_standard_models.cpp) calls `pCell->ingest_cell(pTarget)`, `pCell->attack_cell(pTarget, dt)`, and `pCell->fuse_cell(pTarget)` directly from a per-cell loop. All three mutate the *target* cell's live state (volume, damage, death flag, position) and are wrapped in `#pragma omp critical` internally — but the guard that decides whether to act (e.g. `pCell_to_eat->phenotype.volume.total < 1e-15`) runs *before* that critical section. Two different cells can both decide to consume, attack, or fuse with the same target in the same step; which one's action actually lands first, and therefore what the target's final state is, depends on thread scheduling. +[`core/PhysiCell_standard_models.cpp`](../core/PhysiCell_standard_models.cpp) calls `pCell->ingest_cell(pTarget)`, `pCell->attack_cell(pTarget, dt)`, and `pCell->fuse_cell(pTarget)` directly from a per-cell loop. All three mutate the *target* cell's live state (volume, damage, death flag, position) and are wrapped in `#pragma omp critical` internally — but the guard that decides whether to act (e.g. `pCell_to_eat->phenotype.volume.total < 1e-15`) runs *before* that critical section. Two different cells can both decide to consume, attack, or fuse with the same target in the same step; which one's action actually lands first, and therefore what the target's final state is, depends on thread scheduling. -Fixed by wrapping the call to `standard_cell_cell_interactions()` in [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp) in `#pragma omp ordered` (loop changed to `schedule(static) ordered`). This also covers the function's early-exit check on `phenotype.death.dead`, which was itself an unsynchronized read of a flag another cell's thread could be writing at the same time. +Fixed by wrapping the call to `standard_cell_cell_interactions()` in [core/PhysiCell_cell_container.cpp](../core/PhysiCell_cell_container.cpp) in `#pragma omp ordered` — via `run_cell_cell_interactions_phase_ordered()`, taken only when `use_counter_based_rng` is enabled; see "Why these fixes are gated behind counter-based RNG" below. This also covers the function's early-exit check on `phenotype.death.dead`, which was itself an unsynchronized read of a flag another cell's thread could be writing at the same time. ### Fix 2: `dynamic_spring_attachments` — attachment-capacity race -[`core/PhysiCell_standard_models.cpp`](core/PhysiCell_standard_models.cpp) checks `pTest->state.spring_attachments.size() < pTest->phenotype.mechanics.maximum_number_of_attachments` to decide whether a neighbor still has room for a new attachment, then calls `attach_cells_as_spring()`. The capacity check is never re-validated once the lock inside `attach_cells_as_spring()` is actually held. Two different cells can both see "room for one more" on the same neighbor at once and both attach, silently pushing that cell past its configured maximum, in a way that depends on scheduling. +[`core/PhysiCell_standard_models.cpp`](../core/PhysiCell_standard_models.cpp) checks `pTest->state.spring_attachments.size() < pTest->phenotype.mechanics.maximum_number_of_attachments` to decide whether a neighbor still has room for a new attachment, then calls `attach_cells_as_spring()`. The capacity check is never re-validated once the lock inside `attach_cells_as_spring()` is actually held. Two different cells can both see "room for one more" on the same neighbor at once and both attach, silently pushing that cell past its configured maximum, in a way that depends on scheduling. -Fixed the same way: the call to `dynamic_spring_attachments()` in [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp) is now inside `#pragma omp ordered`. +Fixed the same way: the call to `dynamic_spring_attachments()` in [core/PhysiCell_cell_container.cpp](../core/PhysiCell_cell_container.cpp) is now inside `#pragma omp ordered` — via `run_spring_attachment_phase_ordered()` — when `use_counter_based_rng` is enabled. ### Fix 3: `max_cell_interactive_distance_in_voxel` — stale/racy neighbor-search bound -Each mechanics voxel keeps a running maximum of `radius * relative_maximum_adhesion_distance` over the cells that have been recorded there, used to decide whether a neighbor search needs to look inside that voxel at all. It's updated with an unprotected read-compare-write (`if (current < new) current = new;`) in `Cell::convert_to_cell_definition()` (fires on a rule-triggered type transformation) and, matching upstream PR [#409](https://github.com/MathCancer/PhysiCell/pull/409), in `Cell::update_voxel_in_container()` (fires on ordinary movement between voxels — without this, the bound goes stale as cells migrate, and a voxel's neighbor search can wrongly skip a cell that has since moved in). Both live in [core/PhysiCell_cell.cpp](core/PhysiCell_cell.cpp). +Each mechanics voxel keeps a running maximum of `radius * relative_maximum_adhesion_distance` over the cells that have been recorded there, used to decide whether a neighbor search needs to look inside that voxel at all. It's updated with an unprotected read-compare-write (`if (current < new) current = new;`) in `Cell::convert_to_cell_definition()` (fires on a rule-triggered type transformation) and, matching upstream PR [#409](https://github.com/MathCancer/PhysiCell/pull/409), in `Cell::update_voxel_in_container()` (fires on ordinary movement between voxels — without this, the bound goes stale as cells migrate, and a voxel's neighbor search can wrongly skip a cell that has since moved in). Both live in [core/PhysiCell_cell.cpp](../core/PhysiCell_cell.cpp). Unlike Fixes 1 and 2, this doesn't need `ordered`: a `max` over a fixed set of candidate values is the same regardless of the order you compare them in, so a plain `#pragma omp critical` around the read-compare-write is enough to make it both race-free and deterministic. ### Fix 4: cell secretion/uptake into shared diffusion voxels -PhysiCell's real per-cell secretion path is `Cell_Container::update_all_cells()`'s first loop → `Secretion::advance()` → `Basic_Agent::simulate_secretion_and_uptake()` (in [BioFVM/BioFVM_basic_agent.cpp](BioFVM/BioFVM_basic_agent.cpp)), which does: +PhysiCell's real per-cell secretion path is `Cell_Container::update_all_cells()`'s first loop → `Secretion::advance()` → `Basic_Agent::simulate_secretion_and_uptake()` (in [BioFVM/BioFVM_basic_agent.cpp](../BioFVM/BioFVM_basic_agent.cpp)), which does: ```cpp (*pS)(current_voxel_index) += cell_source_sink_solver_temp1; @@ -197,11 +213,23 @@ PhysiCell's real per-cell secretion path is `Cell_Container::update_all_cells()` directly on the microenvironment's shared voxel density vector, with no lock at all. Diffusion voxels are typically larger than a single cell, so it's normal for several cells to share one; when two cells in the same voxel are processed by different threads at the same time, this isn't just reordering, it's a genuine lost update — one thread's contribution can be silently overwritten instead of added. -Fixed by wrapping the call to `phenotype.secretion.advance(...)` in [core/PhysiCell_cell_container.cpp](core/PhysiCell_cell_container.cpp) in `#pragma omp ordered`. +Fixed by wrapping the call to `phenotype.secretion.advance(...)` in [core/PhysiCell_cell_container.cpp](../core/PhysiCell_cell_container.cpp) in `#pragma omp ordered` — via `run_secretion_phase_ordered()` — when `use_counter_based_rng` is enabled. + +### Why these fixes are gated behind counter-based RNG + +`#pragma omp ordered` isn't free: it forces the loop's marked region to execute in strict sequence across *every* iteration, not just the ones that actually conflict, which caps how much of that loop can ever be parallelized (this is "DOACROSS" parallelism in the literature, as opposed to the fully-independent "DOALL" case a plain `parallel for` handles). Benchmarks on comparable loop-carried-dependency patterns show real, sometimes large, slowdowns from coarse-grained `ordered` use, so applying it unconditionally to every PhysiCell run — including the large fraction of existing projects that have never asked for cross-thread-count reproducibility — was not an acceptable default. + +So Fixes 1, 2, and 4 are gated at their call site in `Cell_Container::update_all_cells()`, not inside a shared function: each is an `if( use_counter_based_rng ) { ...call a small ordered-only function... } else { ...original inline loop... }` written directly where the phase runs. The `if` branch calls one of three small, ordered-only functions added to the anonymous namespace in [core/PhysiCell_cell_container.cpp](../core/PhysiCell_cell_container.cpp): `run_secretion_phase_ordered()`, `run_spring_attachment_phase_ordered()`, and `run_cell_cell_interactions_phase_ordered()`. `use_counter_based_rng` is checked once per phase per step, not once per cell. + +The `else` branch of each does *not* call `activate_random_context()`/`clear_random_context()`, even though the equivalent (non-split) loops elsewhere in `update_all_cells()` do. That's not an inconsistency: those other loops keep the guarded call so one loop body works for both RNG modes, avoiding a needless second copy. But inside an `else` block that only runs when `use_counter_based_rng` is already `false`, that same guarded call is provably dead code — the internal `if( use_counter_based_rng )` check inside it can never be true there, so calling it can never do anything. Once the racy phases already need two separate loop bodies (for the ordered-vs-not structural difference), there's no cost to leaving those calls out of the `else` branch specifically, and doing so makes that branch true original behavior — no context bookkeeping at all — rather than "calls that happen to no-op." + +This keeps the diff against the pre-fix code additive with one small subtraction: the three new ordered functions, each racy loop wrapped in an `if`/`else`, and the now-pointless `activate_random_context()`/`clear_random_context()` calls dropped from the three `else` branches. Nothing else in `update_all_cells()` — intracellular update, phenotype update, division/death, `evaluate_interactions`, `custom_cell_rule`, `update_velocity`, spring-force application, position update — was touched. + +This is a deliberate, explicit trade-off, not an oversight: **with `counter_based_rng` disabled (the default), the three races described in Fixes 1, 2, and 4 are still present.** They are not reproducibility-only concerns — the secretion race in particular is a genuine lost-update on shared substrate mass, independent of RNG mode entirely, and predates the counter-based RNG feature. Turning on `counter_based_rng` is currently the only way to get both deterministic RNG *and* these race-free interaction/attachment/secretion phases. Fix 3 (`max_cell_interactive_distance_in_voxel`) is not gated — a `critical`-protected `max` update is cheap enough that there was no reason to make it conditional. ## How to enable counter-based RNG -Counter-based RNG is opt-in and **off by default**. If `` is not set at all, the simulation uses the legacy thread-local generator (`PhysiCell_settings.use_counter_based_rng = false` in [modules/PhysiCell_settings.h](modules/PhysiCell_settings.h)), which is not guaranteed to be reproducible across different thread counts. +Counter-based RNG is opt-in and **off by default**. If `` is not set at all, the simulation uses the legacy thread-local generator (`PhysiCell_settings.use_counter_based_rng = false` in [modules/PhysiCell_settings.h](../modules/PhysiCell_settings.h)), which is not guaranteed to be reproducible across different thread counts. To enable it, add (or set) an `` element inside `` in the config XML: @@ -212,7 +240,7 @@ To enable it, add (or set) an `` element inside `` in the con ``` -Accepted values, parsed in [modules/PhysiCell_settings.cpp](modules/PhysiCell_settings.cpp): +Accepted values, parsed in [modules/PhysiCell_settings.cpp](../modules/PhysiCell_settings.cpp): - `counter_based`, `counter`, or `philox` → enables the deterministic Philox-based path (`use_counter_based_rng = true`). - `legacy`, `thread_local`, or `mt19937` → explicitly selects the old, non-deterministic thread-local generator (`use_counter_based_rng = false`). @@ -222,7 +250,7 @@ An unrecognized value prints an error and leaves the setting unchanged. ## Testing a new model for thread reproducibility -[beta/test_thread_repro.py](beta/test_thread_repro.py) runs an already-built executable twice, once per thread count, and diffs the output. It does not build the project for you — build first (`make`), then run: +[beta/test_thread_repro.py](../beta/test_thread_repro.py) runs an already-built executable twice, once per thread count, and diffs the output. It does not build the project for you — build first (`make`), then run: ```bash python beta/test_thread_repro.py @@ -241,7 +269,7 @@ python beta/test_thread_repro.py project config/PhysiCell_settings.xml 1440 1 4 - `rng_mode`: written into each run's copy of `` before running (omit, or pass `off`/`none`/`default`, to leave the config's existing setting untouched). - `work_dir`: where the two runs' output is written (each invocation creates a timestamped subfolder inside it). -The script runs both thread counts, then diffs the final XML/SVG snapshots with [beta/test_diff_svg.py](beta/test_diff_svg.py) and reports the first mismatch, if any. +The script runs both thread counts, then diffs the final XML/SVG snapshots with [beta/test_diff_svg.py](../beta/test_diff_svg.py) and reports the first mismatch, if any. Because these races can be intermittent — several of the fixes above only failed some fraction of the time — a single passing run is not strong evidence on its own. Repeat the same command several times (5–10x) before trusting a "pass," especially for a new or modified model. @@ -251,4 +279,6 @@ The new system tests pass with 1 thread and 4 threads for the template project. The four cross-cell-mutation fixes above were found and validated by testing the template project with multiple cell behaviors activated (phagocytosis, attack, transform, and others) using `beta/test_thread_repro.py` across single and multiple threads. Repeated trials failed consistently before Fix 4 (the secretion race) was in place; after all four fixes were applied, repeated trials passed bit-exact. -Known open item: a project-level (not core-engine) bug remains in the worm sample's custom `contact_function` ([custom_modules/custom.cpp](custom_modules/custom.cpp) when the worm project is loaded), which writes into an attached cell's `custom_data` from within a parallel loop without going through the engine's RNG/ordering machinery. That is a project-code issue rather than a PhysiCell core issue and is tracked separately. +After gating Fixes 1, 2, and 4 behind `use_counter_based_rng`, both branches of each rewritten phase were re-checked: with `rng_mode=counter_based`, repeated trials still pass bit-exact; with `rng_mode=legacy`, both thread counts still run to completion without errors (confirming the unordered branch itself is correct), and are not expected to, and do not, match each other — legacy mode was never meant to guarantee cross-thread-count reproducibility, gated fixes or not. + +Known open item: a project-level (not core-engine) bug remains in the worm sample's custom `contact_function` ([custom_modules/custom.cpp](../custom_modules/custom.cpp) when the worm project is loaded), which writes into an attached cell's `custom_data` from within a parallel loop without going through the engine's RNG/ordering machinery. That is a project-code issue rather than a PhysiCell core issue and is tracked separately.