64f2a6721d8bc5f984d14e615d1b5f7bd0cf3fba
[charm.git] / 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 /*@}*/