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, "Logging specific to SURF (lagrange)");
  22 XBT_LOG_NEW_SUBCATEGORY(surf_lagrange_dichotomy, surf, "Logging specific to SURF (lagrange dichotomy)");
  23
  24 /*
  25  * Local prototypes to implement the lagrangian optimization with optimal step, also called dichotomy.
  26  */
  27 //solves the proportional fairness using a lagrange optimizition with dichotomy step
  28 void   lagrange_solve       (lmm_system_t sys);
  29 //computes the value of the dichotomy using a initial values, init, with a specific variable or constraint
  30 double dichotomy(double init, double diff(double, void*), void *var_cnst, double min_error);
  31 //computes the value of the differential of variable param_var applied to mu
  32 double partial_diff_mu      (double mu, void * param_var);
  33 //computes the value of the differential of constraint param_cnst applied to lambda
  34 double partial_diff_lambda  (double lambda, void * param_cnst);
  35 //auxiliar function to compute the partial_diff
  36 double diff_aux(lmm_variable_t var, double x);
  37
  38
  39 static int __check_kkt(xbt_swag_t cnst_list, xbt_swag_t var_list,int warn)
  40 {
  41   xbt_swag_t elem_list  = NULL;
  42   lmm_element_t elem    = NULL;
  43   lmm_constraint_t cnst = NULL;
  44   lmm_variable_t var    = NULL;
  45
  46   double tmp;
  47
  48   //verify the KKT property for each link
  49   xbt_swag_foreach(cnst, cnst_list){
  50     tmp = 0;
  51     elem_list = &(cnst->element_set);
  52     xbt_swag_foreach(elem, elem_list) {
  53       var = elem->variable;
  54       if(var->weight<=0) continue;
  55       tmp += var->value;
  56     }
  57
  58     if(double_positive(tmp - cnst->bound)) {
  59       if(warn) WARN3("The link (%p) is over-used. Expected less than %f and got %f", cnst, cnst->bound, tmp);
  60       return 0;
  61     }
  62     DEBUG3("Checking KKT for constraint (%p): sat = %f, lambda = %f ",cnst, tmp - cnst->bound, cnst->lambda);
  63
  64 /*     if(!((fabs(tmp - cnst->bound)<MAXMIN_PRECISION && cnst->lambda>=MAXMIN_PRECISION) || */
  65 /*       (fabs(tmp - cnst->bound)>=MAXMIN_PRECISION && cnst->lambda<MAXMIN_PRECISION))) { */
  66 /*       if(warn) WARN1("The KKT condition is not verified for cnst %p...", cnst); */
  67 /*       return 0; */
  68 /*     } */
  69   }
  70
  71   //verify the KKT property of each flow
  72   xbt_swag_foreach(var, var_list){
  73     if(var->bound < 0 || var->weight <= 0) continue;
  74     DEBUG3("Checking KKT for variable (%p): sat = %f mu = %f",var, var->value - var->bound,var->mu);
  75
  76     if(double_positive(var->value - var->bound)) {
  77       if(warn) WARN3("The variable (%p) is too large. Expected less than %f and got %f", var, var->bound, var->value);
  78       return 0;
  79     }
  80
  81 /*     if(!((fabs(var->value - var->bound)<MAXMIN_PRECISION && var->mu>=MAXMIN_PRECISION) || */
  82 /*       (fabs(var->value - var->bound)>=MAXMIN_PRECISION && var->mu<MAXMIN_PRECISION))) { */
  83 /*       if(warn) WARN1("The KKT condition is not verified for var %p...",var); */
  84 /*       return 0; */
  85 /*     } */
  86   }
  87   return 1;
  88 }
  89
  90 void lagrange_solve(lmm_system_t sys)
  91 {
  92   /*
  93    * Lagrange Variables.
  94    */
  95   int max_iterations= 10000;
  96   double epsilon_min_error  = 1e-6;
  97   double dichotomy_min_error = 1e-8;
  98   double overall_error = 1;
  99
 100   /*
 101    * Variables to manipulate the data structure proposed to model the maxmin
 102    * fairness. See docummentation for more details.
 103    */
 104   xbt_swag_t cnst_list  = NULL;
 105   lmm_constraint_t cnst = NULL;
 106
 107   xbt_swag_t var_list   = NULL;
 108   lmm_variable_t var    = NULL;
 109
 110   /*
 111    * Auxiliar variables.
 112    */
 113   int iteration=0;
 114   double tmp=0;
 115   int i;
 116
 117
 118   DEBUG0("Iterative method configuration snapshot =====>");
 119   DEBUG1("#### Maximum number of iterations       : %d", max_iterations);
 120   DEBUG1("#### Minimum error tolerated            : %e", epsilon_min_error);
 121   DEBUG1("#### Minimum error tolerated (dichotomy) : %e", 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_error > epsilon_min_error && iteration < max_iterations){
 160
 161     iteration++;
 162     DEBUG1("************** ITERATION %d **************", iteration);
 163
 164     /*
 165      * Compute the value of mu_i
 166      */
 167     //forall mu_i in mu_1, mu_2, ..., mu_n
 168     xbt_swag_foreach(var, var_list) {
 169       if((var->bound >= 0) && (var->weight > 0) ){
 170         var->new_mu = dichotomy(var->mu, partial_diff_mu, var, dichotomy_min_error);
 171         if(var->new_mu < 0) var->new_mu = 0;
 172         DEBUG3("====> var->mu (%p) : %g -> %g", var, var->mu, var->new_mu);
 173         var->mu = var->new_mu;
 174       }
 175     }
 176
 177     /*
 178      * Compute the value of lambda_i
 179      */
 180     //forall lambda_i in lambda_1, lambda_2, ..., lambda_n
 181     xbt_swag_foreach(cnst, cnst_list) {
 182       cnst->new_lambda = dichotomy(cnst->lambda, partial_diff_lambda, cnst, dichotomy_min_error);
 183       DEBUG2("====> cnst->lambda (%p) = %e", cnst, cnst->new_lambda);
 184       cnst->lambda = cnst->new_lambda;
 185     }
 186
 187     /*
 188      * Now computes the values of each variable (\rho) based on
 189      * the values of \lambda and \mu.
 190      */
 191     overall_error=0;
 192     xbt_swag_foreach(var, var_list) {
 193       if(var->weight <=0)
 194         var->value = 0.0;
 195       else {
 196         //compute sigma_i + mu_i
 197         tmp = 0;
 198         for(i=0; i<var->cnsts_number; i++){
 199           tmp += (var->cnsts[i].constraint)->lambda;
 200         }
 201         if(var->bound > 0)
 202           tmp+=var->mu;
 203         DEBUG3("\t Working on var (%p). cost = %e; Df = %e", var, tmp, var->df);
 204
 205         //uses the partial differential inverse function
 206         tmp = var->func_fpi(var, tmp);
 207
 208         //computes de overall_error using normalized value
 209         if(overall_error <  (fabs(var->value - tmp)/tmp) ){
 210           overall_error = (fabs(var->value - tmp)/tmp);
 211         }
 212
 213         var->value = tmp;
 214       }
 215       DEBUG3("======> value of var (%p)  = %e, overall_error = %e", var, var->value, overall_error);
 216     }
 217
 218     if(!__check_kkt(cnst_list,var_list,0)) overall_error=1.0;
 219     DEBUG2("Iteration %d: Overall_error : %f",iteration,overall_error);
 220   }
 221
 222
 223   __check_kkt(cnst_list,var_list,1);
 224
 225   if(overall_error <= epsilon_min_error){
 226     DEBUG1("The method converges in %d iterations.", iteration);
 227   }
 228   if(iteration>= max_iterations) {
 229     WARN1("Method reach %d iterations, which is the maximum number of iterations allowed.", iteration);
 230   }
 231   INFO1("Method converged after %d iterations", iteration);
 232
 233   if(XBT_LOG_ISENABLED(surf_lagrange, xbt_log_priority_debug)) {
 234     lmm_print(sys);
 235   }
 236 }
 237
 238 /*
 239  * Returns a double value corresponding to the result of a dichotomy proccess with
 240  * respect to a given variable/constraint (\mu in the case of a variable or \lambda in
 241  * case of a constraint) and a initial value init.
 242  *
 243  * @param init initial value for \mu or \lambda
 244  * @param diff a function that computes the differential of with respect a \mu or \lambda
 245  * @param var_cnst a pointer to a variable or constraint
 246  * @param min_erro a minimun error tolerated
 247  *
 248  * @return a double correponding to the result of the dichotomyal process
 249  */
 250 double dichotomy(double init, double diff(double, void*), void *var_cnst, double min_error){
 251   double min, max;
 252   double overall_error;
 253   double middle;
 254   double min_diff, max_diff, middle_diff;
 255   double diff_0=0.0;
 256   min = max = init;
 257
 258   if(init == 0){
 259     min = max = 1;
 260   }
 261
 262   min_diff = max_diff = middle_diff = 0.0;
 263   overall_error = 1;
 264
 265   if((diff_0=diff(0.0, var_cnst)) >= 0){
 266     CDEBUG1(surf_lagrange_dichotomy,"====> returning 0.0 (diff = %e)", diff(0.0, var_cnst));
 267     return 0.0;
 268   }
 269
 270   CDEBUG1(surf_lagrange_dichotomy,"====> not detected positive diff in 0 (%e)",diff_0);
 271
 272   while(overall_error > min_error){
 273
 274     min_diff = diff(min, var_cnst);
 275     max_diff = diff(max, var_cnst);
 276
 277     CDEBUG2(surf_lagrange_dichotomy,"DICHOTOMY ===> min = %1.20f , max = %1.20f", min, max);
 278     CDEBUG2(surf_lagrange_dichotomy,"DICHOTOMY ===> diffmin = %1.20f , diffmax = %1.20f", min_diff, max_diff);
 279
 280     if( min_diff > 0 && max_diff > 0 ){
 281       if(min == max){
 282         CDEBUG0(surf_lagrange_dichotomy,"Decreasing min");
 283         min = min / 2.0;
 284       }else{
 285         CDEBUG0(surf_lagrange_dichotomy,"Decreasing max");
 286         max = min;
 287       }
 288     }else if( min_diff < 0 && max_diff < 0 ){
 289       if(min == max){
 290         CDEBUG0(surf_lagrange_dichotomy,"Increasing max");
 291         max = max * 2.0;
 292       }else{
 293         CDEBUG0(surf_lagrange_dichotomy,"Increasing min");
 294         min = max;
 295       }
 296     }else if( min_diff < 0 && max_diff > 0 ){
 297       middle = (max + min)/2.0;
 298       middle_diff = diff(middle, var_cnst);
 299
 300       if(max != 0.0 && min != 0.0){
 301         overall_error = fabs(min - max)/max;
 302       }
 303
 304       if( middle_diff < 0 ){
 305         min = middle;
 306       }else if( middle_diff > 0 ){
 307         max = middle;
 308       }else{
 309         CWARN0(surf_lagrange_dichotomy,"Found an optimal solution with 0 error!");
 310         overall_error = 0;
 311         return middle;
 312       }
 313
 314     }else if(min_diff == 0){
 315       return min;
 316     }else if(max_diff == 0){
 317       return max;
 318     }else if(min_diff > 0 && max_diff < 0){
 319       CWARN0(surf_lagrange_dichotomy,"The impossible happened, partial_diff(min) > 0 && partial_diff(max) < 0");
 320     }else {
 321       CWARN0(surf_lagrange_dichotomy,"diffmin or diffmax are something I don't know, taking no action.");
 322     }
 323   }
 324
 325
 326   CDEBUG1(surf_lagrange_dichotomy,"====> returning %e", (min+max)/2.0);
 327   return ((min+max)/2.0);
 328 }
 329
 330 /*
 331  *
 332  */
 333 double partial_diff_mu(double mu, void *param_var){
 334   double mu_partial=0.0;
 335   double sigma_mu=0.0;
 336   lmm_variable_t var = (lmm_variable_t)param_var;
 337   int i;
 338
 339   //compute sigma_i
 340   for(i=0; i<var->cnsts_number; i++)
 341     sigma_mu += (var->cnsts[i].constraint)->lambda;
 342
 343   //compute sigma_i + mu_i
 344   sigma_mu += mu;
 345
 346   //use auxiliar function passing (sigma_i + mu_i)
 347   mu_partial = diff_aux(var, sigma_mu) ;
 348
 349   //add the RTT limit
 350   mu_partial += var->bound;
 351
 352   return mu_partial;
 353 }
 354
 355 /*
 356  *
 357  */
 358 double partial_diff_lambda(double lambda, void *param_cnst){
 359
 360   int i;
 361   xbt_swag_t elem_list = NULL;
 362   lmm_element_t elem = NULL;
 363   lmm_variable_t var = NULL;
 364   lmm_constraint_t cnst= (lmm_constraint_t) param_cnst;
 365   double lambda_partial=0.0;
 366   double sigma_i=0.0;
 367
 368   elem_list = &(cnst->element_set);
 369
 370   DEBUG1("Computting diff of cnst (%p)", cnst);
 371
 372   xbt_swag_foreach(elem, elem_list) {
 373     var = elem->variable;
 374     if(var->weight<=0) continue;
 375
 376     //initilize de sumation variable
 377     sigma_i = 0.0;
 378
 379     //compute sigma_i of variable var
 380     for(i=0; i<var->cnsts_number; i++){
 381       sigma_i += (var->cnsts[i].constraint)->lambda;
 382     }
 383
 384     //add mu_i if this flow has a RTT constraint associated
 385     if(var->bound > 0) sigma_i += var->mu;
 386
 387     //replace value of cnst->lambda by the value of parameter lambda
 388     sigma_i = (sigma_i - cnst->lambda) + lambda;
 389
 390     //use the auxiliar function passing (\sigma_i + \mu_i)
 391     lambda_partial += diff_aux(var, sigma_i);
 392   }
 393
 394
 395   lambda_partial += cnst->bound;
 396
 397
 398   CDEBUG1(surf_lagrange_dichotomy,"returnning = %1.20f", lambda_partial);
 399
 400   return lambda_partial;
 401 }
 402
 403
 404 double diff_aux(lmm_variable_t var, double x){
 405   double tmp_fp, tmp_fpi, tmp_fpip, result;
 406
 407   xbt_assert0(var->func_fp, "Initialize the protocol functions first create variables before.");
 408
 409   tmp_fp = var->func_fp(var, x);
 410   tmp_fpi = var->func_fpi(var, x);
 411   tmp_fpip = var->func_fpip(var, x);
 412
 413   result = tmp_fpip*(var->func_fp(var, tmp_fpi));
 414
 415   result = result - tmp_fpi;
 416
 417   result = result - (tmp_fpip * x);
 418
 419   CDEBUG2(surf_lagrange_dichotomy,"diff_aux(%1.20f) = %1.20f", x, result);
 420   return result;
 421 }
 422
 423
 424
 425
 426
 427
 428
 429