src/surf/lagrange.c

   1 /*      $Id$     */
   2
   3 /* Copyright (c) 2007 Arnaud Legrand, Pedro Velho. All rights reserved.     */
   4
   5 /* This program is free software; you can redistribute it and/or modify it
   6  * under the terms of the license (GNU LGPL) which comes with this package. */
   7
   8 /*
   9  * Modelling the proportional fairness using the Lagrange Optimization
  10  * Approach. For a detailed description see:
  11  * "ssh://username@scm.gforge.inria.fr/svn/memo/people/pvelho/lagrange/ppf.ps".
  12  */
  13 #include "xbt/log.h"
  14 #include "xbt/sysdep.h"
  15 #include "xbt/mallocator.h"
  16 #include "maxmin_private.h"
  17
  18 #include <stdlib.h>
  19 #ifndef MATH
  20 #include <math.h>
  21 #endif
  22
  23 #define LAMBDA_STEP 0.01
  24
  25
  26 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_lagrange, surf,
  27                                 "Logging specific to SURF (lagrange)");
  28
  29 XBT_LOG_NEW_SUBCATEGORY(surf_writelambda, surf,
  30                                 "Generates the lambda.in file. WARNING: the size of this file might be a few GBs.");
  31
  32 void lagrange_solve(lmm_system_t sys);
  33
  34 void lagrange_solve(lmm_system_t sys)
  35 {
  36   /*
  37    * Lagrange Variables.
  38    */
  39   int max_iterations= 1000000;
  40   double epsilon_min_error = 0.00001;
  41   double overall_error = 1;
  42   double sigma_step = LAMBDA_STEP;
  43   //double capacity_error=0, bound_error=0;
  44   int watch_out = 0;
  45
  46   /*
  47    * Variables to manipulate the data structure proposed to model the maxmin
  48    * fairness. See docummentation for more details.
  49    */
  50   xbt_swag_t elem_list = NULL;
  51   //lmm_element_t elem = NULL;
  52   lmm_element_t elem1 = NULL;
  53
  54
  55   xbt_swag_t cnst_list = NULL;
  56   //lmm_constraint_t cnst = NULL;
  57   lmm_constraint_t cnst1 = NULL;
  58   //lmm_constraint_t cnst2 = NULL;
  59
  60
  61   xbt_swag_t var_list = NULL;
  62   lmm_variable_t var1 = NULL;
  63   lmm_variable_t var2 = NULL;
  64
  65   /*
  66    * Auxiliar variables.
  67    */
  68   int iteration=0;
  69   double mu_partial=0;
  70   double lambda_partial=0;
  71   double tmp=0;
  72   int i,j;
  73   FILE *gnuplot_file=NULL;
  74   //char print_buf[1024];
  75   //char *trace_buf=xbt_malloc0(sizeof(char));
  76   //double sum;
  77
  78
  79   DEBUG0("Iterative method configuration snapshot =====>");
  80   DEBUG1("#### Maximum number of iterations : %d", max_iterations);
  81   DEBUG1("#### Minimum error tolerated      : %e", epsilon_min_error);
  82   DEBUG1("#### Step                         : %e", sigma_step);
  83
  84
  85   if ( !(sys->modified))
  86     return;
  87
  88   /*
  89    * Initialize the var list variable with only the active variables.
  90    * Associate an index in the swag variables. Initialize mu.
  91    */
  92   var_list = &(sys->variable_set);
  93   i=0;
  94   xbt_swag_foreach(var1, var_list) {
  95     if((var1->bound < 0.0) || (var1->weight <= 0.0)){
  96       DEBUG1("#### NOTE var1(%d) is a boundless variable", i);
  97       var1->mu = -1.0;
  98     } else{
  99       var1->mu = 1.0;
 100       var1->new_mu = 2.0;
 101     }
 102     DEBUG2("#### var1(%d)->mu:  %e", i, var1->mu);
 103     DEBUG2("#### var1(%d)->weight: %e", i, var1->weight);
 104     i++;
 105   }
 106
 107   /*
 108    * Initialize lambda.
 109    */
 110   cnst_list=&(sys->active_constraint_set);
 111   xbt_swag_foreach(cnst1, cnst_list) {
 112     cnst1->lambda = 1.0;
 113     cnst1->new_lambda = 2.0;
 114     DEBUG2("#### cnst1(%p)->lambda:  %e", cnst1, cnst1->lambda);
 115   }
 116
 117   if(XBT_LOG_ISENABLED(surf_writelambda, xbt_log_priority_debug)) {
 118     gnuplot_file = fopen("lambda.in", "w");
 119     fprintf(gnuplot_file, "# iteration    lambda1  lambda2 lambda3 ... lambdaP\n");
 120   }
 121
 122
 123   /*
 124    * While doesn't reach a minimun error or a number maximum of iterations.
 125    */
 126   while(overall_error > epsilon_min_error && iteration < max_iterations){
 127     iteration++;
 128
 129     /*                        d Dual
 130      * Compute the value of ----------- (\lambda^k, \mu^k) this portion
 131      *                       d \mu_i^k
 132      * of code depends on function f(x).
 133      */
 134     var_list = &(sys->variable_set);
 135     xbt_swag_foreach(var1, var_list) {
 136       mu_partial = 0;
 137       if((var1->bound >= 0) && (var1->weight > 0) ){
 138         //for each link with capacity cnsts[i] that uses flow of variable var1 do
 139         for(i=0; i<var1->cnsts_number; i++)
 140           mu_partial += (var1->cnsts[i].constraint)->lambda + var1->mu;
 141
 142         mu_partial = -1.0 / mu_partial + var1->bound;
 143         var1->new_mu = var1->mu - sigma_step * mu_partial;
 144
 145         if(var1->new_mu < 0){
 146           var1->new_mu = 0;
 147         }
 148       }
 149     }
 150
 151
 152     /*                         d Dual
 153      * Compute the value of ------------- (\lambda^k, \mu^k) this portion
 154      *                      d \lambda_i^k
 155      * of code depends on function f(x).
 156      */
 157     j=0;
 158     if(XBT_LOG_ISENABLED(surf_writelambda, xbt_log_priority_debug)) {
 159       fprintf(gnuplot_file, "\n%d",iteration);
 160     }
 161     xbt_swag_foreach(cnst1, cnst_list) {
 162       j++;
 163
 164       lambda_partial = 0;
 165
 166       elem_list = &(cnst1->element_set);
 167       watch_out=0;
 168       xbt_swag_foreach(elem1, elem_list) {
 169
 170         var2 = elem1->variable;
 171
 172         if(var2->weight<=0) continue;
 173
 174         tmp = 0;
 175
 176         for(i=0; i<var2->cnsts_number; i++){
 177           tmp += (var2->cnsts[i].constraint)->lambda;
 178         }
 179         if(var2->bound > 0)
 180           tmp += var2->mu;
 181
 182
 183         if(tmp==0) break;
 184
 185         if (tmp==cnst1->lambda)
 186           watch_out=1;
 187         lambda_partial += (-1.0 / tmp);
 188       }
 189
 190       if(tmp == 0)
 191         cnst1->new_lambda = LAMBDA_STEP;
 192       else {
 193         lambda_partial += cnst1->bound;
 194         if(watch_out && (lambda_partial>0)) {
 195           /* INFO6("Watch Out (%d) %p! lambda_partial: %e; lambda : %e ; (%e %e) \n",iteration, cnst1,  */
 196           /*            lambda_partial, cnst1->lambda, cnst1->lambda / 2, */
 197           /*            cnst1->lambda - sigma_step * lambda_partial); */
 198
 199           if(cnst1->lambda < 0) WARN2("Value of cnst1->lambda(%p) = %e < 0", cnst1, cnst1->lambda);
 200           if((cnst1->lambda - sigma_step * lambda_partial) < 0) WARN1("Value of lambda_new = %e < 0", (cnst1->lambda - sigma_step * lambda_partial));
 201
 202           if(cnst1->lambda - sigma_step * lambda_partial < cnst1->lambda / 2)
 203             cnst1->new_lambda = cnst1->lambda / 2;
 204           else
 205             cnst1->new_lambda = cnst1->lambda - sigma_step * lambda_partial;
 206         } else
 207           cnst1->new_lambda = cnst1->lambda - sigma_step * lambda_partial;
 208         if(cnst1->new_lambda < 0){
 209           cnst1->new_lambda = 0;
 210         }
 211       }
 212
 213       if(XBT_LOG_ISENABLED(surf_writelambda, xbt_log_priority_debug)) {
 214         fprintf(gnuplot_file, "  %e", cnst1->lambda);
 215       }
 216
 217     }
 218
 219
 220     /*
 221      * Now computes the values of each variable (\rho) based on
 222      * the values of \lambda and \mu.
 223      */
 224     overall_error=0;
 225     xbt_swag_foreach(var1, var_list) {
 226       if(var1->weight <=0)
 227         var1->value = 0.0;
 228       else {
 229         tmp = 0;
 230         for(i=0; i<var1->cnsts_number; i++){
 231           tmp += (var1->cnsts[i].constraint)->lambda;
 232           if(var1->bound > 0)
 233             tmp+=var1->mu;
 234         }
 235
 236         //computes de overall_error
 237         if(overall_error < fabs(var1->value - 1.0/tmp)){
 238           overall_error = fabs(var1->value - 1.0/tmp);
 239         }
 240
 241         var1->value = 1.0 / tmp;
 242       }
 243
 244     }
 245
 246
 247     /* Updating lambda's and mu's */
 248     xbt_swag_foreach(var1, var_list)
 249       if(!((var1->bound > 0.0) || (var1->weight <= 0.0)))
 250         var1->mu = var1->new_mu;
 251
 252
 253     xbt_swag_foreach(cnst1, cnst_list)
 254       cnst1->lambda = cnst1->new_lambda;
 255   }
 256
 257
 258
 259
 260   //verify the KKT property
 261   xbt_swag_foreach(cnst1, cnst_list){
 262     tmp = 0;
 263     elem_list = &(cnst1->element_set);
 264     xbt_swag_foreach(elem1, elem_list) {
 265       var1 = elem1->variable;
 266       if(var1->weight<=0) continue;
 267       tmp += var1->value;
 268     }
 269
 270     tmp = tmp - cnst1->bound;
 271
 272
 273     if(tmp != 0 ||  cnst1->lambda != 0){
 274       WARN4("The link %s(%p) doesn't match the KKT property, value expected (=0) got (lambda=%e) (sum_rho=%e)", (char *)cnst1->id, cnst1, cnst1->lambda, tmp);
 275     }
 276
 277   }
 278
 279
 280   xbt_swag_foreach(var1, var_list){
 281     if(var1->bound <= 0 || var1->weight <= 0) continue;
 282     tmp = 0;
 283     tmp = (var1->value - var1->bound);
 284
 285
 286     if(tmp != 0 ||  var1->mu != 0){
 287       WARN4("The flow %s(%p) doesn't match the KKT property, value expected (=0) got (lambda=%e) (sum_rho=%e)", (char *)var1->id, var1, var1->mu, tmp);
 288     }
 289
 290   }
 291
 292
 293
 294
 295   if(overall_error <= epsilon_min_error){
 296     DEBUG1("The method converge in %d iterations.", iteration);
 297   }else{
 298     WARN1("Method reach %d iterations, which is the maxmimun number of iterations allowed.", iteration);
 299   }
 300
 301
 302   if(XBT_LOG_ISENABLED(surf_writelambda, xbt_log_priority_debug)) {
 303     fclose(gnuplot_file);
 304   }
 305
 306
 307
 308
 309
 310 /*   /\* */
 311 /*    * Now computes the values of each variable (\rho) based on */
 312 /*    * the values of \lambda and \mu. */
 313 /*    *\/ */
 314 /*   var_list = &(sys->variable_set); */
 315 /*   xbt_swag_foreach(var1, var_list) { */
 316 /*     tmp = 0; */
 317 /*     for(i=0; i<var1->cnsts_number; i++){ */
 318 /*       elem1 = &(var1->cnsts[i]); */
 319 /*       tmp += (elem1->constraint)->lambda + var1->mu; */
 320 /*     } */
 321 /*     var1->weight = 1 / tmp; */
 322
 323 /*     DEBUG2("var1->weight (id=%s) : %e", (char *)var1->id, var1->weight); */
 324 /*   } */
 325
 326
 327
 328
 329 }