Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

refactor update_all_cells #338

Closed
wants to merge 5 commits into from
Closed

refactor update_all_cells #338

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
7d53c63
refactor update_all_cells
drbergman Nov 24, 2024
2e8fb1f
fix logic in timing of phenotype/mechanics updates
drbergman Nov 26, 2024
7ac3e39
remove change to detach
drbergman Nov 27, 2024
9abcf80
edit whitespaces
drbergman Nov 27, 2024
686f16b
rename var to time_since_last_phenotype
drbergman Dec 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
91 changes: 38 additions & 53 deletions core/PhysiCell_cell_container.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,14 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me
{
// secretions and uptakes. Syncing with BioFVM is automated.

static double time_since_last_phenotype = phenotype_dt_;
static double time_since_last_mechanics = mechanics_dt_;
static double phenotype_threshold = phenotype_dt_ - 0.5 * diffusion_dt_;
static double mechanics_threshold = mechanics_dt_ - 0.5 * diffusion_dt_;

bool time_for_phenotype = time_since_last_phenotype > phenotype_threshold;
bool time_for_mechanics = time_since_last_mechanics > mechanics_threshold;

#pragma omp parallel for
for( int i=0; i < (*all_cells).size(); i++ )
{
Expand All @@ -133,18 +141,12 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me
}
}

//if it is the time for running cell cycle, do it!
double time_since_last_cycle= t- last_cell_cycle_time;

static double phenotype_dt_tolerance = 0.001 * phenotype_dt_;
static double mechanics_dt_tolerance = 0.001 * mechanics_dt_;

// intracellular update. called for every diffusion_dt, but actually depends on the intracellular_dt of each cell (as it can be noisy)

#pragma omp parallel for
for( int i=0; i < (*all_cells).size(); i++ )
{
if( (*all_cells)[i]->is_out_of_domain == false && initialzed ) {
if( (*all_cells)[i]->is_out_of_domain == false ) {

if( (*all_cells)[i]->phenotype.intracellular != NULL && (*all_cells)[i]->phenotype.intracellular->need_update())
{
Expand All @@ -159,54 +161,58 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me
}
}

if( time_since_last_cycle > phenotype_dt_ - 0.5 * diffusion_dt_ || !initialzed )
if( time_for_phenotype )
{
// Reset the max_radius in each voxel. It will be filled in set_total_volume
// It might be better if we calculate it before mechanics each time
// std::fill(max_cell_interactive_distance_in_voxel.begin(), max_cell_interactive_distance_in_voxel.end(), 0.0);

if(!initialzed)
{
time_since_last_cycle = phenotype_dt_;
}

// new as of 1.2.1 -- bundles cell phenotype parameter update, volume update, geometry update,
// checking for death, and advancing the cell cycle. Not motility, though. (that's in mechanics)
#pragma omp parallel for
for( int i=0; i < (*all_cells).size(); i++ )
{
if( (*all_cells)[i]->is_out_of_domain == false )
{
(*all_cells)[i]->advance_bundled_phenotype_functions( time_since_last_cycle );
(*all_cells)[i]->advance_bundled_phenotype_functions(time_since_last_phenotype);
}
}

// process divides / removes
for( int i=0; i < cells_ready_to_divide.size(); i++ )
{
cells_ready_to_divide[i]->divide();
}
// remove cells here (and not only below after interactions) in case phenotype and mechanics time steps don't line up
for( int i=0; i < cells_ready_to_die.size(); i++ )
{
cells_ready_to_die[i]->die();
{
cells_ready_to_die[i]->die();
}
num_divisions_in_current_step+= cells_ready_to_divide.size();
num_deaths_in_current_step+= cells_ready_to_die.size();

cells_ready_to_die.clear();
cells_ready_to_divide.clear();
last_cell_cycle_time= t;
}

double time_since_last_mechanics= t- last_mechanics_time;

// if( time_since_last_mechanics>= mechanics_dt || !initialzed)
if( time_since_last_mechanics > mechanics_dt_ - 0.5 * diffusion_dt_ || !initialzed )
time_since_last_phenotype = diffusion_dt_; // setting it for next cycle
}
else
{ time_since_last_phenotype += diffusion_dt_; }

if( time_for_mechanics )
{
if(!initialzed)
// new March 2022:
// run standard interactions (phagocytosis, attack, fusion) here
#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);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We need to discuss this change in call order further, as it will impact phagocytosis, attack, and fusion. I noticed it also relates to other PRs (#340 and #341). This change might introduce too many side effects.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, very happy to discuss more. By side effects, do you mean unintended behavior and/or errors? Or do you mean changes from how it ran before?

If we find this to be a long discussion, I would prefer to get a quick fix for this bug while we keep working on a longterm fix. Perhaps something like (the highly inefficient) #339?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I mean unintended behavior (like in PR #340 ), as the attach and detach dynamics were working fine in the previous release. I feel that by making this call order change, we might introduce these issues and possibly others we haven’t considered yet. Does that make sense?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, that makes sense. To be clear, #340 and #341 are built off the dev branch / previous release (1.14.0), not this fix. Their point of origin is how we implemented changes in 1.14.0. This fix does not change the existence of those behaviors.

}
for( int i=0; i < cells_ready_to_die.size(); i++ )
{
time_since_last_mechanics = mechanics_dt_;
cells_ready_to_die[i]->die();
}
cells_ready_to_die.clear();

// new February 2018
// if we need gradients, compute them
Expand Down Expand Up @@ -271,30 +277,6 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me
}
}

