Go to the documentation of this file.00001
00002
00003
00004
00005 #include <assert.h>
00006
00007 #include <string>
00008 #include <fstream>
00009 #include <sstream>
00010 #include <cmath>
00011
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
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
00034
00035
00036
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
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
00116
00117
00118
00119
00120 }
00121
00122
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
00165 currentBucket.setCumulatedProtection (yj);
00166
00173
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
00195 Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00196 currentBucket.setCumulatedProtection (iCabinCapacity);
00197
00203
00204 const double ynm1 = ioBucketHolder.getPreviousCumulatedProtection ();
00205
00206 if (ynm1 < iCabinCapacity) {
00207
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
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
00259
00260
00261 BasChronometer lDrawBasChronometer;
00262 BasChronometer lSortBasChronometer;
00263 BasChronometer lBVPCalculationBasChronometer;
00264
00265
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
00297 std::cout << "K" << j << " = " << Kj << std::endl;
00298
00300 Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00301 Bucket& nextBucket = ioBucketHolder.getNextBucket();
00302
00303
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
00349
00350
00351
00352
00353 }
00354
00355 const double lDrawTimeValue = lDrawBasChronometer.elapsed();
00356 ioStudyStatManager.addMeasure ("DrawTime", lDrawTimeValue);
00357
00358
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
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
00407 currentBucket.setCumulatedProtection (yj);
00408
00415
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
00442 Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00443 currentBucket.setCumulatedProtection (iCabinCapacity);
00444
00450
00451 const double ynm1 = ioBucketHolder.getPreviousCumulatedProtection ();
00452
00453 if (ynm1 < iCabinCapacity) {
00454
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
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
00519 Bucket& lFirstBucket = ioBucketHolder.getCurrentBucket();
00520
00521 GeneratedDemandVector_T lPartialSumVector =
00522 lFirstBucket.getGeneratedDemandVector ();
00523
00524
00525 std::sort (lPartialSumVector.begin(), lPartialSumVector.end(),
00526 std::greater<double>());
00527
00528
00529 const unsigned int K = lPartialSumVector.size();
00530
00531
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
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
00560 std::cout << "l" << j << " = " << lj << std::endl;
00561
00562
00563
00564
00565
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
00580 currentBucket.setCumulatedProtection (yj);
00581
00583 Kj = lj;
00584 lPartialSumVector.resize (Kj);
00585
00586
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
00596 std::sort (lPartialSumVector.begin(), lPartialSumVector.end(),
00597 std::greater<double>());
00598 }
00599
00600
00601 Bucket& currentBucket = ioBucketHolder.getCurrentBucket();
00602 currentBucket.setCumulatedProtection (iCabinCapacity);
00603
00608 ioBucketHolder.recalculate ();
00609 }
00610
00611 }