diff --git a/BioFVM/BioFVM_microenvironment.cpp b/BioFVM/BioFVM_microenvironment.cpp index 84c975ed6..70e434004 100644 --- a/BioFVM/BioFVM_microenvironment.cpp +++ b/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 split_substrate_csv_row(const std::string &line) +{ + std::vector 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,,..." 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,77 +1497,99 @@ 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 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 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 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 voxel_set = {}; // set to check that no voxel value is set twice - + std::vector 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(); @@ -1543,37 +1597,69 @@ void load_initial_conditions_from_csv(std::string filename) 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) +void get_row_from_substrate_initial_condition_csv(std::vector &voxel_is_set, const std::string &line, const std::vector &substrate_indices, const unsigned int line_number) { - static bool warning_issued = false; - std::vector 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 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 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 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; + } } + }; diff --git a/BioFVM/BioFVM_microenvironment.h b/BioFVM/BioFVM_microenvironment.h index 6fec93728..98759ec56 100644 --- a/BioFVM/BioFVM_microenvironment.h +++ b/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 &voxel_set, const std::string line, const std::vector substrate_indices, const bool header_provided); +void get_row_from_substrate_initial_condition_csv(std::vector &voxel_is_set, const std::string &line, const std::vector &substrate_indices, const unsigned int line_number); }; #endif diff --git a/BioFVM/BioFVM_vector.cpp b/BioFVM/BioFVM_vector.cpp index e12a3672d..6a183b1db 100644 --- a/BioFVM/BioFVM_vector.cpp +++ b/BioFVM/BioFVM_vector.cpp @@ -47,6 +47,8 @@ */ #include "BioFVM_vector.h" +#include +#include /* some global BioFVM strings */ @@ -376,6 +378,53 @@ void csv_to_vector( const char* buffer , std::vector& 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& 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::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& vect ) { static int datum_size = 16; // format = %.7e, 1 (sign) + 1 (lead) + 1 (decimal) + 7 (figs) + 2 (e, sign) + 3 (exponent) + 1 (delimiter) = 16 diff --git a/BioFVM/BioFVM_vector.h b/BioFVM/BioFVM_vector.h index e3c36d449..1b9507ae5 100644 --- a/BioFVM/BioFVM_vector.h +++ b/BioFVM/BioFVM_vector.h @@ -129,7 +129,11 @@ void double_axpy_div( std::vector* y, std::vector& 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& vect ); +void csv_to_vector( const char* buffer , std::vector& 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& vect); + char* vector_to_csv( const std::vector& vect ); void vector_to_csv_safe( const std::vector& vect , char*& buffer ); void vector_to_csv( const std::vector& vect , char*& buffer );