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