src/core/biolccc.cpp
00001 #include <iostream>
00002 #include <vector>
00003 #include <cmath>
00004 #include <algorithm>
00005 #include "biolccc.h"
00006 
00007 namespace BioLCCC
00008 {
00009 
00010 // Auxiliary functions that shouldn't be exposed to a user.
00011 namespace
00012 {
00013 double calculateKd(const std::vector<ChemicalGroup> &parsedSequence,
00014                    const double secondSolventConcentration,
00015                    const ChemicalBasis &chemBasis,
00016                    const double columnPoreSize,
00017                    const double columnRelativeStrength,
00018                    const double temperature) throw(BioLCCCException)
00019 {
00020     // assymetricCalculations shows whether the Kd for the reversed molecule 
00021     // will differ. It happens when a molecule cannot be divided into an integer
00022     // number of Kuhn segments.
00023     bool assymetricCalculations = 
00024         (fmod(chemBasis.monomerLength() * parsedSequence.size(), 
00025               chemBasis.kuhnLength()) != 0);
00026     // Choosing the appropriate polymerModel.
00027     if (chemBasis.polymerModel()==CHAIN)
00028     {
00029         double Kd = calculateKdChain(parsedSequence,
00030                                     secondSolventConcentration, 
00031                                     chemBasis, columnPoreSize,
00032                                     columnRelativeStrength, temperature);
00033         if (assymetricCalculations) 
00034         {
00035             std::vector<ChemicalGroup> revParsedSequence = 
00036                 parsedSequence;
00037             std::reverse(revParsedSequence.begin(),
00038                          revParsedSequence.end());
00039             Kd = (Kd + calculateKdChain(revParsedSequence,
00040                                        secondSolventConcentration,
00041                                        chemBasis,
00042                                        columnPoreSize,
00043                                        columnRelativeStrength, 
00044                                        temperature)) / 2.0 ;
00045         }
00046         return Kd;
00047     }
00048     else if (chemBasis.polymerModel() == ROD)
00049     {
00050         double Kd = calculateKdRod(parsedSequence,
00051                                    secondSolventConcentration, 
00052                                    chemBasis, columnPoreSize,
00053                                    columnRelativeStrength, temperature);
00054         if (assymetricCalculations) 
00055         {
00056             std::vector<ChemicalGroup> revParsedSequence = 
00057                 parsedSequence;
00058             std::reverse(revParsedSequence.begin(),
00059                          revParsedSequence.end());
00060             Kd = (Kd + calculateKdRod(revParsedSequence,
00061                                       secondSolventConcentration,
00062                                       chemBasis,
00063                                       columnPoreSize,
00064                                       columnRelativeStrength, 
00065                                       temperature)) / 2.0 ;
00066         }
00067         return Kd;
00068     }
00069     else
00070     {
00071         throw BioLCCCException("Model error.");
00072     }
00073 }
00074 
00075 class KdCalculator
00076 {
00077 public:
00078     KdCalculator(const std::vector<ChemicalGroup> &parsedSequence,
00079                  const ChemicalBasis &chemBasis, 
00080                  const double columnPoreSize,
00081                  const double columnRelativeStrength,
00082                  const double temperature,
00083                  const int numInterpolationPoints) :
00084                      mParsedSequence(parsedSequence),
00085                      mChemicalBasis(chemBasis),
00086                      mColumnPoreSize(columnPoreSize),
00087                      mColumnRelativeStrength(columnRelativeStrength),
00088                      mTemperature(temperature),
00089                      mN(numInterpolationPoints)
00090     {
00091         if (mN > 0)
00092         {
00093             mSecondSolventConcentrations = new double[mN];
00094             mLogKds = new double[mN];
00095             // The number of extra points in the terminal segments.
00096             // This points significantly increase the accuracy of spline
00097             // interpolation.
00098             int NETP = 1;
00099             for (int i=0; i < mN; i++) 
00100             {
00101                 if (i <= NETP)
00102                 {
00103                     mSecondSolventConcentrations[i] =
00104                         i * 100.0 / (mN - 2.0 * NETP - 1.0) / (NETP + 1.0);
00105                 }
00106                 else if (i > (mN - NETP - 2))
00107                 {
00108                     mSecondSolventConcentrations[i] = 
00109                         ((mN - 2.0 * NETP - 2.0)
00110                             + (i - mN + NETP + 2.0) / (NETP + 1.0))
00111                          * 100.0 / (mN - 2.0 * NETP - 1.0);
00112                 }
00113                 else
00114                 {
00115                     mSecondSolventConcentrations[i] =
00116                         (i - NETP) * 100.0 / (mN - 2.0 * NETP - 1.0);
00117                 }
00118                 mLogKds[i] = log(calculateKd(mParsedSequence,
00119                     mSecondSolventConcentrations[i],
00120                     mChemicalBasis, 
00121                     mColumnPoreSize,
00122                     mColumnRelativeStrength,
00123                     mTemperature));
00124             }
00125 
00126             mSecondDers = new double[mN];
00127             fitSpline(mSecondSolventConcentrations, mLogKds,
00128                 mN, mSecondDers);
00129         }
00130     }
00131 
00132     ~KdCalculator()
00133     {
00134         if (mN > 0)
00135         {
00136             delete[] mSecondSolventConcentrations;
00137             delete[] mLogKds;
00138             delete[] mSecondDers;
00139         }
00140     }
00141 
00142     double operator()(double secondSolventConcentration) 
00143         throw (BioLCCCException)
00144     {
00145         if (mN == 0) 
00146         {
00147             return calculateKd(mParsedSequence,
00148                                secondSolventConcentration,
00149                                mChemicalBasis, 
00150                                mColumnPoreSize,
00151                                mColumnRelativeStrength,
00152                                mTemperature);
00153         }
00154         else 
00155         {
00156             return exp(calculateSpline(mSecondSolventConcentrations, mLogKds,
00157                 mSecondDers, mN, 
00158                 secondSolventConcentration));
00159         }
00160     }
00161 
00162 private:
00163     const std::vector<ChemicalGroup> &mParsedSequence;
00164     const ChemicalBasis & mChemicalBasis;
00165     const double mColumnPoreSize;
00166     const double mColumnRelativeStrength;
00167     const double mTemperature;
00168     const int mN;
00169     double * mSecondSolventConcentrations;
00170     double * mLogKds;
00171     double * mSecondDers;
00172 };
00173 
00174 double calculateRT(const std::vector<ChemicalGroup> &parsedSequence,
00175                    const ChemicalBasis &chemBasis,
00176                    const ChromoConditions &conditions,
00177                    const int numInterpolationPoints,
00178                    const bool continueGradient,
00179                    const bool backwardCompatibility
00180                    ) throw(BioLCCCException)
00181 {
00182     if (numInterpolationPoints < 0)
00183     {
00184         throw BioLCCCException(
00185             "The number of interpolation points must be non-negative.");
00186     }
00187 
00188     KdCalculator kdCalculator(parsedSequence, chemBasis,
00189                               conditions.columnPoreSize(),
00190                               conditions.columnRelativeStrength(),
00191                               conditions.temperature(),
00192                               numInterpolationPoints);
00193 
00194     double RT = 0.0;
00195     // Use simplified expression for isocratic elution.
00196     if (conditions.SSConcentrations().size() == 1)
00197     {
00198         RT = kdCalculator(conditions.SSConcentrations()[0]) 
00199              * conditions.columnPoreVolume() / conditions.flowRate();
00200     }
00201     else
00202     {
00203         // The part of a column passed by molecules. When it exceeds 1.0,
00204         // the analyte elutes from the column.
00205         double S = 0.0;
00206         double dS = 0.0;
00207         int j = 0;
00208         double currentSSConcentration = 0.0;
00209         while (S < 1.0)
00210         {
00211             j++;
00212             if (j < conditions.SSConcentrations().size())
00213             {
00214                 currentSSConcentration = conditions.SSConcentrations()[j];
00215             }
00216             else
00217             {
00218                 // If continue gradient then use the slope of the last section.
00219                 if (continueGradient)
00220                 {
00221                     currentSSConcentration += 
00222                         conditions.SSConcentrations().back()
00223                         - *(conditions.SSConcentrations().end() - 2);
00224                 }
00225                 else
00226                 {
00227                     break;
00228                 }
00229             }
00230             dS = conditions.dV() 
00231                  / kdCalculator(currentSSConcentration) 
00232                  / conditions.columnPoreVolume();
00233             S += dS;
00234         }
00235 
00236         RT = j * conditions.dV() / conditions.flowRate();
00237         // Correction for the discreteness of integration.
00238         if ((!backwardCompatibility) && (S > 1.0))
00239         {
00240             RT -= (S - 1.0) / dS * conditions.dV() / conditions.flowRate();
00241         }
00242     }
00243 
00244     RT += conditions.delayTime();
00245     // Correction for V0.
00246     RT += conditions.columnInterstitialVolume() / conditions.flowRate();
00247     return RT;
00248 }
00249 }
00250 
00251 double calculateRT(const std::string &sequence,
00252                    const ChemicalBasis &chemBasis,
00253                    const ChromoConditions &conditions,
00254                    const int numInterpolationPoints,
00255                    const bool continueGradient,
00256                    const bool backwardCompatibility) 
00257                    throw(BioLCCCException)
00258 {
00259     std::vector<ChemicalGroup> parsedSequence = 
00260         parseSequence(sequence, chemBasis);
00261     return calculateRT(parsedSequence,
00262                        chemBasis,
00263                        conditions,
00264                        numInterpolationPoints,
00265                        continueGradient,
00266                        backwardCompatibility);
00267 }
00268 
00269 double calculateKd(const std::string &sequence,
00270                    const double secondSolventConcentration,
00271                    const ChemicalBasis & chemBasis,
00272                    const double columnPoreSize,
00273                    const double columnRelativeStrength,
00274                    const double temperature)
00275                    throw(BioLCCCException)
00276 {
00277     return calculateKd(parseSequence(sequence, chemBasis),
00278                        secondSolventConcentration,
00279                        chemBasis,
00280                        columnPoreSize,
00281                        columnRelativeStrength,
00282                        temperature);
00283 }
00284 
00285 double calculateAverageMass(const std::string &sequence,
00286                             const ChemicalBasis &chemBasis)
00287                             throw(BioLCCCException)
00288 {
00289     std::vector<ChemicalGroup> parsedSequence =
00290         parseSequence(sequence, chemBasis);
00291     double peptideAverageMass = 0;
00292     for (std::vector<ChemicalGroup>::const_iterator i =
00293                 parsedSequence.begin();
00294             i < parsedSequence.end();
00295             i++)
00296     {
00297         peptideAverageMass += i -> averageMass();
00298     }
00299 
00300     return peptideAverageMass;
00301 }
00302 
00303 double calculateMonoisotopicMass(const std::string &sequence,
00304                                  const ChemicalBasis &chemBasis)
00305                                  throw(BioLCCCException)
00306 {
00307     std::vector<ChemicalGroup> parsedSequence =
00308         parseSequence(sequence, chemBasis);
00309     double peptideMonoisotopicMass = 0;
00310     for (std::vector<ChemicalGroup>::const_iterator i =
00311                 parsedSequence.begin();
00312             i < parsedSequence.end();
00313             i++)
00314     {
00315         peptideMonoisotopicMass += i -> monoisotopicMass();
00316     }
00317 
00318     return peptideMonoisotopicMass;
00319 }
00320 
00321 }
 All Classes Namespaces Functions Variables Enumerations Enumerator