Conversation
also make substrates.csv in template project. disable in xml
…-and-ecm-ics hide argument_parser in core
avoid dividing by cell volume
- about midway between 1.14.0 and 1.14.1
- need to go back and decide how to handle the "old" rules - right now, just ignoring them (and breaking old projects in the process)
- also some cleaning - make sure to look at previous commit for notes
- only need a signal reference for AbstractHill signals - make these fall into new class RelativeSignals - others fall into AbsoluteSignals - example showing off many of these new rules
- these allow for memory-ful behaviors where signals lead to either accumulation or attenuation of the behavior over time - basically, if a behavior is ruled, then they it can either be set, accumulated, or attenuated - set is the default behavior and what has happened before - accumulation and attenuation will move away from the base behavior if they process a positive signal - otherwise they relax to the base - the top-level signal for these behaviors are still mediators - but now the mediator sets the rate of accumulation or attenuation - a positive rate means the behavior is moving away from the base - a negative rate means the behavior is relaxing to the base - users will likely want the base behavior to be negative so that in the absence of increasing signals the behavior will relax to the base
…athCancer#60) get_parameter_value() and set_parameter_value() index species_result_column_index with operator[], the inserting form. An unknown species name silently resolves to column 0, so a typo -- or a model whose column map was never populated -- reads, and writes, a different species than the caller named. A read also mutates the map. Both accessors now look the species up with find() and, on a miss, report and exit. Exiting rather than warning is deliberate: there is no value to return and nothing sensible to write, so continuing would feed a fabricated number into the model and every step after it. validate_SBML_species() already exits for exactly this condition at setup; it just cannot cover these calls, because it only validates the mappings declared in the XML and custom code calls the accessors directly with its own strings. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
The block tested phenotype.cell_interactions.pAttackTarget for NULL and then loaded it again into a local. Those are two separate loads of a field another thread can clear: remove_all_attackers() nulls an attacker's pAttackTarget when its target is eaten, fused or lysed, and that runs from standard_cell_cell_interactions inside the mechanics parallel-for. If the clear lands between the test and the reload, the local is NULL and attack_cell() dereferences it -- it only guards against attacking itself. Take one snapshot and use it throughout. Follow-up to MathCancer#59, which introduced the cross-thread write. Mirror of MathCancer#422's 6c79881. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
ingest_cell, fuse_cell and lyse_cell each flag their victim for removal and then immediately tear down its attachments, springs and attack links, from inside the mechanics parallel-for. That mutates other cells' state from a thread that does not own them. remove_all_spring_attachments walks state.spring_attachments by index with no lock, while attach_cell_as_spring and detach_cell_as_spring edit that same vector under the unnamed critical. detach is swap-and-pop, so a removal below the walker's cursor moves the tail element into a slot already passed and it never gets detached -- leaving a neighbour holding a pointer to a cell about to be freed. The work was already redundant: flag_for_removal() queues the cell and the serial cells_ready_to_die -> die() -> delete_cell() pass does the same steps. Dropping the in-loop copies removes every unlocked walk from a parallel region. Mirror of the same commit on MathCancer#422. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
…meters (MathCancer#56) * Give each cell definition its own fixed_duration and death-phase parameters Fixes MathCancer#199. Cycle_Model objects are shared by pointer between cell definitions: apoptosis and necrosis are globals, and every definition registers the same address via add_death_model( rate, &apoptosis, ... ). Cell_Definition's copy constructor copies the pointer, not the model. fixed_duration lived in Phase_Link, inside that shared object, while the XML parser wrote it per cell definition. So each definition overwrote the flag for every other one and the last definition parsed decided for all of them. A model asking for deterministic apoptosis on most of its cell types silently got the stochastic branch on all of them, with no warning and nothing in the saved output to reveal it. Migrate the flag from the graph structure to the param structure, as Paul suggested on the issue: Cycle_Data gains fixed_durations, indexed exactly like transition_rates, with fixed_duration() / exit_fixed_duration() accessors mirroring transition_rate() / exit_rate(). Cycle_Model::advance_model now reads the per-cell copy. Phase_Link::fixed_duration is removed; the six standard models declare their defaults through the model's own Cycle_Data instead, and those defaults still reach cells that do not override them in XML. That alone is not enough for the death models. There is no per-definition Cycle_Data for a death model -- Death holds rates, models and parameters, none of which carry cycle parameters -- so the parser wrote durations to the shared models[i]->data, and Cycle::sync_to_cycle_model() copies that same shared data over the cell at start_death(). Relocating the flag would have moved it from one shared object to another. This is also why rheiland observed on the issue that the <duration> value is flattened along with the flag. So Death gains model_data, one Cycle_Data per death model, seeded from the model at add_death_model(), written by the parser, and applied at start_death() through a new two-argument Cycle::sync_to_cycle_model( cm, cd ). The single-argument form delegates to it and keeps its exact previous behaviour. Verified on a 7-cell-type model where six definitions ask for <phase_durations fixed_duration="true"> and the last asks for a rate. Before: all seven report fixed=0. After: six report 1 and the last reports 0, while all seven still share one Cycle_Model object, so the phase graph stays shared and Phenotype copies stay cheap. A definition with no apoptosis <phase_durations> at all still inherits the standard model's default. Note: this changes results for any model that mixes <phase_durations> and <phase_transition_rates> across its cell definitions. Those models were not getting the death timing their XML asked for; they will now. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Apply the per-definition death parameters on the runtime death path too Cell::advance_bundled_phenotype_functions syncs the cycle to the death model directly when check_for_death fires, rather than going through start_death. That site was still using the single-argument sync_to_cycle_model, so cells dying by death_rate -- the common case -- kept getting parameters from the shared Cycle_Model and MathCancer#199 was unfixed on that path. Also from review: hoist the transition_rates / fixed_durations resize out of the inner phase-link loop in Cycle_Data::sync_to_cycle_model, and take the Cycle_Data by const reference in the new Cycle overload. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Report per-definition death parameters in the cell definition summary display_cell_definitions() read the death phase durations out of death.models[k]->data, the shared Cycle_Model. The parser no longer writes that object, so after the rest of this PR the summary printed the compiled-in standard-model defaults for every definition regardless of the XML -- and that summary is exactly what a user reads to check whether MathCancer#199 is fixed. Take the parameters from death.model_data[k] instead. Phases and phase links still come from the shared model, which is correct: only the parameters moved. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
…ancer#61) Every object rule names only its .cpp, so make never rebuilds an object when a header it includes changes. Editing a header and running make relinks the stale objects and reports success. That is merely confusing when the change is behavioural. It is dangerous when the header changes a type's layout: the objects that did get rebuilt disagree with the ones that did not about member offsets and sizeof, the link succeeds because mangled names encode types rather than layouts, and the result is a binary that reads fields at the wrong addresses. The crash looks like a bug in the code under test. -MMD makes the compiler emit a .d file beside each .o listing the headers that object actually included; -MP adds a dummy target for each so that deleting or renaming a header does not break the build. The .d files are read back in at the bottom of each Makefile, which is where they have to be: each one declares rules for its own .o, so including them earlier would replace `all` as the default goal. Applied to all 44 Makefiles, clean removes the generated .d files, and they are gitignored. Mirror of MathCancer#426 for my-physicell. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
…Cancer#64) asymmetric_division_function tolerates a small excess when checking that the probabilities sum to at most 1, because probabilities meant to sum to exactly 1 land a hair over it in double precision. That tolerance was hard-coded at 1e-12, which is the right default but not right for every model: how much slack is needed depends on how many probabilities are summed and how they are produced. A model whose probabilities come out of rules rather than literals can exceed 1 by far more than 1e-12 -- excesses around 8e-3 have been observed in practice -- and such a model currently exits with no way to say that is acceptable. Adds <asymmetric_division_probability_tolerance> to the <options> block. Absent, the default 1e-12 applies and nothing is printed, so existing configs are unaffected. Present, the value is used and echoed at startup. A negative value is rejected at parse time rather than silently inverting both comparisons. The name says which quantity is toleranced: it is slack on the probability SUM, not on division timing or on the division itself. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
…cer#65) dynamic_spring_attachments() runs inside the mechanics parallel-for and walks pCell->state.spring_attachments with no lock. That vector is not private to this thread: another thread running the same function for a neighbouring cell calls attach_cells_as_spring(), which reaches into THIS cell's spring_attachments and push_back()s. Those writers serialize on the unnamed critical inside Cell::attach_cell_as_spring() and Cell::detach_cell_as_spring(); the read did not take it. A concurrent push_back that reallocates leaves the loop indexing a freed buffer, and the next iteration dereferences whatever it finds. The result is an intermittent SIGSEGV in dynamic_spring_attachments, inside the OpenMP outlined region, faulting on a garbage address. It is thread-count dependent and therefore invisible on single-threaded runs, which is why it survives casual testing. Copy the vector under the same critical the writers use, then walk the copy. The lock is deliberately not held across the loop body: detach_cells_as_spring() takes the same unnamed critical, and OpenMP criticals are not reentrant. Measured on a model that reproduces the fault reliably, 60 runs per arm at 4 threads: 10 crashes without the change, 0 with it. Cost is +2.4% wall time (paired median ratio 1.0237, slower in 9 of 9 paired blocks). This addresses the crash only. The function still reads .size() outside the lock in the attachment half, so it is not yet data-race-free. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
* Adjust asymmetric division probabilities whenever they sum past 1 asymmetric_division_function entered its renormalization block on total > 1.0 + tolerance, so raising the tolerance did not only relax when an overshoot is an error -- it also decided whether the probabilities were rewritten at all. That block sets the symmetric division probability, so gating it on a user-settable threshold let the tolerance change the model rather than its strictness. With a tolerance of 0.5 and probabilities of 0.3 for (stem, stem) and 1.0 for (stem, progenitor_1), the total of 1.3 fell inside the old gate, so nothing was adjusted and the draw ran against a distribution summing to 1.3. Because select_daughter_types draws from [0,1), (stem, stem) kept its 0.3 instead of being taken down to 0, and progenitor_1 was starved to 0.7. Running the asymmetric division sample that way ends with 464 cells, 137 of them stem; adjusting whenever the total passes 1 ends with 16 cells and a stem population that stays at 1, which is what those probabilities describe. Adjust on total > 1.0 and leave the tolerance its one job: deciding whether the resulting symmetric division probability is negative enough to be an error rather than round-off. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Add a weighted mode for asymmetric division The extended asymmetric division values are probabilities: they must sum to at most 1, and whatever they leave short of 1 is the chance of dividing symmetrically. That constraint is awkward for the models this feature is aimed at. The intended use is rule-driven values, and a Hypothesis Grammar rule cannot see what the other daughter-type pairs currently evaluate to, so keeping a sum at or below 1 is not something a rule set can enforce. The configurable tolerance buys slack for round-off, but cannot rescue a model whose values legitimately sum to 2 or 3. Adds <asymmetric_division_mode> to the <options> block, accepting "probabilities" (the default, and what an absent element means) or "weights". Under weights the values are relative weights, normalized by their own sum at each division, so their scale stops mattering and each rule can be written on its own. A string rather than a bool so a typo errors out instead of quietly reading as false and running the other model. The mode is model-wide rather than per cell definition: mixing them would make one rule line mean a probability for one source type and a weight for another, with nothing in the line to say which. Weights carry no implicit symmetric-division remainder -- symmetric division needs its own (type,type) weight -- except that an all-zero total, which has no normalized distribution, is defined as symmetric division. That needs no special case: leaving the draw scale at 1.0 makes select_daughter_types fall through to the parent and daughter types unchanged. The draw itself is one line. select_daughter_types takes a total_weight that scales the random draw; probability mode passes 1.0, which is exact in double precision, so no existing model's RNG stream moves. Confirmed by running the asymmetric division sample before and after: 487 .mat files byte-identical, and 244 SVG and 243 XML identical once wall-clock timestamps are stripped. Setting <asymmetric_division_probability_tolerance> together with weights is rejected at parse time. The tolerance bounds how far probabilities may sum past 1, and weights are normalized rather than bounded, so there is nothing for it to do; erroring beats letting one of the two silently win. An explicit "probabilities" mode alongside a tolerance stays legal, since that is unambiguous. The extended_asym_div sample gains a second config and rules file expressing the same model in weights, plus a README covering both modes and the trap that separates them: probabilities take the whole overshoot out of the symmetric entry, while weights scale every entry proportionally. The weights rules saturate at 1.0 where the probabilities rules use 0.5, chosen so the normalized weights reproduce the probabilities exactly. Verified single-threaded: across all 121 snapshots the two runs agree on every cell, and the only recorded difference is the asymmetric division values. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
…ncer#73) get_row_from_dirichlet_condition_csv took line and substrate_indices by value, copying a string and a vector for every row of the file. Take both by const reference, matching the substrate reader. The "wrong number of density data" warning printed number_of_voxels() where it means number_of_densities(); Found is a count of substrate columns, and the matlab loader's equivalent message already uses the density count. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
…quire well-formed rows (MathCancer#66) (MathCancer#69) * Bounds-check substrate IC csv rows before indexing them get_row_from_substrate_initial_condition_csv indexed data[0..2] and data[ci + 3] without consulting data.size(). Both reads are out of bounds for a malformed row: - substrate_csv_to_vector always emits a final field, so a blank line parses to a single 0.0 and data[1]/data[2] read past the end. - substrate_indices is sized from the header's column count, so a row with fewer columns than the header runs data[ci + 3] off the end. Reproduced with libc++ bounds checking on a 4-voxel microenvironment: a trailing blank line, a whitespace-only line, a CRLF blank line, a 2-column row, and a row supplying 1 of 3 header substrates all abort with "vector[] index out of bounds" before this change and are handled cleanly after it. Note these reads are not new here -- development crashes on the same five inputs. Its coverage check ran after the read loop, so it never guarded them. Blank lines are now skipped, matching load_cells_csv_v1 and process_csv_v2_line, and a short row is a hard error rather than a silent partial write. Ordering the guards ahead of the existing warning also stops data.size() - 3 from underflowing on a short row. Also correct that warning: it reported number_of_voxels() where it means number_of_densities(), and the block comment still described the one-row-per-voxel, no-header format this branch replaced. * Revive the headerless csv path and require well-formed rows The headerless branch of load_initial_conditions_from_csv could never have worked. Two bugs had to be fixed together to see either one: - "if (i<3) {continue;}" jumped past the "i++" below it, so i never advanced and substrate_indices came out empty. - the reopened "std::ifstream file(filename, ...)" shadowed the enclosing stream, which had just been closed, so the row loop read a closed stream and processed no rows. -Wshadow flags this. A headerless csv therefore loaded nothing. On development that at least failed loudly, because the deleted coverage check saw voxel_set.size() == 0 and exited; without it the run continued silently on the config file's uniform initial conditions. Reading the first row's column count and rewinding the stream fixes both, and the row loop now counts lines so every diagnostic can name the row it rejected. Rows are now held to being well formed rather than parsed as far as they go: - a field must be empty or a complete finite number. strtod's endptr was discarded, so "NA" and "1.5abc" silently became 0 and "inf" became a density. - a row must have exactly the column count the header (or the first row) established. Extra columns were silently dropped. - x, y and z must all be present, so ",,,," no longer resolves to the origin. - a position must lie inside the domain, since nearest_voxel_index clamps and would otherwise snap a typo onto an edge voxel. - a header may not name the same substrate twice, and must name at least one. - an empty file is an error rather than a header sniff on nothing. An omitted entry now travels as NaN instead of 0, which is what lets the row loop tell "the user left this blank" from "the user asked for zero". It still resolves to 0, as before; the distinction only makes the checks above possible. The header sniff also no longer indexes line.c_str()[2] and [4] without knowing they exist -- a first line of "x" read past the end. It splits the row and compares fields instead, which also makes the sniff whitespace tolerant. Verified against real BioFVM on a 2x2x1 microenvironment with libc++ bounds checking: 36 cases covering header/headerless, subsets, reordering, whitespace, CRLF, blank lines, and every rejection above. All 36 behave as intended with no out-of-bounds reads. Headerless files load correctly for the first time. --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
…er#74) template_BM and tutorial both wrote if( !load_PhysiCell_config_file(); ) so neither sample compiled at all: main.cpp:92:44: error: expected primary-expression before ')' token g++ also warned about the accidental init-statement in a selection statement (-Wc++17-extensions), which goes away with the semicolon. cancer_invasion already has the intended form, and is what these two are matched to. Introduced in 7eba010 when these samples moved to the argument_parser plus load_PhysiCell_config_file() API. Not present on upstream development, which still uses the older XML_status form, so this is fork-only and needs no upstream PR. These are the only two occurrences in the repository. This is why "Testing PhysiBoSS Template" and "Testing PhysiBoSS Tutorial" failed on Ubuntu, Windows and MacOS 14 on every PR, along with the windows/macos_step0b template_BM build jobs. Both samples now build, link and run: template_BM completes with the shipped config, and tutorial completes against config/cell_cycle/PhysiCell_settings.xml (it ships no top-level config, only the three scenario subdirectories). Note this does not address the other failure in those jobs, "Makefile.maboss:193: BooleanNetwork_256n.o Error 127", which is a missing command in the MaBoSS build environment rather than anything in this tree. MaBoSS built here without trouble, so that one is CI-specific. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
MathCancer#69 gave the substrate reader a parser that returns the 1-based index of the first malformed field and yields NaN for an omitted entry. That NaN sentinel is exactly what the parallel is_missing vector existed to carry, so the dirichlet reader can use the same parser and both dirichlet_csv_to_vector and the store_dirichlet_csv_entry helper of the previous revision of this branch go away. Rather than copy the header handling a third time, the two loaders now share it: resolve_substrate_csv_columns reads the optional "x,y,z,<name>,..." header, applies the headerless "first n densities" convention, and leaves the stream on the first data row. The two readers differ only in what a parsed row does to the microenvironment -- the substrate reader resolves an omitted entry to 0, the dirichlet reader leaves that voxel-substrate pair alone -- and in the word their diagnostics use. Sharing it also carries MathCancer#69's fixes onto the dirichlet path, which had the same two bugs its substrate counterpart did: - "if (i<3) {continue;}" jumped past the "i++" below it, so substrate_indices came out empty on a headerless file. - the reopened "std::ifstream file(filename, ...)" shadowed the enclosing stream, which had just been closed, so the row loop read a closed stream. A headerless dcs.csv therefore set no dirichlet conditions at all and said nothing about it. Verified on the dirichlet_from_file sample: before this change a headerless dcs.csv produced a byte-identical final microenvironment to disabling the file entirely; after it, it matches the headered file. Blank lines are skipped rather than read as a row of zeroes, a row whose column count disagrees with the header is a hard error, and every diagnostic names the line and column it rejected. The out-of-bounds write in dirichlet_csv_to_vector, which had no bound on ind, disappears with the function. The stock sample is unaffected: its dcs.csv leans on empty fields meaning "leave this pair alone" on nearly every row, and it produces byte-identical initial and final microenvironments before and after. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
…r#75) * Seed the RNG from OS entropy and widen the seed to 64 bits `SeedRandom(void)` took the seed straight from `system_clock::now().time_since_epoch().count()`, truncated into an `unsigned int`. The clock advances in ~1 microsecond steps (measured with homebrew g++-15 on macOS: 200k back-to-back `now()` calls returned only 10059 distinct counts), so processes launched together routinely read the same count and share a seed. Running 8 template sims simultaneously, 2 of 40 batches contained a duplicate `random_seed.txt`. `SeedRandom(void)` now mixes `std::random_device` with the clock and the process id. The clock and pid keep the seed varying on toolchains where `random_device` is a deterministic stub (older MinGW), and the pid separates processes that read the clock in the same tick. The seed chain is also widened to 64 bits: `physicell_random_seed`, `physicell_random_seeds`, `SeedRandom`'s argument, and the config parse (`std::stoull`). The generator is `std::mt19937_64`, so the old 32-bit seed both threw away entropy and silently truncated configured seeds at or above 2^32 (`<random_seed>4294967296</random_seed>` was recorded as `0`). `std::seed_seq` consumes 32-bit words, so per-thread seeds derived from a seed above 2^32 now contribute a second (high) word. Seeds that fit in 32 bits still contribute one word each, exactly as before, so fixed-seed runs are bit-identical to the previous behavior -- verified for seeds 0 and 42 at 1 and 6 threads against a binary built from the parent commit. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Derive per-thread RNG seeds with splitmix64 and seed every thread that draws Follow-up to the 64-bit seed change, dropping the requirement that fixed-seed runs reproduce previous releases. Three separate defects: 1. Per-thread seeds were 32-bit and thread 0's was the raw base seed. They came from `std::seed_seq`, whose input and output are both 32-bit words, so a 64-bit generator was started from a 32-bit value on every thread but one, and `<random_seed>0</random_seed>` started thread 0's mt19937_64 from 0. Seeds now come from splitmix64: thread i's seed is `mix64(base + (i+1) * 0x9E3779B97F4A7C15)`, the closed form of splitmix64's (i+1)-th output. Full 64 bits, a bijection so seeds never collide, ~32 of 64 bits flipped between adjacent base seeds or adjacent threads, and dependent only on (base seed, thread index) -- not on the thread count, as the `seed_seq` version was. 2. `NormalRandom()`, `UniformInt()` and `UniformRandom_old_not_thread_safe()` read `physicell_PRNG_generator` without the lazy seeding that `UniformRandom()` did. A thread whose first random number came from one of them drew from a default-constructed `std::mt19937_64` -- seed 5489 on every thread. Measured on the parent commit with 6 threads: `NormalRandom` and `UniformInt` returned one value shared by 5 of the 6 threads. This reached PhysiBoSS (`UniformInt()` seeds each MaBoSS engine, `LogNormalRandom` sets `time_to_update`) and PhysiMeSS fibre setup. All four entry points now seed first. 3. Re-seeding never reached worker threads. `local_pnrg_setup_done` was set once and never cleared, so a second `SeedRandom()` -- the episode sample project, or a `random_seed` user parameter -- left every worker thread on its old stream. Measured on the parent commit: `SeedRandom(42)` twice in one process gave different results the second time on 5 of 6 threads, so `random_seed.txt` did not actually pin the run. A generation counter now makes each thread reseed on its next draw. `SeedRandom(void)` absorbs its three entropy sources through `mix64` one at a time instead of a flat xor, so a small process id cannot cancel against the clock's low bits. Fixed-seed runs are reproducible run to run (1 and 6 threads, and after a re-seed), and no longer match previous releases -- that is the trade this commit makes deliberately. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Drop thread_local from the uniform distributions; fold the seed install together Two follow-ups from review. The per-call `seed_thread_rng_if_needed()` is not the cost it looks like: `UniformRandom()` already did exactly this check inline (`if( local_pnrg_setup_done == false )`), so the added cost there is a second integer read. Measured, best of 5 x 200M calls, 1 thread: before this branch after this commit UniformRandom() 10.10 ns 10.31 ns 6.96 ns UniformInt() 5.85 ns 10.83 ns 6.94 ns NormalRandom(0,1) 17.71 ns 19.65 ns 19.43 ns The real cost was the `thread_local` distribution objects, not the check. `std::uniform_real_distribution` and `std::uniform_int_distribution` hold only their parameters and their `reset()` is a no-op, so a plain local is free to construct and produces an identical stream -- verified on four fixed-seed configs. It also removes the shared `static` int distribution in `UniformInt()`, which every thread was mutating through `operator()`. `std::normal_distribution` really is stateful (it caches the second Box-Muller variate in `_M_saved`), so `NormalRandom()` keeps building one per call exactly as it always has. End to end, heterogeneity-sample at one thread, seed 0, best of 2: 76.86 s before this branch, 76.79 s with it, 74.13 s with this commit. Second, `install_thread_seeds()` now fills `physicell_random_seeds` and stamps the new generation in one place, so a new set of seeds and the stamp that makes threads pick it up cannot drift apart. Kept as a function rather than a struct because `physicell_random_seed` and `physicell_random_seeds` are `extern` in PhysiCell_utilities.h; moving them into a type would break user code that touches them. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Put the RNG seed state in a class and derive thread seeds on demand Replaces the loose seed globals with an RNG_State holding the base seed and the generation stamp, mutable only through reseed(). Three things this buys: - The base seed can no longer be set without reseeding. `physicell_random_seed` was a writable extern; assigning it compiled, changed nothing (the generators were already seeded and nothing told them to look again), and left random_seed.txt describing a run that did not happen. The same class of bug as the stale worker threads fixed two commits ago, now unrepresentable. - `physicell_random_seeds` is gone. Thread seeds are a closed-form function of (base seed, thread index), so caching them in a vector bought nothing and cost a bounds check plus a fallback path on the seeding route, for threads beyond the configured count. seed_for_thread() is defined for every index, so that case is no longer special. - The header exports two read-only accessors, get_random_seed() and get_thread_random_seed(int), instead of two mutable globals. Behavior is unchanged: identical output on four fixed-seed configs, and get_thread_random_seed(3) returns the same value the cached vector held. No measurable cost. Alternating rounds against the previous commit, best of 5 x 200M calls, one thread, UniformRandom(): 8.49 / 8.89 / 8.59 ns before, 8.84 / 8.79 / 8.47 ns after. Also corrects the numbers in the previous commit message, which compared runs made minutes apart under different machine load. Measured back to back, and stable across rounds, UniformRandom() is 12.3-12.9 ns before this PR and 8.5-8.8 ns after -- a ~32% improvement, the same ratio that message reported from less careful absolute numbers. UniformInt() goes 7.2 -> 8.2-8.8 ns and NormalRandom() 21.9-22.5 -> 23.7 ns; both gained the seeding check they should always have had. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Seed every thread when the seed is installed, not lazily on first draw reseed() now sets the base seed and seeds every thread's generator in one parallel region, so the link between the recorded seed and the streams threads draw from is a single visible statement. That removes the generation stamp, the thread_local copy of it, seed_thread_rng_if_needed(), and its four call sites. The lazy check existed to cover a thread whose index exceeds <omp_num_threads>. That was robustness this code never had: the original indexed `physicell_random_seeds[omp_get_thread_num()]` into a vector sized to the configured thread count, so such a thread read past the end -- undefined behavior, which is why it produced 7 of 8 and 14 of 16 distinct streams rather than failing cleanly. Nothing in PhysiCell, its addons, or the sample projects uses a num_threads clause or changes the thread count after startup; every main() calls omp_set_num_threads(PhysiCell_settings.omp_num_threads) once. The eager version gives 5 of 8 and 6 of 16 in that same test: still wrong, but now defined rather than an out-of-bounds read. Everything this PR is actually about is unaffected, because every thread is seeded before it can draw: NormalRandom() and UniformInt() give six distinct values across six threads on a first call, and reseeding to the same value twice reproduces. Output is identical to the previous commit on four fixed-seed configs. Best of 5 x 200M calls, one thread, alternating rounds: lazy check eager UniformRandom() 6.85-6.92 5.55 ns UniformInt() 6.85-6.89 5.30 ns NormalRandom(0,1) 19.13-19.24 17.35-17.51 ns End to end the difference does not clear run-to-run noise (heterogeneity-sample at one thread: 73.4-75.4 s with the check, 72.7-74.2 s without), so this is for the simpler code, with the speed as a bonus. 65 lines of code added over my-physicell, down from 74. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
…dels - mesenchymal_tumor_class1_class2 no longer returns to epithelial_tumor_class1: its base transformation rate to that type goes from 1e-2 to 0 in the HTAN config, the IMC template and all 48 IMC ROI configs. Class I+II now returns to epithelial only as class I+II, like every other class. - HTAN config: T-cell base speed 0.649298757545331 -> 0.6665 (CD4_Tcell, CD8_Tcell, Treg, CD8_exhausted), the value set with the corrected motility fit in eed649e, which updated the IMC config only. - IMC template and ROI configs: CD4_Tcell, CD8_Tcell and Treg base cycle rate 0 -> 0.00000351636, the HTAN value; CD8_exhausted stays 0 as in HTAN. Only these values change; every other byte of the 50 files is untouched. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Fix the class I+II leak and align T-cell speed and division across models
save/full_data/interval and save/SVG/interval go from 30 to 60 min in the HTAN config, the IMC template and all 48 IMC ROI configs. A 7-day run now writes 169 saves instead of 337, halving its output. Nothing else changes. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Save full data and SVG snapshots every 60 min instead of 30 save/full_data/interval and save/SVG/interval go from 30 to 60 min in the HTAN config, the IMC template and all 48 IMC ROI configs. A 7-day run now writes 169 saves instead of 337, halving its output. Nothing else changes. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
- user_projects/ root: a project saved without a name (Makefile, main.cpp, config/, custom_modules/, VERSION.txt). - config/ics: make-load copies and older IC generations. PCMM passes every run's ICs with -i, and the IMC scripts now read the tracked copy under user_projects/antigen_presentation/config/ics. - HTAN project config: copies of sample-project files, other studies' ICs and settings, and a config/Makefile copy. - IMC project config: rule variants, the old ics/substrate folder, Xenium ICs, two JHH387ROI1 duct-filler variants, and config/Makefile and VERSION.txt copies. Nothing in the runs or the prep scripts reads these files: regenerating the IMC run inputs gives byte-identical output. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
- Replace absolute laptop paths with paths relative to PhysiCell/, the convention the 48 IMC ROI configs already use. - Set omp_num_threads to 1 in both templates and the 48 ROI configs, as in the run configs. - IMC template: name its output folder outputs/template and note that it is the stage-3 template (generate_roi_configs.py overwrites the domain, output folder, ICs and rules folder per ROI; its volumes are the global defaults). - HTAN template: output folder output, an IC file that exists in the project, no reference to the disabled ECM file, and the spatial_config_index user parameter from the run config. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
- HTAN project: main.cpp, Makefile and custom.h now match data/inputs/custom_codes, so it builds the extended rules engine with the argument parser, as the runs do. - custom.cpp (both projects): keep initialize_default_cell_definition(), which PhysiCell needs, correct its comment, and move the secretion sync that sat after the comment onto its own commented-out line. - custom.cpp and custom.h: comment out setup_tissue_domain(), which read the domain bounds into locals and never used them. - VERSION.txt: 1.14.2-drbergman-2.5.1, as in PhysiCell/VERSION.txt. A build from these files is byte-identical to one from the run custom code, and 6-hour test runs give byte-identical outputs before and after. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Cleanup/user projects
…ple projects Merge the 1.14.2-drbergman-3.2.1 release. Everything outside user_projects/ and .gitignore now matches the release exactly, apart from the sample projects: - Engine (core, modules, BioFVM), root Makefile and VERSION.txt from 3.2.1: per-thread seeds mixed with the process id, per-cell-type death-phase flags, a symmetric attack link, thread-safety fixes and stricter CSV parsing. - addons/ as released, restoring PhysiECM and PhysiPKPD, which the 2.5.1 merge left out; unit_tests, tests and beta take the release versions. - sample_projects, sample_projects_intracellular and sample_projects_physipkpd removed: the model and PCMM runs do not use them. - user_projects: VERSION.txt set to 1.14.2-drbergman-3.2.1. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Update/drbergman 3.2.1
Eight rows in both user projects' cell_rules.csv: a tumour cell touching a CAF or apCAF gains class II (class1 -> class1_class2, no class -> class2) in both lineages. Their rate is 0, so runs are unchanged; the main repo's runners switch them on with a rules variation. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Halves the output per run (about 0.67 TB instead of 1.3 TB for a full set of 534 runs), to leave room for replicates. Applied to both user-project templates and the 48 IMC ROI configs, which setup_imc_spatial_pcmm.py turns into the IMC run config. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
4 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Add CAF-contact MHC-II rules (off by default) and save every 2 hours
Checks:
🤖 Generated with Claude Code