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