#include <CInterpolation.h>

Collaboration diagram for CInterpolation:

Public Member Functions
	CInterpolation ()

	CInterpolation (size_t stateNum, size_t reactionNum)

CStateRecord *	getInterpolationState ()

void	recordReset ()

void	recordState (const C_FLOAT64 time, const C_FLOAT64 *y)

	~CInterpolation ()

Private Member Functions
void	calculateFiringTime ()

void	calculateState ()

C_FLOAT64	interpolation (C_FLOAT64 x)

Private Attributes
size_t	mBeginRecordIndex

size_t	mHaveRecordNum

CStateRecord *	mpInterpolationState

CStateRecord *	mpState

C_FLOAT64 *	mpX

C_FLOAT64 *	mpY

bool	mShiftBeginIndex

size_t	mStateNum

Friends
void	printInterpolation (const CInterpolation &interp)

Detailed Description

Definition at line 102 of file CInterpolation.h.

Constructor & Destructor Documentation

CInterpolation::CInterpolation ( )

Default Constructor

Definition at line 95 of file CInterpolation.cpp.

                               :
   mpX(NULL),
   mpY(NULL),
   mStateNum(0),
   mBeginRecordIndex(0),
   mpState(NULL),
   mHaveRecordNum(0),
   mpInterpolationState(NULL),
   mShiftBeginIndex(false)
 {}

CInterpolation::CInterpolation	(	size_t	stateNum,
		size_t	reactionNum
	)

Assignment Constructor

Definition at line 106 of file CInterpolation.cpp.

References mpInterpolationState, and mpState.

                                                                 :
   mpX(new C_FLOAT64[stateNum]),
   mpY(new C_FLOAT64[stateNum]),
   mStateNum(stateNum),
   mBeginRecordIndex(0),
   mHaveRecordNum(0),
   mShiftBeginIndex(false)
 {
   mpInterpolationState = new CStateRecord(speciesNum);
 
   mpState = static_cast< CStateRecord* >
             (operator new[](sizeof(CStateRecord) * stateNum));
 
   for (size_t i = 0; i < stateNum; i++)
     new(&mpState[i]) CStateRecord(speciesNum);
 
   return;
 }

CInterpolation::~CInterpolation ( )

Destructor

Definition at line 125 of file CInterpolation.cpp.

References mpInterpolationState, mpState, mpX, mpY, mStateNum, pdelete, pdeletev, and CStateRecord::~CStateRecord().

 {
   pdeletev(mpX);
   pdeletev(mpY);
 
   pdelete(mpInterpolationState);
 
   for (int i = mStateNum - 1; i >= 0; i--)
     {
       mpState[i].~CStateRecord();
     }
 
   if (mpState != NULL)
     operator delete[](mpState);
 }

Member Function Documentation

void CInterpolation::calculateFiringTime ( )

private

Calculate time when slow reaction fires

Definition at line 232 of file CInterpolation.cpp.

References C_FLOAT64, CStateRecord::getArray(), CStateRecord::getSlowPropensitySum(), CStateRecord::getTime(), interpolation(), mBeginRecordIndex, mHaveRecordNum, mpInterpolationState, mpState, mpX, mpY, and mStateNum.

Referenced by getInterpolationState().

 {
   size_t i = mBeginRecordIndex;
 
   //(1)Load time into mpY
   //(2)Load sum of slow reaction propensities
   //   into mpX
 
   for (size_t j = 0; j < mHaveRecordNum; j++)
     {
       mpY[j] = mpState[i].getTime();
       mpX[j] = mpState[i].getSlowPropensitySum();
 
       if (++i > (mStateNum - 1))
         i -= mStateNum;
     }
 
   //(3)do interpolation at propensity=0
   C_FLOAT64 *pArray = mpInterpolationState->getArray();
   pArray[0] = interpolation(0.0);
 }

void CInterpolation::calculateState ( )

private

Calculate Model State at event time

Definition at line 254 of file CInterpolation.cpp.

References C_FLOAT64, CStateRecord::getArray(), CStateRecord::getArrayLen(), CStateRecord::getTime(), interpolation(), mBeginRecordIndex, mHaveRecordNum, mpInterpolationState, mpState, mpX, mpY, and mStateNum.

Referenced by getInterpolationState().

 {
   C_FLOAT64 *pInterArray = mpInterpolationState->getArray();
   C_FLOAT64 **ppArray = new C_FLOAT64*[mHaveRecordNum];
   C_FLOAT64 time = mpInterpolationState->getTime();
 
   //(1)Load time into array mpX
   //(2)Load pointers of objects CStateRecord into ppArray
   size_t i = mBeginRecordIndex;
 
   for (size_t j = 0; j < mHaveRecordNum; j++)
     {
       ppArray[j] = mpState[i].getArray();
       mpX[j] = mpState[i].getTime();
 
       if (++i > mStateNum - 1)
         i -= mStateNum;
     }
 
   //(3)Load propensities into mpY
   //(4)Do interpolation
   size_t speciesNum = (mpInterpolationState->getArrayLen()) - 2;
 
   for (size_t species = 1; species <= speciesNum; species++)
     {
       for (size_t j = 0; j < mHaveRecordNum; j++)
         mpY[j] = ppArray[j][species];
 
       pInterArray[species] = interpolation(time);
     }
 }