// new March 2022:
// run standard interactions (phagocytosis, attack, fusion) here
#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);
}
// 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 )
{
/*
std::cout << "\tClearing dummy cells from phagocytosis and fusion events ... " << std::endl;
std::cout << "\t\tClearing " << cells_ready_to_die.size() << " cells ... " << std::endl;
// there might be a lot of "dummy" cells ready for removal. Let's do it.
*/
for( int i=0; i < cells_ready_to_die.size(); i++ )
{ cells_ready_to_die[i]->die(); }
cells_ready_to_die.clear();
}


// update positions

#pragma omp parallel for
Expand All @@ -305,16 +287,19 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me
{ pC->update_position(time_since_last_mechanics); }
}

// When somebody reviews this code, let's add proper braces for clarity!!!

// Update cell indices in the container
for( int i=0; i < (*all_cells).size(); i++ )
{
if(!(*all_cells)[i]->is_out_of_domain && (*all_cells)[i]->is_movable)
{
(*all_cells)[i]->update_voxel_in_container();
last_mechanics_time=t;
}
}
time_since_last_mechanics = diffusion_dt_; // setting it for next cycle
}
else
{ time_since_last_mechanics += diffusion_dt_; }

initialzed=true;
return;
}

Expand Down
1 change: 0 additions & 1 deletion core/PhysiCell_cell_container.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,6 @@ class Cell_Container : public BioFVM::Agent_Container
std::vector<Cell*> cells_ready_to_divide; // the index of agents ready to divide
std::vector<Cell*> cells_ready_to_die;
int boundary_condition_for_pushed_out_agents; // what to do with pushed out cells
bool initialzed = false;

public:
BioFVM::Cartesian_Mesh underlying_mesh;
Expand Down
48 changes: 24 additions & 24 deletions core/PhysiCell_standard_models.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1317,36 +1317,36 @@ void standard_cell_cell_interactions( Cell* pCell, Phenotype& phenotype, double

// move effector attack here.

if( pCell->phenotype.cell_interactions.pAttackTarget != NULL )
{
Cell* pTarget = pCell->phenotype.cell_interactions.pAttackTarget;
if( pCell->phenotype.cell_interactions.pAttackTarget != NULL )
{
Cell* pTarget = pCell->phenotype.cell_interactions.pAttackTarget;

pCell->attack_cell(pTarget,dt);
attacked = true; // attacked at least one cell in this time step
pCell->attack_cell(pTarget,dt);
attacked = true; // attacked at least one cell in this time step

// attack_cell
// attack_cell

// probability of ending attack
// end attack if target is dead
probability = dt / (1e-15 + pCell->phenotype.cell_interactions.attack_duration);
// probability of ending attack
// end attack if target is dead
probability = dt / (1e-15 + pCell->phenotype.cell_interactions.attack_duration);


if( UniformRandom() < probability || pTarget->phenotype.death.dead )
{
/*
std::cout << "********* ********* ******** attack done **** " << PhysiCell_globals.current_time << " "
<< probability << " "
<< "attack time: " << pTarget->state.total_attack_time << " "
<< "damage: " << pTarget->phenotype.cell_integrity.damage << " "
<< "dead? " << (int) pTarget->phenotype.death.dead << " "
<< "damage delivered: " << pCell->phenotype.cell_interactions.total_damage_delivered << std::endl;
*/

detach_cells_as_spring(pCell,pTarget);

pCell->phenotype.cell_interactions.pAttackTarget = NULL;
}
if( UniformRandom() < probability || pTarget->phenotype.death.dead )
{
/*
std::cout << "********* ********* ******** attack done **** " << PhysiCell_globals.current_time << " "
<< probability << " "
<< "attack time: " << pTarget->state.total_attack_time << " "
<< "damage: " << pTarget->phenotype.cell_integrity.damage << " "
<< "dead? " << (int) pTarget->phenotype.death.dead << " "
<< "damage delivered: " << pCell->phenotype.cell_interactions.total_damage_delivered << std::endl;
*/

detach_cells_as_spring(pCell,pTarget);

pCell->phenotype.cell_interactions.pAttackTarget = NULL;
}
}


// move decision if to end attack here.
Expand Down
Loading