Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
bdffdc4d418223d5443d168d82ccef5e3fa4c46f
[simgrid.git] / src / surf / lagrange.c
1 /*      $Id$     */
2 /* Copyright (c) 2007 Arnaud Legrand, Pedro Velho. All rights reserved.     */
3 /* This program is free software; you can redistribute it and/or modify it
4  * under the terms of the license (GNU LGPL) which comes with this package. */
5 /*
6  * Modelling the proportional fairness using the Lagrange Optimization 
7  * Approach. For a detailed description see:
8  * "ssh://username@scm.gforge.inria.fr/svn/memo/people/pvelho/lagrange/ppf.ps".
9  */
10 #include "xbt/log.h"
11 #include "xbt/sysdep.h"
12 #include "xbt/mallocator.h"
13 #include "maxmin_private.h"
14
15 #include <stdlib.h>
16 #ifndef MATH
17 #include <math.h>
18 #endif
19
20
21 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_lagrange, surf, "Logging specific to SURF (lagrange)");
22
23
24 /*
25  * Local prototypes to implement the lagrangian optimization with optimal step, also called dicotomi.
26  */
27 //solves the proportional fairness using a lagrange optimizition with dicotomi step
28 void   lagrange_solve       (lmm_system_t sys);
29 //computes the value of the dicotomi using a initial values, init, with a specific variable or constraint
30 double dicotomi(double init, double diff(double, void*), void *var_cnst, double min_error);
31 //computes the value of the differential of variable param_var applied to mu  
32 double partial_diff_mu      (double mu, void * param_var);
33 //computes the value of the differential of constraint param_cnst applied to lambda  
34 double partial_diff_lambda  (double lambda, void * param_cnst);
35
36
37
38 void lagrange_solve(lmm_system_t sys)
39 {
40   /*
41    * Lagrange Variables.
42    */
43   int max_iterations= 10000;
44   double epsilon_min_error = 1e-4;
45   double dicotomi_min_error = 1e-8;
46   double overall_error = 1;
47
48
49   /*
50    * Variables to manipulate the data structure proposed to model the maxmin
51    * fairness. See docummentation for more details.
52    */
53   xbt_swag_t elem_list  = NULL;
54   lmm_element_t elem    = NULL;
55
56   xbt_swag_t cnst_list  = NULL;
57   lmm_constraint_t cnst = NULL;
58   
59   xbt_swag_t var_list   = NULL;
60   lmm_variable_t var    = NULL;
61
62
63   /*
64    * Auxiliar variables.
65    */
66   int iteration=0;
67   double tmp=0;
68   int i;
69    
70
71   DEBUG0("Iterative method configuration snapshot =====>");
72   DEBUG1("#### Maximum number of iterations       : %d", max_iterations);
73   DEBUG1("#### Minimum error tolerated            : %e", epsilon_min_error);  
74   DEBUG1("#### Minimum error tolerated (dicotomi) : %e", dicotomi_min_error);
75
76   if ( !(sys->modified))
77     return;
78
79   /* 
80    * Initialize the var list variable with only the active variables. 
81    * Associate an index in the swag variables. Initialize mu.
82    */
83   var_list = &(sys->variable_set);
84   i=0;
85   xbt_swag_foreach(var, var_list) {    
86     if((var->bound > 0.0) || (var->weight <= 0.0)){
87       DEBUG1("#### NOTE var(%d) is a boundless variable", i);
88       var->mu = -1.0;
89     } else{ 
90       var->mu =   1.0;
91       var->new_mu = 2.0;
92     }
93     DEBUG2("#### var(%d)->mu :  %e", i, var->mu);
94     DEBUG2("#### var(%d)->weight: %e", i, var->weight);
95     i++;
96   }
97
98   /* 
99    * Initialize lambda.
100    */
101   cnst_list=&(sys->active_constraint_set); 
102   xbt_swag_foreach(cnst, cnst_list){
103     cnst->lambda = 1.0;
104     cnst->new_lambda = 2.0;
105     DEBUG2("#### cnst(%p)->lambda :  %e", cnst, cnst->lambda);
106   }
107   
108   /*
109    * While doesn't reach a minimun error or a number maximum of iterations.
110    */
111   while(overall_error > epsilon_min_error && iteration < max_iterations){
112    
113     iteration++;
114     DEBUG1("************** ITERATION %d **************", iteration);    
115
116     /*                       
117      * Compute the value of mu_i
118      */
119     //forall mu_i in mu_1, mu_2, ..., mu_n
120     xbt_swag_foreach(var, var_list) {
121       if((var->bound >= 0) && (var->weight > 0) ){
122         var->new_mu = dicotomi(var->mu, partial_diff_mu, var, dicotomi_min_error);
123         if(var->new_mu < 0) var->new_mu = 0;
124       }
125     }
126
127     /*
128      * Compute the value of lambda_i
129      */
130     //forall lambda_i in lambda_1, lambda_2, ..., lambda_n
131     xbt_swag_foreach(cnst, cnst_list) {
132       cnst->new_lambda = dicotomi(cnst->lambda, partial_diff_lambda, cnst, dicotomi_min_error);
133       DEBUG2("====> cnst->lambda (%p) = %e", cnst, cnst->new_lambda);      
134     }
135
136     /*                       
137      * Update values of mu and lambda
138      */
139     //forall mu_i in mu_1, mu_2, ..., mu_n
140     xbt_swag_foreach(var, var_list) {
141       var->mu = var->new_mu ;
142     }
143   
144     //forall lambda_i in lambda_1, lambda_2, ..., lambda_n
145     xbt_swag_foreach(cnst, cnst_list) {
146       cnst->lambda = cnst->new_lambda;
147     }
148
149     /*
150      * Now computes the values of each variable (\rho) based on
151      * the values of \lambda and \mu.
152      */
153     overall_error=0;   
154     xbt_swag_foreach(var, var_list) {
155       if(var->weight <=0) 
156         var->value = 0.0;
157       else {
158         tmp = 0;
159         for(i=0; i<var->cnsts_number; i++){
160           tmp += (var->cnsts[i].constraint)->lambda;
161           if(var->bound > 0) 
162             tmp+=var->mu;
163         }
164
165         if(tmp == 0.0)
166           WARN0("CAUTION: division by 0.0");
167
168         //computes de overall_error
169         if(overall_error < fabs(var->value - 1.0/tmp)){
170           overall_error = fabs(var->value - 1.0/tmp);
171         }
172         var->value = 1.0 / tmp;
173       }
174       DEBUG4("======> value of var %s (%p)  = %e, overall_error = %e", (char *)var->id, var, var->value, overall_error);       
175     }
176   }
177
178
179   //verify the KKT property for each link
180   xbt_swag_foreach(cnst, cnst_list){
181     tmp = 0;
182     elem_list = &(cnst->element_set);
183     xbt_swag_foreach(elem, elem_list) {
184       var = elem->variable;
185       if(var->weight<=0) continue;
186       tmp += var->value;
187     }
188   
189     tmp = tmp - cnst->bound;
190
191     if(tmp > epsilon_min_error){
192       WARN4("The link %s(%p) doesn't match the KKT property, expected less than %e and got %e", (char *)cnst->id, cnst, epsilon_min_error, tmp);
193     }
194   
195   }
196
197   
198   //verify the KKT property of each flow
199   xbt_swag_foreach(var, var_list){
200     if(var->bound <= 0 || var->weight <= 0) continue;
201     tmp = 0;
202     tmp = (var->value - var->bound);
203
204     
205     if(tmp != 0 ||  var->mu != 0){
206       WARN4("The flow %s(%p) doesn't match the KKT property, value expected (=0) got (lambda=%e) (sum_rho=%e)", (char *)var->id, var, var->mu, tmp);
207     }
208
209   }
210
211
212   if(overall_error <= epsilon_min_error){
213     DEBUG1("The method converge in %d iterations.", iteration);
214   }else{
215     WARN1("Method reach %d iterations, which is the maxmimun number of iterations allowed.", iteration);
216   }
217 }
218
219 /*
220  * Returns a double value corresponding to the result of a dicotomi proccess with
221  * respect to a given variable/constraint (\mu in the case of a variable or \lambda in
222  * case of a constraint) and a initial value init. 
223  *
224  * @param init initial value for \mu or \lambda
225  * @param diff a function that computes the differential of with respect a \mu or \lambda
226  * @param var_cnst a pointer to a variable or constraint 
227  * @param min_erro a minimun error tolerated
228  *
229  * @return a double correponding to the result of the dicotomial process
230  */
231 double dicotomi(double init, double diff(double, void*), void *var_cnst, double min_error){
232   double min, max;
233   double overall_error;
234   double middle;
235   double min_diff, max_diff, middle_diff;
236   
237   min = max = init;
238   min_diff = max_diff = middle_diff = 0.0;
239   overall_error = 1;
240
241   if(diff(0.0, var_cnst) > 0){
242     DEBUG1("====> returning 0.0 (diff = %e)", diff(0.0, var_cnst));
243     return 0.0;
244   }
245
246   while(overall_error > min_error){
247     min_diff = diff(min, var_cnst);
248     max_diff = diff(max, var_cnst);
249
250     if( min_diff > 0 && max_diff > 0 ){
251       if(min == max){
252         min = min / 2.0;
253       }else{
254         max = min;
255       }
256     }else if( min_diff < 0 && max_diff < 0 ){
257       if(min == max){
258         max = max * 2.0;
259       }else{
260         min = max;
261       }
262     }else if( min_diff < 0 && max_diff > 0 ){
263       middle = (max + min)/2.0;
264       middle_diff = diff(middle, var_cnst);
265       overall_error = fabs(min - max);
266
267       if( middle_diff < 0 ){
268         min = middle;
269       }else if( middle_diff > 0 ){
270         max = middle;
271       }else{
272         WARN0("Found an optimal solution with 0 error!");
273         overall_error = 0;
274         return middle;
275       }
276
277     }else if(min_diff == 0){
278       return min;
279     }else if(max_diff == 0){
280       return max;
281     }else if(min_diff > 0 && max_diff < 0){
282       WARN0("The impossible happened, partial_diff(min) > 0 && partial_diff(max) < 0");
283     }
284   }
285
286
287   DEBUG1("====> returning %e", (min+max)/2.0);
288   return ((min+max)/2.0);
289 }
290
291 /*
292  *
293  */
294 double partial_diff_mu(double mu, void *param_var){
295   double mu_partial=0.0;
296   lmm_variable_t var = (lmm_variable_t)param_var;
297   int i;
298
299   //for each link with capacity cnsts[i] that uses flow of variable var do
300   for(i=0; i<var->cnsts_number; i++)
301     mu_partial += (var->cnsts[i].constraint)->lambda + mu;
302   
303   mu_partial = (-1.0/mu_partial) + var->bound;
304
305   return mu_partial;
306 }
307
308 /*
309  *
310  */
311 double partial_diff_lambda(double lambda, void *param_cnst){
312
313   double tmp=0.0;
314   int i;
315   xbt_swag_t elem_list = NULL;
316   lmm_element_t elem = NULL;
317   lmm_variable_t var = NULL;
318   lmm_constraint_t cnst= (lmm_constraint_t) param_cnst;
319   double lambda_partial=0.0;
320
321
322   elem_list = &(cnst->element_set);
323
324
325   DEBUG2("Computting diff of cnst (%p) %s", cnst, (char *)cnst->id);
326   
327   xbt_swag_foreach(elem, elem_list) {
328     var = elem->variable;
329     if(var->weight<=0) continue;
330     
331     tmp = 0;
332
333     DEBUG2("===> Variable (%p) %s", var, (char *)var->id);
334
335     for(i=0; i<var->cnsts_number; i++){
336       tmp += (var->cnsts[i].constraint)->lambda;
337       DEBUG1("======> lambda %e + ", (var->cnsts[i].constraint)->lambda);
338     }
339         
340     if(var->bound > 0)
341       tmp += var->mu;
342     
343
344     DEBUG2("======> lambda - %e + %e ", cnst->lambda, lambda);
345
346     tmp = tmp - cnst->lambda + lambda;
347     
348     //avoid a disaster value of lambda
349     //if(tmp==0) tmp = 10e-8;
350     
351     lambda_partial += (-1.0/tmp);
352
353     DEBUG1("======> %e ", (-1.0/tmp));
354   }
355
356   lambda_partial += cnst->bound;
357
358   DEBUG1("===> %e ", lambda_partial);
359
360   return lambda_partial;
361 }
362   
363