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