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= 10;
  39   double epsilon_min_error = 1e-10;
  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 _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 = var1->mu;
 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->mu = max;
 159
 160         if(var1->mu < 0){
 161           var1->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
 175       DEBUG2("cnst1 (id=%s) (%p)", (char *)cnst1->id, cnst1);
 176
 177       //begin dicotomi
 178       i=0;
 179       overall_error = 1;
 180       min = max = cnst1->lambda;
 181       while(overall_error > epsilon_min_error){
 182         i++;
 183
 184
 185         //      DEBUG4("====> Dicotomi debug. [%e, %e], D(min,max) = [%e, %e]", min, max, partial_diff_lambda(min, cnst1), partial_diff_lambda(max, cnst1));
 186
 187         if( partial_diff_lambda(min, cnst1) > 0 && partial_diff_lambda(max, cnst1) > 0 ){
 188           if(min == max){
 189             min = min / 2.0;
 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             max = max * 2.0;
 196           }else{
 197             min = max;
 198           }
 199         }else if( partial_diff_lambda(min,cnst1) < 0 && partial_diff_lambda(max,cnst1) > 0 ){
 200           middle = (max + min)/2.0;
 201
 202
 203           //DEBUG2("Ideal state reached middle = %e, D(fabs(min-max)/2.0) = %e", middle, partial_diff_lambda(middle, cnst1));
 204           if( partial_diff_lambda(middle, cnst1) < 0 ){
 205             min = middle;
 206           }else if( partial_diff_lambda(middle, cnst1) > 0 ){
 207             max = middle;
 208           }else{
 209             WARN0("Found an optimal solution with 0 error!");
 210             overall_error = 0;
 211           }
 212           overall_error = fabs(min - max);
 213         }else{
 214           WARN0("The impossible happened, partial_diff(min) >0 && partial_diff(max) < 0");
 215         }
 216       }
 217
 218
 219       DEBUG1("Number of iteration in the dicotomi %d", i);
 220
 221       cnst1->lambda = min;
 222
 223       if(cnst1->lambda < 0){
 224         cnst1->lambda = 0;
 225       }
 226     }
 227
 228
 229     /*
 230      * Now computes the values of each variable (\rho) based on
 231      * the values of \lambda and \mu.
 232      */
 233     overall_error=0;
 234     DEBUG1("Iteration %d ", iteration);
 235     xbt_swag_foreach(var1, var_list) {
 236       if(var1->weight <=0)
 237         var1->value = 0.0;
 238       else {
 239         tmp = 0;
 240         for(i=0; i<var1->cnsts_number; i++){
 241           tmp += (var1->cnsts[i].constraint)->lambda;
 242           if(var1->bound > 0)
 243             tmp+=var1->mu;
 244         }
 245
 246         //computes de overall_error
 247         if(overall_error < fabs(var1->value - 1.0/tmp)){
 248           overall_error = fabs(var1->value - 1.0/tmp);
 249         }
 250
 251         var1->value = 1.0 / tmp;
 252       }
 253
 254
 255       DEBUG2("======> value of var1 (%p)  = %e", var1, var1->value);
 256     }
 257   }
 258
 259
 260
 261
 262
 263   //verify the KKT property
 264   xbt_swag_foreach(cnst1, cnst_list){
 265     tmp = 0;
 266     elem_list = &(cnst1->element_set);
 267     xbt_swag_foreach(elem1, elem_list) {
 268       var1 = elem1->variable;
 269       if(var1->weight<=0) continue;
 270       tmp += var1->value;
 271     }
 272
 273     tmp = tmp - cnst1->bound;
 274
 275
 276     if(tmp != 0 ||  cnst1->lambda != 0){
 277       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);
 278     }
 279
 280   }
 281
 282
 283   xbt_swag_foreach(var1, var_list){
 284     if(var1->bound <= 0 || var1->weight <= 0) continue;
 285     tmp = 0;
 286     tmp = (var1->value - var1->bound);
 287
 288
 289     if(tmp != 0 ||  var1->mu != 0){
 290       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);
 291     }
 292
 293   }
 294
 295
 296   if(overall_error <= epsilon_min_error){
 297     DEBUG1("The method converge in %d iterations.", iteration);
 298   }else{
 299     WARN1("Method reach %d iterations, which is the maxmimun number of iterations allowed.", iteration);
 300   }
 301
 302
 303
 304 }
 305
 306
 307 double dicotomi(double init, void *diff(double, void*), void *var_cnst){
 308   double min, max;
 309   double overall_error;
 310
 311   min = max = init;
 312   overall_error = 1;
 313
 314   while(overall_error > epsilon_min_error){
 315     if( diff(min, var_cnst) > 0 && diff(max, var_cnst) > 0 ){
 316       if(min == max){
 317         min = min / 2.0;
 318       }else{
 319         max = min;
 320       }
 321     }else if( diff(min, var_cnst) < 0 && diff(max, var_cnst) < 0 ){
 322       if(min == max){
 323         max = max * 2.0;
 324       }else{
 325         min = max;
 326       }
 327     }else if( diff(min, var_cnst) < 0 && diff(max, var_cnst) > 0 ){
 328       middle = (max + min)/2.0;
 329
 330       if( diff(middle, var_cnst) < 0 ){
 331         min = middle;
 332       }else if( diff(middle, var_cnst) > 0 ){
 333         max = middle;
 334       }else{
 335         WARN0("Found an optimal solution with 0 error!");
 336         overall_error = 0;
 337       }
 338       overall_error = fabs(min - max);
 339     }else{
 340       WARN0("The impossible happened, partial_diff(min) >0 && partial_diff(max) < 0");
 341     }
 342   }
 343 }
 344
 345 double partial_diff_mu(double mu, lmm_variable_t var1){
 346   double mu_partial=0.0;
 347   int i;
 348
 349   //for each link with capacity cnsts[i] that uses flow of variable var1 do
 350   for(i=0; i<var1->cnsts_number; i++)
 351     mu_partial += (var1->cnsts[i].constraint)->lambda + mu;
 352
 353   mu_partial = (-1.0/mu_partial) + var1->bound;
 354
 355   return mu_partial;
 356 }
 357
 358
 359 double partial_diff_lambda(double lambda, lmm_constraint_t cnst1){
 360
 361   double tmp=0.0;
 362   int i;
 363   double lambda_partial=0.0;
 364   xbt_swag_t elem_list = NULL;
 365   lmm_element_t elem1 = NULL;
 366   lmm_variable_t var1 = NULL;
 367
 368
 369   elem_list = &(cnst1->element_set);
 370
 371   xbt_swag_foreach(elem1, elem_list) {
 372     var1 = elem1->variable;
 373     if(var1->weight<=0) continue;
 374
 375     tmp = 0;
 376     for(i=0; i<var1->cnsts_number; i++){
 377       tmp += (var1->cnsts[i].constraint)->lambda;
 378     }
 379
 380     if(var1->bound > 0)
 381       tmp += var1->mu;
 382
 383     tmp = tmp - cnst1->lambda + lambda;
 384
 385     //un peux du bricolage pour evite la catastrophe
 386     if(tmp==0) lambda_partial = 10e-8;
 387
 388     lambda_partial += (-1.0 /tmp);
 389   }
 390
 391
 392   lambda_partial += cnst1->bound;
 393
 394   //DEBUG3("Partial diff lambda result cnst1 %s (%p) : %e", (char *)cnst1->id, cnst1, lambda_partial);
 395   return lambda_partial;
 396 }
 397