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

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

Disabling garbage collection in PySP during instance construction.

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