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

Last change on this file since 2446 was 2446, checked in by jwatson, 9 years ago

Added --output-scenario-tree-solution option to PySP runph script, which will:
1) Create a solution from the node averages in a scenario tree -and-
2) Output the full solution in scenario tree format (which includes the leaf nodes).

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