source: coopr.pysp/trunk/coopr/pysp/ph.py @ 2201

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

Update to Coopr to account for changes in PyUtilib? package names.

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