source: coopr.pysp/stable/2.3/coopr/pysp/ph.py @ 2317

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

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

........

r2246 | wehart | 2010-02-01 21:10:48 -0700 (Mon, 01 Feb 2010) | 2 lines


Tagging coopr.pysp release 2.2

........

r2247 | wehart | 2010-02-01 21:46:08 -0700 (Mon, 01 Feb 2010) | 2 lines


Documentation fix.

........

r2248 | jwatson | 2010-02-03 08:44:59 -0700 (Wed, 03 Feb 2010) | 3 lines


Adding os.path.expanduser wrappers around all directory/filenames, to facilitate correct processing of ~ characters.

........

r2249 | jwatson | 2010-02-03 09:28:25 -0700 (Wed, 03 Feb 2010) | 3 lines


Misc fixes.

........

r2254 | jwatson | 2010-02-03 21:09:38 -0700 (Wed, 03 Feb 2010) | 3 lines


Changed PyomoModelData? call from add_data_mumble() to add().

........

r2255 | jwatson | 2010-02-03 21:38:35 -0700 (Wed, 03 Feb 2010) | 3 lines


Re-factoring of PH options parser code to accomodate MRP work being done by DLW.

........

r2256 | jwatson | 2010-02-03 22:18:11 -0700 (Wed, 03 Feb 2010) | 3 lines


Major speed improvements in the EF writer by avoiding Python deep-copes - saves a few orders of magnitude of run-time.

........

r2260 | jwatson | 2010-02-04 21:44:29 -0700 (Thu, 04 Feb 2010) | 3 lines


Refactoring of ph initialization routines to support sampling and bundling.

........

r2261 | jwatson | 2010-02-04 21:46:37 -0700 (Thu, 04 Feb 2010) | 3 lines


Renaming ph_script module to phinit, which is more accurate with the newly factored, library-like functionality.

........

r2262 | jwatson | 2010-02-04 21:54:50 -0700 (Thu, 04 Feb 2010) | 3 lines


Minor improvement to phinit functionality.

........

r2267 | jwatson | 2010-02-05 15:09:02 -0700 (Fri, 05 Feb 2010) | 3 lines


Initial PySP unit tests!!!

........

r2268 | wehart | 2010-02-05 15:25:52 -0700 (Fri, 05 Feb 2010) | 2 lines


A fix to the tests.

........

r2269 | jwatson | 2010-02-05 15:46:35 -0700 (Fri, 05 Feb 2010) | 3 lines


Added SIZES3 PySP test and added absolute paths to "runph" script.

........

r2270 | jwatson | 2010-02-05 16:01:08 -0700 (Fri, 05 Feb 2010) | 3 lines


PySP tests need absolute output paths!

........

r2271 | jwatson | 2010-02-05 16:20:01 -0700 (Fri, 05 Feb 2010) | 3 lines


Making unit tests for PySP compatible with coverage utilities. Farmer examples work, SIZES3 doesn't for some reason.

........

r2272 | jwatson | 2010-02-05 16:59:46 -0700 (Fri, 05 Feb 2010) | 3 lines


Various fixes to PySP unit tests. Changing name of testphextension to examplephextension - with the test prefix, coverage tests import the module, which causes all kinds of issues.

........

r2275 | jwatson | 2010-02-06 13:45:11 -0700 (Sat, 06 Feb 2010) | 3 lines


Testing improvements. From lpython, the tests run individually just fine. In aggregate, only the first run passes - independent of what test that might be! Something to stare at later....

........

r2277 | jwatson | 2010-02-06 23:30:13 -0700 (Sat, 06 Feb 2010) | 3 lines


Fixed problem with runph --profile option, broken by my recent factoring of phinit.py.

........

r2288 | jwatson | 2010-02-08 13:32:41 -0700 (Mon, 08 Feb 2010) | 3 lines


Significant initialization speed reductions in the WW PH extension for PySP.

........

r2290 | jwatson | 2010-02-08 19:22:32 -0700 (Mon, 08 Feb 2010) | 1 line


Initial commit of multi-stage capacity expansion problem in PySP

........

r2291 | jwatson | 2010-02-08 19:23:24 -0700 (Mon, 08 Feb 2010) | 1 line


Miscellaneous fix to ef writer involving indexed cost variables

........

r2296 | jwatson | 2010-02-09 18:56:48 -0700 (Tue, 09 Feb 2010) | 3 lines


Removing monster-sized LP file from the PySP forestry examples directory.

........

r2297 | jwatson | 2010-02-09 18:59:05 -0700 (Tue, 09 Feb 2010) | 3 lines


The extensive forms in the PySP forestry example were massive - and are now gone.

........

