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

Last change on this file since 2406 was 2406, checked in by jwatson, 11 years ago

More tweaks to the PH solver manager/servers.

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