Logo AND Algorithmique Numérique Distribuée

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