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

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

Drive-by cleanup of PH option parsing, implement logging facility, mild rework of imports

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