RMOL Logo  0.25.3
C++ library of Revenue Management and Optimisation classes and functions
MCOptimiser.cpp
Go to the documentation of this file.
00001 // //////////////////////////////////////////////////////////////////////
00002 // Import section
00003 // //////////////////////////////////////////////////////////////////////
00004 // STL
00005 #include <cassert>
00006 #include <string>
00007 #include <sstream>
00008 #include <algorithm>
00009 #include <cmath>
00010 // StdAir
00011 #include <stdair/stdair_basic_types.hpp>
00012 #include <stdair/bom/BomManager.hpp>
00013 #include <stdair/bom/LegCabin.hpp>
00014 #include <stdair/bom/SegmentCabin.hpp>
00015 #include <stdair/bom/BookingClass.hpp>
00016 #include <stdair/bom/VirtualClassStruct.hpp>
00017 #include <stdair/service/Logger.hpp>
00018 
00019 #include <stdair/basic/RandomGeneration.hpp>
00020 #include <stdair/basic/BasConst_General.hpp>
00021 // RMOL
00022 #include <rmol/bom/MCOptimiser.hpp>
00023 
00024 namespace RMOL {
00025 
00026   // // //////////////////////////////////////////////////////////////////////
00027   void MCOptimiser::
00028   optimalOptimisationByMCIntegration (stdair::LegCabin& ioLegCabin) { 
00029     // Retrieve the segment-cabin
00030     const stdair::SegmentCabinList_T lSegmentCabinList =
00031       stdair::BomManager::getList<stdair::SegmentCabin> (ioLegCabin);
00032     stdair::SegmentCabinList_T::const_iterator itSC = lSegmentCabinList.begin();
00033     assert (itSC != lSegmentCabinList.end());
00034     const stdair::SegmentCabin* lSegmentCabin_ptr = *itSC;
00035     assert (lSegmentCabin_ptr != NULL);
00036     
00037     // Retrieve the class list.
00038     const stdair::BookingClassList_T lBookingClassList =
00039       stdair::BomManager::getList<stdair::BookingClass> (*lSegmentCabin_ptr);
00040 
00041     // Retrieve the remaining cabin capacity.
00042     const stdair::Availability_T& lCap = ioLegCabin.getAvailabilityPool();
00043     const int lCapacity = static_cast<const int> (lCap);
00044     const stdair::UnsignedIndex_T lCapacityIndex =
00045       static_cast<const stdair::UnsignedIndex_T> ((lCapacity+abs(lCapacity))/2);
00046     
00047     // Retrieve the virtual class list.
00048     stdair::VirtualClassList_T& lVCList = ioLegCabin.getVirtualClassList();
00049 
00050     // Parse the virtual class list and compute the protection levels.
00051     stdair::VirtualClassList_T::iterator itCurrentVC = lVCList.begin();
00052     assert (itCurrentVC != lVCList.end());
00053     stdair::VirtualClassList_T::iterator itNextVC = itCurrentVC; ++itNextVC;
00054 
00055     // Initialise  the partial sum holder with the demand sample of the first
00056     // virtual class.
00057     stdair::VirtualClassStruct& lFirstVC = *itCurrentVC;
00058     stdair::GeneratedDemandVector_T lPartialSumHolder =
00059       lFirstVC.getGeneratedDemandVector();      
00060 
00061     // Initialise the booking limit for the first class, which is equal to
00062     // the remaining capacity.
00063     lFirstVC.setCumulatedBookingLimit (lCap);
00064 
00065     // Initialise bid price vector with the first element (index 0) equal to
00066     // the highest yield.
00067     ioLegCabin.emptyBidPriceVector();
00068     stdair::BidPriceVector_T& lBPV = ioLegCabin.getBidPriceVector();
00069     //const stdair::Yield_T& y1 = lFirstVC.getYield ();
00070     //lBPV.push_back (y1);
00071     stdair::UnsignedIndex_T idx = 1;    
00072     
00073     for (; itNextVC != lVCList.end(); ++itCurrentVC, ++itNextVC) {
00074       // Get the yields of the two classes.
00075       stdair::VirtualClassStruct& lCurrentVC = *itCurrentVC;
00076       stdair::VirtualClassStruct& lNextVC = *itNextVC;
00077       const stdair::Yield_T& yj = lCurrentVC.getYield ();
00078       const stdair::Yield_T& yj1 = lNextVC.getYield ();
00079 
00080       // Consistency check: the yield/price of a higher class/bucket 
00081       // (with the j index lower) must be higher.
00082       assert (yj > yj1);
00083 
00084       // Sort the partial sum holder.
00085       std::sort (lPartialSumHolder.begin(), lPartialSumHolder.end());
00086       const stdair::UnsignedIndex_T K = lPartialSumHolder.size ();
00087 
00088       // Compute the optimal index lj = floor {[y(j)-y(j+1)]/y(j) . K}
00089       const double ljdouble = std::floor (K * (yj - yj1) / yj);
00090       stdair::UnsignedIndex_T lj =
00091         static_cast<stdair::UnsignedIndex_T> (ljdouble);
00092       
00093       // Consistency check. 
00094       assert (lj >= 1 && lj < K);
00095 
00096       //  The optimal protection: p(j) = 1/2 [S(j,lj) + S(j, lj+1)]
00097       const double sjl = lPartialSumHolder.at (lj - 1);
00098       const double sjlp1 = lPartialSumHolder.at (lj + 1 - 1);
00099       const double pj = (sjl + sjlp1) / 2;
00100 
00101       // Set the cumulated protection level for the current class.
00102       lCurrentVC.setCumulatedProtection (pj);
00103       // Set the cumulated booking limit for the next class.
00104       lNextVC.setCumulatedBookingLimit (lCap - pj);
00105 
00110       const stdair::UnsignedIndex_T pjint = static_cast<const int> (pj);
00111       stdair::GeneratedDemandVector_T::iterator itLowerBound =
00112         lPartialSumHolder.begin();
00113       for (; idx <= pjint && idx <= lCapacityIndex; ++idx) {
00114         itLowerBound =
00115           std::lower_bound (itLowerBound, lPartialSumHolder.end(), idx);
00116         const stdair::UnsignedIndex_T pos =
00117           itLowerBound - lPartialSumHolder.begin();
00118 
00119         const stdair::BidPrice_T lBP = yj * (K - pos) / K;
00120         lBPV.push_back (lBP);
00121       }
00122 
00123       // Update the partial sum holder.
00124       const stdair::GeneratedDemandVector_T& lNextPSH =
00125         lNextVC.getGeneratedDemandVector();
00126       assert (K <= lNextPSH.size());
00127       for (stdair::UnsignedIndex_T i = 0; i < K - lj; ++i) {
00128         lPartialSumHolder.at(i) = lPartialSumHolder.at(i + lj) + lNextPSH.at(i);
00129       }
00130       lPartialSumHolder.resize (K - lj);
00131     }
00132     
00137     stdair::VirtualClassStruct& lLastVC = *itCurrentVC;
00138     const stdair::Yield_T& yn = lLastVC.getYield();
00139     stdair::GeneratedDemandVector_T::iterator itLowerBound =
00140       lPartialSumHolder.begin();
00141     for (; idx <= lCapacityIndex; ++idx) {
00142       itLowerBound =
00143         std::lower_bound (itLowerBound, lPartialSumHolder.end(), idx);
00144       const stdair::UnsignedIndex_T pos =
00145         itLowerBound - lPartialSumHolder.begin();
00146       const stdair::UnsignedIndex_T K = lPartialSumHolder.size();
00147       const stdair::BidPrice_T lBP = yn * (K - pos) / K;
00148       lBPV.push_back (lBP);
00149     }
00150   }
00151 
00152   // ///////////////////////////////////////////////////////////////////
00153   stdair::GeneratedDemandVector_T MCOptimiser::
00154   generateDemandVector (const stdair::MeanValue_T& iMean,
00155                         const stdair::StdDevValue_T& iStdDev,
00156                         const unsigned int& K) {
00157     stdair::GeneratedDemandVector_T oDemandVector;
00158     if (iStdDev > 0) {
00159       stdair::RandomGeneration lGenerator (stdair::DEFAULT_RANDOM_SEED);
00160       for (unsigned int i = 0; i < K; ++i) {
00161         stdair::RealNumber_T lDemandSample =
00162           lGenerator.generateNormal (iMean, iStdDev);
00163         oDemandVector.push_back (lDemandSample);        
00164       }
00165     } else {
00166       for (unsigned int i = 0; i < K; ++i) {
00167         oDemandVector.push_back (iMean);
00168       }
00169     }
00170     return oDemandVector;
00171   }
00172 
00173   // /////////////////////////////////////////////////////////////////////////
00174   void MCOptimiser::
00175   optimisationByMCIntegration (stdair::LegCabin& ioLegCabin) {
00176     // Number of MC samples
00177     unsigned int K = 100000;
00178 
00179     const stdair::YieldLevelDemandMap_T& lYieldDemandMap =
00180       ioLegCabin.getYieldLevelDemandMap();
00181     assert (!lYieldDemandMap.empty());
00182 
00183     std::ostringstream oStr;
00184     oStr << "Yield list ";
00185     for (stdair::YieldLevelDemandMap_T::const_iterator itYD =
00186            lYieldDemandMap.begin();
00187          itYD != lYieldDemandMap.end(); ++itYD) {
00188       const stdair::Yield_T& y = itYD->first;
00189       oStr << y << " ";
00190     }
00191 
00192     STDAIR_LOG_DEBUG (oStr.str());
00193     ioLegCabin.emptyBidPriceVector();
00194     stdair::BidPriceVector_T& lBidPriceVector =
00195       ioLegCabin.getBidPriceVector();
00196     const stdair::Availability_T& lAvailabilityPool =
00197       ioLegCabin.getAvailabilityPool();
00198     // Initialise the minimal bid price to 1.0 (just to avoid problems
00199     // of division by zero).
00200     const stdair::BidPrice_T& lMinBP = 1.0;
00201 
00202     stdair::YieldLevelDemandMap_T::const_reverse_iterator itCurrentYD =
00203       lYieldDemandMap.rbegin();
00204     stdair::YieldLevelDemandMap_T::const_reverse_iterator itNextYD = itCurrentYD;
00205     ++itNextYD;
00206     
00207     // Initialise the partial sum holder
00208     stdair::MeanStdDevPair_T lMeanStdDevPair = itCurrentYD->second;
00209     stdair::GeneratedDemandVector_T lPartialSumHolder =
00210       generateDemandVector(lMeanStdDevPair.first, lMeanStdDevPair.second, K);
00211     
00212     stdair::UnsignedIndex_T idx = 1;
00213     for (; itNextYD!=lYieldDemandMap.rend(); ++itCurrentYD, ++itNextYD) {
00214       const stdair::Yield_T& yj = itCurrentYD->first;
00215       const stdair::Yield_T& yj1 = itNextYD->first;      
00216       // Consistency check: the yield/price of a higher class/bucket 
00217       // (with the j index lower) must be higher.
00218       assert (yj > yj1);
00219       // Sort the partial sum holder.
00220       std::sort (lPartialSumHolder.begin(), lPartialSumHolder.end());
00221       // STDAIR_LOG_DEBUG ("Partial sums : max = " << lPartialSumHolder.back()
00222       //                   << " min = " << lPartialSumHolder.front());
00223       K = lPartialSumHolder.size ();
00224       // Compute the optimal index lj = floor {[y(j)-y(j+1)]/y(j) . K}
00225       const double ljdouble = std::floor (K * (yj - yj1) / yj);
00226       stdair::UnsignedIndex_T lj =
00227         static_cast<stdair::UnsignedIndex_T> (ljdouble);
00228       // Consistency check. 
00229       assert (lj >= 1 && lj < K);
00230       //  The optimal protection: p(j) = 1/2 [S(j,lj) + S(j, lj+1)]
00231       const double sjl = lPartialSumHolder.at (lj - 1);
00232       const double sjlp1 = lPartialSumHolder.at (lj + 1 - 1);
00233       const double pj = (sjl + sjlp1) / 2;
00238       const stdair::UnsignedIndex_T pjint = static_cast<const int> (pj);
00239       stdair::GeneratedDemandVector_T::iterator itLowerBound =
00240         lPartialSumHolder.begin();
00241       for (; idx <= pjint && idx <= lAvailabilityPool; ++idx) {
00242         itLowerBound =
00243           std::lower_bound (itLowerBound, lPartialSumHolder.end(), idx);
00244         const stdair::UnsignedIndex_T pos =
00245           itLowerBound - lPartialSumHolder.begin();        
00246 
00247         const stdair::BidPrice_T lBP = yj * (K - pos) / K;
00248         lBidPriceVector.push_back (lBP);
00249       }
00250       // Update the partial sum holder.
00251       lMeanStdDevPair = itNextYD->second;
00252       const stdair::GeneratedDemandVector_T& lNextDV =
00253         generateDemandVector (lMeanStdDevPair.first,
00254                               lMeanStdDevPair.second, K - lj);
00255       for (stdair::UnsignedIndex_T i = 0; i < K - lj; ++i) {
00256         lPartialSumHolder.at(i) = lPartialSumHolder.at(i + lj) + lNextDV.at(i);
00257       }
00258       lPartialSumHolder.resize (K - lj);      
00259     }
00268     // STDAIR_LOG_DEBUG ("Partial sums : max = " << lPartialSumHolder.back()
00269     //                   << " min = " << lPartialSumHolder.front());
00270     
00271     std::sort (lPartialSumHolder.begin(), lPartialSumHolder.end());
00272     const stdair::Yield_T& yn = itCurrentYD->first;
00273     stdair::GeneratedDemandVector_T::iterator itLowerBound =
00274       lPartialSumHolder.begin();
00275     K = lPartialSumHolder.size();
00276     
00277     bool lMinBPReached = false;
00278     for (; idx <= lAvailabilityPool; ++idx) {
00279       itLowerBound =
00280         std::lower_bound (itLowerBound, lPartialSumHolder.end(), idx);
00281       
00282       if (!lMinBPReached) {
00283         const stdair::UnsignedIndex_T pos =
00284           itLowerBound - lPartialSumHolder.begin();      
00285         stdair::BidPrice_T lBP = yn * (K - pos) / K;
00286         
00287         if (lBP < lMinBP) {
00288           lBP = lMinBP; lMinBPReached = true;
00289         }       
00290         
00291         lBidPriceVector.push_back (lBP);
00292         
00293       } else {        
00294         lBidPriceVector.push_back (lMinBP);
00295       }      
00296       
00297     }
00298     
00299     // Updating the bid price values
00300     ioLegCabin.updatePreviousBidPrice();
00301     ioLegCabin.setCurrentBidPrice (lBidPriceVector.back());
00302 
00303     // Compute and display the bid price variation after optimisation
00304     const stdair::BidPrice_T lPreviousBP = ioLegCabin.getPreviousBidPrice();
00305     stdair::BidPrice_T lNewBP = ioLegCabin.getCurrentBidPrice();
00306     // Check
00307     assert (lPreviousBP != 0);
00308     stdair::BidPrice_T lBidPriceDelta = lNewBP - lPreviousBP;
00309         
00310     double lBidPriceVariation = 100*lBidPriceDelta/lPreviousBP;
00311     
00312     STDAIR_LOG_DEBUG ("Bid price: previous value " << lPreviousBP
00313                       << ", new value " << lNewBP
00314                       << ", variation " << lBidPriceVariation << " %"
00315                       <<", BPV size " << lBidPriceVector.size());
00316   }
00317   
00318 }