source: coopr.pysp/stable/coopr/pysp/ph.py @ 3185

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

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

........

r3120 | jwatson | 2010-10-18 21:11:21 -0600 (Mon, 18 Oct 2010) | 3 lines


Adding PySP options for linearizing expressions.

........

r3134 | jwatson | 2010-10-21 15:45:49 -0600 (Thu, 21 Oct 2010) | 3 lines


Bug fix associated with linear expressions.

........

r3137 | khunter | 2010-10-22 10:19:44 -0600 (Fri, 22 Oct 2010) | 2 lines


Add name as contributor.

........

r3138 | jwatson | 2010-10-22 14:11:43 -0600 (Fri, 22 Oct 2010) | 3 lines


Added "--simplify-expressions" option to runph, in order to eliminate the memory and time costs associated with simplifying expressions (e.g., in formulation of the PH objective) for which simpification is very unlikely to help. By default, expression simplification is disabled. This may cause issues with certain writers, e.g., NL - which is why I have retained the option.

........

r3139 | jwatson | 2010-10-22 15:20:49 -0600 (Fri, 22 Oct 2010) | 3 lines


When forming the PH linear terms for all variables and the proximal terms for binary variables, use value() to immediately extract the fixed value of the underlying parameters. This speeds up the associated expression tree construction time, which was non-trivial (5% of the total run-time, for example). Because we re-form the expressions each iteration, this is OK.

........

r3143 | jwatson | 2010-10-22 22:07:07 -0600 (Fri, 22 Oct 2010) | 3 lines


Modifying more PH parameters to be of the "nochecking" variety. Also slightly improved the performance of the update_variable_statistics PH function.

........

r3145 | jwatson | 2010-10-22 22:21:06 -0600 (Fri, 22 Oct 2010) | 3 lines


Update of PySP baseline output files, due to small numerical differences (4th significant digit) associated with some code mods I recently performed in the course of optimization.

........

r3146 | jwatson | 2010-10-23 13:27:38 -0600 (Sat, 23 Oct 2010) | 3 lines


Speed improvements in weight update.

........

r3147 | jwatson | 2010-10-23 13:38:52 -0600 (Sat, 23 Oct 2010) | 3 lines


More speed improvements to PH weight computation procedure.

........

r3152 | jwatson | 2010-10-24 13:20:03 -0600 (Sun, 24 Oct 2010) | 3 lines


Further code optimizations to PH weight and statistic update routines.

........

r3155 | jwatson | 2010-10-24 22:15:29 -0600 (Sun, 24 Oct 2010) | 3 lines


Added --linearize-expression and associated processing to runef PySP script.

........

r3156 | jwatson | 2010-10-24 22:19:44 -0600 (Sun, 24 Oct 2010) | 3 lines


Correcting shift from "integer" to "general" variable blocks when writing the LP file for an extensive form.

........

r3165 | jwatson | 2010-10-26 18:46:24 -0600 (Tue, 26 Oct 2010) | 3 lines


Update of PySP CHANGELOG for 2.4 release.

........

r3168 | jwatson | 2010-10-26 20:12:23 -0600 (Tue, 26 Oct 2010) | 3 lines


Completing implementation of optional processing of expression simplification in PySP.

........

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