CStateRecord * CInterpolation::getInterpolationState ( )

Retrieve the model State at the slow reaction firing time by call calculateFiringTime() and calculateState()

Returns: CStateRecord mInterpolationState

Definition at line 286 of file CInterpolation.cpp.

References calculateFiringTime(), calculateState(), and mpInterpolationState.

Referenced by CHybridMethodODE45::doInverseInterpolation().

 {
   calculateFiringTime();
   calculateState();
   return mpInterpolationState;
 }

C_FLOAT64 CInterpolation::interpolation ( C_FLOAT64 x )

private

Calculate Interpolation Value

Parameters

C_FLOAT64 x

Definition at line 184 of file CInterpolation.cpp.

References C_FLOAT64, mHaveRecordNum, mpX, and mpY.

Referenced by calculateFiringTime(), and calculateState().

 {
   //Before use this function, data have been stored in
   //arraies mpX, mpY. Lengths of array is given by
   //mHaveRecordNum
 
   if ((x > mpX[mHaveRecordNum - 1])
       || (x < mpX[0])) //If x is within the range
     //Show Error !!!
     std::cout << "x " << x << " is out of range" << std::endl;
 
   if (mHaveRecordNum == 1)
     return mpY[0]; //only one point recorded
 
   size_t i, j;
 
   for (i = 0; i < mHaveRecordNum; i++) //x is in *mpX
     {
       if (x == mpX[i])
         return mpY[i];
     }
 
   //Calculate interpolation value
   C_FLOAT64 numer, tmpNumer = 1;
   C_FLOAT64 denom;
 
   for (i = 0; i < mHaveRecordNum; i++)
     tmpNumer *= x - mpX[i];
 
   C_FLOAT64 y = 0;
 
   for (i = 0; i < mHaveRecordNum; i++)
     {
       numer = tmpNumer / (x - mpX[i]);
       denom = 1;
 
       for (j = 0; j < mHaveRecordNum; j++)
         {
           if (i != j)
             denom *= mpX[i] - mpX[j];
         }
 
       y += mpY[i] * numer / denom;
     }
 
   return y;
 }

void CInterpolation::recordReset ( )

Reset the record

Definition at line 177 of file CInterpolation.cpp.

References mBeginRecordIndex, mHaveRecordNum, and mShiftBeginIndex.

Referenced by CHybridMethodODE45::doInverseInterpolation().

 {
   mBeginRecordIndex = 0;
   mHaveRecordNum    = 0;
   mShiftBeginIndex  = false;
 }

void CInterpolation::recordState	(	const C_FLOAT64	time,
		const C_FLOAT64 *	y
	)

Load the New State const C_FLOAT64 time, const unsigned C_INT32 len const C_FLOAT54 *y

Definition at line 141 of file CInterpolation.cpp.

References C_FLOAT64, CStateRecord::getArray(), CStateRecord::getArrayLen(), mBeginRecordIndex, mHaveRecordNum, mpState, mShiftBeginIndex, and mStateNum.

Referenced by CHybridMethodODE45::doInverseInterpolation().

 {
   mHaveRecordNum++;
 
   if (mHaveRecordNum > mStateNum)
     {
       mHaveRecordNum = mStateNum;
 
       if (!mShiftBeginIndex)
         mShiftBeginIndex = true;
     }
 
   if (mShiftBeginIndex)
     (++mBeginRecordIndex) %= mStateNum;
 
   size_t newStateIndex = mBeginRecordIndex
                          + mHaveRecordNum - 1;
 
   newStateIndex %= mStateNum;
   //std::cout << "newStateIndex: " << newStateIndex << std::endl;
 
   C_FLOAT64 *pArray = mpState[newStateIndex].getArray();
   const size_t size = mpState[newStateIndex].getArrayLen();
 
   pArray[0] = time;
 
   for (size_t i = 1; i < size; i++, y++)
     {
       ++pArray;
       *pArray = *y;
     }
 
   return;
 }

Friends And Related Function Documentation

void printInterpolation ( const CInterpolation & interp )

friend

Friend Function for Debug

Definition at line 294 of file CInterpolation.cpp.

 {
   std::cout << "State Num: " << interp.mStateNum << std::endl;
   std::cout << "Record Num: " << interp.mHaveRecordNum << std::endl;
   std::cout << "Begin Record Index: " << interp.mBeginRecordIndex << std::endl;
   std::cout << "Shift or not: ";
 
   if (interp.mShiftBeginIndex)
     std::cout << "Yes" << std::endl;
   else
     std::cout << "No" << std::endl;
 
   int site;
 
   for (size_t i = 0; i < interp.mHaveRecordNum; i++)
     {
       std::cout << "Record #:" << i << std::endl;
       site = (interp.mBeginRecordIndex + i) % interp.mStateNum;
       printRecord(interp.mpState[site]);
     }
 
   std::cout << "Finish, aha" << std::endl;
   return;
 }