timeseries.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002   file : $URL: https://frepple.svn.sourceforge.net/svnroot/frepple/trunk/modules/forecast/forecastsolver.cpp $
00003   version : $LastChangedRevision: 1108 $  $LastChangedBy: jdetaeye $
00004   date : $LastChangedDate: 2009-12-06 18:54:18 +0100 (Sun, 06 Dec 2009) $
00005  ***************************************************************************/
00006 
00007 /***************************************************************************
00008  *                                                                         *
00009  * Copyright (C) 2007 by Johan De Taeye                                    *
00010  *                                                                         *
00011  * This library is free software; you can redistribute it and/or modify it *
00012  * under the terms of the GNU Lesser General Public License as published   *
00013  * by the Free Software Foundation; either version 2.1 of the License, or  *
00014  * (at your option) any later version.                                     *
00015  *                                                                         *
00016  * This library is distributed in the hope that it will be useful,         *
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of          *
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser *
00019  * General Public License for more details.                                *
00020  *                                                                         *
00021  * You should have received a copy of the GNU Lesser General Public        *
00022  * License along with this library; if not, write to the Free Software     *
00023  * Foundation Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,*
00024  * USA                                                                     *
00025  *                                                                         *
00026  ***************************************************************************/
00027 
00028 #include "forecast.h"
00029 
00030 namespace module_forecast
00031 {
00032 
00033 #define ACCURACY 0.01
00034 
00035 void Forecast::generateFutureValues(
00036   const double history[], unsigned int historycount,
00037   const Date buckets[], unsigned int bucketcount,
00038   bool debug)
00039 {
00040   // Validate the input
00041   if (!history || !buckets)
00042     throw RuntimeException("Null argument to forecast function");
00043   if (bucketcount < 2)
00044     throw DataException("Need at least 2 forecast dates");
00045 
00046   // Create a list of forecasting methods.
00047   // We create the forecasting objects in stack memory for best performance.
00048   MovingAverage moving_avg;
00049   SingleExponential single_exp;
00050   DoubleExponential double_exp;
00051   int numberOfMethods = 3;
00052   ForecastMethod* methods[3];
00053   methods[0] = &single_exp;
00054   methods[1] = &double_exp;
00055   methods[2] = &moving_avg;
00056 
00057   // Initialize a vector with the mad weights
00058   double *madWeight = new double[historycount+1];
00059   madWeight[historycount] = madWeight[historycount-1] = 1.0;
00060   for (int i=historycount-2; i>=0; --i)
00061     madWeight[i] = madWeight[i+1] * Forecast::getForecastMadAlfa();
00062 
00063   // Evaluate each forecast method
00064   double best_error = DBL_MAX;
00065   int best_method = -1;
00066   double error;
00067   try
00068   {
00069     for (int i=0; i<numberOfMethods; ++i)
00070     {
00071       error = methods[i]->generateForecast(this, history, historycount, madWeight, debug);
00072       if (error<best_error)
00073       {
00074         best_error = error;
00075         best_method = i;
00076       }
00077     }
00078   }
00079   catch (...)
00080   {
00081     delete madWeight;
00082     throw;
00083   }
00084   delete madWeight;
00085 
00086   // Apply the most appropriate forecasting method
00087   if (best_method >= 0)
00088   {
00089     if (debug)
00090       logger << getName() << ": chosen method : " << methods[best_method]->getName() << endl;
00091     methods[best_method]->applyForecast(this, buckets, bucketcount, debug);
00092   }
00093 }
00094 
00095 
00096 //
00097 // MOVING AVERAGE FORECAST
00098 //
00099 
00100 
00101 unsigned int Forecast::MovingAverage::defaultbuckets = 5;
00102 
00103 
00104 double Forecast::MovingAverage::generateForecast
00105   (Forecast* fcst, const double history[], unsigned int count, const double madWeight[], bool debug)
00106 {
00107   double error_mad = 0.0;
00108   double sum = 0.0;
00109 
00110   // Calculate the forecast and forecast error.
00111   for (unsigned int i = 1; i <= count; ++i)
00112   {
00113     sum += history[i-1];
00114     if (i>buckets)
00115     {
00116       sum -= history[i-1-buckets];
00117       avg = sum / buckets;
00118     }
00119     else
00120       avg = sum / i;
00121     if (i >= fcst->getForecastSkip())
00122       error_mad += fabs(avg - history[i]) * madWeight[i];
00123   }
00124 
00125   // Echo the result
00126   if (debug)
00127     logger << (fcst ? fcst->getName() : "") << ": moving average : "
00128       << "mad " << error_mad
00129       << ", forecast " << avg << endl;
00130   return error_mad;
00131 }
00132 
00133 
00134 void Forecast::MovingAverage::applyForecast
00135   (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00136 {
00137   // Loop over all buckets and set the forecast to a constant value
00138   if (avg < 0) return;
00139   for (unsigned int i = 1; i < bucketcount; ++i)
00140     forecast->setTotalQuantity(
00141       DateRange(buckets[i-1], buckets[i]),
00142       avg
00143       );
00144 }
00145 
00146 
00147 //
00148 // SINGLE EXPONENTIAL FORECAST
00149 //
00150 
00151 
00152 double Forecast::SingleExponential::initial_alfa = 0.2;
00153 double Forecast::SingleExponential::min_alfa = 0.03;
00154 double Forecast::SingleExponential::max_alfa = 1.0;
00155 
00156 
00157 double Forecast::SingleExponential::generateForecast
00158   (Forecast* fcst, const double history[], unsigned int count, const double madWeight[], bool debug)
00159 {
00160   // Verify whether this is a valid forecast method.
00161   //   - We need at least 5 buckets after the warmup period.
00162   if (count < fcst->getForecastSkip() + 5)
00163     return DBL_MAX;
00164 
00165   unsigned int iteration = 1;
00166   bool upperboundarytested = false;
00167   bool lowerboundarytested = false;
00168   double error_mad = 0.0, delta, df_dalfa_i, sum_11, sum_12;
00169   double best_error_mad = DBL_MAX, best_alfa = initial_alfa, best_f_i = 0.0;
00170   for (; iteration <= Forecast::getForecastIterations(); ++iteration)
00171   {
00172     // Initialize the iteration
00173     f_i = history[0];
00174     df_dalfa_i = sum_11 = sum_12 = error_mad = 0.0;
00175 
00176     // Calculate the forecast and forecast error.
00177     // We also compute the sums required for the Marquardt optimization.
00178     for (unsigned long i = 1; i <= count; ++i)
00179     {
00180       df_dalfa_i = history[i-1] - f_i + (1 - alfa) * df_dalfa_i;
00181       f_i = history[i-1] * alfa + (1 - alfa) * f_i;
00182       if (i == count) break;
00183       sum_12 += df_dalfa_i * (history[i] - f_i) * madWeight[i];
00184       sum_11 += df_dalfa_i * df_dalfa_i * madWeight[i];
00185       if (i >= fcst->getForecastSkip())
00186         error_mad += fabs(f_i - history[i]) * madWeight[i];
00187     }
00188 
00189     // Better than earlier iterations?
00190     if (error_mad < best_error_mad)
00191     {
00192       best_error_mad = error_mad;
00193       best_alfa = alfa;
00194       best_f_i = f_i;
00195     }
00196 
00197     // Add Levenberg - Marquardt damping factor
00198     if (fabs(sum_11 + error_mad / iteration) > ROUNDING_ERROR)
00199       sum_11 += error_mad / iteration;
00200 
00201     // Calculate a delta for the alfa parameter
00202     if (fabs(sum_11) < ROUNDING_ERROR) break;
00203     delta = sum_12 / sum_11;
00204 
00205     // Stop when we are close enough and have tried hard enough
00206     if (fabs(delta) < ACCURACY && iteration > 3) break;
00207 
00208     // New alfa
00209     alfa += delta;
00210 
00211     // Stop when we repeatedly bounce against the limits.
00212     // Testing a limits once is allowed.
00213     if (alfa > max_alfa)
00214     {
00215       alfa = max_alfa;
00216       if (upperboundarytested) break;
00217       upperboundarytested = true;
00218     }
00219     else if (alfa < min_alfa)
00220     {
00221       alfa = min_alfa;
00222       if (lowerboundarytested) break;
00223       lowerboundarytested = true;
00224     }
00225   }
00226 
00227   // Keep the best result
00228   f_i = best_f_i;
00229 
00230   // Echo the result
00231   if (debug)
00232     logger << (fcst ? fcst->getName() : "") << ": single exponential : "
00233       << "alfa " << best_alfa
00234       << ", mad " << best_error_mad
00235       << ", " << iteration << " iterations"
00236       << ", forecast " << f_i << endl;
00237   return best_error_mad;
00238 }
00239 
00240 
00241 void Forecast::SingleExponential::applyForecast
00242   (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00243 {
00244   // Loop over all buckets and set the forecast to a constant value
00245   if (f_i < 0) return;
00246   for (unsigned int i = 1; i < bucketcount; ++i)
00247     forecast->setTotalQuantity(
00248       DateRange(buckets[i-1], buckets[i]),
00249       f_i
00250       );
00251 }
00252 
00253 
00254 //
00255 // DOUBLE EXPONENTIAL FORECAST
00256 //
00257 
00258 
00259 double Forecast::DoubleExponential::initial_alfa = 0.2;
00260 double Forecast::DoubleExponential::min_alfa = 0.02;
00261 double Forecast::DoubleExponential::max_alfa = 1.0;
00262 double Forecast::DoubleExponential::initial_gamma = 0.2;
00263 double Forecast::DoubleExponential::min_gamma = 0.05;
00264 double Forecast::DoubleExponential::max_gamma = 1.0;
00265 
00266 
00267 double Forecast::DoubleExponential::generateForecast
00268   (Forecast* fcst, const double history[], unsigned int count, const double madWeight[], bool debug)
00269 {
00270   // Verify whether this is a valid forecast method.
00271   //   - We need at least 5 buckets after the warmup period.
00272   if (count < fcst->getForecastSkip() + 5)
00273     return DBL_MAX;
00274 
00275   // Define variables
00276   double error_mad = 0.0, delta_alfa, delta_gamma, determinant;
00277   double constant_i_prev, trend_i_prev, d_constant_d_gamma_prev,
00278     d_constant_d_alfa_prev, d_constant_d_alfa, d_constant_d_gamma,  
00279     d_trend_d_alfa, d_trend_d_gamma, d_forecast_d_alfa, d_forecast_d_gamma,
00280     sum11, sum12, sum22, sum13, sum23;
00281   double best_error_mad = DBL_MAX, best_alfa = initial_alfa,
00282     best_gamma = initial_gamma, best_constant_i = 0.0, best_trend_i = 0.0;
00283 
00284   // Iterations
00285   unsigned int iteration = 1, boundarytested = 0;
00286   for (; iteration <= Forecast::getForecastIterations(); ++iteration)
00287   {
00288     // Initialize the iteration
00289     constant_i = history[0];
00290     trend_i = history[1] - history[0];
00291     error_mad = sum11 = sum12 = sum22 = sum13 = sum23 = 0.0;
00292     d_constant_d_alfa = d_constant_d_gamma = d_trend_d_alfa = d_trend_d_gamma = 0.0;
00293     d_forecast_d_alfa = d_forecast_d_gamma = 0.0;
00294 
00295     // Calculate the forecast and forecast error.
00296     // We also compute the sums required for the Marquardt optimization.
00297     for (unsigned long i = 1; i <= count; ++i)
00298     {
00299       constant_i_prev = constant_i;
00300       trend_i_prev = trend_i;
00301       constant_i = history[i-1] * alfa + (1 - alfa) * (constant_i_prev + trend_i_prev);
00302       trend_i = gamma * (constant_i - constant_i_prev) + (1 - gamma) * trend_i_prev;
00303       if (i == count) break;
00304       d_constant_d_gamma_prev = d_constant_d_gamma;
00305       d_constant_d_alfa_prev = d_constant_d_alfa;
00306       d_constant_d_alfa = history[i-1] - constant_i_prev - trend_i_prev
00307         + (1 - alfa) * d_forecast_d_alfa;
00308       d_constant_d_gamma = (1 - alfa) * d_forecast_d_gamma;
00309       d_trend_d_alfa = gamma * (d_constant_d_alfa - d_constant_d_alfa_prev)
00310         + (1 - gamma) * d_trend_d_alfa;
00311       d_trend_d_gamma = constant_i - constant_i_prev - trend_i_prev
00312         + gamma * (d_constant_d_gamma - d_constant_d_gamma_prev)
00313         + (1 - gamma) * d_trend_d_gamma;
00314       d_forecast_d_alfa = d_constant_d_alfa + d_trend_d_alfa;
00315       d_forecast_d_gamma = d_constant_d_gamma + d_trend_d_gamma;
00316       sum11 += madWeight[i] * d_forecast_d_alfa * d_forecast_d_alfa;
00317       sum12 += madWeight[i] * d_forecast_d_alfa * d_forecast_d_gamma;
00318       sum22 += madWeight[i] * d_forecast_d_gamma * d_forecast_d_gamma;
00319       sum13 += madWeight[i] * d_forecast_d_alfa * (history[i] - constant_i - trend_i);
00320       sum23 += madWeight[i] * d_forecast_d_gamma * (history[i] - constant_i - trend_i);
00321       if (i >= fcst->getForecastSkip()) // Don't measure during the warmup period
00322         error_mad += fabs(constant_i + trend_i - history[i]) * madWeight[i];
00323     }
00324 
00325     // Better than earlier iterations?
00326     if (error_mad < best_error_mad)
00327     {
00328       best_error_mad = error_mad;
00329       best_alfa = alfa;
00330       best_gamma = gamma;
00331       best_constant_i = constant_i;
00332       best_trend_i = trend_i;
00333     }
00334 
00335     // Add Levenberg - Marquardt damping factor
00336     sum11 += error_mad / iteration;
00337     sum22 += error_mad / iteration;
00338 
00339     // Calculate a delta for the alfa and gamma parameters
00340     determinant = sum11 * sum22 - sum12 * sum12;
00341     if (fabs(determinant) < ROUNDING_ERROR)
00342     {
00343       // Almost singular matrix. Try without the damping factor.
00344       sum11 -= error_mad / iteration;
00345       sum22 -= error_mad / iteration;
00346       determinant = sum11 * sum22 - sum12 * sum12;
00347     }
00348     if (fabs(determinant) < ROUNDING_ERROR) break; // Still singular
00349     delta_alfa = (sum13 * sum22 - sum23 * sum12) / determinant;
00350     delta_gamma = (sum23 * sum11 - sum13 * sum12) / determinant;
00351 
00352     // Stop when we are close enough and have tried hard enough
00353     if (fabs(delta_alfa) + fabs(delta_gamma) < ACCURACY && iteration > 3)
00354       break;
00355 
00356     // New values for the next iteration
00357     alfa += delta_alfa;
00358     gamma += delta_gamma;
00359 
00360     // Limit the parameters in their allowed range.
00361     if (alfa > max_alfa)
00362       alfa = max_alfa;
00363     else if (alfa < min_alfa)
00364       alfa = min_alfa;
00365     if (gamma > max_gamma)
00366       gamma = max_gamma;
00367     else if (gamma < min_gamma)
00368       gamma = min_gamma;
00369 
00370     // Verify repeated running with both parameters at the boundary
00371     if ((gamma == min_gamma || gamma == max_gamma)
00372       && (alfa == min_alfa || alfa == max_alfa))
00373     {
00374       if (boundarytested++ > 2) break;
00375     }
00376   }
00377 
00378   // Keep the best result
00379   constant_i = best_constant_i;
00380   trend_i = best_trend_i;
00381 
00382   // Echo the result
00383   if (debug)
00384     logger << (fcst ? fcst->getName() : "") << ": double exponential : "
00385       << "alfa " << best_alfa
00386       << ", gamma " << best_gamma
00387       << ", mad " << best_error_mad
00388       << ", " << iteration << " iterations"
00389       << ", constant " << constant_i
00390       << ", trend " << trend_i
00391       << ", forecast " << (trend_i + constant_i) << endl;
00392   return best_error_mad;
00393 }
00394 
00395 
00396 void Forecast::DoubleExponential::applyForecast
00397   (Forecast* forecast, const Date buckets[], unsigned int bucketcount, bool debug)
00398 {
00399   // Loop over all buckets and set the forecast to a linearly changing value
00400   for (unsigned int i = 1; i < bucketcount; ++i)
00401     if (constant_i + i * trend_i > 0)
00402       forecast->setTotalQuantity(
00403         DateRange(buckets[i-1], buckets[i]),
00404         constant_i + i * trend_i
00405         );
00406 }
00407 
00408 
00409 }       // end namespace

Generated on 16 Apr 2010 for frePPLe by  doxygen 1.6.1