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

#  _________________________________________________________________________
#
#  Coopr: A COmmon Optimization Python Repository
#  Copyright (c) 2008 Sandia Corporation.
#  This software is distributed under the BSD License.
#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
#  the U.S. Government retains certain rights in this software.
#  For more information, see the Coopr README.txt file.
#  _________________________________________________________________________

import sys
import types
from coopr.pyomo import *
from coopr.pyomo.base.expr import *
import copy
import os.path
import traceback
import copy
from coopr.opt import SolverResults,SolverStatus
from coopr.opt.base import SolverFactory
from coopr.opt.parallel import SolverManagerFactory
import time
import types
import pickle
import gc
from math import fabs, log, exp

from scenariotree import *
from phutils import *
from phobjective import *

from pyutilib.component.core import ExtensionPoint

from coopr.pysp.phextension import IPHExtension

class ProgressiveHedging(object):

   #
   # a utility intended for folks who are brave enough to script rho setting in a python file.
   #

   def setRhoAllScenarios(self, variable_value, rho_expression):

      variable_name = None
      variable_index = None

      if isVariableNameIndexed(variable_value.name) is True:

         variable_name, variable_index = extractVariableNameAndIndex(variable_value.name)

      else:

         variable_name = variable_value.name
         variable_index = None

      new_rho_value = None
      if isinstance(rho_expression, float):
         new_rho_value = rho_expression
      else:
         new_rho_value = rho_expression()

      if self._verbose is True:
         print "Setting rho="+str(new_rho_value)+" for variable="+variable_value.name

      for instance_name, instance in self._instances.items():

         rho_param = getattr(instance, "PHRHO_"+variable_name)
         rho_param[variable_index] = new_rho_value

   #
   # a utility intended for folks who are brave enough to script variable bounds setting in a python file.
   #

   def setVariableBoundsAllScenarios(self, variable_name, variable_index, lower_bound, upper_bound):

      if isinstance(lower_bound, float) is False:
         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)

      if isinstance(upper_bound, float) is False:
         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)

      for instance_name, instance in self._instances.items():

         variable = getattr(instance, variable_name)
         variable[variable_index].setlb(lower_bound)
         variable[variable_index].setub(upper_bound)

   #
   # a utility intended for folks who are brave enough to script variable bounds setting in a python file.
   # same functionality as above, but applied to all indicies of the variable, in all scenarios.
   #

   def setVariableBoundsAllIndicesAllScenarios(self, variable_name, lower_bound, upper_bound):

      if isinstance(lower_bound, float) is False:
         raise ValueError, "Lower bound supplied to PH method setVariableBoundsAllIndiciesAllScenarios for variable="+variable_name+" must be a constant; value supplied="+str(lower_bound)

      if isinstance(upper_bound, float) is False:
         raise ValueError, "Upper bound supplied to PH method setVariableBoundsAllIndicesAllScenarios for variable="+variable_name+" must be a constant; value supplied="+str(upper_bound)

      for instance_name, instance in self._instances.items():

         variable = getattr(instance, variable_name)
         for index in variable:
            variable[index].setlb(lower_bound)
            variable[index].setub(upper_bound)

   #
   # checkpoint the current PH state via pickle'ing. the input iteration count
   # simply serves as a tag to create the output file name. everything with the
   # exception of the _ph_plugins, _solver_manager, and _solver attributes are
   # pickled. currently, plugins fail in the pickle process, which is fine as
   # JPW doesn't think you want to pickle plugins (particularly the solver and
   # solver manager) anyway. For example, you might want to change those later,
   # after restoration - and the PH state is independent of how scenario
   # sub-problems are solved.
   #

   def checkpoint(self, iteration_count):

      checkpoint_filename = "checkpoint."+str(iteration_count)

      tmp_ph_plugins = self._ph_plugins
      tmp_solver_manager = self._solver_manager
      tmp_solver = self._solver

      self._ph_plugins = None
      self._solver_manager = None
      self._solver = None

      checkpoint_file = open(checkpoint_filename, "w")
      pickle.dump(self,checkpoint_file)
      checkpoint_file.close()

      self._ph_plugins = tmp_ph_plugins
      self._solver_manager = tmp_solver_manager
      self._solver = tmp_solver

      print "Checkpoint written to file="+checkpoint_filename

   #
   # a simple utility to count the number of continuous and discrete variables in a set of instances.
   # unused variables are ignored, and counts include all active indices. returns a pair - num-discrete,
   # num-continuous.
   #

   def compute_variable_counts(self):

      num_continuous_vars = 0
      num_discrete_vars = 0

      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage

         for tree_node in stage._tree_nodes:

            for (variable, index_template) in stage._variables:

               variable_indices = tree_node._variable_indices[variable.name]

               for index in variable_indices:

                  is_used = True # until proven otherwise
                  for scenario in tree_node._scenarios:
                     instance = self._instances[scenario._name]
                     if getattr(instance,variable.name)[index].status == VarStatus.unused:
                        is_used = False

                  if is_used is True:

                     if isinstance(variable.domain, IntegerSet) or isinstance(variable.domain, BooleanSet):
                        num_discrete_vars = num_discrete_vars + 1
                     else:
                        num_continuous_vars = num_continuous_vars + 1

      return (num_discrete_vars, num_continuous_vars)

   #
   # ditto above, but count the number of fixed discrete and continuous variables.
   # important: once a variable (value) is fixed, it is flagged as unused in the
   # course of presolve - because it is no longer referenced. this makes sense,
   # of course; it's just something to watch for. this is an obvious assumption
   # that we won't be fixing unused variables, which should not be an issue.
   #

   def compute_fixed_variable_counts(self):

      num_fixed_continuous_vars = 0
      num_fixed_discrete_vars = 0

      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage

         for tree_node in stage._tree_nodes:

            for (variable, index_template) in stage._variables:

               variable_indices = tree_node._variable_indices[variable.name]

               for index in variable_indices:

                  # implicit assumption is that if a variable value is fixed in one
                  # scenario, it is fixed in all scenarios.

                  is_fixed = False # until proven otherwise
                  for scenario in tree_node._scenarios:
                     instance = self._instances[scenario._name]
                     var_value = getattr(instance,variable.name)[index]
                     if var_value.fixed is True:
                        is_fixed = True

                  if is_fixed is True:

                     if isinstance(variable.domain, IntegerSet) or isinstance(variable.domain, BooleanSet):
                        num_fixed_discrete_vars = num_fixed_discrete_vars + 1
                     else:
                        num_fixed_continuous_vars = num_fixed_continuous_vars + 1

      return (num_fixed_discrete_vars, num_fixed_continuous_vars)

   # when the quadratic penalty terms are approximated via piecewise linear segments,
   # we end up (necessarily) "littering" the scenario instances with extra constraints.
   # these need to and should be cleaned up after PH, for purposes of post-PH manipulation,
   # e.g., writing the extensive form.
   def _cleanup_scenario_instances(self):

      for instance_name, instance in self._instances.items():

         for constraint_name in self._instance_augmented_attributes[instance_name]:

            instance._clear_attribute(constraint_name)

         # if you don't pre-solve, the name collections won't be updated.
         instance.preprocess()

   # create PH weight and xbar vectors, on a per-scenario basis, for each variable that is not in the
   # final stage, i.e., for all variables that are being blended by PH. the parameters are created
   # in the space of each scenario instance, so that they can be directly and automatically
   # incorporated into the (appropriately modified) objective function.

   def _create_ph_scenario_parameters(self):

      for (instance_name, instance) in self._instances.items():

         new_penalty_variable_names = create_ph_parameters(instance, self._scenario_tree, self._rho, self._linearize_nonbinary_penalty_terms)
         self._instance_augmented_attributes[instance_name].extend(new_penalty_variable_names)

   # a simple utility to extract the first-stage cost statistics, e.g., min, average, and max.
   def _extract_first_stage_cost_statistics(self):

      maximum_value = 0.0
      minimum_value = 0.0
      sum_values = 0.0
      num_values = 0
      first_time = True

      first_stage = self._scenario_tree._stages[0]
      (cost_variable, cost_variable_index) = first_stage._cost_variable
      for scenario_name, instance in self._instances.items():
         this_value = getattr(instance, cost_variable.name)[cost_variable_index].value
         if this_value is not None: # None means not reported by the solver.
            num_values += 1
            sum_values += this_value
            if first_time is True:
               first_time = False
               maximum_value = this_value
               minimum_value = this_value
            else:
               if this_value > maximum_value:
                  maximum_value = this_value
               if this_value < minimum_value:
                  minimum_value = this_value

      if num_values > 0:
         sum_values = sum_values / num_values

      return minimum_value, sum_values, maximum_value

   # a utility to transmit - across the PH solver manager - the current weights
   # and averages for each of my problem instances. used to set up iteration K solves.
   def _transmit_weights_and_averages(self):

      for scenario_name, scenario_instance in self._instances.items():

         weights_to_transmit = []
         averages_to_transmit = []

         for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.

            for (variable, index_template) in stage._variables:

               variable_name = variable.name
               weight_parameter_name = "PHWEIGHT_"+variable_name
               weights_to_transmit.append(getattr(scenario_instance, weight_parameter_name))
               average_parameter_name = "PHAVG_"+variable_name
               averages_to_transmit.append(getattr(scenario_instance, average_parameter_name))

         self._solver_manager.transmit_weights_and_averages(scenario_instance, weights_to_transmit, averages_to_transmit)

   #
   # a utility to transmit - across the PH solver manager - the current rho values
   # for each of my problem instances. mainly after PH iteration 0 is complete,
   # in preparation for the iteration K solves.
   #

   def _transmit_rhos(self):

      for scenario_name, scenario_instance in self._instances.items():

         rhos_to_transmit = []

         for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no rhos to worry about.

            for (variable, index_template) in stage._variables:

               variable_name = variable.name
               rho_parameter_name = "PHRHO_"+variable_name
               rhos_to_transmit.append(getattr(scenario_instance, rho_parameter_name))

         self._solver_manager.transmit_rhos(scenario_instance, rhos_to_transmit)

   #
   # a utility to transmit - across the PH solver manager - the current scenario
   # tree node statistics to each of my problem instances. done prior to each
   # PH iteration k.
   #

   def _transmit_tree_node_statistics(self):

      for scenario_name, scenario_instance in self._instances.items():

         tree_node_minimums = {}
         tree_node_maximums = {}

         scenario = self._scenario_tree._scenario_map[scenario_name]

         for tree_node in scenario._node_list:

            tree_node_minimums[tree_node._name] = tree_node._minimums
            tree_node_maximums[tree_node._name] = tree_node._maximums

         self._solver_manager.transmit_tree_node_statistics(scenario_instance, tree_node_minimums, tree_node_maximums)

   #
   # a utility to enable - across the PH solver manager - weighted penalty objectives.
   #

   def _enable_ph_objectives(self):

      for scenario_name, scenario_instance in self._instances.items():

         self._solver_manager.enable_ph_objective(scenario_instance)


   """ Constructor
       Arguments:
          max_iterations        the maximum number of iterations to run PH (>= 0). defaults to 0.
          rho                   the global rho value (> 0). defaults to 0.
          rho_setter            an optional name of a python file used to set particular variable rho values.
          solver                the solver type that PH uses to solve scenario sub-problems. defaults to "cplex".
          solver_manager        the solver manager type that coordinates scenario sub-problem solves. defaults to "serial".
          keep_solver_files     do I keep intermediate solver files around (for debugging)? defaults to False.
          output_solver_log     do I dump the solver log (as it is being generated) to the screen? defaults to False.
          output_solver_results do I output (for debugging) the detailed solver results, including solutions, for scenario solves? defaults to False.
          verbose               does the PH object stream debug/status output? defaults to False.
          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).
          checkpoint_interval   how many iterations between writing a checkpoint file containing the entire PH state? defaults to 0, indicating never.

   """
   def __init__(self, *args, **kwds):

      # PH configuration parameters
      self._rho = 0.0 # a default, global values for rhos.
      self._rho_setter = None # filename for the modeler to set rho on a per-variable or per-scenario basis.
      self._bounds_setter = None # filename for the modeler to set rho on a per-variable basis, after all scenarios are available.
      self._max_iterations = 0

      # PH reporting parameters
      self._verbose = False # do I flood the screen with status output?
      self._report_solutions = False # do I report solutions after each PH iteration?
      self._report_weights = False # do I report PH weights prior to each PH iteration?
      self._report_only_statistics = False # do I report only variable statistics when outputting solutions and weights?
      self._output_continuous_variable_stats = True # when in verbose mode, do I output weights/averages for continuous variables?
      self._output_solver_results = False
      self._output_times = False
      self._output_scenario_tree_solution = False

      # PH run-time variables
      self._current_iteration = 0 # the 'k'
      self._xbar = {} # current per-variable averages. maps (node_id, variable_name) -> value
      self._initialized = False # am I ready to call "solve"? Set to True by the initialize() method.

      # PH solver information / objects.
      self._solver_type = "cplex"
      self._solver_manager_type = "serial" # serial or pyro are the options currently available

      self._solver = None
      self._solver_manager = None

      self._keep_solver_files = False
      self._output_solver_log = False

      # PH convergence computer/updater.
      self._converger = None

      # PH history
      self._solutions = {}

      # the checkpoint interval - expensive operation, but worth it for big models.
      # 0 indicates don't checkpoint.
      self._checkpoint_interval = 0

      # all information related to the scenario tree (implicit and explicit).
      self._model = None # not instantiated
      self._model_instance = None # instantiated

      self._scenario_tree = None

      self._scenario_data_directory = "" # this the prefix for all scenario data
      self._instances = {} # maps scenario name to the corresponding model instance
      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.

      # for various reasons (mainly hacks at this point), it's good to know whether we're minimizing or maximizing.
      self._is_minimizing = None

      # global handle to ph extension plugins
      self._ph_plugins = ExtensionPoint(IPHExtension)

      # PH timing statistics - relative to last invocation.
      self._init_start_time = None # for initialization() method
      self._init_end_time = None
      self._solve_start_time = None # for solve() method
      self._solve_end_time = None
      self._cumulative_solve_time = None # seconds, over course of solve()
      self._cumulative_xbar_time = None # seconds, over course of update_xbars()
      self._cumulative_weight_time = None # seconds, over course of update_weights()

      # do I disable warm-start for scenario sub-problem solves during PH iterations >= 1?
      self._disable_warmstarts = False

      # do I drop proximal (quadratic penalty) terms from the weighted objective functions?
      self._drop_proximal_terms = False

      # do I linearize the quadratic penalty term for continuous variables via a
      # piecewise linear approximation? the default should always be 0 (off), as the
      # user should be aware when they are forcing an approximation.
      self._linearize_nonbinary_penalty_terms = 0

      # the breakpoint distribution strategy employed when linearizing. 0 implies uniform
      # distribution between the variable lower and upper bounds.
      self._breakpoint_strategy = 0

      # do I retain quadratic objective terms associated with binary variables? in general,
      # there is no good reason to not linearize, but just in case, I introduced the option.
      self._retain_quadratic_binary_terms = False

      # PH default tolerances - for use in fixing and testing equality across scenarios,
      # and other stuff.
      self._integer_tolerance = 0.00001

      # PH maintains a mipgap that is applied to each scenario solve that is performed.
      # this attribute can be changed by PH extensions, and the change will be applied
      # on all subsequent solves - until it is modified again. the default is None,
      # indicating unassigned.
      self._mipgap = None

      # should PH simplify expressions after creating them? an example includes
      # the objective function expression, either with or without linearization.
      self._simplify_expressions = False

      # we only store these temporarily...
      scenario_solver_options = None

      # process the keyword options
      for key in kwds.keys():
         if key == "max_iterations":
            self._max_iterations = kwds[key]
         elif key == "rho":
            self._rho = kwds[key]
         elif key == "rho_setter":
            self._rho_setter = kwds[key]
         elif key == "bounds_setter":
            self._bounds_setter = kwds[key]
         elif key == "solver":
            self._solver_type = kwds[key]
         elif key == "solver_manager":
            self._solver_manager_type = kwds[key]
         elif key == "scenario_solver_options":
            scenario_solver_options = kwds[key]
         elif key == "scenario_mipgap":
            self._mipgap = kwds[key]
         elif key == "keep_solver_files":
            self._keep_solver_files = kwds[key]
         elif key == "output_solver_results":
            self._output_solver_results = kwds[key]
         elif key == "output_solver_log":
            self._output_solver_log = kwds[key]
         elif key == "verbose":
            self._verbose = kwds[key]
         elif key == "report_solutions":
            self._report_solutions = kwds[key]
         elif key == "report_weights":
            self._report_weights = kwds[key]
         elif key == "report_only_statistics":
            self._report_only_statistics = kwds[key]
         elif key == "output_times":
            self._output_times = kwds[key]
         elif key == "disable_warmstarts":
            self._disable_warmstarts = kwds[key]
         elif key == "drop_proximal_terms":
            self._drop_proximal_terms = kwds[key]
         elif key == "retain_quadratic_binary_terms":
            self._retain_quadratic_binary_terms = kwds[key]
         elif key == "linearize_nonbinary_penalty_terms":
            self._linearize_nonbinary_penalty_terms = kwds[key]
         elif key == "breakpoint_strategy":
            self._breakpoint_strategy = kwds[key]
         elif key == "checkpoint_interval":
            self._checkpoint_interval = kwds[key]
         elif key == "output_scenario_tree_solution":
            self._output_scenario_tree_solution = kwds[key]
         elif key == "simplify_expressions":
            self._simplify_expressions = kwds[key]
         else:
            print "Unknown option=" + key + " specified in call to PH constructor"

      # validate all "atomic" options (those that can be validated independently)
      if self._max_iterations < 0:
         raise ValueError, "Maximum number of PH iterations must be non-negative; value specified=" + `self._max_iterations`
      if self._rho <= 0.0:
         raise ValueError, "Value of the rho parameter in PH must be non-zero positive; value specified=" + `self._rho`
      if (self._mipgap is not None) and ((self._mipgap < 0.0) or (self._mipgap > 1.0)):
         raise ValueError, "Value of the mipgap parameter in PH must be on the unit interval; value specified=" + `self._mipgap`

      # validate the linearization (number of pieces) and breakpoint distribution parameters.
      if self._linearize_nonbinary_penalty_terms < 0:
         raise ValueError, "Value of linearization parameter for nonbinary penalty terms must be non-negative; value specified=" + `self._linearize_nonbinary_penalty_terms`
      if self._breakpoint_strategy < 0:
         raise ValueError, "Value of the breakpoint distribution strategy parameter must be non-negative; value specified=" + str(self._breakpoint_strategy)
      if self._breakpoint_strategy > 3:
         raise ValueError, "Unknown breakpoint distribution strategy specified - valid values are between 0 and 2, inclusive; value specified=" + str(self._breakpoint_strategy)

      # validate rho setter file if specified.
      if self._rho_setter is not None:
         if os.path.exists(self._rho_setter) is False:
            raise ValueError, "The rho setter script file="+self._rho_setter+" does not exist"

      # validate bounds setter file if specified.
      if self._bounds_setter is not None:
         if os.path.exists(self._bounds_setter) is False:
            raise ValueError, "The bounds setter script file="+self._bounds_setter+" does not exist"

      # validate the checkpoint interval.
      if self._checkpoint_interval < 0:
         raise ValueError, "A negative checkpoint interval with value="+str(self._checkpoint_interval)+" was specified in call to PH constructor"

      # construct the sub-problem solver.
      if self._verbose is True:
         print "Constructing solver type="+self._solver_type
      self._solver = SolverFactory(self._solver_type)
      if self._solver == None:
         raise ValueError, "Unknown solver type=" + self._solver_type + " specified in call to PH constructor"
      if self._keep_solver_files is True:
         self._solver.keepFiles = True
      if len(scenario_solver_options) > 0:
         if self._verbose is True:
            print "Initializing scenario sub-problem solver with options="+str(scenario_solver_options)
         self._solver.set_options("".join(scenario_solver_options))
      if self._output_times is True:
         self._solver._report_timing = True

      # construct the solver manager.
      if self._verbose is True:
         print "Constructing solver manager of type="+self._solver_manager_type
      self._solver_manager = SolverManagerFactory(self._solver_manager_type)
      if self._solver_manager is None:
         raise ValueError, "Failed to create solver manager of type="+self._solver_manager_type+" specified in call to PH constructor"

      # a set of all valid PH iteration indicies is generally useful for plug-ins, so create it here.
      self._iteration_index_set = Set(name="PHIterations")
      for i in range(0,self._max_iterations + 1):
         self._iteration_index_set.add(i)

      # spit out parameterization if verbosity is enabled
      if self._verbose is True:
         print "PH solver configuration: "
         print "   Max iterations=" + `self._max_iterations`
         print "   Default global rho=" + `self._rho`
         if self._rho_setter is not None:
            print "   Rho initialization file=" + self._rho_setter
         if self._bounds_setter is not None:
            print "   Variable bounds initialization file=" + self._bounds_setter
         print "   Sub-problem solver type=" + `self._solver_type`
         print "   Solver manager type=" + `self._solver_manager_type`
         print "   Keep solver files? " + str(self._keep_solver_files)
         print "   Output solver results? " + str(self._output_solver_results)
         print "   Output solver log? " + str(self._output_solver_log)
         print "   Output times? " + str(self._output_times)
         print "   Checkpoint interval="+str(self._checkpoint_interval)

   """ Initialize PH with model and scenario data, in preparation for solve().
       Constructs and reads instances.
   """
   def initialize(self, scenario_data_directory_name=".", model=None, model_instance=None, scenario_tree=None, converger=None, linearize=False):

      self._init_start_time = time.time()

      if self._verbose is True:
         print "Initializing PH"
         print "   Scenario data directory=" + scenario_data_directory_name

      if not os.path.exists(scenario_data_directory_name):
         raise ValueError, "Scenario data directory=" + scenario_data_directory_name + " either does not exist or cannot be read"

      self._scenario_data_directory_name = scenario_data_directory_name

      # IMPT: The input model should be an *instance*, as it is very useful (critical!) to know
      #       the dimensions of sets, be able to store suffixes on variable values, etc.
      if model is None:
         raise ValueError, "A model must be supplied to the PH initialize() method"

      if scenario_tree is None:
         raise ValueError, "A scenario tree must be supplied to the PH initialize() method"

      if converger is None:
         raise ValueError, "A convergence computer must be supplied to the PH initialize() method"

      self._model = model
      self._model_instance = model_instance
      self._scenario_tree = scenario_tree
      self._converger = converger

      model_objective = model.active_components(Objective)
      self._is_minimizing = (model_objective[ model_objective.keys()[0] ].sense == minimize)

      self._converger.reset()

      # construct the instances for each scenario.
      #
      # garbage collection noticeably slows down PH when dealing with
      # large numbers of scenarios. disable prior to instance construction,
      # and then re-enable. there isn't much collection to do as instances
      # are constructed.
      re_enable_gc = gc.isenabled()
      gc.disable()

      if self._verbose is True:
         if self._scenario_tree._scenario_based_data == 1:
            print "Scenario-based instance initialization enabled"
         else:
            print "Node-based instance initialization enabled"

      for scenario in self._scenario_tree._scenarios:

         scenario_instance = construct_scenario_instance(self._scenario_tree,
                                                         self._scenario_data_directory_name,
                                                         scenario._name,
                                                         self._model,
                                                         self._verbose,
                                                         preprocess=False,
                                                         linearize=linearize,
                                                         simplify=self._simplify_expressions)

         # IMPT: disable canonical representation construction for ASL solvers.
         #       this is a hack, in that we need to address encodings and
         #       the like at a more general level.
         if self._solver_type == "asl":
            scenario_instance.skip_canonical_repn = True
         else:
            scenario_instance.preprocess()

         self._instances[scenario._name] = scenario_instance
         self._instances[scenario._name].name = scenario._name
         self._instance_augmented_attributes[scenario._name] = []

      # perform a single pass of garbage collection and re-enable automatic collection.
      if re_enable_gc is True:
         gc.collect()
         gc.enable()

      # with the scenario instances now available, have the scenario tree compute the
      # variable match indices at each node.
      self._scenario_tree.defineVariableIndexSets(self._instances)

      # let plugins know if they care - this callback point allows
      # users to create/modify the original scenario instances and/or
      # the scenario tree prior to creating PH-related parameters,
      # variables, and the like.
      for plugin in self._ph_plugins:
         plugin.post_instance_creation(self)

      # create ph-specific parameters (weights, xbar, etc.) for each instance.

      if self._verbose is True:
         print "Creating weight, average, and rho parameter vectors for scenario instances"

      self._create_ph_scenario_parameters()

      # if specified, run the user script to initialize variable rhos at their whim.
      if self._rho_setter is not None:
         print "Executing user rho set script from filename="+self._rho_setter
         execfile(self._rho_setter)

      # with the instances created, run the user script to initialize variable bounds.
      if self._bounds_setter is not None:
         print "Executing user variable bounds set script from filename=", self._bounds_setter
         execfile(self._bounds_setter)

      # create parameters to store variable statistics (of general utility) at each node in the scenario tree.

      if self._verbose is True:
         print "Creating variable statistic (min/avg/max) parameter vectors for scenario tree nodes"

      # do this for all stages, simply for completeness, i.e., to create a fully populated scenario tree.
      for stage in self._scenario_tree._stages:

         # first, gather all unique variables referenced in this stage
         # this "gather" step is currently required because we're being lazy
         # in terms of index management in the scenario tree - which
         # should really be done in terms of sets of indices.

         stage_variables = {}
         for (reference_variable, index_template) in stage._variables:
            if reference_variable.name not in stage_variables.keys():
               stage_variables[reference_variable.name] = reference_variable

         # next, create min/avg/max parameters for each variable in the corresponding tree node.
         # NOTE: the parameter names below could really be empty, as they are never referenced
         #       explicitly.
         for (variable_name, reference_variable) in stage_variables.items():
            
            for tree_node in stage._tree_nodes:

               new_variable_index = getattr(self._instances[tree_node._scenarios[0]._name], variable_name)._index

               new_min_index = new_variable_index
               new_min_parameter_name = "NODEMIN_"+reference_variable.name
               # this bit of ugliness is due to Pyomo not correctly handling the Param construction
               # case when the supplied index set consists strictly of None, i.e., the source variable
               # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
               new_min_parameter = None
               if (len(new_min_index) is 1) and (None in new_min_index):
                  new_min_parameter = Param(name=new_min_parameter_name, mutable=True, nochecking=True)
               else:
                  new_min_parameter = Param(new_min_index, name=new_min_parameter_name, mutable=True, nochecking=True)
               for index in new_min_index:
                  new_min_parameter[index] = 0.0
               tree_node._minimums[reference_variable.name] = new_min_parameter

               new_avg_index = new_variable_index
               new_avg_parameter_name = "NODEAVG_"+reference_variable.name
               new_avg_parameter = None
               if (len(new_avg_index) is 1) and (None in new_avg_index):
                  new_avg_parameter = Param(name=new_avg_parameter_name, mutable=True, nochecking=True)
               else:
                  new_avg_parameter = Param(new_avg_index, name=new_avg_parameter_name, mutable=True, nochecking=True)
               for index in new_avg_index:
                  new_avg_parameter[index] = 0.0
               tree_node._averages[reference_variable.name] = new_avg_parameter

               new_max_index = new_variable_index
               new_max_parameter_name = "NODEMAX_"+reference_variable.name
               new_max_parameter = None
               if (len(new_max_index) is 1) and (None in new_max_index):
                  new_max_parameter = Param(name=new_max_parameter_name, mutable=True, nochecking=True)
               else:
                  new_max_parameter = Param(new_max_index, name=new_max_parameter_name, mutable=True, nochecking=True)
               for index in new_max_index:
                  new_max_parameter[index] = 0.0
               tree_node._maximums[reference_variable.name] = new_max_parameter

      # the objective functions are modified throughout the course of PH iterations.
      # save the original, as a baseline to modify in subsequent iterations. reserve
      # the original objectives, for subsequent modification.
      self._original_objective_expression = {}
      for instance_name, instance in self._instances.items():
         objective_name = instance.active_components(Objective).keys()[0]
         expr = instance.active_components(Objective)[objective_name]._data[None].expr
         if isinstance(expr, Expression) is False:
            expr = _IdentityExpression(expr)
         self._original_objective_expression[instance_name] = expr

      # cache the number of discrete and continuous variables in the master instance. this value
      # is of general use, e.g., in the converger classes and in plugins.
      (self._total_discrete_vars,self._total_continuous_vars) = self.compute_variable_counts()
      if self._verbose is True:
         print "Total number of discrete instance variables="+str(self._total_discrete_vars)
         print "Total number of continuous instance variables="+str(self._total_continuous_vars)

      # track the total number of fixed variables of each category at the end of each PH iteration.
      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()

      # indicate that we're ready to run.
      self._initialized = True

      if self._verbose is True:
         print "PH successfully created model instances for all scenarios"

      self._init_end_time = time.time()

      if self._verbose is True:
         print "PH is successfully initialized"
         if self._output_times is True:
            print "Initialization time=%8.2f seconds" % (self._init_end_time - self._init_start_time)

      # let plugins know if they care.
      if self._verbose is True:
         print "Initializing PH plugins"
      for plugin in self._ph_plugins:
         plugin.post_ph_initialization(self)
      if self._verbose is True:
         print "PH plugin initialization complete"

   """ Perform the non-weighted scenario solves and form the initial w and xbars.
   """
   def iteration_0_solve(self):

      if self._verbose is True:
         print "------------------------------------------------"
         print "Starting PH iteration 0 solves"

      self._current_iteration = 0

      solve_start_time = time.time()

      # STEP 0: set up all global solver options.
      self._solver.mipgap = self._mipgap

      # STEP 1: queue up the solves for all scenario sub-problems and
      #         grab all the action handles for the subsequent barrier sync.

      action_handles = []
      scenario_action_handle_map = {} # maps scenario names to action handles
      action_handle_scenario_map = {} # maps action handles to scenario names

      for scenario in self._scenario_tree._scenarios:

         instance = self._instances[scenario._name]

         if self._verbose is True:
            print "Queuing solve for scenario=" + scenario._name

         # IMPT: You have to re-presolve if approximating continuous variable penalty terms with a
         #       piecewise linear function. otherwise, the newly introduced variables won't be flagged
         #       as unused (as is correct for iteration 0), and the output writer will crater.
         # IMPT: I decided to presolve unconditionally, as PH extensions can add arbitrary components
         #       to the base scenario instances - and the variable values/etc. need to be collectged.
         instance.preprocess()

         # there's nothing to warm-start from in iteration 0, so don't include the keyword in the solve call.
         # the reason you don't want to include it is that some solvers don't know how to handle the keyword
         # at all (despite it being false). you might want to solve iteration 0 solves using some other solver.

         new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
         scenario_action_handle_map[scenario._name] = new_action_handle
         action_handle_scenario_map[new_action_handle] = scenario._name

         action_handles.append(new_action_handle)

      # STEP 2: loop for the solver results, reading them and loading
      #         them into instances as they are available.

      if self._verbose is True:
         print "Waiting for scenario sub-problem solves"

      num_results_so_far = 0

      while (num_results_so_far < len(self._scenario_tree._scenarios)):

         action_handle = self._solver_manager.wait_any()
         results = self._solver_manager.get_results(action_handle)
         scenario_name = action_handle_scenario_map[action_handle]
         instance = self._instances[scenario_name]

         if self._verbose is True:
            print "Results obtained for scenario="+scenario_name

         if len(results.solution) == 0:
            results.write(num=1)
            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"

         if self._output_solver_results is True:
            print "Results for scenario=",scenario_name
            results.write(num=1)

         start_time = time.time()
         instance.load(results)
         end_time = time.time()
         if self._output_times is True:
            print "Time loading results into instance="+str(end_time-start_time)+" seconds"

         if self._verbose is True:
            print "Successfully loaded solution for scenario="+scenario_name

         num_results_so_far = num_results_so_far + 1

      if self._verbose is True:
         print "Scenario sub-problem solves completed"

      solve_end_time = time.time()
      self._cumulative_solve_time += (solve_end_time - solve_start_time)

      if self._output_times is True:
         print "Aggregate sub-problem solve time this iteration=%8.2f" % (solve_end_time - solve_start_time)

      if self._verbose is True:
         print "Successfully completed PH iteration 0 solves - solution statistics:"
         print "         Scenario              Objective                  Value"
         for scenario in self._scenario_tree._scenarios:
            instance = self._instances[scenario._name]
            for objective_name in instance.active_components(Objective):
               objective = instance.active_components(Objective)[objective_name]
               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
         print "------------------------------------------------"

   #
   # recompute the averages, minimum, and maximum statistics for all variables to be blended by PH, i.e.,
   # not appearing in the final stage. technically speaking, the min/max aren't required by PH, but they
   # are used often enough to warrant their computation and it's basically free if you're computing the
   # average.
   #
   def update_variable_statistics(self):

      start_time = time.time()

      # NOTE: the following code has some optimizations that are not normally recommended, in particular
      #       the direct access and manipulation of parameters via the .value attribute instead of the 
      #       user-level-preferred value() method. this is justifiable in this particular instance
      #       because we are creating the PH parameters (and therefore can manipulate them safely), and
      #       this routine takes a non-trivial amount of the overall run-time.

      # compute statistics over all stages, even the last. this is necessary in order to
      # successfully snapshot a scenario tree solution from the average values.
      for stage in self._scenario_tree._stages:

         for tree_node in stage._tree_nodes:

            for (variable, index_template) in stage._variables:

               variable_name = variable.name

               variable_indices = tree_node._variable_indices[variable_name]

               avg_parameter_name = "PHAVG_"+variable_name

               tree_node_var_mins = tree_node._minimums[variable_name]
               tree_node_var_avgs = tree_node._averages[variable_name]
               tree_node_var_maxs = tree_node._maximums[variable_name]

               scenario_variables = []
               for scenario in tree_node._scenarios:
                  instance = self._instances[scenario._name]
                  scenario_variables.append(getattr(instance, variable_name))

               for index in variable_indices:

                  min = float("inf")
                  avg = 0.0
                  max = float("-inf")
                  node_probability = 0.0

                  is_used = True # until proven otherwise
                  for scenario_variable in scenario_variables:

                     if scenario_variable[index].status == VarStatus.unused:
                        is_used = False
                     else:
                        node_probability += scenario._probability
                        var_value = scenario_variable[index].value
                        if var_value < min:
                           min = var_value
                        avg += (scenario._probability * var_value)
                        if var_value > max:
                           max = var_value

                  if is_used is True:

                     tree_node_var_mins[index].value = min
                     tree_node_var_avgs[index].value = avg / node_probability
                     tree_node_var_maxs[index].value = max

                     # distribute the newly computed average to the xbar variable in
                     # each instance/scenario associated with this node. only do this
                     # if the variable is used!
                     for scenario in tree_node._scenarios:
                        instance = self._instances[scenario._name]
                        try:
                           avg_parameter = getattr(instance, avg_parameter_name)
                           avg_parameter[index].value = avg / node_probability
                        except:
                           pass

      end_time = time.time()
      self._cumulative_xbar_time += (end_time - start_time)

   def update_weights(self):

      # because the weight updates rely on the xbars, and the xbars are node-based,
      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
      start_time = time.time()

      # NOTE: the following code has some optimizations that are not normally recommended, in particular
      #       the direct access and manipulation of parameters via the .value attribute instead of the 
      #       user-level-preferred value() method. this is justifiable in this particular instance
      #       because we are creating the PH parameters (and therefore can manipulate them safely), and
      #       this routine takes a non-trivial amount of the overall run-time.

      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.

         for tree_node in stage._tree_nodes:

            for (variable, index_template) in stage._variables:

               variable_name = variable.name
               blend_parameter_name = "PHBLEND_"+variable_name
               weight_parameter_name = "PHWEIGHT_"+variable_name
               rho_parameter_name = "PHRHO_"+variable_name

               variable_indices = tree_node._variable_indices[variable_name]

               for index in variable_indices:

                  tree_node_average = tree_node._averages[variable_name][index].value

                  for scenario in tree_node._scenarios:

                     instance = self._instances[scenario._name]

                     if getattr(instance,variable_name)[index].status != VarStatus.unused:

                        weight_parameter = getattr(instance, weight_parameter_name)

                        # we are currently not updating weights if blending is disabled for a variable.
                        # this is done on the premise that unless you are actively trying to move
                        # the variable toward the mean, the weights will blow up and be huge by the
                        # time that blending is activated.
                        variable_blend_indicator = getattr(instance, blend_parameter_name)[index].value

                        # get the weight and rho parameters for this variable/index combination.
                        rho_value = getattr(instance, rho_parameter_name)[index].value
                        current_variable_weight = weight_parameter[index].value

                        # if I'm maximizing, invert value prior to adding (hack to implement negatives).
                        # probably fixed in Pyomo at this point - I just haven't checked in a long while.
                        if self._is_minimizing is False:
                           current_variable_weight = (-current_variable_weight)
                        current_variable_value = getattr(instance, variable_name)[index]()
                        new_variable_weight = current_variable_weight + variable_blend_indicator * rho_value * (current_variable_value - tree_node_average)
                        # I have the correct updated value, so now invert if maximizing.
                        if self._is_minimizing is False:
                           new_variable_weight = (-new_variable_weight)
                        weight_parameter[index].value = new_variable_weight

      # we shouldn't have to re-simplify the expression, as we aren't adding any constant-variable terms - just modifying parameters.

      end_time = time.time()
      self._cumulative_weight_time += (end_time - start_time)

   def form_iteration_k_objectives(self):

      for instance_name, instance in self._instances.items():

         new_attrs = form_ph_objective(instance_name, \
                                       instance, \
                                       self._original_objective_expression[instance_name], \
                                       self._scenario_tree, \
                                       self._linearize_nonbinary_penalty_terms, \
                                       self._drop_proximal_terms, \
                                       self._retain_quadratic_binary_terms, \
                                       self._breakpoint_strategy, \
                                       self._integer_tolerance, \
                                       self._simplify_expressions)
         self._instance_augmented_attributes[instance_name].extend(new_attrs)

   def iteration_k_solve(self):

      if self._verbose is True:
         print "------------------------------------------------"
         print "Starting PH iteration " + str(self._current_iteration) + " solves"

      # cache the objective values generated by PH for output at the end of this function.
      ph_objective_values = {}

      solve_start_time = time.time()

      # STEP -1: if using a PH solver manager, propagate current weights/averages to the appropriate solver servers.
      #          ditto the tree node statistics, which are necessary if linearizing (so an optimization could be
      #          performed here).
      # NOTE: We aren't currently propagating rhos, as they generally don't change - we need to
      #       have a flag, though, indicating whether the rhos have changed, so they can be
      #       transmitted if needed.
      if self._solver_manager_type == "ph":
         self._transmit_weights_and_averages()
         self._transmit_tree_node_statistics()

      # STEP 0: set up all global solver options.
      self._solver.mipgap = self._mipgap

      # STEP 1: queue up the solves for all scenario sub-problems and
      #         grab all of the action handles for the subsequent barrier sync.

      action_handles = []
      scenario_action_handle_map = {} # maps scenario names to action handles
      action_handle_scenario_map = {} # maps action handles to scenario names

      for scenario in self._scenario_tree._scenarios:

         instance = self._instances[scenario._name]

         if self._verbose is True:
            print "Queuing solve for scenario=" + scenario._name

         # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
         # don't do this, you won't see any chance in the output files as you vary the problem parameters!
         # ditto for instance fixing!
         instance.preprocess()

         # once past iteration 0, there is always a feasible solution from which to warm-start.
         # however, you might want to disable warm-start when the solver is behaving badly (which does happen).
         new_action_handle = None
         if (self._disable_warmstarts is False) and (self._solver.warm_start_capable() is True):
            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
         else:
            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)

         scenario_action_handle_map[scenario._name] = new_action_handle
         action_handle_scenario_map[new_action_handle] = scenario._name

         action_handles.append(new_action_handle)

      # STEP 2: loop for the solver results, reading them and loading
      #         them into instances as they are available.
      if self._verbose is True:
         print "Waiting for scenario sub-problem solves"

      num_results_so_far = 0

      while (num_results_so_far < len(self._scenario_tree._scenarios)):

         action_handle = self._solver_manager.wait_any()
         results = self._solver_manager.get_results(action_handle)
         scenario_name = action_handle_scenario_map[action_handle]
         instance = self._instances[scenario_name]

         if self._verbose is True:
            print "Results obtained for scenario="+scenario_name

         if len(results.solution) == 0:
            results.write(num=1)
            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"

         if self._output_solver_results is True:
            print "Results for scenario=",scenario_name
            results.write(num=1)

         start_time = time.time()
         instance.load(results)
         end_time = time.time()
         if self._output_times is True:
            print "Time loading results into instance="+str(end_time-start_time)+" seconds"

         if self._verbose is True:
            print "Successfully loaded solution for scenario="+scenario_name

         # we're assuming there is a single solution.
         # the "value" attribute is a pre-defined feature of any solution - it is relative to whatever
         # objective was selected during optimization, which of course should be the PH objective.
         try:
            ph_objective_values[instance.name] = float(results.solution(0).objective['f'].value)
         except AttributeError:
            # some solvers (e.g., through the SOL interface) don't report objective function values.
            ph_objective_values[instance.name] = 0.0

         num_results_so_far = num_results_so_far + 1

      if self._verbose is True:
         print "Scenario sub-problem solves completed"

      solve_end_time = time.time()
      self._cumulative_solve_time += (solve_end_time - solve_start_time)

      if self._output_times is True:
         print "Aggregate sub-problem solve time this iteration=%8.2f" % (solve_end_time - solve_start_time)

      if self._verbose is True:
         print "Successfully completed PH iteration " + str(self._current_iteration) + " solves - solution statistics:"
         print "  Scenario             PH Objective             Cost Objective"
         for scenario in self._scenario_tree._scenarios:
            instance = self._instances[scenario._name]
            for objective_name in instance.active_components(Objective):
               objective = instance.active_components(Objective)[objective_name]
               print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], self._scenario_tree.compute_scenario_cost(instance))

   def solve(self):

      self._solve_start_time = time.time()
      self._cumulative_solve_time = 0.0
      self._cumulative_xbar_time = 0.0
      self._cumulative_weight_time = 0.0

      print "Starting PH"

      if self._initialized == False:
         raise RuntimeError, "PH is not initialized - cannot invoke solve() method"

      # garbage collection noticeably slows down PH when dealing with
      # large numbers of scenarios. fortunately, there are well-defined
      # points at which garbage collection makes sense (and there isn't a
      # lot of collection to do). namely, after each PH iteration.
      re_enable_gc = gc.isenabled()
      gc.disable()

      print "Initiating PH iteration=" + `self._current_iteration`

      self.iteration_0_solve()

      # update variable statistics prior to any output.
      self.update_variable_statistics()

      if (self._verbose is True) or (self._report_solutions is True):
         print "Variable values following scenario solves:"
         self.pprint(False, False, True, False, output_only_statistics=self._report_only_statistics)

      # let plugins know if they care.
      for plugin in self._ph_plugins:
         plugin.post_iteration_0_solves(self)

      # update the fixed variable statistics.
      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()

      if self._verbose is True:
         print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
         print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"

      # always output the convergence metric and first-stage cost statistics, to give a sense of progress.
      self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
      first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()
      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)

      self.update_weights()

      # let plugins know if they care.
      for plugin in self._ph_plugins:
         plugin.post_iteration_0(self)

      # if using a PH solver server, trasnsmit the rhos prior to the iteration
      # k solve sequence. for now, we are assuming that the rhos don't change
      # on a per-iteration basis, but that assumption can be easily relaxed.
      # it is important to do this after the plugins have a chance to do their
      # computation.
      if self._solver_manager_type == "ph":
         self._transmit_rhos()
         self._enable_ph_objectives()

      # checkpoint if it's time - which it always is after iteration 0,
      # if the interval is >= 1!
      if (self._checkpoint_interval > 0):
         self.checkpoint(0)

      # garbage-collect if it wasn't disabled entirely.
      if re_enable_gc is True:
         gc.collect()

      # there is an upper bound on the number of iterations to execute -
      # the actual bound depends on the converger supplied by the user.
      for i in range(1, self._max_iterations+1):

         self._current_iteration = self._current_iteration + 1

         print "Initiating PH iteration=" + `self._current_iteration`

         if (self._verbose is True) or (self._report_weights is True):
            print "Variable averages and weights prior to scenario solves:"
            self.pprint(True, True, False, False, output_only_statistics=self._report_only_statistics)

         # with the introduction of piecewise linearization, the form of the
         # penalty-weighted objective is no longer fixed. thus, we need to
         # create the objectives each PH iteration. we could be fancy in
         # the case of quadratic objectives (eliminating redundancy), but
         # if we decide to do this, see the "IMPT" note below.
         self.form_iteration_k_objectives()

         # IMPT: for efficiency purposes, the PH objective formation calls above create expression trees by
         #       evaluating the parameters involved in either quadratic or linear expressions directly, e.g.,
         #       via value(). as a consequence, if the user modifies PH weights or rhos in the plugin 
         #       callback below, it will *not* be propagated to the objective expression.

         # let plugins know if they care.
         for plugin in self._ph_plugins:
            plugin.pre_iteration_k_solves(self)

         # do the actual solves.
         self.iteration_k_solve()

         # update variable statistics prior to any output.
         self.update_variable_statistics()

         if (self._verbose is True) or (self._report_solutions is True):
            print "Variable values following scenario solves:"
            self.pprint(False, False, True, False, output_only_statistics=self._report_only_statistics)

         # we don't technically have to do this at the last iteration,
         # but with checkpointing and re-starts, you're never sure
         # when you're executing the last iteration.
         self.update_weights()

         # let plugins know if they care.
         for plugin in self._ph_plugins:
            plugin.post_iteration_k_solves(self)

         # update the fixed variable statistics.
         (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()

         if self._verbose is True:
            print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
            print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"

         # let plugins know if they care.
         for plugin in self._ph_plugins:
            plugin.post_iteration_k(self)

         # at this point, all the real work of an iteration is complete.

         # checkpoint if it's time.
         if (self._checkpoint_interval > 0) and (i % self._checkpoint_interval is 0):
            self.checkpoint(i)

         # check for early termination.
         self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
         first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()
         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)

         if self._converger.isConverged(self) is True:
            if self._total_discrete_vars == 0:
               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)
            else:
               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)+" or all discrete variables are fixed"
            break

         # if we're terminating due to exceeding the maximum iteration count, print a message
         # indicating so - otherwise, you get a quiet, information-free output trace.
         if i == self._max_iterations:
            print "Halting PH - reached maximal iteration count="+str(self._max_iterations)

         # garbage-collect if it wasn't disabled entirely.
         if re_enable_gc is True:
            gc.collect()

      # re-enable the normal garbage collection mode.
      if re_enable_gc is True:
         gc.enable()

      if self._verbose is True:
         print "Number of discrete variables fixed before final plugin calls="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
         print "Number of continuous variables fixed before final plugin calls="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"

      # let plugins know if they care. do this before
      # the final solution / statistics output, as the plugins
      # might do some final tweaking.
      for plugin in self._ph_plugins:
         plugin.post_ph_execution(self)

      # update the fixed variable statistics - the plugins might have done something.
      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()

      self._solve_end_time = time.time()

      print "PH complete"

      print "Convergence history:"
      self._converger.pprint()

      print "Final number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
      print "Final number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"

      # populate the scenario tree solution from the instances - to ensure consistent state
      # across the scenario tree instance and the scenario instances.
      self._scenario_tree.snapshotSolutionFromInstances(self._instances)      

      print "Final variable values:"
      self.pprint(False, False, True, True, output_only_statistics=self._report_only_statistics)

      print "Final costs:"
      self._scenario_tree.pprintCosts(self._instances)

      if self._output_scenario_tree_solution is True:
         self._scenario_tree.snapshotSolutionFromAverages(self._instances)
         print "Final solution (scenario tree format):"
         self._scenario_tree.pprintSolution()

      if (self._verbose is True) and (self._output_times is True):
         print "Overall run-time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)

      # cleanup the scenario instances for post-processing - ideally, we want to leave them in
      # their original state, minus all the PH-specific stuff. we don't do all cleanup (leaving
      # things like rhos, etc), but we do clean up constraints, as that really hoses up the ef writer.
      self._cleanup_scenario_instances()

   #
   # prints a summary of all collected time statistics
   #
   def print_time_stats(self):

      print "PH run-time statistics (user):"

      print "Initialization time=  %8.2f seconds" % (self._init_end_time - self._init_start_time)
      print "Overall solve time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
      print "Scenario solve time=  %8.2f seconds" % self._cumulative_solve_time
      print "Average update time=  %8.2f seconds" % self._cumulative_xbar_time
      print "Weight update time=   %8.2f seconds" % self._cumulative_weight_time

   #
   # a utility to determine whether to output weight / average / etc. information for
   # a variable/node combination. when the printing is moved into a callback/plugin,
   # this routine will go there. for now, we don't dive down into the node resolution -
   # just the variable/stage.
   #
   def should_print(self, stage, variable):

      if self._output_continuous_variable_stats is False:

         variable_type = variable.domain

         if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):

            return False

      return True

   #
   # pretty-prints the state of the current variable averages, weights, and values.
   # inputs are booleans indicating which components should be output.
   #
   def pprint(self, output_averages, output_weights, output_values, output_fixed, output_only_statistics=False):

      if self._initialized is False:
         raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"

      # print tree nodes and associated variable/xbar/ph information in stage-order
      # we don't blend in the last stage, so we don't current care about printing the associated information.
      for stage in self._scenario_tree._stages[:-1]:

         print "\tStage=" + str(stage._name)

         num_outputs_this_stage = 0 # tracks the number of outputs on a per-index basis.

         for (variable, index_template) in stage._variables:

            variable_name = variable.name

            if self.should_print(stage, variable) is True:

               num_outputs_this_variable = 0 # track, so we don't output the variable names unless there is an entry to report.

               for tree_node in stage._tree_nodes:

                  variable_indices = tree_node._variable_indices[variable_name]               

                  for index in variable_indices:

                     weight_parameter_name = "PHWEIGHT_"+variable_name

                     num_outputs_this_index = 0 # track, so we don't output the variable index more than once.

                     # determine if the variable/index pair is used across the set of scenarios (technically,
                     # it should be good enough to check one scenario). ditto for "fixed" status. fixed does
                     # imply unused (see note below), but we care about the fixed status when outputting
                     # final solutions.

                     is_used = True # should be consistent across scenarios, so one "unused" flags as invalid.
                     is_fixed = False

                     for scenario in tree_node._scenarios:
                        instance = self._instances[scenario._name]
                        variable_value = getattr(instance,variable_name)[index]
                        if variable_value.status == VarStatus.unused:
                           is_used = False
                        if variable_value.fixed is True:
                           is_fixed = True

                     # IMPT: this is far from obvious, but variables that are fixed will - because
                     #       presolve will identify them as constants and eliminate them from all
                     #       expressions - be flagged as "unused" and therefore not output.

                     if ((output_fixed is True) and (is_fixed is True)) or (is_used is True):

                           minimum_value = None
                           maximum_value = None

                           if index is None:
                              minimum_value = tree_node._minimums[variable_name]
                              maximum_value = tree_node._maximums[variable_name]
                           else:
                              minimum_value = value(tree_node._minimums[variable_name][index])
                              maximum_value = value(tree_node._maximums[variable_name][index])

                           # there really isn't a need to output variables whose
                           # values are equal to 0 across-the-board. and there is
                           # good reason not to, i.e., the volume of output.
                           if (fabs(minimum_value) > self._integer_tolerance) or \
                              (fabs(maximum_value) > self._integer_tolerance):

                              num_outputs_this_stage = num_outputs_this_stage + 1
                              num_outputs_this_variable = num_outputs_this_variable + 1
                              num_outputs_this_index = num_outputs_this_index + 1

                              if num_outputs_this_variable == 1:
                                 print "\t\tVariable=" + variable_name

                              if num_outputs_this_index == 1:
                                 if index is not None:
                                    print "\t\t\tIndex:", indexToString(index),

                              if len(stage._tree_nodes) > 1:
                                 print ""
                                 print "\t\t\t\tTree Node="+tree_node._name,
                              if output_only_statistics is False:
                                 print "\t\t (Scenarios: ",
                                 for scenario in tree_node._scenarios:
                                    print scenario._name," ",
                                    if scenario == tree_node._scenarios[-1]:
                                       print ")"

                              if output_values is True:
                                 average_value = value(tree_node._averages[variable_name][index])
                                 if output_only_statistics is False:
                                    print "\t\t\t\tValues: ",
                                 for scenario in tree_node._scenarios:
                                    instance = self._instances[scenario._name]
                                    this_value = getattr(instance,variable_name)[index].value
                                    if output_only_statistics is False:
                                       print "%12.4f" % this_value,
                                    if scenario == tree_node._scenarios[-1]:
                                       if output_only_statistics is True:
                                          # there technically isn't any good reason not to always report
                                          # the min and max; the only reason we're not doing this currently
                                          # is to avoid updating our regression test baseline output.
                                          print "    Min=%12.4f" % (minimum_value),
                                          print "    Avg=%12.4f" % (average_value),
                                          print "    Max=%12.4f" % (maximum_value),
                                       else:
                                          print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
                                          print "    Avg=%12.4f" % (average_value),
                                       print ""
                              if output_weights:
                                 print "\t\t\t\tWeights: ",
                                 for scenario in tree_node._scenarios:
                                    instance = self._instances[scenario._name]
                                    print "%12.4f" % value(getattr(instance,weight_parameter_name)[index]),
                                    if scenario == tree_node._scenarios[-1]:
                                       print ""

                              if output_averages:
                                 print "\t\t\t\tAverage: %12.4f" % (tree_node._averages[variable_name][index].value)

         if num_outputs_this_stage == 0:
            print "\t\tNo non-converged variables in stage"

         # cost variables aren't blended, so go through the gory computation of min/max/avg.
         # we currently always print these.
         cost_variable_name = stage._cost_variable[0].name
         cost_variable_index = stage._cost_variable[1]
         if cost_variable_index is None:
            print "\t\tCost Variable=" + cost_variable_name
         else:
            print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)
         for tree_node in stage._tree_nodes:
            print "\t\t\tTree Node=" + tree_node._name,
            if output_only_statistics is False:
               print "\t\t (Scenarios: ",
               for scenario in tree_node._scenarios:
                  print scenario._name," ",
                  if scenario == tree_node._scenarios[-1]:
                     print ")"
            maximum_value = 0.0
            minimum_value = 0.0
            sum_values = 0.0
            num_values = 0
            first_time = True
            if output_only_statistics is False:
               print "\t\t\tValues: ",
            else:
               print "\t\t\t",
            for scenario in tree_node._scenarios:
                instance = self._instances[scenario._name]
                this_value = getattr(instance,cost_variable_name)[cost_variable_index].value
                if output_only_statistics is False:
                   if this_value is not None:
                      print "%12.4f" % this_value,
                   else:
                      # this is a hack, in case the stage cost variables are not returned. ipopt
                      # does this occasionally, for example, if stage cost variables are constrained
                      # to a constant value (and consequently preprocessed out).
                      print "%12s" % "Not Rprted",
                if this_value is not None:
                   num_values += 1
                   sum_values += this_value
                   if first_time is True:
                      first_time = False
                      maximum_value = this_value
                      minimum_value = this_value
                   else:
                      if this_value > maximum_value:
                         maximum_value = this_value
                      if this_value < minimum_value:
                         minimum_value = this_value
                if scenario == tree_node._scenarios[-1]:
                   if num_values > 0:
                      if output_only_statistics is True:
                         print "    Min=%12.4f" % (minimum_value),
                         print "    Avg=%12.4f" % (sum_values/num_values),
                         print "    Max=%12.4f" % (maximum_value),
                      else:
                         print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
                         print "    Avg=%12.4f" % (sum_values/num_values),
                   print ""