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