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

Last change on this file since 3478 was 3478, checked in by khunter, 9 years ago

NFC: EOL whitespace removal

File size: 75.4 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         # important - mark the scenario instances as "concrete" - PH augments
663         # the instances in various ways, including via variable, constraint, and
664         # objective modifications when dealing with quadratic proximal and the weight
665         # penalty terms.
666         scenario_instance.concrete_mode()
667
668         # IMPT: disable canonical representation construction for ASL solvers.
669         #       this is a hack, in that we need to address encodings and
670         #       the like at a more general level.
671         if self._solver_type == "asl":
672            scenario_instance.skip_canonical_repn = True
673         else:
674            scenario_instance.preprocess()
675
676         self._instances[scenario._name] = scenario_instance
677         self._instances[scenario._name].name = scenario._name
678         self._instance_augmented_attributes[scenario._name] = []
679
680      # perform a single pass of garbage collection and re-enable automatic collection.
681      if re_enable_gc is True:
682         gc.collect()
683         gc.enable()
684
685      # with the scenario instances now available, have the scenario tree compute the
686      # variable match indices at each node.
687      self._scenario_tree.defineVariableIndexSets(self._instances)
688
689      # let plugins know if they care - this callback point allows
690      # users to create/modify the original scenario instances and/or
691      # the scenario tree prior to creating PH-related parameters,
692      # variables, and the like.
693      for plugin in self._ph_plugins:
694         plugin.post_instance_creation(self)
695
696      # create ph-specific parameters (weights, xbar, etc.) for each instance.
697
698      if self._verbose is True:
699         print "Creating weight, average, and rho parameter vectors for scenario instances"
700
701      self._create_ph_scenario_parameters()
702
703      # if specified, run the user script to initialize variable rhos at their whim.
704      if self._rho_setter is not None:
705         print "Executing user rho set script from filename="+self._rho_setter
706         execfile(self._rho_setter)
707
708      # with the instances created, run the user script to initialize variable bounds.
709      if self._bounds_setter is not None:
710         print "Executing user variable bounds set script from filename=", self._bounds_setter
711         execfile(self._bounds_setter)
712
713      # create parameters to store variable statistics (of general utility) at each node in the scenario tree.
714
715      if self._verbose is True:
716         print "Creating variable statistic (min/avg/max) parameter vectors for scenario tree nodes"
717
718      # do this for all stages, simply for completeness, i.e., to create a fully populated scenario tree.
719      for stage in self._scenario_tree._stages:
720
721         # first, gather all unique variables referenced in this stage
722         # this "gather" step is currently required because we're being lazy
723         # in terms of index management in the scenario tree - which
724         # should really be done in terms of sets of indices.
725
726         stage_variables = {}
727         for (reference_variable, index_template) in stage._variables:
728            if reference_variable.name not in stage_variables.keys():
729               stage_variables[reference_variable.name] = reference_variable
730
731         # next, create min/avg/max parameters for each variable in the corresponding tree node.
732         # NOTE: the parameter names below could really be empty, as they are never referenced
733         #       explicitly.
734         for (variable_name, reference_variable) in stage_variables.items():
735
736            for tree_node in stage._tree_nodes:
737
738               new_variable_index = getattr(self._instances[tree_node._scenarios[0]._name], variable_name)._index
739
740               new_min_index = new_variable_index
741               new_min_parameter_name = "NODEMIN_"+reference_variable.name
742               # this bit of ugliness is due to Pyomo not correctly handling the Param construction
743               # case when the supplied index set consists strictly of None, i.e., the source variable
744               # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
745               new_min_parameter = None
746               if (len(new_min_index) is 1) and (None in new_min_index):
747                  new_min_parameter = Param(name=new_min_parameter_name, nochecking=True)
748               else:
749                  new_min_parameter = Param(new_min_index, name=new_min_parameter_name, nochecking=True)
750               for index in new_min_index:
751                  new_min_parameter[index] = 0.0
752               tree_node._minimums[reference_variable.name] = new_min_parameter
753
754               new_avg_index = new_variable_index
755               new_avg_parameter_name = "NODEAVG_"+reference_variable.name
756               new_avg_parameter = None
757               if (len(new_avg_index) is 1) and (None in new_avg_index):
758                  new_avg_parameter = Param(name=new_avg_parameter_name, nochecking=True)
759               else:
760                  new_avg_parameter = Param(new_avg_index, name=new_avg_parameter_name, nochecking=True)
761               for index in new_avg_index:
762                  new_avg_parameter[index] = 0.0
763               tree_node._averages[reference_variable.name] = new_avg_parameter
764
765               new_max_index = new_variable_index
766               new_max_parameter_name = "NODEMAX_"+reference_variable.name
767               new_max_parameter = None
768               if (len(new_max_index) is 1) and (None in new_max_index):
769                  new_max_parameter = Param(name=new_max_parameter_name, nochecking=True)
770               else:
771                  new_max_parameter = Param(new_max_index, name=new_max_parameter_name, nochecking=True)
772               for index in new_max_index:
773                  new_max_parameter[index] = 0.0
774               tree_node._maximums[reference_variable.name] = new_max_parameter
775
776      # the objective functions are modified throughout the course of PH iterations.
777      # save the original, as a baseline to modify in subsequent iterations. reserve
778      # the original objectives, for subsequent modification.
779      self._original_objective_expression = {}
780      for instance_name, instance in self._instances.items():
781         objective_name = instance.active_components(Objective).keys()[0]
782         expr = instance.active_components(Objective)[objective_name]._data[None].expr
783         if isinstance(expr, Expression) is False:
784            expr = _IdentityExpression(expr)
785         self._original_objective_expression[instance_name] = expr
786
787      # cache the number of discrete and continuous variables in the master instance. this value
788      # is of general use, e.g., in the converger classes and in plugins.
789      (self._total_discrete_vars,self._total_continuous_vars) = self.compute_variable_counts()
790      if self._verbose is True:
791         print "Total number of discrete instance variables="+str(self._total_discrete_vars)
792         print "Total number of continuous instance variables="+str(self._total_continuous_vars)
793
794      # track the total number of fixed variables of each category at the end of each PH iteration.
795      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
796
797      # indicate that we're ready to run.
798      self._initialized = True
799
800      if self._verbose is True:
801         print "PH successfully created model instances for all scenarios"
802
803      self._init_end_time = time.time()
804
805      if self._verbose is True:
806         print "PH is successfully initialized"
807         if self._output_times is True:
808            print "Initialization time=%8.2f seconds" % (self._init_end_time - self._init_start_time)
809
810      # let plugins know if they care.
811      if self._verbose is True:
812         print "Initializing PH plugins"
813      for plugin in self._ph_plugins:
814         plugin.post_ph_initialization(self)
815      if self._verbose is True:
816         print "PH plugin initialization complete"
817
818   """ Perform the non-weighted scenario solves and form the initial w and xbars.
819   """
820   def iteration_0_solve(self):
821
822      if self._verbose is True:
823         print "------------------------------------------------"
824         print "Starting PH iteration 0 solves"
825
826      self._current_iteration = 0
827
828      solve_start_time = time.time()
829
830      # STEP 0: set up all global solver options.
831      self._solver.mipgap = self._mipgap
832
833      # STEP 1: queue up the solves for all scenario sub-problems and
834      #         grab all the action handles for the subsequent barrier sync.
835
836      action_handles = []
837      scenario_action_handle_map = {} # maps scenario names to action handles
838      action_handle_scenario_map = {} # maps action handles to scenario names
839
840      for scenario in self._scenario_tree._scenarios:
841
842         instance = self._instances[scenario._name]
843
844         if self._verbose is True:
845            print "Queuing solve for scenario=" + scenario._name
846
847         # IMPT: You have to re-presolve if approximating continuous variable penalty terms with a
848         #       piecewise linear function. otherwise, the newly introduced variables won't be flagged
849         #       as unused (as is correct for iteration 0), and the output writer will crater.
850         # IMPT: I decided to presolve unconditionally, as PH extensions can add arbitrary components
851         #       to the base scenario instances - and the variable values/etc. need to be collectged.
852         instance.preprocess()
853
854         # there's nothing to warm-start from in iteration 0, so don't include the keyword in the solve call.
855         # the reason you don't want to include it is that some solvers don't know how to handle the keyword
856         # at all (despite it being false). you might want to solve iteration 0 solves using some other solver.
857
858         new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
859         scenario_action_handle_map[scenario._name] = new_action_handle
860         action_handle_scenario_map[new_action_handle] = scenario._name
861
862         action_handles.append(new_action_handle)
863
864      # STEP 2: loop for the solver results, reading them and loading
865      #         them into instances as they are available.
866
867      if self._verbose is True:
868         print "Waiting for scenario sub-problem solves"
869
870      num_results_so_far = 0
871
872      while (num_results_so_far < len(self._scenario_tree._scenarios)):
873
874         action_handle = self._solver_manager.wait_any()
875         results = self._solver_manager.get_results(action_handle)
876         scenario_name = action_handle_scenario_map[action_handle]
877         instance = self._instances[scenario_name]
878
879         if self._verbose is True:
880            print "Results obtained for scenario="+scenario_name
881
882         if len(results.solution) == 0:
883            results.write(num=1)
884            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
885
886         if self._output_solver_results is True:
887            print "Results for scenario=",scenario_name
888            results.write(num=1)
889
890         start_time = time.time()
891         instance.load(results)
892         end_time = time.time()
893         if self._output_times is True:
894            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
895
896         if self._verbose is True:
897            print "Successfully loaded solution for scenario="+scenario_name
898
899         num_results_so_far = num_results_so_far + 1
900
901      if self._verbose is True:
902         print "Scenario sub-problem solves completed"
903
904      solve_end_time = time.time()
905      self._cumulative_solve_time += (solve_end_time - solve_start_time)
906
907      if self._output_times is True:
908         print "Aggregate sub-problem solve time this iteration=%8.2f" % (solve_end_time - solve_start_time)
909
910      if self._verbose is True:
911         print "Successfully completed PH iteration 0 solves - solution statistics:"
912         print "         Scenario              Objective                  Value"
913         for scenario in self._scenario_tree._scenarios:
914            instance = self._instances[scenario._name]
915            for objective_name in instance.active_components(Objective):
916               objective = instance.active_components(Objective)[objective_name]
917               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
918         print "------------------------------------------------"
919
920   #
921   # recompute the averages, minimum, and maximum statistics for all variables to be blended by PH, i.e.,
922   # not appearing in the final stage. technically speaking, the min/max aren't required by PH, but they
923   # are used often enough to warrant their computation and it's basically free if you're computing the
924   # average.
925   #
926   def update_variable_statistics(self):
927
928      start_time = time.time()
929
930      # NOTE: the following code has some optimizations that are not normally recommended, in particular
931      #       the direct access and manipulation of parameters via the .value attribute instead of the
932      #       user-level-preferred value() method. this is justifiable in this particular instance
933      #       because we are creating the PH parameters (and therefore can manipulate them safely), and
934      #       this routine takes a non-trivial amount of the overall run-time.
935
936      # compute statistics over all stages, even the last. this is necessary in order to
937      # successfully snapshot a scenario tree solution from the average values.
938      for stage in self._scenario_tree._stages:
939
940         for tree_node in stage._tree_nodes:
941
942            for (variable, index_template) in stage._variables:
943
944               variable_name = variable.name
945
946               variable_indices = tree_node._variable_indices[variable_name]
947
948               avg_parameter_name = "PHAVG_"+variable_name
949
950               tree_node_var_mins = tree_node._minimums[variable_name]
951               tree_node_var_avgs = tree_node._averages[variable_name]
952               tree_node_var_maxs = tree_node._maximums[variable_name]
953
954               scenario_variables = []
955               for scenario in tree_node._scenarios:
956                  instance = self._instances[scenario._name]
957                  scenario_variables.append(getattr(instance, variable_name))
958
959               for index in variable_indices:
960
961                  min = float("inf")
962                  avg = 0.0
963                  max = float("-inf")
964                  node_probability = 0.0
965
966                  is_used = True # until proven otherwise
967                  for scenario_variable in scenario_variables:
968
969                     if scenario_variable[index].status == VarStatus.unused:
970                        is_used = False
971                     else:
972                        node_probability += scenario._probability
973                        var_value = scenario_variable[index].value
974                        if var_value < min:
975                           min = var_value
976                        avg += (scenario._probability * var_value)
977                        if var_value > max:
978                           max = var_value
979
980                  if is_used is True:
981
982                     tree_node_var_mins[index].value = min
983                     tree_node_var_avgs[index].value = avg / node_probability
984                     tree_node_var_maxs[index].value = max
985
986                     # distribute the newly computed average to the xbar variable in
987                     # each instance/scenario associated with this node. only do this
988                     # if the variable is used!
989                     for scenario in tree_node._scenarios:
990                        instance = self._instances[scenario._name]
991                        try:
992                           avg_parameter = getattr(instance, avg_parameter_name)
993                           avg_parameter[index].value = avg / node_probability
994                        except:
995                           pass
996
997      end_time = time.time()
998      self._cumulative_xbar_time += (end_time - start_time)
999
1000   def update_weights(self):
1001
1002      # because the weight updates rely on the xbars, and the xbars are node-based,
1003      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
1004      start_time = time.time()
1005
1006      # NOTE: the following code has some optimizations that are not normally recommended, in particular
1007      #       the direct access and manipulation of parameters via the .value attribute instead of the
1008      #       user-level-preferred value() method. this is justifiable in this particular instance
1009      #       because we are creating the PH parameters (and therefore can manipulate them safely), and
1010      #       this routine takes a non-trivial amount of the overall run-time.
1011
1012      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.
1013
1014         for tree_node in stage._tree_nodes:
1015
1016            for (variable, index_template) in stage._variables:
1017
1018               variable_name = variable.name
1019               blend_parameter_name = "PHBLEND_"+variable_name
1020               weight_parameter_name = "PHWEIGHT_"+variable_name
1021               rho_parameter_name = "PHRHO_"+variable_name
1022
1023               variable_indices = tree_node._variable_indices[variable_name]
1024
1025               for index in variable_indices:
1026
1027                  tree_node_average = tree_node._averages[variable_name][index].value
1028
1029                  for scenario in tree_node._scenarios:
1030
1031                     instance = self._instances[scenario._name]
1032
1033                     if getattr(instance,variable_name)[index].status != VarStatus.unused:
1034
1035                        weight_parameter = getattr(instance, weight_parameter_name)
1036
1037                        # we are currently not updating weights if blending is disabled for a variable.
1038                        # this is done on the premise that unless you are actively trying to move
1039                        # the variable toward the mean, the weights will blow up and be huge by the
1040                        # time that blending is activated.
1041                        variable_blend_indicator = getattr(instance, blend_parameter_name)[index].value
1042
1043                        # get the weight and rho parameters for this variable/index combination.
1044                        rho_value = getattr(instance, rho_parameter_name)[index].value
1045                        current_variable_weight = weight_parameter[index].value
1046
1047                        # if I'm maximizing, invert value prior to adding (hack to implement negatives).
1048                        # probably fixed in Pyomo at this point - I just haven't checked in a long while.
1049                        if self._is_minimizing is False:
1050                           current_variable_weight = (-current_variable_weight)
1051                        current_variable_value = getattr(instance, variable_name)[index]()
1052                        new_variable_weight = current_variable_weight + variable_blend_indicator * rho_value * (current_variable_value - tree_node_average)
1053                        # I have the correct updated value, so now invert if maximizing.
1054                        if self._is_minimizing is False:
1055                           new_variable_weight = (-new_variable_weight)
1056                        weight_parameter[index].value = new_variable_weight
1057
1058      # we shouldn't have to re-simplify the expression, as we aren't adding any constant-variable terms - just modifying parameters.
1059
1060      end_time = time.time()
1061      self._cumulative_weight_time += (end_time - start_time)
1062
1063   def form_iteration_k_objectives(self):
1064
1065      for instance_name, instance in self._instances.items():
1066
1067         new_attrs = form_ph_objective(instance_name, \
1068                                       instance, \
1069                                       self._original_objective_expression[instance_name], \
1070                                       self._scenario_tree, \
1071                                       self._linearize_nonbinary_penalty_terms, \
1072                                       self._drop_proximal_terms, \
1073                                       self._retain_quadratic_binary_terms, \
1074                                       self._breakpoint_strategy, \
1075                                       self._integer_tolerance, \
1076                                       self._simplify_expressions)
1077         self._instance_augmented_attributes[instance_name].extend(new_attrs)
1078
1079   def iteration_k_solve(self):
1080
1081      if self._verbose is True:
1082         print "------------------------------------------------"
1083         print "Starting PH iteration " + str(self._current_iteration) + " solves"
1084
1085      # cache the objective values generated by PH for output at the end of this function.
1086      ph_objective_values = {}
1087
1088      solve_start_time = time.time()
1089
1090      # STEP -1: if using a PH solver manager, propagate current weights/averages to the appropriate solver servers.
1091      #          ditto the tree node statistics, which are necessary if linearizing (so an optimization could be
1092      #          performed here).
1093      # NOTE: We aren't currently propagating rhos, as they generally don't change - we need to
1094      #       have a flag, though, indicating whether the rhos have changed, so they can be
1095      #       transmitted if needed.
1096      if self._solver_manager_type == "ph":
1097         self._transmit_weights_and_averages()
1098         self._transmit_tree_node_statistics()
1099
1100      # STEP 0: set up all global solver options.
1101      self._solver.mipgap = self._mipgap
1102
1103      # STEP 1: queue up the solves for all scenario sub-problems and
1104      #         grab all of the action handles for the subsequent barrier sync.
1105
1106      action_handles = []
1107      scenario_action_handle_map = {} # maps scenario names to action handles
1108      action_handle_scenario_map = {} # maps action handles to scenario names
1109
1110      for scenario in self._scenario_tree._scenarios:
1111
1112         instance = self._instances[scenario._name]
1113
1114         if self._verbose is True:
1115            print "Queuing solve for scenario=" + scenario._name
1116
1117         # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
1118         # don't do this, you won't see any chance in the output files as you vary the problem parameters!
1119         # ditto for instance fixing!
1120         instance.preprocess()
1121
1122         # once past iteration 0, there is always a feasible solution from which to warm-start.
1123         # however, you might want to disable warm-start when the solver is behaving badly (which does happen).
1124         new_action_handle = None
1125         if (self._disable_warmstarts is False) and (self._solver.warm_start_capable() is True):
1126            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
1127         else:
1128            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
1129
1130         scenario_action_handle_map[scenario._name] = new_action_handle
1131         action_handle_scenario_map[new_action_handle] = scenario._name
1132
1133         action_handles.append(new_action_handle)
1134
1135      # STEP 2: loop for the solver results, reading them and loading
1136      #         them into instances as they are available.
1137      if self._verbose is True:
1138         print "Waiting for scenario sub-problem solves"
1139
1140      num_results_so_far = 0
1141
1142      while (num_results_so_far < len(self._scenario_tree._scenarios)):
1143
1144         action_handle = self._solver_manager.wait_any()
1145         results = self._solver_manager.get_results(action_handle)
1146         scenario_name = action_handle_scenario_map[action_handle]
1147         instance = self._instances[scenario_name]
1148
1149         if self._verbose is True:
1150            print "Results obtained for scenario="+scenario_name
1151
1152         if len(results.solution) == 0:
1153            results.write(num=1)
1154            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
1155
1156         if self._output_solver_results is True:
1157            print "Results for scenario=",scenario_name
1158            results.write(num=1)
1159
1160         start_time = time.time()
1161         instance.load(results)
1162         end_time = time.time()
1163         if self._output_times is True:
1164            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
1165
1166         if self._verbose is True:
1167            print "Successfully loaded solution for scenario="+scenario_name
1168
1169         # we're assuming there is a single solution.
1170         # the "value" attribute is a pre-defined feature of any solution - it is relative to whatever
1171         # objective was selected during optimization, which of course should be the PH objective.
1172         try:
1173            ph_objective_values[instance.name] = float(results.solution(0).objective['f'].value)
1174         except AttributeError:
1175            # some solvers (e.g., through the SOL interface) don't report objective function values.
1176            ph_objective_values[instance.name] = 0.0
1177
1178         num_results_so_far = num_results_so_far + 1
1179
1180      if self._verbose is True:
1181         print "Scenario sub-problem solves completed"
1182
1183      solve_end_time = time.time()
1184      self._cumulative_solve_time += (solve_end_time - solve_start_time)
1185
1186      if self._output_times is True:
1187         print "Aggregate sub-problem solve time this iteration=%8.2f" % (solve_end_time - solve_start_time)
1188
1189      if self._verbose is True:
1190         print "Successfully completed PH iteration " + str(self._current_iteration) + " solves - solution statistics:"
1191         print "  Scenario             PH Objective             Cost Objective"
1192         for scenario in self._scenario_tree._scenarios:
1193            instance = self._instances[scenario._name]
1194            for objective_name in instance.active_components(Objective):
1195               objective = instance.active_components(Objective)[objective_name]
1196               print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], self._scenario_tree.compute_scenario_cost(instance))
1197
1198   def solve(self):
1199
1200      self._solve_start_time = time.time()
1201      self._cumulative_solve_time = 0.0
1202      self._cumulative_xbar_time = 0.0
1203      self._cumulative_weight_time = 0.0
1204
1205      print "Starting PH"
1206
1207      if self._initialized == False:
1208         raise RuntimeError, "PH is not initialized - cannot invoke solve() method"
1209
1210      # garbage collection noticeably slows down PH when dealing with
1211      # large numbers of scenarios. fortunately, there are well-defined
1212      # points at which garbage collection makes sense (and there isn't a
1213      # lot of collection to do). namely, after each PH iteration.
1214      re_enable_gc = gc.isenabled()
1215      gc.disable()
1216
1217      print "Initiating PH iteration=" + `self._current_iteration`
1218
1219      self.iteration_0_solve()
1220
1221      # update variable statistics prior to any output.
1222      self.update_variable_statistics()
1223
1224      if (self._verbose is True) or (self._report_solutions is True):
1225         print "Variable values following scenario solves:"
1226         self.pprint(False, False, True, False, output_only_statistics=self._report_only_statistics)
1227
1228      # let plugins know if they care.
1229      for plugin in self._ph_plugins:
1230         plugin.post_iteration_0_solves(self)
1231
1232      # update the fixed variable statistics.
1233      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
1234
1235      if self._verbose is True:
1236         print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1237         print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1238
1239      # always output the convergence metric and first-stage cost statistics, to give a sense of progress.
1240      self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1241      first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()
1242      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)
1243
1244      self.update_weights()
1245
1246      # let plugins know if they care.
1247      for plugin in self._ph_plugins:
1248         plugin.post_iteration_0(self)
1249
1250      # if using a PH solver server, trasnsmit the rhos prior to the iteration
1251      # k solve sequence. for now, we are assuming that the rhos don't change
1252      # on a per-iteration basis, but that assumption can be easily relaxed.
1253      # it is important to do this after the plugins have a chance to do their
1254      # computation.
1255      if self._solver_manager_type == "ph":
1256         self._transmit_rhos()
1257         self._enable_ph_objectives()
1258
1259      # checkpoint if it's time - which it always is after iteration 0,
1260      # if the interval is >= 1!
1261      if (self._checkpoint_interval > 0):
1262         self.checkpoint(0)
1263
1264      # garbage-collect if it wasn't disabled entirely.
1265      if re_enable_gc is True:
1266         gc.collect()
1267
1268      # there is an upper bound on the number of iterations to execute -
1269      # the actual bound depends on the converger supplied by the user.
1270      for i in range(1, self._max_iterations+1):
1271
1272         self._current_iteration = self._current_iteration + 1
1273
1274         print "Initiating PH iteration=" + `self._current_iteration`
1275
1276         if (self._verbose is True) or (self._report_weights is True):
1277            print "Variable averages and weights prior to scenario solves:"
1278            self.pprint(True, True, False, False, output_only_statistics=self._report_only_statistics)
1279
1280         # with the introduction of piecewise linearization, the form of the
1281         # penalty-weighted objective is no longer fixed. thus, we need to
1282         # create the objectives each PH iteration. we could be fancy in
1283         # the case of quadratic objectives (eliminating redundancy), but
1284         # if we decide to do this, see the "IMPT" note below.
1285         self.form_iteration_k_objectives()
1286
1287         # IMPT: for efficiency purposes, the PH objective formation calls above create expression trees by
1288         #       evaluating the parameters involved in either quadratic or linear expressions directly, e.g.,
1289         #       via value(). as a consequence, if the user modifies PH weights or rhos in the plugin
1290         #       callback below, it will *not* be propagated to the objective expression.
1291
1292         # let plugins know if they care.
1293         for plugin in self._ph_plugins:
1294            plugin.pre_iteration_k_solves(self)
1295
1296         # do the actual solves.
1297         self.iteration_k_solve()
1298
1299         # update variable statistics prior to any output.
1300         self.update_variable_statistics()
1301
1302         if (self._verbose is True) or (self._report_solutions is True):
1303            print "Variable values following scenario solves:"
1304            self.pprint(False, False, True, False, output_only_statistics=self._report_only_statistics)
1305
1306         # we don't technically have to do this at the last iteration,
1307         # but with checkpointing and re-starts, you're never sure
1308         # when you're executing the last iteration.
1309         self.update_weights()
1310
1311         # let plugins know if they care.
1312         for plugin in self._ph_plugins:
1313            plugin.post_iteration_k_solves(self)
1314
1315         # update the fixed variable statistics.
1316         (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
1317
1318         if self._verbose is True:
1319            print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1320            print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1321
1322         # let plugins know if they care.
1323         for plugin in self._ph_plugins:
1324            plugin.post_iteration_k(self)
1325
1326         # at this point, all the real work of an iteration is complete.
1327
1328         # checkpoint if it's time.
1329         if (self._checkpoint_interval > 0) and (i % self._checkpoint_interval is 0):
1330            self.checkpoint(i)
1331
1332         # check for early termination.
1333         self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1334         first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()
1335         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)
1336
1337         if self._converger.isConverged(self) is True:
1338            if self._total_discrete_vars == 0:
1339               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)
1340            else:
1341               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)+" or all discrete variables are fixed"
1342            break
1343
1344         # if we're terminating due to exceeding the maximum iteration count, print a message
1345         # indicating so - otherwise, you get a quiet, information-free output trace.
1346         if i == self._max_iterations:
1347            print "Halting PH - reached maximal iteration count="+str(self._max_iterations)
1348
1349         # garbage-collect if it wasn't disabled entirely.
1350         if re_enable_gc is True:
1351            gc.collect()
1352
1353      # re-enable the normal garbage collection mode.
1354      if re_enable_gc is True:
1355         gc.enable()
1356
1357      if self._verbose is True:
1358         print "Number of discrete variables fixed before final plugin calls="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1359         print "Number of continuous variables fixed before final plugin calls="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1360
1361      # let plugins know if they care. do this before
1362      # the final solution / statistics output, as the plugins
1363      # might do some final tweaking.
1364      for plugin in self._ph_plugins:
1365         plugin.post_ph_execution(self)
1366
1367      # update the fixed variable statistics - the plugins might have done something.
1368      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
1369
1370      self._solve_end_time = time.time()
1371
1372      print "PH complete"
1373
1374      print "Convergence history:"
1375      self._converger.pprint()
1376
1377      print "Final number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1378      print "Final number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1379
1380      # populate the scenario tree solution from the instances - to ensure consistent state
1381      # across the scenario tree instance and the scenario instances.
1382      self._scenario_tree.snapshotSolutionFromInstances(self._instances)
1383
1384      print "Final variable values:"
1385      self.pprint(False, False, True, True, output_only_statistics=self._report_only_statistics)
1386
1387      print "Final costs:"
1388      self._scenario_tree.pprintCosts(self._instances)
1389
1390      if self._output_scenario_tree_solution is True:
1391         self._scenario_tree.snapshotSolutionFromAverages(self._instances)
1392         print "Final solution (scenario tree format):"
1393         self._scenario_tree.pprintSolution()
1394
1395      if (self._verbose is True) and (self._output_times is True):
1396         print "Overall run-time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1397
1398      # cleanup the scenario instances for post-processing - ideally, we want to leave them in
1399      # their original state, minus all the PH-specific stuff. we don't do all cleanup (leaving
1400      # things like rhos, etc), but we do clean up constraints, as that really hoses up the ef writer.
1401      self._cleanup_scenario_instances()
1402
1403   #
1404   # prints a summary of all collected time statistics
1405   #
1406   def print_time_stats(self):
1407
1408      print "PH run-time statistics (user):"
1409
1410      print "Initialization time=  %8.2f seconds" % (self._init_end_time - self._init_start_time)
1411      print "Overall solve time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1412      print "Scenario solve time=  %8.2f seconds" % self._cumulative_solve_time
1413      print "Average update time=  %8.2f seconds" % self._cumulative_xbar_time
1414      print "Weight update time=   %8.2f seconds" % self._cumulative_weight_time
1415
1416   #
1417   # a utility to determine whether to output weight / average / etc. information for
1418   # a variable/node combination. when the printing is moved into a callback/plugin,
1419   # this routine will go there. for now, we don't dive down into the node resolution -
1420   # just the variable/stage.
1421   #
1422   def should_print(self, stage, variable):
1423
1424      if self._output_continuous_variable_stats is False:
1425
1426         variable_type = variable.domain
1427
1428         if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):
1429
1430            return False
1431
1432      return True
1433
1434   #
1435   # pretty-prints the state of the current variable averages, weights, and values.
1436   # inputs are booleans indicating which components should be output.
1437   #
1438   def pprint(self, output_averages, output_weights, output_values, output_fixed, output_only_statistics=False):
1439
1440      if self._initialized is False:
1441         raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"
1442
1443      # print tree nodes and associated variable/xbar/ph information in stage-order
1444      # we don't blend in the last stage, so we don't current care about printing the associated information.
1445      for stage in self._scenario_tree._stages[:-1]:
1446
1447         print "\tStage=" + str(stage._name)
1448
1449         num_outputs_this_stage = 0 # tracks the number of outputs on a per-index basis.
1450
1451         for (variable, index_template) in stage._variables:
1452
1453            variable_name = variable.name
1454
1455            if self.should_print(stage, variable) is True:
1456
1457               num_outputs_this_variable = 0 # track, so we don't output the variable names unless there is an entry to report.
1458
1459               for tree_node in stage._tree_nodes:
1460
1461                  variable_indices = tree_node._variable_indices[variable_name]
1462
1463                  for index in variable_indices:
1464
1465                     weight_parameter_name = "PHWEIGHT_"+variable_name
1466
1467                     num_outputs_this_index = 0 # track, so we don't output the variable index more than once.
1468
1469                     # determine if the variable/index pair is used across the set of scenarios (technically,
1470                     # it should be good enough to check one scenario). ditto for "fixed" status. fixed does
1471                     # imply unused (see note below), but we care about the fixed status when outputting
1472                     # final solutions.
1473
1474                     is_used = True # should be consistent across scenarios, so one "unused" flags as invalid.
1475                     is_fixed = False
1476
1477                     for scenario in tree_node._scenarios:
1478                        instance = self._instances[scenario._name]
1479                        variable_value = getattr(instance,variable_name)[index]
1480                        if variable_value.status == VarStatus.unused:
1481                           is_used = False
1482                        if variable_value.fixed is True:
1483                           is_fixed = True
1484
1485                     # IMPT: this is far from obvious, but variables that are fixed will - because
1486                     #       presolve will identify them as constants and eliminate them from all
1487                     #       expressions - be flagged as "unused" and therefore not output.
1488
1489                     if ((output_fixed is True) and (is_fixed is True)) or (is_used is True):
1490
1491                           minimum_value = None
1492                           maximum_value = None
1493
1494                           if index is None:
1495                              minimum_value = tree_node._minimums[variable_name]
1496                              maximum_value = tree_node._maximums[variable_name]
1497                           else:
1498                              minimum_value = value(tree_node._minimums[variable_name][index])
1499                              maximum_value = value(tree_node._maximums[variable_name][index])
1500
1501                           # there really isn't a need to output variables whose
1502                           # values are equal to 0 across-the-board. and there is
1503                           # good reason not to, i.e., the volume of output.
1504                           if (fabs(minimum_value) > self._integer_tolerance) or \
1505                              (fabs(maximum_value) > self._integer_tolerance):
1506
1507                              num_outputs_this_stage = num_outputs_this_stage + 1
1508                              num_outputs_this_variable = num_outputs_this_variable + 1
1509                              num_outputs_this_index = num_outputs_this_index + 1
1510
1511                              if num_outputs_this_variable == 1:
1512                                 print "\t\tVariable=" + variable_name
1513
1514                              if num_outputs_this_index == 1:
1515                                 if index is not None:
1516                                    print "\t\t\tIndex:", indexToString(index),
1517
1518                              if len(stage._tree_nodes) > 1:
1519                                 print ""
1520                                 print "\t\t\t\tTree Node="+tree_node._name,
1521                              if output_only_statistics is False:
1522                                 print "\t\t (Scenarios: ",
1523                                 for scenario in tree_node._scenarios:
1524                                    print scenario._name," ",
1525                                    if scenario == tree_node._scenarios[-1]:
1526                                       print ")"
1527
1528                              if output_values is True:
1529                                 average_value = value(tree_node._averages[variable_name][index])
1530                                 if output_only_statistics is False:
1531                                    print "\t\t\t\tValues: ",
1532                                 for scenario in tree_node._scenarios:
1533                                    instance = self._instances[scenario._name]
1534                                    this_value = getattr(instance,variable_name)[index].value
1535                                    if output_only_statistics is False:
1536                                       print "%12.4f" % this_value,
1537                                    if scenario == tree_node._scenarios[-1]:
1538                                       if output_only_statistics is True:
1539                                          # there technically isn't any good reason not to always report
1540                                          # the min and max; the only reason we're not doing this currently
1541                                          # is to avoid updating our regression test baseline output.
1542                                          print "    Min=%12.4f" % (minimum_value),
1543                                          print "    Avg=%12.4f" % (average_value),
1544                                          print "    Max=%12.4f" % (maximum_value),
1545                                       else:
1546                                          print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1547                                          print "    Avg=%12.4f" % (average_value),
1548                                       print ""
1549                              if output_weights:
1550                                 print "\t\t\t\tWeights: ",
1551                                 for scenario in tree_node._scenarios:
1552                                    instance = self._instances[scenario._name]
1553                                    print "%12.4f" % value(getattr(instance,weight_parameter_name)[index]),
1554                                    if scenario == tree_node._scenarios[-1]:
1555                                       print ""
1556
1557                              if output_averages:
1558                                 print "\t\t\t\tAverage: %12.4f" % (tree_node._averages[variable_name][index].value)
1559
1560         if num_outputs_this_stage == 0:
1561            print "\t\tNo non-converged variables in stage"
1562
1563         # cost variables aren't blended, so go through the gory computation of min/max/avg.
1564         # we currently always print these.
1565         cost_variable_name = stage._cost_variable[0].name
1566         cost_variable_index = stage._cost_variable[1]
1567         if cost_variable_index is None:
1568            print "\t\tCost Variable=" + cost_variable_name
1569         else:
1570            print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)
1571         for tree_node in stage._tree_nodes:
1572            print "\t\t\tTree Node=" + tree_node._name,
1573            if output_only_statistics is False:
1574               print "\t\t (Scenarios: ",
1575               for scenario in tree_node._scenarios:
1576                  print scenario._name," ",
1577                  if scenario == tree_node._scenarios[-1]:
1578                     print ")"
1579            maximum_value = 0.0
1580            minimum_value = 0.0
1581            sum_values = 0.0
1582            num_values = 0
1583            first_time = True
1584            if output_only_statistics is False:
1585               print "\t\t\tValues: ",
1586            else:
1587               print "\t\t\t",
1588            for scenario in tree_node._scenarios:
1589                instance = self._instances[scenario._name]
1590                this_value = getattr(instance,cost_variable_name)[cost_variable_index].value
1591                if output_only_statistics is False:
1592                   if this_value is not None:
1593                      print "%12.4f" % this_value,
1594                   else:
1595                      # this is a hack, in case the stage cost variables are not returned. ipopt
1596                      # does this occasionally, for example, if stage cost variables are constrained
1597                      # to a constant value (and consequently preprocessed out).
1598                      print "%12s" % "Not Rprted",
1599                if this_value is not None:
1600                   num_values += 1
1601                   sum_values += this_value
1602                   if first_time is True:
1603                      first_time = False
1604                      maximum_value = this_value
1605                      minimum_value = this_value
1606                   else:
1607                      if this_value > maximum_value:
1608                         maximum_value = this_value
1609                      if this_value < minimum_value:
1610                         minimum_value = this_value
1611                if scenario == tree_node._scenarios[-1]:
1612                   if num_values > 0:
1613                      if output_only_statistics is True:
1614                         print "    Min=%12.4f" % (minimum_value),
1615                         print "    Avg=%12.4f" % (sum_values/num_values),
1616                         print "    Max=%12.4f" % (maximum_value),
1617                      else:
1618                         print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1619                         print "    Avg=%12.4f" % (sum_values/num_values),
1620                   print ""
Note: See TracBrowser for help on using the repository browser.