src/ck-ldb/CentralPredictor.C

   1 /**
   2  * \addtogroup CkLdb
   3 */
   4 /*@{*/
   5
   6 #include <math.h>
   7 #include <charm++.h>
   8 #include "CentralLB.h"
   9
  10 void CentralLB::staticPredictorOn(void* data, void *model)
  11 {
  12   CentralLB *me = (CentralLB*)(data);
  13   me->predictorOn((LBPredictorFunction*)model);
  14 }
  15
  16 void CentralLB::staticPredictorOnWin(void* data, void *model, int wind)
  17 {
  18   CentralLB *me = (CentralLB*)(data);
  19   me->predictorOn((LBPredictorFunction*)model, wind);
  20 }
  21
  22 void CentralLB::staticPredictorOff(void* data)
  23 {
  24   CentralLB *me = (CentralLB*)(data);
  25   me->predictorOff();
  26 }
  27
  28 void CentralLB::staticChangePredictor(void* data, void *model)
  29 {
  30   CentralLB *me = (CentralLB*)(data);
  31   me->changePredictor((LBPredictorFunction*)model);
  32 }
  33
  34 #define MAX_CHISQ_ITER 10000
  35 #define SMALL_NUMBER   0.00001    // to avoid singular matrix in gaussj
  36
  37 void gaussj(double **a, double *b, int n) {
  38 #if CMK_LBDB_ON
  39   int i,j,k;
  40   int irow, icol;
  41   double big, dum, pivinv;
  42   // arrays for bookkeeping on the pivoting
  43   int *indxc, *indxr, *ipiv;
  44
  45   indxc = new int[n];
  46   indxr = new int[n];
  47   ipiv = new int[n];
  48   for (j=0; j<n; ++j) ipiv[j]=0;
  49   // main loop over the columns to be reduced
  50   for (i=0; i<n; ++i) {
  51     big = 0;
  52     // outer loop of the search for a pivot element
  53     for (j=0; j<n; ++j)
  54       if (ipiv[j] != 1)
  55         for (k=0; k<n; ++k) {
  56           if (ipiv[k] == 0 && fabs(a[j][k]) >= big) {
  57             big = fabs(a[j][k]);
  58             irow=j;
  59             icol=k;
  60           }
  61         }
  62     ++(ipiv[icol]);
  63
  64     if (irow != icol) {
  65       for (j=0; j<n; ++j) {dum=a[irow][j]; a[irow][j]=a[icol][j]; a[icol][j]=dum;}
  66       dum=b[irow]; b[irow]=b[icol]; b[icol]=dum;
  67     }
  68     // we are now ready to divide the pivot row by the pivot element, located at irow, icol
  69     indxr[i]=irow;
  70     indxc[i]=icol;
  71     if (a[icol][icol] == 0) {
  72       a[icol][icol] = SMALL_NUMBER;
  73       CkPrintf("LB: Singular Matrix\n");
  74     }
  75     pivinv = 1.0/a[icol][icol];
  76     a[icol][icol] = 1;
  77     for (j=0; j<n; ++j) a[icol][j] *= pivinv;
  78     b[icol] *= pivinv;
  79     for (j=0; j<n; ++j)
  80       if (j != icol) {
  81         dum = a[j][icol];
  82         a[j][icol] = 0;
  83         for (k=0; k<n; ++k) a[j][k] -= a[icol][k]*dum;
  84         b[j] -= b[icol]*dum;
  85       }
  86   }
  87   // unscramble the matrix
  88   for (i=n-1; i>=0; --i) {
  89     if (indxr[i] != indxc[i])
  90       for (j=0; j<n; ++j) {dum=a[j][indxr[i]]; a[j][indxr[i]]=a[j][indxc[i]]; a[j][indxc[i]]=dum;}
  91   }
  92   delete[] indxr;
  93   delete[] indxc;
  94   delete[] ipiv;
  95
  96 #endif
  97 }
  98
  99 void Marquardt_coefficients(double *x, double *y, double *param, double **alpha, double *beta, double &chisq, LBPredictorFunction *predict) {
 100 #if CMK_LBDB_ON
 101   int i,j,k,l,m;
 102   double ymod, dy;
 103   double *dyda = new double[predict->num_params];
 104
 105   for (i=0; i<predict->num_params; ++i) {
 106     for (j=0; j<=i; ++j) alpha[i][j] = 0;
 107     beta[i]=0;
 108   }
 109   chisq = 0;
 110
 111   // summation loop over all data
 112   for (i=0; i<predict->num_params; ++i) {
 113     predict->function(x[i], param, ymod, dyda);
 114     dy = y[i] - ymod;
 115     for (j=0, l=0; l<predict->num_params; ++l) {
 116       for (k=0, m=0; m<l+1; ++m) {
 117         alpha[j][k++] += dyda[l]*dyda[m];
 118       }
 119       beta[j++] += dy*dyda[l];
 120     }
 121     chisq += dy*dy;
 122   }
 123
 124   // fill the symmetric side
 125   for (j=1; j<predict->num_params; ++j) {
 126     for (k=0; k<j; ++k) alpha[k][j] = alpha[j][k];
 127   }
 128
 129   delete[] dyda;
 130 #endif
 131 }
 132
 133 bool Marquardt_solver(CentralLB::FutureModel *mod, int object) {
 134 #if CMK_LBDB_ON
 135   double chisq, ochisq;
 136   double lambda = 0.001;
 137   int i,j;
 138   int iterations=0;
 139   bool allow_stop = false;
 140
 141   double *oneda = new double[mod->predictor->num_params];
 142   double *atry = new double[mod->predictor->num_params];
 143   double *beta = new double[mod->predictor->num_params];
 144   double *da = new double[mod->predictor->num_params];
 145   double **covar = new double*[mod->predictor->num_params];
 146   double **alpha = new double*[mod->predictor->num_params];
 147   double *x = new double[mod->cur_stats-1];
 148   double *y = new double[mod->cur_stats-1];
 149   double **temp = new double*[mod->predictor->num_params];
 150
 151   for (i=0; i<mod->predictor->num_params; ++i) {
 152     alpha[i] = new double[mod->predictor->num_params];
 153     covar[i] = new double[mod->predictor->num_params];
 154     temp[i] = new double[mod->predictor->num_params];
 155     atry[i] = mod->parameters[object][i];
 156   }
 157   for (i=0; i<mod->cur_stats-2; ++i) {
 158     x[i] = mod->collection[i].objData[object].wallTime;
 159     y[i] = mod->collection[i+1].objData[object].wallTime;
 160   }
 161
 162   Marquardt_coefficients(x,y,mod->parameters[object],alpha,beta,chisq,mod->predictor);
 163   ochisq = chisq;
 164
 165   while (chisq > 0.01 || !allow_stop) {
 166     if (++iterations > MAX_CHISQ_ITER) {
 167       // something wrong!!!
 168       return false;
 169     }
 170     // alter linearized fitting matrix, by augmenting diagonal elements
 171     for (i=0; i<mod->predictor->num_params; ++i) {
 172       for (j=0; j<mod->predictor->num_params; ++j) covar[i][j] = alpha[i][j];
 173       covar[i][i] = alpha[i][i] * (1 + lambda);
 174       for (j=0; j<mod->predictor->num_params; ++j) temp[i][j] = covar[i][j];
 175       oneda[i] = beta[i];
 176     }
 177
 178     // matrix solution
 179     gaussj(temp, oneda, mod->predictor->num_params);
 180     for (i=0; i<mod->predictor->num_params; ++i) {
 181       for (j=0; j<mod->predictor->num_params; ++j) covar[i][j] = temp[i][j];
 182       da[i] = oneda[i];
 183     }
 184
 185     // did the trial succeed?
 186     for (i=0, j=0; j<mod->predictor->num_params; ++j) atry[j] = mod->parameters[object][j] + da[i++];
 187     Marquardt_coefficients(x,y,atry,covar,da,chisq,mod->predictor);
 188     if (chisq < ochisq) {  // success, accept the new solution
 189       lambda *= 0.1;
 190       allow_stop = true;
 191       for (i=0; i<mod->predictor->num_params; ++i) {
 192         for (j=0; j<mod->predictor->num_params; ++j) alpha[i][j] = covar[i][j];
 193         beta[i] = da[i];
 194         mod->parameters[object][i] = atry[i];
 195       }
 196     } else {  // failure, increase lamda
 197       lambda *= 10;
 198       allow_stop = false;
 199     }
 200     ochisq = chisq;
 201   }
 202   for (i=0; i<mod->predictor->num_params; ++i) {
 203     delete[] alpha[i];
 204     delete[] covar[i];
 205     delete[] temp[i];
 206   }
 207   delete[] oneda;
 208   delete[] atry;
 209   delete[] beta;
 210   delete[] da;
 211   delete[] covar;
 212   delete[] alpha;
 213   delete[] x;
 214   delete[] y;
 215   delete[] temp;
 216
 217 #endif
 218   return true;
 219 }
 220
 221 // routine that update LDStats given a predictor model
 222 void CentralLB::FuturePredictor(BaseLB::LDStats* stats) {
 223 #if CMK_LBDB_ON
 224   bool model_done;
 225   int i;
 226
 227   if (predicted_model->cur_stats < _lb_predict_delay) {
 228     // not yet ready to create the model, just store the relevant statistic
 229     predicted_model->collection[predicted_model->start_stats].objData.resize(stats->n_objs);
 230     predicted_model->collection[predicted_model->start_stats].commData.resize(stats->n_comm);
 231     predicted_model->collection[predicted_model->start_stats].n_objs = stats->n_objs;
 232     predicted_model->collection[predicted_model->start_stats].n_migrateobjs = stats->n_migrateobjs;
 233     predicted_model->collection[predicted_model->start_stats].n_comm = stats->n_comm;
 234     for (i=0; i<stats->n_objs; ++i)
 235       predicted_model->collection[predicted_model->start_stats].objData[i] = stats->objData[i];
 236     for (i=0; i<stats->n_comm; ++i)
 237       predicted_model->collection[predicted_model->start_stats].commData[i] = stats->commData[i];
 238     ++predicted_model->cur_stats;
 239     ++predicted_model->start_stats;
 240
 241   } else {
 242
 243     if (predicted_model->parameters == NULL) {     // time to create the new prediction model
 244       // allocate parameters
 245       predicted_model->model_valid = new bool[stats->n_objs];
 246       predicted_model->parameters = new double*[stats->n_objs];
 247       for (i=0; i<stats->n_objs; ++i) predicted_model->parameters[i] = new double[predicted_model->predictor->num_params];
 248       for (i=0; i<stats->n_objs; ++i) {
 249         // initialization
 250         predicted_model->predictor->initialize_params(predicted_model->parameters[i]);
 251         predicted_model->predictor->print(predicted_model->parameters[i]);
 252
 253         model_done = Marquardt_solver(predicted_model, i);
 254         // always initialize to false for conservativity
 255         predicted_model->model_valid[i] = false;
 256         CkPrintf("LB: Model for object %d %s\n",i,model_done?"found":"not found");
 257         predicted_model->predictor->print(predicted_model->parameters[i]);
 258       }
 259
 260       if (predicted_model->model_valid) {
 261         CkPrintf("LB: New model completely constructed\n");
 262       } else {
 263         CkPrintf("LB: Construction of new model failed\n");
 264       }
 265
 266     } else {     // model already constructed, update it
 267
 268       double *error_model = new double[stats->n_objs];
 269       double *error_default = new double[stats->n_objs];
 270
 271       CkPrintf("Error in estimation:\n");
 272       for (i=0; i<stats->n_objs; ++i) {
 273         error_model[i] = stats->objData[i].wallTime-predicted_model->predictor->predict(predicted_model->collection[(predicted_model->start_stats-1)%predicted_model->n_stats].objData[i].wallTime,predicted_model->parameters[i]);
 274         error_default[i] = stats->objData[i].wallTime-predicted_model->collection[(predicted_model->start_stats-1)%predicted_model->n_stats].objData[i].wallTime;
 275         CkPrintf("object %d: real time=%f, model error=%f, default error=%f\n",i,stats->objData[i].wallTime,error_model[i],error_default[i]);
 276       }
 277
 278       // save statistics in the last position
 279       if (predicted_model->start_stats >= predicted_model->n_stats) predicted_model->start_stats -= predicted_model->n_stats;
 280       if (predicted_model->cur_stats < predicted_model->n_stats) ++predicted_model->cur_stats;
 281
 282       predicted_model->collection[predicted_model->start_stats].objData.resize(stats->n_objs);
 283       predicted_model->collection[predicted_model->start_stats].commData.resize(stats->n_comm);
 284
 285       predicted_model->collection[predicted_model->start_stats].n_objs = stats->n_objs;
 286       predicted_model->collection[predicted_model->start_stats].n_migrateobjs = stats->n_migrateobjs;
 287       predicted_model->collection[predicted_model->start_stats].n_comm = stats->n_comm;
 288       for (i=0; i<stats->n_objs; ++i)
 289         predicted_model->collection[predicted_model->start_stats].objData[i] = stats->objData[i];
 290       for (i=0; i<stats->n_comm; ++i)
 291         predicted_model->collection[predicted_model->start_stats].commData[i] = stats->commData[i];
 292       ++predicted_model->start_stats;
 293
 294       // check if model is ok
 295       // the check can be performed even if the model is not valid since it will
 296       // releave which objects are wrongly updated and will try to fix them
 297
 298       // the update of the model is done if the model does not approximate
 299       // sufficiently well the underlining function or if the time-invariante
 300       // approach is performing better
 301       for (i=0; i<stats->n_objs; ++i) {
 302         //if (fabs(error_model[i]) > 0.2*stats->objData[i].wallTime || fabs(error_model[i]) > fabs(error_default[i])) {
 303         if (fabs(error_model[i]) > fabs(error_default[i])) {  // no absolute error check
 304           predicted_model->model_valid[i] = false;
 305           // model wrong, rebuild it now
 306           predicted_model->predictor->initialize_params(predicted_model->parameters[i]);
 307           model_done = Marquardt_solver(predicted_model, i);
 308           CkPrintf("LB: Updated model for object %d %s",i,model_done?"success":"failed. ");
 309           predicted_model->predictor->print(predicted_model->parameters[i]);
 310         }
 311         if (fabs(error_model[i]) < fabs(error_default[i])) predicted_model->model_valid[i] = true;
 312       }
 313
 314     }
 315
 316     // use the model to update statistics
 317     double *param;
 318     for (int i=0; i<stats->n_objs; ++i) {
 319       if (predicted_model->model_valid[i]) {
 320         param = predicted_model->parameters[i];
 321         stats->objData[i].wallTime = predicted_model->predictor->predict(stats->objData[i].wallTime, param);
 322 #if CMK_LB_CPUTIMER
 323         stats->objData[i].cpuTime = predicted_model->predictor->predict(stats->objData[i].cpuTime, param);
 324 #endif
 325       }
 326     }
 327
 328   }
 329
 330 #endif
 331 }
 332
 333 /*@}*/