improve substrate csv initial conditions: omit voxels, fix out-of-bounds reads, revive headerless path

BioFVM/BioFVM_microenvironment.cpp

@@ -1450,12 +1450,44 @@ void load_initial_conditions_from_matlab(std::string filename)
 	return;
 }
 
+// split a csv row into trimmed fields. A row of n commas always yields n+1 fields, matching how
+// substrate_csv_to_vector counts them, so a header and its data rows are measured the same way.
+static std::vector<std::string> split_substrate_csv_row(const std::string &line)
+{
+	std::vector<std::string> fields;
+	std::size_t start = 0;
+	while (true)
+	{
+		std::size_t comma = line.find(',', start);
+		std::size_t end = (comma == std::string::npos) ? line.size() : comma;
+		std::size_t first = line.find_first_not_of(" \t\r", start);
+		std::string field; // stays empty when the field is empty or all whitespace
+		if (first != std::string::npos && first < end)
+		{
+			std::size_t last = line.find_last_not_of(" \t\r", end - 1);
+			field = line.substr(first, last - first + 1);
+		}
+		fields.push_back(field);
+		if (comma == std::string::npos)
+		{ return fields; }
+		start = comma + 1;
+	}
+}
+
+static bool is_csv_label(const std::string &field, const char lower, const char upper)
+{ return field.size() == 1 && (field[0] == lower || field[0] == upper); }
+
 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.
+	// Each row locates one voxel and sets substrate densities in it:
+	//     [x coord, y coord, z coord, value, value, ...]
+	// An optional header row "x,y,z,<substrate name>,..." names the substrates being set, which may be
+	// any subset of the densities, in any order. Without a header, the columns after x,y,z are taken to
+	// be the first n densities in index order.
+	// Rows need not cover every voxel: a voxel with no row keeps the initial condition from the config file.
+	// Within a row, an entry may be omitted by leaving the field empty (e.g. "0,0,0,,3.5"), which sets 0.
+	// Every row must otherwise be well formed -- the same column count as the header (or as the first row),
+	// a position inside the domain, and a finite number in every field that is not empty.
 
 	// open file 
 	std::ifstream file( filename, std::ios::in );
@@ -1465,115 +1497,169 @@ void load_initial_conditions_from_csv(std::string filename)
 		exit(-1);
 	}
 
-	// determine if header row exists 
 	std::string line; 
-	std::getline( file , line );
+	if( !std::getline( file , line ) )
+	{
+		std::cout << "ERROR: " << filename << " is empty." << std::endl
+				  << "\tIt must contain at least one row of substrate initial conditions." << std::endl;
+		file.close();
+		exit(-1);
+	}
 	trim_cr(line);
-	char c = line.c_str()[0];
+
+	// determine if header row exists 
+	std::vector<std::string> first_row = split_substrate_csv_row(line);
+	bool header_provided = is_csv_label(first_row[0], 'x', 'X'); // split always returns at least one field
 	std::vector<int> substrate_indices;
