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

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

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

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