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 }
Generated on Fri Jul 30 21:53:39 2010 for RMOL by Doxygen 1.6.1