---

ChemicalTimeSeriesData.cpp

Go to the documentation of this file.

// ChemicalTimeSeriesData.cpp: implementation of the ChemicalTimeSeriesData class.
//

#include "ChemicalTimeSeriesData.h"

// Construction/Destruction

// the single-argument constructor + AttachNewTimeSeries combination can
// eventually replace the two-argument constructor
// KSB 10/11/00

ChemicalTimeSeriesData::ChemicalTimeSeriesData(ReactionNetwork *reactionNetwork)
{
        this->nChem = reactionNetwork->GetNumberOfChemicals();
        chemData.resize(0);
        for(int i = 0; i < nChem; i++)
        {
                std::vector<DataPoint *> localVector;
                chemData.push_back(localVector);
        }

}

// Chemical time series data is read in from a file, with the format
//              time            Chemical Name           value           error
// You must put something for the "error" column even if there are 
// no errors in your data (i.e. if it is fake test data); errors of
// 1.0 are fine.  Data does not have to be ordered in any way; entries
// in the above format should be input with no separating whitespace
//              The constructor has to take a "ReactionNetwork" argument 
// because otherwise it can't look up where the chemicals are

ChemicalTimeSeriesData::ChemicalTimeSeriesData(ReactionNetwork *reactionNetwork, std::string fileName)
{
        this->nChem = reactionNetwork->GetNumberOfChemicals();

        // used for reading in the data
        int whereAmI;
        int i;
        char lineBuf[1000];
        double time;
        char name[100];
        double value;
        double error;

        this->nChem = reactionNetwork->GetNumberOfChemicals();
        chemData.resize(0);
        for(i = 0; i < nChem; i++)
        {
                std::vector<DataPoint *> localVector;
                chemData.push_back(localVector);
        }

        // open the file for reading
        ifstream dataFile(fileName.c_str());
        assert(dataFile);
        
        // put the data line by line into the linebuffer and then parse it
        // for right now, this is really slow 
        // We could (and should) speed this up by allowing
        // a search of the chemicals vector to find the chemical number based
        // on the name
        while(dataFile.getline(lineBuf,sizeof(lineBuf)))
        {
                if(lineBuf[0] == '%') {continue;}
                std::istringstream dataStream(lineBuf);
                dataStream >> time >> name >> value >> error;
                
                // find the number the network has assigned to this chemical
                std::string nameString(name);
                for(whereAmI = 0; whereAmI < nChem; whereAmI++)
                {       
                        bool isEqual = ( nameString == (reactionNetwork->GetChemical(whereAmI)->GetName()) );
                        if(isEqual) {break;}    
                }
                // create a data point and put it in the right place
                DataPoint *dp = new DataPoint(time,value,error);
                chemData[whereAmI].push_back(dp);
        }
        // sort the chemData
        SortChemTimeSeriesData();
        // create the time vector and synchronize it to the data
        TimeVectorSynchronize();
        
        // set the number of time series we have data for
        int counter = 0;
        for(i = 0; i < nChem; i++)
        {
                if(chemData[i].size() > 0) {counter++;}
        }
        nTimeSeries = counter;
        
        return;
}

// This reads a new time series from fileName and adds it to the cTSD.
// As of now, there is no safety mechanism to ensure you don't contra-
// dict existing data!  You should be careful that each new file adds
// data for chemicals for which you have no previous data (for now).
// - K.S.B.

void ChemicalTimeSeriesData::AttachNewTimeSeries(ReactionNetwork *reactionNetwork, std::string fileName)
{
        this->nChem = reactionNetwork->GetNumberOfChemicals();

        // used for reading in the data
        int whereAmI;
        char lineBuf[1000];
        double time;
        char name[100];
        double value;
        double error;


        // open the file for reading
        ifstream dataFile(fileName.c_str());
        assert(dataFile);
        
        // put the data line by line into the linebuffer and then parse it
        while(dataFile.getline(lineBuf,sizeof(lineBuf)))
        {
                if(lineBuf[0] == '%') {continue;}
                std::istringstream dataStream(lineBuf);
                dataStream >> time >> name >> value >> error;
                // find the number the network has assigned to this chemical
                std::string nameString(name);
                for(whereAmI = 0; whereAmI < nChem; whereAmI++)
                {       
                        bool isEqual = ( nameString == (reactionNetwork->GetChemical(whereAmI)->GetName()) );
                        if(isEqual) {break;}    
                }
                // create a data point and put it in the right place
                //DataPoint *dp = new DataPoint(time,value,error);
                chemData[whereAmI].push_back(new DataPoint(time,value,error));
        }
        // sort the chemData
        SortChemTimeSeriesData();
        // recreate and synchronize the time vector
        TimeVectorSynchronize();
        
        // set the number of time series 
        // you can't just increment nTimeSeries because the new files
        // you add can have any number of new time series in them
        int counter = 0;
        for(int i = 0; i < nChem; i++)
        {
                //if( chemData[i] != (std::vector<DataPoint *>)NULL ) {counter++;}
                //if( chemData[i] != NULL ) {counter++;}
                if(chemData[i].size() > 0) {counter++;}
        }
        nTimeSeries = counter;

        
        return;
}

