RMOL Logo Get Revenue Management Optimisation Library at SourceForge.net. Fast, secure and Free Open Source software downloads

MCOptimiser.cpp

Go to the documentation of this file.
00001 // //////////////////////////////////////////////////////////////////////
00002 // Import section
00003 // //////////////////////////////////////////////////////////////////////
00004 // C
00005 #include <assert.h>
00006 // STL
00007 #include <string>
00008 #include <fstream>
00009 #include <sstream>
00010 #include <cmath>
00011 // RMOL
00012 #include <rmol/basic/BasChronometer.hpp>
00013 #include <rmol/bom/StudyStatManager.hpp>
00014 #include <rmol/bom/VariateList.hpp>
00015 #include <rmol/bom/Gaussian.hpp>
00016 #include <rmol/bom/Bucket.hpp>
00017 #include <rmol/bom/BucketHolder.hpp>
00018 #include <rmol/bom/PartialSumHolder.hpp>
00019 #include <rmol/bom/PartialSumHolderHolder.hpp>
00020 //#include <rmol/bom/Resource.hpp>
00021 #include <rmol/bom/MCOptimiser.hpp>
00022 #include <rmol/service/Logger.hpp>
00023 
00024 namespace RMOL {
00025 
00026   // //////////////////////////////////////////////////////////////////////
00027   void MCOptimiser::
00028   optimalOptimisationByMCIntegration(const int K, 
00029                                      const ResourceCapacity_T iCabinCapacity,
00030                                      BucketHolder& ioBucketHolder,
00031                                      PartialSumHolderHolder& ioPSHolderHolder,
00032                                      BidPriceVector_T& ioBidPriceVector){
00033     // Retrieve the BucketHolder
00034     // BucketHolder& ioBucketHolder = ioResource.getBucketHolder();
00035 
00036     // Number of classes/buckets: n
00037     const short nbOfClasses = ioBucketHolder.getSize();
00038 
00046     ioPSHolderHolder.begin();
00047     PartialSumHolder& firstPartialSumHolder = 
00048       ioPSHolderHolder.getCurrentPartialSumHolder();
00049     firstPartialSumHolder.initSize (K);
00050 
00051     for (int k=1 ; k <= K; k++) {
00052       firstPartialSumHolder.addPartialSum (0.0);
00053     }
00054 
00060     ioBucketHolder.begin();
00061     ioPSHolderHolder.iterate();
00062     int Kj = K;
00063     int lj = 0;
00064     const int cabinCapacityInt = static_cast<int> (iCabinCapacity);
00065     for (short j = 1 ; j <= nbOfClasses - 1; 
00066          ++j, ioBucketHolder.iterate(), ioPSHolderHolder.iterate()) {
00068       Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00069       Bucket& nextBucket = ioBucketHolder.getNextBucket();
00070 
00071       // STEP 1.
00076       const FldDistributionParameters& aDistribParams = 
00077         currentBucket.getDistributionParameters();
00078       const Gaussian gaussianDemandGenerator (aDistribParams);
00079 
00090       VariateList_T aVariateList;
00091 
00092       PartialSumHolder& previousPartialSumList = 
00093         ioPSHolderHolder.getPreviousPartialSumHolder();
00094       PartialSumHolder& currentPartialSumList = 
00095         ioPSHolderHolder.getCurrentPartialSumHolder();
00096       currentPartialSumList.initSize (Kj);
00097       for (int k=1; k <= Kj; k++) {
00098         const double djk = gaussianDemandGenerator.generateVariate();
00099         aVariateList.push_back (djk);
00100 
00110         const double spjm1lpk = 
00111           previousPartialSumList.getPartialSum (lj + k - 1);
00112         const double sjk = spjm1lpk + djk;
00113         currentPartialSumList.addPartialSum (sjk);
00114 
00115         /* DEBUG
00116            RMOL_LOG_DEBUG ("d(" << j << ", " << k << "); " << djk 
00117            << "; S'(" << j-1 << ", " << lj+k << "); " << spjm1lpk
00118            << "; S(" << j << ", " << k << "); " << sjk);
00119         */
00120       }
00121 
00122       // STEP 2.
00126       currentPartialSumList.sort ();
00127 
00129       const double pj = currentBucket.getAverageYield();
00130       const double pj1 = nextBucket.getAverageYield();
00131 
00134       assert (pj > pj1);
00135 
00140       const double ljdouble = std::floor (Kj * (pj - pj1) / pj);
00141       lj = static_cast<int> (ljdouble);
00142 
00149       assert (lj >= 1 && lj < Kj);
00150 
00155       const double sjl = currentPartialSumList.getPartialSum (lj - 1);
00156       const double sjlp1 = currentPartialSumList.getPartialSum (lj + 1 - 1);
00157       const double yj = (sjl + sjlp1) / 2;
00158 
00164       // Set the cumulated protection for Bucket(j) (j ranging from 1 to n-1)
00165       currentBucket.setCumulatedProtection (yj);
00166 
00173       // Get the previous cumulated protection y(j-1).
00174       const double yjm1 = ioBucketHolder.getPreviousCumulatedProtection ();
00175       const int yjm1int = static_cast<int> (yjm1);
00176       const int yjint = static_cast<int> (yj);
00177       int currentIndex = 0;
00178 
00179       for (int x = yjm1int + 1; 
00180            x <= yjint && x <= cabinCapacityInt; ++x) {
00181         currentIndex = currentPartialSumList.getLowerBound (x);
00182         const double bidPriceAtX = pj * (Kj - currentIndex) / Kj;
00183         ioBidPriceVector.push_back (bidPriceAtX);
00184       }   
00185 
00187       Kj = Kj - lj;
00188 
00192     }
00193 
00194     // Set the protection of Bucket(n) to be equal to the capacity
00195     Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00196     currentBucket.setCumulatedProtection (iCabinCapacity);
00197 
00203     // Get the previous cumulated protection y(n-1).
00204     const double ynm1 = ioBucketHolder.getPreviousCumulatedProtection ();
00205 
00206     if (ynm1 < iCabinCapacity) {
00207       // STEP 1.
00208       const FldDistributionParameters& aDistribParams = 
00209         currentBucket.getDistributionParameters();
00210       const Gaussian gaussianDemandGenerator (aDistribParams);
00211       
00212       VariateList_T aVariateList;
00213       
00214       PartialSumHolder& previousPartialSumList = 
00215         ioPSHolderHolder.getPreviousPartialSumHolder();
00216       PartialSumHolder& currentPartialSumList = 
00217         ioPSHolderHolder.getCurrentPartialSumHolder();
00218       currentPartialSumList.initSize (Kj);
00219       
00220       for (int k = 1; k <= Kj; ++k) {
00221         const double djk = gaussianDemandGenerator.generateVariate();
00222         aVariateList.push_back (djk);
00223         
00224         const double spjm1lpk = 
00225           previousPartialSumList.getPartialSum (lj + k - 1);
00226         const double sjk = spjm1lpk + djk;
00227         currentPartialSumList.addPartialSum (sjk);
00228       }
00229       
00230       // STEP 2.
00231       currentPartialSumList.sort ();
00232       
00233       const int ynm1int = static_cast<int> (ynm1);
00234       const double pn = currentBucket.getAverageYield();
00235       int currentIndex = 0;
00236       for (int x = ynm1int + 1; x <= cabinCapacityInt; ++x) {
00237         currentIndex = currentPartialSumList.getLowerBound (x);
00238         const double bidPriceAtX = pn * (Kj - currentIndex) / Kj;
00239         ioBidPriceVector.push_back (bidPriceAtX);
00240       } 
00241     }
00242     
00247     ioBucketHolder.recalculate ();
00248   }
00249 
00250   // //////////////////////////////////////////////////////////////////////
00251   void MCOptimiser::
00252   optimalOptimisationByMCIntegration(const int K, 
00253                                      const ResourceCapacity_T iCabinCapacity,
00254                                      BucketHolder& ioBucketHolder,
00255                                      PartialSumHolderHolder& ioPSHolderHolder,
00256                                      BidPriceVector_T& ioBidPriceVector,
00257                                      StudyStatManager& ioStudyStatManager){
00258     // Retrieve the BucketHolder
00259     // BucketHolder& ioBucketHolder = ioResource.getBucketHolder();
00260 
00261     BasChronometer lDrawBasChronometer;
00262     BasChronometer lSortBasChronometer;
00263     BasChronometer lBVPCalculationBasChronometer;
00264 
00265     // Number of classes/buckets: n
00266     const short nbOfClasses = ioBucketHolder.getSize();
00267 
00275     ioPSHolderHolder.begin();
00276     PartialSumHolder& firstPartialSumHolder = 
00277       ioPSHolderHolder.getCurrentPartialSumHolder();
00278     firstPartialSumHolder.initSize (K);
00279 
00280     for (int k=1 ; k <= K; k++) {
00281       firstPartialSumHolder.addPartialSum (0.0);
00282     }
00283 
00289     ioBucketHolder.begin();
00290     ioPSHolderHolder.iterate();
00291     int Kj = K;
00292     int lj = 0;
00293     const int cabinCapacityInt = static_cast<int> (iCabinCapacity);
00294     for (short j = 1 ; j <= nbOfClasses - 1; 
00295          ++j, ioBucketHolder.iterate(), ioPSHolderHolder.iterate()) {
00296       // DEBUG
00297       std::cout << "K" << j << " = " << Kj << std::endl;
00298 
00300       Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00301       Bucket& nextBucket = ioBucketHolder.getNextBucket();
00302 
00303       // STEP 1.
00308       const FldDistributionParameters& aDistribParams = 
00309         currentBucket.getDistributionParameters();
00310       const Gaussian gaussianDemandGenerator (aDistribParams);
00311 
00322       VariateList_T aVariateList;
00323 
00324       PartialSumHolder& previousPartialSumList = 
00325         ioPSHolderHolder.getPreviousPartialSumHolder();
00326       PartialSumHolder& currentPartialSumList = 
00327         ioPSHolderHolder.getCurrentPartialSumHolder();
00328       currentPartialSumList.initSize (Kj);
00329       lDrawBasChronometer.start();
00330       for (int k=1; k <= Kj; k++) {
00331         const double djk = gaussianDemandGenerator.generateVariate();
00332         aVariateList.push_back (djk);
00333 
00343         const double spjm1lpk = 
00344           previousPartialSumList.getPartialSum (lj + k - 1);
00345         const double sjk = spjm1lpk + djk;
00346         currentPartialSumList.addPartialSum (sjk);
00347 
00348         /* DEBUG
00349            RMOL_LOG_DEBUG ("d(" << j << ", " << k << "); " << djk 
00350            << "; S'(" << j-1 << ", " << lj+k << "); " << spjm1lpk
00351            << "; S(" << j << ", " << k << "); " << sjk);
00352         */
00353       }
00354 
00355       const double lDrawTimeValue = lDrawBasChronometer.elapsed();
00356       ioStudyStatManager.addMeasure ("DrawTime", lDrawTimeValue);
00357       
00358       // STEP 2.
00362       lSortBasChronometer.start();
00363       currentPartialSumList.sort ();
00364       const double lSortTimeValue = lSortBasChronometer.elapsed();
00365       ioStudyStatManager.addMeasure ("SortTime", lSortTimeValue);
00366       
00368       const double pj = currentBucket.getAverageYield();
00369       const double pj1 = nextBucket.getAverageYield();
00370 
00373       assert (pj > pj1);
00374 
00379       const double ljdouble = std::floor (Kj * (pj - pj1) / pj);
00380       lj = static_cast<int> (ljdouble);
00381       
00382       // DEBUG
00383       std::cout << "l" << j << " = " << lj << std::endl;
00384 
00391       assert (lj >= 1 && lj < Kj);
00392 
00397       const double sjl = currentPartialSumList.getPartialSum (lj - 1);
00398       const double sjlp1 = currentPartialSumList.getPartialSum (lj + 1 - 1);
00399       const double yj = (sjl + sjlp1) / 2;
00400 
00406       // Set the cumulated protection for Bucket(j) (j ranging from 1 to n-1)
00407       currentBucket.setCumulatedProtection (yj);
00408 
00415       // Get the previous cumulated protection y(j-1).
00416       const double yjm1 = ioBucketHolder.getPreviousCumulatedProtection ();
00417       const int yjm1int = static_cast<int> (yjm1);
00418       const int yjint = static_cast<int> (yj);
00419       int currentIndex = 0;
00420 
00421       lBVPCalculationBasChronometer.start();
00422       for (int x = yjm1int + 1; 
00423            x <= yjint && x <= cabinCapacityInt; ++x) {
00424         currentIndex = currentPartialSumList.getLowerBound (x);
00425         const double bidPriceAtX = pj * (Kj - currentIndex) / Kj;
00426         ioBidPriceVector.push_back (bidPriceAtX);
00427       }   
00428       const double lBVPCalculationTimeValue =
00429         lBVPCalculationBasChronometer.elapsed();
00430       ioStudyStatManager.addMeasure ("BVPCalculationTime",
00431                                      lBVPCalculationTimeValue);
00432 
00434       Kj = Kj - lj;
00435 
00439     }
00440 
00441     // Set the protection of Bucket(n) to be equal to the capacity
00442     Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00443     currentBucket.setCumulatedProtection (iCabinCapacity);
00444 
00450     // Get the previous cumulated protection y(n-1).
00451     const double ynm1 = ioBucketHolder.getPreviousCumulatedProtection ();
00452 
00453     if (ynm1 < iCabinCapacity) {
00454       // STEP 1.
00455       const FldDistributionParameters& aDistribParams = 
00456         currentBucket.getDistributionParameters();
00457       const Gaussian gaussianDemandGenerator (aDistribParams);
00458       
00459       VariateList_T aVariateList;
00460       
00461       PartialSumHolder& previousPartialSumList = 
00462         ioPSHolderHolder.getPreviousPartialSumHolder();
00463       PartialSumHolder& currentPartialSumList = 
00464         ioPSHolderHolder.getCurrentPartialSumHolder();
00465       currentPartialSumList.initSize (Kj);
00466       
00467       lDrawBasChronometer.start();
00468       for (int k = 1; k <= Kj; ++k) {
00469         const double djk = gaussianDemandGenerator.generateVariate();
00470         aVariateList.push_back (djk);
00471         
00472         const double spjm1lpk = 
00473           previousPartialSumList.getPartialSum (lj + k - 1);
00474         const double sjk = spjm1lpk + djk;
00475         currentPartialSumList.addPartialSum (sjk);
00476       }
00477       
00478       const double lDrawTimeValue = lDrawBasChronometer.elapsed();
00479       ioStudyStatManager.addMeasure ("DrawTime", lDrawTimeValue);
00480       
00481       // STEP 2.
00482       lSortBasChronometer.start();
00483       currentPartialSumList.sort ();
00484       const double lSortTimeValue = lSortBasChronometer.elapsed();
00485       ioStudyStatManager.addMeasure ("SortTime", lSortTimeValue);
00486       
00487       const int ynm1int = static_cast<int> (ynm1);
00488       const double pn = currentBucket.getAverageYield();
00489       int currentIndex = 0;
00490       
00491       lBVPCalculationBasChronometer.start();
00492       for (int x = ynm1int + 1; x <= cabinCapacityInt; ++x) {
00493         currentIndex = currentPartialSumList.getLowerBound (x);
00494         const double bidPriceAtX = pn * (Kj - currentIndex) / Kj;
00495         ioBidPriceVector.push_back (bidPriceAtX);
00496       }
00497       const double lBVPCalculationTimeValue =
00498         lBVPCalculationBasChronometer.elapsed();
00499       ioStudyStatManager.addMeasure ("BVPCalculationTime",
00500                                      lBVPCalculationTimeValue);
00501     }
00502     
00507     ioBucketHolder.recalculate ();
00508   }
00509 
00510   // ////////////////////////////////////////////////////////////////////////
00511   void MCOptimiser::
00512   legOptimisationByMC (const ResourceCapacity_T iCabinCapacity,
00513                        BucketHolder& ioBucketHolder,
00514                        BidPriceVector_T& ioBidPriceVector) {
00515 
00516     ioBucketHolder.begin();
00517 
00518     // Get the first bucket (the one with the highest average yield).
00519     Bucket& lFirstBucket = ioBucketHolder.getCurrentBucket();
00520 
00521     GeneratedDemandVector_T lPartialSumVector =
00522       lFirstBucket.getGeneratedDemandVector ();
00523 
00524     // Sort the vector from high to low.
00525     std::sort (lPartialSumVector.begin(), lPartialSumVector.end(),
00526                std::greater<double>());
00527 
00528     // Get the number of draws (K).
00529     const unsigned int K = lPartialSumVector.size();
00530 
00531     // Number of classes/buckets: n
00532     const short nbOfClasses = ioBucketHolder.getSize();
00533     
00539     unsigned int Kj = K;
00540     const int cabinCapacityInt = static_cast<int> (iCabinCapacity);
00541     for (short j = 1 ; j <= nbOfClasses - 1; ++j, ioBucketHolder.iterate()) {
00542       // DEBUG
00543       std::cout << "K" << j << " = " << Kj << std::endl;
00544       
00546       Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00547       Bucket& nextBucket = ioBucketHolder.getNextBucket();
00548 
00550       const double pj = currentBucket.getAverageYield();
00551       const double pj1 = nextBucket.getAverageYield();
00552 
00557       const unsigned int lj = Kj - std::floor (Kj * (pj - pj1) / pj);
00558 
00559       // DEBUG
00560       std::cout << "l" << j << " = " << lj << std::endl;
00561 
00562       /*
00563       std::cout << "p(j+1) = " << pj1 << std::endl
00564                 << "p(j) = " << pj << std::endl
00565                 << "Kj = " << Kj << std::endl;
00566       */
00567       
00569       assert (lj >= 1 && lj < Kj);
00570 
00575       const double sjl = lPartialSumVector.at (lj - 1);
00576       const double sjlp1 = lPartialSumVector.at (lj + 1 - 1);
00577       const double yj = (sjl + sjlp1) / 2;
00578 
00579       // Set the cumulated protection for Bucket(j) (j ranging from 1 to n-1)
00580       currentBucket.setCumulatedProtection (yj);
00581 
00583       Kj = lj;
00584       lPartialSumVector.resize (Kj);
00585 
00586       // Generated demand of the (j+1)th bucket for the next iteration.
00587       const GeneratedDemandVector_T& lNextGeneratedDemandVector =
00588         nextBucket.getGeneratedDemandVector ();
00589 
00590       for (unsigned int i = 0; i < Kj; ++i) {
00591         const double lGeneratedDemand = lNextGeneratedDemandVector.at(i);
00592         lPartialSumVector.at(i) += lGeneratedDemand;
00593       }
00594 
00595       // Sort the vector from high to low.
00596       std::sort (lPartialSumVector.begin(), lPartialSumVector.end(),
00597                  std::greater<double>());
00598     }
00599 
00600     // Set the protection of Bucket(n) to be equal to the capacity
00601     Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00602     currentBucket.setCumulatedProtection (iCabinCapacity);
00603 
00608     ioBucketHolder.recalculate ();
00609   }
00610 
00611 }
SourceForge Logo

Generated on Sat Jun 6 13:49:02 2009 for RMOL by Doxygen 1.5.7.1