10 void CentralLB::staticPredictorOn(void* data, void *model)
12 CentralLB *me = (CentralLB*)(data);
13 me->predictorOn((LBPredictorFunction*)model);
16 void CentralLB::staticPredictorOnWin(void* data, void *model, int wind)
18 CentralLB *me = (CentralLB*)(data);
19 me->predictorOn((LBPredictorFunction*)model, wind);
22 void CentralLB::staticPredictorOff(void* data)
24 CentralLB *me = (CentralLB*)(data);
28 void CentralLB::staticChangePredictor(void* data, void *model)
30 CentralLB *me = (CentralLB*)(data);
31 me->changePredictor((LBPredictorFunction*)model);
34 #define MAX_CHISQ_ITER 10000
35 #define SMALL_NUMBER 0.00001 // to avoid singular matrix in gaussj
37 void gaussj(double **a, double *b, int n) {
41 double big, dum, pivinv;
42 // arrays for bookkeeping on the pivoting
43 int *indxc, *indxr, *ipiv;
48 for (j=0; j<n; ++j) ipiv[j]=0;
49 // main loop over the columns to be reduced
52 // outer loop of the search for a pivot element
56 if (ipiv[k] == 0 && fabs(a[j][k]) >= big) {
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;
68 // we are now ready to divide the pivot row by the pivot element, located at irow, icol
71 if (a[icol][icol] == 0) {
72 a[icol][icol] = SMALL_NUMBER;
73 CkPrintf("LB: Singular Matrix\n");
75 pivinv = 1.0/a[icol][icol];
77 for (j=0; j<n; ++j) a[icol][j] *= pivinv;
83 for (k=0; k<n; ++k) a[j][k] -= a[icol][k]*dum;
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;}
99 void Marquardt_coefficients(double *x, double *y, double *param, double **alpha, double *beta, double &chisq, LBPredictorFunction *predict) {
103 double *dyda = new double[predict->num_params];
105 for (i=0; i<predict->num_params; ++i) {
106 for (j=0; j<=i; ++j) alpha[i][j] = 0;
111 // summation loop over all data
112 for (i=0; i<predict->num_params; ++i) {
113 predict->function(x[i], param, ymod, dyda);
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];
119 beta[j++] += dy*dyda[l];
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];
133 bool Marquardt_solver(CentralLB::FutureModel *mod, int object) {
135 double chisq, ochisq;
136 double lambda = 0.001;
139 bool allow_stop = false;
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];
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];
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;
162 Marquardt_coefficients(x,y,mod->parameters[object],alpha,beta,chisq,mod->predictor);
165 while (chisq > 0.01 || !allow_stop) {
166 if (++iterations > MAX_CHISQ_ITER) {
167 // something wrong!!!
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];
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];
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
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];
194 mod->parameters[object][i] = atry[i];
196 } else { // failure, increase lamda
202 for (i=0; i<mod->predictor->num_params; ++i) {
221 // routine that update LDStats given a predictor model
222 void CentralLB::FuturePredictor(BaseLB::LDStats* stats) {
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;
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) {
250 predicted_model->predictor->initialize_params(predicted_model->parameters[i]);
251 predicted_model->predictor->print(predicted_model->parameters[i]);
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]);
260 if (predicted_model->model_valid) {
261 CkPrintf("LB: New model completely constructed\n");
263 CkPrintf("LB: Construction of new model failed\n");
266 } else { // model already constructed, update it
268 double *error_model = new double[stats->n_objs];
269 double *error_default = new double[stats->n_objs];
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]);
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;
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);
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;
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
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]);
311 if (fabs(error_model[i]) < fabs(error_default[i])) predicted_model->model_valid[i] = true;
316 // use the model to update statistics
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);
323 stats->objData[i].cpuTime = predicted_model->predictor->predict(stats->objData[i].cpuTime, param);