---

ReactionNetwork.cpp

Go to the documentation of this file.

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

#ifdef _WINDOWS
#include "stdafx.h"
#endif

#include "ReactionNetwork.h"
#include "Reaction.h"

// Construction/Destruction

ReactionNetwork::ReactionNetwork()
{
        // sort the list of algebraic indices - this way we don't 
        // have to jump around in the chemical vector
        std::sort(algebraicIndices.begin(),algebraicIndices.end());
}

ReactionNetwork::~ReactionNetwork()
{
        int i;
        // first delete chemicals and rate constants
        for (i = 0; i < chemicals.size(); i++) delete chemicals[i];
        for (i = 0; i < rateConstants.size(); i++) delete rateConstants[i];
        // now delete reactions
        for (i = 0; i < reactions.size(); i++) delete reactions[i];
        // clean up vector of algebraic indices
        algebraicIndices.erase(algebraicIndices.begin(),algebraicIndices.end());
}

int ReactionNetwork::GetNumberOfChemicals() const
{
        return chemicals.size();
}       

int ReactionNetwork::GetNumberOfAlgebraicChemicals() const
{
        return algebraicIndices.size();
}

int ReactionNetwork::GetNumberOfReactions() const
{
        return reactions.size();
}

int ReactionNetwork::GetNumberOfRateConstants() const
{
        return rateConstants.size();
}       

void ReactionNetwork::ChemicalReset()
{
        for(int i = 0; i < GetNumberOfChemicals(); i++)
        {
                GetChemical(i)->Reset();
        }
}

void ReactionNetwork::RateConstantReset()
{
        for(int i = 0; i < GetNumberOfRateConstants(); i++)
        {
                GetRateConstant(i)->Reset();
        }
}

Chemical *ReactionNetwork::GetChemical(int chemicalNumber)
{
        return chemicals[chemicalNumber];
}       

RateConstant *ReactionNetwork::GetRateConstant(int rateConstantNumber)
{
        return rateConstants[rateConstantNumber];
}       

Reaction *ReactionNetwork::GetReaction(int reactionNumber)
{
        return reactions[reactionNumber];
}

std::vector<int> *ReactionNetwork::GetAlgebraicIndices()
{
        return &algebraicIndices;
}

std::vector<Reaction*> *ReactionNetwork::GetReactions()
{
        return &reactions;
}

std::vector<Chemical*> *ReactionNetwork::GetChemicals()
{
        return &chemicals;
}

std::vector<RateConstant*> *ReactionNetwork::GetRateConstants()
{
        return &rateConstants;
}

std::vector<double> *ReactionNetwork::GetReactionRates()
{
        // ensure synchrony with the reactions
        if(reactionRates.size() != this->GetNumberOfReactions())
        {
                // clear out the old vector
                reactionRates.erase(reactionRates.begin(),reactionRates.end());
                // fill the vector with the current rates
                for(int i = 0; i < this->GetNumberOfReactions(); i++)
                {
                        reactionRates.push_back(reactions[i]->GetRate());
                }
        }
        else
        {
                for(int i = 0; i < this->GetNumberOfReactions(); i++)
                {
                        reactionRates[i] = reactions[i]->GetRate();
                }
        }
        
        // return the vector
        return &reactionRates;
}

void ReactionNetwork::SetAmount(int chemicalNumber, double amount)
{
        this->GetChemical(chemicalNumber)->SetAmount(amount);
}

void ReactionNetwork::AddAlgebraicIndex(Chemical *chemical)
{
        int index = chemical->GetChemicalNumber();
        algebraicIndices.push_back(index);
}

void ReactionNetwork::RemoveReaction(int rxnNumber)
{
        // first delete the reaction so memory is properly freed
        delete reactions[rxnNumber];
        // now remove from the vector
        reactions.erase(reactions.begin() + rxnNumber);
}

