00001 // ////////////////////////////////////////////////////////////////////// 00002 // Import section 00003 // ////////////////////////////////////////////////////////////////////// 00004 // STL 00005 #include <cassert> 00006 #include <sstream> 00007 #include <vector> 00008 #include <cmath> 00009 // Boost Math 00010 #include <boost/math/distributions/normal.hpp> 00011 // StdAir 00012 #include <stdair/bom/LegCabin.hpp> 00013 #include <stdair/bom/VirtualClassStruct.hpp> 00014 #include <stdair/service/Logger.hpp> 00015 // RMOL 00016 #include <rmol/basic/BasConst_General.hpp> 00017 #include <rmol/bom/DPOptimiser.hpp> 00018 00019 namespace RMOL { 00020 00021 // //////////////////////////////////////////////////////////////////// 00022 void DPOptimiser::optimalOptimisationByDP (stdair::LegCabin& ioLegCabin) { 00023 // // Number of classes/buckets: n 00024 // const short nbOfClasses = ioBucketHolder.getSize(); 00025 00026 // // Number of values of x to compute for each Vj(x). 00027 // const int maxValue = static_cast<int> (iCabinCapacity * DEFAULT_PRECISION); 00028 00029 // // Vector of the Expected Maximal Revenue (Vj). 00030 // std::vector< std::vector<double> > MERVectorHolder; 00031 00032 // // Vector of V_0(x). 00033 // std::vector<double> initialMERVector (maxValue+1, 0.0); 00034 // MERVectorHolder.push_back (initialMERVector); 00035 00036 // // Current cumulative protection level (y_j * DEFAULT_PRECISION). 00037 // // Initialise with y_0 = 0. 00038 // int currentProtection = 0; 00039 00040 // int currentBucketIndex = 1; 00041 // ioBucketHolder.begin(); 00042 00043 // while (currentProtection < maxValue && currentBucketIndex < nbOfClasses) { 00044 // //while (currentBucketIndex == 1) { 00045 // bool protectionChanged = false; 00046 // double nextProtection = 0.0; 00047 // std::vector<double> currentMERVector; 00048 // // double testGradient = 10000; 00049 00050 // Bucket& currentBucket = ioBucketHolder.getCurrentBucket(); 00051 // const double meanDemand = currentBucket.getMean(); 00052 // const double SDDemand = currentBucket.getStandardDeviation(); 00053 // const double currentYield = currentBucket.getAverageYield(); 00054 // const double errorFactor = 1.0; 00055 00056 // Bucket& nextBucket = ioBucketHolder.getNextBucket(); 00057 // const double nextYield = nextBucket.getAverageYield(); 00058 00059 // // For x <= currentProtection (y_(j-1)), V_j(x) = V_(j-1)(x). 00060 // for (int x = 0; x <= currentProtection; ++x) { 00061 // const double MERValue = MERVectorHolder.at(currentBucketIndex-1).at(x); 00062 // currentMERVector.push_back (MERValue); 00063 // } 00064 00065 // // 00066 // boost::math::normal lNormalDistribution (meanDemand, SDDemand); 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 // boost::math::pdf (lNormalDistribution, s/DEFAULT_PRECISION); 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 // boost::math::cdf (boost::math::complement (lNormalDistribution, 00081 // s/DEFAULT_PRECISION)); 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 = 00092 // - 0.5 * meanDemand * meanDemand / (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 = std::exp (power2); 00099 00100 // const double cdfValue0 = 00101 // boost::math::cdf (boost::math::complement (lNormalDistribution, 00102 // 0.0)); 00103 // const double cdfValue1 = 00104 // boost::math::cdf(boost::math::complement(lNormalDistribution, 00105 // lowerBound/DEFAULT_PRECISION)); 00106 // const double integralResult1 = currentYield 00107 // * ((e1 - e2) * SDDemand / sqrt (2 * 3.14159265) 00108 // + meanDemand * (cdfValue0 - cdfValue1)); 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 00145 // const double secondElement = constCoefOfSecondElement 00146 // * boost::math::cdf(boost::math::complement(lNormalDistribution, 00147 // lowerBound/DEFAULT_PRECISION)); 00148 00149 // const double MERValue = (firstElement + secondElement) / errorFactor; 00150 00151 00152 // assert (currentMERVector.size() > 0); 00153 // const double lastMERValue = currentMERVector.back(); 00154 00155 // const double currentGradient = 00156 // (MERValue - lastMERValue) * DEFAULT_PRECISION; 00157 00158 // //assert (currentGradient >= 0); 00159 // if (currentGradient < -0) { 00160 // std::ostringstream ostr; 00161 // ostr << currentGradient << std::endl 00162 // << "x = " << x << std::endl 00163 // << "class: " << currentBucketIndex << std::endl; 00164 // STDAIR_LOG_DEBUG (ostr.str()); 00165 // } 00166 00167 // /* 00168 // assert (currentGradient <= testGradient); 00169 // testGradient = currentGradient; 00170 // */ 00171 // if (protectionChanged == false && currentGradient <= nextYield) { 00172 // nextProtection = x - 1; 00173 // protectionChanged = true; 00174 // } 00175 00176 // if (protectionChanged == true && currentGradient > nextYield) { 00177 // protectionChanged = false; 00178 // } 00179 00180 // if (protectionChanged == false && x == maxValue) { 00181 // nextProtection = maxValue; 00182 // } 00183 00184 // currentMERVector.push_back (MERValue); 00185 // } 00186 00187 // // DEBUG 00188 // STDAIR_LOG_DEBUG ("Vmaxindex = " << currentMERVector.back()); 00189 00190 // MERVectorHolder.push_back (currentMERVector); 00191 00192 // const double realProtection = nextProtection / DEFAULT_PRECISION; 00193 // const double bookingLimit = iCabinCapacity - realProtection; 00194 00195 // currentBucket.setCumulatedProtection (realProtection); 00196 // nextBucket.setCumulatedBookingLimit (bookingLimit); 00197 00198 // currentProtection = static_cast<int> (std::floor (nextProtection)); 00199 00200 // ioBucketHolder.iterate(); 00201 // ++currentBucketIndex; 00202 // } 00203 } 00204 00205 }