diff --git a/.gitignore b/.gitignore index 765ccdc8c..63546c554 100644 --- a/.gitignore +++ b/.gitignore @@ -21,3 +21,4 @@ user_projects/* !user_projects/empty.txt Studio.zip /studio +.vscode \ No newline at end of file diff --git a/Makefile b/Makefile index 5d256d1c3..acd1d9f23 100644 --- a/Makefile +++ b/Makefile @@ -74,7 +74,7 @@ list-projects: @echo "Sample projects: template biorobots-sample cancer-biorobots-sample cancer-immune-sample" @echo " celltypes3-sample heterogeneity-sample pred-prey-farmer virus-macrophage-sample" @echo " worm-sample interaction-sample mechano-sample rules-sample physimess-sample custom-division-sample" - @echo " asymmetric-division-sample immune-function-sample episode-sample" + @echo " asymmetric-division-sample extended-asym-div-sample immune-function-sample episode-sample" @echo "" @echo "Sample intracellular projects: template_BM ode-energy-sample physiboss-cell-lines-sample" @echo " cancer-metabolism-sample physiboss-tutorial physiboss-tutorial-invasion" @@ -178,6 +178,7 @@ mechano-sample: cp ./sample_projects/mechano/main.cpp ./main.cpp cp Makefile Makefile-backup cp ./sample_projects/mechano/Makefile . + cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml cp ./sample_projects/mechano/config/* ./config/ rules-sample: @@ -231,8 +232,18 @@ episode-sample: cp ./sample_projects/episode/main.cpp ./main.cpp cp Makefile Makefile-backup cp ./sample_projects/episode/Makefile . + cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml cp -r ./sample_projects/episode/config/* ./config +extended-asym-div-sample: + cp -r ./sample_projects/extended_asym_div/custom_modules/* ./custom_modules/ + touch main.cpp && cp main.cpp main-backup.cpp + cp ./sample_projects/extended_asym_div/main.cpp ./main.cpp + cp Makefile Makefile-backup + cp ./sample_projects/extended_asym_div/Makefile . + cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml + cp -r ./sample_projects/extended_asym_div/config/* ./config/ + # ---- intracellular projects ode-energy-sample: cp ./sample_projects_intracellular/ode/ode_energy/custom_modules/* ./custom_modules/ diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp index a14c70b99..6cad413fb 100644 --- a/core/PhysiCell_cell.cpp +++ b/core/PhysiCell_cell.cpp @@ -173,7 +173,6 @@ Cell_Definition::Cell_Definition() // are appropriately sized. Same on motiltiy. phenotype.cell_interactions.sync_to_cell_definitions(); phenotype.cell_transformations.sync_to_cell_definitions(); - phenotype.cycle.asymmetric_division.sync_to_cell_definitions(); phenotype.motility.sync_to_current_microenvironment(); phenotype.mechanics.sync_to_cell_definitions(); @@ -2073,8 +2072,6 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node // transformation pCD->phenotype.cell_transformations.transformation_rates.assign(number_of_cell_defs,0.0); - // asymmetric division - pCD->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities.assign(number_of_cell_defs,0.0); } else { @@ -2112,7 +2109,6 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node // this requires that prebuild_cell_definition_index_maps was already run pCD->phenotype.cell_interactions.sync_to_cell_definitions(); pCD->phenotype.cell_transformations.sync_to_cell_definitions(); - pCD->phenotype.cycle.asymmetric_division.sync_to_cell_definitions(); pCD->phenotype.mechanics.sync_to_cell_definitions(); // set the reference phenotype @@ -2268,8 +2264,13 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node node = node.child( "standard_asymmetric_division" ); if( node && node.attribute("enabled").as_bool() ) { + pugi::xml_node node_EAD = cd_node.child( "phenotype" ).child( "cycle" ).child( "extended_asymmetric_division" ); + if( node_EAD && node_EAD.attribute("enabled").as_bool() ) + { + std::cerr << "Error: Both standard and extended asymmetric division enabled for cell type " + pCD->name << std::endl; + exit(-1); + } Asymmetric_Division *pAD = &(pCD->phenotype.cycle.asymmetric_division); - // asymmetric division rates pugi::xml_node node_adp = node.child( "asymmetric_division_probability"); while (node_adp) @@ -2285,7 +2286,7 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node int target_index = search->second; double asymmetric_division_probability = xml_get_my_double_value(node_adp); - pAD->asymmetric_division_probabilities[target_index] = asymmetric_division_probability; + pAD->set_asymmetric_division_probability(pCD->type, target_index, asymmetric_division_probability); } else { @@ -2296,14 +2297,50 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node } node_adp = node_adp.next_sibling("asymmetric_division_probability"); } - std::cout << "Asymmetric division probabilities for " << pCD->name << ": "; - for (int i = 0; i < pAD->asymmetric_division_probabilities.size(); i++) + pCD->functions.cell_division_function = asymmetric_division_function; + } + + node = cd_node.child( "phenotype" ); + node = node.child( "cycle" ); + node = node.child( "extended_asymmetric_division" ); + if( node && node.attribute("enabled").as_bool() ) + { + Asymmetric_Division *pAD = &(pCD->phenotype.cycle.asymmetric_division); + + // asymmetric division rates + pugi::xml_node node_eadp = node.child( "extended_asymmetric_division_probability" ); + + while (node_eadp) { - std::cout << pAD->asymmetric_division_probabilities[i] << " "; + // get the name of the target cell type + std::string first_target_name = node_eadp.attribute("name1").value(); + std::string second_target_name = node_eadp.attribute("name2").value(); + // now find its index + auto first_search = cell_definition_indices_by_name.find(first_target_name); + auto second_search = cell_definition_indices_by_name.find(second_target_name); + + // safety first! + if( first_search == cell_definition_indices_by_name.end() || second_search == cell_definition_indices_by_name.end() ) + { + std::cout << "Error: When processing the " << pCD->name << " cell definition: " << std::endl + << "\tCould not find cell type " << first_target_name << " or " << second_target_name << " for asymmetric division." << std::endl + << "\tRemove this cell type from the extended asymmetric division probabilities!" << std::endl << std::endl; + exit(-1); + } + + // if the target is found, set the appropriate rate + int first_target_index = first_search->second; + int second_target_index = second_search->second; + + double extended_asymmetric_division_probability = xml_get_my_double_value(node_eadp); + pAD->set_asymmetric_division_probability(first_target_index, second_target_index, extended_asymmetric_division_probability); + + node_eadp = node_eadp.next_sibling("extended_asymmetric_division_probability"); } std::cout << std::endl; - pCD->functions.cell_division_function = standard_asymmetric_division_function; + pCD->functions.cell_division_function = asymmetric_division_function; } + } // here's what it ***should*** do: diff --git a/core/PhysiCell_phenotype.cpp b/core/PhysiCell_phenotype.cpp index 4d03249be..8f6a53743 100644 --- a/core/PhysiCell_phenotype.cpp +++ b/core/PhysiCell_phenotype.cpp @@ -1358,36 +1358,93 @@ double& Cell_Transformations::transformation_rate( std::string type_name ) return transformation_rates[n]; } -Asymmetric_Division::Asymmetric_Division() +void Asymmetric_Division::set_asymmetric_division_probability(std::pair types, double probability) { - asymmetric_division_probabilities = {0.0}; + asymmetric_division_probabilities[types] = probability; } -void Asymmetric_Division::sync_to_cell_definitions() +void Asymmetric_Division::set_asymmetric_division_probability(int upper_triangular_index, double probability) +{ + std::pair types = extended_asym_index_to_upper_triangle(upper_triangular_index); + set_asymmetric_division_probability(types, probability); +} + +void Asymmetric_Division::set_asymmetric_division_probability(int type_1, int type_2, double probability) +{ + std::pair types = std::make_pair(type_1, type_2); + set_asymmetric_division_probability(types, probability); +} + +void Asymmetric_Division::set_asymmetric_division_probability(std::string type_name_1, std::string type_name_2, double probability) { extern std::unordered_map cell_definition_indices_by_name; - int number_of_cell_defs = cell_definition_indices_by_name.size(); - - if( asymmetric_division_probabilities.size() != number_of_cell_defs ) - { asymmetric_division_probabilities.resize( number_of_cell_defs, 0.0); } - - return; + int n = cell_definition_indices_by_name[type_name_1]; + int m = cell_definition_indices_by_name[type_name_2]; + set_asymmetric_division_probability(n, m, probability); +} + +double Asymmetric_Division::asymmetric_division_probability( std::pair types ) +{ + if (asymmetric_division_probabilities.find(types) == asymmetric_division_probabilities.end()) + { return 0.0; } // key not found, assume this means 0 probability of this type of division for the cell + return asymmetric_division_probabilities[types]; +} + +double Asymmetric_Division::asymmetric_division_probability( int type_1, int type_2 ) +{ + return asymmetric_division_probability(std::make_pair(type_1, type_2)); +} + +double Asymmetric_Division::asymmetric_division_probability( int upper_triangular_index ) +{ + std::pair types = extended_asym_index_to_upper_triangle(upper_triangular_index); + return asymmetric_division_probability(types); +} + +double Asymmetric_Division::asymmetric_division_probability( std::string type_name_1, std::string type_name_2 ) +{ + extern std::unordered_map cell_definition_indices_by_name; + int n = cell_definition_indices_by_name[type_name_1]; + int m = cell_definition_indices_by_name[type_name_2]; + return asymmetric_division_probability(n, m); +} + +std::pair Asymmetric_Division::select_daughter_types(int type_1, int type_2) +{ + double r = UniformRandom(); + for( auto it = asymmetric_division_probabilities.begin(); it != asymmetric_division_probabilities.end(); ++it ) + { + if( r <= it->second ) + { return it->first; } // return the pair of types + } + return std::make_pair(type_1, type_2); // if we reach here, return the original types } double Asymmetric_Division::probabilities_total( void ) { double total = 0.0; - for( int i=0; i < asymmetric_division_probabilities.size(); i++ ) - { total += asymmetric_division_probabilities[i]; } + for (const auto& pair : asymmetric_division_probabilities) + { total += pair.second; } return total; } -// ease of access -double& Asymmetric_Division::asymmetric_division_probability( std::string type_name ) +std::pair extended_asym_index_to_upper_triangle(int index) +{ + static std::vector< std::pair > pairs_vector = initialize_pairs_vector(); + return pairs_vector[index]; +} + +std::vector< std::pair > initialize_pairs_vector() { extern std::unordered_map cell_definition_indices_by_name; - int n = cell_definition_indices_by_name[type_name]; - return asymmetric_division_probabilities[n]; + std::vector< std::pair > output; + int n = cell_definition_indices_by_name.size(); + for( int i = 0; i < n; i++ ) + { + for( int j = i; j < n; j++ ) + { output.push_back( std::make_pair(i,j) ); } + } + return output; } // beta functionality in 1.10.3 diff --git a/core/PhysiCell_phenotype.h b/core/PhysiCell_phenotype.h index 701d816b5..b82b0bc72 100644 --- a/core/PhysiCell_phenotype.h +++ b/core/PhysiCell_phenotype.h @@ -213,23 +213,50 @@ class Cycle_Model std::ostream& display( std::ostream& os ); // done }; +// A hash function for pairs of ints as keys for extended_asymmetric_division_probabilities. +// Should follow upper-triangular variant of Cantor set to prevent collisions. +struct pair_hash { + + std::size_t operator () (const std::pair& pair) const + { + int lower = std::min(pair.first, pair.second); + int upper = std::max(pair.first, pair.second); + int UT = upper * (upper + 1) / 2 + lower; + auto hash = std::hash{}(UT); + return hash; + } +}; + +struct equality_function { + bool operator()(const std::pair& lhs, const std::pair& rhs) const + { + return ((lhs.first == rhs.first && lhs.second == rhs.second) || (lhs.first == rhs.second && lhs.second == rhs.first)); + } +}; + class Asymmetric_Division { -private: public: - // rates of asymmetric division into different cell types - std::vector asymmetric_division_probabilities; + std::unordered_map, double, pair_hash, equality_function> asymmetric_division_probabilities; - // initialization - Asymmetric_Division(); // done - void sync_to_cell_definitions(); // done + void set_asymmetric_division_probability(std::pair types, double probability); + void set_asymmetric_division_probability(int upper_triangular_index, double probability); + void set_asymmetric_division_probability(int type_1, int type_2, double probability); + void set_asymmetric_division_probability(std::string type_name_1, std::string type_name_2, double probability); + + double asymmetric_division_probability(std::pair types); + double asymmetric_division_probability(int upper_triangular_index); + double asymmetric_division_probability(int type_1, int type_2); + double asymmetric_division_probability(std::string type_name_1, std::string type_name_2); double probabilities_total(); - // ease of access - double& asymmetric_division_probability( std::string type_name ); // done + std::pair select_daughter_types(int type_1, int type_2); }; +std::pair extended_asym_index_to_upper_triangle(int index); +std::vector> initialize_pairs_vector(void); + class Cycle { private: diff --git a/core/PhysiCell_signal_behavior.cpp b/core/PhysiCell_signal_behavior.cpp index 8780e8d7c..57cfa1315 100644 --- a/core/PhysiCell_signal_behavior.cpp +++ b/core/PhysiCell_signal_behavior.cpp @@ -538,7 +538,7 @@ void setup_signal_behavior_dictionaries( void ) } - // asymmetic division + // asymmetric division (this remains necessary even with extended asym div because this item implicitly uses the current cell type name, meaning the index depends on the cell using this rule) for( int i=0; i < n ; i++ ) { map_index++; @@ -548,6 +548,23 @@ void setup_signal_behavior_dictionaries( void ) int_to_behavior[map_index] = temp; } + // extended asym div + for( int i=0; i < n ; i++ ) + { + for( int j=i; j < n ; j++ ) + { + map_index++; + Cell_Definition* pCD1 = cell_definitions_by_type[i]; + Cell_Definition* pCD2 = cell_definitions_by_type[j]; + std::string signal_str = "extended asymmetric division to " + pCD1->name + " and " + pCD2->name; + behavior_to_int[signal_str] = map_index; + signal_str = "extended asymmetric division to " + pCD2->name + " and " + pCD1->name; + behavior_to_int[signal_str] = map_index; + + int_to_behavior[map_index] = signal_str; + } + } + // custom behaviors for( int nc=0 ; nc < cell_defaults.custom_data.variables.size() ; nc++ ) { @@ -1417,9 +1434,18 @@ void set_behaviors( Cell* pCell , std::vector parameters ) // asymmetric division static int first_asymmetric_division_index = find_behavior_index( "asymmetric division to " + cell_definitions_by_type[0]->name ); - std::copy( parameters.begin()+first_asymmetric_division_index , - parameters.begin()+first_asymmetric_division_index+n , - pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities.begin() ); + for ( int i = 0; i < n ; i++ ) + { + pCell->phenotype.cycle.asymmetric_division.set_asymmetric_division_probability(pCell->type, i, parameters[i + first_asymmetric_division_index]); + } + + // extended asym div + std::cout << "Warning: set_behaviors is overwriting values for standard asymmetric division with extended asymmetric division probabilities!" << std::endl; + static int first_extended_asymmetric_division_index = find_behavior_index( "extended asymmetric division to " + cell_definitions_by_type[0]->name + " and " + cell_definitions_by_type[0]->name ); + for( int ind = 0; ind < (n + 1) * n / 2; ind++ ) + { + pCell->phenotype.cycle.asymmetric_division.set_asymmetric_division_probability(ind, parameters[ind + first_extended_asymmetric_division_index]); + } // custom behaviors static int first_custom_ind = find_behavior_index( "custom 0"); @@ -1626,7 +1652,12 @@ void set_single_behavior( Cell* pCell, int index , double parameter ) // asymmetric division static int first_asymmetric_division_index = find_behavior_index( "asymmetric division to " + cell_definitions_by_type[0]->name ); if( index >= first_asymmetric_division_index && index < first_asymmetric_division_index + n ) - { pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities[index-first_asymmetric_division_index] = parameter; return; } + { pCell->phenotype.cycle.asymmetric_division.set_asymmetric_division_probability(pCell->type, index-first_asymmetric_division_index, parameter); return; } + + // extended asym div + static int first_extended_asymmetric_division_index = find_behavior_index( "extended asymmetric division to " + cell_definitions_by_type[0]->name + " and " + cell_definitions_by_type[0]->name ); + if( index >= first_extended_asymmetric_division_index && index < first_extended_asymmetric_division_index + n * (n + 1) / 2 ) + { pCell->phenotype.cycle.asymmetric_division.set_asymmetric_division_probability(index - first_extended_asymmetric_division_index, parameter); return; } // custom behavior static int first_custom_ind = find_behavior_index( "custom 0"); @@ -1837,9 +1868,17 @@ std::vector get_behaviors( Cell* pCell ) // asymmetric division static int first_asymmetric_division_index = find_behavior_index( "asymmetric division to " + cell_definitions_by_type[0]->name ); - std::copy( pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities.begin(), - pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities.end(), - parameters.begin()+first_asymmetric_division_index ); + for ( int i = 0; i < n ; i++ ) + { + parameters[i + first_asymmetric_division_index] = pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probability(pCell->type, i); + } + + // extended asym div + static int first_extended_asymmetric_division_index = find_behavior_index( "extended asymmetric division to " + cell_definitions_by_type[0]->name + " and " + cell_definitions_by_type[0]->name ); + for( int ind = 0; ind < (n + 1) * n / 2; ind++ ) + { + parameters[ind + first_extended_asymmetric_division_index] = pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probability(ind); + } // custom behavior static int first_custom_ind = find_behavior_index( "custom 0"); @@ -2052,7 +2091,12 @@ double get_single_behavior( Cell* pCell , int index ) // asymmetric division static int first_asymmetric_division_index = find_behavior_index( "asymmetric division to " + cell_definitions_by_type[0]->name ); if( index >= first_asymmetric_division_index && index < first_asymmetric_division_index+n ) - { return pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities[index-first_asymmetric_division_index]; } + { return pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probability(pCell->type, index-first_asymmetric_division_index); } + + // extended asym div + static int first_extended_asymmetric_division_index = find_behavior_index( "extended asymmetric division to " + cell_definitions_by_type[0]->name + " and " + cell_definitions_by_type[0]->name ); + if( index >= first_extended_asymmetric_division_index && index < first_extended_asymmetric_division_index + n * (n + 1) / 2 ) + { return pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probability(index - first_extended_asymmetric_division_index); } // custom behavior static int first_custom_ind = find_behavior_index( "custom 0"); @@ -2293,9 +2337,17 @@ std::vector get_base_behaviors( Cell* pCell ) // asymmetric division static int first_asymmetric_division_index = find_behavior_index( "asymmetric division to " + cell_definitions_by_type[0]->name ); - std::copy( pCD->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities.begin(), - pCD->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities.end(), - parameters.begin()+first_asymmetric_division_index ); + for ( int i = 0; i < n ; i++ ) + { + parameters[i + first_asymmetric_division_index] = pCD->phenotype.cycle.asymmetric_division.asymmetric_division_probability(pCD->type, i); + } + + // extended asym div + static int first_extended_asymmetric_division_index = find_behavior_index( "extended asymmetric division to " + cell_definitions_by_type[0]->name + " and " + cell_definitions_by_type[0]->name ); + for( int ind = 0; ind < (n + 1) * n / 2; ind++ ) + { + parameters[ind + first_extended_asymmetric_division_index] = pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probability(ind); + } // custom behavior static int first_custom_ind = find_behavior_index( "custom 0"); @@ -2511,7 +2563,12 @@ double get_single_base_behavior( Cell* pCell , int index ) // asymmetric division static int first_asymmetric_division_index = find_behavior_index( "asymmetric division to " + cell_definitions_by_type[0]->name ); if( index >= first_asymmetric_division_index && index < first_asymmetric_division_index + n ) - { return pCD->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities[index-first_asymmetric_division_index]; } + { return pCD->phenotype.cycle.asymmetric_division.asymmetric_division_probability(pCell->type, index-first_asymmetric_division_index); } + + // extended asym div + static int first_extended_asymmetric_division_index = find_behavior_index( "extended asymmetric division to " + cell_definitions_by_type[0]->name + " and " + cell_definitions_by_type[0]->name ); + if( index >= first_extended_asymmetric_division_index && index < first_extended_asymmetric_division_index + n * (n + 1) / 2 ) + { return pCD->phenotype.cycle.asymmetric_division.asymmetric_division_probability(index - first_extended_asymmetric_division_index); } // custom behavior static int first_custom_ind = find_behavior_index( "custom 0"); @@ -2733,7 +2790,12 @@ double get_single_base_behavior( Cell_Definition* pCD , int index ) // asymmetric division static int first_asymmetric_division_index = find_behavior_index( "asymmetric division to " + cell_definitions_by_type[0]->name ); if( index >= first_asymmetric_division_index && index < first_asymmetric_division_index + n ) - { return pCD->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities[index-first_asymmetric_division_index]; } + { return pCD->phenotype.cycle.asymmetric_division.asymmetric_division_probability(pCD->type, index-first_asymmetric_division_index); } + + // extended asym div + static int first_extended_asymmetric_division_index = find_behavior_index( "extended asymmetric division to " + cell_definitions_by_type[0]->name + " and " + cell_definitions_by_type[0]->name ); + if( index >= first_extended_asymmetric_division_index && index < first_extended_asymmetric_division_index + n * (n + 1) / 2 ) + { return pCD->phenotype.cycle.asymmetric_division.asymmetric_division_probability(index - first_extended_asymmetric_division_index); } // custom behavior static int first_custom_ind = find_behavior_index( "custom 0"); diff --git a/core/PhysiCell_signal_behavior.h b/core/PhysiCell_signal_behavior.h index 29a7cc128..f6db9641d 100644 --- a/core/PhysiCell_signal_behavior.h +++ b/core/PhysiCell_signal_behavior.h @@ -156,7 +156,7 @@ void set_single_behavior( Cell* pCell, std::string name , double parameter ); // std::vector get_behaviors( Cell* pCell ); // done // get selected current behavior -std::vector get_behaviors( Cell* pCell , std::vector indices ); // doen +std::vector get_behaviors( Cell* pCell , std::vector indices ); // done std::vector get_behaviors( Cell* pCell , std::vector names ); // done // get single current behavior diff --git a/core/PhysiCell_standard_models.cpp b/core/PhysiCell_standard_models.cpp index 48a70ccb8..d8ea0d933 100644 --- a/core/PhysiCell_standard_models.cpp +++ b/core/PhysiCell_standard_models.cpp @@ -779,7 +779,6 @@ void initialize_default_cell_definition( void ) // these for the cell_defaults cell_defaults.phenotype.cell_interactions.sync_to_cell_definitions(); cell_defaults.phenotype.cell_transformations.sync_to_cell_definitions(); - cell_defaults.phenotype.cycle.asymmetric_division.sync_to_cell_definitions(); cell_defaults.phenotype.motility.sync_to_current_microenvironment(); cell_defaults.phenotype.mechanics.sync_to_cell_definitions(); @@ -1373,52 +1372,41 @@ void standard_cell_transformations( Cell* pCell, Phenotype& phenotype, double dt } -void standard_asymmetric_division_function( Cell* pCell_parent, Cell* pCell_daughter ) +void asymmetric_division_function( Cell* pCell_parent, Cell* pCell_daughter ) { - Cell_Definition* pCD_parent = cell_definitions_by_name[pCell_parent->type_name]; + std::string parent_name = pCell_parent->type_name; + int parent_type = pCell_parent->type; + Cell_Definition* pCD_parent = cell_definitions_by_name[parent_name]; double total = pCell_parent->phenotype.cycle.asymmetric_division.probabilities_total(); if (total > 1.0) { - double sym_div_prob = pCell_parent->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities[pCell_parent->type] + 1.0 - total; + double sym_div_prob = pCell_parent->phenotype.cycle.asymmetric_division.asymmetric_division_probability(parent_type, parent_type) + 1.0 - total; if (sym_div_prob < 0.0) { - throw std::runtime_error("Error: Asymmetric division probabilities for " + pCD_parent->name + " sum to greater than 1.0 and cannot be normalized."); - } - pCell_parent->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities[pCell_parent->type] = sym_div_prob; - pCell_daughter->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities[pCell_daughter->type] = sym_div_prob; - } - double r = UniformRandom(); - for( int i=0; i < pCD_parent->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities.size(); i++ ) - { - if( r <= pCell_parent->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities[i] ) - { - if (i != pCell_daughter->type) // only convert if the daughter is not already the correct type - { pCell_daughter->convert_to_cell_definition( *cell_definitions_by_index[i] ); } - return; + std::cerr << "Error: Asymmetric division probabilities for " + pCD_parent->name + " sum to greater than 1.0 and cannot be normalized." << std::endl; + std::cerr << "Adjusted sym_div_prob = " << sym_div_prob << std::endl; + std::cerr << "List of all asym div probabilities:" << std::endl; + for (int i = 0; i < cell_definitions_by_index.size(); i++) + { + for (int j = i; j < cell_definitions_by_index.size(); j++) + { + std::cerr << " - " << cell_definitions_by_index[i]->name << " and " << cell_definitions_by_index[j]->name << ": " << pCell_parent->phenotype.cycle.asymmetric_division.asymmetric_division_probability(i, j) << std::endl; + } + } + exit(-1); } - r -= pCell_parent->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities[i]; + pCell_parent->phenotype.cycle.asymmetric_division.set_asymmetric_division_probability(parent_type, parent_type, sym_div_prob); + pCell_daughter->phenotype.cycle.asymmetric_division.set_asymmetric_division_probability(pCell_daughter->type, pCell_daughter->type, sym_div_prob); } - // if we're here, then do not do asym div + std::pair daughter_types = pCell_parent->phenotype.cycle.asymmetric_division.select_daughter_types(pCell_parent->type, pCell_daughter->type); + + if (daughter_types.first != pCell_parent->type) // only convert if the parent is not already the correct type + { pCell_parent->convert_to_cell_definition( *cell_definitions_by_index[daughter_types.first] ); } + if (daughter_types.second != pCell_daughter->type) // only convert if the daughter is not already the correct type + { pCell_daughter->convert_to_cell_definition( *cell_definitions_by_index[daughter_types.second] ); } return; } -// alternative way to select the index from weights that could be faster (is faster as # cell types --> infinity) -// int select_by_probabilities( const std::vector& probabilities ) -// { -// double r = UniformRandom(); - -// std::vector cumulative_weights(probabilities.size()); -// std::partial_sum(probabilities.begin(), probabilities.end(), cumulative_weights.begin()); - -// // Use binary search to find the index -// auto it = std::upper_bound(cumulative_weights.begin(), cumulative_weights.end(), r); -// int index = std::distance(cumulative_weights.begin(), it); -// if (index >= probabilities.size()) -// { return -1; } - -// return index; -// } - void dynamic_attachments( Cell* pCell , Phenotype& phenotype, double dt ) { // check for detachments diff --git a/core/PhysiCell_standard_models.h b/core/PhysiCell_standard_models.h index c67e2c3f7..fcba08469 100644 --- a/core/PhysiCell_standard_models.h +++ b/core/PhysiCell_standard_models.h @@ -143,7 +143,7 @@ void evaluate_interactions( Cell* pCell, Phenotype& phenotype, double dt ); // automated cell phagocytosis, attack, and fusion void standard_cell_cell_interactions( Cell* pCell, Phenotype& phenotype, double dt ); void standard_cell_transformations( Cell* pCell, Phenotype& phenotype, double dt ); -void standard_asymmetric_division_function( Cell* pCell_parent, Cell* pCell_daughter ); +void asymmetric_division_function( Cell* pCell_parent, Cell* pCell_daughter ); void advanced_chemotaxis_function_normalized( Cell* pCell, Phenotype& phenotype , double dt ); void advanced_chemotaxis_function( Cell* pCell, Phenotype& phenotype , double dt ); diff --git a/modules/PhysiCell_MultiCellDS.cpp b/modules/PhysiCell_MultiCellDS.cpp index 92b06743e..cdb984894 100644 --- a/modules/PhysiCell_MultiCellDS.cpp +++ b/modules/PhysiCell_MultiCellDS.cpp @@ -570,7 +570,7 @@ void add_PhysiCell_cells_to_open_xml_pugi_v2( pugi::xml_document& xml_dom, std:: // asymmetric division // std::vector asymmetric_division_probabilities; // n add_variable_to_labels( data_names, data_units, data_start_indices, data_sizes, - "asymmetric_division_probabilities" , "none" , n ); + "asymmetric_division_probabilities" , "none" , n * (n+1) / 2 ); // cell integrity @@ -983,7 +983,14 @@ void add_PhysiCell_cells_to_open_xml_pugi_v2( pugi::xml_document& xml_dom, std:: // asymmetric division // name = "asymmetric_division_rate"; - std::fwrite( pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probabilities.data() , sizeof(double) , n, fp ); + for ( int i1 = 0; i1 < n; i1++ ) + { + for ( int i2 = i1; i2 < n; i2++ ) + { + double prob = pCell->phenotype.cycle.asymmetric_division.asymmetric_division_probability(i1, i2); + std::fwrite(&prob, sizeof(double), 1, fp); + } + } // cell integrity // name = "damage"; diff --git a/modules/PhysiCell_settings.cpp b/modules/PhysiCell_settings.cpp index 538f6a8ea..09fd77158 100644 --- a/modules/PhysiCell_settings.cpp +++ b/modules/PhysiCell_settings.cpp @@ -357,7 +357,7 @@ bool create_directories(const std::string &path) bool create_directory(const std::string &path) { -#if defined(__MINGW32__) || defined(__MINGW64__) +#if defined(_WIN32) bool success = mkdir(path.c_str()) == 0; #else bool success = mkdir(path.c_str(), 0755) == 0; diff --git a/sample_projects/Makefile-default b/sample_projects/Makefile-default index 5d256d1c3..acd1d9f23 100644 --- a/sample_projects/Makefile-default +++ b/sample_projects/Makefile-default @@ -74,7 +74,7 @@ list-projects: @echo "Sample projects: template biorobots-sample cancer-biorobots-sample cancer-immune-sample" @echo " celltypes3-sample heterogeneity-sample pred-prey-farmer virus-macrophage-sample" @echo " worm-sample interaction-sample mechano-sample rules-sample physimess-sample custom-division-sample" - @echo " asymmetric-division-sample immune-function-sample episode-sample" + @echo " asymmetric-division-sample extended-asym-div-sample immune-function-sample episode-sample" @echo "" @echo "Sample intracellular projects: template_BM ode-energy-sample physiboss-cell-lines-sample" @echo " cancer-metabolism-sample physiboss-tutorial physiboss-tutorial-invasion" @@ -178,6 +178,7 @@ mechano-sample: cp ./sample_projects/mechano/main.cpp ./main.cpp cp Makefile Makefile-backup cp ./sample_projects/mechano/Makefile . + cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml cp ./sample_projects/mechano/config/* ./config/ rules-sample: @@ -231,8 +232,18 @@ episode-sample: cp ./sample_projects/episode/main.cpp ./main.cpp cp Makefile Makefile-backup cp ./sample_projects/episode/Makefile . + cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml cp -r ./sample_projects/episode/config/* ./config +extended-asym-div-sample: + cp -r ./sample_projects/extended_asym_div/custom_modules/* ./custom_modules/ + touch main.cpp && cp main.cpp main-backup.cpp + cp ./sample_projects/extended_asym_div/main.cpp ./main.cpp + cp Makefile Makefile-backup + cp ./sample_projects/extended_asym_div/Makefile . + cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml + cp -r ./sample_projects/extended_asym_div/config/* ./config/ + # ---- intracellular projects ode-energy-sample: cp ./sample_projects_intracellular/ode/ode_energy/custom_modules/* ./custom_modules/ diff --git a/sample_projects/extended_asym_div/Makefile b/sample_projects/extended_asym_div/Makefile new file mode 100644 index 000000000..865e7c5b8 --- /dev/null +++ b/sample_projects/extended_asym_div/Makefile @@ -0,0 +1,327 @@ +VERSION := $(shell grep . VERSION.txt | cut -f1 -d:) +PROGRAM_NAME := project + +CC := g++ +# CC := g++-mp-7 # typical macports compiler name +# CC := g++-7 # typical homebrew compiler name + +# Check for environment definitions of compiler +# e.g., on CC = g++-7 on OSX +ifdef PHYSICELL_CPP + CC := $(PHYSICELL_CPP) +endif + +ifndef STATIC_OPENMP + STATIC_OPENMP = -fopenmp +endif + +ARCH := native # best auto-tuning +# ARCH := core2 # a reasonably safe default for most CPUs since 2007 +# ARCH := corei7 +# ARCH := corei7-avx # earlier i7 +# ARCH := core-avx-i # i7 ivy bridge or newer +# ARCH := core-avx2 # i7 with Haswell or newer +# ARCH := nehalem +# ARCH := westmere +# ARCH := sandybridge # circa 2011 +# ARCH := ivybridge # circa 2012 +# ARCH := haswell # circa 2013 +# ARCH := broadwell # circa 2014 +# ARCH := skylake # circa 2015 +# ARCH := bonnell +# ARCH := silvermont +# ARCH := skylake-avx512 +# ARCH := nocona #64-bit pentium 4 or later + +# CFLAGS := -march=$(ARCH) -Ofast -s -fomit-frame-pointer -mfpmath=both -fopenmp -m64 -std=c++11 +CFLAGS := -march=$(ARCH) -O3 -fomit-frame-pointer -mfpmath=both -fopenmp -m64 -std=c++11 + +ifeq ($(OS),Windows_NT) +else + UNAME_S := $(shell uname -s) + ifeq ($(UNAME_S),Darwin) + UNAME_P := $(shell uname -p) + var := $(shell which $(CC) | xargs file) + ifeq ($(lastword $(var)),arm64) + CFLAGS := -march=$(ARCH) -O3 -fomit-frame-pointer -fopenmp -m64 -std=c++11 + endif + endif +endif + +CFLAGS_LINK := $(shell echo $(CFLAGS) | sed -e "s/-fopenmp//g") +COMPILE_COMMAND := $(CC) $(CFLAGS) $(EXTRA_FLAGS) +LINK_COMMAND := $(CC) $(CFLAGS_LINK) $(EXTRA_FLAGS) + +BioFVM_OBJECTS := BioFVM_vector.o BioFVM_mesh.o BioFVM_microenvironment.o BioFVM_solvers.o BioFVM_matlab.o \ +BioFVM_utilities.o BioFVM_basic_agent.o BioFVM_MultiCellDS.o BioFVM_agent_container.o + +PhysiCell_core_OBJECTS := PhysiCell_phenotype.o PhysiCell_cell_container.o PhysiCell_standard_models.o \ +PhysiCell_cell.o PhysiCell_custom.o PhysiCell_utilities.o PhysiCell_constants.o PhysiCell_basic_signaling.o \ +PhysiCell_signal_behavior.o PhysiCell_rules.o + +PhysiCell_module_OBJECTS := PhysiCell_SVG.o PhysiCell_pathology.o PhysiCell_MultiCellDS.o PhysiCell_various_outputs.o \ +PhysiCell_pugixml.o PhysiCell_settings.o PhysiCell_geometry.o + +# put your custom objects here (they should be in the custom_modules directory) + +PhysiCell_custom_module_OBJECTS := custom.o + +pugixml_OBJECTS := pugixml.o + +PhysiCell_OBJECTS := $(BioFVM_OBJECTS) $(pugixml_OBJECTS) $(PhysiCell_core_OBJECTS) $(PhysiCell_module_OBJECTS) +ALL_OBJECTS := $(PhysiCell_OBJECTS) $(PhysiCell_custom_module_OBJECTS) + +# compile the project + +all: main.cpp $(ALL_OBJECTS) + $(COMPILE_COMMAND) -o $(PROGRAM_NAME) $(ALL_OBJECTS) main.cpp + make name + +static: main.cpp $(ALL_OBJECTS) $(MaBoSS) + $(LINK_COMMAND) $(INC) -o $(PROGRAM_NAME) $(ALL_OBJECTS) main.cpp $(LIB) -static-libgcc -static-libstdc++ $(STATIC_OPENMP) + +name: + @echo "" + @echo "Executable name is" $(PROGRAM_NAME) + @echo "" + +# PhysiCell core components + +PhysiCell_phenotype.o: ./core/PhysiCell_phenotype.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_phenotype.cpp + +PhysiCell_digital_cell_line.o: ./core/PhysiCell_digital_cell_line.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_digital_cell_line.cpp + +PhysiCell_cell.o: ./core/PhysiCell_cell.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_cell.cpp + +PhysiCell_cell_container.o: ./core/PhysiCell_cell_container.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_cell_container.cpp + +PhysiCell_standard_models.o: ./core/PhysiCell_standard_models.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_standard_models.cpp + +PhysiCell_utilities.o: ./core/PhysiCell_utilities.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_utilities.cpp + +PhysiCell_custom.o: ./core/PhysiCell_custom.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_custom.cpp + +PhysiCell_constants.o: ./core/PhysiCell_constants.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_constants.cpp + +PhysiCell_signal_behavior.o: ./core/PhysiCell_signal_behavior.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_signal_behavior.cpp + +PhysiCell_rules.o: ./core/PhysiCell_rules.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_rules.cpp + +# BioFVM core components (needed by PhysiCell) + +BioFVM_vector.o: ./BioFVM/BioFVM_vector.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_vector.cpp + +BioFVM_agent_container.o: ./BioFVM/BioFVM_agent_container.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_agent_container.cpp + +BioFVM_mesh.o: ./BioFVM/BioFVM_mesh.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_mesh.cpp + +BioFVM_microenvironment.o: ./BioFVM/BioFVM_microenvironment.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_microenvironment.cpp + +BioFVM_solvers.o: ./BioFVM/BioFVM_solvers.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_solvers.cpp + +BioFVM_utilities.o: ./BioFVM/BioFVM_utilities.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_utilities.cpp + +BioFVM_basic_agent.o: ./BioFVM/BioFVM_basic_agent.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_basic_agent.cpp + +BioFVM_matlab.o: ./BioFVM/BioFVM_matlab.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_matlab.cpp + +BioFVM_MultiCellDS.o: ./BioFVM/BioFVM_MultiCellDS.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_MultiCellDS.cpp + +pugixml.o: ./BioFVM/pugixml.cpp + $(COMPILE_COMMAND) -c ./BioFVM/pugixml.cpp + +# standard PhysiCell modules + +PhysiCell_SVG.o: ./modules/PhysiCell_SVG.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_SVG.cpp + +PhysiCell_pathology.o: ./modules/PhysiCell_pathology.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_pathology.cpp + +PhysiCell_MultiCellDS.o: ./modules/PhysiCell_MultiCellDS.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_MultiCellDS.cpp + +PhysiCell_various_outputs.o: ./modules/PhysiCell_various_outputs.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_various_outputs.cpp + +PhysiCell_pugixml.o: ./modules/PhysiCell_pugixml.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_pugixml.cpp + +PhysiCell_settings.o: ./modules/PhysiCell_settings.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_settings.cpp + +PhysiCell_basic_signaling.o: ./core/PhysiCell_basic_signaling.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_basic_signaling.cpp + +PhysiCell_geometry.o: ./modules/PhysiCell_geometry.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_geometry.cpp + +# user-defined PhysiCell modules + +custom.o: ./custom_modules/custom.cpp + $(COMPILE_COMMAND) -c ./custom_modules/custom.cpp + +# cleanup + +reset: + rm -f *.cpp + cp ./sample_projects/Makefile-default Makefile + rm -f ./custom_modules/* + touch ./custom_modules/empty.txt + touch ALL_CITATIONS.txt + touch ./core/PhysiCell_cell.cpp + rm ALL_CITATIONS.txt + cp ./config/PhysiCell_settings-backup.xml ./config/PhysiCell_settings.xml + touch ./config/empty.csv + rm -f ./config/*.csv + +clean: + rm -f *.o + rm -f $(PROGRAM_NAME)* + +data-cleanup: + rm -rf ./output + mkdir ./output + touch ./output/empty.txt + +# archival + +checkpoint: + zip -r $$(date +%b_%d_%Y_%H%M).zip Makefile *.cpp *.h config/*.xml custom_modules/* + +zip: + zip -r latest.zip Makefile* *.cpp *.h BioFVM/* config/* core/* custom_modules/* matlab/* modules/* sample_projects/* + cp latest.zip $$(date +%b_%d_%Y_%H%M).zip + cp latest.zip VERSION_$(VERSION).zip + mv *.zip archives/ + +tar: + tar --ignore-failed-read -czf latest.tar Makefile* *.cpp *.h BioFVM/* config/* core/* custom_modules/* matlab/* modules/* sample_projects/* + cp latest.tar $$(date +%b_%d_%Y_%H%M).tar + cp latest.tar VERSION_$(VERSION).tar + mv *.tar archives/ + +unzip: + cp ./archives/latest.zip . + unzip latest.zip + +untar: + cp ./archives/latest.tar . + tar -xzf latest.tar + +# easier animation + +FRAMERATE := 24 +OUTPUT := output + +jpeg: + @magick identify -format "%h" $(OUTPUT)/initial.svg > __H.txt + @magick identify -format "%w" $(OUTPUT)/initial.svg > __W.txt + @expr 2 \* \( $$(grep . __H.txt) / 2 \) > __H1.txt + @expr 2 \* \( $$(grep . __W.txt) / 2 \) > __W1.txt + @echo "$$(grep . __W1.txt)!x$$(grep . __H1.txt)!" > __resize.txt + @magick mogrify -format jpg -resize $$(grep . __resize.txt) $(OUTPUT)/s*.svg + rm -f __H*.txt __W*.txt __resize.txt + +gif: + magick convert $(OUTPUT)/s*.svg $(OUTPUT)/out.gif + +movie: + ffmpeg -r $(FRAMERATE) -f image2 -i $(OUTPUT)/snapshot%08d.jpg -vcodec libx264 -pix_fmt yuv420p -strict -2 -tune animation -crf 15 -acodec none $(OUTPUT)/out.mp4 + +# upgrade rules + +SOURCE := PhysiCell_upgrade.zip +get-upgrade: + @echo $$(curl https://raw.githubusercontent.com/MathCancer/PhysiCell/master/VERSION.txt) > VER.txt + @echo https://github.com/MathCancer/PhysiCell/releases/download/$$(grep . VER.txt)/PhysiCell_V.$$(grep . VER.txt).zip > DL_FILE.txt + rm -f VER.txt + $$(curl -L $$(grep . DL_FILE.txt) --output PhysiCell_upgrade.zip) + rm -f DL_FILE.txt + +PhysiCell_upgrade.zip: + make get-upgrade + +upgrade: $(SOURCE) + unzip $(SOURCE) PhysiCell/VERSION.txt + mv -f PhysiCell/VERSION.txt . + unzip $(SOURCE) PhysiCell/core/* + cp -r PhysiCell/core/* core + unzip $(SOURCE) PhysiCell/modules/* + cp -r PhysiCell/modules/* modules + unzip $(SOURCE) PhysiCell/sample_projects/* + cp -r PhysiCell/sample_projects/* sample_projects + unzip $(SOURCE) PhysiCell/BioFVM/* + cp -r PhysiCell/BioFVM/* BioFVM + unzip $(SOURCE) PhysiCell/documentation/User_Guide.pdf + mv -f PhysiCell/documentation/User_Guide.pdf documentation + rm -f -r PhysiCell + rm -f $(SOURCE) + +# use: make save PROJ=your_project_name +PROJ := my_project + +save: + echo "Saving project as $(PROJ) ... " + mkdir -p ./user_projects + mkdir -p ./user_projects/$(PROJ) + mkdir -p ./user_projects/$(PROJ)/custom_modules + mkdir -p ./user_projects/$(PROJ)/config + cp main.cpp ./user_projects/$(PROJ) + cp Makefile ./user_projects/$(PROJ) + cp VERSION.txt ./user_projects/$(PROJ) + cp ./config/* ./user_projects/$(PROJ)/config + cp ./custom_modules/* ./user_projects/$(PROJ)/custom_modules + +load: + echo "Loading project from $(PROJ) ... " + cp ./user_projects/$(PROJ)/main.cpp . + cp ./user_projects/$(PROJ)/Makefile . + cp ./user_projects/$(PROJ)/config/* ./config/ + cp ./user_projects/$(PROJ)/custom_modules/* ./custom_modules/ + +pack: + @echo " " + @echo "Preparing project $(PROJ) for sharing ... " + @echo " " + cd ./user_projects && zip -r $(PROJ).zip $(PROJ) + @echo " " + @echo "Share ./user_projects/$(PROJ).zip ... " + @echo "Other users can unzip $(PROJ).zip in their ./user_projects, compile, and run." + @echo " " + +unpack: + @echo " " + @echo "Preparing shared project $(PROJ).zip for use ... " + @echo " " + cd ./user_projects && unzip $(PROJ).zip + @echo " " + @echo "Load this project via make load PROJ=$(PROJ) ... " + @echo " " + +list-user-projects: + @echo "user projects::" + @cd ./user_projects && ls -dt1 * | grep . | sed 's!empty.txt!!' + +run: + rm output/* + ./$(PROGRAM_NAME) \ No newline at end of file diff --git a/sample_projects/extended_asym_div/config/PhysiCell_settings.xml b/sample_projects/extended_asym_div/config/PhysiCell_settings.xml new file mode 100644 index 000000000..92067fd67 --- /dev/null +++ b/sample_projects/extended_asym_div/config/PhysiCell_settings.xml @@ -0,0 +1,797 @@ + + + + -500 + 500 + -500 + 500 + -10 + 10 + 20 + 20 + 20 + true + + + + 7200 + min + micron + 0.01 + 0.1 + 6 + + + + 6 + + + + output + + 60 + true + + + 60 + true + + substrate + YlOrRd + 0 + 1 + + + + false + + + + + false + true + false + 0 + + + + + + 100000.0 + 10 + + 0 + 0 + + 0 + 0 + 0 + 0 + 0 + 0 + + + + true + true + + ./config/initial.mat + + + ./config/dirichlet.mat + + + + + + + + + + 0.005 + + + 1.0 + + + + + 0 + + 516 + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + 0.0 + + 0 + 86400 + + + 1.11667e-2 + 8.33333e-4 + 5.33333e-5 + 2.16667e-3 + 0 + 2.0 + + + + + 2494 + 0.75 + 540 + 0.05 + 0.0045 + 0.0055 + 0 + 0 + 2.0 + + + 0.4 + 10.0 + 1.25 + + 1 + + + 1.8 + 15.12 + + 0.01 + 0.0 + 0.0 + 12 + + + 1 + 1 + .5 + + false + true + + false + substrate + 1 + + + false + false + + 0.0 + + + + + + + 0 + 1 + 0 + 0 + + + + 0 + 0 + 0 + + 0 + + + 0 + + 1 + 0.1 + + 0 + + + + + 0 + + + + 0.0 + 0.0 + + + + 1.0 + + + + Volume + 4 + 2 + 100000 + + + apoptosis + 1e-6 + 1e-2 + + + + + + + + 0.0 + + + + + 0 + + 516 + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + 0.0 + + 0 + 86400 + + + 1.11667e-2 + 8.33333e-4 + 5.33333e-5 + 2.16667e-3 + 0 + 2.0 + + + + + 2494 + 0.75 + 540 + 0.05 + 0.0045 + 0.0055 + 0 + 0 + 2.0 + + + 0.4 + 10.0 + 1.25 + + 1 + + + 1.8 + 15.12 + + 0.01 + 0.0 + 0.0 + 12 + + + 1 + 1 + .5 + + false + true + + false + substrate + 1 + + + false + false + + 0.0 + + + + + + + 0 + 1 + 0 + 0 + + + + 0 + 0 + 0 + + 0 + + + 0 + + 1 + 0.1 + + 0 + + + + + 0 + + + + 0.0 + 0.0 + + + + 1.0 + + + + Volume + 4 + 2 + 100000 + + + apoptosis + 1e-6 + 1e-2 + + + + + + + + 0 + + + + + 0 + + 516 + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + 0.0 + + 0 + 86400 + + + 1.11667e-2 + 8.33333e-4 + 5.33333e-5 + 2.16667e-3 + 0 + 2.0 + + + + + 2494 + 0.75 + 540 + 0.05 + 0.0045 + 0.0055 + 0 + 0 + 2.0 + + + 0.4 + 10.0 + 1.25 + + 1 + + + 1.8 + 15.12 + + 0.01 + 0.0 + 0.0 + 12 + + + 1 + 1 + .5 + + false + true + + false + substrate + 1 + + + false + false + + 0.0 + + + + + + + 0 + 1 + 0 + 0 + + + + 0 + 0 + 0 + + 0 + + + 0 + + 1 + 0.1 + + 0 + + + + + 0 + + + + 0.0 + 0.0 + + + + 1.0 + + + + Volume + 4 + 2 + 100000 + + + apoptosis + 1e-6 + 1e-2 + + + + + + + + 0 + + + + + 0 + + 516 + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + 0.0 + + 0 + 86400 + + + 1.11667e-2 + 8.33333e-4 + 5.33333e-5 + 2.16667e-3 + 0 + 2.0 + + + + + 2494 + 0.75 + 540 + 0.05 + 0.0045 + 0.0055 + 0 + 0 + 2.0 + + + 0.4 + 10.0 + 1.25 + + 1 + + + 1.8 + 15.12 + + 0.01 + 0.0 + 0.0 + 12 + + + 1 + 1 + .5 + + false + true + + false + substrate + 1 + + + false + false + + 0.0 + + + + + + + 0 + 1 + 0 + 0 + + + + 0 + 0 + 0 + + 0 + + + 0 + + 1 + 0.1 + + 0 + + + + + 0 + + + + 0.0 + 0.0 + + + + 1.0 + + + + Volume + 4 + 2 + 100000 + + + apoptosis + 1e-6 + 1e-2 + + + + + + + + 0 + + + + + 0 + + 516 + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + 0.0 + + 0 + 86400 + + + 1.11667e-2 + 8.33333e-4 + 5.33333e-5 + 2.16667e-3 + 0 + 2.0 + + + + + 2494 + 0.75 + 540 + 0.05 + 0.0045 + 0.0055 + 0 + 0 + 2.0 + + + 0.4 + 10.0 + 1.25 + + 1 + + + 1.8 + 15.12 + + 0.01 + 0.0 + 0.0 + 12 + + + 1 + 1 + .5 + + false + true + + false + substrate + 1 + + + false + false + + 0.0 + + + + + + + 0 + 1 + 0 + 0 + + + + 0 + 0 + 0 + + 0 + + + 0 + + 1 + 0.1 + + 0 + + + + + 0 + + + + 0.0 + 0.0 + + + + 1.0 + + + + Volume + 4 + 2 + 100000 + + + apoptosis + 1e-6 + 1e-2 + + + + + + + + ./config + cells.csv + + + + + + + ./config + cell_rules.csv + + + + + + + 0 + + \ No newline at end of file diff --git a/sample_projects/extended_asym_div/config/cell_rules.csv b/sample_projects/extended_asym_div/config/cell_rules.csv new file mode 100644 index 000000000..3f54e3355 --- /dev/null +++ b/sample_projects/extended_asym_div/config/cell_rules.csv @@ -0,0 +1,3 @@ +type0,time,increases,extended asymmetric division to type1 and type2,0.5,600,256,0 +type0,time,decreases,extended asymmetric division to type1 and type2,0.0,1200,256,0 +type0,time,increases,extended asymmetric division to type3 and type4,0.5,1200,256,0 \ No newline at end of file diff --git a/sample_projects/extended_asym_div/config/cells.csv b/sample_projects/extended_asym_div/config/cells.csv new file mode 100644 index 000000000..5684724fb --- /dev/null +++ b/sample_projects/extended_asym_div/config/cells.csv @@ -0,0 +1,2 @@ +x,y,z,type,volume,cycle entry,custom:GFP,custom:sample +0,0,0,type0,,,, \ No newline at end of file diff --git a/sample_projects/extended_asym_div/custom_modules/custom.cpp b/sample_projects/extended_asym_div/custom_modules/custom.cpp new file mode 100644 index 000000000..59b23bd93 --- /dev/null +++ b/sample_projects/extended_asym_div/custom_modules/custom.cpp @@ -0,0 +1,214 @@ +/* +############################################################################### +# If you use PhysiCell in your project, please cite PhysiCell and the version # +# number, such as below: # +# # +# We implemented and solved the model using PhysiCell (Version x.y.z) [1]. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# See VERSION.txt or call get_PhysiCell_version() to get the current version # +# x.y.z. Call display_citations() to get detailed information on all cite-# +# able software used in your PhysiCell application. # +# # +# Because PhysiCell extensively uses BioFVM, we suggest you also cite BioFVM # +# as below: # +# # +# We implemented and solved the model using PhysiCell (Version x.y.z) [1], # +# with BioFVM [2] to solve the transport equations. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# [2] A Ghaffarizadeh, SH Friedman, and P Macklin, BioFVM: an efficient para- # +# llelized diffusive transport solver for 3-D biological simulations, # +# Bioinformatics 32(8): 1256-8, 2016. DOI: 10.1093/bioinformatics/btv730 # +# # +############################################################################### +# # +# BSD 3-Clause License (see https://opensource.org/licenses/BSD-3-Clause) # +# # +# Copyright (c) 2015-2021, Paul Macklin and the PhysiCell Project # +# All rights reserved. # +# # +# Redistribution and use in source and binary forms, with or without # +# modification, are permitted provided that the following conditions are met: # +# # +# 1. Redistributions of source code must retain the above copyright notice, # +# this list of conditions and the following disclaimer. # +# # +# 2. 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. # +# # +# 3. Neither the name of the copyright holder 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 HOLDER 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. # +# # +############################################################################### +*/ + +#include "./custom.h" + +void create_cell_types( void ) +{ + // set the random seed + if (parameters.ints.find_index("random_seed") != -1) + { + SeedRandom(parameters.ints("random_seed")); + } + + /* + Put any modifications to default cell definition here if you + want to have "inherited" by other cell types. + + This is a good place to set default functions. + */ + + initialize_default_cell_definition(); + cell_defaults.phenotype.secretion.sync_to_microenvironment( µenvironment ); + + cell_defaults.functions.volume_update_function = standard_volume_update_function; + cell_defaults.functions.update_velocity = standard_update_cell_velocity; + + cell_defaults.functions.update_migration_bias = NULL; + cell_defaults.functions.update_phenotype = NULL; // update_cell_and_death_parameters_O2_based; + cell_defaults.functions.custom_cell_rule = NULL; + cell_defaults.functions.contact_function = NULL; + + cell_defaults.functions.add_cell_basement_membrane_interactions = NULL; + cell_defaults.functions.calculate_distance_to_membrane = NULL; + + /* + This parses the cell definitions in the XML config file. + */ + + initialize_cell_definitions_from_pugixml(); + + /* + This builds the map of cell definitions and summarizes the setup. + */ + + build_cell_definitions_maps(); + + /* + This intializes cell signal and response dictionaries + */ + + setup_signal_behavior_dictionaries(); + + /* + Cell rule definitions + */ + + setup_cell_rules(); + + /* + Put any modifications to individual cell definitions here. + + This is a good place to set custom functions. + */ + + cell_defaults.functions.update_phenotype = phenotype_function; + cell_defaults.functions.custom_cell_rule = custom_function; + cell_defaults.functions.contact_function = contact_function; + + /* + This builds the map of cell definitions and summarizes the setup. + */ + + display_cell_definitions( std::cout ); + + return; +} + +void setup_microenvironment( void ) +{ + // set domain parameters + + // put any custom code to set non-homogeneous initial conditions or + // extra Dirichlet nodes here. + + // initialize BioFVM + + initialize_microenvironment(); + + return; +} + +void setup_tissue( void ) +{ + double Xmin = microenvironment.mesh.bounding_box[0]; + double Ymin = microenvironment.mesh.bounding_box[1]; + double Zmin = microenvironment.mesh.bounding_box[2]; + + double Xmax = microenvironment.mesh.bounding_box[3]; + double Ymax = microenvironment.mesh.bounding_box[4]; + double Zmax = microenvironment.mesh.bounding_box[5]; + + if( default_microenvironment_options.simulate_2D == true ) + { + Zmin = 0.0; + Zmax = 0.0; + } + + double Xrange = Xmax - Xmin; + double Yrange = Ymax - Ymin; + double Zrange = Zmax - Zmin; + + // create some of each type of cell + + Cell* pC; + + for( int k=0; k < cell_definitions_by_index.size() ; k++ ) + { + Cell_Definition* pCD = cell_definitions_by_index[k]; + std::cout << "Placing cells of type " << pCD->name << " ... " << std::endl; + for( int n = 0 ; n < parameters.ints("number_of_cells") ; n++ ) + { + std::vector position = {0,0,0}; + position[0] = Xmin + UniformRandom()*Xrange; + position[1] = Ymin + UniformRandom()*Yrange; + position[2] = Zmin + UniformRandom()*Zrange; + + pC = create_cell( *pCD ); + pC->assign_position( position ); + } + } + std::cout << std::endl; + + // load cells from your CSV file (if enabled) + load_cells_from_pugixml(); + set_parameters_from_distributions(); + + return; +} + +std::vector my_coloring_function( Cell* pCell ) +{ return paint_by_number_cell_coloring(pCell); } + +void phenotype_function( Cell* pCell, Phenotype& phenotype, double dt ) +{ return; } + +void custom_function( Cell* pCell, Phenotype& phenotype , double dt ) +{ return; } + +void contact_function( Cell* pMe, Phenotype& phenoMe , Cell* pOther, Phenotype& phenoOther , double dt ) +{ return; } \ No newline at end of file diff --git a/sample_projects/extended_asym_div/custom_modules/custom.h b/sample_projects/extended_asym_div/custom_modules/custom.h new file mode 100644 index 000000000..0e6df8d02 --- /dev/null +++ b/sample_projects/extended_asym_div/custom_modules/custom.h @@ -0,0 +1,92 @@ +/* +############################################################################### +# If you use PhysiCell in your project, please cite PhysiCell and the version # +# number, such as below: # +# # +# We implemented and solved the model using PhysiCell (Version x.y.z) [1]. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# See VERSION.txt or call get_PhysiCell_version() to get the current version # +# x.y.z. Call display_citations() to get detailed information on all cite-# +# able software used in your PhysiCell application. # +# # +# Because PhysiCell extensively uses BioFVM, we suggest you also cite BioFVM # +# as below: # +# # +# We implemented and solved the model using PhysiCell (Version x.y.z) [1], # +# with BioFVM [2] to solve the transport equations. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# [2] A Ghaffarizadeh, SH Friedman, and P Macklin, BioFVM: an efficient para- # +# llelized diffusive transport solver for 3-D biological simulations, # +# Bioinformatics 32(8): 1256-8, 2016. DOI: 10.1093/bioinformatics/btv730 # +# # +############################################################################### +# # +# BSD 3-Clause License (see https://opensource.org/licenses/BSD-3-Clause) # +# # +# Copyright (c) 2015-2021, Paul Macklin and the PhysiCell Project # +# All rights reserved. # +# # +# Redistribution and use in source and binary forms, with or without # +# modification, are permitted provided that the following conditions are met: # +# # +# 1. Redistributions of source code must retain the above copyright notice, # +# this list of conditions and the following disclaimer. # +# # +# 2. 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. # +# # +# 3. Neither the name of the copyright holder 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 HOLDER 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. # +# # +############################################################################### +*/ + +#include "../core/PhysiCell.h" +#include "../modules/PhysiCell_standard_modules.h" + +using namespace BioFVM; +using namespace PhysiCell; + +// setup functions to help us along + +void create_cell_types( void ); +void setup_tissue( void ); + +// set up the BioFVM microenvironment +void setup_microenvironment( void ); + +// custom pathology coloring function + +std::vector my_coloring_function( Cell* ); + +// custom functions can go here + +void phenotype_function( Cell* pCell, Phenotype& phenotype, double dt ); +void custom_function( Cell* pCell, Phenotype& phenotype , double dt ); + +void contact_function( Cell* pMe, Phenotype& phenoMe , Cell* pOther, Phenotype& phenoOther , double dt ); + diff --git a/sample_projects/extended_asym_div/main.cpp b/sample_projects/extended_asym_div/main.cpp new file mode 100644 index 000000000..bf6e65d15 --- /dev/null +++ b/sample_projects/extended_asym_div/main.cpp @@ -0,0 +1,254 @@ +/* +############################################################################### +# If you use PhysiCell in your project, please cite PhysiCell and the version # +# number, such as below: # +# # +# We implemented and solved the model using PhysiCell (Version x.y.z) [1]. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# See VERSION.txt or call get_PhysiCell_version() to get the current version # +# x.y.z. Call display_citations() to get detailed information on all cite-# +# able software used in your PhysiCell application. # +# # +# Because PhysiCell extensively uses BioFVM, we suggest you also cite BioFVM # +# as below: # +# # +# We implemented and solved the model using PhysiCell (Version x.y.z) [1], # +# with BioFVM [2] to solve the transport equations. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# [2] A Ghaffarizadeh, SH Friedman, and P Macklin, BioFVM: an efficient para- # +# llelized diffusive transport solver for 3-D biological simulations, # +# Bioinformatics 32(8): 1256-8, 2016. DOI: 10.1093/bioinformatics/btv730 # +# # +############################################################################### +# # +# BSD 3-Clause License (see https://opensource.org/licenses/BSD-3-Clause) # +# # +# Copyright (c) 2015-2022, Paul Macklin and the PhysiCell Project # +# All rights reserved. # +# # +# Redistribution and use in source and binary forms, with or without # +# modification, are permitted provided that the following conditions are met: # +# # +# 1. Redistributions of source code must retain the above copyright notice, # +# this list of conditions and the following disclaimer. # +# # +# 2. 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. # +# # +# 3. Neither the name of the copyright holder 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 HOLDER 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. # +# # +############################################################################### +*/ + +#include +#include +#include +#include +#include +#include +#include + +#include "./core/PhysiCell.h" +#include "./modules/PhysiCell_standard_modules.h" + +// put custom code modules here! + +#include "./custom_modules/custom.h" + +using namespace BioFVM; +using namespace PhysiCell; + +int main( int argc, char* argv[] ) +{ + // load and parse settings file(s) + + bool XML_status = false; + char copy_command [1024]; + if( argc > 1 ) + { + XML_status = load_PhysiCell_config_file( argv[1] ); + sprintf( copy_command , "cp %s %s" , argv[1] , PhysiCell_settings.folder.c_str() ); + } + else + { + XML_status = load_PhysiCell_config_file( "./config/PhysiCell_settings.xml" ); + sprintf( copy_command , "cp ./config/PhysiCell_settings.xml %s" , PhysiCell_settings.folder.c_str() ); + } + if( !XML_status ) + { exit(-1); } + + // copy config file to output directry + system( copy_command ); + + // OpenMP setup + omp_set_num_threads(PhysiCell_settings.omp_num_threads); + + // time setup + std::string time_units = "min"; + + /* Microenvironment setup */ + + setup_microenvironment(); // modify this in the custom code + + /* PhysiCell setup */ + + // set mechanics voxel size, and match the data structure to BioFVM + double mechanics_voxel_size = 30; + Cell_Container* cell_container = create_cell_container_for_microenvironment( microenvironment, mechanics_voxel_size ); + + /* Users typically start modifying here. START USERMODS */ + + create_cell_types(); + + setup_tissue(); + + /* Users typically stop modifying here. END USERMODS */ + + // set MultiCellDS save options + + set_save_biofvm_mesh_as_matlab( true ); + set_save_biofvm_data_as_matlab( true ); + set_save_biofvm_cell_data( true ); + set_save_biofvm_cell_data_as_custom_matlab( true ); + + // save a simulation snapshot + + char filename[1024]; + sprintf( filename , "%s/initial" , PhysiCell_settings.folder.c_str() ); + save_PhysiCell_to_MultiCellDS_v2( filename , microenvironment , PhysiCell_globals.current_time ); + + // save a quick SVG cross section through z = 0, after setting its + // length bar to 200 microns + + PhysiCell_SVG_options.length_bar = 200; + + // for simplicity, set a pathology coloring function + + std::vector (*cell_coloring_function)(Cell*) = my_coloring_function; + std::string (*substrate_coloring_function)(double, double, double) = paint_by_density_percentage; + + sprintf( filename , "%s/initial.svg" , PhysiCell_settings.folder.c_str() ); + SVG_plot( filename , microenvironment, 0.0 , PhysiCell_globals.current_time, cell_coloring_function, substrate_coloring_function ); + + sprintf( filename , "%s/legend.svg" , PhysiCell_settings.folder.c_str() ); + create_plot_legend( filename , cell_coloring_function ); + + display_citations(); + + // set the performance timers + + BioFVM::RUNTIME_TIC(); + BioFVM::TIC(); + + std::ofstream report_file; + if( PhysiCell_settings.enable_legacy_saves == true ) + { + sprintf( filename , "%s/simulation_report.txt" , PhysiCell_settings.folder.c_str() ); + + report_file.open(filename); // create the data log file + report_file<<"simulated time\tnum cells\tnum division\tnum death\twall time"< PhysiCell_globals.next_full_save_time - 0.5 * diffusion_dt ) + { + display_simulation_status( std::cout ); + if( PhysiCell_settings.enable_legacy_saves == true ) + { + log_output( PhysiCell_globals.current_time , PhysiCell_globals.full_output_index, microenvironment, report_file); + } + + if( PhysiCell_settings.enable_full_saves == true ) + { + sprintf( filename , "%s/output%08u" , PhysiCell_settings.folder.c_str(), PhysiCell_globals.full_output_index ); + + save_PhysiCell_to_MultiCellDS_v2( filename , microenvironment , PhysiCell_globals.current_time ); + } + + PhysiCell_globals.full_output_index++; + PhysiCell_globals.next_full_save_time += PhysiCell_settings.full_save_interval; + } + + // save SVG plot if it's time + if( PhysiCell_globals.current_time > PhysiCell_globals.next_SVG_save_time - 0.5 * diffusion_dt ) + { + if( PhysiCell_settings.enable_SVG_saves == true ) + { + sprintf( filename , "%s/snapshot%08u.svg" , PhysiCell_settings.folder.c_str() , PhysiCell_globals.SVG_output_index ); + SVG_plot( filename , microenvironment, 0.0 , PhysiCell_globals.current_time, cell_coloring_function, substrate_coloring_function); + + PhysiCell_globals.SVG_output_index++; + PhysiCell_globals.next_SVG_save_time += PhysiCell_settings.SVG_save_interval; + } + } + + // update the microenvironment + microenvironment.simulate_diffusion_decay( diffusion_dt ); + + // run PhysiCell + ((Cell_Container *)microenvironment.agent_container)->update_all_cells( PhysiCell_globals.current_time ); + + /* + Custom add-ons could potentially go here. + */ + + PhysiCell_globals.current_time += diffusion_dt; + } + + if( PhysiCell_settings.enable_legacy_saves == true ) + { + log_output(PhysiCell_globals.current_time, PhysiCell_globals.full_output_index, microenvironment, report_file); + report_file.close(); + } + } + catch( const std::exception& e ) + { // reference to the base of a polymorphic object + std::cout << e.what(); // information from length_error printed + } + + // save a final simulation snapshot + + sprintf( filename , "%s/final" , PhysiCell_settings.folder.c_str() ); + save_PhysiCell_to_MultiCellDS_v2( filename , microenvironment , PhysiCell_globals.current_time ); + + sprintf( filename , "%s/final.svg" , PhysiCell_settings.folder.c_str() ); + SVG_plot(filename, microenvironment, 0.0, PhysiCell_globals.current_time, cell_coloring_function, substrate_coloring_function); + + // timer + + std::cout << std::endl << "Total simulation runtime: " << std::endl; + BioFVM::display_stopwatch_value( std::cout , BioFVM::runtime_stopwatch_value() ); + + return 0; +}