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

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

"Manually" controlling garbage collection in PH - Python does it stupidly and too often. Basically, there isn't any real good reason to do it more frequently than after every PH iteration, which is what I am now doing.

File size: 68.7 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      if self._verbose is True:
617         if self._scenario_tree._scenario_based_data == 1:
618            print "Scenario-based instance initialization enabled"
619         else:
620            print "Node-based instance initialization enabled"
621         
622      for scenario in self._scenario_tree._scenarios:
623
624         scenario_instance = construct_scenario_instance(self._scenario_tree,
625                                                         self._scenario_data_directory_name,
626                                                         scenario._name,
627                                                         self._model,
628                                                         self._verbose)
629
630         self._instances[scenario._name] = scenario_instance
631         self._instances[scenario._name].name = scenario._name
632         self._instance_augmented_attributes[scenario._name] = []
633
634      # let plugins know if they care - this callback point allows
635      # users to create/modify the original scenario instances and/or
636      # the scenario tree prior to creating PH-related parameters,
637      # variables, and the like.
638      for plugin in self._ph_plugins:     
639         plugin.post_instance_creation(self)         
640
641      # create ph-specific parameters (weights, xbar, etc.) for each instance.
642
643      if self._verbose is True:
644         print "Creating weight, average, and rho parameter vectors for scenario instances"
645
646      self._create_ph_scenario_parameters()
647           
648      # if specified, run the user script to initialize variable rhos at their whim.
649      if self._rho_setter is not None:
650         print "Executing user rho set script from filename="+self._rho_setter
651         execfile(self._rho_setter)
652
653      # with the instances created, run the user script to initialize variable bounds.
654      if self._bounds_setter is not None:
655         print "Executing user variable bounds set script from filename=", self._bounds_setter
656         execfile(self._bounds_setter)
657
658      # create parameters to store variable statistics (of general utility) at each node in the scenario tree.
659
660      if self._verbose is True:
661         print "Creating variable statistic (min/avg/max) parameter vectors for scenario tree nodes"
662
663      for stage in self._scenario_tree._stages[:-1]:
664
665         # first, gather all unique variables referenced in this stage
666         # this "gather" step is currently required because we're being lazy
667         # in terms of index management in the scenario tree - which
668         # should really be done in terms of sets of indices.
669
670         stage_variables = {}
671         for (reference_variable, index_template, reference_index) in stage._variables:
672            if reference_variable.name not in stage_variables.keys():
673               stage_variables[reference_variable.name] = reference_variable
674
675         # next, create min/avg/max parameters for each variable in the corresponding tree node.
676         # NOTE: the parameter names below could really be empty, as they are never referenced
677         #       explicitly.
678         for (variable_name, reference_variable) in stage_variables.items():
679            for tree_node in stage._tree_nodes:
680
681               new_min_index = reference_variable._index
682               new_min_parameter_name = "NODEMIN_"+reference_variable.name
683               # this bit of ugliness is due to Pyomo not correctly handling the Param construction
684               # case when the supplied index set consists strictly of None, i.e., the source variable
685               # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
686               new_min_parameter = None
687               if (len(new_min_index) is 1) and (None in new_min_index):
688                  new_min_parameter = Param(name=new_min_parameter_name)
689               else:
690                  new_min_parameter = Param(new_min_index,name=new_min_parameter_name)
691               for index in new_min_index:
692                  new_min_parameter[index] = 0.0
693               tree_node._minimums[reference_variable.name] = new_min_parameter                                   
694
695               new_avg_index = reference_variable._index
696               new_avg_parameter_name = "NODEAVG_"+reference_variable.name
697               new_avg_parameter = None
698               if (len(new_avg_index) is 1) and (None in new_avg_index):
699                  new_avg_parameter = Param(name=new_avg_parameter_name)
700               else:
701                  new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
702               for index in new_avg_index:
703                  new_avg_parameter[index] = 0.0                     
704               tree_node._averages[reference_variable.name] = new_avg_parameter
705
706               new_max_index = reference_variable._index
707               new_max_parameter_name = "NODEMAX_"+reference_variable.name
708               new_max_parameter = None
709               if (len(new_max_index) is 1) and (None in new_max_index):
710                  new_max_parameter = Param(name=new_max_parameter_name)
711               else:
712                  new_max_parameter = Param(new_max_index,name=new_max_parameter_name)
713               for index in new_max_index:
714                  new_max_parameter[index] = 0.0                     
715               tree_node._maximums[reference_variable.name] = new_max_parameter
716
717      # the objective functions are modified throughout the course of PH iterations.
718      # save the original, as a baseline to modify in subsequent iterations. reserve
719      # the original objectives, for subsequent modification.
720      self._original_objective_expression = {}
721      for instance_name, instance in self._instances.items():
722         objective_name = instance.active_components(Objective).keys()[0]
723         expr = instance.active_components(Objective)[objective_name]._data[None].expr
724         if isinstance(expr, Expression) is False:
725            expr = _IdentityExpression(expr)
726         self._original_objective_expression[instance_name] = expr
727
728      # cache the number of discrete and continuous variables in the master instance. this value
729      # is of general use, e.g., in the converger classes and in plugins.
730      (self._total_discrete_vars,self._total_continuous_vars) = self.compute_variable_counts()
731      if self._verbose is True:
732         print "Total number of discrete instance variables="+str(self._total_discrete_vars)
733         print "Total number of continuous instance variables="+str(self._total_continuous_vars)
734
735      # track the total number of fixed variables of each category at the end of each PH iteration.
736      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
737
738      # indicate that we're ready to run.
739      self._initialized = True
740
741      if self._verbose is True:
742         print "PH successfully created model instances for all scenarios"
743
744      self._init_end_time = time.time()
745
746      if self._verbose is True:
747         print "PH is successfully initialized"
748         if self._output_times is True:
749            print "Initialization time=%8.2f seconds" % (self._init_end_time - self._init_start_time)
750
751      # let plugins know if they care.
752      if self._verbose is True:
753         print "Initializing PH plugins"
754      for plugin in self._ph_plugins:
755         plugin.post_ph_initialization(self)
756      if self._verbose is True:
757         print "PH plugin initialization complete"
758
759   """ Perform the non-weighted scenario solves and form the initial w and xbars.
760   """   
761   def iteration_0_solve(self):
762
763      if self._verbose is True:
764         print "------------------------------------------------"
765         print "Starting PH iteration 0 solves"
766
767      self._current_iteration = 0
768
769      solve_start_time = time.time()
770
771      # STEP 0: set up all global solver options.
772      self._solver.mipgap = self._mipgap
773
774      # STEP 1: queue up the solves for all scenario sub-problems and
775      #         grab all the action handles for the subsequent barrier sync.
776
777      action_handles = []
778      scenario_action_handle_map = {} # maps scenario names to action handles
779      action_handle_scenario_map = {} # maps action handles to scenario names
780
781      for scenario in self._scenario_tree._scenarios:
782         
783         instance = self._instances[scenario._name]
784         
785         if self._verbose is True:
786            print "Queuing solve for scenario=" + scenario._name
787
788         # IMPT: You have to re-presolve if approximating continuous variable penalty terms with a
789         #       piecewise linear function. otherwise, the newly introduced variables won't be flagged
790         #       as unused (as is correct for iteration 0), and the output writer will crater.
791         # IMPT: I decided to presolve unconditionally, as PH extensions can add arbitrary components
792         #       to the base scenario instances - and the variable values/etc. need to be collectged.
793         instance.preprocess()
794           
795         # there's nothing to warm-start from in iteration 0, so don't include the keyword in the solve call.
796         # the reason you don't want to include it is that some solvers don't know how to handle the keyword
797         # at all (despite it being false). you might want to solve iteration 0 solves using some other solver.
798
799         new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
800         scenario_action_handle_map[scenario._name] = new_action_handle
801         action_handle_scenario_map[new_action_handle] = scenario._name
802
803         action_handles.append(new_action_handle)
804
805      # STEP 2: loop for the solver results, reading them and loading
806      #         them into instances as they are available.
807
808      if self._verbose is True:
809         print "Waiting for scenario sub-problem solves"
810
811      num_results_so_far = 0
812
813      while (num_results_so_far < len(self._scenario_tree._scenarios)):
814
815         action_handle = self._solver_manager.wait_any()
816         results = self._solver_manager.get_results(action_handle)         
817         scenario_name = action_handle_scenario_map[action_handle]
818         instance = self._instances[scenario_name]         
819
820         if self._verbose is True:
821            print "Results obtained for scenario="+scenario_name
822
823         if len(results.solution) == 0:
824            results.write(num=1)
825            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
826
827         if self._output_solver_results is True:
828            print "Results for scenario=",scenario_name
829            results.write(num=1)
830
831         start_time = time.time()
832         instance.load(results)
833         end_time = time.time()
834         if self._output_times is True:
835            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
836
837         if self._verbose is True:                 
838            print "Successfully loaded solution for scenario="+scenario_name
839
840         num_results_so_far = num_results_so_far + 1
841         
842      if self._verbose is True:     
843         print "Scenario sub-problem solves completed"     
844
845      solve_end_time = time.time()
846      self._cumulative_solve_time += (solve_end_time - solve_start_time)
847
848      if self._output_times is True:
849         print "Aggregate sub-problem solve time this iteration=%8.2f" % (solve_end_time - solve_start_time)
850
851      if self._verbose is True:
852         print "Successfully completed PH iteration 0 solves - solution statistics:"
853         print "         Scenario              Objective                  Value"
854         for scenario in self._scenario_tree._scenarios:
855            instance = self._instances[scenario._name]
856            for objective_name in instance.active_components(Objective):
857               objective = instance.active_components(Objective)[objective_name]
858               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
859         print "------------------------------------------------"
860
861   #
862   # recompute the averages, minimum, and maximum statistics for all variables to be blended by PH, i.e.,
863   # not appearing in the final stage. technically speaking, the min/max aren't required by PH, but they
864   # are used often enough to warrant their computation and it's basically free if you're computing the
865   # average.
866   #
867   def update_variable_statistics(self):
868
869      start_time = time.time()
870     
871      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage
872         
873         for tree_node in stage._tree_nodes:
874               
875            for (variable, index_template, variable_indices) in stage._variables:
876                 
877               variable_name = variable.name
878                     
879               avg_parameter_name = "PHAVG_"+variable_name
880                 
881               for index in variable_indices:
882                  min = float("inf")
883                  avg = 0.0
884                  max = float("-inf")
885                  node_probability = 0.0
886                     
887                  is_used = True # until proven otherwise                     
888                  for scenario in tree_node._scenarios:
889                       
890                     instance = self._instances[scenario._name]
891                       
892                     if getattr(instance,variable_name)[index].status == VarStatus.unused:
893                        is_used = False
894                     else:                       
895                        node_probability += scenario._probability
896                        var_value = getattr(instance, variable.name)[index].value
897                        if var_value < min:
898                           min = var_value
899                        avg += (scenario._probability * var_value)
900                        if var_value > max:
901                           max = var_value
902
903                  if is_used is True:
904
905                     tree_node._minimums[variable.name][index] = min
906                     tree_node._averages[variable.name][index] = avg / node_probability
907                     tree_node._maximums[variable.name][index] = max                       
908
909                     # distribute the newly computed average to the xbar variable in
910                     # each instance/scenario associated with this node. only do this
911                     # if the variable is used!
912                     for scenario in tree_node._scenarios:
913                        instance = self._instances[scenario._name]
914                        avg_parameter = getattr(instance, avg_parameter_name)
915                        avg_parameter[index] = avg / node_probability
916
917      end_time = time.time()
918      self._cumulative_xbar_time += (end_time - start_time)
919
920   def update_weights(self):
921
922      # because the weight updates rely on the xbars, and the xbars are node-based,
923      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
924      start_time = time.time()     
925
926      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.
927         
928         for tree_node in stage._tree_nodes:
929
930            for (variable, index_template, variable_indices) in stage._variables:
931
932               variable_name = variable.name
933               blend_parameter_name = "PHBLEND_"+variable_name
934               weight_parameter_name = "PHWEIGHT_"+variable_name
935               rho_parameter_name = "PHRHO_"+variable_name
936
937               for index in variable_indices:
938
939                  tree_node_average = tree_node._averages[variable.name][index]()
940
941                  for scenario in tree_node._scenarios:
942
943                     instance = self._instances[scenario._name]
944
945                     if getattr(instance,variable.name)[index].status != VarStatus.unused:
946
947                        # we are currently not updating weights if blending is disabled for a variable.
948                        # this is done on the premise that unless you are actively trying to move
949                        # the variable toward the mean, the weights will blow up and be huge by the
950                        # time that blending is activated.
951                        variable_blend_indicator = getattr(instance, blend_parameter_name)[index]()                             
952
953                        # get the weight and rho parameters for this variable/index combination.
954                        rho_value = getattr(instance, rho_parameter_name)[index]()
955                        current_variable_weight = getattr(instance, weight_parameter_name)[index]()
956                                               
957                        # if I'm maximizing, invert value prior to adding (hack to implement negatives).
958                        # probably fixed in Pyomo at this point - I just haven't checked in a long while.
959                        if self._is_minimizing is False:
960                           current_variable_weight = (-current_variable_weight)
961                        current_variable_value = getattr(instance,variable.name)[index]()
962                        new_variable_weight = current_variable_weight + variable_blend_indicator * rho_value * (current_variable_value - tree_node_average)
963                        # I have the correct updated value, so now invert if maximizing.
964                        if self._is_minimizing is False:
965                           new_variable_weight = (-new_variable_weight)
966                        getattr(instance, weight_parameter_name)[index].value = new_variable_weight
967
968      # we shouldn't have to re-simplify the expression, as we aren't adding any constant-variable terms - just modifying parameters.
969
970      end_time = time.time()
971      self._cumulative_weight_time += (end_time - start_time)
972
973   def form_iteration_k_objectives(self):
974
975      for instance_name, instance in self._instances.items():
976
977         new_attrs = form_ph_objective(instance_name, \
978                                       instance, \
979                                       self._original_objective_expression[instance_name], \
980                                       self._scenario_tree, \
981                                       self._linearize_nonbinary_penalty_terms, \
982                                       self._drop_proximal_terms, \
983                                       self._retain_quadratic_binary_terms, \
984                                       self._breakpoint_strategy, \
985                                       self._integer_tolerance)
986         self._instance_augmented_attributes[instance_name].extend(new_attrs)
987
988   def iteration_k_solve(self):
989
990      if self._verbose is True:
991         print "------------------------------------------------"       
992         print "Starting PH iteration " + str(self._current_iteration) + " solves"
993
994      # cache the objective values generated by PH for output at the end of this function.
995      ph_objective_values = {}
996
997      solve_start_time = time.time()
998
999      # STEP -1: if using a PH solver manager, propagate current weights/averages to the appropriate solver servers.
1000      #          ditto the tree node statistics, which are necessary if linearizing (so an optimization could be
1001      #          performed here).
1002      # NOTE: We aren't currently propagating rhos, as they generally don't change - we need to
1003      #       have a flag, though, indicating whether the rhos have changed, so they can be
1004      #       transmitted if needed.
1005      if self._solver_manager_type == "ph":
1006         self._transmit_weights_and_averages()
1007         self._transmit_tree_node_statistics()
1008
1009      # STEP 0: set up all global solver options.
1010      self._solver.mipgap = self._mipgap     
1011
1012      # STEP 1: queue up the solves for all scenario sub-problems and
1013      #         grab all of the action handles for the subsequent barrier sync.
1014
1015      action_handles = []
1016      scenario_action_handle_map = {} # maps scenario names to action handles
1017      action_handle_scenario_map = {} # maps action handles to scenario names
1018
1019      for scenario in self._scenario_tree._scenarios:     
1020
1021         instance = self._instances[scenario._name]
1022
1023         if self._verbose is True:
1024            print "Queuing solve for scenario=" + scenario._name
1025
1026         # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
1027         # don't do this, you won't see any chance in the output files as you vary the problem parameters!
1028         # ditto for instance fixing!
1029         instance.preprocess()
1030
1031         # once past iteration 0, there is always a feasible solution from which to warm-start.
1032         # however, you might want to disable warm-start when the solver is behaving badly (which does happen).
1033         new_action_handle = None
1034         if (self._disable_warmstarts is False) and (self._solver.warm_start_capable() is True):
1035            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
1036         else:
1037            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)           
1038
1039         scenario_action_handle_map[scenario._name] = new_action_handle
1040         action_handle_scenario_map[new_action_handle] = scenario._name         
1041
1042         action_handles.append(new_action_handle)
1043
1044      # STEP 2: loop for the solver results, reading them and loading
1045      #         them into instances as they are available.
1046      if self._verbose is True:
1047         print "Waiting for scenario sub-problem solves"
1048
1049      num_results_so_far = 0
1050
1051      while (num_results_so_far < len(self._scenario_tree._scenarios)):
1052
1053         action_handle = self._solver_manager.wait_any()
1054         results = self._solver_manager.get_results(action_handle)         
1055         scenario_name = action_handle_scenario_map[action_handle]
1056         instance = self._instances[scenario_name]         
1057
1058         if self._verbose is True:
1059            print "Results obtained for scenario="+scenario_name
1060
1061         if len(results.solution) == 0:
1062            results.write(num=1)
1063            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
1064
1065         if self._output_solver_results is True:
1066            print "Results for scenario=",scenario_name
1067            results.write(num=1)
1068
1069         start_time = time.time()
1070         instance.load(results)
1071         end_time = time.time()
1072         if self._output_times is True:
1073            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
1074
1075         if self._verbose is True:                 
1076            print "Successfully loaded solution for scenario="+scenario_name
1077
1078         # we're assuming there is a single solution.
1079         # the "value" attribute is a pre-defined feature of any solution - it is relative to whatever
1080         # objective was selected during optimization, which of course should be the PH objective.
1081         ph_objective_values[instance.name] = float(results.solution(0).objective['f'].value)
1082
1083         num_results_so_far = num_results_so_far + 1
1084         
1085      if self._verbose is True:     
1086         print "Scenario sub-problem solves completed"     
1087
1088      solve_end_time = time.time()
1089      self._cumulative_solve_time += (solve_end_time - solve_start_time)
1090
1091      if self._output_times is True:
1092         print "Aggregate sub-problem solve time this iteration=%8.2f" % (solve_end_time - solve_start_time)
1093
1094      if self._verbose is True:
1095         print "Successfully completed PH iteration " + str(self._current_iteration) + " solves - solution statistics:"
1096         print "  Scenario             PH Objective             Cost Objective"
1097         for scenario in self._scenario_tree._scenarios:
1098            instance = self._instances[scenario._name]
1099            for objective_name in instance.active_components(Objective):
1100               objective = instance.active_components(Objective)[objective_name]
1101               print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], self._scenario_tree.compute_scenario_cost(instance))
1102
1103   def solve(self):
1104
1105      self._solve_start_time = time.time()
1106      self._cumulative_solve_time = 0.0
1107      self._cumulative_xbar_time = 0.0
1108      self._cumulative_weight_time = 0.0
1109
1110      print "Starting PH"
1111
1112      if self._initialized == False:
1113         raise RuntimeError, "PH is not initialized - cannot invoke solve() method"
1114
1115      # garbage collection noticeably slows down PH when dealing with
1116      # large numbers of scenarios. fortunately, there are well-defined
1117      # points at which garbage collection makes sense (and there isn't a
1118      # lot of collection to do). namely, after each PH iteration.
1119      re_enable_gc = gc.isenabled()
1120      gc.disable()
1121
1122      print "Initiating PH iteration=" + `self._current_iteration`
1123
1124      self.iteration_0_solve()
1125
1126      # update variable statistics prior to any output.
1127      self.update_variable_statistics()     
1128
1129      if (self._verbose is True) or (self._report_solutions is True):
1130         print "Variable values following scenario solves:"
1131         self.pprint(False,False,True,False)
1132
1133      # let plugins know if they care.
1134      for plugin in self._ph_plugins:     
1135         plugin.post_iteration_0_solves(self)
1136
1137      # update the fixed variable statistics.
1138      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
1139
1140      if self._verbose is True:
1141         print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1142         print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1143
1144      # always output the convergence metric and first-stage cost statistics, to give a sense of progress.
1145      self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1146      first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()
1147      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)               
1148
1149      self.update_weights()
1150
1151      # let plugins know if they care.
1152      for plugin in self._ph_plugins:     
1153         plugin.post_iteration_0(self)
1154
1155      # if using a PH solver server, trasnsmit the rhos prior to the iteration
1156      # k solve sequence. for now, we are assuming that the rhos don't change
1157      # on a per-iteration basis, but that assumption can be easily relaxed.
1158      # it is important to do this after the plugins have a chance to do their
1159      # computation.
1160      if self._solver_manager_type == "ph":
1161         self._transmit_rhos()
1162         self._enable_ph_objectives()
1163
1164      # checkpoint if it's time - which it always is after iteration 0,
1165      # if the interval is >= 1!
1166      if (self._checkpoint_interval > 0):
1167         self.checkpoint(0)
1168
1169      # garbage-collect if it wasn't disabled entirely.
1170      if re_enable_gc is True:
1171         gc.collect()
1172
1173      # there is an upper bound on the number of iterations to execute -
1174      # the actual bound depends on the converger supplied by the user.
1175      for i in range(1, self._max_iterations+1):
1176
1177         self._current_iteration = self._current_iteration + 1                 
1178
1179         print "Initiating PH iteration=" + `self._current_iteration`         
1180
1181         if (self._verbose is True) or (self._report_weights is True):
1182            print "Variable averages and weights prior to scenario solves:"
1183            self.pprint(True,True,False,False)
1184
1185         # with the introduction of piecewise linearization, the form of the
1186         # penalty-weighted objective is no longer fixed. thus, we need to
1187         # create the objectives each PH iteration.
1188         self.form_iteration_k_objectives()
1189
1190         # let plugins know if they care.
1191         for plugin in self._ph_plugins:
1192            plugin.pre_iteration_k_solves(self)
1193
1194         # do the actual solves.
1195         self.iteration_k_solve()
1196
1197         # update variable statistics prior to any output.
1198         self.update_variable_statistics()
1199         
1200         if (self._verbose is True) or (self._report_solutions is True):
1201            print "Variable values following scenario solves:"
1202            self.pprint(False,False,True,False)
1203
1204         # we don't technically have to do this at the last iteration,
1205         # but with checkpointing and re-starts, you're never sure
1206         # when you're executing the last iteration.
1207         self.update_weights()
1208
1209         # let plugins know if they care.
1210         for plugin in self._ph_plugins:
1211            plugin.post_iteration_k_solves(self)
1212
1213         # update the fixed variable statistics.
1214         (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
1215
1216         if self._verbose is True:
1217            print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1218            print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1219
1220         # let plugins know if they care.
1221         for plugin in self._ph_plugins:
1222            plugin.post_iteration_k(self)
1223
1224         # at this point, all the real work of an iteration is complete.
1225
1226         # checkpoint if it's time.
1227         if (self._checkpoint_interval > 0) and (i % self._checkpoint_interval is 0):
1228            self.checkpoint(i)
1229
1230         # check for early termination.
1231         self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1232         first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()         
1233         print "Convergence metric=%12.4f  First stage cost avg=%12.4f  Max-Min=%8.2f" % (self._converger.lastMetric(), first_stage_avg, first_stage_max-first_stage_min)         
1234
1235         if self._converger.isConverged(self) is True:
1236            if self._total_discrete_vars == 0:
1237               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)
1238            else:
1239               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)+" or all discrete variables are fixed"               
1240            break
1241
1242         # if we're terminating due to exceeding the maximum iteration count, print a message
1243         # indicating so - otherwise, you get a quiet, information-free output trace.
1244         if i == self._max_iterations:
1245            print "Halting PH - reached maximal iteration count="+str(self._max_iterations)
1246
1247         # garbage-collect if it wasn't disabled entirely.
1248         if re_enable_gc is True:
1249            gc.collect()
1250
1251      # re-enable the normal garbage collection mode.
1252      if re_enable_gc is True:
1253         gc.enable()
1254
1255      if self._verbose is True:
1256         print "Number of discrete variables fixed before final plugin calls="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1257         print "Number of continuous variables fixed before final plugin calls="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"           
1258
1259      # let plugins know if they care. do this before
1260      # the final solution / statistics output, as the plugins
1261      # might do some final tweaking.
1262      for plugin in self._ph_plugins:
1263         plugin.post_ph_execution(self)
1264
1265      # update the fixed variable statistics - the plugins might have done something.
1266      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()                       
1267
1268      print "PH complete"
1269
1270      print "Convergence history:"
1271      self._converger.pprint()
1272
1273      print "Final number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1274      print "Final number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"         
1275
1276      print "Final variable values:"
1277      self.pprint(False,False,True,True)         
1278
1279      print "Final costs:"
1280      self._scenario_tree.pprintCosts(self._instances)
1281
1282      self._solve_end_time = time.time()
1283
1284      if (self._verbose is True) and (self._output_times is True):
1285         print "Overall run-time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1286
1287      # cleanup the scenario instances for post-processing - ideally, we want to leave them in
1288      # their original state, minus all the PH-specific stuff. we don't do all cleanup (leaving
1289      # things like rhos, etc), but we do clean up constraints, as that really hoses up the ef writer.
1290      self._cleanup_scenario_instances()
1291
1292   #
1293   # prints a summary of all collected time statistics
1294   #
1295   def print_time_stats(self):
1296
1297      print "PH run-time statistics (user):"
1298
1299      print "Initialization time=  %8.2f seconds" % (self._init_end_time - self._init_start_time)
1300      print "Overall solve time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1301      print "Scenario solve time=  %8.2f seconds" % self._cumulative_solve_time
1302      print "Average update time=  %8.2f seconds" % self._cumulative_xbar_time
1303      print "Weight update time=   %8.2f seconds" % self._cumulative_weight_time
1304
1305   #
1306   # a utility to determine whether to output weight / average / etc. information for
1307   # a variable/node combination. when the printing is moved into a callback/plugin,
1308   # this routine will go there. for now, we don't dive down into the node resolution -
1309   # just the variable/stage.
1310   #
1311   def should_print(self, stage, variable, variable_indices):
1312
1313      if self._output_continuous_variable_stats is False:
1314
1315         variable_type = variable.domain         
1316
1317         if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):
1318
1319            return False
1320
1321      return True
1322     
1323   #
1324   # pretty-prints the state of the current variable averages, weights, and values.
1325   # inputs are booleans indicating which components should be output.
1326   #
1327   def pprint(self, output_averages, output_weights, output_values, output_fixed):
1328
1329      if self._initialized is False:
1330         raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"         
1331     
1332      # print tree nodes and associated variable/xbar/ph information in stage-order
1333      # we don't blend in the last stage, so we don't current care about printing the associated information.
1334      for stage in self._scenario_tree._stages[:-1]:
1335
1336         print "\tStage=" + stage._name
1337
1338         num_outputs_this_stage = 0 # tracks the number of outputs on a per-index basis.
1339
1340         for (variable, index_template, variable_indices) in stage._variables:
1341
1342            variable_name = variable.name
1343
1344            if self.should_print(stage, variable, variable_indices) is True:
1345
1346               num_outputs_this_variable = 0 # track, so we don't output the variable names unless there is an entry to report.
1347
1348               for index in variable_indices:               
1349
1350                  weight_parameter_name = "PHWEIGHT_"+variable_name
1351
1352                  num_outputs_this_index = 0 # track, so we don't output the variable index more than once.
1353
1354                  for tree_node in stage._tree_nodes:                 
1355
1356                     # determine if the variable/index pair is used across the set of scenarios (technically,
1357                     # it should be good enough to check one scenario). ditto for "fixed" status. fixed does
1358                     # imply unused (see note below), but we care about the fixed status when outputting
1359                     # final solutions.
1360
1361                     is_used = True # should be consistent across scenarios, so one "unused" flags as invalid.
1362                     is_fixed = False
1363
1364                     for scenario in tree_node._scenarios:
1365                        instance = self._instances[scenario._name]
1366                        variable_value = getattr(instance,variable_name)[index]
1367                        if variable_value.status == VarStatus.unused:
1368                           is_used = False
1369                        if variable_value.fixed is True:
1370                           is_fixed = True
1371 
1372                     # IMPT: this is far from obvious, but variables that are fixed will - because
1373                     #       presolve will identify them as constants and eliminate them from all
1374                     #       expressions - be flagged as "unused" and therefore not output.
1375
1376                     if ((output_fixed is True) and (is_fixed is True)) or (is_used is True):
1377
1378                           minimum_value = None
1379                           maximum_value = None
1380
1381                           if index is None:
1382                              minimum_value = tree_node._minimums[variable_name]
1383                              maximum_value = tree_node._maximums[variable_name]
1384                           else:
1385                              minimum_value = tree_node._minimums[variable_name][index]()
1386                              maximum_value = tree_node._maximums[variable_name][index]()
1387
1388                           # there really isn't a need to output variables whose
1389                           # values are equal to 0 across-the-board. and there is
1390                           # good reason not to, i.e., the volume of output.
1391                           if (fabs(minimum_value) > self._integer_tolerance) or \
1392                              (fabs(maximum_value) > self._integer_tolerance):
1393
1394                              num_outputs_this_stage = num_outputs_this_stage + 1                           
1395                              num_outputs_this_variable = num_outputs_this_variable + 1
1396                              num_outputs_this_index = num_outputs_this_index + 1
1397
1398                              if num_outputs_this_variable == 1:
1399                                 print "\t\tVariable=" + variable_name
1400
1401                              if num_outputs_this_index == 1:
1402                                 if index is not None:
1403                                    print "\t\t\tIndex:", indexToString(index)                             
1404
1405                              print "\t\t\t\tTree Node="+tree_node._name+"\t\t (Scenarios: ",                             
1406                              for scenario in tree_node._scenarios:
1407                                 print scenario._name," ",
1408                                 if scenario == tree_node._scenarios[-1]:
1409                                    print ")"
1410                           
1411                              if output_values is True:
1412                                 average_value = tree_node._averages[variable_name][index]()
1413                                 print "\t\t\t\tValues: ",                       
1414                                 for scenario in tree_node._scenarios:
1415                                    instance = self._instances[scenario._name]
1416                                    this_value = getattr(instance,variable_name)[index].value
1417                                    print "%12.4f" % this_value,
1418                                    if scenario == tree_node._scenarios[-1]:
1419                                       print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1420                                       print "    Avg=%12.4f" % (average_value),
1421                                       print ""
1422                              if output_weights:
1423                                 print "\t\t\t\tWeights: ",
1424                                 for scenario in tree_node._scenarios:
1425                                    instance = self._instances[scenario._name]
1426                                    print "%12.4f" % getattr(instance,weight_parameter_name)[index].value,
1427                                    if scenario == tree_node._scenarios[-1]:
1428                                       print ""
1429
1430                              if output_averages:
1431                                 print "\t\t\t\tAverage: %12.4f" % (tree_node._averages[variable_name][index].value)
1432
1433         if num_outputs_this_stage == 0:
1434            print "\t\tNo non-converged variables in stage"
1435
1436         # cost variables aren't blended, so go through the gory computation of min/max/avg.
1437         # we currently always print these.
1438         cost_variable_name = stage._cost_variable[0].name
1439         cost_variable_index = stage._cost_variable[1]
1440         if cost_variable_index is None:
1441            print "\t\tCost Variable=" + cost_variable_name
1442         else:
1443            print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)           
1444         for tree_node in stage._tree_nodes:
1445            print "\t\t\tTree Node=" + tree_node._name + "\t\t (Scenarios: ",
1446            for scenario in tree_node._scenarios:
1447               print scenario._name," ",
1448               if scenario == tree_node._scenarios[-1]:
1449                  print ")"
1450            maximum_value = 0.0
1451            minimum_value = 0.0
1452            sum_values = 0.0
1453            num_values = 0
1454            first_time = True
1455            print "\t\t\tValues: ",                       
1456            for scenario in tree_node._scenarios:
1457                instance = self._instances[scenario._name]
1458                this_value = getattr(instance,cost_variable_name)[cost_variable_index].value
1459                print "%12.4f" % this_value,
1460                num_values += 1
1461                sum_values += this_value
1462                if first_time is True:
1463                   first_time = False
1464                   maximum_value = this_value
1465                   minimum_value = this_value
1466                else:
1467                   if this_value > maximum_value:
1468                      maximum_value = this_value
1469                   if this_value < minimum_value:
1470                      minimum_value = this_value
1471                if scenario == tree_node._scenarios[-1]:
1472                   print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1473                   print "    Avg=%12.4f" % (sum_values/num_values),
1474                   print ""
1475           
Note: See TracBrowser for help on using the repository browser.