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

Last change on this file since 3104 was 3096, checked in by jwatson, 11 years ago

Enhancements to PySP to deal with reporting issues caused by the lack of objective function values in SOL files (in general, and those produced by ipopt in particular).

File size: 72.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, 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      # we only store these temporarily...
462      scenario_solver_options = None
463
464      # process the keyword options
465      for key in kwds.keys():
466         if key == "max_iterations":
467            self._max_iterations = kwds[key]
468         elif key == "rho":
469            self._rho = kwds[key]
470         elif key == "rho_setter":
471            self._rho_setter = kwds[key]
472         elif key == "bounds_setter":
473            self._bounds_setter = kwds[key]
474         elif key == "solver":
475            self._solver_type = kwds[key]
476         elif key == "solver_manager":
477            self._solver_manager_type = kwds[key]
478         elif key == "scenario_solver_options":
479            scenario_solver_options = kwds[key]
480         elif key == "scenario_mipgap":
481            self._mipgap = kwds[key]
482         elif key == "keep_solver_files":
483            self._keep_solver_files = kwds[key]
484         elif key == "output_solver_results":
485            self._output_solver_results = kwds[key]
486         elif key == "output_solver_log":
487            self._output_solver_log = kwds[key]
488         elif key == "verbose":
489            self._verbose = kwds[key]
490         elif key == "report_solutions":
491            self._report_solutions = kwds[key]
492         elif key == "report_weights":
493            self._report_weights = kwds[key]
494         elif key == "report_only_statistics":
495            self._report_only_statistics = kwds[key]
496         elif key == "output_times":
497            self._output_times = kwds[key]
498         elif key == "disable_warmstarts":
499            self._disable_warmstarts = kwds[key]
500         elif key == "drop_proximal_terms":
501            self._drop_proximal_terms = kwds[key]
502         elif key == "retain_quadratic_binary_terms":
503            self._retain_quadratic_binary_terms = kwds[key]
504         elif key == "linearize_nonbinary_penalty_terms":
505            self._linearize_nonbinary_penalty_terms = kwds[key]
506         elif key == "breakpoint_strategy":
507            self._breakpoint_strategy = kwds[key]
508         elif key == "checkpoint_interval":
509            self._checkpoint_interval = kwds[key]
510         elif key == "output_scenario_tree_solution":
511            self._output_scenario_tree_solution = kwds[key]
512         else:
513            print "Unknown option=" + key + " specified in call to PH constructor"
514
515      # validate all "atomic" options (those that can be validated independently)
516      if self._max_iterations < 0:
517         raise ValueError, "Maximum number of PH iterations must be non-negative; value specified=" + `self._max_iterations`
518      if self._rho <= 0.0:
519         raise ValueError, "Value of the rho parameter in PH must be non-zero positive; value specified=" + `self._rho`
520      if (self._mipgap is not None) and ((self._mipgap < 0.0) or (self._mipgap > 1.0)):
521         raise ValueError, "Value of the mipgap parameter in PH must be on the unit interval; value specified=" + `self._mipgap`
522
523      # validate the linearization (number of pieces) and breakpoint distribution parameters.
524      if self._linearize_nonbinary_penalty_terms < 0:
525         raise ValueError, "Value of linearization parameter for nonbinary penalty terms must be non-negative; value specified=" + `self._linearize_nonbinary_penalty_terms`
526      if self._breakpoint_strategy < 0:
527         raise ValueError, "Value of the breakpoint distribution strategy parameter must be non-negative; value specified=" + str(self._breakpoint_strategy)
528      if self._breakpoint_strategy > 3:
529         raise ValueError, "Unknown breakpoint distribution strategy specified - valid values are between 0 and 2, inclusive; value specified=" + str(self._breakpoint_strategy)
530
531      # validate rho setter file if specified.
532      if self._rho_setter is not None:
533         if os.path.exists(self._rho_setter) is False:
534            raise ValueError, "The rho setter script file="+self._rho_setter+" does not exist"
535
536      # validate bounds setter file if specified.
537      if self._bounds_setter is not None:
538         if os.path.exists(self._bounds_setter) is False:
539            raise ValueError, "The bounds setter script file="+self._bounds_setter+" does not exist"
540
541      # validate the checkpoint interval.
542      if self._checkpoint_interval < 0:
543         raise ValueError, "A negative checkpoint interval with value="+str(self._checkpoint_interval)+" was specified in call to PH constructor"
544
545      # construct the sub-problem solver.
546      if self._verbose is True:
547         print "Constructing solver type="+self._solver_type
548      self._solver = SolverFactory(self._solver_type)
549      if self._solver == None:
550         raise ValueError, "Unknown solver type=" + self._solver_type + " specified in call to PH constructor"
551      if self._keep_solver_files is True:
552         self._solver.keepFiles = True
553      if len(scenario_solver_options) > 0:
554         if self._verbose is True:
555            print "Initializing scenario sub-problem solver with options="+str(scenario_solver_options)
556         self._solver.set_options("".join(scenario_solver_options))
557      if self._output_times is True:
558         self._solver._report_timing = True
559
560      # construct the solver manager.
561      if self._verbose is True:
562         print "Constructing solver manager of type="+self._solver_manager_type
563      self._solver_manager = SolverManagerFactory(self._solver_manager_type)
564      if self._solver_manager is None:
565         raise ValueError, "Failed to create solver manager of type="+self._solver_manager_type+" specified in call to PH constructor"
566
567      # a set of all valid PH iteration indicies is generally useful for plug-ins, so create it here.
568      self._iteration_index_set = Set(name="PHIterations")
569      for i in range(0,self._max_iterations + 1):
570         self._iteration_index_set.add(i)
571
572      # spit out parameterization if verbosity is enabled
573      if self._verbose is True:
574         print "PH solver configuration: "
575         print "   Max iterations=" + `self._max_iterations`
576         print "   Default global rho=" + `self._rho`
577         if self._rho_setter is not None:
578            print "   Rho initialization file=" + self._rho_setter
579         if self._bounds_setter is not None:
580            print "   Variable bounds initialization file=" + self._bounds_setter
581         print "   Sub-problem solver type=" + `self._solver_type`
582         print "   Solver manager type=" + `self._solver_manager_type`
583         print "   Keep solver files? " + str(self._keep_solver_files)
584         print "   Output solver results? " + str(self._output_solver_results)
585         print "   Output solver log? " + str(self._output_solver_log)
586         print "   Output times? " + str(self._output_times)
587         print "   Checkpoint interval="+str(self._checkpoint_interval)
588
589   """ Initialize PH with model and scenario data, in preparation for solve().
590       Constructs and reads instances.
591   """
592   def initialize(self, scenario_data_directory_name=".", model=None, model_instance=None, scenario_tree=None, converger=None):
593
594      self._init_start_time = time.time()
595
596      if self._verbose is True:
597         print "Initializing PH"
598         print "   Scenario data directory=" + scenario_data_directory_name
599
600      if not os.path.exists(scenario_data_directory_name):
601         raise ValueError, "Scenario data directory=" + scenario_data_directory_name + " either does not exist or cannot be read"
602
603      self._scenario_data_directory_name = scenario_data_directory_name
604
605      # IMPT: The input model should be an *instance*, as it is very useful (critical!) to know
606      #       the dimensions of sets, be able to store suffixes on variable values, etc.
607      if model is None:
608         raise ValueError, "A model must be supplied to the PH initialize() method"
609
610      if scenario_tree is None:
611         raise ValueError, "A scenario tree must be supplied to the PH initialize() method"
612
613      if converger is None:
614         raise ValueError, "A convergence computer must be supplied to the PH initialize() method"
615
616      self._model = model
617      self._model_instance = model_instance
618      self._scenario_tree = scenario_tree
619      self._converger = converger
620
621      model_objective = model.active_components(Objective)
622      self._is_minimizing = (model_objective[ model_objective.keys()[0] ].sense == minimize)
623
624      self._converger.reset()
625
626      # construct the instances for each scenario.
627      #
628      # garbage collection noticeably slows down PH when dealing with
629      # large numbers of scenarios. disable prior to instance construction,
630      # and then re-enable. there isn't much collection to do as instances
631      # are constructed.
632      re_enable_gc = gc.isenabled()
633      gc.disable()
634
635      if self._verbose is True:
636         if self._scenario_tree._scenario_based_data == 1:
637            print "Scenario-based instance initialization enabled"
638         else:
639            print "Node-based instance initialization enabled"
640
641      for scenario in self._scenario_tree._scenarios:
642
643         scenario_instance = construct_scenario_instance(self._scenario_tree,
644                                                         self._scenario_data_directory_name,
645                                                         scenario._name,
646                                                         self._model,
647                                                         self._verbose,
648                                                         preprocess=False)
649
650         # IMPT: disable canonical representation construction for ASL solvers.
651         #       this is a hack, in that we need to address encodings and
652         #       the like at a more general level.
653         if self._solver_type == "asl":
654            scenario_instance.skip_canonical_repn = True
655         else:
656            scenario_instance.preprocess()
657
658         self._instances[scenario._name] = scenario_instance
659         self._instances[scenario._name].name = scenario._name
660         self._instance_augmented_attributes[scenario._name] = []
661
662
663      # perform a single pass of garbage collection and re-enable automatic collection.
664      if re_enable_gc is True:
665         gc.collect()
666         gc.enable()
667
668      # let plugins know if they care - this callback point allows
669      # users to create/modify the original scenario instances and/or
670      # the scenario tree prior to creating PH-related parameters,
671      # variables, and the like.
672      for plugin in self._ph_plugins:
673         plugin.post_instance_creation(self)
674
675      # create ph-specific parameters (weights, xbar, etc.) for each instance.
676
677      if self._verbose is True:
678         print "Creating weight, average, and rho parameter vectors for scenario instances"
679
680      self._create_ph_scenario_parameters()
681
682      # if specified, run the user script to initialize variable rhos at their whim.
683      if self._rho_setter is not None:
684         print "Executing user rho set script from filename="+self._rho_setter
685         execfile(self._rho_setter)
686
687      # with the instances created, run the user script to initialize variable bounds.
688      if self._bounds_setter is not None:
689         print "Executing user variable bounds set script from filename=", self._bounds_setter
690         execfile(self._bounds_setter)
691
692      # create parameters to store variable statistics (of general utility) at each node in the scenario tree.
693
694      if self._verbose is True:
695         print "Creating variable statistic (min/avg/max) parameter vectors for scenario tree nodes"
696
697      # do this for all stages, simply for completeness, i.e., to create a fully populated scenario tree.
698      for stage in self._scenario_tree._stages:
699
700         # first, gather all unique variables referenced in this stage
701         # this "gather" step is currently required because we're being lazy
702         # in terms of index management in the scenario tree - which
703         # should really be done in terms of sets of indices.
704
705         stage_variables = {}
706         for (reference_variable, index_template, reference_index) in stage._variables:
707            if reference_variable.name not in stage_variables.keys():
708               stage_variables[reference_variable.name] = reference_variable
709
710         # next, create min/avg/max parameters for each variable in the corresponding tree node.
711         # NOTE: the parameter names below could really be empty, as they are never referenced
712         #       explicitly.
713         for (variable_name, reference_variable) in stage_variables.items():
714            for tree_node in stage._tree_nodes:
715
716               new_min_index = reference_variable._index
717               new_min_parameter_name = "NODEMIN_"+reference_variable.name
718               # this bit of ugliness is due to Pyomo not correctly handling the Param construction
719               # case when the supplied index set consists strictly of None, i.e., the source variable
720               # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
721               new_min_parameter = None
722               if (len(new_min_index) is 1) and (None in new_min_index):
723                  new_min_parameter = Param(name=new_min_parameter_name, mutable=True)
724               else:
725                  new_min_parameter = Param(new_min_index, name=new_min_parameter_name, mutable=True)
726               for index in new_min_index:
727                  new_min_parameter[index] = 0.0
728               tree_node._minimums[reference_variable.name] = new_min_parameter
729
730               new_avg_index = reference_variable._index
731               new_avg_parameter_name = "NODEAVG_"+reference_variable.name
732               new_avg_parameter = None
733               if (len(new_avg_index) is 1) and (None in new_avg_index):
734                  new_avg_parameter = Param(name=new_avg_parameter_name, mutable=True)
735               else:
736                  new_avg_parameter = Param(new_avg_index, name=new_avg_parameter_name, mutable=True)
737               for index in new_avg_index:
738                  new_avg_parameter[index] = 0.0
739               tree_node._averages[reference_variable.name] = new_avg_parameter
740
741               new_max_index = reference_variable._index
742               new_max_parameter_name = "NODEMAX_"+reference_variable.name
743               new_max_parameter = None
744               if (len(new_max_index) is 1) and (None in new_max_index):
745                  new_max_parameter = Param(name=new_max_parameter_name, mutable=True)
746               else:
747                  new_max_parameter = Param(new_max_index, name=new_max_parameter_name, mutable=True)
748               for index in new_max_index:
749                  new_max_parameter[index] = 0.0
750               tree_node._maximums[reference_variable.name] = new_max_parameter
751
752      # the objective functions are modified throughout the course of PH iterations.
753      # save the original, as a baseline to modify in subsequent iterations. reserve
754      # the original objectives, for subsequent modification.
755      self._original_objective_expression = {}
756      for instance_name, instance in self._instances.items():
757         objective_name = instance.active_components(Objective).keys()[0]
758         expr = instance.active_components(Objective)[objective_name]._data[None].expr
759         if isinstance(expr, Expression) is False:
760            expr = _IdentityExpression(expr)
761         self._original_objective_expression[instance_name] = expr
762
763      # cache the number of discrete and continuous variables in the master instance. this value
764      # is of general use, e.g., in the converger classes and in plugins.
765      (self._total_discrete_vars,self._total_continuous_vars) = self.compute_variable_counts()
766      if self._verbose is True:
767         print "Total number of discrete instance variables="+str(self._total_discrete_vars)
768         print "Total number of continuous instance variables="+str(self._total_continuous_vars)
769
770      # track the total number of fixed variables of each category at the end of each PH iteration.
771      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
772
773      # indicate that we're ready to run.
774      self._initialized = True
775
776      if self._verbose is True:
777         print "PH successfully created model instances for all scenarios"
778
779      self._init_end_time = time.time()
780
781      if self._verbose is True:
782         print "PH is successfully initialized"
783         if self._output_times is True:
784            print "Initialization time=%8.2f seconds" % (self._init_end_time - self._init_start_time)
785
786      # let plugins know if they care.
787      if self._verbose is True:
788         print "Initializing PH plugins"
789      for plugin in self._ph_plugins:
790         plugin.post_ph_initialization(self)
791      if self._verbose is True:
792         print "PH plugin initialization complete"
793
794   """ Perform the non-weighted scenario solves and form the initial w and xbars.
795   """
796   def iteration_0_solve(self):
797
798      if self._verbose is True:
799         print "------------------------------------------------"
800         print "Starting PH iteration 0 solves"
801
802      self._current_iteration = 0
803
804      solve_start_time = time.time()
805
806      # STEP 0: set up all global solver options.
807      self._solver.mipgap = self._mipgap
808
809      # STEP 1: queue up the solves for all scenario sub-problems and
810      #         grab all the action handles for the subsequent barrier sync.
811
812      action_handles = []
813      scenario_action_handle_map = {} # maps scenario names to action handles
814      action_handle_scenario_map = {} # maps action handles to scenario names
815
816      for scenario in self._scenario_tree._scenarios:
817
818         instance = self._instances[scenario._name]
819
820         if self._verbose is True:
821            print "Queuing solve for scenario=" + scenario._name
822
823         # IMPT: You have to re-presolve if approximating continuous variable penalty terms with a
824         #       piecewise linear function. otherwise, the newly introduced variables won't be flagged
825         #       as unused (as is correct for iteration 0), and the output writer will crater.
826         # IMPT: I decided to presolve unconditionally, as PH extensions can add arbitrary components
827         #       to the base scenario instances - and the variable values/etc. need to be collectged.
828         instance.preprocess()
829
830         # there's nothing to warm-start from in iteration 0, so don't include the keyword in the solve call.
831         # the reason you don't want to include it is that some solvers don't know how to handle the keyword
832         # at all (despite it being false). you might want to solve iteration 0 solves using some other solver.
833
834         new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
835         scenario_action_handle_map[scenario._name] = new_action_handle
836         action_handle_scenario_map[new_action_handle] = scenario._name
837
838         action_handles.append(new_action_handle)
839
840      # STEP 2: loop for the solver results, reading them and loading
841      #         them into instances as they are available.
842
843      if self._verbose is True:
844         print "Waiting for scenario sub-problem solves"
845
846      num_results_so_far = 0
847
848      while (num_results_so_far < len(self._scenario_tree._scenarios)):
849
850         action_handle = self._solver_manager.wait_any()
851         results = self._solver_manager.get_results(action_handle)
852         scenario_name = action_handle_scenario_map[action_handle]
853         instance = self._instances[scenario_name]
854
855         if self._verbose is True:
856            print "Results obtained for scenario="+scenario_name
857
858         if len(results.solution) == 0:
859            results.write(num=1)
860            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
861
862         if self._output_solver_results is True:
863            print "Results for scenario=",scenario_name
864            results.write(num=1)
865
866         start_time = time.time()
867         instance.load(results)
868         end_time = time.time()
869         if self._output_times is True:
870            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
871
872         if self._verbose is True:
873            print "Successfully loaded solution for scenario="+scenario_name
874
875         num_results_so_far = num_results_so_far + 1
876
877      if self._verbose is True:
878         print "Scenario sub-problem solves completed"
879
880      solve_end_time = time.time()
881      self._cumulative_solve_time += (solve_end_time - solve_start_time)
882
883      if self._output_times is True:
884         print "Aggregate sub-problem solve time this iteration=%8.2f" % (solve_end_time - solve_start_time)
885
886      if self._verbose is True:
887         print "Successfully completed PH iteration 0 solves - solution statistics:"
888         print "         Scenario              Objective                  Value"
889         for scenario in self._scenario_tree._scenarios:
890            instance = self._instances[scenario._name]
891            for objective_name in instance.active_components(Objective):
892               objective = instance.active_components(Objective)[objective_name]
893               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
894         print "------------------------------------------------"
895
896   #
897   # recompute the averages, minimum, and maximum statistics for all variables to be blended by PH, i.e.,
898   # not appearing in the final stage. technically speaking, the min/max aren't required by PH, but they
899   # are used often enough to warrant their computation and it's basically free if you're computing the
900   # average.
901   #
902   def update_variable_statistics(self):
903
904      start_time = time.time()
905
906      # compute statistics over all stages, even the last. this is necessary in order to
907      # successfully snapshot a scenario tree solution from the average values.
908      for stage in self._scenario_tree._stages:
909
910         for tree_node in stage._tree_nodes:
911
912            for (variable, index_template, variable_indices) in stage._variables:
913
914               variable_name = variable.name
915
916               avg_parameter_name = "PHAVG_"+variable_name
917
918               for index in variable_indices:
919                  min = float("inf")
920                  avg = 0.0
921                  max = float("-inf")
922                  node_probability = 0.0
923
924                  is_used = True # until proven otherwise
925                  for scenario in tree_node._scenarios:
926
927                     instance = self._instances[scenario._name]
928
929                     if getattr(instance,variable_name)[index].status == VarStatus.unused:
930                        is_used = False
931                     else:
932                        node_probability += scenario._probability
933                        var_value = getattr(instance, variable.name)[index].value
934                        if var_value < min:
935                           min = var_value
936                        avg += (scenario._probability * var_value)
937                        if var_value > max:
938                           max = var_value
939
940                  if is_used is True:
941
942                     tree_node._minimums[variable.name][index] = min
943                     tree_node._averages[variable.name][index] = avg / node_probability
944                     tree_node._maximums[variable.name][index] = max
945
946                     # distribute the newly computed average to the xbar variable in
947                     # each instance/scenario associated with this node. only do this
948                     # if the variable is used!
949                     for scenario in tree_node._scenarios:
950                        instance = self._instances[scenario._name]
951                        try:
952                           avg_parameter = getattr(instance, avg_parameter_name)
953                           avg_parameter[index] = avg / node_probability
954                        except:
955                           pass
956
957      end_time = time.time()
958      self._cumulative_xbar_time += (end_time - start_time)
959
960   def update_weights(self):
961
962      # because the weight updates rely on the xbars, and the xbars are node-based,
963      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
964      start_time = time.time()
965
966      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.
967
968         for tree_node in stage._tree_nodes:
969
970            for (variable, index_template, variable_indices) in stage._variables:
971
972               variable_name = variable.name
973               blend_parameter_name = "PHBLEND_"+variable_name
974               weight_parameter_name = "PHWEIGHT_"+variable_name
975               rho_parameter_name = "PHRHO_"+variable_name
976
977               for index in variable_indices:
978
979                  tree_node_average = value(tree_node._averages[variable.name][index])
980
981                  for scenario in tree_node._scenarios:
982
983                     instance = self._instances[scenario._name]
984
985                     if getattr(instance,variable.name)[index].status != VarStatus.unused:
986
987                        # we are currently not updating weights if blending is disabled for a variable.
988                        # this is done on the premise that unless you are actively trying to move
989                        # the variable toward the mean, the weights will blow up and be huge by the
990                        # time that blending is activated.
991                        variable_blend_indicator = value(getattr(instance, blend_parameter_name)[index])
992
993                        # get the weight and rho parameters for this variable/index combination.
994                        rho_value = value(getattr(instance, rho_parameter_name)[index])
995                        current_variable_weight = value(getattr(instance, weight_parameter_name)[index])
996
997                        # if I'm maximizing, invert value prior to adding (hack to implement negatives).
998                        # probably fixed in Pyomo at this point - I just haven't checked in a long while.
999                        if self._is_minimizing is False:
1000                           current_variable_weight = (-current_variable_weight)
1001                        current_variable_value = value(getattr(instance,variable.name)[index])
1002                        new_variable_weight = current_variable_weight + variable_blend_indicator * rho_value * (current_variable_value - tree_node_average)
1003                        # I have the correct updated value, so now invert if maximizing.
1004                        if self._is_minimizing is False:
1005                           new_variable_weight = (-new_variable_weight)
1006                        getattr(instance, weight_parameter_name)[index] = new_variable_weight
1007
1008      # we shouldn't have to re-simplify the expression, as we aren't adding any constant-variable terms - just modifying parameters.
1009
1010      end_time = time.time()
1011      self._cumulative_weight_time += (end_time - start_time)
1012
1013   def form_iteration_k_objectives(self):
1014
1015      for instance_name, instance in self._instances.items():
1016
1017         new_attrs = form_ph_objective(instance_name, \
1018                                       instance, \
1019                                       self._original_objective_expression[instance_name], \
1020                                       self._scenario_tree, \
1021                                       self._linearize_nonbinary_penalty_terms, \
1022                                       self._drop_proximal_terms, \
1023                                       self._retain_quadratic_binary_terms, \
1024                                       self._breakpoint_strategy, \
1025                                       self._integer_tolerance)
1026         self._instance_augmented_attributes[instance_name].extend(new_attrs)
1027
1028   def iteration_k_solve(self):
1029
1030      if self._verbose is True:
1031         print "------------------------------------------------"
1032         print "Starting PH iteration " + str(self._current_iteration) + " solves"
1033
1034      # cache the objective values generated by PH for output at the end of this function.
1035      ph_objective_values = {}
1036
1037      solve_start_time = time.time()
1038
1039      # STEP -1: if using a PH solver manager, propagate current weights/averages to the appropriate solver servers.
1040      #          ditto the tree node statistics, which are necessary if linearizing (so an optimization could be
1041      #          performed here).
1042      # NOTE: We aren't currently propagating rhos, as they generally don't change - we need to
1043      #       have a flag, though, indicating whether the rhos have changed, so they can be
1044      #       transmitted if needed.
1045      if self._solver_manager_type == "ph":
1046         self._transmit_weights_and_averages()
1047         self._transmit_tree_node_statistics()
1048
1049      # STEP 0: set up all global solver options.
1050      self._solver.mipgap = self._mipgap
1051
1052      # STEP 1: queue up the solves for all scenario sub-problems and
1053      #         grab all of the action handles for the subsequent barrier sync.
1054
1055      action_handles = []
1056      scenario_action_handle_map = {} # maps scenario names to action handles
1057      action_handle_scenario_map = {} # maps action handles to scenario names
1058
1059      for scenario in self._scenario_tree._scenarios:
1060
1061         instance = self._instances[scenario._name]
1062
1063         if self._verbose is True:
1064            print "Queuing solve for scenario=" + scenario._name
1065
1066         # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
1067         # don't do this, you won't see any chance in the output files as you vary the problem parameters!
1068         # ditto for instance fixing!
1069         instance.preprocess()
1070
1071         # once past iteration 0, there is always a feasible solution from which to warm-start.
1072         # however, you might want to disable warm-start when the solver is behaving badly (which does happen).
1073         new_action_handle = None
1074         if (self._disable_warmstarts is False) and (self._solver.warm_start_capable() is True):
1075            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
1076         else:
1077            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
1078
1079         scenario_action_handle_map[scenario._name] = new_action_handle
1080         action_handle_scenario_map[new_action_handle] = scenario._name
1081
1082         action_handles.append(new_action_handle)
1083
1084      # STEP 2: loop for the solver results, reading them and loading
1085      #         them into instances as they are available.
1086      if self._verbose is True:
1087         print "Waiting for scenario sub-problem solves"
1088
1089      num_results_so_far = 0
1090
1091      while (num_results_so_far < len(self._scenario_tree._scenarios)):
1092
1093         action_handle = self._solver_manager.wait_any()
1094         results = self._solver_manager.get_results(action_handle)
1095         scenario_name = action_handle_scenario_map[action_handle]
1096         instance = self._instances[scenario_name]
1097
1098         if self._verbose is True:
1099            print "Results obtained for scenario="+scenario_name
1100
1101         if len(results.solution) == 0:
1102            results.write(num=1)
1103            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
1104
1105         if self._output_solver_results is True:
1106            print "Results for scenario=",scenario_name
1107            results.write(num=1)
1108
1109         start_time = time.time()
1110         instance.load(results)
1111         end_time = time.time()
1112         if self._output_times is True:
1113            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
1114
1115         if self._verbose is True:
1116            print "Successfully loaded solution for scenario="+scenario_name
1117
1118         # we're assuming there is a single solution.
1119         # the "value" attribute is a pre-defined feature of any solution - it is relative to whatever
1120         # objective was selected during optimization, which of course should be the PH objective.
1121         try:
1122            ph_objective_values[instance.name] = float(results.solution(0).objective['f'].value)
1123         except AttributeError:
1124            # some solvers (e.g., through the SOL interface) don't report objective function values.
1125            ph_objective_values[instance.name] = 0.0
1126
1127         num_results_so_far = num_results_so_far + 1
1128
1129      if self._verbose is True:
1130         print "Scenario sub-problem solves completed"
1131
1132      solve_end_time = time.time()
1133      self._cumulative_solve_time += (solve_end_time - solve_start_time)
1134
1135      if self._output_times is True:
1136         print "Aggregate sub-problem solve time this iteration=%8.2f" % (solve_end_time - solve_start_time)
1137
1138      if self._verbose is True:
1139         print "Successfully completed PH iteration " + str(self._current_iteration) + " solves - solution statistics:"
1140         print "  Scenario             PH Objective             Cost Objective"
1141         for scenario in self._scenario_tree._scenarios:
1142            instance = self._instances[scenario._name]
1143            for objective_name in instance.active_components(Objective):
1144               objective = instance.active_components(Objective)[objective_name]
1145               print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], self._scenario_tree.compute_scenario_cost(instance))
1146
1147   def solve(self):
1148
1149      self._solve_start_time = time.time()
1150      self._cumulative_solve_time = 0.0
1151      self._cumulative_xbar_time = 0.0
1152      self._cumulative_weight_time = 0.0
1153
1154      print "Starting PH"
1155
1156      if self._initialized == False:
1157         raise RuntimeError, "PH is not initialized - cannot invoke solve() method"
1158
1159      # garbage collection noticeably slows down PH when dealing with
1160      # large numbers of scenarios. fortunately, there are well-defined
1161      # points at which garbage collection makes sense (and there isn't a
1162      # lot of collection to do). namely, after each PH iteration.
1163      re_enable_gc = gc.isenabled()
1164      gc.disable()
1165
1166      print "Initiating PH iteration=" + `self._current_iteration`
1167
1168      self.iteration_0_solve()
1169
1170      # update variable statistics prior to any output.
1171      self.update_variable_statistics()
1172
1173      if (self._verbose is True) or (self._report_solutions is True):
1174         print "Variable values following scenario solves:"
1175         self.pprint(False, False, True, False, output_only_statistics=self._report_only_statistics)
1176
1177      # let plugins know if they care.
1178      for plugin in self._ph_plugins:
1179         plugin.post_iteration_0_solves(self)
1180
1181      # update the fixed variable statistics.
1182      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
1183
1184      if self._verbose is True:
1185         print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1186         print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1187
1188      # always output the convergence metric and first-stage cost statistics, to give a sense of progress.
1189      self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1190      first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()
1191      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)
1192
1193      self.update_weights()
1194
1195      # let plugins know if they care.
1196      for plugin in self._ph_plugins:
1197         plugin.post_iteration_0(self)
1198
1199      # if using a PH solver server, trasnsmit the rhos prior to the iteration
1200      # k solve sequence. for now, we are assuming that the rhos don't change
1201      # on a per-iteration basis, but that assumption can be easily relaxed.
1202      # it is important to do this after the plugins have a chance to do their
1203      # computation.
1204      if self._solver_manager_type == "ph":
1205         self._transmit_rhos()
1206         self._enable_ph_objectives()
1207
1208      # checkpoint if it's time - which it always is after iteration 0,
1209      # if the interval is >= 1!
1210      if (self._checkpoint_interval > 0):
1211         self.checkpoint(0)
1212
1213      # garbage-collect if it wasn't disabled entirely.
1214      if re_enable_gc is True:
1215         gc.collect()
1216
1217      # there is an upper bound on the number of iterations to execute -
1218      # the actual bound depends on the converger supplied by the user.
1219      for i in range(1, self._max_iterations+1):
1220
1221         self._current_iteration = self._current_iteration + 1
1222
1223         print "Initiating PH iteration=" + `self._current_iteration`
1224
1225         if (self._verbose is True) or (self._report_weights is True):
1226            print "Variable averages and weights prior to scenario solves:"
1227            self.pprint(True, True, False, False, output_only_statistics=self._report_only_statistics)
1228
1229         # with the introduction of piecewise linearization, the form of the
1230         # penalty-weighted objective is no longer fixed. thus, we need to
1231         # create the objectives each PH iteration.
1232         self.form_iteration_k_objectives()
1233
1234         # let plugins know if they care.
1235         for plugin in self._ph_plugins:
1236            plugin.pre_iteration_k_solves(self)
1237
1238         # do the actual solves.
1239         self.iteration_k_solve()
1240
1241         # update variable statistics prior to any output.
1242         self.update_variable_statistics()
1243
1244         if (self._verbose is True) or (self._report_solutions is True):
1245            print "Variable values following scenario solves:"
1246            self.pprint(False, False, True, False, output_only_statistics=self._report_only_statistics)
1247
1248         # we don't technically have to do this at the last iteration,
1249         # but with checkpointing and re-starts, you're never sure
1250         # when you're executing the last iteration.
1251         self.update_weights()
1252
1253         # let plugins know if they care.
1254         for plugin in self._ph_plugins:
1255            plugin.post_iteration_k_solves(self)
1256
1257         # update the fixed variable statistics.
1258         (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
1259
1260         if self._verbose is True:
1261            print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1262            print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1263
1264         # let plugins know if they care.
1265         for plugin in self._ph_plugins:
1266            plugin.post_iteration_k(self)
1267
1268         # at this point, all the real work of an iteration is complete.
1269
1270         # checkpoint if it's time.
1271         if (self._checkpoint_interval > 0) and (i % self._checkpoint_interval is 0):
1272            self.checkpoint(i)
1273
1274         # check for early termination.
1275         self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1276         first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()
1277         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)
1278
1279         if self._converger.isConverged(self) is True:
1280            if self._total_discrete_vars == 0:
1281               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)
1282            else:
1283               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)+" or all discrete variables are fixed"
1284            break
1285
1286         # if we're terminating due to exceeding the maximum iteration count, print a message
1287         # indicating so - otherwise, you get a quiet, information-free output trace.
1288         if i == self._max_iterations:
1289            print "Halting PH - reached maximal iteration count="+str(self._max_iterations)
1290
1291         # garbage-collect if it wasn't disabled entirely.
1292         if re_enable_gc is True:
1293            gc.collect()
1294
1295      # re-enable the normal garbage collection mode.
1296      if re_enable_gc is True:
1297         gc.enable()
1298
1299      if self._verbose is True:
1300         print "Number of discrete variables fixed before final plugin calls="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1301         print "Number of continuous variables fixed before final plugin calls="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1302
1303      # let plugins know if they care. do this before
1304      # the final solution / statistics output, as the plugins
1305      # might do some final tweaking.
1306      for plugin in self._ph_plugins:
1307         plugin.post_ph_execution(self)
1308
1309      # update the fixed variable statistics - the plugins might have done something.
1310      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
1311
1312      self._solve_end_time = time.time()
1313
1314      print "PH complete"
1315
1316      print "Convergence history:"
1317      self._converger.pprint()
1318
1319      print "Final number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1320      print "Final number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1321
1322      # populate the scenario tree solution from the instances - to ensure consistent state
1323      # across the scenario tree instance and the scenario instances.
1324      self._scenario_tree.snapshotSolutionFromInstances(self._instances)     
1325
1326      print "Final variable values:"
1327      self.pprint(False, False, True, True, output_only_statistics=self._report_only_statistics)
1328
1329      print "Final costs:"
1330      self._scenario_tree.pprintCosts(self._instances)
1331
1332      if self._output_scenario_tree_solution is True:
1333         self._scenario_tree.snapshotSolutionFromAverages()
1334         print "Final solution (scenario tree format):"
1335         self._scenario_tree.pprintSolution()
1336
1337      if (self._verbose is True) and (self._output_times is True):
1338         print "Overall run-time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1339
1340      # cleanup the scenario instances for post-processing - ideally, we want to leave them in
1341      # their original state, minus all the PH-specific stuff. we don't do all cleanup (leaving
1342      # things like rhos, etc), but we do clean up constraints, as that really hoses up the ef writer.
1343      self._cleanup_scenario_instances()
1344
1345   #
1346   # prints a summary of all collected time statistics
1347   #
1348   def print_time_stats(self):
1349
1350      print "PH run-time statistics (user):"
1351
1352      print "Initialization time=  %8.2f seconds" % (self._init_end_time - self._init_start_time)
1353      print "Overall solve time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1354      print "Scenario solve time=  %8.2f seconds" % self._cumulative_solve_time
1355      print "Average update time=  %8.2f seconds" % self._cumulative_xbar_time
1356      print "Weight update time=   %8.2f seconds" % self._cumulative_weight_time
1357
1358   #
1359   # a utility to determine whether to output weight / average / etc. information for
1360   # a variable/node combination. when the printing is moved into a callback/plugin,
1361   # this routine will go there. for now, we don't dive down into the node resolution -
1362   # just the variable/stage.
1363   #
1364   def should_print(self, stage, variable, variable_indices):
1365
1366      if self._output_continuous_variable_stats is False:
1367
1368         variable_type = variable.domain
1369
1370         if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):
1371
1372            return False
1373
1374      return True
1375
1376   #
1377   # pretty-prints the state of the current variable averages, weights, and values.
1378   # inputs are booleans indicating which components should be output.
1379   #
1380   def pprint(self, output_averages, output_weights, output_values, output_fixed, output_only_statistics=False):
1381
1382      if self._initialized is False:
1383         raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"
1384
1385      # print tree nodes and associated variable/xbar/ph information in stage-order
1386      # we don't blend in the last stage, so we don't current care about printing the associated information.
1387      for stage in self._scenario_tree._stages[:-1]:
1388
1389         print "\tStage=" + str(stage._name)
1390
1391         num_outputs_this_stage = 0 # tracks the number of outputs on a per-index basis.
1392
1393         for (variable, index_template, variable_indices) in stage._variables:
1394
1395            variable_name = variable.name
1396
1397            if self.should_print(stage, variable, variable_indices) is True:
1398
1399               num_outputs_this_variable = 0 # track, so we don't output the variable names unless there is an entry to report.
1400
1401               for index in variable_indices:
1402
1403                  weight_parameter_name = "PHWEIGHT_"+variable_name
1404
1405                  num_outputs_this_index = 0 # track, so we don't output the variable index more than once.
1406
1407                  for tree_node in stage._tree_nodes:
1408
1409                     # determine if the variable/index pair is used across the set of scenarios (technically,
1410                     # it should be good enough to check one scenario). ditto for "fixed" status. fixed does
1411                     # imply unused (see note below), but we care about the fixed status when outputting
1412                     # final solutions.
1413
1414                     is_used = True # should be consistent across scenarios, so one "unused" flags as invalid.
1415                     is_fixed = False
1416
1417                     for scenario in tree_node._scenarios:
1418                        instance = self._instances[scenario._name]
1419                        variable_value = getattr(instance,variable_name)[index]
1420                        if variable_value.status == VarStatus.unused:
1421                           is_used = False
1422                        if variable_value.fixed is True:
1423                           is_fixed = True
1424
1425                     # IMPT: this is far from obvious, but variables that are fixed will - because
1426                     #       presolve will identify them as constants and eliminate them from all
1427                     #       expressions - be flagged as "unused" and therefore not output.
1428
1429                     if ((output_fixed is True) and (is_fixed is True)) or (is_used is True):
1430
1431                           minimum_value = None
1432                           maximum_value = None
1433
1434                           if index is None:
1435                              minimum_value = tree_node._minimums[variable_name]
1436                              maximum_value = tree_node._maximums[variable_name]
1437                           else:
1438                              minimum_value = value(tree_node._minimums[variable_name][index])
1439                              maximum_value = value(tree_node._maximums[variable_name][index])
1440
1441                           # there really isn't a need to output variables whose
1442                           # values are equal to 0 across-the-board. and there is
1443                           # good reason not to, i.e., the volume of output.
1444                           if (fabs(minimum_value) > self._integer_tolerance) or \
1445                              (fabs(maximum_value) > self._integer_tolerance):
1446
1447                              num_outputs_this_stage = num_outputs_this_stage + 1
1448                              num_outputs_this_variable = num_outputs_this_variable + 1
1449                              num_outputs_this_index = num_outputs_this_index + 1
1450
1451                              if num_outputs_this_variable == 1:
1452                                 print "\t\tVariable=" + variable_name
1453
1454                              if num_outputs_this_index == 1:
1455                                 if index is not None:
1456                                    print "\t\t\tIndex:", indexToString(index),
1457
1458                              if len(stage._tree_nodes) > 1:
1459                                 print ""
1460                                 print "\t\t\t\tTree Node="+tree_node._name,
1461                              if output_only_statistics is False:
1462                                 print "\t\t (Scenarios: ",
1463                                 for scenario in tree_node._scenarios:
1464                                    print scenario._name," ",
1465                                    if scenario == tree_node._scenarios[-1]:
1466                                       print ")"
1467
1468                              if output_values is True:
1469                                 average_value = value(tree_node._averages[variable_name][index])
1470                                 if output_only_statistics is False:
1471                                    print "\t\t\t\tValues: ",
1472                                 for scenario in tree_node._scenarios:
1473                                    instance = self._instances[scenario._name]
1474                                    this_value = getattr(instance,variable_name)[index].value
1475                                    if output_only_statistics is False:
1476                                       print "%12.4f" % this_value,
1477                                    if scenario == tree_node._scenarios[-1]:
1478                                       if output_only_statistics is True:
1479                                          # there technically isn't any good reason not to always report
1480                                          # the min and max; the only reason we're not doing this currently
1481                                          # is to avoid updating our regression test baseline output.
1482                                          print "    Min=%12.4f" % (minimum_value),
1483                                          print "    Avg=%12.4f" % (average_value),
1484                                          print "    Max=%12.4f" % (maximum_value),
1485                                       else:
1486                                          print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1487                                          print "    Avg=%12.4f" % (average_value),
1488                                       print ""
1489                              if output_weights:
1490                                 print "\t\t\t\tWeights: ",
1491                                 for scenario in tree_node._scenarios:
1492                                    instance = self._instances[scenario._name]
1493                                    print "%12.4f" % value(getattr(instance,weight_parameter_name)[index]),
1494                                    if scenario == tree_node._scenarios[-1]:
1495                                       print ""
1496
1497                              if output_averages:
1498                                 print "\t\t\t\tAverage: %12.4f" % (tree_node._averages[variable_name][index].value)
1499
1500         if num_outputs_this_stage == 0:
1501            print "\t\tNo non-converged variables in stage"
1502
1503         # cost variables aren't blended, so go through the gory computation of min/max/avg.
1504         # we currently always print these.
1505         cost_variable_name = stage._cost_variable[0].name
1506         cost_variable_index = stage._cost_variable[1]
1507         if cost_variable_index is None:
1508            print "\t\tCost Variable=" + cost_variable_name
1509         else:
1510            print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)
1511         for tree_node in stage._tree_nodes:
1512            print "\t\t\tTree Node=" + tree_node._name,
1513            if output_only_statistics is False:
1514               print "\t\t (Scenarios: ",
1515               for scenario in tree_node._scenarios:
1516                  print scenario._name," ",
1517                  if scenario == tree_node._scenarios[-1]:
1518                     print ")"
1519            maximum_value = 0.0
1520            minimum_value = 0.0
1521            sum_values = 0.0
1522            num_values = 0
1523            first_time = True
1524            if output_only_statistics is False:
1525               print "\t\t\tValues: ",
1526            else:
1527               print "\t\t\t",
1528            for scenario in tree_node._scenarios:
1529                instance = self._instances[scenario._name]
1530                this_value = getattr(instance,cost_variable_name)[cost_variable_index].value
1531                if output_only_statistics is False:
1532                   if this_value is not None:
1533                      print "%12.4f" % this_value,
1534                   else:
1535                      # this is a hack, in case the stage cost variables are not returned. ipopt
1536                      # does this occasionally, for example, if stage cost variables are constrained
1537                      # to a constant value (and consequently preprocessed out).
1538                      print "%12s" % "Not Rprted",
1539                if this_value is not None:
1540                   num_values += 1
1541                   sum_values += this_value
1542                   if first_time is True:
1543                      first_time = False
1544                      maximum_value = this_value
1545                      minimum_value = this_value
1546                   else:
1547                      if this_value > maximum_value:
1548                         maximum_value = this_value
1549                      if this_value < minimum_value:
1550                         minimum_value = this_value
1551                if scenario == tree_node._scenarios[-1]:
1552                   if num_values > 0:
1553                      if output_only_statistics is True:
1554                         print "    Min=%12.4f" % (minimum_value),
1555                         print "    Avg=%12.4f" % (sum_values/num_values),
1556                         print "    Max=%12.4f" % (maximum_value),
1557                      else:
1558                         print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1559                         print "    Avg=%12.4f" % (sum_values/num_values),
1560                   print ""
Note: See TracBrowser for help on using the repository browser.