src/surf/lagrange.cpp

   1 /* Copyright (c) 2007, 2008, 2009, 2010. The SimGrid Team.
   2  * All rights reserved.                                                     */
   3
   4 /* This program is free software; you can redistribute it and/or modify it
   5  * under the terms of the license (GNU LGPL) which comes with this package. */
   6
   7 /*
   8  * Modelling the proportional fairness using the Lagrange Optimization
   9  * Approach. For a detailed description see:
  10  * "ssh://username@scm.gforge.inria.fr/svn/memo/people/pvelho/lagrange/ppf.ps".
  11  */
  12 #include "xbt/log.h"
  13 #include "xbt/sysdep.h"
  14 //#include "maxmin_private.h"
  15 #include "solver.h"
  16 #include "solver.hpp"
  17 #include <stdlib.h>
  18 #ifndef MATH
  19 #include <math.h>
  20 #endif
  21
  22 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_lagrange, surf,
  23                                 "Logging specific to SURF (lagrange)");
  24 XBT_LOG_NEW_SUBCATEGORY(surf_lagrange_dichotomy, surf_lagrange,
  25                         "Logging specific to SURF (lagrange dichotomy)");
  26
  27 #define SHOW_EXPR(expr) XBT_CDEBUG(surf_lagrange,#expr " = %g",expr);
  28
  29 double (*func_f_def) (lmm_variable_t, double);
  30 double (*func_fp_def) (lmm_variable_t, double);
  31 double (*func_fpi_def) (lmm_variable_t, double);
  32
  33 /*
  34  * Local prototypes to implement the lagrangian optimization with optimal step, also called dichotomy.
  35  */
  36 //solves the proportional fairness using a lagrange optimizition with dichotomy step
  37 void lagrange_solve(lmm_system_t sys);
  38 //computes the value of the dichotomy using a initial values, init, with a specific variable or constraint
  39 static double dichotomy(double init, double diff(double, void *),
  40                         void *var_cnst, double min_error);
  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(std::vector<ConstraintPtr> *cnstList, std::vector<VariablePtr> *varList,
  45                             int warn)
  46 {
  47   std::vector<ElementPtr> *elemList = NULL;
  48   lmm_element_t elem = NULL;
  49   lmm_constraint_t cnst = NULL;
  50   lmm_variable_t var = NULL;
  51   std::vector<VariablePtr>::iterator varIt;
  52   std::vector<ElementPtr>::iterator elemIt;
  53   std::vector<ConstraintPtr>::iterator cnstIt;
  54
  55   double tmp;
  56
  57   for (cnstIt=cnstList->begin(); cnstIt!=cnstList->end(); ++cnstIt) {
  58     cnst = (*cnstIt);
  59     tmp = 0;
  60     elemList = &(cnst->m_elementSet);
  61     for (elemIt=elemList->begin(); elemIt!=elemList->end(); ++elemIt) {
  62       var = (*elemIt)->p_variable;
  63       if (var->m_weight <= 0)
  64         continue;
  65       tmp += var->m_value;
  66     }
  67
  68     if (double_positive(tmp - cnst->m_bound)) {
  69       if (warn)
  70         XBT_WARN
  71             ("The link (%p) is over-used. Expected less than %f and got %f",
  72              cnst, cnst->m_bound, tmp);
  73       return 0;
  74     }
  75     XBT_DEBUG
  76         ("Checking feasability for constraint (%p): sat = %f, lambda = %f ",
  77          cnst, tmp - cnst->m_bound, cnst->m_lambda);
  78   }
  79
  80   for (varIt=varList->begin(); varIt!=varList->end(); ++varIt) {
  81     if (!var->m_weight)
  82       break;
  83     if (var->m_bound < 0)
  84       continue;
  85     XBT_DEBUG("Checking feasability for variable (%p): sat = %f mu = %f", var,
  86            var->m_value - var->m_bound, var->m_mu);
  87
  88     if (double_positive(var->m_value - var->m_bound)) {
  89       if (warn)
  90         XBT_WARN
  91             ("The variable (%p) is too large. Expected less than %f and got %f",
  92              var, var->m_bound, var->m_value);
  93       return 0;
  94     }
  95   }
  96   return 1;
  97 }
  98
  99 static double new_value(lmm_variable_t var)
 100 {
 101   double tmp = 0;
 102   int i;
 103   std::vector<ElementPtr>::iterator elemIt;
 104
 105   for (elemIt=var->m_cnsts.begin(); elemIt!=var->m_cnsts.end(); ++elemIt) {
 106     tmp += ((*elemIt)->p_constraint)->m_lambda;
 107   }
 108   if (var->m_bound > 0)
 109     tmp += var->m_mu;
 110   XBT_DEBUG("\t Working on var (%p). cost = %e; Weight = %e", var, tmp,
 111          var->m_weight);
 112   //uses the partial differential inverse function
 113   return var->p_funcFPI(var, tmp);
 114 }
 115
 116 static double new_mu(lmm_variable_t var)
 117 {
 118   double mu_i = 0.0;
 119   double sigma_i = 0.0;
 120   int j;
 121   std::vector<ElementPtr>::iterator elemIt;
 122
 123   for (elemIt=var->m_cnsts.begin(); elemIt!=var->m_cnsts.end(); ++elemIt) {
 124     sigma_i += ((*elemIt)->p_constraint)->m_lambda;
 125   }
 126   mu_i = var->p_funcFP(var, var->m_bound) - sigma_i;
 127   if (mu_i < 0.0)
 128     return 0.0;
 129   return mu_i;
 130 }
 131
 132 static double dual_objective(std::vector<VariablePtr> *varList, std::vector<ConstraintPtr> *cnstList)
 133 {
 134   lmm_constraint_t cnst = NULL;
 135   lmm_variable_t var = NULL;
 136
 137   double obj = 0.0;
 138   std::vector<VariablePtr>::iterator varIt;
 139   std::vector<ElementPtr>::iterator elemIt;
 140   std::vector<ConstraintPtr>::iterator cnstIt;
 141
 142   for (varIt=varList->begin(); varIt!=varList->end(); ++varIt) {
 143     var = (*varIt);
 144     double sigma_i = 0.0;
 145     int j;
 146
 147     if (!var->m_weight)
 148       break;
 149
 150     for (elemIt=var->m_cnsts.begin(); elemIt!=var->m_cnsts.end(); ++elemIt)
 151       sigma_i += ((*elemIt)->p_constraint)->m_lambda;
 152
 153     if (var->m_bound > 0)
 154       sigma_i += var->m_mu;
 155
 156     XBT_DEBUG("var %p : sigma_i = %1.20f", var, sigma_i);
 157
 158     obj += var->p_funcF(var, var->p_funcFPI(var, sigma_i)) -
 159         sigma_i * var->p_funcFPI(var, sigma_i);
 160
 161     if (var->m_bound > 0)
 162       obj += var->m_mu * var->m_bound;
 163   }
 164
 165   for (cnstIt=cnstList->begin(); cnstIt!=cnstList->end(); ++cnstIt)
 166     obj += (*cnstIt)->m_lambda * (*cnstIt)->m_bound;
 167
 168   return obj;
 169 }
 170
 171 void lagrange_solve(lmm_system_t sys)
 172 {
 173   /*
 174    * Lagrange Variables.
 175    */
 176   int max_iterations = 100;
 177   double epsilon_min_error = MAXMIN_PRECISION;
 178   double dichotomy_min_error = 1e-14;
 179   double overall_modification = 1;
 180
 181   /*
 182    * Variables to manipulate the data structure proposed to model the maxmin
 183    * fairness. See docummentation for more details.
 184    */
 185   std::vector<ConstraintPtr> *cnstList = NULL;
 186   std::vector<ConstraintPtr>::iterator cnstIt;
 187   lmm_constraint_t cnst = NULL;
 188
 189   std::vector<VariablePtr> *varList = NULL;
 190   std::vector<VariablePtr>::iterator varIt;
 191   lmm_variable_t var = NULL;
 192
 193   std::vector<ElementPtr>::iterator elemIt;
 194
 195   /*
 196    * Auxiliar variables.
 197    */
 198   int iteration = 0;
 199   double tmp = 0;
 200   int i;
 201   double obj, new_obj;
 202
 203   XBT_DEBUG("Iterative method configuration snapshot =====>");
 204   XBT_DEBUG("#### Maximum number of iterations       : %d", max_iterations);
 205   XBT_DEBUG("#### Minimum error tolerated            : %e",
 206          epsilon_min_error);
 207   XBT_DEBUG("#### Minimum error tolerated (dichotomy) : %e",
 208          dichotomy_min_error);
 209
 210   if (XBT_LOG_ISENABLED(surf_lagrange, xbt_log_priority_debug)) {
 211     lmm_print(sys);
 212   }
 213
 214   if (!(sys->m_modified))
 215     return;
 216
 217   /*
 218    * Initialize lambda.
 219    */
 220   cnstList = &(sys->m_activeConstraintSet);
 221   for (cnstIt=cnstList->begin(); cnstIt!=cnstList->end(); ++cnstIt) {
 222     cnst = *cnstIt;
 223     cnst->m_lambda = 1.0;
 224     cnst->m_newLambda = 2.0;
 225     XBT_DEBUG("#### cnst(%p)->lambda :  %e", cnst, cnst->m_lambda);
 226   }
 227
 228   /*
 229    * Initialize the var list variable with only the active variables.
 230    * Associate an index in the swag variables. Initialize mu.
 231    */
 232   varList = &(sys->m_variableSet);
 233   i = 0;
 234   for (varIt=varList->begin(); varIt!=varList->end(); ++varIt) {
 235     var = *varIt;
 236     if (!var->m_weight)
 237       var->m_value = 0.0;
 238     else {
 239       int nb = 0;
 240       if (var->m_bound < 0.0) {
 241         XBT_DEBUG("#### NOTE var(%d) is a boundless variable", i);
 242         var->m_mu = -1.0;
 243         var->m_value = new_value(var);
 244       } else {
 245         var->m_mu = 1.0;
 246         var->m_newMu = 2.0;
 247         var->m_value = new_value(var);
 248       }
 249       XBT_DEBUG("#### var(%p) ->weight :  %e", var, var->m_weight);
 250       XBT_DEBUG("#### var(%p) ->mu :  %e", var, var->m_mu);
 251       XBT_DEBUG("#### var(%p) ->weight: %e", var, var->m_weight);
 252       XBT_DEBUG("#### var(%p) ->bound: %e", var, var->m_bound);
 253       for (elemIt=var->m_cnsts.begin(); elemIt!=var->m_cnsts.end(); ++elemIt) {
 254         if ((*elemIt)->m_value == 0.0)
 255           nb++;
 256       }
 257       if (nb == var->m_cnsts.size())
 258         var->m_value = 1.0;
 259     }
 260   }
 261
 262   /*
 263    * Compute dual objective.
 264    */
 265   obj = dual_objective(varList, cnstList);
 266
 267   /*
 268    * While doesn't reach a minimun error or a number maximum of iterations.
 269    */
 270   while (overall_modification > epsilon_min_error
 271          && iteration < max_iterations) {
 272 /*     int dual_updated=0; */
 273
 274     iteration++;
 275     XBT_DEBUG("************** ITERATION %d **************", iteration);
 276     XBT_DEBUG("-------------- Gradient Descent ----------");
 277
 278     /*
 279      * Improve the value of mu_i
 280      */
 281     for (varIt=varList->begin(); varIt!=varList->end(); ++varIt) {
 282       var = *varIt;
 283       if (!var->m_weight)
 284         break;
 285       if (var->m_bound >= 0) {
 286         XBT_DEBUG("Working on var (%p)", var);
 287         var->m_newMu = new_mu(var);
 288 /*   dual_updated += (fabs(var->new_mu-var->mu)>dichotomy_min_error); */
 289 /*   XBT_DEBUG("dual_updated (%d) : %1.20f",dual_updated,fabs(var->new_mu-var->mu)); */
 290         XBT_DEBUG("Updating mu : var->mu (%p) : %1.20f -> %1.20f", var,
 291                var->m_mu, var->m_newMu);
 292         var->m_mu = var->m_newMu;
 293
 294         new_obj = dual_objective(varList, cnstList);
 295         XBT_DEBUG("Improvement for Objective (%g -> %g) : %g", obj, new_obj,
 296                obj - new_obj);
 297         xbt_assert(obj - new_obj >= -epsilon_min_error,
 298                     "Our gradient sucks! (%1.20f)", obj - new_obj);
 299         obj = new_obj;
 300       }
 301     }
 302
 303     /*
 304      * Improve the value of lambda_i
 305      */
 306     for (cnstIt=cnstList->begin(); cnstIt!=cnstList->end(); ++cnstIt) {
 307       cnst = *cnstIt;
 308       XBT_DEBUG("Working on cnst (%p)", cnst);
 309       cnst->m_newLambda =
 310           dichotomy(cnst->m_lambda, partial_diff_lambda, cnst,
 311                     dichotomy_min_error);
 312 /*       dual_updated += (fabs(cnst->new_lambda-cnst->lambda)>dichotomy_min_error); */
 313 /*       XBT_DEBUG("dual_updated (%d) : %1.20f",dual_updated,fabs(cnst->new_lambda-cnst->lambda)); */
 314       XBT_DEBUG("Updating lambda : cnst->lambda (%p) : %1.20f -> %1.20f",
 315              cnst, cnst->m_lambda, cnst->m_newLambda);
 316       cnst->m_lambda = cnst->m_newLambda;
 317
 318       new_obj = dual_objective(varList, cnstList);
 319       XBT_DEBUG("Improvement for Objective (%g -> %g) : %g", obj, new_obj,
 320              obj - new_obj);
 321       xbt_assert(obj - new_obj >= -epsilon_min_error,
 322                   "Our gradient sucks! (%1.20f)", obj - new_obj);
 323       obj = new_obj;
 324     }
 325
 326     /*
 327      * Now computes the values of each variable (\rho) based on
 328      * the values of \lambda and \mu.
 329      */
 330     XBT_DEBUG("-------------- Check convergence ----------");
 331     overall_modification = 0;
 332     for (varIt=varList->begin(); varIt!=varList->end(); ++varIt) {
 333       var = *varIt;
 334       if (var->m_weight <= 0)
 335         var->m_value = 0.0;
 336       else {
 337         tmp = new_value(var);
 338
 339         overall_modification =
 340             MAX(overall_modification, fabs(var->m_value - tmp));
 341
 342         var->m_value = tmp;
 343         XBT_DEBUG("New value of var (%p)  = %e, overall_modification = %e",
 344                var, var->m_value, overall_modification);
 345       }
 346     }
 347
 348     XBT_DEBUG("-------------- Check feasability ----------");
 349     if (!__check_feasible(cnstList, varList, 0))
 350       overall_modification = 1.0;
 351     XBT_DEBUG("Iteration %d: overall_modification : %f", iteration,
 352            overall_modification);
 353 /*     if(!dual_updated) { */
 354 /*       XBT_WARN("Could not improve the convergence at iteration %d. Drop it!",iteration); */
 355 /*       break; */
 356 /*     } */
 357   }
 358   __check_feasible(cnstList, varList, 1);
 359
 360   if (overall_modification <= epsilon_min_error) {
 361     XBT_DEBUG("The method converges in %d iterations.", iteration);
 362   }
 363   if (iteration >= max_iterations) {
 364     XBT_DEBUG
 365         ("Method reach %d iterations, which is the maximum number of iterations allowed.",
 366          iteration);
 367   }
 368 /*   XBT_INFO("Method converged after %d iterations", iteration); */
 369
 370   if (XBT_LOG_ISENABLED(surf_lagrange, xbt_log_priority_debug)) {
 371     lmm_print(sys);
 372   }
 373 }
 374
 375 /*
 376  * Returns a double value corresponding to the result of a dichotomy proccess with
 377  * respect to a given variable/constraint (\mu in the case of a variable or \lambda in
 378  * case of a constraint) and a initial value init.
 379  *
 380  * @param init initial value for \mu or \lambda
 381  * @param diff a function that computes the differential of with respect a \mu or \lambda
 382  * @param var_cnst a pointer to a variable or constraint
 383  * @param min_erro a minimun error tolerated
 384  *
 385  * @return a double correponding to the result of the dichotomyal process
 386  */
 387 static double dichotomy(double init, double diff(double, void *),
 388                         void *var_cnst, double min_error)
 389 {
 390   #ifdef TOREPAIR
 391   double min, max;
 392   double overall_error;
 393   double middle;
 394   double min_diff, max_diff, middle_diff;
 395   double diff_0 = 0.0;
 396   min = max = init;
 397
 398   XBT_IN();
 399
 400   if (init == 0.0) {
 401     min = max = 0.5;
 402   }
 403
 404   min_diff = max_diff = middle_diff = 0.0;
 405   overall_error = 1;
 406
 407   if ((diff_0 = diff(1e-16, var_cnst)) >= 0) {
 408     XBT_CDEBUG(surf_lagrange_dichotomy, "returning 0.0 (diff = %e)", diff_0);
 409     XBT_OUT();
 410     return 0.0;
 411   }
 412
 413   min_diff = diff(min, var_cnst);
 414   max_diff = diff(max, var_cnst);
 415
 416   while (overall_error > min_error) {
 417     XBT_CDEBUG(surf_lagrange_dichotomy,
 418             "[min, max] = [%1.20f, %1.20f] || diffmin, diffmax = %1.20f, %1.20f",
 419             min, max, min_diff, max_diff);
 420
 421     if (min_diff > 0 && max_diff > 0) {
 422       if (min == max) {
 423         XBT_CDEBUG(surf_lagrange_dichotomy, "Decreasing min");
 424         min = min / 2.0;
 425         min_diff = diff(min, var_cnst);
 426       } else {
 427         XBT_CDEBUG(surf_lagrange_dichotomy, "Decreasing max");
 428         max = min;
 429         max_diff = min_diff;
 430       }
 431     } else if (min_diff < 0 && max_diff < 0) {
 432       if (min == max) {
 433         XBT_CDEBUG(surf_lagrange_dichotomy, "Increasing max");
 434         max = max * 2.0;
 435         max_diff = diff(max, var_cnst);
 436       } else {
 437         XBT_CDEBUG(surf_lagrange_dichotomy, "Increasing min");
 438         min = max;
 439         min_diff = max_diff;
 440       }
 441     } else if (min_diff < 0 && max_diff > 0) {
 442       middle = (max + min) / 2.0;
 443       XBT_CDEBUG(surf_lagrange_dichotomy, "Trying (max+min)/2 : %1.20f",
 444               middle);
 445
 446       if ((min == middle) || (max == middle)) {
 447         XBT_CWARN(surf_lagrange_dichotomy,
 448                "Cannot improve the convergence! min=max=middle=%1.20f, diff = %1.20f."
 449                " Reaching the 'double' limits. Maybe scaling your function would help ([%1.20f,%1.20f]).",
 450                min, max - min, min_diff, max_diff);
 451         break;
 452       }
 453       middle_diff = diff(middle, var_cnst);
 454
 455       if (middle_diff < 0) {
 456         XBT_CDEBUG(surf_lagrange_dichotomy, "Increasing min");
 457         min = middle;
 458         overall_error = max_diff - middle_diff;
 459         min_diff = middle_diff;
 460 /*   SHOW_EXPR(overall_error); */
 461       } else if (middle_diff > 0) {
 462         XBT_CDEBUG(surf_lagrange_dichotomy, "Decreasing max");
 463         max = middle;
 464         overall_error = max_diff - middle_diff;
 465         max_diff = middle_diff;
 466 /*   SHOW_EXPR(overall_error); */
 467       } else {
 468         overall_error = 0;
 469 /*   SHOW_EXPR(overall_error); */
 470       }
 471     } else if (min_diff == 0) {
 472       max = min;
 473       overall_error = 0;
 474 /*       SHOW_EXPR(overall_error); */
 475     } else if (max_diff == 0) {
 476       min = max;
 477       overall_error = 0;
 478 /*       SHOW_EXPR(overall_error); */
 479     } else if (min_diff > 0 && max_diff < 0) {
 480       XBT_CWARN(surf_lagrange_dichotomy,
 481              "The impossible happened, partial_diff(min) > 0 && partial_diff(max) < 0");
 482       xbt_abort();
 483     } else {
 484       XBT_CWARN(surf_lagrange_dichotomy,
 485              "diffmin (%1.20f) or diffmax (%1.20f) are something I don't know, taking no action.",
 486              min_diff, max_diff);
 487       xbt_abort();
 488     }
 489   }
 490
 491   XBT_CDEBUG(surf_lagrange_dichotomy, "returning %e", (min + max) / 2.0);
 492   XBT_OUT();
 493   return ((min + max) / 2.0);
 494   #endif
 495 }
 496
 497 static double partial_diff_lambda(double lambda, void *param_cnst)
 498 {
 499   #ifdef TOREPAIR
 500   int j;
 501   xbt_swag_t elem_list = NULL;
 502   lmm_element_t elem = NULL;
 503   lmm_variable_t var = NULL;
 504   lmm_constraint_t cnst = (lmm_constraint_t) param_cnst;
 505   double diff = 0.0;
 506   double sigma_i = 0.0;
 507
 508   XBT_IN();
 509   elem_list = &(cnst->element_set);
 510
 511   XBT_CDEBUG(surf_lagrange_dichotomy, "Computing diff of cnst (%p)", cnst);
 512
 513   xbt_swag_foreach(elem, elem_list) {
 514     var = elem->variable;
 515     if (var->weight <= 0)
 516       continue;
 517
 518     XBT_CDEBUG(surf_lagrange_dichotomy, "Computing sigma_i for var (%p)",
 519             var);
 520     // Initialize the summation variable
 521     sigma_i = 0.0;
 522
 523     // Compute sigma_i
 524     for (j = 0; j < var->cnsts_number; j++) {
 525       sigma_i += (var->cnsts[j].constraint)->lambda;
 526     }
 527
 528     //add mu_i if this flow has a RTT constraint associated
 529     if (var->bound > 0)
 530       sigma_i += var->mu;
 531
 532     //replace value of cnst->lambda by the value of parameter lambda
 533     sigma_i = (sigma_i - cnst->lambda) + lambda;
 534
 535     diff += -var->func_fpi(var, sigma_i);
 536   }
 537
 538
 539   diff += cnst->bound;
 540
 541   XBT_CDEBUG(surf_lagrange_dichotomy,
 542           "d D/d lambda for cnst (%p) at %1.20f = %1.20f", cnst, lambda,
 543           diff);
 544   XBT_OUT();
 545   return diff;
 546   #endif
 547 }
 548
 549 /** \brief Attribute the value bound to var->bound.
 550  *
 551  *  \param func_fpi  inverse of the partial differential of f (f prime inverse, (f')^{-1})
 552  *
 553  *  Set default functions to the ones passed as parameters. This is a polimorfism in C pure, enjoy the roots of programming.
 554  *
 555  */
 556 void lmm_set_default_protocol_function(double (*func_f) (lmm_variable_t var, double x),
 557                                        double (*func_fp) (lmm_variable_t var, double x),
 558                                        double (*func_fpi) (lmm_variable_t var, double x))
 559 {
 560   func_f_def = func_f;
 561   func_fp_def = func_fp;
 562   func_fpi_def = func_fpi;
 563 }
 564
 565
 566 /**************** Vegas and Reno functions *************************/
 567 /*
 568  * NOTE for Reno: all functions consider the network
 569  * coeficient (alpha) equal to 1.
 570  */
 571
 572 /*
 573  * For Vegas: $f(x) = \alpha D_f\ln(x)$
 574  * Therefore: $fp(x) = \frac{\alpha D_f}{x}$
 575  * Therefore: $fpi(x) = \frac{\alpha D_f}{x}$
 576  */
 577 #define VEGAS_SCALING 1000.0
 578
 579 double func_vegas_f(lmm_variable_t var, double x)
 580 {
 581   #ifdef TOREPAIR
 582   xbt_assert(x > 0.0, "Don't call me with stupid values! (%1.20f)", x);
 583   return VEGAS_SCALING * var->weight * log(x);
 584   #endif
 585 }
 586
 587 double func_vegas_fp(lmm_variable_t var, double x)
 588 {
 589   #ifdef TOREPAIR
 590   xbt_assert(x > 0.0, "Don't call me with stupid values! (%1.20f)", x);
 591   return VEGAS_SCALING * var->weight / x;
 592   #endif
 593 }
 594
 595 double func_vegas_fpi(lmm_variable_t var, double x)
 596 {
 597   #ifdef TOREPAIR
 598   xbt_assert(x > 0.0, "Don't call me with stupid values! (%1.20f)", x);
 599   return var->weight / (x / VEGAS_SCALING);
 600   #endif
 601 }
 602
 603 /*
 604  * For Reno:  $f(x) = \frac{\sqrt{3/2}}{D_f} atan(\sqrt{3/2}D_f x)$
 605  * Therefore: $fp(x)  = \frac{3}{3 D_f^2 x^2+2}$
 606  * Therefore: $fpi(x)  = \sqrt{\frac{1}{{D_f}^2 x} - \frac{2}{3{D_f}^2}}$
 607  */
 608 #define RENO_SCALING 1.0
 609 double func_reno_f(lmm_variable_t var, double x)
 610 {
 611   xbt_assert(var->m_weight > 0.0, "Don't call me with stupid values!");
 612
 613   return RENO_SCALING * sqrt(3.0 / 2.0) / var->m_weight *
 614       atan(sqrt(3.0 / 2.0) * var->m_weight * x);
 615 }
 616
 617 double func_reno_fp(lmm_variable_t var, double x)
 618 {
 619   return RENO_SCALING * 3.0 / (3.0 * var->m_weight * var->m_weight * x * x +
 620                                2.0);
 621 }
 622
 623 double func_reno_fpi(lmm_variable_t var, double x)
 624 {
 625   double res_fpi;
 626
 627   xbt_assert(var->m_weight > 0.0, "Don't call me with stupid values!");
 628   xbt_assert(x > 0.0, "Don't call me with stupid values!");
 629
 630   res_fpi =
 631       1.0 / (var->m_weight * var->m_weight * (x / RENO_SCALING)) -
 632       2.0 / (3.0 * var->m_weight * var->m_weight);
 633   if (res_fpi <= 0.0)
 634     return 0.0;
 635 //   xbt_assert(res_fpi>0.0,"Don't call me with stupid values!");
 636   return sqrt(res_fpi);
 637 }
 638
 639
 640 /* Implementing new Reno-2
 641  * For Reno-2:  $f(x)   = U_f(x_f) = \frac{{2}{D_f}}*ln(2+x*D_f)$
 642  * Therefore:   $fp(x)  = 2/(Weight*x + 2)
 643  * Therefore:   $fpi(x) = (2*Weight)/x - 4
 644  */
 645 #define RENO2_SCALING 1.0
 646 double func_reno2_f(lmm_variable_t var, double x)
 647 {
 648   xbt_assert(var->m_weight > 0.0, "Don't call me with stupid values!");
 649   return RENO2_SCALING * (1.0 / var->m_weight) * log((x * var->m_weight) /
 650                                                    (2.0 * x * var->m_weight +
 651                                                     3.0));
 652 }
 653
 654 double func_reno2_fp(lmm_variable_t var, double x)
 655 {
 656   return RENO2_SCALING * 3.0 / (var->m_weight * x *
 657                                 (2.0 * var->m_weight * x + 3.0));
 658 }
 659
 660 double func_reno2_fpi(lmm_variable_t var, double x)
 661 {
 662   double res_fpi;
 663   double tmp;
 664
 665   xbt_assert(x > 0.0, "Don't call me with stupid values!");
 666   tmp = x * var->m_weight * var->m_weight;
 667   res_fpi = tmp * (9.0 * x + 24.0);
 668
 669   if (res_fpi <= 0.0)
 670     return 0.0;
 671
 672   res_fpi = RENO2_SCALING * (-3.0 * tmp + sqrt(res_fpi)) / (4.0 * tmp);
 673   return res_fpi;
 674 }