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