Feature extended asym div

.gitignore

@@ -21,3 +21,4 @@ user_projects/*
 !user_projects/empty.txt
 Studio.zip
 /studio
+.vscode

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/

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: 

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<int, int> 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<int, int> 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<int, int> 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<std::string,int> 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<int, int> 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<int, int> 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<std::string,int> 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<int, int> 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<int, int> extended_asym_index_to_upper_triangle(int index)
+{
+	static std::vector< std::pair<int, int> > pairs_vector = initialize_pairs_vector();
+	return pairs_vector[index];
+}
+
+std::vector< std::pair<int, int> > initialize_pairs_vector()
 {
 	extern std::unordered_map<std::string,int> cell_definition_indices_by_name; 
-	int n = cell_definition_indices_by_name[type_name]; 
-	return asymmetric_division_probabilities[n]; 
+	std::vector< std::pair<int, int> > 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 

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<int, int>& 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<int>{}(UT);
+		return hash; 
+	}
+};
+
+struct equality_function {
+	bool operator()(const std::pair<int, int>& lhs, const std::pair<int, int>& 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<double> asymmetric_division_probabilities; 
+	std::unordered_map<std::pair<int, int>, 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<int, int> 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<int, int> 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<int, int> select_daughter_types(int type_1, int type_2);
 };
 
+std::pair<int, int> extended_asym_index_to_upper_triangle(int index);
+std::vector<std::pair<int, int>> initialize_pairs_vector(void);
+
 class Cycle
 {
  private: