From 5ee31233fd724f8f939c025ceeea9c3f08249e9c Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Mon, 21 Aug 2023 19:21:01 -0400 Subject: [PATCH] added ability to initialize BioFVM substrates with initial conditions specified in a .mat file CSVs (w/o headers) can be used for substrate ics. Substrate ICs from csv with headers. Now can include a header row in the CSV defining substrate ICs so that only a subset of substrates can be included for specification. If no header row is given, then it assumes the first n substrates are specified, which could be problematic? It issues a warning about that though. --- BioFVM/BioFVM_matlab.cpp | 3 +- BioFVM/BioFVM_microenvironment.cpp | 206 ++++++++++++++++++++++++++++- BioFVM/BioFVM_microenvironment.h | 8 ++ modules/PhysiCell_settings.cpp | 15 ++- 4 files changed, 223 insertions(+), 9 deletions(-) diff --git a/BioFVM/BioFVM_matlab.cpp b/BioFVM/BioFVM_matlab.cpp index ed056f45f..f64f72936 100644 --- a/BioFVM/BioFVM_matlab.cpp +++ b/BioFVM/BioFVM_matlab.cpp @@ -108,7 +108,8 @@ std::vector< std::vector > read_matlab( std::string filename ) type_data_format > 5 || // unknown format type_matrix_type != 0 ) // want full matrices, not sparse { - std::cout << "Error reading file " << filename << ": I can't read this format yet!" << std::endl; + std::cout << "Error reading file " << filename << ": I can't read this format yet!" << std::endl + << "\t Make sure you save the .mat file in '-v4'" << std::endl; return output; } diff --git a/BioFVM/BioFVM_microenvironment.cpp b/BioFVM/BioFVM_microenvironment.cpp index 1c330870f..a0e245ea9 100644 --- a/BioFVM/BioFVM_microenvironment.cpp +++ b/BioFVM/BioFVM_microenvironment.cpp @@ -1319,9 +1319,29 @@ void initialize_microenvironment( void ) default_microenvironment_options.initial_condition_vector = default_microenvironment_options.Dirichlet_condition_vector; } - // set the initial condition - for( unsigned int n=0; n < microenvironment.number_of_voxels() ; n++ ) - { microenvironment.density_vector(n) = default_microenvironment_options.initial_condition_vector; } + // set the initial condition + + if (default_microenvironment_options.initial_condition_from_file_enabled) + { + if (default_microenvironment_options.initial_condition_file_type=="matlab") + { + load_initial_conditions_from_matlab(default_microenvironment_options.initial_condition_file); + } + else if (default_microenvironment_options.initial_condition_file_type=="csv" || default_microenvironment_options.initial_condition_file_type=="CSV") + { + load_initial_conditions_from_csv(default_microenvironment_options.initial_condition_file); + } + else // eventually can add other file types + { + std::cout << "ERROR : Load BioFVM initial conditions from " << default_microenvironment_options.initial_condition_file_type << " not yet supported." << std::endl; + exit(-1); + } + } + else // do what was done before + { + for (unsigned int n = 0; n < microenvironment.number_of_voxels(); n++) + { microenvironment.density_vector(n) = default_microenvironment_options.initial_condition_vector; } + } // now, figure out which sides have BCs (for at least one substrate): @@ -1531,14 +1551,188 @@ void initialize_microenvironment( void ) // April 2023: no longer necessary after flipping our approach and doing an "additive" instead of "subtractive" DCs handling. I.e., we assume DC activation is false by default; make true on-demand. - // // set the Dirichlet condition activation vector to match the microenvironment options + // // set the Dirichlet condition activation vector to match the microenvironment options // for( int i=0 ; i < default_microenvironment_options.Dirichlet_activation_vector.size(); i++ ) // { - // microenvironment.set_substrate_dirichlet_activation( i , default_microenvironment_options.Dirichlet_activation_vector[i] ); + // microenvironment.set_substrate_dirichlet_activation( i , default_microenvironment_options.Dirichlet_activation_vector[i] ); // } + + microenvironment.display_information(std::cout); + return; +} + +void load_initial_conditions_from_matlab(std::string filename) +{ + // The .mat file needs to contain exactly one matrix where each row corresponds to a single voxel. + // Each row is a vector of values as follows: [x coord, y coord, z coord, substrate id 0 value, substrate id 1 value, ...] + // Thus, your matrix should be of size #voxels x (3 + #densities) (rows x columns) + // M will be read in such that each element of M is one of these rows + std::vector< std::vector > M = read_matlab( filename ); + int num_mat_voxels = M.size(); + int num_mat_features = M[1].size(); + if (num_mat_voxels != microenvironment.number_of_voxels()) + { + std::cout << "ERROR : Wrong number of voxels supplied in the .mat file specifying BioFVM initial conditions." << std::endl + << "\tExpected: " << microenvironment.number_of_voxels() << std::endl + << "\tFound: " << num_mat_voxels << std::endl + << "\tRemember, save your matrix as a #voxels x (3 + #densities) matrix." << std::endl; + exit(-1); + } + if (num_mat_features != (microenvironment.number_of_densities() + 3)) + { + std::cout << "ERROR : Wrong number of density values supplied in the .mat file specifying BioFVM initial conditions." << std::endl + << "\tExpected: " << microenvironment.number_of_densities() << std::endl + << "\tFound: " << num_mat_features - 3 << std::endl + << "\tRemember, save your matrix as a #voxels x (3 + #densities) matrix." << std::endl; + exit(-1); + } + std::vector voxel_set = {}; // set to check that no voxel value is set twice + int voxel_ind; + std::vector position; + for (unsigned int i = 0; i < num_mat_voxels; i++) + { + position = {M[i][0], M[i][1], M[i][2]}; + voxel_ind = microenvironment.mesh.nearest_voxel_index(position); + for (unsigned int ci = 0; ci < microenvironment.number_of_densities(); ci++) + { + microenvironment.density_vector(voxel_ind)[ci] = M[i][ci+3]; + } + for (unsigned int j = 0; j < i; j++) + { + if (voxel_ind==voxel_set[j]) + { + std::cout << "ERROR : the matlab-supplied initial conditions for BioFVM repeat the same voxel. Fix the .mat file and try again." << std::endl + << "\tPosition that was repeated: " << position << std::endl; + exit(-1); + } + } + voxel_set.push_back(voxel_ind); + } + + return; +} + +void load_initial_conditions_from_csv(std::string filename) +{ + // The .csv file needs to contain one row per voxel. + // Each row is a vector of values as follows: [x coord, y coord, z coord, substrate id 0 value, substrate id 1 value, ...] + // Thus, your table should be of size #voxels x (3 + #densities) (rows x columns) + // Do not include a header row. + + // open file + std::ifstream file( filename, std::ios::in ); + if( !file ) + { + std::cout << "ERROR: " << filename << " not found during cell loading. Quitting." << std::endl; + exit(-1); + } + + // determine if header row exists + std::string line; + std::getline( file , line ); + char c = line.c_str()[0]; + std::vector substrate_indices; + bool header_provided = false; + if( c == 'X' || c == 'x' ) + { + // do not support this with a header yet + if ((line.c_str()[2] != 'Y' && line.c_str()[2] != 'y') || (line.c_str()[4] != 'Z' && line.c_str()[4] != 'z')) + { + std::cout << "ERROR: Header row starts with x but then not y,z? What is this? Exiting now." << std::endl; + file.close(); + exit(-1); + } + std::vector< std::string> column_names; // this will include x,y,z (so make sure to skip those below) + std::stringstream stream(line); + std::string field; + + while (std::getline(stream, field, ',')) + { + column_names.push_back(field); + } + for (int i = 3; i voxel_set = {}; // set to check that no voxel value is set twice + + while (std::getline(file, line)) + { + get_row_from_substrate_initial_condition_csv(voxel_set, line, substrate_indices, header_provided); + } + + if (voxel_set.size() != microenvironment.number_of_voxels()) + { + std::cout << "ERROR : Wrong number of voxels supplied in the .csv file specifying BioFVM initial conditions." << std::endl + << "\tExpected: " << microenvironment.number_of_voxels() << std::endl + << "\tFound: " << voxel_set.size() << std::endl + << "\tRemember, your table should have dimensions #voxels x (3 + #densities)." << std::endl; + exit(-1); + } + + file.close(); - microenvironment.display_information( std::cout ); return; } +void get_row_from_substrate_initial_condition_csv(std::vector &voxel_set, const std::string line, const std::vector substrate_indices, const bool header_provided) +{ + static bool warning_issued = false; + std::vector data; + csv_to_vector(line.c_str(), data); + + if (!(warning_issued) && !(header_provided) && (data.size() != (microenvironment.number_of_densities() + 3))) + { + std::cout << "WARNING: Wrong number of density values supplied in the .csv file specifying BioFVM initial conditions." << std::endl + << "\tExpected: " << microenvironment.number_of_voxels() << std::endl + << "\tFound: " << data.size() - 3 << std::endl + << "\tRemember, save your csv with columns as: x, y, z, substrate_0, substrate_1,...." << std::endl + << "\tThis could also be resolved by including a header row \"x,y,z,[substrate_i0,substrate_i1]\"" << std::endl; + warning_issued = true; + } + + std::vector position = {data[0], data[1], data[2]}; + int voxel_ind = microenvironment.mesh.nearest_voxel_index(position); + for (unsigned int ci = 0; ci < substrate_indices.size(); ci++) // column index, counting from the first substrate (or just the index of the vector substrate_indices) + { + microenvironment.density_vector(voxel_ind)[substrate_indices[ci]] = data[ci + 3]; + } + for (unsigned int j = 0; j < voxel_set.size(); j++) + { + if (voxel_ind == voxel_set[j]) + { + std::cout << "ERROR : the csv-supplied initial conditions for BioFVM repeat the same voxel. Fix the .csv file and try again." << std::endl + << "\tPosition that was repeated: " << position << std::endl; + exit(-1); + } + } + voxel_set.push_back(voxel_ind); +} }; diff --git a/BioFVM/BioFVM_microenvironment.h b/BioFVM/BioFVM_microenvironment.h index 4896b44ad..8b894e880 100644 --- a/BioFVM/BioFVM_microenvironment.h +++ b/BioFVM/BioFVM_microenvironment.h @@ -49,6 +49,7 @@ #ifndef __BioFVM_microenvironment_h__ #define __BioFVM_microenvironment_h__ +#include #include "BioFVM_mesh.h" #include "BioFVM_agent_container.h" #include "BioFVM_MultiCellDS.h" @@ -339,6 +340,10 @@ class Microenvironment_Options std::vector Dirichlet_zmax_values; std::vector initial_condition_vector; + + bool initial_condition_from_file_enabled; + std::string initial_condition_file_type; + std::string initial_condition_file; bool simulate_2D; std::vector X_range; @@ -359,6 +364,9 @@ extern Microenvironment microenvironment; void initialize_microenvironment( void ); +void load_initial_conditions_from_matlab( std::string filename ); +void load_initial_conditions_from_csv( std::string filename ); +void get_row_from_substrate_initial_condition_csv(std::vector &voxel_set, const std::string line, const std::vector substrate_indices, const bool header_provided); }; #endif diff --git a/modules/PhysiCell_settings.cpp b/modules/PhysiCell_settings.cpp index 5dfec6a68..0b47d73dd 100644 --- a/modules/PhysiCell_settings.cpp +++ b/modules/PhysiCell_settings.cpp @@ -911,8 +911,19 @@ bool setup_microenvironment_from_XML( pugi::xml_node root_node ) // track internalized substrates in each agent? default_microenvironment_options.track_internalized_substrates_in_each_agent - = xml_get_bool_value( node, "track_internalized_substrates_in_each_agent" ); - + = xml_get_bool_value( node, "track_internalized_substrates_in_each_agent" ); + + node = xml_find_node(node, "initial_condition"); + if (node) + { + default_microenvironment_options.initial_condition_from_file_enabled = node.attribute("enabled").as_bool(); + if (default_microenvironment_options.initial_condition_from_file_enabled) + { + default_microenvironment_options.initial_condition_file_type = node.attribute("type").as_string(); + default_microenvironment_options.initial_condition_file = xml_get_string_value(node, "filename"); + } + } + // not yet supported : read initial conditions /* // read in initial conditions from an external file