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

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

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

........

r3201 | jwatson | 2010-10-29 13:18:17 -0600 (Fri, 29 Oct 2010) | 3 lines


Inverting order of .dat files in PySP when loading from a node representation - now root-to-leaf, instead of leaf-to-root. This allows for deeper-in-the-tree nodes to over-write parameter values defined higher in the tree, which is a more "expected" behavior than the converse. The real answer is to throw an exception if a parameter is re-defined, but we're not there yet.

........

r3217 | jwatson | 2010-11-05 11:29:42 -0600 (Fri, 05 Nov 2010) | 3 lines


Various updates to support heteogeneous index sets in PH for different nodes in the scenario tree - more work / testing remains.

........

r3218 | jwatson | 2010-11-05 12:35:51 -0600 (Fri, 05 Nov 2010) | 3 lines


More changes associated with generalizing the PySP index structures from per-stage to per-node.

........

r3220 | jwatson | 2010-11-05 20:28:29 -0600 (Fri, 05 Nov 2010) | 3 lines


Various fixes to the WW PH extension, bringing it in compliance to previous commit changes.

........

r3221 | jwatson | 2010-11-05 20:43:40 -0600 (Fri, 05 Nov 2010) | 1 line


Removing some older PySP test problems from the repository

........

r3222 | jwatson | 2010-11-05 20:55:15 -0600 (Fri, 05 Nov 2010) | 1 line


Moving PySP forestry example to local sandbox, to streamline distribution.

........

r3261 | jwatson | 2010-11-29 15:26:46 -0700 (Mon, 29 Nov 2010) | 9 lines


Adding two options to the runef and runph pysp scripts, to facilitate scenario downsampling - the case where you have a big tree, but you don't want to use it all.


The options are:
--scenario-tree-downsample-fraction=X
--scenario-tree-random-seed


The options are fairly self-explanatory - the only possible nuance is that the downsample fraction is the fraction of scenarios retained.

........

r3263 | jwatson | 2010-12-01 10:47:21 -0700 (Wed, 01 Dec 2010) | 3 lines


Corrected issue with cvar generation introduced a while back.

........

r3264 | jwatson | 2010-12-01 11:16:01 -0700 (Wed, 01 Dec 2010) | 3 lines


Adding PySP extensive form tests.

........

r3265 | wehart | 2010-12-01 13:19:56 -0700 (Wed, 01 Dec 2010) | 2 lines


Auxmenting the filter

........

r3266 | wehart | 2010-12-01 13:58:59 -0700 (Wed, 01 Dec 2010) | 2 lines


Adding further diagnostics to the filter.

........

r3267 | wehart | 2010-12-01 14:12:15 -0700 (Wed, 01 Dec 2010) | 2 lines


Another attempt to fix this filter...

........

r3271 | wehart | 2010-12-01 15:51:23 -0700 (Wed, 01 Dec 2010) | 4 lines


Simplifying filter.


The error was really a Python 2.5 portability issue in EF. :(

........

r3275 | jwatson | 2010-12-02 19:35:52 -0700 (Thu, 02 Dec 2010) | 3 lines


Fixing Python 2.5-related issue with string.translate in phutils.py - using None as the translation table is not allowed in Python 2.5.

........

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