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. */
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".
11 #include "xbt/sysdep.h"
12 #include "xbt/mallocator.h"
13 #include "maxmin_private.h"
21 XBT_LOG_NEW_DEFAULT_SUBCATEGORY(surf_lagrange, surf,
22 "Logging specific to SURF (lagrange)");
23 XBT_LOG_NEW_SUBCATEGORY(surf_lagrange_dichotomy, surf,
24 "Logging specific to SURF (lagrange dichotomy)");
27 * Local prototypes to implement the lagrangian optimization with optimal step, also called dichotomy.
29 //solves the proportional fairness using a lagrange optimizition with dichotomy step
30 void lagrange_solve(lmm_system_t sys);
31 //computes the value of the dichotomy using a initial values, init, with a specific variable or constraint
32 double dichotomy(double init, double diff(double, void *), void *var_cnst,
34 //computes the value of the differential of variable param_var applied to mu
35 double partial_diff_mu(double mu, void *param_var);
36 //computes the value of the differential of constraint param_cnst applied to lambda
37 double partial_diff_lambda(double lambda, void *param_cnst);
38 //auxiliar function to compute the partial_diff
39 double diff_aux(lmm_variable_t var, double x);
42 static int __check_kkt(xbt_swag_t cnst_list, xbt_swag_t var_list, int warn)
44 xbt_swag_t elem_list = NULL;
45 lmm_element_t elem = NULL;
46 lmm_constraint_t cnst = NULL;
47 lmm_variable_t var = NULL;
51 //verify the KKT property for each link
52 xbt_swag_foreach(cnst, cnst_list) {
54 elem_list = &(cnst->element_set);
55 xbt_swag_foreach(elem, elem_list) {
62 if (double_positive(tmp - cnst->bound)) {
65 ("The link (%p) is over-used. Expected less than %f and got %f",
66 cnst, cnst->bound, tmp);
69 DEBUG3("Checking KKT for constraint (%p): sat = %f, lambda = %f ",
70 cnst, tmp - cnst->bound, cnst->lambda);
72 /* if(!((fabs(tmp - cnst->bound)<MAXMIN_PRECISION && cnst->lambda>=MAXMIN_PRECISION) || */
73 /* (fabs(tmp - cnst->bound)>=MAXMIN_PRECISION && cnst->lambda<MAXMIN_PRECISION))) { */
74 /* if(warn) WARN1("The KKT condition is not verified for cnst %p...", cnst); */
79 //verify the KKT property of each flow
80 xbt_swag_foreach(var, var_list) {
81 if (var->bound < 0 || var->weight <= 0)
83 DEBUG3("Checking KKT for variable (%p): sat = %f mu = %f", var,
84 var->value - var->bound, var->mu);
86 if (double_positive(var->value - var->bound)) {
89 ("The variable (%p) is too large. Expected less than %f and got %f",
90 var, var->bound, var->value);
94 /* if(!((fabs(var->value - var->bound)<MAXMIN_PRECISION && var->mu>=MAXMIN_PRECISION) || */
95 /* (fabs(var->value - var->bound)>=MAXMIN_PRECISION && var->mu<MAXMIN_PRECISION))) { */
96 /* if(warn) WARN1("The KKT condition is not verified for var %p...",var); */
103 void lagrange_solve(lmm_system_t sys)
106 * Lagrange Variables.
108 int max_iterations = 10000;
109 double epsilon_min_error = 1e-6;
110 double dichotomy_min_error = 1e-8;
111 double overall_error = 1;
114 * Variables to manipulate the data structure proposed to model the maxmin
115 * fairness. See docummentation for more details.
117 xbt_swag_t cnst_list = NULL;
118 lmm_constraint_t cnst = NULL;
120 xbt_swag_t var_list = NULL;
121 lmm_variable_t var = NULL;
124 * Auxiliar variables.
131 DEBUG0("Iterative method configuration snapshot =====>");
132 DEBUG1("#### Maximum number of iterations : %d", max_iterations);
133 DEBUG1("#### Minimum error tolerated : %e",
135 DEBUG1("#### Minimum error tolerated (dichotomy) : %e",
136 dichotomy_min_error);
138 if (!(sys->modified))
142 * Initialize the var list variable with only the active variables.
143 * Associate an index in the swag variables. Initialize mu.
145 var_list = &(sys->variable_set);
147 xbt_swag_foreach(var, var_list) {
148 if ((var->bound < 0.0) || (var->weight <= 0.0)) {
149 DEBUG1("#### NOTE var(%d) is a boundless (or inactive) variable", i);
155 DEBUG3("#### var(%d) %p ->mu : %e", i, var, var->mu);
156 DEBUG3("#### var(%d) %p ->weight: %e", i, var, var->weight);
157 DEBUG3("#### var(%d) %p ->bound: %e", i, var, var->bound);
164 cnst_list = &(sys->active_constraint_set);
165 xbt_swag_foreach(cnst, cnst_list) {
167 cnst->new_lambda = 2.0;
168 DEBUG2("#### cnst(%p)->lambda : %e", cnst, cnst->lambda);
172 * While doesn't reach a minimun error or a number maximum of iterations.
174 while (overall_error > epsilon_min_error && iteration < max_iterations) {
177 DEBUG1("************** ITERATION %d **************", iteration);
180 * Compute the value of mu_i
182 //forall mu_i in mu_1, mu_2, ..., mu_n
183 xbt_swag_foreach(var, var_list) {
184 if ((var->bound >= 0) && (var->weight > 0)) {
185 DEBUG1("====> Working on var (%p)", var);
187 dichotomy(var->mu, partial_diff_mu, var, dichotomy_min_error);
190 DEBUG3("====> var->mu (%p) : %g -> %g", var, var->mu, var->new_mu);
191 var->mu = var->new_mu;
196 * Compute the value of lambda_i
198 //forall lambda_i in lambda_1, lambda_2, ..., lambda_n
199 xbt_swag_foreach(cnst, cnst_list) {
200 DEBUG1("====> Working on cnst (%p)", cnst);
202 dichotomy(cnst->lambda, partial_diff_lambda, cnst,
203 dichotomy_min_error);
204 DEBUG2("====> cnst->lambda (%p) = %e", cnst, cnst->new_lambda);
205 cnst->lambda = cnst->new_lambda;
209 * Now computes the values of each variable (\rho) based on
210 * the values of \lambda and \mu.
213 xbt_swag_foreach(var, var_list) {
214 if (var->weight <= 0)
217 //compute sigma_i + mu_i
219 for (i = 0; i < var->cnsts_number; i++) {
220 tmp += (var->cnsts[i].constraint)->lambda;
224 DEBUG3("\t Working on var (%p). cost = %e; Df = %e", var, tmp,
227 //uses the partial differential inverse function
228 tmp = var->func_fpi(var, tmp);
230 //computes de overall_error using normalized value
231 if (overall_error < (fabs(var->value - tmp) / tmp)) {
232 overall_error = (fabs(var->value - tmp) / tmp);
237 DEBUG3("======> value of var (%p) = %e, overall_error = %e", var,
238 var->value, overall_error);
241 if (!__check_kkt(cnst_list, var_list, 0))
243 DEBUG2("Iteration %d: Overall_error : %f", iteration, overall_error);
247 __check_kkt(cnst_list, var_list, 1);
249 if (overall_error <= epsilon_min_error) {
250 DEBUG1("The method converges in %d iterations.", iteration);
252 if (iteration >= max_iterations) {
254 ("Method reach %d iterations, which is the maximum number of iterations allowed.",
257 /* INFO1("Method converged after %d iterations", iteration); */
259 if (XBT_LOG_ISENABLED(surf_lagrange, xbt_log_priority_debug)) {
265 * Returns a double value corresponding to the result of a dichotomy proccess with
266 * respect to a given variable/constraint (\mu in the case of a variable or \lambda in
267 * case of a constraint) and a initial value init.
269 * @param init initial value for \mu or \lambda
270 * @param diff a function that computes the differential of with respect a \mu or \lambda
271 * @param var_cnst a pointer to a variable or constraint
272 * @param min_erro a minimun error tolerated
274 * @return a double correponding to the result of the dichotomyal process
276 double dichotomy(double init, double diff(double, void *), void *var_cnst,
280 double overall_error;
282 double min_diff, max_diff, middle_diff;
292 min_diff = max_diff = middle_diff = 0.0;
295 if ((diff_0 = diff(1e-16, var_cnst)) >= 0) {
296 CDEBUG1(surf_lagrange_dichotomy, "====> returning 0.0 (diff = %e)",
301 CDEBUG1(surf_lagrange_dichotomy,
302 "====> not detected positive diff in 0 (%e)", diff_0);
304 while (overall_error > min_error) {
306 min_diff = diff(min, var_cnst);
307 max_diff = diff(max, var_cnst);
309 CDEBUG2(surf_lagrange_dichotomy,
310 "DICHOTOMY ===> min = %1.20f , max = %1.20f", min, max);
311 CDEBUG2(surf_lagrange_dichotomy,
312 "DICHOTOMY ===> diffmin = %1.20f , diffmax = %1.20f", min_diff,
315 if (min_diff > 0 && max_diff > 0) {
317 CDEBUG0(surf_lagrange_dichotomy, "Decreasing min");
320 CDEBUG0(surf_lagrange_dichotomy, "Decreasing max");
323 } else if (min_diff < 0 && max_diff < 0) {
325 CDEBUG0(surf_lagrange_dichotomy, "Increasing max");
328 CDEBUG0(surf_lagrange_dichotomy, "Increasing min");
331 } else if (min_diff < 0 && max_diff > 0) {
332 middle = (max + min) / 2.0;
333 middle_diff = diff(middle, var_cnst);
335 if (max != 0.0 && min != 0.0) {
336 overall_error = fabs(min - max) / max;
339 if (middle_diff < 0) {
341 } else if (middle_diff > 0) {
344 CWARN0(surf_lagrange_dichotomy,
345 "Found an optimal solution with 0 error!");
350 } else if (min_diff == 0) {
352 } else if (max_diff == 0) {
354 } else if (min_diff > 0 && max_diff < 0) {
355 CWARN0(surf_lagrange_dichotomy,
356 "The impossible happened, partial_diff(min) > 0 && partial_diff(max) < 0");
358 CWARN2(surf_lagrange_dichotomy,
359 "diffmin (%1.20f) or diffmax (%1.20f) are something I don't know, taking no action.",
367 CDEBUG1(surf_lagrange_dichotomy, "====> returning %e",
369 return ((min + max) / 2.0);
375 double partial_diff_mu(double mu, void *param_var)
377 double mu_partial = 0.0;
378 double sigma_mu = 0.0;
379 lmm_variable_t var = (lmm_variable_t) param_var;
383 for (i = 0; i < var->cnsts_number; i++)
384 sigma_mu += (var->cnsts[i].constraint)->lambda;
386 //compute sigma_i + mu_i
389 //use auxiliar function passing (sigma_i + mu_i)
390 mu_partial = diff_aux(var, sigma_mu);
393 mu_partial += var->bound;
402 double partial_diff_lambda(double lambda, void *param_cnst)
406 xbt_swag_t elem_list = NULL;
407 lmm_element_t elem = NULL;
408 lmm_variable_t var = NULL;
409 lmm_constraint_t cnst = (lmm_constraint_t) param_cnst;
410 double lambda_partial = 0.0;
411 double sigma_i = 0.0;
414 elem_list = &(cnst->element_set);
416 CDEBUG1(surf_lagrange_dichotomy,"Computting diff of cnst (%p)", cnst);
418 xbt_swag_foreach(elem, elem_list) {
419 var = elem->variable;
420 if (var->weight <= 0)
423 //initilize de sumation variable
426 //compute sigma_i of variable var
427 for (i = 0; i < var->cnsts_number; i++) {
428 sigma_i += (var->cnsts[i].constraint)->lambda;
431 //add mu_i if this flow has a RTT constraint associated
435 //replace value of cnst->lambda by the value of parameter lambda
436 sigma_i = (sigma_i - cnst->lambda) + lambda;
438 //use the auxiliar function passing (\sigma_i + \mu_i)
439 lambda_partial += diff_aux(var, sigma_i);
443 lambda_partial += cnst->bound;
446 CDEBUG1(surf_lagrange_dichotomy, "returning = %1.20f", lambda_partial);
449 return lambda_partial;
453 double diff_aux(lmm_variable_t var, double x)
455 double tmp_fpi, result;
457 XBT_IN2("(var (%p), x (%1.20f))", var, x);
458 xbt_assert0(var->func_fp,
459 "Initialize the protocol functions first create variables before.");
461 tmp_fpi = var->func_fpi(var, x);
464 CDEBUG1(surf_lagrange_dichotomy, "returning %1.20f", result);
470 /**************** Vegas and Reno functions *************************/
472 * NOTE for Reno: all functions consider the network
473 * coeficient (alpha) equal to 1.
477 * For Vegas f: $\alpha_f d_f \log\left(x_f\right)$
479 double func_vegas_f(lmm_variable_t var, double x){
480 return var->df * log(x);
484 * For Vegas fp: $\frac{\alpha D_f}{x}$
486 double func_vegas_fp(lmm_variable_t var, double x){
487 //avoid a disaster value - c'est du bricolage mais ca marche
488 /* if(x == 0) x = 10e-8; */
493 * For Vegas fpi: $\frac{\alpha D_f}{x}$
495 double func_vegas_fpi(lmm_variable_t var, double x){
496 //avoid a disaster value - c'est du bricolage mais ca marche
497 /* if(x == 0) x = 10e-8; */
502 * For Vegas fpip: $-\frac{\alpha D_f}{x^2}$
504 double func_vegas_fpip(lmm_variable_t var, double x){
505 //avoid a disaster value - c'est du bricolage mais ca marche
506 /* if(x == 0) x = 10e-8; */
507 return -( var->df/(x*x) ) ;
512 * For Reno f: $\frac{\sqrt{\frac{3}{2}}}{D_f} \arctan\left(\sqrt{\frac{3}{2}}x_f D_f\right)$
514 double func_reno_f(lmm_variable_t var, double x){
515 xbt_assert0(var->df>0.0,"Don't call me with stupid values!");
516 // \sqrt{3/2} = 0.8164965808
517 return (0.8164965808 / var->df) * atan( (0.8164965808 / var->df)*x );
521 * For Reno fp: $\frac{3}{3 {D_f}^2 x^2 + 2}$
523 double func_reno_fp(lmm_variable_t var, double x){
524 return 3 / (3*var->df*var->df*x*x + 2);
528 * For Reno fpi: $\sqrt{\frac{1}{{D_f}^2 x} - \frac{2}{3{D_f}^2}}$
530 double func_reno_fpi(lmm_variable_t var, double x){
533 xbt_assert0(var->df>0.0,"Don't call me with stupid values!");
534 xbt_assert0(x>0.0,"Don't call me with stupid values!");
536 res_fpi = 1/(var->df*var->df*x) - 2/(3*var->df*var->df);
537 if(res_fpi<=0.0) return 0.0;
538 xbt_assert0(res_fpi>0.0,"Don't call me with stupid values!");
539 return sqrt(res_fpi);
543 * For Reno fpip: $-\frac{1}{2 {D_f}^2 x^2\sqrt{\frac{1}{{D_f}^2 x} - \frac{2}{3{D_f}^2}}}$
545 double func_reno_fpip(lmm_variable_t var, double x){
547 double critical_test;
549 xbt_assert0(var->df>0.0,"Don't call me with stupid values!");
550 xbt_assert0(x>0.0,"Don't call me with stupid values!");
552 res_fpip = 1/(var->df*var->df*x) - 2/(3*var->df*var->df);
553 xbt_assert0(res_fpip>0.0,"Don't call me with stupid values!");
554 critical_test = (2*var->df*var->df*x*x*sqrt(res_fpip));
556 return -(1.0/critical_test);