r2298 | jwatson | 2010-02-09 19:03:09 -0700 (Tue, 09 Feb 2010) | 1 line


Removing output logs for PySP network flow example

........

r2299 | jwatson | 2010-02-09 19:04:12 -0700 (Tue, 09 Feb 2010) | 1 line


Removing a big network flow EF

........

r2300 | jwatson | 2010-02-09 19:05:08 -0700 (Tue, 09 Feb 2010) | 3 lines


Removing PySP cap example EF to free up space.

........

r2301 | jwatson | 2010-02-09 19:14:26 -0700 (Tue, 09 Feb 2010) | 1 line


Performance improvements to PH obtained by processing scenario sub-problem results as they come in, instead of waiting for them after a solver barrier sync

........

File size: 87.7 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(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      scenario_action_handle_map = {} # maps scenario names to action handles
1012      action_handle_scenario_map = {} # maps action handles to scenario names
1013
1014      for scenario in self._scenario_tree._scenarios:
1015         
1016         instance = self._instances[scenario._name]
1017         
1018         if self._verbose is True:
1019            print "Queuing solve for scenario=" + scenario._name
1020
1021         # IMPT: You have to re-presolve if approximating continuous variable penalty terms with a
1022         #       piecewise linear function. otherwise, the newly introduced variables won't be flagged
1023         #       as unused (as is correct for iteration 0), and the output writer will crater.
1024         # IMPT: I decided to presolve unconditionally, as PH extensions can add arbitrary components
1025         #       to the base scenario instances - and the variable values/etc. need to be collectged.
1026         instance.presolve()
1027           
1028         # there's nothing to warm-start from in iteration 0, so don't include the keyword in the solve call.
1029         # the reason you don't want to include it is that some solvers don't know how to handle the keyword
1030         # at all (despite it being false). you might want to solve iteration 0 solves using some other solver.
1031
1032         new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
1033         scenario_action_handle_map[scenario._name] = new_action_handle
1034         action_handle_scenario_map[new_action_handle] = scenario._name
1035
1036         action_handles.append(new_action_handle)
1037
1038      # STEP 2: loop for the solver results, reading them and loading
1039      #         them into instances as they are available.
1040
1041      if self._verbose is True:
1042         print "Waiting for scenario sub-problem solves"
1043
1044      num_results_so_far = 0
1045
1046      while (num_results_so_far < len(self._scenario_tree._scenarios)):
1047
1048         action_handle = self._solver_manager.wait_any()
1049         results = self._solver_manager.get_results(action_handle)         
1050         scenario_name = action_handle_scenario_map[action_handle]
1051         instance = self._instances[scenario_name]         
1052
1053         if self._verbose is True:
1054            print "Results obtained for scenario="+scenario_name
1055
1056         if len(results.solution) == 0:
1057            results.write(num=1)
1058            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
1059
1060         if self._output_solver_results is True:
1061            print "Results for scenario=",scenario_name
1062            results.write(num=1)
1063
1064         start_time = time.time()
1065         instance.load(results)
1066         end_time = time.time()
1067         if self._output_times is True:
1068            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
1069
1070         if self._verbose is True:                 
1071            print "Successfully loaded solution for scenario="+scenario_name
1072
1073         num_results_so_far = num_results_so_far + 1
1074         
1075      if self._verbose is True:     
1076         print "Scenario sub-problem solves completed"     
1077
1078      solve_end_time = time.time()
1079      self._cumulative_solve_time += (solve_end_time - solve_start_time)
1080
1081      if self._output_times is True:
1082         print "Aggregate sub-problem solve time this iteration=%8.2f" % (solve_end_time - solve_start_time)
1083
1084      if self._verbose is True:
1085         print "Successfully completed PH iteration 0 solves - solution statistics:"
1086         print "         Scenario              Objective                  Value"
1087         for scenario in self._scenario_tree._scenarios:
1088            instance = self._instances[scenario._name]
1089            for objective_name in instance.active_components(Objective):
1090               objective = instance.active_components(Objective)[objective_name]
1091               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
1092         print "------------------------------------------------"
1093
1094   #
1095   # recompute the averages, minimum, and maximum statistics for all variables to be blended by PH, i.e.,
1096   # not appearing in the final stage. technically speaking, the min/max aren't required by PH, but they
1097   # are used often enough to warrant their computation and it's basically free if you're computing the
1098   # average.
1099   #
1100   def update_variable_statistics(self):
1101
1102      start_time = time.time()
1103     
1104      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage
1105         
1106         for tree_node in stage._tree_nodes:
1107               
1108            for (variable, index_template, variable_indices) in stage._variables:
1109                 
1110               variable_name = variable.name
1111                     
1112               avg_parameter_name = "PHAVG_"+variable_name
1113                 
1114               for index in variable_indices:
1115                  min = float("inf")
1116                  avg = 0.0
1117                  max = float("-inf")
1118                  node_probability = 0.0
1119                     
1120                  is_used = True # until proven otherwise                     
1121                  for scenario in tree_node._scenarios:
1122                       
1123                     instance = self._instances[scenario._name]
1124                       
1125                     if getattr(instance,variable_name)[index].status == VarStatus.unused:
1126                        is_used = False
1127                     else:                       
1128                        node_probability += scenario._probability
1129                        var_value = getattr(instance, variable.name)[index].value
1130                        if var_value < min:
1131                           min = var_value
1132                        avg += (scenario._probability * var_value)
1133                        if var_value > max:
1134                           max = var_value
1135
1136                  if is_used is True:
1137
1138                     tree_node._minimums[variable.name][index] = min
1139                     tree_node._averages[variable.name][index] = avg / node_probability
1140                     tree_node._maximums[variable.name][index] = max                       
1141
1142                     # distribute the newly computed average to the xbar variable in
1143                     # each instance/scenario associated with this node. only do this
1144                     # if the variable is used!
1145                     for scenario in tree_node._scenarios:
1146                        instance = self._instances[scenario._name]
1147                        avg_parameter = getattr(instance, avg_parameter_name)
1148                        avg_parameter[index] = avg / node_probability
1149
1150      end_time = time.time()
1151      self._cumulative_xbar_time += (end_time - start_time)
1152
1153   def update_weights(self):
1154
1155      # because the weight updates rely on the xbars, and the xbars are node-based,
1156      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
1157      start_time = time.time()     
1158
1159      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.
1160         
1161         for tree_node in stage._tree_nodes:
1162
1163            for (variable, index_template, variable_indices) in stage._variables:
1164
1165               variable_name = variable.name
1166               blend_parameter_name = "PHBLEND_"+variable_name
1167               weight_parameter_name = "PHWEIGHT_"+variable_name
1168               rho_parameter_name = "PHRHO_"+variable_name
1169
1170               for index in variable_indices:
1171
1172                  tree_node_average = tree_node._averages[variable.name][index]()
1173
1174                  for scenario in tree_node._scenarios:
1175
1176                     instance = self._instances[scenario._name]
1177
1178                     if getattr(instance,variable.name)[index].status != VarStatus.unused:
1179
1180                        # we are currently not updating weights if blending is disabled for a variable.
1181                        # this is done on the premise that unless you are actively trying to move
1182                        # the variable toward the mean, the weights will blow up and be huge by the
1183                        # time that blending is activated.
1184                        variable_blend_indicator = getattr(instance, blend_parameter_name)[index]()                             
1185
1186                        # get the weight and rho parameters for this variable/index combination.
1187                        rho_value = getattr(instance, rho_parameter_name)[index]()
1188                        current_variable_weight = getattr(instance, weight_parameter_name)[index]()
1189                                               
1190                        # if I'm maximizing, invert value prior to adding (hack to implement negatives).
1191                        # probably fixed in Pyomo at this point - I just haven't checked in a long while.
1192                        if self._is_minimizing is False:
1193                           current_variable_weight = (-current_variable_weight)
1194                        current_variable_value = getattr(instance,variable.name)[index]()
1195                        new_variable_weight = current_variable_weight + variable_blend_indicator * rho_value * (current_variable_value - tree_node_average)
1196                        # I have the correct updated value, so now invert if maximizing.
1197                        if self._is_minimizing is False:
1198                           new_variable_weight = (-new_variable_weight)
1199                        getattr(instance, weight_parameter_name)[index].value = new_variable_weight
1200
1201      # we shouldn't have to re-simplify the expression, as we aren't adding any constant-variable terms - just modifying parameters.
1202
1203      end_time = time.time()
1204      self._cumulative_weight_time += (end_time - start_time)
1205
1206   def form_iteration_k_objectives(self):
1207
1208      # for each blended variable (i.e., those not appearing in the final stage),
1209      # add the linear and quadratic penalty terms to the objective.
1210      for instance_name, instance in self._instances.items():
1211         
1212         objective_name = instance.active_components(Objective).keys()[0]         
1213         objective = instance.active_components(Objective)[objective_name]
1214         # clone the objective, because we're going to augment (repeatedly) the original objective.
1215         objective_expression = self._original_objective_expression[instance_name].clone() 
1216         # the quadratic expression is really treated as just a list - eventually should be treated as a full expression.
1217         quad_expression = 0.0
1218         
1219         for stage in self._scenario_tree._stages[:-1]: # skip the last stage, as no blending occurs
1220
1221            variable_tree_node = None
1222            for node in stage._tree_nodes:
1223               for scenario in node._scenarios:
1224                  if scenario._name == instance_name:
1225                     variable_tree_node = node
1226                     break
1227
1228            for (reference_variable, index_template, variable_indices) in stage._variables:
1229
1230               variable_name = reference_variable.name
1231               variable_type = reference_variable.domain
1232
1233               w_parameter_name = "PHWEIGHT_"+variable_name
1234               w_parameter = instance.active_components(Param)[w_parameter_name]
1235                 
1236               average_parameter_name = "PHAVG_"+variable_name
1237               average_parameter = instance.active_components(Param)[average_parameter_name]
1238
1239               rho_parameter_name = "PHRHO_"+variable_name
1240               rho_parameter = instance.active_components(Param)[rho_parameter_name]
1241
1242               blend_parameter_name = "PHBLEND_"+variable_name
1243               blend_parameter = instance.active_components(Param)[blend_parameter_name]
1244
1245               node_min_parameter = variable_tree_node._minimums[variable_name]
1246               node_max_parameter = variable_tree_node._maximums[variable_name]               
1247
1248               quad_penalty_term_variable = None
1249               if self._linearize_nonbinary_penalty_terms > 0:
1250                  quad_penalty_term_variable_name = "PHQUADPENALTY_"+variable_name
1251                  quad_penalty_term_variable = instance.active_components(Var)[quad_penalty_term_variable_name]
1252
1253               instance_variable = instance.active_components(Var)[variable_name]
1254
1255               for index in variable_indices:
1256
1257                  if (instance_variable[index].status is not VarStatus.unused) and (instance_variable[index].fixed is False):
1258
1259                     # add the linear (w-weighted) term is a consistent fashion, independent of variable type.
1260                     # if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes for reasons currently unknown.
1261                     objective_expression += (w_parameter[index] * instance_variable[index])
1262
1263                     # there are some cases in which a user may want to avoid the proximal term completely.
1264                     # it's of course only a good idea when there are at least bounds (both lb and ub) on
1265                     # the variables to be blended.
1266                     if self._drop_proximal_terms is False:
1267
1268                        # deal with binaries
1269                        if isinstance(variable_type, BooleanSet) is True:
1270
1271                           if self._retain_quadratic_binary_terms is False:
1272                              # this rather ugly form of the linearized quadratic expression term is required
1273                              # due to a pyomo bug - the form (rho/2) * (x+y+z) chokes in presolve when distributing
1274                              # over the sum.
1275                              new_term = (blend_parameter[index] * rho_parameter[index] / 2.0 * instance_variable[index]) - \
1276                                         (blend_parameter[index] * rho_parameter[index] * average_parameter[index] * instance_variable[index]) + \
1277                                         (blend_parameter[index] * rho_parameter[index] / 2.0 * average_parameter[index] * average_parameter[index])                             
1278                              if objective.sense is minimize:
1279                                 objective_expression += new_term
1280                              else:
1281                                 objective_expression -= new_term                                 
1282                           else:
1283                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)
1284
1285                        # deal with everything else
1286                        else:
1287
1288                           if self._linearize_nonbinary_penalty_terms > 0:
1289
1290                              # the variables are easy - just a single entry.
1291                              if objective.sense is minimize:
1292                                 objective_expression += (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
1293                              else:
1294                                 objective_expression -= (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
1295                                 
1296                              # the constraints are somewhat nastier.
1297
1298                              # TBD - DEFINE CONSTRAINTS ON-THE-FLY?? (INDIVIDUALLY NAMED FOR NOW - CREATE INDEX SETS!) - OR A LEAST AN INDEX SET PER "PIECE"
1299                              xavg = average_parameter[index]
1300                              x = instance_variable[index]
1301
1302                              lb = None
1303                              ub = None
1304
1305                              if x.lb is None:
1306                                 raise ValueError, "No lower bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
1307                              else:
1308                                 lb = x.lb()
1309
1310                              if x.ub is None:
1311                                 raise ValueError, "No upper bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
1312                              else:
1313                                 ub = x.ub()
1314
1315                              node_min = node_min_parameter[index]()
1316                              node_max = node_max_parameter[index]()
1317
1318                              # compute the breakpoint sequence according to the specified strategy.
1319                              breakpoints = []
1320                              if self._breakpoint_strategy == 0:
1321                                 breakpoints = self.compute_uniform_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
1322                              elif self._breakpoint_strategy == 1:
1323                                 breakpoints = self.compute_uniform_between_nodestat_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
1324                              elif self._breakpoint_strategy == 2:
1325                                 breakpoints = self.compute_uniform_between_woodruff_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
1326                              elif self._breakpoint_strategy == 3:
1327                                 breakpoints = self.compute_exponential_from_mean_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)                                                                                                   
1328                              else:
1329                                 raise ValueError, "A breakpoint distribution strategy="+str(self._breakpoint_strategy)+" is currently not supported within PH!"
1330
1331                              for i in range(0,len(breakpoints)-1):
1332
1333                                 this_lb = breakpoints[i]
1334                                 this_ub = breakpoints[i+1]
1335
1336                                 piece_constraint_name = "QUAD_PENALTY_PIECE_"+str(i)+"_"+variable_name+str(index)
1337                                 if hasattr(instance, piece_constraint_name) is False:
1338                                    # this is the first time the constraint is being added - add it to the list of PH-specific constraints for this instance.
1339                                    self._instance_augmented_attributes[instance_name].append(piece_constraint_name)                                                                           
1340                                 piece_constraint = Constraint(name=piece_constraint_name)
1341                                 piece_constraint.model = instance
1342                                 piece_expression = self._create_piecewise_constraint_expression(this_lb, this_ub, x, xavg, quad_penalty_term_variable[index])
1343                                 piece_constraint.add(None, piece_expression)
1344                                 setattr(instance, piece_constraint_name, piece_constraint)                                   
1345
1346                           else:
1347
1348                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)                           
1349                   
1350         # strictly speaking, this probably isn't necessary - parameter coefficients won't get
1351         # pre-processed out of the expression tree. however, if the under-the-hood should change,
1352         # we'll be covered.
1353         objective_expression.simplify(instance)
1354         instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression
1355         # if we are linearizing everything, then nothing will appear in the quadratic expression -
1356         # don't add the empty "0.0" expression to the objective. otherwise, the output file won't
1357         # be properly generated.
1358         if quad_expression != 0.0:
1359           instance.active_components(Objective)[objective_name]._quad_subexpr = quad_expression
1360
1361   def iteration_k_solve(self):
1362
1363      if self._verbose is True:
1364         print "------------------------------------------------"       
1365         print "Starting PH iteration " + str(self._current_iteration) + " solves"
1366
1367      # cache the objective values generated by PH for output at the end of this function.
1368      ph_objective_values = {}
1369
1370      solve_start_time = time.time()
1371
1372      # STEP 0: set up all global solver options.
1373      self._solver.mipgap = self._mipgap     
1374
1375      # STEP 1: queue up the solves for all scenario sub-problems and
1376      #         grab all of the action handles for the subsequent barrier sync.
1377
1378      action_handles = []
1379      scenario_action_handle_map = {} # maps scenario names to action handles
1380      action_handle_scenario_map = {} # maps action handles to scenario names
1381
1382      for scenario in self._scenario_tree._scenarios:     
1383
1384         instance = self._instances[scenario._name]
1385
1386         if self._verbose is True:
1387            print "Queuing solve for scenario=" + scenario._name
1388
1389         # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
1390         # don't do this, you won't see any chance in the output files as you vary the problem parameters!
1391         # ditto for instance fixing!
1392         instance.presolve()
1393
1394         # once past iteration 0, there is always a feasible solution from which to warm-start.
1395         # however, you might want to disable warm-start when the solver is behaving badly (which does happen).
1396         new_action_handle = None
1397         if (self._disable_warmstarts is False) and (self._solver.warm_start_capable() is True):
1398            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
1399         else:
1400            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)           
1401
1402         scenario_action_handle_map[scenario._name] = new_action_handle
1403         action_handle_scenario_map[new_action_handle] = scenario._name         
1404
1405         action_handles.append(new_action_handle)
1406
1407      # STEP 2: loop for the solver results, reading them and loading
1408      #         them into instances as they are available.
1409      if self._verbose is True:
1410         print "Waiting for scenario sub-problem solves"
1411
1412      num_results_so_far = 0
1413
1414      while (num_results_so_far < len(self._scenario_tree._scenarios)):
1415
1416         action_handle = self._solver_manager.wait_any()
1417         results = self._solver_manager.get_results(action_handle)         
1418         scenario_name = action_handle_scenario_map[action_handle]
1419         instance = self._instances[scenario_name]         
1420
1421         if self._verbose is True:
1422            print "Results obtained for scenario="+scenario_name
1423
1424         if len(results.solution) == 0:
1425            results.write(num=1)
1426            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
1427
1428         if self._output_solver_results is True:
1429            print "Results for scenario=",scenario_name
1430            results.write(num=1)
1431
1432         start_time = time.time()
1433         instance.load(results)
1434         end_time = time.time()
1435         if self._output_times is True:
1436            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
1437
1438         if self._verbose is True:                 
1439            print "Successfully loaded solution for scenario="+scenario_name
1440
1441         # we're assuming there is a single solution.
1442         # the "value" attribute is a pre-defined feature of any solution - it is relative to whatever
1443         # objective was selected during optimization, which of course should be the PH objective.
1444         ph_objective_values[instance.name] = float(results.solution(0).objective['f'].value)
1445
1446         num_results_so_far = num_results_so_far + 1
1447         
1448      if self._verbose is True:     
1449         print "Scenario sub-problem solves completed"     
1450
1451      solve_end_time = time.time()
1452      self._cumulative_solve_time += (solve_end_time - solve_start_time)
1453
1454      if self._output_times is True:
1455         print "Aggregate sub-problem solve time this iteration=%8.2f" % (solve_end_time - solve_start_time)
1456
1457      if self._verbose is True:
1458         print "Successfully completed PH iteration " + str(self._current_iteration) + " solves - solution statistics:"
1459         print "  Scenario             PH Objective             Cost Objective"
1460         for scenario in self._scenario_tree._scenarios:
1461            instance = self._instances[scenario._name]
1462            for objective_name in instance.active_components(Objective):
1463               objective = instance.active_components(Objective)[objective_name]
1464               print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], 0.0)
1465
1466   def solve(self):
1467
1468      self._solve_start_time = time.time()
1469      self._cumulative_solve_time = 0.0
1470      self._cumulative_xbar_time = 0.0
1471      self._cumulative_weight_time = 0.0
1472
1473      print "Starting PH"
1474
1475      if self._initialized == False:
1476         raise RuntimeError, "PH is not initialized - cannot invoke solve() method"
1477
1478      print "Initiating PH iteration=" + `self._current_iteration`
1479
1480      self.iteration_0_solve()
1481
1482      # update variable statistics prior to any output.
1483      self.update_variable_statistics()     
1484
1485      if (self._verbose is True) or (self._report_solutions is True):
1486         print "Variable values following scenario solves:"
1487         self.pprint(False,False,True,False)
1488
1489      # let plugins know if they care.
1490      for plugin in self._ph_plugins:     
1491         plugin.post_iteration_0_solves(self)
1492
1493      # update the fixed variable statistics.
1494      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
1495
1496      if self._verbose is True:
1497         print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1498         print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1499
1500      # always output the convergence metric and first-stage cost statistics, to give a sense of progress.
1501      self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1502      first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()
1503      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)               
1504
1505      self.update_weights()
1506
1507      # let plugins know if they care.
1508      for plugin in self._ph_plugins:     
1509         plugin.post_iteration_0(self)
1510
1511      # checkpoint if it's time - which it always is after iteration 0,
1512      # if the interval is >= 1!
1513      if (self._checkpoint_interval > 0):
1514         self.checkpoint(0)
1515
1516      # there is an upper bound on the number of iterations to execute -
1517      # the actual bound depends on the converger supplied by the user.
1518      for i in range(1, self._max_iterations+1):
1519
1520         self._current_iteration = self._current_iteration + 1                 
1521
1522         print "Initiating PH iteration=" + `self._current_iteration`         
1523
1524         if (self._verbose is True) or (self._report_weights is True):
1525            print "Variable averages and weights prior to scenario solves:"
1526            self.pprint(True,True,False,False)
1527
1528         # with the introduction of piecewise linearization, the form of the
1529         # penalty-weighted objective is no longer fixed. thus, we need to
1530         # create the objectives each PH iteration.
1531         self.form_iteration_k_objectives()
1532
1533         # let plugins know if they care.
1534         for plugin in self._ph_plugins:
1535            plugin.pre_iteration_k_solves(self)
1536
1537         # do the actual solves.
1538         self.iteration_k_solve()
1539
1540         # update variable statistics prior to any output.
1541         self.update_variable_statistics()
1542         
1543         if (self._verbose is True) or (self._report_solutions is True):
1544            print "Variable values following scenario solves:"
1545            self.pprint(False,False,True,False)
1546
1547         # we don't technically have to do this at the last iteration,
1548         # but with checkpointing and re-starts, you're never sure
1549         # when you're executing the last iteration.
1550         self.update_weights()
1551
1552         # let plugins know if they care.
1553         for plugin in self._ph_plugins:
1554            plugin.post_iteration_k_solves(self)
1555
1556         # update the fixed variable statistics.
1557         (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
1558
1559         if self._verbose is True:
1560            print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1561            print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1562
1563         # let plugins know if they care.
1564         for plugin in self._ph_plugins:
1565            plugin.post_iteration_k(self)
1566
1567         # at this point, all the real work of an iteration is complete.
1568
1569         # checkpoint if it's time.
1570         if (self._checkpoint_interval > 0) and (i % self._checkpoint_interval is 0):
1571            self.checkpoint(i)
1572
1573         # check for early termination.
1574         self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1575         first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()         
1576         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)         
1577
1578         if self._converger.isConverged(self) is True:
1579            if self._total_discrete_vars == 0:
1580               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)
1581            else:
1582               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)+" or all discrete variables are fixed"               
1583            break
1584
1585         # if we're terminating due to exceeding the maximum iteration count, print a message
1586         # indicating so - otherwise, you get a quiet, information-free output trace.
1587         if i == self._max_iterations:
1588            print "Halting PH - reached maximal iteration count="+str(self._max_iterations)
1589
1590      if self._verbose is True:
1591         print "Number of discrete variables fixed before final plugin calls="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1592         print "Number of continuous variables fixed before final plugin calls="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"           
1593
1594      # let plugins know if they care. do this before
1595      # the final solution / statistics output, as the plugins
1596      # might do some final tweaking.
1597      for plugin in self._ph_plugins:
1598         plugin.post_ph_execution(self)
1599
1600      # update the fixed variable statistics - the plugins might have done something.
1601      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()                       
1602
1603      print "PH complete"
1604
1605      print "Convergence history:"
1606      self._converger.pprint()
1607
1608      print "Final number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1609      print "Final number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"         
1610
1611      print "Final variable values:"
1612      self.pprint(False,False,True,True)         
1613
1614      print "Final costs:"
1615      self._scenario_tree.pprintCosts(self._instances)
1616
1617      self._solve_end_time = time.time()
1618
1619      if (self._verbose is True) and (self._output_times is True):
1620         print "Overall run-time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1621
1622      # cleanup the scenario instances for post-processing - ideally, we want to leave them in
1623      # their original state, minus all the PH-specific stuff. we don't do all cleanup (leaving
1624      # things like rhos, etc), but we do clean up constraints, as that really hoses up the ef writer.
1625      self._cleanup_scenario_instances()
1626
1627   #
1628   # prints a summary of all collected time statistics
1629   #
1630   def print_time_stats(self):
1631
1632      print "PH run-time statistics (user):"
1633
1634      print "Initialization time=  %8.2f seconds" % (self._init_end_time - self._init_start_time)
1635      print "Overall solve time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1636      print "Scenario solve time=  %8.2f seconds" % self._cumulative_solve_time
1637      print "Average update time=  %8.2f seconds" % self._cumulative_xbar_time
1638      print "Weight update time=   %8.2f seconds" % self._cumulative_weight_time
1639
1640   #
1641   # a utility to determine whether to output weight / average / etc. information for
1642   # a variable/node combination. when the printing is moved into a callback/plugin,
1643   # this routine will go there. for now, we don't dive down into the node resolution -
1644   # just the variable/stage.
1645   #
1646   def should_print(self, stage, variable, variable_indices):
1647
1648      if self._output_continuous_variable_stats is False:
1649
1650         variable_type = variable.domain         
1651
1652         if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):
1653
1654            return False
1655
1656      return True
1657     
1658   #
1659   # pretty-prints the state of the current variable averages, weights, and values.
1660   # inputs are booleans indicating which components should be output.
1661   #
1662   def pprint(self, output_averages, output_weights, output_values, output_fixed):
1663
1664      if self._initialized is False:
1665         raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"         
1666     
1667      # print tree nodes and associated variable/xbar/ph information in stage-order
1668      # we don't blend in the last stage, so we don't current care about printing the associated information.
1669      for stage in self._scenario_tree._stages[:-1]:
1670
1671         print "\tStage=" + stage._name
1672
1673         num_outputs_this_stage = 0 # tracks the number of outputs on a per-index basis.
1674
1675         for (variable, index_template, variable_indices) in stage._variables:
1676
1677            variable_name = variable.name
1678
1679            if self.should_print(stage, variable, variable_indices) is True:
1680
1681               num_outputs_this_variable = 0 # track, so we don't output the variable names unless there is an entry to report.
1682
1683               for index in variable_indices:               
1684
1685                  weight_parameter_name = "PHWEIGHT_"+variable_name
1686
1687                  num_outputs_this_index = 0 # track, so we don't output the variable index more than once.
1688
1689                  for tree_node in stage._tree_nodes:                 
1690
1691                     # determine if the variable/index pair is used across the set of scenarios (technically,
1692                     # it should be good enough to check one scenario). ditto for "fixed" status. fixed does
1693                     # imply unused (see note below), but we care about the fixed status when outputting
1694                     # final solutions.
1695
1696                     is_used = True # should be consistent across scenarios, so one "unused" flags as invalid.
1697                     is_fixed = False
1698
1699                     for scenario in tree_node._scenarios:
1700                        instance = self._instances[scenario._name]
1701                        variable_value = getattr(instance,variable_name)[index]
1702                        if variable_value.status == VarStatus.unused:
1703                           is_used = False
1704                        if variable_value.fixed is True:
1705                           is_fixed = True
1706 
1707                     # IMPT: this is far from obvious, but variables that are fixed will - because
1708                     #       presolve will identify them as constants and eliminate them from all
1709                     #       expressions - be flagged as "unused" and therefore not output.
1710
1711                     if ((output_fixed is True) and (is_fixed is True)) or (is_used is True):
1712
1713                           minimum_value = None
1714                           maximum_value = None
1715
1716                           if index is None:
1717                              minimum_value = tree_node._minimums[variable_name]
1718                              maximum_value = tree_node._maximums[variable_name]
1719                           else:
1720                              minimum_value = tree_node._minimums[variable_name][index]()
1721                              maximum_value = tree_node._maximums[variable_name][index]()
1722
1723                           # there really isn't a need to output variables whose
1724                           # values are equal to 0 across-the-board. and there is
1725                           # good reason not to, i.e., the volume of output.
1726                           if (fabs(minimum_value) > self._integer_tolerance) or \
1727                              (fabs(maximum_value) > self._integer_tolerance):
1728
1729                              num_outputs_this_stage = num_outputs_this_stage + 1                           
1730                              num_outputs_this_variable = num_outputs_this_variable + 1
1731                              num_outputs_this_index = num_outputs_this_index + 1
1732
1733                              if num_outputs_this_variable == 1:
1734                                 print "\t\tVariable=" + variable_name
1735
1736                              if num_outputs_this_index == 1:
1737                                 if index is not None:
1738                                    print "\t\t\tIndex:", indexToString(index)                             
1739
1740                              print "\t\t\t\tTree Node="+tree_node._name+"\t\t (Scenarios: ",                             
1741                              for scenario in tree_node._scenarios:
1742                                 print scenario._name," ",
1743                                 if scenario == tree_node._scenarios[-1]:
1744                                    print ")"
1745                           
1746                              if output_values is True:
1747                                 average_value = tree_node._averages[variable_name][index]()
1748                                 print "\t\t\t\tValues: ",                       
1749                                 for scenario in tree_node._scenarios:
1750                                    instance = self._instances[scenario._name]
1751                                    this_value = getattr(instance,variable_name)[index].value
1752                                    print "%12.4f" % this_value,
1753                                    if scenario == tree_node._scenarios[-1]:
1754                                       print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1755                                       print "    Avg=%12.4f" % (average_value),
1756                                       print ""
1757                              if output_weights:
1758                                 print "\t\t\t\tWeights: ",
1759                                 for scenario in tree_node._scenarios:
1760                                    instance = self._instances[scenario._name]
1761                                    print "%12.4f" % getattr(instance,weight_parameter_name)[index].value,
1762                                    if scenario == tree_node._scenarios[-1]:
1763                                       print ""
1764
1765                              if output_averages:
1766                                 print "\t\t\t\tAverage: %12.4f" % (tree_node._averages[variable_name][index].value)
1767
1768         if num_outputs_this_stage == 0:
1769            print "\t\tNo non-converged variables in stage"
1770
1771         # cost variables aren't blended, so go through the gory computation of min/max/avg.
1772         # we currently always print these.
1773         cost_variable_name = stage._cost_variable[0].name
1774         cost_variable_index = stage._cost_variable[1]
1775         if cost_variable_index is None:
1776            print "\t\tCost Variable=" + cost_variable_name
1777         else:
1778            print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)           
1779         for tree_node in stage._tree_nodes:
1780            print "\t\t\tTree Node=" + tree_node._name + "\t\t (Scenarios: ",
1781            for scenario in tree_node._scenarios:
1782               print scenario._name," ",
1783               if scenario == tree_node._scenarios[-1]:
1784                  print ")"
1785            maximum_value = 0.0
1786            minimum_value = 0.0
1787            sum_values = 0.0
1788            num_values = 0
1789            first_time = True
1790            print "\t\t\tValues: ",                       
1791            for scenario in tree_node._scenarios:
1792                instance = self._instances[scenario._name]
1793                this_value = getattr(instance,cost_variable_name)[cost_variable_index].value
1794                print "%12.4f" % this_value,
1795                num_values += 1
1796                sum_values += this_value
1797                if first_time is True:
1798                   first_time = False
1799                   maximum_value = this_value
1800                   minimum_value = this_value
1801                else:
1802                   if this_value > maximum_value:
1803                      maximum_value = this_value
1804                   if this_value < minimum_value:
1805                      minimum_value = this_value
1806                if scenario == tree_node._scenarios[-1]:
1807                   print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1808                   print "    Avg=%12.4f" % (sum_values/num_values),
1809                   print ""
1810           
Note: See TracBrowser for help on using the repository browser.