source: coopr.pysp/stable/2.1/coopr/pysp/ph.py @ 2068

Last change on this file since 2068 was 2068, checked in by wehart, 10 years ago

Merged revisions 1952-2067 via svnmerge from
https://software.sandia.gov/svn/public/coopr/coopr.pysp/trunk

........

r1956 | jwatson | 2009-12-02 17:56:53 -0700 (Wed, 02 Dec 2009) | 3 lines


Added --scenario-solver-options and --ef-solver-options options to the "runph" script.

........

r1957 | dlwoodr | 2009-12-03 14:17:35 -0700 (Thu, 03 Dec 2009) | 2 lines


Documentation updates for pysp

........

r1974 | wehart | 2009-12-06 17:20:56 -0700 (Sun, 06 Dec 2009) | 2 lines


Updating PyPI categories

........

r1978 | jwatson | 2009-12-10 21:29:33 -0700 (Thu, 10 Dec 2009) | 3 lines


Eliminated exception-handling logic when loading user-defined extension modules in PH. The range of exceptions is too large, and for debugging purposes, it is more useful to see the raw trace output.

........

r1979 | jwatson | 2009-12-10 22:23:17 -0700 (Thu, 10 Dec 2009) | 5 lines


Biggest enhancement: The efwriter command-line script now has the option to output a CVaR-weighted objective term. Not extensively tested, but behaves sane on a number of small test cases.


Improved exception handling and error diagnostics in both the runph and efwriter scripts.

........

r1985 | jwatson | 2009-12-12 10:45:17 -0700 (Sat, 12 Dec 2009) | 3 lines


Modified PH to only use warm-starts if a solver has the capability!

........

r1998 | jwatson | 2009-12-13 15:17:58 -0700 (Sun, 13 Dec 2009) | 3 lines


Changed references to _component to active_component.

........

r2026 | wehart | 2009-12-21 23:27:06 -0700 (Mon, 21 Dec 2009) | 2 lines


Attempting to update PH. I'm not sure if this works, since I don't know how to test PH.

........

r2029 | jwatson | 2009-12-22 09:52:21 -0700 (Tue, 22 Dec 2009) | 3 lines


Some fixes to the ef writer based on Bill's recent changes to _initialize_constraint.

........

r2035 | jwatson | 2009-12-22 21:10:32 -0700 (Tue, 22 Dec 2009) | 3 lines


Added --scenario-mipgap option to PH script. Added _mipgap attribute to PH object. This attribute is mirrored to the solver plugin at the initiation of each iteration, after any PH extensions have the opportunity to provide a new value to the attribute. It is currently made use of by the WW PH extension.

........

r2037 | dlwoodr | 2009-12-23 14:38:43 -0700 (Wed, 23 Dec 2009) | 2 lines


add this example from Fernando

........

r2038 | dlwoodr | 2009-12-23 14:46:56 -0700 (Wed, 23 Dec 2009) | 3 lines


finish the job: we can now duplicate Fernando's example

........

r2039 | jwatson | 2009-12-23 15:13:54 -0700 (Wed, 23 Dec 2009) | 3 lines


Missed fix with new constraint initialization syntax in PH linearization.

........

r2058 | jwatson | 2009-12-29 10:57:58 -0700 (Tue, 29 Dec 2009) | 3 lines


Missed some _initialize_constraint function calls within the PySP EF writer during the recent switch to the corresponding "add" method.

........

r2059 | jwatson | 2009-12-29 10:58:34 -0700 (Tue, 29 Dec 2009) | 3 lines


Enabling garbage collection by default in PH.

........

r2060 | jwatson | 2009-12-29 10:59:05 -0700 (Tue, 29 Dec 2009) | 3 lines


Elimnating some debug output.

........

r2061 | jwatson | 2009-12-29 11:07:47 -0700 (Tue, 29 Dec 2009) | 3 lines


Fixing some option documentation in PH.

........

r2062 | jwatson | 2009-12-29 12:00:37 -0700 (Tue, 29 Dec 2009) | 3 lines


Added ef-mipgap option to PH scripts.

........

