stim/biomodels/cellset.h

#ifndef STIM_CELLSET_H
#define STIM_CELLSET_H
#include <stim/math/vec3.h>
#include <vector>
#include <list>
#include <unordered_map>
#include <fstream>
namespace stim{
class cellset{
private:
	static const char delim = ' ';
protected:
	std::vector<double*> cells;							//vector storing field data for each cell
	std::unordered_map<std::string, size_t> fields;		//unordered map storing field->index information for each field
	size_t ip[3];										//hard code to position indices (for speed)
	void init(){
	}
	/// Initialize fields for a standard cell set containing three points and a radius
	void init_p3(){
		fields.insert(std::pair<std::string, size_t>("x", 0));
		ip[0] = 0;
		fields.insert(std::pair<std::string, size_t>("y", 1));
		ip[1] = 1;
		fields.insert(std::pair<std::string, size_t>("z", 2));
		ip[2] = 2;
	}
	void init_p3r(){
		init_p3();
		fields.insert(std::pair<std::string, size_t>("radius", 3));
		ip[3] = 3;
	}
public:
	/// Constructor - create an empty cell set
	cellset(){
		init();
	}
	/// Constructor - load a cellset from a file
	cellset(std::string filename){
		init();											//initialize an empty cellset
		load(filename);									//load the cellset from an existing file
	}
	
	/// Loads a cellset from a file
	void load(std::string filename){
		std::ifstream infile(filename);
		std::string header;								//allocate space for the file header
		std::getline(infile, header);					//get the file header
		// break the header into fields
		std::stringstream ss(header);					//create a string stream
		std::string field;								//store a single field name
		size_t i = 0;									//current field index
		while (std::getline(ss, field, delim)) {		//split the header into individual fields
			std::pair<std::string, size_t> p(field, i);	//create a pair associating the header name with the index
			fields.insert(p);							//insert the pair into the fields map
			i++;										//increment the data index
		}
		size_t nfields = fields.size();					//store the number of fields for each cell
		//load each cell and all associated fields
		std::string cell_line;							//string holds all information for a cell
		std::list<std::string> cell_list;				//list will be temporary storage for the cell fields
		while(std::getline(infile, cell_line)){			//for each cell entry
			cell_list.push_back(cell_line);				//push the cell entry into the list
		}
		//convert the list into actual data
		size_t ncells = cell_list.size();				//count the number of cells
		cells.resize(ncells);							//allocate enough space in the array to store all cells
		for(size_t c = 0; c < ncells; c++){				//for each cell entry in the list
			cells[c] = (double*) malloc(sizeof(double) * nfields);	//allocate enough space for each field
			std::stringstream fss(cell_list.front());	//turn the string representing the cell list into a stringstream
			for(size_t f = 0; f < nfields; f++){		//for each field
				fss>>cells[c][f];						//load the field
			}
			cell_list.pop_front();						//pop the read string off of the front of the list
		}
		infile.close();									//close the input file
		ip[0] = fields["x"];							//hard code the position indices for speed
		ip[1] = fields["y"];							//	this assumes all cells have positions
		ip[2] = fields["z"];
	}
	/// Return the value a specified field for a cell
	/// @param c is the cell index
	/// @param f is the field
	double value(size_t c, std::string f){
		size_t idx = fields[f];
		return cells[c][idx];
	}
	/// returns an ID used to look up a field
	bool exists(std::string f){
		std::unordered_map<std::string, size_t>::iterator iter = fields.find(f);
		if(iter == fields.end()) return false;
		else return true;
	}
	/// Return the position of cell [i]
	stim::vec3<double> p(size_t i){
		stim::vec3<double> pos(cells[i][ip[0]], cells[i][ip[1]], cells[i][ip[2]]);
		return pos;
	}
	/// Return the number of cells in the set
	size_t size(){
		return cells.size();
	}
	/// Return the maximum value of a field in this cell set
	double maximum(std::string field){
		size_t idx = fields[field];						//get the field index
		size_t ncells = cells.size();					//get the total number of cells
		double maxval, val;								//stores the current and maximum values
		for(size_t c = 0; c < ncells; c++){				//for each cell
			val = cells[c][idx];						//get the field value for this cell
			if(c == 0) maxval = val;					//if this is the first cell, just assign the maximum
			else if(val > maxval) maxval = val;			//	otherwise text for the size of val and assign it as appropriate
		}
		return maxval;
	}
	/// Return the maximum value of a field in this cell set
	double minimum(std::string field){
		size_t idx = fields[field];						//get the field index
		size_t ncells = cells.size();					//get the total number of cells
		double minval, val;								//stores the current and maximum values
		for(size_t c = 0; c < ncells; c++){				//for each cell
			val = cells[c][idx];						//get the field value for this cell
			if(c == 0) minval = val;					//if this is the first cell, just assign the maximum
			else if(val < minval) minval = val;			//	otherwise text for the size of val and assign it as appropriate
		}
		return minval;
	}
	/// adds a cell to the cell set
	void add(double x, double y, double z, double r = 0){
		if(cells.size() == 0){							//if the cell set is empty
			if(r == 0)									//if the radius is zero
				init_p3();								//initialize without a radius
			else
				init_p3r();								//otherwise initialize with a radius
		}
		size_t nf = fields.size();						//get the number of fields
		size_t bytes = sizeof(double) * nf;				//get the size of a cell field
		double* newcell = (double*) malloc(bytes);		//allocate memory for the new cell
		memset(newcell, 0, bytes);						//initialize all fields to zero
		newcell[ip[0]] = x;
		newcell[ip[1]] = y;
		newcell[ip[2]] = z;
		if(r != 0){										//if the user specifies a radius
			size_t ir = fields["radius"];				//get the index for the radius field
			newcell[ir] = r;							//add the radius to the field
		}
		cells.push_back(newcell);						//push the new memory entry into the cell array
	}
void save(std::string filename){
		
		size_t ncells = cells.size();	// get the number of cells
		std::ofstream file(filename);	//open a file to store the cell's coordinates
		if (file.is_open()) {
		
			file << "x y z radius\n";	//write the file header
			for (size_t c=0; c < ncells; c++){   //for each cell
				if (cells[c][ip[3]] != NULL)	//check if for the current cell, radius has been assigned
					file << cells[c][ip[0]]  << delim << cells[c][ip[1]] << delim << cells[c][ip[2]] << delim << cells[c][ip[3]] << '\n' ;
				else							//check if for the current cell, radius has not been assigned, set radius to zero
					file << cells[c][ip[0]]  << delim << cells[c][ip[1]] << delim << cells[c][ip[2]] << delim << 0 << '\n' ;
			}
			
		}
			file.close();
	}
	
};		//end class cellset
};		//end namespace stim
#endif