src/surf/lagrangedico.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
  24 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_lagrange, surf,
  25                                 "Logging specific to SURF (lagrange)");
  26
  27 XBT_LOG_NEW_SUBCATEGORY(surf_writelambda, surf,
  28                                 "Generates the lambda.in file. WARNING: the size of this file might be a few GBs.");
  29
  30 void lagrange_dicotomi_solve(lmm_system_t sys);
  31
  32 double partial_diff_mu(double mu, lmm_variable_t var1);
  33 double partial_diff_lambda(double lambda, lmm_constraint_t cnst1);
  34
  35 void lagrange_dicotomi_solve(lmm_system_t sys)
  36 {
  37   /*
  38    * Lagrange Variables.
  39    */
  40   int max_iterations= 1000000;
  41   double epsilon_min_error = 0.00001;
  42   double overall_error = 1;
  43   double min, max, middle;
  44
  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 elem1 = NULL;
  52
  53
  54   xbt_swag_t cnst_list = NULL;
  55   lmm_constraint_t cnst1 = NULL;
  56
  57   xbt_swag_t var_list = NULL;
  58   lmm_variable_t var1 = NULL;
  59
  60
  61   /*
  62    * Auxiliar variables.
  63    */
  64   int iteration=0;
  65   double tmp=0;
  66   int i;
  67   FILE *gnuplot_file=NULL;
  68
  69
  70   DEBUG0("Iterative method configuration snapshot =====>");
  71   DEBUG1("#### Maximum number of iterations : %d", max_iterations);
  72   DEBUG1("#### Minimum error tolerated      : %e", epsilon_min_error);
  73
  74
  75   if ( !(sys->modified))
  76     return;
  77
  78   /*
  79    * Initialize the var list variable with only the active variables.
  80    * Associate an index in the swag variables. Initialize mu.
  81    */
  82   var_list = &(sys->variable_set);
  83   i=0;
  84   xbt_swag_foreach(var1, var_list) {
  85     if((var1->bound > 0.0) || (var1->weight <= 0.0)){
  86       DEBUG1("#### NOTE var1(%d) is a boundless variable", i);
  87       var1->mu = -1.0;
  88     } else{
  89       var1->mu =   1.0;
  90       var1->new_mu = 2.0;
  91     }
  92     DEBUG2("#### var1(%d)->mu :  %e", i, var1->mu);
  93     DEBUG2("#### var1(%d)->weight: %e", i, var1->weight);
  94     i++;
  95   }
  96
  97   /*
  98    * Initialize lambda.
  99    */
 100   cnst_list=&(sys->active_constraint_set);
 101   xbt_swag_foreach(cnst1, cnst_list) {
 102     cnst1->lambda = 1.0;
 103     cnst1->new_lambda = 2.0;
 104     DEBUG2("#### cnst1(%p)->lambda :  %e", cnst1, cnst1->lambda);
 105   }
 106
 107
 108
 109   /*
 110    * While doesn't reach a minimun error or a number maximum of iterations.
 111    */
 112   while(overall_error > epsilon_min_error && iteration < max_iterations){
 113
 114     iteration++;
 115
 116
 117     /*                        d Dual
 118      * Compute the value of ----------- (\lambda^k, \mu^k) this portion
 119      *                       d \mu_i^k
 120      * of code depends on function f(x).
 121      */
 122     var_list = &(sys->variable_set);
 123     //forall mu_i in mu_1, mu_2, ..., mu_n
 124     xbt_swag_foreach(var1, var_list) {
 125       if((var1->bound >= 0) && (var1->weight > 0) ){
 126         //for each link with capacity cnsts[i] that uses flow of variable var1 do
 127         //begin dicotomi
 128         min = max = 1.0;
 129         overall_error = 1;
 130         while(overall_error < epsilon_min_error){
 131           if( partial_diff_mu(min, var1)>0 && partial_diff_mu(max, var1)>0 ){
 132             if(min == max){
 133               min = min / 2;
 134             }else{
 135               max = min;
 136             }
 137           }else if( partial_diff_mu(min, var1)<0 && partial_diff_mu(max, var1)<0 ){
 138             if(min == max){
 139               max = max * 2;
 140             }else{
 141             max = min;
 142             }
 143           }else if( partial_diff_mu(min,var1)<0 && partial_diff_mu(max,var1) > 0 ){
 144             if(min == max){
 145               middle =  partial_diff_mu((fabs(min - max)/2), var1);
 146               if( middle > 0 ){
 147                 max = (fabs(min - max)/2);
 148               }else if( middle < 0 ){
 149                 min = (fabs(min - max)/2);
 150               }else{
 151                 WARN0("Found an optimal solution with 0 error!");
 152                 overall_error = 0;
 153               }
 154               overall_error = fabs(min - max);
 155             }
 156           }else{
 157             WARN0("The impossible happened, partial_diff(min) >0 && partial_diff(max) < 0");
 158           }
 159         }
 160
 161         var1->new_mu = max;
 162
 163         if(var1->new_mu < 0){
 164           var1->new_mu = 0;
 165         }
 166       }
 167     }
 168
 169
 170     /*                         d Dual
 171      * Compute the value of ------------- (\lambda^k, \mu^k) this portion
 172      *                      d \lambda_i^k
 173      * of code depends on function f(x).
 174      */
 175     xbt_swag_foreach(cnst1, cnst_list) {
 176
 177       //begin dicotomi
 178       overall_error = 1;
 179       min = max = 1.0;
 180       while(overall_error < epsilon_min_error){
 181         if( partial_diff_lambda(min, cnst1) > 0 && partial_diff_lambda(max, cnst1) > 0 ){
 182           if(min == max){
 183             min = min / 2;
 184           }else{
 185             max = min;
 186           }
 187         }else if( partial_diff_lambda(min, cnst1) < 0 && partial_diff_lambda(max, cnst1) < 0 ){
 188           if(min == max){
 189             max = max * 2;
 190           }else{
 191             max = min;
 192           }
 193         }else if( partial_diff_lambda(min,cnst1) < 0 && partial_diff_lambda(max,cnst1) > 0 ){
 194           if(min == max){
 195             middle = partial_diff_lambda((fabs(min - max)/2), cnst1);
 196             if(  middle > 0 ){
 197               max = (fabs(min - max)/2);
 198             }else if( middle < 0 ){
 199               min = (fabs(min - max)/2);
 200             }else{
 201               WARN0("Found an optimal solution with 0 error!");
 202               overall_error = 0;
 203             }
 204             overall_error = fabs(min - max);
 205           }
 206         }else{
 207           WARN0("The impossible happened, partial_diff(min) >0 && partial_diff(max) < 0");
 208         }
 209       }
 210
 211       var1->new_mu = max;
 212
 213       cnst1->new_lambda = cnst1->lambda;
 214
 215       if(cnst1->new_lambda < 0){
 216         cnst1->new_lambda = 0;
 217         }
 218     }
 219
 220
 221     /*
 222      * Now computes the values of each variable (\rho) based on
 223      * the values of \lambda and \mu.
 224      */
 225     overall_error=0;
 226     xbt_swag_foreach(var1, var_list) {
 227       if(var1->weight <=0)
 228         var1->value = 0.0;
 229       else {
 230         tmp = 0;
 231         for(i=0; i<var1->cnsts_number; i++){
 232           tmp += (var1->cnsts[i].constraint)->lambda;
 233           if(var1->bound > 0)
 234             tmp+=var1->mu;
 235         }
 236
 237         //computes de overall_error
 238         if(overall_error < fabs(var1->value - 1.0/tmp)){
 239           overall_error = fabs(var1->value - 1.0/tmp);
 240         }
 241
 242         var1->value = 1.0 / tmp;
 243       }
 244
 245     }
 246
 247
 248     /* Updating lambda's and mu's */
 249     xbt_swag_foreach(var1, var_list)
 250       if(!((var1->bound > 0.0) || (var1->weight <= 0.0)))
 251         var1->mu = var1->new_mu;
 252
 253
 254     xbt_swag_foreach(cnst1, cnst_list)
 255       cnst1->lambda = cnst1->new_lambda;
 256   }
 257
 258
 259   //verify the KKT property
 260   xbt_swag_foreach(cnst1, cnst_list){
 261     tmp = 0;
 262     elem_list = &(cnst1->element_set);
 263     xbt_swag_foreach(elem1, elem_list) {
 264       var1 = elem1->variable;
 265       if(var1->weight<=0) continue;
 266       tmp += var1->value;
 267     }
 268
 269     tmp = tmp - cnst1->bound;
 270
 271
 272     if(tmp != 0 ||  cnst1->lambda != 0){
 273       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);
 274     }
 275
 276   }
 277
 278
 279   xbt_swag_foreach(var1, var_list){
 280     if(var1->bound <= 0 || var1->weight <= 0) continue;
 281     tmp = 0;
 282     tmp = (var1->value - var1->bound);
 283
 284
 285     if(tmp != 0 ||  var1->mu != 0){
 286       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);
 287     }
 288
 289   }
 290
 291
 292   if(overall_error <= epsilon_min_error){
 293     DEBUG1("The method converge in %d iterations.", iteration);
 294   }else{
 295     WARN1("Method reach %d iterations, which is the maxmimun number of iterations allowed.", iteration);
 296   }
 297
 298
 299   if(XBT_LOG_ISENABLED(surf_writelambda, xbt_log_priority_debug)) {
 300     fclose(gnuplot_file);
 301   }
 302
 303
 304 }
 305
 306
 307
 308 double partial_diff_mu(double mu, lmm_variable_t var1){
 309   double mu_partial=0.0;
 310   int i;
 311
 312   //for each link with capacity cnsts[i] that uses flow of variable var1 do
 313   for(i=0; i<var1->cnsts_number; i++)
 314     mu_partial += (var1->cnsts[i].constraint)->lambda + mu;
 315
 316   mu_partial = (-1.0/mu_partial) + var1->bound;
 317
 318   return mu_partial;
 319 }
 320
 321
 322 double partial_diff_lambda(double lambda, lmm_constraint_t cnst1){
 323
 324   double tmp=0.0;
 325   int i;
 326   double lambda_partial=0.0;
 327   xbt_swag_t elem_list = NULL;
 328   lmm_element_t elem1 = NULL;
 329   lmm_variable_t var1 = NULL;
 330
 331
 332   elem_list = &(cnst1->element_set);
 333
 334   xbt_swag_foreach(elem1, elem_list) {
 335     var1 = elem1->variable;
 336     if(var1->weight<=0) continue;
 337
 338     tmp = 0;
 339     for(i=0; i<var1->cnsts_number; i++){
 340       tmp += (var1->cnsts[i].constraint)->lambda;
 341     }
 342
 343     if(var1->bound > 0)
 344       tmp += var1->mu;
 345
 346     if(tmp==0) lambda_partial = 10e-8;
 347     lambda_partial += (-1.0 / (tmp - 3*cnst1->lambda + 3*cnst1->lambda));
 348   }
 349
 350   return lambda_partial;
 351 }
 352