Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
Debug cosmetics.
[simgrid.git] / src / surf / lagrange.c
1 /*      $Id$     */
2 /* Copyright (c) 2007 Arnaud Legrand, Pedro Velho. All rights reserved.     */
3 /* This program is free software; you can redistribute it and/or modify it
4  * under the terms of the license (GNU LGPL) which comes with this package. */
5 /*
6  * Modelling the proportional fairness using the Lagrange Optimization 
7  * Approach. For a detailed description see:
8  * "ssh://username@scm.gforge.inria.fr/svn/memo/people/pvelho/lagrange/ppf.ps".
9  */
10 #include "xbt/log.h"
11 #include "xbt/sysdep.h"
12 #include "xbt/mallocator.h"
13 #include "maxmin_private.h"
14
15 #include <stdlib.h>
16 #ifndef MATH
17 #include <math.h>
18 #endif
19
20
21 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_lagrange, surf,
22                                 "Logging specific to SURF (lagrange)");
23 XBT_LOG_NEW_SUBCATEGORY(surf_lagrange_dichotomy, surf,
24                         "Logging specific to SURF (lagrange dichotomy)");
25
26 /*
27  * Local prototypes to implement the lagrangian optimization with optimal step, also called dichotomy.
28  */
29 //solves the proportional fairness using a lagrange optimizition with dichotomy step
30 void lagrange_solve(lmm_system_t sys);
31 //computes the value of the dichotomy using a initial values, init, with a specific variable or constraint
32 double dichotomy(double init, double diff(double, void *), void *var_cnst,
33                  double min_error);
34 //computes the value of the differential of variable param_var applied to mu  
35 double partial_diff_mu(double mu, void *param_var);
36 //computes the value of the differential of constraint param_cnst applied to lambda  
37 double partial_diff_lambda(double lambda, void *param_cnst);
38 //auxiliar function to compute the partial_diff
39 double diff_aux(lmm_variable_t var, double x);
40
41
42 static int __check_kkt(xbt_swag_t cnst_list, xbt_swag_t var_list, int warn)
43 {
44   xbt_swag_t elem_list = NULL;
45   lmm_element_t elem = NULL;
46   lmm_constraint_t cnst = NULL;
47   lmm_variable_t var = NULL;
48
49   double tmp;
50
51   //verify the KKT property for each link
52   xbt_swag_foreach(cnst, cnst_list) {
53     tmp = 0;
54     elem_list = &(cnst->element_set);
55     xbt_swag_foreach(elem, elem_list) {
56       var = elem->variable;
57       if (var->weight <= 0)
58         continue;
59       tmp += var->value;
60     }
61
62     if (double_positive(tmp - cnst->bound)) {
63       if (warn)
64         WARN3
65             ("The link (%p) is over-used. Expected less than %f and got %f",
66              cnst, cnst->bound, tmp);
67       return 0;
68     }
69     DEBUG3("Checking KKT for constraint (%p): sat = %f, lambda = %f ",
70            cnst, tmp - cnst->bound, cnst->lambda);
71
72 /*     if(!((fabs(tmp - cnst->bound)<MAXMIN_PRECISION && cnst->lambda>=MAXMIN_PRECISION) || */
73 /*       (fabs(tmp - cnst->bound)>=MAXMIN_PRECISION && cnst->lambda<MAXMIN_PRECISION))) { */
74 /*       if(warn) WARN1("The KKT condition is not verified for cnst %p...", cnst); */
75 /*       return 0; */
76 /*     } */
77   }
78
79   //verify the KKT property of each flow
80   xbt_swag_foreach(var, var_list) {
81     if (var->bound < 0 || var->weight <= 0)
82       continue;
83     DEBUG3("Checking KKT for variable (%p): sat = %f mu = %f", var,
84            var->value - var->bound, var->mu);
85
86     if (double_positive(var->value - var->bound)) {
87       if (warn)
88         WARN3
89             ("The variable (%p) is too large. Expected less than %f and got %f",
90              var, var->bound, var->value);
91       return 0;
92     }
93
94 /*     if(!((fabs(var->value - var->bound)<MAXMIN_PRECISION && var->mu>=MAXMIN_PRECISION) || */
95 /*       (fabs(var->value - var->bound)>=MAXMIN_PRECISION && var->mu<MAXMIN_PRECISION))) { */
96 /*       if(warn) WARN1("The KKT condition is not verified for var %p...",var); */
97 /*       return 0; */
98 /*     } */
99   }
100   return 1;
101 }
102
103 void lagrange_solve(lmm_system_t sys)
104 {
105   /*
106    * Lagrange Variables.
107    */
108   int max_iterations = 100;
109   double epsilon_min_error = 1e-6;
110   double dichotomy_min_error = 1e-20;
111   double overall_error = 1;
112
113   /*
114    * Variables to manipulate the data structure proposed to model the maxmin
115    * fairness. See docummentation for more details.
116    */
117   xbt_swag_t cnst_list = NULL;
118   lmm_constraint_t cnst = NULL;
119
120   xbt_swag_t var_list = NULL;
121   lmm_variable_t var = NULL;
122
123   /*
124    * Auxiliar variables.
125    */
126   int iteration = 0;
127   double tmp = 0;
128   int i;
129
130
131   DEBUG0("Iterative method configuration snapshot =====>");
132   DEBUG1("#### Maximum number of iterations       : %d", max_iterations);
133   DEBUG1("#### Minimum error tolerated            : %e",
134          epsilon_min_error);
135   DEBUG1("#### Minimum error tolerated (dichotomy) : %e",
136          dichotomy_min_error);
137
138   if (!(sys->modified))
139     return;
140
141   /* 
142    * Initialize the var list variable with only the active variables. 
143    * Associate an index in the swag variables. Initialize mu.
144    */
145   var_list = &(sys->variable_set);
146   i = 0;
147   xbt_swag_foreach(var, var_list) {
148     if ((var->bound < 0.0) || (var->weight <= 0.0)) {
149       DEBUG1("#### NOTE var(%d) is a boundless (or inactive) variable", i);
150       var->mu = -1.0;
151     } else {
152       var->mu = 1.0;
153       var->new_mu = 2.0;
154     }
155     DEBUG3("#### var(%d) %p ->mu :  %e", i, var, var->mu);
156     DEBUG3("#### var(%d) %p ->weight: %e", i, var, var->weight);
157     DEBUG3("#### var(%d) %p ->bound: %e", i, var, var->bound);
158     i++;
159   }
160
161   /* 
162    * Initialize lambda.
163    */
164   cnst_list = &(sys->active_constraint_set);
165   xbt_swag_foreach(cnst, cnst_list) {
166     cnst->lambda = 1.0;
167     cnst->new_lambda = 2.0;
168     DEBUG2("#### cnst(%p)->lambda :  %e", cnst, cnst->lambda);
169   }
170
171   /*
172    * While doesn't reach a minimun error or a number maximum of iterations.
173    */
174   while (overall_error > epsilon_min_error && iteration < max_iterations) {
175     int dual_updated=0;
176
177     iteration++;
178     DEBUG1("************** ITERATION %d **************", iteration);
179     DEBUG0("-------------- Gradient Descent ----------");
180     /*                       
181      * Compute the value of mu_i
182      */
183     //forall mu_i in mu_1, mu_2, ..., mu_n
184     xbt_swag_foreach(var, var_list) {
185       if ((var->bound >= 0) && (var->weight > 0)) {
186         DEBUG1("Working on var (%p)", var);
187         var->new_mu =
188             dichotomy(var->mu, partial_diff_mu, var, dichotomy_min_error);
189         dual_updated += (fabs(var->new_mu-var->mu)>dichotomy_min_error);
190         DEBUG2("dual_updated (%d) : %1.20f",dual_updated,fabs(var->new_mu-var->mu));
191         DEBUG3("Updating mu : var->mu (%p) : %1.20f -> %1.20f", var, var->mu, var->new_mu);
192         var->mu = var->new_mu;
193       }
194     }
195
196     /*
197      * Compute the value of lambda_i
198      */
199     //forall lambda_i in lambda_1, lambda_2, ..., lambda_n
200     xbt_swag_foreach(cnst, cnst_list) {
201       DEBUG1("Working on cnst (%p)", cnst);
202       cnst->new_lambda =
203           dichotomy(cnst->lambda, partial_diff_lambda, cnst,
204                     dichotomy_min_error);
205       dual_updated += (fabs(cnst->new_lambda-cnst->lambda)>dichotomy_min_error);
206       DEBUG2("dual_updated (%d) : %1.20f",dual_updated,fabs(cnst->new_lambda-cnst->lambda));
207       DEBUG3("Updating lambda : cnst->lambda (%p) : %1.20f -> %1.20f", cnst, cnst->lambda, cnst->new_lambda);
208       cnst->lambda = cnst->new_lambda;
209     }
210
211     /*
212      * Now computes the values of each variable (\rho) based on
213      * the values of \lambda and \mu.
214      */
215     DEBUG0("-------------- Check convergence ----------");
216     overall_error = 0;
217     xbt_swag_foreach(var, var_list) {
218       if (var->weight <= 0)
219         var->value = 0.0;
220       else {
221         //compute sigma_i + mu_i
222         tmp = 0;
223         for (i = 0; i < var->cnsts_number; i++) {
224           tmp += (var->cnsts[i].constraint)->lambda;
225         }
226         if (var->bound > 0)
227           tmp += var->mu;
228         DEBUG3("\t Working on var (%p). cost = %e; Df = %e", var, tmp,
229                var->df);
230
231         //uses the partial differential inverse function
232         tmp = var->func_fpi(var, tmp);
233
234         //computes de overall_error using normalized value
235         if (overall_error < (fabs(var->value - tmp))) {
236           overall_error = (fabs(var->value - tmp));
237         }
238
239         var->value = tmp;
240         DEBUG3("New value of var (%p)  = %e, overall_error = %e", var,
241                var->value, overall_error);
242       }
243     }
244
245     if (!__check_kkt(cnst_list, var_list, 0))
246       overall_error = 1.0;
247     DEBUG2("Iteration %d: Overall_error : %f", iteration, overall_error);
248     if(!dual_updated) {
249       WARN1("Could not improve the convergence at iteration %d. Drop it!",iteration);
250       break;
251     }
252   }
253
254
255   __check_kkt(cnst_list, var_list, 1);
256
257   if (overall_error <= epsilon_min_error) {
258     DEBUG1("The method converges in %d iterations.", iteration);
259   }
260   if (iteration >= max_iterations) {
261     WARN1
262         ("Method reach %d iterations, which is the maximum number of iterations allowed.",
263          iteration);
264   }
265 /*   INFO1("Method converged after %d iterations", iteration); */
266
267   if (XBT_LOG_ISENABLED(surf_lagrange, xbt_log_priority_debug)) {
268     lmm_print(sys);
269   }
270 }
271
272 /*
273  * Returns a double value corresponding to the result of a dichotomy proccess with
274  * respect to a given variable/constraint (\mu in the case of a variable or \lambda in
275  * case of a constraint) and a initial value init. 
276  *
277  * @param init initial value for \mu or \lambda
278  * @param diff a function that computes the differential of with respect a \mu or \lambda
279  * @param var_cnst a pointer to a variable or constraint 
280  * @param min_erro a minimun error tolerated
281  *
282  * @return a double correponding to the result of the dichotomyal process
283  */
284 double dichotomy(double init, double diff(double, void *), void *var_cnst,
285                  double min_error)
286 {
287   double min, max;
288   double overall_error;
289   double middle;
290   double min_diff, max_diff, middle_diff;
291   double diff_0 = 0.0;
292   min = max = init;
293
294   XBT_IN;
295
296   if (init == 0.0) {
297     min = max = 0.5;
298   }
299
300   min_diff = max_diff = middle_diff = 0.0;
301   overall_error = 1;
302
303   if ((diff_0 = diff(1e-16, var_cnst)) >= 0) {
304     CDEBUG1(surf_lagrange_dichotomy, "returning 0.0 (diff = %e)",
305             diff_0);
306     XBT_OUT;
307     return 0.0;
308   }
309
310   min_diff = diff(min, var_cnst);
311   max_diff = diff(max, var_cnst);
312
313   while (overall_error > min_error) {
314     CDEBUG4(surf_lagrange_dichotomy,
315             "min, max = %1.20f, %1.20f || diffmin, diffmax %1.20f, %1.20f", min, max,
316             min_diff,max_diff);
317
318     if (min_diff > 0 && max_diff > 0) {
319       if (min == max) {
320         CDEBUG0(surf_lagrange_dichotomy, "Decreasing min");
321         min = min / 2.0;
322         min_diff = diff(min, var_cnst);
323       } else {
324         CDEBUG0(surf_lagrange_dichotomy, "Decreasing max");
325         max = min;
326         max_diff = min_diff;
327
328       }
329     } else if (min_diff < 0 && max_diff < 0) {
330       if (min == max) {
331         CDEBUG0(surf_lagrange_dichotomy, "Increasing max");
332         max = max * 2.0;
333         max_diff = diff(max, var_cnst);
334       } else {
335         CDEBUG0(surf_lagrange_dichotomy, "Increasing min");
336         min = max;
337         min_diff = max_diff;
338       }
339     } else if (min_diff < 0 && max_diff > 0) {
340       middle = (max + min) / 2.0;
341       CDEBUG1(surf_lagrange_dichotomy, "Trying (max+min)/2 : %1.20f",middle);
342
343       if((min==middle) || (max==middle)) {
344         WARN0("Cannot improve the convergence!");
345         break;
346       }
347       middle_diff = diff(middle, var_cnst);
348
349       if (middle_diff < 0) {
350         CDEBUG0(surf_lagrange_dichotomy, "Increasing min");
351         min = middle;
352         min_diff = middle_diff;
353       } else if (middle_diff > 0) {
354         CDEBUG0(surf_lagrange_dichotomy, "Decreasing max");
355         max = middle;
356         max_diff = middle_diff;
357       } else {
358         overall_error = 0;
359       }
360     } else if (min_diff == 0) {
361       max=min;
362       overall_error = 0;
363     } else if (max_diff == 0) {
364       min=max;
365       overall_error = 0;
366     } else if (min_diff > 0 && max_diff < 0) {
367       CWARN0(surf_lagrange_dichotomy,
368              "The impossible happened, partial_diff(min) > 0 && partial_diff(max) < 0");
369       abort();
370     } else {
371       CWARN2(surf_lagrange_dichotomy,
372              "diffmin (%1.20f) or diffmax (%1.20f) are something I don't know, taking no action.",
373              min_diff, max_diff);
374       abort();
375     }
376   }
377
378   CDEBUG1(surf_lagrange_dichotomy, "returning %e", (min + max) / 2.0);
379   XBT_OUT;
380   return ((min + max) / 2.0);
381 }
382
383 /*
384  *
385  */
386 double partial_diff_mu(double mu, void *param_var)
387 {
388   double mu_partial = 0.0;
389   double sigma_mu = 0.0;
390   lmm_variable_t var = (lmm_variable_t) param_var;
391   int i;
392   XBT_IN;
393   //compute sigma_i
394   for (i = 0; i < var->cnsts_number; i++)
395     sigma_mu += (var->cnsts[i].constraint)->lambda;
396
397   //compute sigma_i + mu_i
398   sigma_mu += mu;
399
400   //use auxiliar function passing (sigma_i + mu_i)
401   mu_partial = diff_aux(var, sigma_mu);
402
403   //add the RTT limit
404   mu_partial += var->bound;
405
406   XBT_OUT;
407   return mu_partial;
408 }
409
410 /*
411  *
412  */
413 double partial_diff_lambda(double lambda, void *param_cnst)
414 {
415
416   int i;
417   xbt_swag_t elem_list = NULL;
418   lmm_element_t elem = NULL;
419   lmm_variable_t var = NULL;
420   lmm_constraint_t cnst = (lmm_constraint_t) param_cnst;
421   double lambda_partial = 0.0;
422   double sigma_i = 0.0;
423
424   XBT_IN;
425   elem_list = &(cnst->element_set);
426
427   CDEBUG1(surf_lagrange_dichotomy,"Computting diff of cnst (%p)", cnst);
428
429   xbt_swag_foreach(elem, elem_list) {
430     var = elem->variable;
431     if (var->weight <= 0)
432       continue;
433
434     //initilize de sumation variable
435     sigma_i = 0.0;
436
437     //compute sigma_i of variable var
438     for (i = 0; i < var->cnsts_number; i++) {
439       sigma_i += (var->cnsts[i].constraint)->lambda;
440     }
441
442     //add mu_i if this flow has a RTT constraint associated
443     if (var->bound > 0)
444       sigma_i += var->mu;
445
446     //replace value of cnst->lambda by the value of parameter lambda
447     sigma_i = (sigma_i - cnst->lambda) + lambda;
448
449     //use the auxiliar function passing (\sigma_i + \mu_i)
450     lambda_partial += diff_aux(var, sigma_i);
451   }
452
453
454   lambda_partial += cnst->bound;
455
456
457 /*   CDEBUG1(surf_lagrange_dichotomy, "returning = %1.20f", lambda_partial); */
458
459   XBT_OUT;
460   return lambda_partial;
461 }
462
463
464 double diff_aux(lmm_variable_t var, double x)
465 {
466   double tmp_fpi, result;
467
468   XBT_IN2("(var (%p), x (%1.20f))", var, x);
469   xbt_assert0(var->func_fp,
470               "Initialize the protocol functions first create variables before.");
471
472   tmp_fpi = var->func_fpi(var, x);
473   result = - tmp_fpi;
474
475 /*   CDEBUG1(surf_lagrange_dichotomy, "returning %1.20f", result); */
476   XBT_OUT;
477   return result;
478 }
479
480
481 /**************** Vegas and Reno functions *************************/
482 /*
483  * NOTE for Reno: all functions consider the network
484  * coeficient (alpha) equal to 1.
485  */
486
487 /*
488  * For Vegas f: $\alpha_f d_f \log\left(x_f\right)$
489  */
490 #define VEGAS_SCALING 1000.0
491 double func_vegas_f(lmm_variable_t var, double x){
492   return VEGAS_SCALING*var->df * log(x);
493 }
494
495 /*
496  * For Vegas fp: $\frac{\alpha D_f}{x}$
497  */
498 double func_vegas_fp(lmm_variable_t var, double x){
499   //avoid a disaster value - c'est du bricolage mais ca marche
500 /*   if(x == 0) x = 10e-8; */
501   return VEGAS_SCALING*var->df/x;
502 }
503
504 /*
505  * For Vegas fpi: $\frac{\alpha D_f}{x}$
506  */
507 double func_vegas_fpi(lmm_variable_t var, double x){
508   //avoid a disaster value - c'est du bricolage mais ca marche
509 /*   if(x == 0) x = 10e-8; */
510   return VEGAS_SCALING*var->df/x;
511 }
512
513 /*
514  * For Vegas fpip: $-\frac{\alpha D_f}{x^2}$
515  */
516 double func_vegas_fpip(lmm_variable_t var, double x){
517   //avoid a disaster value - c'est du bricolage mais ca marche
518 /*   if(x == 0) x = 10e-8; */
519   return -( VEGAS_SCALING*var->df/(x*x) ) ;
520 }
521
522
523 /*
524  * For Reno f: $\frac{\sqrt{\frac{3}{2}}}{D_f} \arctan\left(\sqrt{\frac{3}{2}}x_f D_f\right)$
525  */
526 double func_reno_f(lmm_variable_t var, double x){
527   xbt_assert0(var->df>0.0,"Don't call me with stupid values!");
528   // \sqrt{3/2} = 0.8164965808
529   return (0.8164965808 / var->df) * atan( (0.8164965808 / var->df)*x );
530 }
531
532 /*
533  * For Reno fp: $\frac{3}{3 {D_f}^2 x^2 + 2}$
534  */
535 double func_reno_fp(lmm_variable_t var, double x){
536   return 3 / (3*var->df*var->df*x*x + 2);
537 }
538
539 /*
540  * For Reno fpi: $\sqrt{\frac{1}{{D_f}^2 x} - \frac{2}{3{D_f}^2}}$
541  */
542 double func_reno_fpi(lmm_variable_t var, double x){
543   double res_fpi; 
544
545   xbt_assert0(var->df>0.0,"Don't call me with stupid values!");
546   xbt_assert0(x>0.0,"Don't call me with stupid values!");
547
548   res_fpi = 1/(var->df*var->df*x) - 2/(3*var->df*var->df);
549   if(res_fpi<=0.0) return 0.0;
550   xbt_assert0(res_fpi>0.0,"Don't call me with stupid values!");
551   return sqrt(res_fpi);
552 }
553
554 /*
555  * For Reno fpip:  $-\frac{1}{2 {D_f}^2 x^2\sqrt{\frac{1}{{D_f}^2 x} - \frac{2}{3{D_f}^2}}}$
556  */
557 double func_reno_fpip(lmm_variable_t var, double x){
558   double res_fpip; 
559   double critical_test;
560
561   xbt_assert0(var->df>0.0,"Don't call me with stupid values!");
562   xbt_assert0(x>0.0,"Don't call me with stupid values!");
563
564   res_fpip = 1/(var->df*var->df*x) - 2/(3*var->df*var->df);
565   xbt_assert0(res_fpip>0.0,"Don't call me with stupid values!");
566   critical_test = (2*var->df*var->df*x*x*sqrt(res_fpip));
567
568   return -(1.0/critical_test);
569 }