00001 // ////////////////////////////////////////////////////////////////////// 00002 // Import section 00003 // ////////////////////////////////////////////////////////////////////// 00004 // GSL Random Number Generation (GSL Reference Manual, version 1.7, Chapter 19) 00005 #include <gsl/gsl_cdf.h> 00006 #include <gsl/gsl_randist.h> 00007 // C 00008 #include <assert.h> 00009 // STL 00010 #include <iostream> 00011 #include <cmath> 00012 #include <vector> 00013 // RMOL 00014 #include <rmol/basic/BasConst_General.hpp> 00015 #include <rmol/bom/DPOptimiser.hpp> 00016 #include <rmol/bom/Bucket.hpp> 00017 #include <rmol/bom/BucketHolder.hpp> 00018 00019 namespace RMOL { 00020 00021 // //////////////////////////////////////////////////////////////////// 00022 void DPOptimiser:: 00023 optimalOptimisationByDP (const ResourceCapacity_T iCabinCapacity, 00024 BucketHolder& ioBucketHolder, 00025 BidPriceVector_T& ioBidPriceVector) { 00026 // Number of classes/buckets: n 00027 const short nbOfClasses = ioBucketHolder.getSize(); 00028 00029 // Number of values of x to compute for each Vj(x). 00030 const int maxValue = static_cast<int> (iCabinCapacity * DEFAULT_PRECISION); 00031 00032 // Vector of the Expected Maximal Revenue (Vj). 00033 std::vector< std::vector<double> > MERVectorHolder; 00034 00035 // Vector of V_0(x). 00036 std::vector<double> initialMERVector (maxValue+1, 0.0); 00037 MERVectorHolder.push_back (initialMERVector); 00038 00039 // Current cumulative protection level (y_j * DEFAULT_PRECISION). 00040 // Initialise with y_0 = 0. 00041 int currentProtection = 0; 00042 00043 int currentBucketIndex = 1; 00044 ioBucketHolder.begin(); 00045 00046 while (currentProtection < maxValue && currentBucketIndex < nbOfClasses) { 00047 //while (currentBucketIndex == 1) { 00048 bool protectionChanged = false; 00049 double nextProtection = 0.0; 00050 std::vector<double> currentMERVector; 00051 // double testGradient = 10000; 00052 00053 Bucket& currentBucket = ioBucketHolder.getCurrentBucket(); 00054 const double meanDemand = currentBucket.getMean(); 00055 const double SDDemand = currentBucket.getStandardDeviation(); 00056 const double currentYield = currentBucket.getAverageYield(); 00057 const double errorFactor = 1;//gsl_cdf_gaussian_Q (-meanDemand, SDDemand); 00058 00059 Bucket& nextBucket = ioBucketHolder.getNextBucket(); 00060 const double nextYield = nextBucket.getAverageYield(); 00061 00062 // For x <= currentProtection (y_(j-1)), V_j(x) = V_(j-1)(x). 00063 for (int x = 0; x <= currentProtection; ++x) { 00064 const double MERValue = MERVectorHolder.at(currentBucketIndex-1).at(x); 00065 currentMERVector.push_back(MERValue); 00066 } 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 gsl_ran_gaussian_pdf (s/DEFAULT_PRECISION - meanDemand, SDDemand); 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 cdfGaussianQ (s/DEFAULT_PRECISION - meanDemand, SDDemand); 00081 cdfVector.push_back(cdfValue); 00082 } 00083 00084 // Compute V_j(x) for x > currentProtection (y_(j-1)). 00085 for (int x = currentProtection + 1; x <= maxValue; ++x) { 00086 const double lowerBound = static_cast<double> (x - currentProtection); 00087 00088 // Compute the first integral in the V_j(x) formulation (see 00089 // the memo of Jerome Contant). 00090 const double power1 = - 0.5 * meanDemand * meanDemand / 00091 (SDDemand * SDDemand); 00092 const double e1 = exp (power1); 00093 const double power2 = 00094 - 0.5 * (lowerBound / DEFAULT_PRECISION - meanDemand) * 00095 (lowerBound / DEFAULT_PRECISION - meanDemand) / 00096 (SDDemand * SDDemand); 00097 const double e2 = exp (power2); 00098 /* 00099 const double integralResult1 = currentYield * 00100 ((e1 - e2) * SDDemand / sqrt (2 * 3.14159265) + 00101 meanDemand * gsl_cdf_gaussian_Q (-meanDemand, SDDemand) - 00102 meanDemand * gsl_cdf_gaussian_Q (lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand)); 00103 */ 00104 const double integralResult1 = currentYield * 00105 ((e1 - e2) * SDDemand / sqrt (2 * 3.14159265) + 00106 meanDemand * cdfGaussianQ (-meanDemand, SDDemand) - 00107 meanDemand * cdfGaussianQ (lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand)); 00108 00109 double integralResult2 = 0.0; 00110 00111 for (int s = 0; s < lowerBound; ++s) { 00112 const double partialResult = 00113 2 * MERVectorHolder.at(currentBucketIndex-1).at(x-s) * 00114 pdfVector.at(s); 00115 00116 integralResult2 += partialResult; 00117 } 00118 integralResult2 -= MERVectorHolder.at(currentBucketIndex-1).at(x) * 00119 pdfVector.at(0); 00120 00121 const int intLowerBound = static_cast<int>(lowerBound); 00122 integralResult2 += 00123 MERVectorHolder.at(currentBucketIndex-1).at(x - intLowerBound) * 00124 pdfVector.at(intLowerBound); 00125 00126 integralResult2 /= 2 * DEFAULT_PRECISION; 00127 /* 00128 for (int s = 0; s < lowerBound; ++s) { 00129 const double partialResult = 00130 (MERVectorHolder.at(currentBucketIndex-1).at(x-s) + 00131 MERVectorHolder.at(currentBucketIndex-1).at(x-s-1)) * 00132 (cdfVector.at(s+1) - cdfVector.at(s)) / 2; 00133 integralResult2 += partialResult; 00134 } 00135 */ 00136 const double firstElement = integralResult1 + integralResult2; 00137 00138 // Compute the second integral in the V_j(x) formulation (see 00139 // the memo of Jerome Contant). 00140 const double constCoefOfSecondElement = 00141 currentYield * lowerBound / DEFAULT_PRECISION + 00142 MERVectorHolder.at(currentBucketIndex-1).at(currentProtection); 00143 const double secondElement = constCoefOfSecondElement * 00144 //gsl_cdf_gaussian_Q(lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand); 00145 cdfGaussianQ (lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand); 00146 const double MERValue = (firstElement + secondElement) / errorFactor; 00147 00148 00149 assert (currentMERVector.size() > 0); 00150 const double lastMERValue = currentMERVector.back(); 00151 00152 const double currentGradient = 00153 (MERValue - lastMERValue) * DEFAULT_PRECISION; 00154 00155 //assert (currentGradient >= 0); 00156 if (currentGradient < -0) { 00157 std::cout << currentGradient << std::endl 00158 << "x = " << x << std::endl 00159 << "class: " << currentBucketIndex << std::endl; 00160 } 00161 00162 /* 00163 assert (currentGradient <= testGradient); 00164 testGradient = currentGradient; 00165 */ 00166 if (!protectionChanged && currentGradient <= nextYield) { 00167 nextProtection = x - 1; 00168 protectionChanged = true; 00169 } 00170 00171 if (protectionChanged && currentGradient > nextYield) { 00172 protectionChanged = false; 00173 } 00174 00175 if (!protectionChanged && x == maxValue) { 00176 nextProtection = maxValue; 00177 } 00178 00179 currentMERVector.push_back (MERValue); 00180 } 00181 std::cout << "Vmaxindex = " << currentMERVector.back() << std::endl; 00182 00183 MERVectorHolder.push_back (currentMERVector); 00184 00185 const double realProtection = nextProtection / DEFAULT_PRECISION; 00186 const double bookingLimit = iCabinCapacity - realProtection; 00187 00188 currentBucket.setCumulatedProtection (realProtection); 00189 nextBucket.setCumulatedBookingLimit (bookingLimit); 00190 00191 currentProtection = static_cast<int> (std::floor (nextProtection)); 00192 00193 ioBucketHolder.iterate(); 00194 ++currentBucketIndex; 00195 } 00196 00197 } 00198 00199 // //////////////////////////////////////////////////////////////////// 00200 double DPOptimiser::cdfGaussianQ (const double c, const double sd) { 00201 const double power = - c * c * 0.625 / (sd * sd); 00202 const double e = sqrt (1-exp(power)); 00203 double result; 00204 if (c >= 0) { 00205 result = 0.5 * (1 - e); 00206 } 00207 else { 00208 result = 0.5 * (1 + e); 00209 } 00210 return result; 00211 } 00212 00213 }
Generated on Sun Jun 14 23:33:59 2009 for RMOL by Doxygen 1.5.8