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

Last change on this file since 3152 was 3152, checked in by jwatson, 10 years ago

Further code optimizations to PH weight and statistic update routines.

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