-	bool header_provided = false;
-	if( c == 'X' || c == 'x' )
+
+	if( header_provided )
 	{ 
-		// 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'))
+		if( first_row.size() < 3 || !is_csv_label(first_row[1], 'y', 'Y') || !is_csv_label(first_row[2], 'z', '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, ','))
+		if( first_row.size() < 4 )
 		{
-			column_names.push_back(field);
+			std::cout << "ERROR: The header row of " << filename << " names no substrates." << std::endl
+					  << "\tExpected a row like \"x,y,z,[substrate_i0,substrate_i1]\"" << std::endl;
+			file.close();
+			exit(-1);
 		}
-		for (int i = 3; i<column_names.size(); i++) // skip x,y,z by starting at 3, not 0
+		for (unsigned int i = 3; i<first_row.size(); i++) // skip x,y,z by starting at 3, not 0
 		{
-			int substrate_index = microenvironment.find_density_index(column_names[i]);
+			int substrate_index = microenvironment.find_density_index(first_row[i]);
 			if (substrate_index == -1)
 			{
-				std::cout << "ERROR: Substrate " << column_names[i] << " not found in the BioFVM microenvironment. Exiting now." << std::endl;
+				std::cout << "ERROR: Substrate " << first_row[i] << " not found in the BioFVM microenvironment. Exiting now." << std::endl;
 				file.close();
 				exit(-1);
 			}
-			substrate_indices.push_back(microenvironment.find_density_index(column_names[i]));
+			for (unsigned int j = 0; j < substrate_indices.size(); j++)
+			{
+				if (substrate_indices[j] == substrate_index) // two columns writing one density would silently race
+				{
+					std::cout << "ERROR: Substrate " << first_row[i] << " appears more than once in the header row of " << filename << ". Exiting now." << std::endl;
+					file.close();
+					exit(-1);
+				}
+			}
+			substrate_indices.push_back(substrate_index);
 		}
-		header_provided = true;
 	}
 	else // no column labels given; just assume that the first n substrates are supplied (n = # columns after x,y,z)
 	{
-		std::stringstream stream(line);
-		std::string field;
-		int i = 0;
-		while (std::getline(stream, field, ','))
+		if( first_row.size() < 4 )
+		{
+			std::cout << "ERROR: The first row of " << filename << " has " << first_row.size() << " column(s)." << std::endl
+					  << "\tExpected at least 4: x, y, z and one substrate value." << std::endl;
+			file.close();
+			exit(-1);
+		}
+		unsigned int number_supplied = first_row.size() - 3; // the substrate index is the column index - 3
+		if( number_supplied > microenvironment.number_of_densities() )
 		{
-			if (i<3) {continue;} // skip (x,y,z)
-			substrate_indices.push_back(i-3); // the substrate index is the column index - 3 (since x,y,z are the first 3 columns)
-			i++;
+			std::cout << "ERROR: The first row of " << filename << " supplies " << number_supplied << " density values," << std::endl
+					  << "\tbut the BioFVM microenvironment has only " << microenvironment.number_of_densities() << "." << std::endl
+					  << "\tRemember, save your csv with columns as: x, y, z, substrate_0, substrate_1,...." << std::endl;
+			file.close();
+			exit(-1);
 		}
-		// in this case, we want to read this first line, so close the file and re-open so that we start with this line
-		file.close();
-		std::ifstream file(filename, std::ios::in);
-		std::getline(file, line);
-		trim_cr(line);
+		if( number_supplied != microenvironment.number_of_densities() )
+		{
+			std::cout << "WARNING: " << filename << " supplies " << number_supplied << " of the "
+					  << microenvironment.number_of_densities() << " substrate densities," << std::endl
+					  << "\tso the first " << number_supplied << " are assumed." << std::endl
+					  << "\tThis could be resolved by including a header row \"x,y,z,[substrate_i0,substrate_i1]\"" << std::endl;
+		}
+		for( unsigned int i = 0; i < number_supplied; i++ )
+		{ substrate_indices.push_back(i); }
+
+		// the first row is data, not labels, so rewind and let the loop below read it again
+		file.clear();
+		file.seekg(0, std::ios::beg);
 	}
 
 	std::cout << "Loading substrate initial conditions from CSV file " << filename << " ... " << std::endl;
-	std::vector<int> voxel_set = {}; // set to check that no voxel value is set twice
-
+	std::vector<bool> voxel_is_set(microenvironment.number_of_voxels(), false); // set to check that no voxel value is set twice
+
+	unsigned int line_number = header_provided ? 1 : 0; // only a header leaves the first line already consumed
 	while (std::getline(file, line))
 	{
+		line_number++;
 		trim_cr(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);
+		get_row_from_substrate_initial_condition_csv(voxel_is_set, line, substrate_indices, line_number);
 	}
 
 	file.close(); 	
 
 	return;
 }
 
-void get_row_from_substrate_initial_condition_csv(std::vector<int> &voxel_set, const std::string line, const std::vector<int> substrate_indices, const bool header_provided)
+void get_row_from_substrate_initial_condition_csv(std::vector<bool> &voxel_is_set, const std::string &line, const std::vector<int> &substrate_indices, const unsigned int line_number)
 {
-	static bool warning_issued = false;
-	std::vector<double> data;
-	csv_to_vector(line.c_str(), data);
+	if (line.find_first_not_of(" \t\r") == std::string::npos)
+	{ return; } // skip blank lines
 
-	if (!(warning_issued) && !(header_provided) && (data.size() != (microenvironment.number_of_densities() + 3)))
+	std::vector<double> data;
+	unsigned int bad_field = substrate_csv_to_vector(line.c_str(), data);
+	if (bad_field != 0)
 	{
-		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::cout << "ERROR : Column " << bad_field << " of line " << line_number << " of the .csv file specifying BioFVM initial conditions is not a number." << std::endl
+				  << "\tEvery column must hold a finite number, or nothing at all to omit a substrate value." << std::endl
+				  << "\tOffending row: " << line << std::endl;
+		exit(-1);
 	}
 
-	std::vector<double> 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)
+	// holding every row to the expected column count is also what keeps every data[...] access below in bounds
+	if (data.size() != substrate_indices.size() + 3)
 	{
-		microenvironment.density_vector(voxel_ind)[substrate_indices[ci]] = data[ci + 3];
+		std::cout << "ERROR : Line " << line_number << " of the .csv file specifying BioFVM initial conditions has the wrong number of columns." << std::endl
+				  << "\tExpected: " << substrate_indices.size() + 3 << " (x, y, z and " << substrate_indices.size() << " substrate value(s))" << std::endl
+				  << "\tFound: " << data.size() << std::endl
+				  << "\tTo omit a value, leave the field empty but keep the comma (e.g. x,y,z,,3.5)." << std::endl
+				  << "\tOffending row: " << line << std::endl;
+		exit(-1);
 	}
-	for (unsigned int j = 0; j < voxel_set.size(); j++)
+
+	for (unsigned int i = 0; i < 3; i++) // a position cannot be omitted the way a density value can
 	{
-		if (voxel_ind == voxel_set[j])
+		if (std::isnan(data[i]))
 		{
-			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;
+			std::cout << "ERROR : Line " << line_number << " of the .csv file specifying BioFVM initial conditions omits an x, y or z coordinate." << std::endl
+					  << "\tOffending row: " << line << std::endl;
 			exit(-1);
 		}
 	}
-	voxel_set.push_back(voxel_ind);
+
+	std::vector<double> position = {data[0], data[1], data[2]};
+	if (!microenvironment.mesh.is_position_valid(position[0], position[1], position[2])) // otherwise it would silently snap to an edge voxel
+	{
+		std::cout << "ERROR : Line " << line_number << " of the .csv file specifying BioFVM initial conditions lies outside the microenvironment domain." << std::endl
+				  << "\tPosition: " << position << std::endl
+				  << "\tDomain: x in [" << microenvironment.mesh.bounding_box[0] << ", " << microenvironment.mesh.bounding_box[3] << "]"
+				  << ", y in [" << microenvironment.mesh.bounding_box[1] << ", " << microenvironment.mesh.bounding_box[4] << "]"
+				  << ", z in [" << microenvironment.mesh.bounding_box[2] << ", " << microenvironment.mesh.bounding_box[5] << "]" << std::endl;
+		exit(-1);
+	}
+
+	int voxel_ind = microenvironment.mesh.nearest_voxel_index(position);
+	if (voxel_is_set[voxel_ind])
+	{
+		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
+				  << "\tRepeated on line: " << line_number << std::endl;
+		exit(-1);
+	}
+	voxel_is_set[voxel_ind] = true;
+	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)
+	{
+		double value = data[ci + 3];
+		if (std::isnan(value))
+		{ value = 0.0; } // an entry the row omitted
+		microenvironment.density_vector(voxel_ind)[substrate_indices[ci]] = value;
+	}
 }