void ReactionNetwork::DumpTeXForm(std::string fileName)
{
    int i;
        int nChem = chemicals.size();
        int nReac = reactions.size();
        int cNum = 0;
        std::string sign;
        
        std::string junk;
        junk = "-";

        int nTerms = this->GetNTermsPerLine();
        
        // keeps track of how many terms have been written so linebreaks can be put in
        std::vector<int> viTerms;
#ifdef _WIN32
//      std::vector<strstream *> vpStreams;
        std::vector<std::ostringstream *> vpStreams;
#else
        std::vector<std::ostringstream *> vpStreams;
#endif  // _WIN32

        for(i = 0; i < nChem; i++)
        {
                viTerms.push_back(0);
#ifdef _WIN32
//              vpStreams.push_back(new strstream);
                vpStreams.push_back(new std::ostringstream);
#else
                vpStreams.push_back(new std::ostringstream);
#endif
        }

        fstream *pTeXFile = new fstream(fileName.c_str(),ios::out);

        *pTeXFile << "\\documentclass[12pt]{article}" << endl;
        *pTeXFile << "\\topmargin 0.0cm" << endl;
        *pTeXFile << "\\oddsidemargin 0.2cm" << endl;
        *pTeXFile << "\\textwidth 16cm" << endl;
        *pTeXFile << "\\textheight 21cm" << endl;
        *pTeXFile << "\\footskip 1.0cm" << endl;
        *pTeXFile << "\\begin{document}" << endl << endl;
        *pTeXFile << "\\begin{eqnarray*}" << endl;
        
        
        // build up the text RHSs the same way that the RHSs are numerically computed
        for (int rNum = 0; rNum < nReac; rNum++)
        {
                std::vector<Chemical *> *localChemicals = reactions[rNum]->GetChemicals();
                int nLC = localChemicals->size();
                
                for (int lcNum = 0; lcNum < nLC; lcNum++)
                {
                        cNum = localChemicals->at(lcNum)->GetChemicalNumber();
                        int nChanged = reactions[rNum]->GetNumberChangedByReaction(lcNum);
                        // next lines are for formatting
                        if(nChanged == 0)
                        {
                                // don't print anything
                        }
                        else 
                        {
                                if(nChanged > 0){sign = "+";}
                                else {sign = "-";}
                                
                                if(abs(nChanged) == 1)
                                {
                                        if((viTerms[cNum] % nTerms == 0) && viTerms[cNum] > 0)
                                        {
                                                *vpStreams[cNum] << "\\\\ \n & &";
                                        }
                                        // don't print the 1, just the sign 
                                        *vpStreams[cNum] << sign.c_str() << " " << reactions[rNum]->GetTeXForm().c_str() << " ";
                                }
                                else
                                {
                                        if((viTerms[cNum] % nTerms == 0) && viTerms[cNum] > 0)
                                        {
                                                *vpStreams[cNum] << "\\\\ \n & &";
                                        }
                                        // print the sign and then abs(nChanged)
                                        *vpStreams[cNum] << sign.c_str() << " " << abs(nChanged) << reactions[rNum]->GetTeXForm().c_str() << " ";
                                }
                                viTerms[cNum]++;
                        }
                }
        }
        

        
        // now loop over chemicals and write everything
        int counter = 0;
        for(int cLoop = 0;  cLoop < nChem; cLoop++)
        {       
                // if this is an algebraic, ignore it
                if((counter < algebraicIndices.size()) && (cLoop == algebraicIndices[counter]))
                {
                        counter++;
                        continue;
                }
                
                *pTeXFile << "{d \\left[ " << chemicals[cLoop]->GetTeXName().c_str() << "\\right]";
                *pTeXFile << " \\over {dt}} & = & ";
                
                if(viTerms[cLoop] == 0)
                {
                        *pTeXFile << "0 \\\\" << endl;
                }
//              *vpStreams[cLoop] << std::ends;
                
                *pTeXFile << vpStreams[cLoop]->str(); 
                *pTeXFile << "\\\\" << endl;
        }
        
        
        *pTeXFile << "\\end{eqnarray*}" << endl << endl;
        *pTeXFile << "\\end{document}" << endl;
        
        for(i = 0; i < vpStreams.size(); i++)
        {
//              vpStreams[i]->rdbuf()->freeze(0);
                delete vpStreams[i];
        }
        delete pTeXFile;
}

void ReactionNetwork::DumpPetriNetDotFile(std::string fileName)
{
        int i;
        int nChem = chemicals.size();
        int nReac = reactions.size();
        int cNum = 0;
        int lcNum = 0;
        
        // open file for writing
        fstream *pDotFile = new fstream(fileName.c_str(),ios::out);

        *pDotFile << "digraph G {" << endl;
        if(this->GetLRGraph())
        {
                *pDotFile << "\t rankdir = LR;" << endl;
        }
        *pDotFile << "\t center = true;" << endl;
        
        // first list all the chemicals (circles) and reactions (boxes)
        *pDotFile << "\t node [shape = circle];" << endl;
        for(i = 0; i < nChem; i++)
        {
                *pDotFile << "\t\t " << "\"" << chemicals[i]->GetName().c_str() << "\";" << endl;
        }
        *pDotFile << "\t node [shape = box, height = 0.75, width = 0.25];" << endl;
        for(i = 0; i < nReac; i++)
        {
                *pDotFile << "\t\t " << "\"" << reactions[i]->GetName().c_str() << "\";" <<  endl;
        }
        
        // now connect up the circles and boxes - stoichiometry gives the direction of
        // the arrows
        for (int rNum = 0; rNum < nReac; rNum++)
        {
                std::vector<Chemical *> *localChemicals = reactions[rNum]->GetChemicals();
                int nLC = localChemicals->size();

                // loop over the chemicals involved in this reaction
                for(lcNum = 0; lcNum < nLC; lcNum++)
                {
                        cNum = localChemicals->at(lcNum)->GetChemicalNumber();
                        int nChanged = reactions[rNum]->GetNumberChangedByReaction(lcNum);
                        
                        if(nChanged == 0) // catalytic, arrow goes to and from reaction
                        {
                                *pDotFile << "\t\t " << "\"" << chemicals[cNum]->GetName().c_str() << "\"" << " -> ";
                                *pDotFile << "\"" << reactions[rNum]->GetName().c_str() << "\"" << ";" << endl;
                                *pDotFile << "\t\t " << "\"" << reactions[rNum]->GetName().c_str() << "\"" << " -> ";
                                *pDotFile << "\"" << chemicals[cNum]->GetName().c_str() << "\"" << ";" << endl;
                        }
                        else if(nChanged > 0) // from transition to place
                        {
                                *pDotFile << "\t\t " << "\"" << reactions[rNum]->GetName().c_str() << "\"" << " -> ";
                                *pDotFile << "\"" << chemicals[cNum]->GetName().c_str() << "\"" << ";" << endl;
                        }
                        else // from place to transition
                        {
                                *pDotFile << "\t\t " << "\"" << chemicals[cNum]->GetName().c_str() << "\"" << " -> ";
                                *pDotFile << "\"" << reactions[rNum]->GetName().c_str() << "\"" << ";" << endl;
                        }
                }
                
        }

        // write closing bracket
        *pDotFile << "}" << endl;

        delete pDotFile;
}

---