File size: 80.2 KB
Line 
1#  _________________________________________________________________________
2#
3#  Coopr: A COmmon Optimization Python Repository
4#  Copyright (c) 2008 Sandia Corporation.
5#  This software is distributed under the BSD License.
6#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
7#  the U.S. Government retains certain rights in this software.
8#  For more information, see the Coopr README.txt file.
9#  _________________________________________________________________________
10
11import sys
12import types
13from coopr.pyomo import *
14import copy
15import os.path
16import traceback
17import copy
18from coopr.opt import SolverResults,SolverStatus
19from coopr.opt.base import SolverFactory
20from coopr.opt.parallel import SolverManagerFactory
21import time
22import types
23import pickle
24from math import fabs, log, exp
25
26from scenariotree import *
27from phutils import *
28
29from pyutilib.plugin.core import ExtensionPoint
30
31from coopr.pysp.phextension import IPHExtension
32
33class ProgressiveHedging(object):
34
35   #
36   # routine to compute linearization breakpoints uniformly between the bounds and the mean.
37   #
38   
39   def compute_uniform_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints_per_side):
40     
41      breakpoints = []
42
43      # add the lower bound - the first breakpoint.
44      breakpoints.append(lb)
45
46      # determine the breakpoints to the left of the mean.
47      left_step = (xavg - lb) / num_breakpoints_per_side
48      current_x = lb
49      for i in range(1,num_breakpoints_per_side+1):
50         this_lb = current_x
51         this_ub = current_x+left_step
52         if (fabs(this_lb-lb) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
53            breakpoints.append(this_lb)
54         current_x += left_step           
55
56      # add the mean - it's always a breakpoint. unless!
57      # the lb or ub and the avg are the same.
58      if (fabs(lb-xavg) > self._integer_tolerance) and (fabs(ub-xavg) > self._integer_tolerance):
59         breakpoints.append(xavg)
60
61      # determine the breakpoints to the right of the mean.
62      right_step = (ub - xavg) / num_breakpoints_per_side
63      current_x = xavg
64      for i in range(1,num_breakpoints_per_side+1):
65         this_lb = current_x
66         this_ub = current_x+right_step
67         if (fabs(this_ub-xavg) > self._integer_tolerance) and (fabs(this_ub-ub) > self._integer_tolerance):         
68            breakpoints.append(this_ub)
69         current_x += right_step
70
71      # add the upper bound - the last breakpoint.
72      # the upper bound should always be different than the lower bound (I say with some
73      # hesitation - it really depends on what plugins are doing to modify the bounds dynamically).
74      breakpoints.append(ub)
75
76      return breakpoints
77
78   #
79   # routine to compute linearization breakpoints uniformly between the current node min/max bounds.
80   #   
81
82   def compute_uniform_between_nodestat_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints):
83
84      breakpoints = []
85
86      # add the lower bound - the first breakpoint.
87      breakpoints.append(lb)
88
89      # add the node-min - the second breakpoint. but only if it is different than the lower bound and the mean.
90      if (fabs(node_min-lb) > self._integer_tolerance) and (fabs(node_min-xavg) > self._integer_tolerance):     
91         breakpoints.append(node_min)
92
93      step = (node_max - node_min) / num_breakpoints
94      current_x = node_min
95      for i in range(1,num_breakpoints+1):
96         this_lb = current_x
97         this_ub = current_x+step
98         if (fabs(this_lb-node_min) > self._integer_tolerance) and (fabs(this_lb-node_max) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
99            breakpoints.append(this_lb)
100         current_x += step           
101
102      # add the node-max - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
103      if (fabs(node_max-ub) > self._integer_tolerance) and (fabs(node_max-xavg) > self._integer_tolerance):           
104         breakpoints.append(node_max)     
105
106      # add the upper bound - the last breakpoint.
107      breakpoints.append(ub)
108
109      # add the mean - it's always a breakpoint. unless!
110      # it happens to be equal to (within tolerance) the lower or upper bounds.
111      # sort to insert it in the correct place.
112      if (fabs(xavg - lb) > self._integer_tolerance) and (fabs(xavg - ub) > self._integer_tolerance):
113         breakpoints.append(xavg)
114      breakpoints.sort()
115
116      return breakpoints
117
118   #
119   # routine to compute linearization breakpoints using "Woodruff" relaxation of the compute_uniform_between_nodestat_breakpoints.
120   #   
121
122   def compute_uniform_between_woodruff_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints):
123
124      breakpoints = []
125
126      # add the lower bound - the first breakpoint.
127      breakpoints.append(lb)
128
129      # be either three "differences" from the mean, or else "halfway to the bound", whichever is closer to the mean.
130      left = max(xavg - 3.0 * (xavg - node_min), xavg - 0.5 * (xavg - lb))
131      right = min(xavg + 3.0 * (node_max - xavg), xavg + 0.5 * (ub - xavg))     
132
133      # add the left bound - the second breakpoint. but only if it is different than the lower bound and the mean.
134      if (fabs(left-lb) > self._integer_tolerance) and (fabs(left-xavg) > self._integer_tolerance):     
135         breakpoints.append(left)
136
137      step = (right - left) / num_breakpoints
138      current_x = left
139      for i in range(1,num_breakpoints+1):
140         this_lb = current_x
141         this_ub = current_x+step
142         if (fabs(this_lb-left) > self._integer_tolerance) and (fabs(this_lb-right) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
143            breakpoints.append(this_lb)
144         current_x += step           
145
146      # add the right bound - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
147      if (fabs(right-ub) > self._integer_tolerance) and (fabs(right-xavg) > self._integer_tolerance):           
148         breakpoints.append(right)     
149
150      # add the upper bound - the last breakpoint.
151      breakpoints.append(ub)
152
153      # add the mean - it's always a breakpoint.
154      # sort to insert it in the correct place.
155      breakpoints.append(xavg)
156      breakpoints.sort()
157
158      return breakpoints
159
160   #
161   # routine to compute linearization breakpoints based on an exponential distribution from the mean in each direction.
162   #   
163
164   def compute_exponential_from_mean_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints_per_side):
165
166      breakpoints = []
167
168      # add the lower bound - the first breakpoint.
169      breakpoints.append(lb)
170
171      # determine the breakpoints to the left of the mean.
172      left_delta = xavg - lb
173      base = exp(log(left_delta) / num_breakpoints_per_side)
174      current_offset = base
175      for i in range(1,num_breakpoints_per_side+1):
176         current_x = xavg - current_offset
177         if (fabs(current_x-lb) > self._integer_tolerance) and (fabs(current_x-xavg) > self._integer_tolerance):
178            breakpoints.append(current_x)
179         current_offset *= base
180
181      # add the mean - it's always a breakpoint.
182      breakpoints.append(xavg)
183
184      # determine the breakpoints to the right of the mean.
185      right_delta = ub - xavg
186      base = exp(log(right_delta) / num_breakpoints_per_side)
187      current_offset = base
188      for i in range(1,num_breakpoints_per_side+1):
189         current_x = xavg + current_offset
190         if (fabs(current_x-xavg) > self._integer_tolerance) and (fabs(current_x-ub) > self._integer_tolerance):         
191            breakpoints.append(current_x)
192         current_offset *= base
193
194      # add the upper bound - the last breakpoint.
195      breakpoints.append(ub)         
196
197      return breakpoints     
198
199   #
200   # a utility intended for folks who are brave enough to script rho setting in a python file.
201   #
202
203   def setRhoAllScenarios(self, variable_value, rho_expression):
204
205      variable_name = None
206      variable_index = None
207
208      if isVariableNameIndexed(variable_value.name) is True:
209
210         variable_name, variable_index = extractVariableNameAndIndex(variable_value.name)
211
212      else:
213
214         variable_name = variable_value.name
215         variable_index = None
216     
217      new_rho_value = rho_expression()
218
219      if self._verbose is True:
220         print "Setting rho="+str(new_rho_value)+" for variable="+variable_value.name
221
222      for instance_name, instance in self._instances.items():
223
224         rho_param = getattr(instance, "PHRHO_"+variable_name)
225         rho_param[variable_index] = new_rho_value
226
227   #
228   # a utility intended for folks who are brave enough to script variable bounds setting in a python file.
229   #
230
231   def setVariableBoundsAllScenarios(self, variable_name, variable_index, lower_bound, upper_bound):
232
233      if isinstance(lower_bound, float) is False:
234         raise ValueError, "Lower bound supplied to PH method setVariableBoundsAllScenarios for variable="+variable_name+indexToString(variable_index)+" must be a constant; value supplied="+str(lower_bound)
235
236      if isinstance(upper_bound, float) is False:
237         raise ValueError, "Upper bound supplied to PH method setVariableBoundsAllScenarios for variable="+variable_name+indexToString(variable_index)+" must be a constant; value supplied="+str(upper_bound)
238
239      for instance_name, instance in self._instances.items():
240
241         variable = getattr(instance, variable_name)
242         variable[variable_index].setlb(lower_bound)
243         variable[variable_index].setub(upper_bound)
244
245   #
246   # a utility intended for folks who are brave enough to script variable bounds setting in a python file.
247   # same functionality as above, but applied to all indicies of the variable, in all scenarios.
248   #
249
250   def setVariableBoundsAllIndicesAllScenarios(self, variable_name, lower_bound, upper_bound):
251
252      if isinstance(lower_bound, float) is False:
253         raise ValueError, "Lower bound supplied to PH method setVariableBoundsAllIndiciesAllScenarios for variable="+variable_name+" must be a constant; value supplied="+str(lower_bound)
254
255      if isinstance(upper_bound, float) is False:
256         raise ValueError, "Upper bound supplied to PH method setVariableBoundsAllIndicesAllScenarios for variable="+variable_name+" must be a constant; value supplied="+str(upper_bound)
257
258      for instance_name, instance in self._instances.items():
259
260         variable = getattr(instance, variable_name)
261         for index in variable:
262            variable[index].setlb(lower_bound)
263            variable[index].setub(upper_bound)
264
265   #
266   # checkpoint the current PH state via pickle'ing. the input iteration count
267   # simply serves as a tag to create the output file name. everything with the
268   # exception of the _ph_plugins, _solver_manager, and _solver attributes are
269   # pickled. currently, plugins fail in the pickle process, which is fine as
270   # JPW doesn't think you want to pickle plugins (particularly the solver and
271   # solver manager) anyway. For example, you might want to change those later,
272   # after restoration - and the PH state is independent of how scenario
273   # sub-problems are solved.
274   #
275
276   def checkpoint(self, iteration_count):
277
278      checkpoint_filename = "checkpoint."+str(iteration_count)
279
280      tmp_ph_plugins = self._ph_plugins
281      tmp_solver_manager = self._solver_manager
282      tmp_solver = self._solver
283
284      self._ph_plugins = None
285      self._solver_manager = None
286      self._solver = None
287     
288      checkpoint_file = open(checkpoint_filename, "w")
289      pickle.dump(self,checkpoint_file)
290      checkpoint_file.close()
291
292      self._ph_plugins = tmp_ph_plugins
293      self._solver_manager = tmp_solver_manager
294      self._solver = tmp_solver
295
296      print "Checkpoint written to file="+checkpoint_filename
297   
298   #
299   # a simple utility to count the number of continuous and discrete variables in a set of instances.
300   # unused variables are ignored, and counts include all active indices. returns a pair - num-discrete,
301   # num-continuous.
302   #
303
304   def compute_variable_counts(self):
305
306      num_continuous_vars = 0
307      num_discrete_vars = 0
308     
309      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage
310
311         for tree_node in stage._tree_nodes:
312
313            for (variable, index_template, variable_indices) in stage._variables:
314
315               for index in variable_indices:
316
317                  is_used = True # until proven otherwise                     
318                  for scenario in tree_node._scenarios:
319                     instance = self._instances[scenario._name]
320                     if getattr(instance,variable.name)[index].status == VarStatus.unused:
321                        is_used = False
322
323                  if is_used is True:                       
324
325                     if isinstance(variable.domain, IntegerSet) or isinstance(variable.domain, BooleanSet):
326                        num_discrete_vars = num_discrete_vars + 1
327                     else:
328                        num_continuous_vars = num_continuous_vars + 1
329
330      return (num_discrete_vars, num_continuous_vars)
331
332   #
333   # ditto above, but count the number of fixed discrete and continuous variables.
334   # important: once a variable (value) is fixed, it is flagged as unused in the
335   # course of presolve - because it is no longer referenced. this makes sense,
336   # of course; it's just something to watch for. this is an obvious assumption
337   # that we won't be fixing unused variables, which should not be an issue.
338   #
339
340   def compute_fixed_variable_counts(self):
341
342      num_fixed_continuous_vars = 0
343      num_fixed_discrete_vars = 0
344     
345      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage
346         
347         for tree_node in stage._tree_nodes:
348
349            for (variable, index_template, variable_indices) in stage._variables:
350
351               for index in variable_indices:
352
353                  # implicit assumption is that if a variable value is fixed in one
354                  # scenario, it is fixed in all scenarios.
355
356                  is_fixed = False # until proven otherwise
357                  for scenario in tree_node._scenarios:
358                     instance = self._instances[scenario._name]
359                     var_value = getattr(instance,variable.name)[index]
360                     if var_value.fixed is True:
361                        is_fixed = True
362
363                  if is_fixed is True:
364
365                     if isinstance(variable.domain, IntegerSet) or isinstance(variable.domain, BooleanSet):
366                        num_fixed_discrete_vars = num_fixed_discrete_vars + 1
367                     else:
368                        num_fixed_continuous_vars = num_fixed_continuous_vars + 1                           
369
370      return (num_fixed_discrete_vars, num_fixed_continuous_vars)
371
372   # a utility to create piece-wise linear constraint expressions for a given variable, for
373   # use in constructing the augmented (penalized) PH objective. lb and ub are the bounds
374   # on this piece, variable is the actual instance variable, and average is the instance
375   # parameter specifying the average of this variable across instances sharing passing
376   # through a common tree node. lb and ub are floats.
377   def _create_piecewise_constraint_expression(self, lb, ub, instance_variable, variable_average, quad_variable):
378
379      penalty_at_lb = (lb - variable_average()) * (lb - variable_average())
380      penalty_at_ub = (ub - variable_average()) * (ub - variable_average())
381      slope = (penalty_at_ub - penalty_at_lb) / (ub - lb)
382      intercept = penalty_at_lb - slope * lb
383      expression = (0.0, quad_variable - slope * instance_variable - intercept, None)
384
385      return expression
386
387   # when the quadratic penalty terms are approximated via piecewise linear segments,
388   # we end up (necessarily) "littering" the scenario instances with extra constraints.
389   # these need to and should be cleaned up after PH, for purposes of post-PH manipulation,
390   # e.g., writing the extensive form.
391   def _cleanup_scenario_instances(self):
392
393      for instance_name, instance in self._instances.items():
394
395         for constraint_name in self._instance_augmented_attributes[instance_name]:
396
397            instance._clear_attribute(constraint_name)
398
399         # if you don't pre-solve, the name collections won't be updated.
400         instance.presolve()
401
402   # create PH weight and xbar vectors, on a per-scenario basis, for each variable that is not in the
403   # final stage, i.e., for all variables that are being blended by PH. the parameters are created
404   # in the space of each scenario instance, so that they can be directly and automatically
405   # incorporated into the (appropriately modified) objective function.
406
407   def _create_ph_scenario_parameters(self):
408
409      for (instance_name, instance) in self._instances.items():
410
411         # first, gather all unique variables referenced in any stage
412         # other than the last, independent of specific indices. this
413         # "gather" step is currently required because we're being lazy
414         # in terms of index management in the scenario tree - which
415         # should really be done in terms of sets of indices.
416         # NOTE: technically, the "instance variables" aren't really references
417         # to the variable in the instance - instead, the reference model. this
418         # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
419         instance_variables = {}
420         
421         for stage in self._scenario_tree._stages[:-1]:
422           
423            for (reference_variable, index_template, reference_indices) in stage._variables:
424                 
425               if reference_variable.name not in instance_variables.keys():
426                     
427                  instance_variables[reference_variable.name] = reference_variable
428
429         # for each blended variable, create a corresponding ph weight and average parameter in the instance.
430         # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor
431         # in the grand scheme of things.
432
433         for (variable_name, reference_variable) in instance_variables.items():
434
435            # PH WEIGHT
436
437            new_w_index = reference_variable._index
438            new_w_parameter_name = "PHWEIGHT_"+reference_variable.name
439            new_w_parameter = Param(new_w_index,name=new_w_parameter_name)
440            setattr(instance,new_w_parameter_name,new_w_parameter)
441
442            # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference
443            # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not
444            # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the
445            # issue in any case for now.
446            for index in new_w_index:
447               new_w_parameter[index] = 0.0
448
449            # PH AVG
450
451            new_avg_index = reference_variable._index
452            new_avg_parameter_name = "PHAVG_"+reference_variable.name
453            new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
454            setattr(instance,new_avg_parameter_name,new_avg_parameter)
455
456            for index in new_avg_index:
457               new_avg_parameter[index] = 0.0
458
459            # PH RHO
460
461            new_rho_index = reference_variable._index
462            new_rho_parameter_name = "PHRHO_"+reference_variable.name
463            new_rho_parameter = Param(new_rho_index,name=new_rho_parameter_name)
464            setattr(instance,new_rho_parameter_name,new_rho_parameter)
465
466            for index in new_rho_index:
467               new_rho_parameter[index] = self._rho
468
469            # PH LINEARIZED PENALTY TERM
470
471            if self._linearize_nonbinary_penalty_terms > 0:
472
473               new_penalty_term_variable_index = reference_variable._index
474               new_penalty_term_variable_name = "PHQUADPENALTY_"+reference_variable.name
475               # this is a quadratic penalty term - the lower bound is 0!
476               new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None))
477               new_penalty_term_variable.construct()
478               setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable)
479               self._instance_augmented_attributes[instance_name].append(new_penalty_term_variable_name)
480
481            # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY.
482
483            # also controls whether weight updates proceed at any iteration.
484
485            new_blend_index = reference_variable._index
486            new_blend_parameter_name = "PHBLEND_"+reference_variable.name
487            new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary)
488            setattr(instance,new_blend_parameter_name,new_blend_parameter)
489
490            for index in new_blend_index:
491               new_blend_parameter[index] = 1
492     
493
494   """ Constructor
495       Arguments:
496          max_iterations        the maximum number of iterations to run PH (>= 0). defaults to 0.
497          rho                   the global rho value (> 0). defaults to 0.
498          rho_setter            an optional name of a python file used to set particular variable rho values.
499          solver                the solver type that PH uses to solve scenario sub-problems. defaults to "cplex".
500          solver_manager        the solver manager type that coordinates scenario sub-problem solves. defaults to "serial".
501          keep_solver_files     do I keep intermediate solver files around (for debugging)? defaults to False.
502          output_solver_log     do I dump the solver log (as it is being generated) to the screen? defaults to False.
503          output_solver_results do I output (for debugging) the detailed solver results, including solutions, for scenario solves? defaults to False.
504          verbose               does the PH object stream debug/status output? defaults to False.
505          output_times          do I output timing statistics? defaults to False (e.g., useful in the case where you want to regression test against baseline output).
506          checkpoint_interval   how many iterations between writing a checkpoint file containing the entire PH state? defaults to 0, indicating never.
507
508   """
509   def __init__(self, *args, **kwds):
510
511      # PH configuration parameters
512      self._rho = 0.0 # a default, global values for rhos.
513      self._rho_setter = None # filename for the modeler to set rho on a per-variable or per-scenario basis.
514      self._bounds_setter = None # filename for the modeler to set rho on a per-variable basis, after all scenarios are available.
515      self._max_iterations = 0
516
517      # PH reporting parameters
518      self._verbose = False # do I flood the screen with status output?
519      self._output_continuous_variable_stats = True # when in verbose mode, do I output weights/averages for continuous variables?
520      self._output_solver_results = False
521      self._output_times = False
522
523      # PH run-time variables
524      self._current_iteration = 0 # the 'k'
525      self._xbar = {} # current per-variable averages. maps (node_id, variable_name) -> value
526      self._initialized = False # am I ready to call "solve"? Set to True by the initialize() method.
527
528      # PH solver information / objects.
529      self._solver_type = "cplex"
530      self._solver_manager_type = "serial" # serial or pyro are the options currently available
531     
532      self._solver = None # will eventually be unnecessary once Bill eliminates the need for a solver in the solver manager constructor.
533      self._solver_manager = None 
534     
535      self._keep_solver_files = False
536      self._output_solver_log = False
537
538      # PH convergence computer/updater.
539      self._converger = None
540
541      # PH history
542      self._solutions = {}
543
544      # the checkpoint interval - expensive operation, but worth it for big models.
545      # 0 indicates don't checkpoint.
546      self._checkpoint_interval = 0
547
548      # all information related to the scenario tree (implicit and explicit).
549      self._model = None # not instantiated
550      self._model_instance = None # instantiated
551
552      self._scenario_tree = None
553     
554      self._scenario_data_directory = "" # this the prefix for all scenario data
555      self._instances = {} # maps scenario name to the corresponding model instance
556      self._instance_augmented_attributes = {} # maps scenario name to a list of the constraints added (e.g., for piecewise linear approximation) to the instance by PH.
557
558      # for various reasons (mainly hacks at this point), it's good to know whether we're minimizing or maximizing.
559      self._is_minimizing = None
560
561      # global handle to ph extension plugins
562      self._ph_plugins = ExtensionPoint(IPHExtension)
563
564      # PH timing statistics - relative to last invocation.
565      self._init_start_time = None # for initialization() method
566      self._init_end_time = None
567      self._solve_start_time = None # for solve() method
568      self._solve_end_time = None
569      self._cumulative_solve_time = None # seconds, over course of solve()
570      self._cumulative_xbar_time = None # seconds, over course of update_xbars()
571      self._cumulative_weight_time = None # seconds, over course of update_weights()
572
573      # do I disable warm-start for scenario sub-problem solves during PH iterations >= 1?
574      self._disable_warmstarts = False
575
576      # do I drop proximal (quadratic penalty) terms from the weighted objective functions?
577      self._drop_proximal_terms = False
578
579      # do I linearize the quadratic penalty term for continuous variables via a
580      # piecewise linear approximation? the default should always be 0 (off), as the
581      # user should be aware when they are forcing an approximation.
582      self._linearize_nonbinary_penalty_terms = 0
583
584      # the breakpoint distribution strategy employed when linearizing. 0 implies uniform
585      # distribution between the variable lower and upper bounds.
586      self._breakpoint_strategy = 0
587
588      # do I retain quadratic objective terms associated with binary variables? in general,
589      # there is no good reason to not linearize, but just in case, I introduced the option.
590      self._retain_quadratic_binary_terms = False
591
592      # PH default tolerances - for use in fixing and testing equality across scenarios,
593      # and other stuff.
594      self._integer_tolerance = 0.00001
595
596      # PH maintains a mipgap that is applied to each scenario solve that is performed.
597      # this attribute can be changed by PH extensions, and the change will be applied
598      # on all subsequent solves - until it is modified again. the default is None,
599      # indicating unassigned.
600      self._mipgap = None
601
602      # we only store these temporarily...
603      scenario_solver_options = None
604
605      # process the keyword options
606      for key in kwds.keys():
607         if key == "max_iterations":
608            self._max_iterations = kwds[key]
609         elif key == "rho":
610            self._rho = kwds[key]
611         elif key == "rho_setter":
612            self._rho_setter = kwds[key]
613         elif key == "bounds_setter":
614            self._bounds_setter = kwds[key]                       
615         elif key == "solver":
616            self._solver_type = kwds[key]
617         elif key == "solver_manager":
618            self._solver_manager_type = kwds[key]
619         elif key == "scenario_solver_options":
620            scenario_solver_options = kwds[key]
621         elif key == "scenario_mipgap":
622            self._mipgap = kwds[key]           
623         elif key == "keep_solver_files":
624            self._keep_solver_files = kwds[key]
625         elif key == "output_solver_results":
626            self._output_solver_results = kwds[key]
627         elif key == "output_solver_log":
628            self._output_solver_log = kwds[key]
629         elif key == "verbose":
630            self._verbose = kwds[key]
631         elif key == "output_times":
632            self._output_times = kwds[key]
633         elif key == "disable_warmstarts":
634            self._disable_warmstarts = kwds[key]
635         elif key == "drop_proximal_terms":
636            self._drop_proximal_terms = kwds[key]
637         elif key == "retain_quadratic_binary_terms":
638            self._retain_quadratic_binary_terms = kwds[key]
639         elif key == "linearize_nonbinary_penalty_terms":
640            self._linearize_nonbinary_penalty_terms = kwds[key]
641         elif key == "breakpoint_strategy":
642            self._breakpoint_strategy = kwds[key]
643         elif key == "checkpoint_interval":
644            self._checkpoint_interval = kwds[key]           
645         else:
646            print "Unknown option=" + key + " specified in call to PH constructor"
647
648      # validate all "atomic" options (those that can be validated independently)
649      if self._max_iterations < 0:
650         raise ValueError, "Maximum number of PH iterations must be non-negative; value specified=" + `self._max_iterations`
651      if self._rho <= 0.0:
652         raise ValueError, "Value of the rho parameter in PH must be non-zero positive; value specified=" + `self._rho`
653      if (self._mipgap is not None) and ((self._mipgap < 0.0) or (self._mipgap > 1.0)):
654         raise ValueError, "Value of the mipgap parameter in PH must be on the unit interval; value specified=" + `self._mipgap`
655
656      # validate the linearization (number of pieces) and breakpoint distribution parameters.
657      if self._linearize_nonbinary_penalty_terms < 0:
658         raise ValueError, "Value of linearization parameter for nonbinary penalty terms must be non-negative; value specified=" + `self._linearize_nonbinary_penalty_terms`         
659      if self._breakpoint_strategy < 0:
660         raise ValueError, "Value of the breakpoint distribution strategy parameter must be non-negative; value specified=" + str(self._breakpoint_strategy)
661      if self._breakpoint_strategy > 3:
662         raise ValueError, "Unknown breakpoint distribution strategy specified - valid values are between 0 and 2, inclusive; value specified=" + str(self._breakpoint_strategy)
663   
664      # validate rho setter file if specified.
665      if self._rho_setter is not None:
666         if os.path.exists(self._rho_setter) is False:
667            raise ValueError, "The rho setter script file="+self._rho_setter+" does not exist"
668
669      # validate bounds setter file if specified.
670      if self._bounds_setter is not None:
671         if os.path.exists(self._bounds_setter) is False:
672            raise ValueError, "The bounds setter script file="+self._bounds_setter+" does not exist"     
673
674      # validate the checkpoint interval.
675      if self._checkpoint_interval < 0:
676         raise ValueError, "A negative checkpoint interval with value="+str(self._checkpoint_interval)+" was specified in call to PH constructor"
677
678      # construct the sub-problem solver.
679      if self._verbose is True:
680         print "Constructing solver type="+self._solver_type         
681      self._solver = SolverFactory(self._solver_type)
682      if self._solver == None:
683         raise ValueError, "Unknown solver type=" + self._solver_type + " specified in call to PH constructor"
684      if self._keep_solver_files is True:
685         self._solver.keepFiles = True
686      if len(scenario_solver_options) > 0:
687         if self._verbose is True:
688            print "Initializing scenario sub-problem solver with options="+str(scenario_solver_options)
689         self._solver.set_options("".join(scenario_solver_options))
690
691      # construct the solver manager.
692      if self._verbose is True:
693         print "Constructing solver manager of type="+self._solver_manager_type
694      self._solver_manager = SolverManagerFactory(self._solver_manager_type)
695      if self._solver_manager is None:
696         raise ValueError, "Failed to create solver manager of type="+self._solver_manager_type+" specified in call to PH constructor"
697
698      # a set of all valid PH iteration indicies is generally useful for plug-ins, so create it here.
699      self._iteration_index_set = Set(name="PHIterations")
700      for i in range(0,self._max_iterations + 1):
701         self._iteration_index_set.add(i)     
702
703      # spit out parameterization if verbosity is enabled
704      if self._verbose is True:
705         print "PH solver configuration: "
706         print "   Max iterations=" + `self._max_iterations`
707         print "   Default global rho=" + `self._rho`
708         if self._rho_setter is not None:
709            print "   Rho initialization file=" + self._rho_setter
710         if self._bounds_setter is not None:
711            print "   Variable bounds initialization file=" + self._bounds_setter           
712         print "   Sub-problem solver type=" + `self._solver_type`
713         print "   Solver manager type=" + `self._solver_manager_type`
714         print "   Keep solver files? " + str(self._keep_solver_files)
715         print "   Output solver results? " + str(self._output_solver_results)
716         print "   Output solver log? " + str(self._output_solver_log)
717         print "   Output times? " + str(self._output_times)
718         print "   Checkpoint interval="+str(self._checkpoint_interval)
719
720   """ Initialize PH with model and scenario data, in preparation for solve().
721       Constructs and reads instances.
722   """
723   def initialize(self, scenario_data_directory_name=".", model=None, model_instance=None, scenario_tree=None, converger=None):
724
725      self._init_start_time = time.time()
726
727      if self._verbose is True:
728         print "Initializing PH"
729         print "   Scenario data directory=" + scenario_data_directory_name
730
731      if not os.path.exists(scenario_data_directory_name):
732         raise ValueError, "Scenario data directory=" + scenario_data_directory_name + " either does not exist or cannot be read"
733
734      self._scenario_data_directory_name = scenario_data_directory_name
735
736      # IMPT: The input model should be an *instance*, as it is very useful (critical!) to know
737      #       the dimensions of sets, be able to store suffixes on variable values, etc.
738      if model is None:
739         raise ValueError, "A model must be supplied to the PH initialize() method"
740
741      if scenario_tree is None:
742         raise ValueError, "A scenario tree must be supplied to the PH initialize() method"
743
744      if converger is None:
745         raise ValueError, "A convergence computer must be supplied to the PH initialize() method"
746
747      self._model = model
748      self._model_instance = model_instance
749      self._scenario_tree = scenario_tree
750      self._converger = converger
751
752      model_objective = model.active_components(Objective)
753      self._is_minimizing = (model_objective[ model_objective.keys()[0] ].sense == minimize)
754
755      self._converger.reset()
756
757      # construct the instances for each scenario
758      if self._verbose is True:
759         if self._scenario_tree._scenario_based_data == 1:
760            print "Scenario-based instance initialization enabled"
761         else:
762            print "Node-based instance initialization enabled"
763         
764      for scenario in self._scenario_tree._scenarios:
765
766         scenario_instance = None
767
768         if self._verbose is True:
769            print "Creating instance for scenario=" + scenario._name
770
771         try:
772            if self._scenario_tree._scenario_based_data == 1:
773               scenario_data_filename = self._scenario_data_directory_name + os.sep + scenario._name + ".dat"
774               if self._verbose is True:
775                  print "Data for scenario=" + scenario._name + " loads from file=" + scenario_data_filename
776               scenario_instance = (self._model).create(scenario_data_filename)
777            else:
778               scenario_instance = self._model.clone()
779               scenario_data = ModelData()
780               current_node = scenario._leaf_node
781               while current_node is not None:
782                  node_data_filename = self._scenario_data_directory_name + os.sep + current_node._name + ".dat"
783                  if self._verbose is True:
784                     print "Node data for scenario=" + scenario._name + " partially loading from file=" + node_data_filename
785                  scenario_data.add_data_file(node_data_filename)
786                  current_node = current_node._parent
787               scenario_data.read(model=scenario_instance)
788               scenario_instance._load_model_data(scenario_data)
789               scenario_instance.presolve()
790         except:
791            print "Encountered exception in model instance creation - traceback:"
792            traceback.print_exc()
793            raise RuntimeError, "Failed to create model instance for scenario=" + scenario._name
794
795         self._instances[scenario._name] = scenario_instance
796         self._instances[scenario._name].name = scenario._name
797         self._instance_augmented_attributes[scenario._name] = []
798
799      # create ph-specific parameters (weights, xbar, etc.) for each instance.
800
801      if self._verbose is True:
802         print "Creating weight, average, and rho parameter vectors for scenario instances"
803
804      self._create_ph_scenario_parameters()
805           
806      # if specified, run the user script to initialize variable rhos at their whim.
807      if self._rho_setter is not None:
808         print "Executing user rho set script from filename=", self._rho_setter
809         execfile(self._rho_setter)
810
811      # with the instances created, run the user script to initialize variable bounds.
812      if self._bounds_setter is not None:
813         print "Executing user variable bounds set script from filename=", self._bounds_setter
814         execfile(self._bounds_setter)
815
816      # create parameters to store variable statistics (of general utility) at each node in the scenario tree.
817
818      if self._verbose is True:
819         print "Creating variable statistic (min/avg/max) parameter vectors for scenario tree nodes"
820
821      for stage in self._scenario_tree._stages[:-1]:
822
823         # first, gather all unique variables referenced in this stage
824         # this "gather" step is currently required because we're being lazy
825         # in terms of index management in the scenario tree - which
826         # should really be done in terms of sets of indices.
827
828         stage_variables = {}
829         for (reference_variable, index_template, reference_index) in stage._variables:
830            if reference_variable.name not in stage_variables.keys():
831               stage_variables[reference_variable.name] = reference_variable
832
833         # next, create min/avg/max parameters for each variable in the corresponding tree node.
834         # NOTE: the parameter names below could really be empty, as they are never referenced
835         #       explicitly.
836         for (variable_name, reference_variable) in stage_variables.items():
837            for tree_node in stage._tree_nodes:
838
839               new_min_index = reference_variable._index
840               new_min_parameter_name = "NODEMIN_"+reference_variable.name
841               new_min_parameter = Param(new_min_index,name=new_min_parameter_name)
842               for index in new_min_index:
843                  new_min_parameter[index] = 0.0
844               tree_node._minimums[reference_variable.name] = new_min_parameter                                   
845
846               new_avg_index = reference_variable._index
847               new_avg_parameter_name = "NODEAVG_"+reference_variable.name
848               new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
849               for index in new_avg_index:
850                  new_avg_parameter[index] = 0.0
851               tree_node._averages[reference_variable.name] = new_avg_parameter
852
853               new_max_index = reference_variable._index
854               new_max_parameter_name = "NODEMAX_"+reference_variable.name
855               new_max_parameter = Param(new_max_index,name=new_max_parameter_name)
856               for index in new_max_index:
857                  new_max_parameter[index] = 0.0
858               tree_node._maximums[reference_variable.name] = new_max_parameter
859
860      # the objective functions are modified throughout the course of PH iterations.
861      # save the original, as a baseline to modify in subsequent iterations. reserve
862      # the original objectives, for subsequent modification.
863      self._original_objective_expression = {}
864      for instance_name, instance in self._instances.items():
865         objective_name = instance.active_components(Objective).keys()[0]
866         self._original_objective_expression[instance_name] = instance.active_components(Objective)[objective_name]._data[None].expr
867
868      # cache the number of discrete and continuous variables in the master instance. this value
869      # is of general use, e.g., in the converger classes and in plugins.
870      (self._total_discrete_vars,self._total_continuous_vars) = self.compute_variable_counts()
871      if self._verbose is True:
872         print "Total number of discrete instance variables="+str(self._total_discrete_vars)
873         print "Total number of continuous instance variables="+str(self._total_continuous_vars)
874
875      # track the total number of fixed variables of each category at the end of each PH iteration.
876      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
877
878      # indicate that we're ready to run.
879      self._initialized = True
880
881      if self._verbose is True:
882         print "PH successfully created model instances for all scenarios"
883
884      self._init_end_time = time.time()
885
886      if self._verbose is True:
887         print "PH is successfully initialized"
888         if self._output_times is True:
889            print "Initialization time=%8.2f seconds" % (self._init_end_time - self._init_start_time)
890
891      # let plugins know if they care.
892      if self._verbose is True:
893         print "Initializing PH plugins"
894      for plugin in self._ph_plugins:
895         plugin.post_ph_initialization(self)
896      if self._verbose is True:
897         print "PH plugin initialization complete"
898
899   """ Perform the non-weighted scenario solves and form the initial w and xbars.
900   """   
901   def iteration_0_solve(self):
902
903      if self._verbose is True:
904         print "------------------------------------------------"
905         print "Starting PH iteration 0 solves"
906
907      self._current_iteration = 0
908
909      solve_start_time = time.time()
910
911      # STEP 0: set up all global solver options.
912      self._solver.mipgap = self._mipgap
913
914      # STEP 1: queue up the solves for all scenario sub-problems and
915      #         grab all the action handles for the subsequent barrier sync.
916
917      action_handles = []
918      action_handle_instance_map = {}
919
920      for scenario in self._scenario_tree._scenarios:
921         
922         instance = self._instances[scenario._name]
923         
924         if self._verbose is True:
925            print "Queuing solve for scenario=" + scenario._name
926
927         # IMPT: You have to re-presolve if approximating continuous variable penalty terms with a
928         #       piecewise linear function. otherwise, the newly introduced variables won't be flagged
929         #       as unused (as is correct for iteration 0), and the output writer will crater.
930         if self._linearize_nonbinary_penalty_terms > 0:
931           instance.presolve()           
932           
933         # there's nothing to warm-start from in iteration 0, so don't include the keyword in the solve call.
934         # the reason you don't want to include it is that some solvers don't know how to handle the keyword
935         # at all (despite it being false). you might want to solve iteration 0 solves using some other solver.
936
937         new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
938         action_handle_instance_map[scenario._name] = new_action_handle
939
940         action_handles.append(new_action_handle)
941
942      # STEP 2: barrier sync for all scenario sub-problem solves.
943      if self._verbose is True:
944         print "Waiting for scenario sub-problem solves"
945      self._solver_manager.wait_all(action_handles)
946      if self._verbose is True:     
947         print "Scenario sub-problem solves completed"     
948
949      solve_end_time = time.time()
950      self._cumulative_solve_time += (solve_end_time - solve_start_time)
951
952      if self._output_times is True:
953         print "Aggregate sub-problem solve time=%8.2f" % (solve_end_time - solve_start_time)
954
955      # STEP 3: Load the results!
956      for scenario_name, action_handle in action_handle_instance_map.items():
957
958         if self._verbose is True:         
959            print "Successfully processed results for scenario="+scenario_name
960
961         instance = self._instances[scenario_name]
962         results = self._solver_manager.get_results(action_handle)
963
964         if len(results.solution) == 0:
965            results.write(num=1)
966            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
967
968         if self._output_solver_results is True:
969            print "Results for scenario=",scenario_name
970            results.write(num=1)           
971
972         instance.load(results)
973
974         if self._verbose is True:                 
975            print "Successfully loaded solution for scenario="+scenario_name
976
977      if self._verbose is True:
978         print "Successfully completed PH iteration 0 solves - solution statistics:"
979         print "         Scenario              Objective                  Value"
980         for scenario in self._scenario_tree._scenarios:
981            instance = self._instances[scenario._name]
982            for objective_name in instance.active_components(Objective):
983               objective = instance.active_components(Objective)[objective_name]
984               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
985         print "------------------------------------------------"
986
987   #
988   # recompute the averages, minimum, and maximum statistics for all variables to be blended by PH, i.e.,
989   # not appearing in the final stage. technically speaking, the min/max aren't required by PH, but they
990   # are used often enough to warrant their computation and it's basically free if you're computing the
991   # average.
992   #
993   def update_variable_statistics(self):
994
995      start_time = time.time()
996     
997      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage
998         
999         for tree_node in stage._tree_nodes:
1000               
1001            for (variable, index_template, variable_indices) in stage._variables:
1002                 
1003               variable_name = variable.name
1004                     
1005               avg_parameter_name = "PHAVG_"+variable_name
1006                 
1007               for index in variable_indices:
1008                  min = float("inf")
1009                  avg = 0.0
1010                  max = float("-inf")
1011                  node_probability = 0.0
1012                     
1013                  is_used = True # until proven otherwise                     
1014                  for scenario in tree_node._scenarios:
1015                       
1016                     instance = self._instances[scenario._name]
1017                       
1018                     if getattr(instance,variable_name)[index].status == VarStatus.unused:
1019                        is_used = False
1020                     else:                       
1021                        node_probability += scenario._probability
1022                        var_value = getattr(instance, variable.name)[index].value
1023                        if var_value < min:
1024                           min = var_value
1025                        avg += (scenario._probability * var_value)
1026                        if var_value > max:
1027                           max = var_value
1028
1029                  if is_used is True:
1030                     tree_node._minimums[variable.name][index] = min
1031                     tree_node._averages[variable.name][index] = avg / node_probability
1032                     tree_node._maximums[variable.name][index] = max                       
1033
1034                     # distribute the newly computed average to the xbar variable in
1035                     # each instance/scenario associated with this node. only do this
1036                     # if the variable is used!
1037                     for scenario in tree_node._scenarios:
1038                        instance = self._instances[scenario._name]
1039                        avg_parameter = getattr(instance, avg_parameter_name)
1040                        avg_parameter[index] = avg / node_probability
1041
1042      end_time = time.time()
1043      self._cumulative_xbar_time += (end_time - start_time)
1044
1045   def update_weights(self):
1046
1047      # because the weight updates rely on the xbars, and the xbars are node-based,
1048      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
1049      start_time = time.time()     
1050
1051      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.
1052         
1053         for tree_node in stage._tree_nodes:
1054
1055            for (variable, index_template, variable_indices) in stage._variables:
1056
1057               variable_name = variable.name
1058               blend_parameter_name = "PHBLEND_"+variable_name
1059               weight_parameter_name = "PHWEIGHT_"+variable_name
1060               rho_parameter_name = "PHRHO_"+variable_name
1061
1062               for index in variable_indices:
1063
1064                  tree_node_average = tree_node._averages[variable.name][index]()
1065
1066                  for scenario in tree_node._scenarios:
1067
1068                     instance = self._instances[scenario._name]
1069
1070                     if getattr(instance,variable.name)[index].status != VarStatus.unused:
1071
1072                        # we are currently not updating weights if blending is disabled for a variable.
1073                        # this is done on the premise that unless you are actively trying to move
1074                        # the variable toward the mean, the weights will blow up and be huge by the
1075                        # time that blending is activated.
1076                        variable_blend_indicator = getattr(instance, blend_parameter_name)[index]()                             
1077
1078                        # get the weight and rho parameters for this variable/index combination.
1079                        rho_value = getattr(instance, rho_parameter_name)[index]()
1080                        current_variable_weight = getattr(instance, weight_parameter_name)[index]()
1081                                               
1082                        # if I'm maximizing, invert value prior to adding (hack to implement negatives).
1083                        # probably fixed in Pyomo at this point - I just haven't checked in a long while.
1084                        if self._is_minimizing is False:
1085                           current_variable_weight = (-current_variable_weight)
1086                        current_variable_value = getattr(instance,variable.name)[index]()
1087                        new_variable_weight = current_variable_weight + variable_blend_indicator * rho_value * (current_variable_value - tree_node_average)
1088                        # I have the correct updated value, so now invert if maximizing.
1089                        if self._is_minimizing is False:
1090                           new_variable_weight = (-new_variable_weight)
1091                        getattr(instance, weight_parameter_name)[index].value = new_variable_weight
1092
1093      # we shouldn't have to re-simplify the expression, as we aren't adding any constant-variable terms - just modifying parameters.
1094
1095      end_time = time.time()
1096      self._cumulative_weight_time += (end_time - start_time)
1097
1098   def form_iteration_k_objectives(self):
1099
1100      # for each blended variable (i.e., those not appearing in the final stage),
1101      # add the linear and quadratic penalty terms to the objective.
1102      for instance_name, instance in self._instances.items():
1103         
1104         objective_name = instance.active_components(Objective).keys()[0]         
1105         objective = instance.active_components(Objective)[objective_name]
1106         # clone the objective, because we're going to augment (repeatedly) the original objective.
1107         objective_expression = self._original_objective_expression[instance_name].clone() 
1108         # the quadratic expression is really treated as just a list - eventually should be treated as a full expression.
1109         quad_expression = 0.0
1110         
1111         for stage in self._scenario_tree._stages[:-1]: # skip the last stage, as no blending occurs
1112
1113            variable_tree_node = None
1114            for node in stage._tree_nodes:
1115               for scenario in node._scenarios:
1116                  if scenario._name == instance_name:
1117                     variable_tree_node = node
1118                     break
1119
1120            for (reference_variable, index_template, variable_indices) in stage._variables:
1121
1122               variable_name = reference_variable.name
1123               variable_type = reference_variable.domain
1124
1125               w_parameter_name = "PHWEIGHT_"+variable_name
1126               w_parameter = instance.active_components(Param)[w_parameter_name]
1127                 
1128               average_parameter_name = "PHAVG_"+variable_name
1129               average_parameter = instance.active_components(Param)[average_parameter_name]
1130
1131               rho_parameter_name = "PHRHO_"+variable_name
1132               rho_parameter = instance.active_components(Param)[rho_parameter_name]
1133
1134               blend_parameter_name = "PHBLEND_"+variable_name
1135               blend_parameter = instance.active_components(Param)[blend_parameter_name]
1136
1137               node_min_parameter = variable_tree_node._minimums[variable_name]
1138               node_max_parameter = variable_tree_node._maximums[variable_name]               
1139
1140               quad_penalty_term_variable = None
1141               if self._linearize_nonbinary_penalty_terms > 0:
1142                  quad_penalty_term_variable_name = "PHQUADPENALTY_"+variable_name
1143                  quad_penalty_term_variable = instance.active_components(Var)[quad_penalty_term_variable_name]
1144
1145               instance_variable = instance.active_components(Var)[variable_name]
1146
1147               for index in variable_indices:
1148
1149                  if (instance_variable[index].status is not VarStatus.unused) and (instance_variable[index].fixed is False):
1150
1151                     # add the linear (w-weighted) term is a consistent fashion, independent of variable type.
1152                     # if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes for reasons currently unknown.
1153                     objective_expression += (w_parameter[index] * instance_variable[index])
1154
1155                     # there are some cases in which a user may want to avoid the proximal term completely.
1156                     # it's of course only a good idea when there are at least bounds (both lb and ub) on
1157                     # the variables to be blended.
1158                     if self._drop_proximal_terms is False:
1159
1160                        # deal with binaries
1161                        if isinstance(variable_type, BooleanSet) is True:
1162
1163                           if self._retain_quadratic_binary_terms is False:
1164                              # this rather ugly form of the linearized quadratic expression term is required
1165                              # due to a pyomo bug - the form (rho/2) * (x+y+z) chokes in presolve when distributing
1166                              # over the sum.
1167                              new_term = (blend_parameter[index] * rho_parameter[index] / 2.0 * instance_variable[index]) - \
1168                                         (blend_parameter[index] * rho_parameter[index] * average_parameter[index] * instance_variable[index]) + \
1169                                         (blend_parameter[index] * rho_parameter[index] / 2.0 * average_parameter[index] * average_parameter[index])                             
1170                              if objective.sense is minimize:
1171                                 objective_expression += new_term
1172                              else:
1173                                 objective_expression -= new_term                                 
1174                           else:
1175                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)
1176
1177                        # deal with everything else
1178                        else:
1179
1180                           if self._linearize_nonbinary_penalty_terms > 0:
1181
1182                              # the variables are easy - just a single entry.
1183                              if objective.sense is minimize:
1184                                 objective_expression += (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
1185                              else:
1186                                 objective_expression -= (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
1187                                 
1188                              # the constraints are somewhat nastier.
1189
1190                              # TBD - DEFINE CONSTRAINTS ON-THE-FLY?? (INDIVIDUALLY NAMED FOR NOW - CREATE INDEX SETS!) - OR A LEAST AN INDEX SET PER "PIECE"
1191                              xavg = average_parameter[index]
1192                              x = instance_variable[index]
1193
1194                              lb = None
1195                              ub = None
1196
1197                              if x.lb is None:
1198                                 raise ValueError, "No lower bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
1199                              else:
1200                                 lb = x.lb()
1201
1202                              if x.ub is None:
1203                                 raise ValueError, "No upper bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
1204                              else:
1205                                 ub = x.ub()
1206
1207                              if x.lb == x.ub:
1208                                 print "***WARNING - LB EQUALS UB"
1209                                 print "VALUE=",x.lb
1210                                 print "VARIABLE="+variable_name+indexToString(index)
1211
1212                              node_min = node_min_parameter[index]()
1213                              node_max = node_max_parameter[index]()
1214
1215                              # compute the breakpoint sequence according to the specified strategy.
1216                              breakpoints = []
1217                              if self._breakpoint_strategy == 0:
1218                                 breakpoints = self.compute_uniform_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
1219                              elif self._breakpoint_strategy == 1:
1220                                 breakpoints = self.compute_uniform_between_nodestat_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
1221                              elif self._breakpoint_strategy == 2:
1222                                 breakpoints = self.compute_uniform_between_woodruff_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
1223                              elif self._breakpoint_strategy == 3:
1224                                 breakpoints = self.compute_exponential_from_mean_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)                                                                                                   
1225                              else:
1226                                 raise ValueError, "A breakpoint distribution strategy="+str(self._breakpoint_strategy)+" is currently not supported within PH!"
1227
1228                              for i in range(0,len(breakpoints)-1):
1229
1230                                 this_lb = breakpoints[i]
1231                                 this_ub = breakpoints[i+1]
1232
1233                                 piece_constraint_name = "QUAD_PENALTY_PIECE_"+str(i)+"_"+variable_name+str(index)
1234                                 if hasattr(instance, piece_constraint_name) is False:
1235                                    # this is the first time the constraint is being added - add it to the list of PH-specific constraints for this instance.
1236                                    self._instance_augmented_attributes[instance_name].append(piece_constraint_name)                                                                           
1237                                 piece_constraint = Constraint(name=piece_constraint_name)
1238                                 piece_constraint.model = instance
1239                                 piece_expression = self._create_piecewise_constraint_expression(this_lb, this_ub, x, xavg, quad_penalty_term_variable[index])
1240                                 piece_constraint.add(None, piece_expression)
1241                                 setattr(instance, piece_constraint_name, piece_constraint)                                   
1242
1243                           else:
1244
1245                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)                           
1246                   
1247         # strictly speaking, this probably isn't necessary - parameter coefficients won't get
1248         # pre-processed out of the expression tree. however, if the under-the-hood should change,
1249         # we'll be covered.
1250         objective_expression.simplify(instance)
1251         instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression
1252         # if we are linearizing everything, then nothing will appear in the quadratic expression -
1253         # don't add the empty "0.0" expression to the objective. otherwise, the output file won't
1254         # be properly generated.
1255         if quad_expression != 0.0:
1256           instance.active_components(Objective)[objective_name]._quad_subexpr = quad_expression
1257
1258   def iteration_k_solve(self):
1259
1260     if self._verbose is True:
1261        print "------------------------------------------------"       
1262        print "Starting PH iteration " + str(self._current_iteration) + " solves"
1263
1264     # cache the objective values generated by PH for output at the end of this function.
1265     ph_objective_values = {}
1266
1267     solve_start_time = time.time()
1268
1269     # STEP 0: set up all global solver options.
1270     self._solver.mipgap = self._mipgap     
1271
1272     # STEP 1: queue up the solves for all scenario sub-problems and
1273     #         grab all of the action handles for the subsequent barrier sync.
1274
1275     action_handles = []
1276     action_handle_instance_map = {}
1277
1278     for scenario in self._scenario_tree._scenarios:     
1279
1280        instance = self._instances[scenario._name]
1281
1282        if self._verbose is True:
1283           print "Queuing solve for scenario=" + scenario._name
1284
1285        # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
1286        # don't do this, you won't see any chance in the output files as you vary the problem parameters!
1287        # ditto for instance fixing!
1288        instance.presolve()
1289
1290        # once past iteration 0, there is always a feasible solution from which to warm-start.
1291        # however, you might want to disable warm-start when the solver is behaving badly (which does happen).
1292        new_action_handle = None
1293        if (self._disable_warmstarts is False) and (self._solver.warm_start_capable() is True):
1294           new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
1295        else:
1296           new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)           
1297
1298        action_handle_instance_map[scenario._name] = new_action_handle
1299
1300        action_handles.append(new_action_handle)
1301
1302     # STEP 2: barrier sync for all scenario sub-problem solves.
1303     if self._verbose is True:           
1304        print "Waiting for scenario sub-problem solves"
1305     self._solver_manager.wait_all(action_handles)
1306     if self._verbose is True:               
1307        print "Scenario sub-problem solves completed"
1308       
1309     solve_end_time = time.time()
1310     self._cumulative_solve_time += (solve_end_time - solve_start_time)
1311
1312     if self._output_times is True:
1313        print "Aggregate sub-problem solve time=%8.2f" % (solve_end_time - solve_start_time)
1314
1315     # STEP 3: Load the results!
1316     for scenario_name, action_handle in action_handle_instance_map.items():
1317
1318        if self._verbose is True:         
1319           print "Successfully processed results for scenario="+scenario_name
1320
1321        instance = self._instances[scenario_name]
1322        results = self._solver_manager.get_results(action_handle)
1323
1324        if len(results.solution) == 0:
1325           raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
1326
1327        if self._output_solver_results is True:
1328           print "Results for scenario=", scenario_name
1329           results.write(num=1)           
1330
1331        instance.load(results)
1332
1333        if self._verbose is True:                 
1334           print "Successfully loaded solution for scenario="+scenario_name
1335
1336        # we're assuming there is a single solution.
1337        # the "value" attribute is a pre-defined feature of any solution - it is relative to whatever
1338        # objective was selected during optimization, which of course should be the PH objective.
1339        ph_objective_values[instance.name] = float(results.solution(0).objective['f'].value)
1340
1341     if self._verbose is True:
1342        print "Successfully completed PH iteration " + str(self._current_iteration) + " solves - solution statistics:"
1343        print "  Scenario             PH Objective             Cost Objective"
1344        for scenario in self._scenario_tree._scenarios:
1345           instance = self._instances[scenario._name]
1346           for objective_name in instance.active_components(Objective):
1347              objective = instance.active_components(Objective)[objective_name]
1348              print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], 0.0)
1349
1350   def solve(self):
1351
1352      self._solve_start_time = time.time()
1353      self._cumulative_solve_time = 0.0
1354      self._cumulative_xbar_time = 0.0
1355      self._cumulative_weight_time = 0.0
1356
1357      print "Starting PH"
1358
1359      if self._initialized == False:
1360         raise RuntimeError, "PH is not initialized - cannot invoke solve() method"
1361
1362      print "Initiating PH iteration=" + `self._current_iteration`
1363
1364      self.iteration_0_solve()
1365
1366      # update variable statistics prior to any output.
1367      self.update_variable_statistics()     
1368
1369      if self._verbose is True:
1370         print "Variable values following scenario solves:"
1371         self.pprint(False,False,True,False)
1372
1373      # let plugins know if they care.
1374      for plugin in self._ph_plugins:     
1375         plugin.post_iteration_0_solves(self)
1376
1377      # update the fixed variable statistics.
1378      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
1379
1380      if self._verbose is True:
1381         print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1382         print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1383
1384      self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1385      print "Convergence metric=%12.4f" % self._converger.lastMetric()
1386
1387      self.update_weights()
1388
1389      # let plugins know if they care.
1390      for plugin in self._ph_plugins:     
1391         plugin.post_iteration_0(self)
1392
1393      # checkpoint if it's time - which it always is after iteration 0,
1394      # if the interval is >= 1!
1395      if (self._checkpoint_interval > 0):
1396         self.checkpoint(0)
1397
1398      # there is an upper bound on the number of iterations to execute -
1399      # the actual bound depends on the converger supplied by the user.
1400      for i in range(1, self._max_iterations+1):
1401
1402         self._current_iteration = self._current_iteration + 1                 
1403
1404         print "Initiating PH iteration=" + `self._current_iteration`         
1405
1406         if self._verbose is True:
1407            print "Variable averages and weights prior to scenario solves:"
1408            self.pprint(True,True,False,False)
1409
1410         # with the introduction of piecewise linearization, the form of the
1411         # penalty-weighted objective is no longer fixed. thus, we need to
1412         # create the objectives each PH iteration.
1413         self.form_iteration_k_objectives()           
1414
1415         self.iteration_k_solve()
1416
1417         # update variable statistics prior to any output.
1418         self.update_variable_statistics()
1419         
1420         if self._verbose is True:
1421            print "Variable values following scenario solves:"
1422            self.pprint(False,False,True,False)
1423
1424         # we don't technically have to do this at the last iteration,
1425         # but with checkpointing and re-starts, you're never sure
1426         # when you're executing the last iteration.
1427         self.update_weights()
1428
1429         # let plugins know if they care.
1430         for plugin in self._ph_plugins:
1431            plugin.post_iteration_k_solves(self)
1432
1433         # update the fixed variable statistics.
1434         (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
1435
1436         if self._verbose is True:
1437            print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1438            print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1439
1440         # let plugins know if they care.
1441         for plugin in self._ph_plugins:
1442            plugin.post_iteration_k(self)
1443
1444         # at this point, all the real work of an iteration is complete.
1445
1446         # checkpoint if it's time.
1447         if (self._checkpoint_interval > 0) and (i % self._checkpoint_interval is 0):
1448            self.checkpoint(i)
1449
1450         # check for early termination.
1451         self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1452         print "Convergence metric=%12.4f" % self._converger.lastMetric()
1453
1454         if self._converger.isConverged(self) is True:
1455            if self._total_discrete_vars == 0:
1456               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)
1457            else:
1458               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)+" or all discrete variables are fixed"               
1459            break
1460
1461         # if we're terminating due to exceeding the maximum iteration count, print a message
1462         # indicating so - otherwise, you get a quiet, information-free output trace.
1463         if i == self._max_iterations:
1464            print "Halting PH - reached maximal iteration count="+str(self._max_iterations)
1465
1466      if self._verbose is True:
1467         print "Number of discrete variables fixed before final plugin calls="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1468         print "Number of continuous variables fixed before final plugin calls="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"           
1469
1470      # let plugins know if they care. do this before
1471      # the final solution / statistics output, as the plugins
1472      # might do some final tweaking.
1473      for plugin in self._ph_plugins:
1474         plugin.post_ph_execution(self)
1475
1476      # update the fixed variable statistics - the plugins might have done something.
1477      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()                       
1478
1479      print "PH complete"
1480
1481      print "Convergence history:"
1482      self._converger.pprint()
1483
1484      print "Final number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1485      print "Final number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"         
1486
1487      print "Final variable values:"
1488      self.pprint(False,False,True,True)         
1489
1490      print "Final costs:"
1491      self._scenario_tree.pprintCosts(self._instances)
1492
1493      self._solve_end_time = time.time()
1494
1495      if (self._verbose is True) and (self._output_times is True):
1496         print "Overall run-time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1497
1498      # cleanup the scenario instances for post-processing - ideally, we want to leave them in
1499      # their original state, minus all the PH-specific stuff. we don't do all cleanup (leaving
1500      # things like rhos, etc), but we do clean up constraints, as that really hoses up the ef writer.
1501      self._cleanup_scenario_instances()
1502
1503   #
1504   # prints a summary of all collected time statistics
1505   #
1506   def print_time_stats(self):
1507
1508      print "PH run-time statistics (user):"
1509
1510      print "Initialization time=  %8.2f seconds" % (self._init_end_time - self._init_start_time)
1511      print "Overall solve time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1512      print "Scenario solve time=  %8.2f seconds" % self._cumulative_solve_time
1513      print "Average update time=  %8.2f seconds" % self._cumulative_xbar_time
1514      print "Weight update time=   %8.2f seconds" % self._cumulative_weight_time
1515
1516   #
1517   # a utility to determine whether to output weight / average / etc. information for
1518   # a variable/node combination. when the printing is moved into a callback/plugin,
1519   # this routine will go there. for now, we don't dive down into the node resolution -
1520   # just the variable/stage.
1521   #
1522   def should_print(self, stage, variable, variable_indices):
1523
1524      if self._output_continuous_variable_stats is False:
1525
1526         variable_type = variable.domain         
1527
1528         if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):
1529
1530            return False
1531
1532      return True
1533     
1534   #
1535   # pretty-prints the state of the current variable averages, weights, and values.
1536   # inputs are booleans indicating which components should be output.
1537   #
1538   def pprint(self, output_averages, output_weights, output_values, output_fixed):
1539
1540      if self._initialized is False:
1541         raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"         
1542     
1543      # print tree nodes and associated variable/xbar/ph information in stage-order
1544      # we don't blend in the last stage, so we don't current care about printing the associated information.
1545      for stage in self._scenario_tree._stages[:-1]:
1546
1547         print "\tStage=" + stage._name
1548
1549         num_outputs_this_stage = 0 # tracks the number of outputs on a per-index basis.
1550
1551         for (variable, index_template, variable_indices) in stage._variables:
1552
1553            variable_name = variable.name
1554
1555            if self.should_print(stage, variable, variable_indices) is True:
1556
1557               num_outputs_this_variable = 0 # track, so we don't output the variable names unless there is an entry to report.
1558
1559               for index in variable_indices:               
1560
1561                  weight_parameter_name = "PHWEIGHT_"+variable_name
1562
1563                  num_outputs_this_index = 0 # track, so we don't output the variable index more than once.
1564
1565                  for tree_node in stage._tree_nodes:                 
1566
1567                     # determine if the variable/index pair is used across the set of scenarios (technically,
1568                     # it should be good enough to check one scenario). ditto for "fixed" status. fixed does
1569                     # imply unused (see note below), but we care about the fixed status when outputting
1570                     # final solutions.
1571
1572                     is_used = True # should be consistent across scenarios, so one "unused" flags as invalid.
1573                     is_fixed = False
1574
1575                     for scenario in tree_node._scenarios:
1576                        instance = self._instances[scenario._name]
1577                        variable_value = getattr(instance,variable_name)[index]
1578                        if variable_value.status == VarStatus.unused:
1579                           is_used = False
1580                        if variable_value.fixed is True:
1581                           is_fixed = True
1582 
1583                     # IMPT: this is far from obvious, but variables that are fixed will - because
1584                     #       presolve will identify them as constants and eliminate them from all
1585                     #       expressions - be flagged as "unused" and therefore not output.
1586
1587                     if ((output_fixed is True) and (is_fixed is True)) or (is_used is True):
1588
1589                           minimum_value = tree_node._minimums[variable_name][index]()
1590                           maximum_value = tree_node._maximums[variable_name][index]()
1591
1592                           num_outputs_this_stage = num_outputs_this_stage + 1                           
1593                           num_outputs_this_variable = num_outputs_this_variable + 1
1594                           num_outputs_this_index = num_outputs_this_index + 1
1595
1596                           if num_outputs_this_variable == 1:
1597                              print "\t\tVariable=",variable_name
1598
1599                           if num_outputs_this_index == 1:
1600                              print "\t\t\tIndex:", indexToString(index)                             
1601
1602                           print "\t\t\t\tTree Node=",tree_node._name,"\t\t (Scenarios: ",                             
1603                           for scenario in tree_node._scenarios:
1604                              print scenario._name," ",
1605                              if scenario == tree_node._scenarios[-1]:
1606                                 print ")"
1607                           
1608                           if output_values is True:
1609                              average_value = tree_node._averages[variable_name][index]()
1610                              print "\t\t\t\tValues: ",                       
1611                              for scenario in tree_node._scenarios:
1612                                 instance = self._instances[scenario._name]
1613                                 this_value = getattr(instance,variable_name)[index].value
1614                                 print "%12.4f" % this_value,
1615                                 if scenario == tree_node._scenarios[-1]:
1616                                    print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1617                                    print "    Avg=%12.4f" % (average_value),
1618                                    print ""
1619                           if output_weights:
1620                              print "\t\t\t\tWeights: ",
1621                              for scenario in tree_node._scenarios:
1622                                 instance = self._instances[scenario._name]
1623                                 print "%12.4f" % getattr(instance,weight_parameter_name)[index].value,
1624                                 if scenario == tree_node._scenarios[-1]:
1625                                    print ""
1626
1627                           if output_averages:
1628                              print "\t\t\t\tAverage: %12.4f" % (tree_node._averages[variable_name][index].value)
1629
1630         if num_outputs_this_stage == 0:
1631            print "\t\tNo non-converged variables in stage"
1632
1633         # cost variables aren't blended, so go through the gory computation of min/max/avg.
1634         # we currently always print these.
1635         cost_variable_name = stage._cost_variable[0].name
1636         cost_variable_index = stage._cost_variable[1]
1637         if cost_variable_index is None:
1638            print "\t\tCost Variable=" + cost_variable_name
1639         else:
1640            print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)           
1641         for tree_node in stage._tree_nodes:
1642            print "\t\t\tTree Node=" + tree_node._name + "\t\t (Scenarios: ",
1643            for scenario in tree_node._scenarios:
1644               print scenario._name," ",
1645               if scenario == tree_node._scenarios[-1]:
1646                  print ")"
1647            maximum_value = 0.0
1648            minimum_value = 0.0
1649            sum_values = 0.0
1650            num_values = 0
1651            first_time = True
1652            print "\t\t\tValues: ",                       
1653            for scenario in tree_node._scenarios:
1654                instance = self._instances[scenario._name]
1655                this_value = getattr(instance,cost_variable_name)[cost_variable_index].value
1656                print "%12.4f" % this_value,
1657                num_values += 1
1658                sum_values += this_value
1659                if first_time is True:
1660                   first_time = False
1661                   maximum_value = this_value
1662                   minimum_value = this_value
1663                else:
1664                   if this_value > maximum_value:
1665                      maximum_value = this_value
1666                   if this_value < minimum_value:
1667                      minimum_value = this_value
1668                if scenario == tree_node._scenarios[-1]:
1669                   print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1670                   print "    Avg=%12.4f" % (sum_values/num_values),
1671                   print ""
1672           
Note: See TracBrowser for help on using the repository browser.