+
 };

BioFVM/BioFVM_microenvironment.h

@@ -367,7 +367,7 @@ void set_microenvironment_initial_condition( 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<int> &voxel_set, const std::string line, const std::vector<int> substrate_indices, const bool header_provided);
+void get_row_from_substrate_initial_condition_csv(std::vector<bool> &voxel_is_set, const std::string &line, const std::vector<int> &substrate_indices, const unsigned int line_number);
 };
 
 #endif

BioFVM/BioFVM_vector.cpp

@@ -47,6 +47,8 @@
 */
 
 #include "BioFVM_vector.h" 
+#include <cctype>
+#include <limits>
 
 /* some global BioFVM strings */ 
 
@@ -376,6 +378,53 @@ void csv_to_vector( const char* buffer , std::vector<double>& vect )
 	return; 
 }
 
+// Parse one row of a substrate initial-condition csv into vect. A row of n commas always yields
+// n+1 values, so that a caller can hold every row to the same column count.
+//
+// A field that is empty or all whitespace yields NaN, which lets a caller tell an entry the row
+// omitted from one it set to zero. Every other field must parse completely as a finite number, so
+// that a typo cannot pass silently as a value: "1.5abc", "NA" and "inf" are all rejected.
+//
+// Returns the 1-based index of the first field that is neither, or 0 if the whole row is well formed.
+unsigned int substrate_csv_to_vector(const char* buffer, std::vector<double>& vect)
+{
+	vect.clear();
+
+	const char* start = buffer;
+	unsigned int field_number = 0;
+
+	while (true)
+	{
+		const char* end = start;
+		while (*end != '\0' && *end != ',')
+		{ end++; }
+		field_number++;
+
+		// trim the field, so that " 1.5 " reads as 1.5 and " " reads as omitted
+		const char* first = start;
+		const char* last = end;
+		while (first < last && isspace((unsigned char) *first))
+		{ first++; }
+		while (last > first && isspace((unsigned char) *(last - 1)))
+		{ last--; }
+
+		if (first == last) // the row omitted this entry
+		{ vect.push_back(std::numeric_limits<double>::quiet_NaN()); }
+		else
+		{
+			char* parse_end;
+			double value = strtod(first, &parse_end);
+			if (parse_end != last || !std::isfinite(value)) // trailing garbage, not a number, or not finite
+			{ return field_number; }
+			vect.push_back(value);
+		}
+
+		if (*end == '\0')
+		{ return 0; }
+		start = end + 1;
+	}
+}
+
 char* vector_to_csv( const std::vector<double>& vect )
 { 
 	static int datum_size = 16;  // format = %.7e, 1 (sign) + 1 (lead) + 1 (decimal) + 7 (figs) + 2 (e, sign) + 3 (exponent) + 1 (delimiter) = 16

BioFVM/BioFVM_vector.h

@@ -129,7 +129,11 @@ void double_axpy_div( std::vector<double>* y, std::vector<double>& a1 , std::vec
 
 // turn a delimited character array (e.g., csv) into a vector of doubles
 
-void csv_to_vector( const char* buffer , std::vector<double>& vect ); 
+void csv_to_vector( const char* buffer , std::vector<double>& vect );
+// returns the 1-based index of the first malformed field, or 0 if the row is well formed;
+// an empty field yields NaN so that an omitted entry is distinguishable from a zero
+unsigned int substrate_csv_to_vector(const char* buffer, std::vector<double>& vect);
+
 char* vector_to_csv( const std::vector<double>& vect );
 void vector_to_csv_safe( const std::vector<double>& vect , char*& buffer );
 void vector_to_csv( const std::vector<double>& vect , char*& buffer );