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 }