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 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 double partial_diff_mu(double mu, void *param_var);
  35 //computes the value of the differential of constraint param_cnst applied to lambda
  36 double partial_diff_lambda(double lambda, void *param_cnst);
  37 //auxiliar function to compute the partial_diff
  38 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 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         CWARN2(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.",
 330                min, max-min);
 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-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-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 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 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 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 1000.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