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

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

Disabling garbage collection in PySP during instance construction.

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