void ChemicalTimeSeriesData::TimeVectorSynchronize()
{
        // erase the old time vector
        timeVector.erase(timeVector.begin(),timeVector.end());
        // run through the data and create the time vector
        for(int i = 0; i < nChem; i++)
        {
                if(chemData[i].size() > 0)
                {
                        for(int pos = 0; pos < chemData[i].size(); pos++)
                        {
                                double time = (chemData[i])[pos]->GetTime();
                                TimeVectorPoint *tp = new TimeVectorPoint(time,pos,i);
                                timeVector.push_back(tp);
                        }
                }
        }
        // sort the newly created time vector
        SortTimeVector();
}


ChemicalTimeSeriesData::~ChemicalTimeSeriesData()
{
        int i,j;
        // cleanup
        for(i = 0; i < nChem; i++)
        {
                for(j = 0; j < chemData[i].size(); j++)
                {
                        delete chemData[i][j];
                }
        }
        for(i = 0; i < timeVector.size(); i++)
        {
                delete timeVector[i];
        }
}

std::vector<DataPoint *> ChemicalTimeSeriesData::GetTimeSeries(int chemIndex)
{
        return chemData[chemIndex];
}

void ChemicalTimeSeriesData::SortChemTimeSeriesData()
{
        // use STL sort with a comparator that accepts pointers
        DataPoint::DataPointComparator comparator;
        for(int i = 0; i < nChem; i++)
        {
                if( chemData[i].size() > 0 )
                {
                        std::sort(chemData[i].begin(),chemData[i].end(),comparator);
                }
        }
}

void ChemicalTimeSeriesData::SortTimeVector()
{
        // use STL sort with the correct comparator
        TimeVectorPoint::TimeVectorPointComparator comparator;
        std::sort(timeVector.begin(),timeVector.end(),comparator);
}

int ChemicalTimeSeriesData::GetNDataPoints()
{
        return timeVector.size();
}

int ChemicalTimeSeriesData::GetNTimeSeries()
{
        return nTimeSeries;
}

int ChemicalTimeSeriesData::GetNChemicals()
{
        return nChem;
}

double ChemicalTimeSeriesData::GetNextDataTime(double time)
{
        // STL algorithm
        // define an iterator for the timeVector
        std::vector<TimeVectorPoint *>::iterator it;
        // the comparator is defined in terms of TimeVectorPoint *, so
        // that's the value we need (it's a dummy in terms of everything
        // except time
        TimeVectorPoint *tp = new TimeVectorPoint(time,0,0);
        // comparator
        TimeVectorPoint::TimeVectorPointComparator comparator;
        // get the iterator
        it = std::lower_bound(timeVector.begin(),timeVector.end(),tp,comparator); 
        // iterator points to a data point; if it's at the end the comparison
        // time is off-the-end of the vector
        if( it != timeVector.end() )
        {
                return (*it)->GetTime();
        }
        else
        {
                // throwing an exception is just temporary
                throw std::runtime_error("Trial value is past final data point");
        }
        
        // cleanup
        delete tp;
//      delete it;
}

// The function GetLatestDataTime() sorts the timeVector and returns
// the latest time at which any chemical data is present.  This is of
// use in deciding how long the ReactionMover moves the network. 
// - K.S.B.

double ChemicalTimeSeriesData::GetLatestDataTime()
{
        SortTimeVector();
        double latestTime = (timeVector[timeVector.size() - 1])->GetTime();
        return latestTime;
}

---