RMOL Logo  0.25.2
C++ library of Revenue Management and Optimisation classes and functions
DPOptimiser.cpp
Go to the documentation of this file.
00001 // //////////////////////////////////////////////////////////////////////
00002 // Import section
00003 // //////////////////////////////////////////////////////////////////////
00004 // STL
00005 #include <cassert>
00006 #include <sstream>
00007 #include <vector>
00008 #include <cmath>
00009 // Boost Math
00010 #include <boost/math/distributions/normal.hpp>
00011 // StdAir
00012 #include <stdair/bom/LegCabin.hpp>
00013 #include <stdair/bom/VirtualClassStruct.hpp>
00014 #include <stdair/service/Logger.hpp>
00015 // RMOL
00016 #include <rmol/basic/BasConst_General.hpp>
00017 #include <rmol/bom/DPOptimiser.hpp>
00018 
00019 namespace RMOL {
00020   
00021   // ////////////////////////////////////////////////////////////////////
00022   void DPOptimiser::optimalOptimisationByDP (stdair::LegCabin& ioLegCabin) {
00023     // // Number of classes/buckets: n
00024     // const short nbOfClasses = ioBucketHolder.getSize();
00025 
00026     // // Number of values of x to compute for each Vj(x).
00027     // const int maxValue = static_cast<int> (iCabinCapacity * DEFAULT_PRECISION);
00028 
00029     // // Vector of the Expected Maximal Revenue (Vj).
00030     // std::vector< std::vector<double> > MERVectorHolder;
00031 
00032     // // Vector of V_0(x).
00033     // std::vector<double> initialMERVector (maxValue+1, 0.0);
00034     // MERVectorHolder.push_back (initialMERVector);
00035 
00036     // // Current cumulative protection level (y_j * DEFAULT_PRECISION).
00037     // // Initialise with y_0 = 0.
00038     // int currentProtection = 0;
00039 
00040     // int currentBucketIndex = 1;
00041     // ioBucketHolder.begin();
00042     
00043     // while (currentProtection < maxValue && currentBucketIndex < nbOfClasses) {
00044     // //while (currentBucketIndex == 1) {
00045     //   bool protectionChanged = false;
00046     //   double nextProtection = 0.0;
00047     //   std::vector<double> currentMERVector;
00048     //   // double testGradient = 10000;
00049       
00050     //   Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00051     //   const double meanDemand = currentBucket.getMean();
00052     //   const double SDDemand = currentBucket.getStandardDeviation();
00053     //   const double currentYield = currentBucket.getAverageYield();
00054     //   const double errorFactor = 1.0;
00055 
00056     //   Bucket& nextBucket = ioBucketHolder.getNextBucket();
00057     //   const double nextYield = nextBucket.getAverageYield();
00058       
00059     //   // For x <= currentProtection (y_(j-1)), V_j(x) = V_(j-1)(x).
00060     //   for (int x = 0; x <= currentProtection; ++x) {
00061     //     const double MERValue = MERVectorHolder.at(currentBucketIndex-1).at(x);
00062     //     currentMERVector.push_back (MERValue);
00063     //   }
00064 
00065     //   //
00066     //   boost::math::normal lNormalDistribution (meanDemand, SDDemand);
00067       
00068     //   // Vector of gaussian pdf values.
00069     //   std::vector<double> pdfVector;
00070     //   for (int s = 0; s <= maxValue - currentProtection; ++s) {
00071     //     const double pdfValue =
00072     //       boost::math::pdf (lNormalDistribution, s/DEFAULT_PRECISION);
00073     //     pdfVector.push_back (pdfValue);
00074     //   }
00075 
00076     //   // Vector of gaussian cdf values.
00077     //   std::vector<double> cdfVector;
00078     //   for (int s = 0; s <= maxValue - currentProtection; ++s) {
00079     //     const double cdfValue =
00080     //       boost::math::cdf (boost::math::complement (lNormalDistribution,
00081     //                                                  s/DEFAULT_PRECISION));
00082     //     cdfVector.push_back (cdfValue);
00083     //   }
00084       
00085     //   // Compute V_j(x) for x > currentProtection (y_(j-1)).
00086     //   for (int x = currentProtection + 1; x <= maxValue; ++x) {
00087     //     const double lowerBound = static_cast<double> (x - currentProtection);
00088         
00089     //     // Compute the first integral in the V_j(x) formulation (see
00090     //     // the memo of Jerome Contant).
00091     //     const double power1 =
00092     //       - 0.5 * meanDemand * meanDemand / (SDDemand * SDDemand);
00093     //     const double e1 = std::exp (power1);
00094     //     const double power2 = 
00095     //       - 0.5 * (lowerBound / DEFAULT_PRECISION - meanDemand)
00096     //       * (lowerBound / DEFAULT_PRECISION - meanDemand)
00097     //       / (SDDemand * SDDemand);
00098     //     const double e2 = std::exp (power2);
00099 
00100     //     const double cdfValue0 =
00101     //       boost::math::cdf (boost::math::complement (lNormalDistribution,
00102     //                                                  0.0));
00103     //     const double cdfValue1 =
00104     //       boost::math::cdf(boost::math::complement(lNormalDistribution,
00105     //                                                lowerBound/DEFAULT_PRECISION));
00106     //     const double integralResult1 = currentYield
00107     //       * ((e1 - e2) * SDDemand / sqrt (2 * 3.14159265)
00108     //          + meanDemand * (cdfValue0 - cdfValue1));
00109         
00110     //     double integralResult2 = 0.0;
00111         
00112     //     for (int s = 0; s < lowerBound; ++s) {
00113     //       const double partialResult =
00114     //         2 * MERVectorHolder.at(currentBucketIndex-1).at(x-s)
00115     //         * pdfVector.at(s);
00116           
00117     //       integralResult2 += partialResult;
00118     //     }
00119     //     integralResult2 -= MERVectorHolder.at(currentBucketIndex-1).at(x) *
00120     //       pdfVector.at(0);
00121         
00122     //     const int intLowerBound = static_cast<int>(lowerBound);
00123     //     integralResult2 += 
00124     //       MERVectorHolder.at(currentBucketIndex-1).at(x - intLowerBound) *
00125     //       pdfVector.at(intLowerBound);
00126         
00127     //     integralResult2 /= 2 * DEFAULT_PRECISION;
00128     //     /*
00129     //     for (int s = 0; s < lowerBound; ++s) {
00130     //       const double partialResult =
00131     //         (MERVectorHolder.at(currentBucketIndex-1).at(x-s) +
00132     //          MERVectorHolder.at(currentBucketIndex-1).at(x-s-1)) *
00133     //         (cdfVector.at(s+1) - cdfVector.at(s)) / 2;
00134     //       integralResult2 += partialResult;
00135     //     }
00136     //     */
00137     //     const double firstElement = integralResult1 + integralResult2;
00138         
00139     //     // Compute the second integral in the V_j(x) formulation (see
00140     //     // the memo of Jerome Contant).
00141     //     const double constCoefOfSecondElement =
00142     //       currentYield * lowerBound / DEFAULT_PRECISION
00143     //       + MERVectorHolder.at(currentBucketIndex-1).at(currentProtection);
00144 
00145     //     const double secondElement = constCoefOfSecondElement
00146     //       * boost::math::cdf(boost::math::complement(lNormalDistribution,
00147     //                                                  lowerBound/DEFAULT_PRECISION));
00148 
00149     //     const double MERValue = (firstElement + secondElement) / errorFactor;
00150 
00151         
00152     //     assert (currentMERVector.size() > 0);
00153     //     const double lastMERValue = currentMERVector.back();
00154 
00155     //     const double currentGradient = 
00156     //       (MERValue - lastMERValue) * DEFAULT_PRECISION;
00157 
00158     //     //assert (currentGradient >= 0);
00159     //     if (currentGradient < -0) {
00160     //       std::ostringstream ostr;
00161     //       ostr << currentGradient  << std::endl
00162     //            << "x = " << x << std::endl
00163     //            << "class: " << currentBucketIndex << std::endl;
00164     //       STDAIR_LOG_DEBUG (ostr.str());
00165     //     }
00166           
00167     //     /*
00168     //     assert (currentGradient <= testGradient);
00169     //     testGradient = currentGradient;
00170     //     */
00171     //     if (protectionChanged == false && currentGradient <= nextYield) {
00172     //       nextProtection = x - 1;
00173     //       protectionChanged = true;
00174     //     }
00175 
00176     //      if (protectionChanged == true && currentGradient > nextYield) {
00177     //       protectionChanged = false;
00178     //     }
00179         
00180     //     if (protectionChanged == false && x == maxValue) {
00181     //       nextProtection = maxValue;
00182     //     }
00183         
00184     //     currentMERVector.push_back (MERValue); 
00185     //   }
00186 
00187     //   // DEBUG
00188     //   STDAIR_LOG_DEBUG ("Vmaxindex = " << currentMERVector.back());
00189         
00190     //   MERVectorHolder.push_back (currentMERVector);
00191      
00192     //   const double realProtection = nextProtection / DEFAULT_PRECISION;
00193     //   const double bookingLimit = iCabinCapacity - realProtection;
00194 
00195     //   currentBucket.setCumulatedProtection (realProtection);
00196     //   nextBucket.setCumulatedBookingLimit (bookingLimit);
00197 
00198     //   currentProtection = static_cast<int> (std::floor (nextProtection));
00199       
00200     //   ioBucketHolder.iterate();
00201     //   ++currentBucketIndex;
00202     // }
00203   }
00204 
00205 }