Menu
▾
▴
[Rmol-svn-commit] SF.net SVN: rmol:[338] trunk/rmol/rmol
|
From: <den...@us...> - 2010-09-27 11:36:32
|
Revision: 338
http://rmol.svn.sourceforge.net/rmol/?rev=338&view=rev
Author: denis_arnaud
Date: 2010-09-27 11:36:24 +0000 (Mon, 27 Sep 2010)
Log Message:
-----------
[Dev] Replaced calls to the GSL (GNU Scientific Library) by Boost.Math.
Modified Paths:
--------------
trunk/rmol/rmol/bom/DPOptimiser.cpp
trunk/rmol/rmol/bom/EmsrUtils.cpp
trunk/rmol/rmol/bom/ExpectationMaximization.cpp
trunk/rmol/rmol/bom/Gaussian.cpp
trunk/rmol/rmol/bom/Gaussian.hpp
trunk/rmol/rmol/bom/HistoricalBookingHolder.cpp
trunk/rmol/rmol/bom/HistoricalBookingHolderHolder.cpp
trunk/rmol/rmol/bom/HistoricalDataHolder.cpp
trunk/rmol/rmol/bom/HistoricalDataHolderHolder.cpp
trunk/rmol/rmol/bom/MCOptimiser.cpp
trunk/rmol/rmol/bom/Makefile.am
trunk/rmol/rmol/bom/Overbooking.cpp
trunk/rmol/rmol/core/Makefile.am
trunk/rmol/rmol/service/RMOL_ServiceContext.cpp
Modified: trunk/rmol/rmol/bom/DPOptimiser.cpp
===================================================================
--- trunk/rmol/rmol/bom/DPOptimiser.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/DPOptimiser.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -1,14 +1,13 @@
// //////////////////////////////////////////////////////////////////////
// Import section
// //////////////////////////////////////////////////////////////////////
-// GSL Random Number Generation (GSL Reference Manual, version 1.7, Chapter 19)
-#include <gsl/gsl_cdf.h>
-#include <gsl/gsl_randist.h>
// STL
#include <cassert>
#include <sstream>
#include <vector>
#include <cmath>
+// Boost Math
+#include <boost/math/distributions/normal.hpp>
// StdAir
#include <stdair/service/Logger.hpp>
// RMOL
@@ -55,7 +54,7 @@
const double meanDemand = currentBucket.getMean();
const double SDDemand = currentBucket.getStandardDeviation();
const double currentYield = currentBucket.getAverageYield();
- const double errorFactor = 1;//gsl_cdf_gaussian_Q (-meanDemand, SDDemand);
+ const double errorFactor = 1.0;
Bucket& nextBucket = ioBucketHolder.getNextBucket();
const double nextYield = nextBucket.getAverageYield();
@@ -63,23 +62,27 @@
// For x <= currentProtection (y_(j-1)), V_j(x) = V_(j-1)(x).
for (int x = 0; x <= currentProtection; ++x) {
const double MERValue = MERVectorHolder.at(currentBucketIndex-1).at(x);
- currentMERVector.push_back(MERValue);
+ currentMERVector.push_back (MERValue);
}
+
+ //
+ boost::math::normal lNormalDistribution (meanDemand, SDDemand);
// Vector of gaussian pdf values.
std::vector<double> pdfVector;
for (int s = 0; s <= maxValue - currentProtection; ++s) {
const double pdfValue =
- gsl_ran_gaussian_pdf (s/DEFAULT_PRECISION - meanDemand, SDDemand);
- pdfVector.push_back(pdfValue);
+ boost::math::pdf (lNormalDistribution, s/DEFAULT_PRECISION);
+ pdfVector.push_back (pdfValue);
}
// Vector of gaussian cdf values.
std::vector<double> cdfVector;
for (int s = 0; s <= maxValue - currentProtection; ++s) {
const double cdfValue =
- cdfGaussianQ (s/DEFAULT_PRECISION - meanDemand, SDDemand);
- cdfVector.push_back(cdfValue);
+ boost::math::cdf (boost::math::complement (lNormalDistribution,
+ s/DEFAULT_PRECISION));
+ cdfVector.push_back (cdfValue);
}
// Compute V_j(x) for x > currentProtection (y_(j-1)).
@@ -88,31 +91,31 @@
// Compute the first integral in the V_j(x) formulation (see
// the memo of Jerome Contant).
- const double power1 = - 0.5 * meanDemand * meanDemand /
- (SDDemand * SDDemand);
+ const double power1 =
+ - 0.5 * meanDemand * meanDemand / (SDDemand * SDDemand);
const double e1 = std::exp (power1);
const double power2 =
- - 0.5 * (lowerBound / DEFAULT_PRECISION - meanDemand) *
- (lowerBound / DEFAULT_PRECISION - meanDemand) /
- (SDDemand * SDDemand);
- const double e2 = exp (power2);
- /*
- const double integralResult1 = currentYield *
- ((e1 - e2) * SDDemand / sqrt (2 * 3.14159265) +
- meanDemand * gsl_cdf_gaussian_Q (-meanDemand, SDDemand) -
- meanDemand * gsl_cdf_gaussian_Q (lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand));
- */
- const double integralResult1 = currentYield *
- ((e1 - e2) * SDDemand / sqrt (2 * 3.14159265) +
- meanDemand * cdfGaussianQ (-meanDemand, SDDemand) -
- meanDemand * cdfGaussianQ (lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand));
+ - 0.5 * (lowerBound / DEFAULT_PRECISION - meanDemand)
+ * (lowerBound / DEFAULT_PRECISION - meanDemand)
+ / (SDDemand * SDDemand);
+ const double e2 = std::exp (power2);
+
+ const double cdfValue0 =
+ boost::math::cdf (boost::math::complement (lNormalDistribution,
+ 0.0));
+ const double cdfValue1 =
+ boost::math::cdf(boost::math::complement(lNormalDistribution,
+ lowerBound/DEFAULT_PRECISION));
+ const double integralResult1 = currentYield
+ * ((e1 - e2) * SDDemand / sqrt (2 * 3.14159265)
+ + meanDemand * (cdfValue0 - cdfValue1));
double integralResult2 = 0.0;
for (int s = 0; s < lowerBound; ++s) {
const double partialResult =
- 2 * MERVectorHolder.at(currentBucketIndex-1).at(x-s) *
- pdfVector.at(s);
+ 2 * MERVectorHolder.at(currentBucketIndex-1).at(x-s)
+ * pdfVector.at(s);
integralResult2 += partialResult;
}
@@ -139,11 +142,13 @@
// Compute the second integral in the V_j(x) formulation (see
// the memo of Jerome Contant).
const double constCoefOfSecondElement =
- currentYield * lowerBound / DEFAULT_PRECISION +
- MERVectorHolder.at(currentBucketIndex-1).at(currentProtection);
- const double secondElement = constCoefOfSecondElement *
- //gsl_cdf_gaussian_Q(lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand);
- cdfGaussianQ (lowerBound / DEFAULT_PRECISION - meanDemand, SDDemand);
+ currentYield * lowerBound / DEFAULT_PRECISION
+ + MERVectorHolder.at(currentBucketIndex-1).at(currentProtection);
+
+ const double secondElement = constCoefOfSecondElement
+ * boost::math::cdf(boost::math::complement(lNormalDistribution,
+ lowerBound/DEFAULT_PRECISION));
+
const double MERValue = (firstElement + secondElement) / errorFactor;
@@ -200,20 +205,4 @@
}
}
- // ////////////////////////////////////////////////////////////////////
- double DPOptimiser::cdfGaussianQ (const double c, const double sd) {
- const double power = - c * c * 0.625 / (sd * sd);
- const double e = std::sqrt (1 - std::exp (power));
- double result = 0.0;
-
- if (c >= 0) {
- result = 0.5 * (1 - e);
-
- } else {
- result = 0.5 * (1 + e);
- }
-
- return result;
- }
-
}
Modified: trunk/rmol/rmol/bom/EmsrUtils.cpp
===================================================================
--- trunk/rmol/rmol/bom/EmsrUtils.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/EmsrUtils.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -1,11 +1,11 @@
// //////////////////////////////////////////////////////////////////////
// Import section
// //////////////////////////////////////////////////////////////////////
-// GSL Random Number Generation (GSL Reference Manual, version 1.7, Chapter 19)
-#include <gsl/gsl_cdf.h>
-// C
-#include <math.h>
-#include <assert.h>
+// STL
+#include <cassert>
+#include <cmath>
+// Boost Math
+#include <boost/math/distributions/normal.hpp>
// RMOL
#include <rmol/bom/EmsrUtils.hpp>
#include <rmol/bom/Bucket.hpp>
@@ -53,17 +53,24 @@
const double lNextYield = ioNextBucket.getAverageYield();
assert (lAggreatedYield != 0);
- // Compute the protection for the aggregated class/bucket
- const double lProtection =
- lMean + gsl_cdf_gaussian_Qinv (lNextYield/lAggreatedYield, lSD);
+ // Compute the yield ratio between the higher bucket and the current one
+ const double lYieldRatio = lNextYield / lAggreatedYield;
+
+ /** Compute the protection for the aggregated class/bucket.
+ <br>Note: The inverse cdf is the quantile function (see also
+ http://en.wikipedia.org/wiki/Quantile_function). */
+ boost::math::normal lNormalDistribution (lMean, lSD);
+ const double lProtection =
+ boost::math::quantile (boost::math::complement (lNormalDistribution,
+ lYieldRatio));
return lProtection;
}
// ////////////////////////////////////////////////////////////////////
- const double EmsrUtils::computeProtectionLevelwithSellup (Bucket& iHigherBucket,
- Bucket& iBucket,
- double iSellupFactor){
+ const double EmsrUtils::
+ computeProtectionLevelwithSellup (Bucket& iHigherBucket, Bucket& iBucket,
+ double iSellupFactor) {
// Retrieve the mean and the standard deviation of the higher
// class(es)/bucket(s) depending EMSR-a or EMSR-b
// and the average yield of each input classes/buckets
@@ -74,25 +81,34 @@
assert (lHigherBucketYield > DEFAULT_EPSILON);
assert (1-iSellupFactor > DEFAULT_EPSILON);
- // compute the protection level for the higher class/bucket
- const double lProtectionLevel =
- lMean +
- gsl_cdf_gaussian_Pinv((lHigherBucketYield-lBucketYield)/
- (lHigherBucketYield*(1-iSellupFactor)),lSD);
-
+ // Compute the yield ratio
+ const double lYieldRatio = (lHigherBucketYield - lBucketYield)
+ / (lHigherBucketYield * (1 - iSellupFactor));
+
+ /** Compute the protection for the for the higher class/bucket.
+ <br>Note: The inverse cdf is the quantile function (see also
+ http://en.wikipedia.org/wiki/Quantile_function). */
+ boost::math::normal lNormalDistribution (lMean, lSD);
+ const double lProtectionLevel = boost::math::quantile (lNormalDistribution,
+ lYieldRatio);
+
return lProtectionLevel;
}
// ////////////////////////////////////////////////////////////////////
- const double EmsrUtils::computeEmsrValue (double iCapacity, Bucket& ioBucket) {
- // Retrive the average yield, mean and standard deviation of the
+ const double EmsrUtils::computeEmsrValue (double iCapacity,
+ Bucket& ioBucket) {
+ // Retrieve the average yield, mean and standard deviation of the
// demand of the class/bucket.
const double lMean = ioBucket.getMean();
const double lSD = ioBucket.getStandardDeviation();
const double lYield = ioBucket.getAverageYield();
// Compute the EMSR value = lYield * Pr (demand >= iCapacity).
- const double emsrValue = lYield * gsl_cdf_gaussian_Q(iCapacity-lMean, lSD);
+ boost::math::normal lNormalDistribution (lMean, lSD);
+ const double emsrValue =
+ lYield * boost::math::cdf (boost::math::complement (lNormalDistribution,
+ iCapacity));
return emsrValue;
}
Modified: trunk/rmol/rmol/bom/ExpectationMaximization.cpp
===================================================================
--- trunk/rmol/rmol/bom/ExpectationMaximization.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/ExpectationMaximization.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -1,11 +1,11 @@
// //////////////////////////////////////////////////////////////////////
// Import section
// //////////////////////////////////////////////////////////////////////
-// GSL
-#include <gsl/gsl_cdf.h>
// STL
#include <cassert>
#include <cmath>
+// Boost Math
+#include <boost/math/distributions/normal.hpp>
// StdAir
#include <stdair/service/Logger.hpp>
// RMOL
@@ -125,18 +125,21 @@
*/
double e, d1, d2;
const double lerror = kthCensoredData - iMean;
- d1 = gsl_cdf_gaussian_Q (lerror, iSD);
- e = - (lerror) * (lerror) * 0.5 / (iSD * iSD);
- d2 = exp(e) * iSD / sqrt(2 * 3.14159265);
+
+ //
+ boost::math::normal lNormalDistribution (iMean, iSD);
+ d1 = boost::math::cdf (boost::math::complement (lNormalDistribution,
+ kthCensoredData));
+ e = -lerror*lerror * 0.5 / (iSD*iSD);
+ d2 = exp(e) * iSD / sqrt (2 * 3.14159265);
if (d1 < DEFAULT_EPSILON) {
- ioUnconstrainedDataHolder.push_back(kthCensoredData);
+ ioUnconstrainedDataHolder.push_back (kthCensoredData);
+
+ } else {
+ ioUnconstrainedDataHolder.push_back (iMean + d2/d1);
}
- else {
- ioUnconstrainedDataHolder.push_back(iMean + d2/d1);
- }
}
}
- // //////////////////////////////////////////////////////////////////////
}
Modified: trunk/rmol/rmol/bom/Gaussian.cpp
===================================================================
--- trunk/rmol/rmol/bom/Gaussian.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/Gaussian.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -2,76 +2,46 @@
// //////////////////////////////////////////////////////////////////////
// Import section
// //////////////////////////////////////////////////////////////////////
-// GSL Random Number Distributions (GSL Reference Manual, version 1.7,
-// Chapter 19)
-#include <gsl/gsl_randist.h>
// RMOL
#include <rmol/bom/Gaussian.hpp>
namespace RMOL {
// //////////////////////////////////////////////////////////////////////
- Gaussian::Gaussian () :
- _rngTypePtr (gsl_rng_default), _rngPtr (NULL),
- _params (FldDistributionParameters()) {
- initRandomGenerator();
+ Gaussian::Gaussian (const FldDistributionParameters& iParams)
+ : _seed (42u), _generator (42u),
+ _normalDistribution (iParams.getMean(), iParams.getStandardDeviation()),
+ _normalDistributionGenerator (_generator, _normalDistribution) {
+ init();
}
// //////////////////////////////////////////////////////////////////////
- Gaussian::Gaussian (const Gaussian& iGaussian) :
- _rngTypePtr (gsl_rng_default), _rngPtr (NULL),
- _params (iGaussian.getDistributionParameters()) {
- initRandomGenerator();
+ void Gaussian::init() {
}
// //////////////////////////////////////////////////////////////////////
- Gaussian::Gaussian (const FldDistributionParameters& iParams) :
- _rngTypePtr (gsl_rng_default), _rngPtr (NULL),
- _params (FldDistributionParameters (iParams.getMean(),
- iParams.getStandardDeviation())) {
- initRandomGenerator();
- }
-
- // //////////////////////////////////////////////////////////////////////
Gaussian::~Gaussian() {
// Release the memory for the random generator
- gsl_rng_free (_rngPtr);
}
// //////////////////////////////////////////////////////////////////////
double Gaussian::getMean() const {
- return _params.getMean();
+ return _normalDistribution.mean();
}
// //////////////////////////////////////////////////////////////////////
double Gaussian::getStandardDeviation() const {
- return _params.getStandardDeviation();
+ return _normalDistribution.sigma();
}
-
// //////////////////////////////////////////////////////////////////////
double Gaussian::getVariance() const {
- return _params.getVariance();
+ return (_normalDistribution.sigma() * _normalDistribution.sigma());
}
// //////////////////////////////////////////////////////////////////////
- void Gaussian::initRandomGenerator () {
- // Initialise the Random Generator
- gsl_rng_env_setup ();
-
- // Allocate the memory for the random generator
- _rngPtr = gsl_rng_alloc (_rngTypePtr);
- }
-
- // //////////////////////////////////////////////////////////////////////
- double Gaussian::generateVariate () const {
-
- const double mean = getMean();
- const double standardDeviation = getStandardDeviation();
-
- double result = gsl_ran_gaussian (_rngPtr, standardDeviation);
- result += mean;
-
+ double Gaussian::generateVariate() {
+ const double result = _normalDistributionGenerator();
return result;
}
Modified: trunk/rmol/rmol/bom/Gaussian.hpp
===================================================================
--- trunk/rmol/rmol/bom/Gaussian.hpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/Gaussian.hpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -1,36 +1,61 @@
-#ifndef __RMOL_GAUSSIAN_HPP
-#define __RMOL_GAUSSIAN_HPP
+#ifndef __RMOL_BOM_GAUSSIAN_HPP
+#define __RMOL_BOM_GAUSSIAN_HPP
// //////////////////////////////////////////////////////////////////////
// Import section
// //////////////////////////////////////////////////////////////////////
-// GSL Random Number Generation (GSL Reference Manual, version 1.7, Chapter 17)
-#include <gsl/gsl_rng.h>
+// Boost Random
+#include <boost/random/linear_congruential.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/variate_generator.hpp>
// RMOL
#include <rmol/field/FldDistributionParameters.hpp>
namespace RMOL {
- /** Gaussian Distribution-based Utilities. */
+ // ////////// Type definitions //////////
+ /** Type definition for the random generator seed.
+ <br>That seed must be unsigned, otherwise the wrong overload may be
+ selected when using mt19937 as the base_generator_type. */
+ typedef unsigned int random_generator_seed_type;
+
+
+ /** Wrapper around a random generator following a normal distribution. */
class Gaussian {
+ private:
+ // ////////// Type definitions //////////
+ /** Type definition for a random number generator base (mt19937). */
+ typedef boost::minstd_rand base_generator_type;
+
+ /** Type definiton for the normal distribution (characteristics). */
+ typedef boost::normal_distribution<> normal_dist_type;
+
+ /** Type definition for the normal distribution random generator. */
+ typedef boost::variate_generator<base_generator_type&,
+ normal_dist_type> normal_gen_type;
+
public:
- /** Constructors. */
+ /** Constructor with mean and sigma (standard deviation) for
+ the normal Distribution.
+ <br>See also
+ http://www.boost.org/doc/libs/1_44_0/doc/html/boost/normal_distribution.html */
+ Gaussian (const FldDistributionParameters&);
+
+ /** Destructor. */
+ ~Gaussian();
+
+ private:
+ /** Default constructors.
+ <br>They are kept private so that the class can be instantiated only
+ with the public constructor. */
Gaussian ();
Gaussian (const Gaussian&);
- /** Constructor with mean and standard deviation of
- the Gaussian Distribution. */
- Gaussian (const FldDistributionParameters&);
- /** Destructors. */
- virtual ~Gaussian();
-
- // Getters
- /** Getter for the parameters for the Gaussian distribution (i.e.,
- mean and standard deviation). */
- const FldDistributionParameters& getDistributionParameters() const {
- return _params;
- }
+ /** Initialise the random generator. */
+ void init();
+ public:
+ // /////////////////// Getters //////////////////
/** Getter for the mean value. */
double getMean() const;
@@ -43,20 +68,21 @@
/** Generate a Gaussian random variate (following the Gaussian
distribution). */
- double generateVariate () const;
+ double generateVariate();
private:
- // Wrapper on GSL Random Generators (type and generator)
- const gsl_rng_type* _rngTypePtr;
- gsl_rng* _rngPtr;
+ /** Seed of the random generator. */
+ const random_generator_seed_type _seed;
- // Gaussian distribution characteristics
- FldDistributionParameters _params;
+ /** Base for the random generator. */
+ base_generator_type _generator;
- private:
- // Initialise the Random Generator
- void initRandomGenerator();
+ /** Normal distribution. */
+ normal_dist_type _normalDistribution;
+
+ /** Random generator for the normal distribution. */
+ normal_gen_type _normalDistributionGenerator;
};
}
-#endif // __RMOL_GAUSSIAN_HPP
+#endif // __RMOL_BOM_GAUSSIAN_HPP
Modified: trunk/rmol/rmol/bom/HistoricalBookingHolder.cpp
===================================================================
--- trunk/rmol/rmol/bom/HistoricalBookingHolder.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/HistoricalBookingHolder.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -1,18 +1,14 @@
// //////////////////////////////////////////////////////////////////////
// Import section
// //////////////////////////////////////////////////////////////////////
-// GSL Random Number Generation
-// (GSL Reference Manual, version 1.7, Chapter 19)
-#include <gsl/gsl_cdf.h>
-#include <gsl/gsl_randist.h>
-// C
-#include <assert.h>
-#include <math.h>
// STL
+#include <cassert>
#include <iostream>
#include <iomanip>
#include <cmath>
-// RMU
+// Boost Math
+#include <boost/math/distributions/normal.hpp>
+// RMOL
#include <rmol/bom/HistoricalBooking.hpp>
#include <rmol/bom/HistoricalBookingHolder.hpp>
@@ -211,13 +207,15 @@
d1 = 0.5 * (1 + s);
}
*/
+
+ //
+ boost::math::normal lNormalDistribution (iMean, iSD);
+ d1 = boost::math::cdf (boost::math::complement (lBooking, _capacity));
- d1 = gsl_cdf_gaussian_Q (lBooking - iMean, iSD);
-
e = - (lBooking - iMean) * (lBooking - iMean) * 0.5 / (iSD * iSD);
e = exp (e);
- d2 = e * iSD / sqrt(2 * 3.14159265);
+ d2 = e * iSD / sqrt (2 * 3.14159265);
// std::cout << "d1, d2 = " << d1 << " " << d2 << std::endl;
Modified: trunk/rmol/rmol/bom/HistoricalBookingHolderHolder.cpp
===================================================================
--- trunk/rmol/rmol/bom/HistoricalBookingHolderHolder.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/HistoricalBookingHolderHolder.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -1,17 +1,14 @@
// //////////////////////////////////////////////////////////////////////
// Import section
// //////////////////////////////////////////////////////////////////////
-// GSL Random Number Generation (GSL Reference Manual, version 1.7, Chapter 19)
-#include <gsl/gsl_cdf.h>
-#include <gsl/gsl_randist.h>
-// C
-#include <assert.h>
-#include <math.h>
// STL
+#include <cassert>
#include <iostream>
#include <iomanip>
#include <cmath>
-// RMU
+// Boost Math
+#include <boost/math/distributions/normal.hpp>
+// RMOL
#include <rmol/bom/HistoricalBookingHolder.hpp>
#include <rmol/bom/HistoricalBookingHolderHolder.hpp>
@@ -210,9 +207,9 @@
// d1 = 0.5 * (1 + s);
// }
// */
-
-// d1 = gsl_cdf_gaussian_Q (lBooking - iMean, iSD);
-
+//
+// boost::math::normal lNormalDistribution (iMean, iSD);
+// d1 = boost::math::cdf (boost::math::complement (lBooking, _capacity));
// e = - (lBooking - iMean) * (lBooking - iMean) * 0.5 / (iSD * iSD);
// e = exp (e);
Modified: trunk/rmol/rmol/bom/HistoricalDataHolder.cpp
===================================================================
--- trunk/rmol/rmol/bom/HistoricalDataHolder.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/HistoricalDataHolder.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -1,17 +1,11 @@
// //////////////////////////////////////////////////////////////////////
// Import section
// //////////////////////////////////////////////////////////////////////
-// GSL Random Number Generation
-// (GSL Reference Manual, version 1.7, Chapter 19)
-#include <gsl/gsl_cdf.h>
-#include <gsl/gsl_randist.h>
-// C
-#include <assert.h>
-#include <math.h>
// STL
+#include <cassert>
#include <iostream>
//#include <iomanip>
-//#include <cmath>
+#include <cmath>
// RMU
#include <rmol/bom/HistoricalDataHolder.hpp>
Modified: trunk/rmol/rmol/bom/HistoricalDataHolderHolder.cpp
===================================================================
--- trunk/rmol/rmol/bom/HistoricalDataHolderHolder.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/HistoricalDataHolderHolder.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -1,18 +1,12 @@
// //////////////////////////////////////////////////////////////////////
// Import section
// //////////////////////////////////////////////////////////////////////
-// GSL Random Number Generation
-// (GSL Reference Manual, version 1.7, Chapter 19)
-// #include <gsl/gsl_cdf.h>
-// #include <gsl/gsl_randist.h>
-// C
-// #include <assert.h>
-// #include <math.h>
// STL
+ #include <cassert>
#include <iostream>
-#include <iomanip>
-// #include <cmath>
-// RMU
+//#include <iomanip>
+//#include <cmath>
+// RMOL
#include <rmol/bom/HistoricalDataHolder.hpp>
#include <rmol/bom/HistoricalDataHolderHolder.hpp>
Modified: trunk/rmol/rmol/bom/MCOptimiser.cpp
===================================================================
--- trunk/rmol/rmol/bom/MCOptimiser.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/MCOptimiser.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -4,9 +4,7 @@
// STL
#include <cassert>
#include <string>
-#include <fstream>
#include <sstream>
-#include <cmath>
// StdAir
#include <stdair/basic/BasChronometer.hpp>
#include <stdair/service/Logger.hpp>
@@ -75,7 +73,7 @@
*/
const FldDistributionParameters& aDistribParams =
currentBucket.getDistributionParameters();
- const Gaussian gaussianDemandGenerator (aDistribParams);
+ Gaussian gaussianDemandGenerator (aDistribParams);
/** DEBUG
STDAIR_LOG_DEBUG ("[" << j << "]: " << Kj << " values with N ( "
@@ -208,7 +206,7 @@
// STEP 1.
const FldDistributionParameters& aDistribParams =
currentBucket.getDistributionParameters();
- const Gaussian gaussianDemandGenerator (aDistribParams);
+ Gaussian gaussianDemandGenerator (aDistribParams);
VariateList_T aVariateList;
@@ -308,7 +306,7 @@
*/
const FldDistributionParameters& aDistribParams =
currentBucket.getDistributionParameters();
- const Gaussian gaussianDemandGenerator (aDistribParams);
+ Gaussian gaussianDemandGenerator (aDistribParams);
/** DEBUG
STDAIR_LOG_DEBUG ("[" << j << "]: " << Kj << " values with N ( "
@@ -452,7 +450,7 @@
// STEP 1.
const FldDistributionParameters& aDistribParams =
currentBucket.getDistributionParameters();
- const Gaussian gaussianDemandGenerator (aDistribParams);
+ Gaussian gaussianDemandGenerator (aDistribParams);
VariateList_T aVariateList;
Modified: trunk/rmol/rmol/bom/Makefile.am
===================================================================
--- trunk/rmol/rmol/bom/Makefile.am 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/Makefile.am 2010-09-27 11:36:24 UTC (rev 338)
@@ -6,8 +6,8 @@
noinst_LTLIBRARIES = librmolbom.la
librmolbom_la_SOURCES = $(rmol_bom_h_sources) $(rmol_bom_cc_sources)
-librmolbom_la_CXXFLAGS = $(GSL_CFLAGS) $(BOOST_CFLAGS) $(STDAIR_CFLAGS)
-librmolbom_la_LDFLAGS = $(GSL_LIBS)
+librmolbom_la_CXXFLAGS = $(BOOST_CFLAGS) $(STDAIR_CFLAGS)
+librmolbom_la_LDFLAGS =
#
#pkgincludedir = $(includedir)/@PACKAGE@/bom
Modified: trunk/rmol/rmol/bom/Overbooking.cpp
===================================================================
--- trunk/rmol/rmol/bom/Overbooking.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/bom/Overbooking.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -1,16 +1,14 @@
// //////////////////////////////////////////////////////////////////////
// Import section
// //////////////////////////////////////////////////////////////////////
-// C
-#include <assert.h>
// STL
-#include <iostream>
-// GSL
-#include <gsl/gsl_sys.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_sf.h>
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_cdf.h>
+#include <cassert>
+// Boost Math
+#include <boost/math/special_functions/pow.hpp>
+#include <boost/math/special_functions/powm1.hpp>
+#include <boost/math/special_functions/round.hpp>
+#include <boost/math/distributions/binomial.hpp>
+#include <boost/math/distributions/normal.hpp>
// RMOL
#include <rmol/basic/BasConst_Overbooking.hpp>
#include <rmol/bom/Overbooking.hpp>
@@ -335,9 +333,13 @@
const double lDemandStdDev =
_demandDistributionParameters.getStandardDeviation();
+ boost::math::normal lNormalDistribution (lDemandMean,
+ lDemandStdDev);
+
// Algorithm
double pNormal =
- gsl_cdf_gaussian_Q (_capacity - lDemandMean, lDemandStdDev);
+ boost::math::cdf (boost::math::complement (lNormalDistribution,
+ _capacity));
// double pNormal = probabilityNormal (CAPACITY, DEMAND_AVERAGE,
// DEMAND_DEVIATION);
@@ -348,20 +350,22 @@
}
assert (PRICE_OVER_DENIED_COST < 1 - pNormal);
- pNormal = gsl_cdf_gaussian_Q (_capacity - 1 - lDemandMean,
- lDemandStdDev);
+ pNormal =
+ boost::math::cdf (boost::math::complement (lNormalDistribution,
+ _capacity - 1));
// pNormal = probabilityNormal (CAPACITY-1, DEMAND_AVERAGE,
// DEMAND_DEVIATION);
const double lNoShowMean = _noShowDistributionParameters.getMean();
-
- double lProbability = (1 - pNormal)
- * gsl_sf_pow_int (lNoShowMean, _capacity);
+ double lProbability =
+ (1 - pNormal) * (boost::math::powm1 (lNoShowMean, _capacity) + 1);
+
int counter = 1;
- while ((lProbability < PRICE_OVER_DENIED_COST) && (counter < 100)) {
- counter++;
+ while (lProbability < PRICE_OVER_DENIED_COST && counter < 100) {
+ ++counter;
+
if (resultBookingLimit >= MAX_BOOKING_LIMIT) {
return resultBookingLimit;
@@ -371,7 +375,9 @@
const unsigned int b = static_cast<unsigned int> (resultBookingLimit);
- pNormal = gsl_cdf_gaussian_Q (b + 1 - lDemandMean, lDemandStdDev);
+ pNormal =
+ boost::math::cdf (boost::math::complement (lNormalDistribution,
+ b + 1));
// pNormal = probabilityNormal (b+1, DEMAND_AVERAGE, DEMAND_DEVIATION);
@@ -380,51 +386,30 @@
// assert (b >= 0); (derived from the other two)
- // The cumulated probability that the number of shows exceeds
- // the leg/cabin physical capacity
- lProbability += (1+pNormal)
- * lNoShowMean * gsl_ran_binomial_pdf (_capacity-1, lNoShowMean, b);
-
+ /**
+ The probability that the number of shows exceeds the leg/cabin
+ physical capacity follows a binomial distribution, with the following
+ parameters:<br>
+ <ul>
+ <li>"Probability of success (== no-show)": p == lNoShowMean</li>
+ <li>"Number of trials (== bookings)": n == b</li>
+ <li>"Number of successes (no-shows)": k == _capacity - 1</li>
+ </ul> */
+ boost::math::binomial lBinomialDistribution (lNoShowMean, b);
+ lProbability += (1 + pNormal) * lNoShowMean
+ * boost::math::pdf (lBinomialDistribution, _capacity-1);
}
return resultBookingLimit;
}
- // Private useful functions, some of them are to be replaced by gsl
- // library's functions
+ // Private useful functions
// //////////////////////////////////////////////////////////////////////
- //to be deleted, already available in gsl lib
- //cumulated binomial probability
- double Overbooking::probabilityNormal (const int b, const double average,
- const double sdeviation) const {
- double resultCdf = 0.0;
- for (int i=0; i < b+1 ; i++) {
- resultCdf += 1 /( sqrt(2 * M_PI) )
- * gsl_sf_exp (-1* (gsl_sf_pow_int(i-average, 2)) / (2*sdeviation));
- }
- return resultCdf;
- }
-
- // //////////////////////////////////////////////////////////////////////
-
- // to be deleted, already available in gsl lib, gsl_ran_binomial_pdf(unsigned int k, double SHOW_RATE, unsigned int b)
- double Overbooking::binomialProbability_F_b_s (const double iShowRate,
- const int b,
- const int k) const {
- double Factorials = gsl_sf_fact(b)/gsl_sf_fact(k)/gsl_sf_fact(b-k);
- double f = Factorials*
- gsl_sf_pow_int(iShowRate, k)*gsl_sf_pow_int(1-iShowRate, b-k) ;
-
- return f;
- }
-
// pSHOW_RATEbability to deny one or more services
// //////////////////////////////////////////////////////////////////////
-
-
////////////// Private functions for Servive level policies/////////////
// service level 1 : cumulated probability of having no-shows
@@ -432,10 +417,21 @@
// //////////////////////////////////////////////////////////////////////
double Overbooking::serviceLevel1 (const double iShowRate,
- const int b, const int C) const {
- double resultProbability = 1.0;
- for (int i = 1; i<=C; i++) {
- resultProbability -= gsl_ran_binomial_pdf (i, iShowRate, b);}
+ const int b, const int iCapacity) const {
+
+ /**
+ The probability that the number of shows exceeds the leg/cabin
+ physical capacity follows a binomial distribution, with the following
+ parameters:<br>
+ <ul>
+ <li>"Probability of success (== no-show)": p == iShowRate</li>
+ <li>"Number of trials (== bookings)": n == b</li>
+ <li>"Number of successes (no-shows)": k == iCapacity</li>
+ </ul> */
+ boost::math::binomial lBinomialDistribution (iShowRate, b);
+ const double resultProbability =
+ boost::math::cdf (lBinomialDistribution, iCapacity);
+
return resultProbability;
}
Modified: trunk/rmol/rmol/core/Makefile.am
===================================================================
--- trunk/rmol/rmol/core/Makefile.am 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/core/Makefile.am 2010-09-27 11:36:24 UTC (rev 338)
@@ -22,6 +22,6 @@
$(top_builddir)/rmol/service/librmolsvc.la
librmol_la_LDFLAGS = \
$(BOOST_DATE_TIME_LIB) $(BOOST_PROGRAM_OPTIONS_LIB) \
- $(BOOST_FILESYSTEM_LIB) $(STDAIR_LIBS) $(GSL_LIBS) \
+ $(BOOST_FILESYSTEM_LIB) $(STDAIR_LIBS) \
-version-info $(GENERIC_LIBRARY_VERSION)
Modified: trunk/rmol/rmol/service/RMOL_ServiceContext.cpp
===================================================================
--- trunk/rmol/rmol/service/RMOL_ServiceContext.cpp 2010-09-13 18:56:50 UTC (rev 337)
+++ trunk/rmol/rmol/service/RMOL_ServiceContext.cpp 2010-09-27 11:36:24 UTC (rev 338)
@@ -120,11 +120,11 @@
lDemandVector.reserve (K);
const FldDistributionParameters aDistributionParam =
FldDistributionParameters (iMean, iDeviation);
- const Gaussian gaussianDemandGenerator (aDistributionParam);
+ Gaussian gaussianDemandGenerator (aDistributionParam);
// Generate K numbers
for (int i = 0; i < K; ++i) {
- const double lGeneratedDemand = gaussianDemandGenerator.generateVariate ();
+ const double lGeneratedDemand= gaussianDemandGenerator.generateVariate ();
lDemandVector.push_back (lGeneratedDemand);
}
@@ -137,17 +137,26 @@
GeneratedDemandVector_T* ioSecondVector) {
if (ioFirstVector == NULL || ioSecondVector == NULL) {
return NULL;
+
} else {
+ assert (ioFirstVector != NULL);
+ assert (ioSecondVector != NULL);
+
const unsigned int K = ioFirstVector->size();
assert (K == ioSecondVector->size());
+
_generatedDemandVectorHolder.push_back (DEFAULT_GENERATED_DEMAND_VECTOR);
GeneratedDemandVectorHolder_T::reverse_iterator itLastVector =
_generatedDemandVectorHolder.rbegin();
+
GeneratedDemandVector_T& lDemandVector = *itLastVector;
lDemandVector.reserve (K);
- for (unsigned int i = 0; i < K; ++i) {
- const double lGeneratedDemand =
- ioFirstVector->at(i) + ioSecondVector->at(i);
+ GeneratedDemandVector_T::const_iterator itFirst = ioFirstVector->begin();
+ GeneratedDemandVector_T::const_iterator itSecond= ioSecondVector->begin();
+ for ( ; itFirst != ioFirstVector->end(); ++itFirst, ++itSecond) {
+ const double& lFirst = *itFirst;
+ const double& lSecond = *itSecond;
+ const double lGeneratedDemand = lFirst + lSecond;
lDemandVector.push_back (lGeneratedDemand);
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|