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

Last change on this file since 1768 was 1768, checked in by wehart, 11 years ago

Rework of Coopr to use the new PyUtilib? package decomposition.

NOTE: to use Coopr with this update, we need to work with a new version of coopr_install.

File size: 59.9 KB
Line 
1#  _________________________________________________________________________
2#
3#  Coopr: A COmmon Optimization Python Repository
4#  Copyright (c) 2008 Sandia Corporation.
5#  This software is distributed under the BSD License.
6#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
7#  the U.S. Government retains certain rights in this software.
8#  For more information, see the Coopr README.txt file.
9#  _________________________________________________________________________
10
11import sys
12import types
13from coopr.pyomo import *
14import copy
15import os.path
16import traceback
17import copy
18from coopr.opt import SolverResults,SolverStatus
19from coopr.opt.base import SolverFactory
20from coopr.opt.parallel import SolverManagerFactory
21import time
22import types
23import pickle
24
25from scenariotree import *
26from phutils import *
27
28from pyutilib.plugin.core import ExtensionPoint
29
30from coopr.pysp.phextension import IPHExtension
31
32class ProgressiveHedging(object):
33
34   #
35   # a utility intended for folks who are brave enough to script rho setting in a python file.
36   #
37
38   def setRhoAllScenarios(self, variable_value, rho_expression):
39
40      variable_name = None
41      variable_index = None
42
43      if isVariableNameIndexed(variable_value.name) is True:
44
45         variable_name, variable_index = extractVariableNameAndIndex(variable_value.name)
46
47      else:
48
49         variable_name = variable_value.name
50         variable_index = None
51     
52      new_rho_value = rho_expression()
53
54      if self._verbose is True:
55         print "Setting rho="+str(new_rho_value)+" for variable="+variable_value.name
56
57      for instance_name, instance in self._instances.items():
58
59         rho_param = getattr(instance, "PHRHO_"+variable_name)
60         rho_param[variable_index] = new_rho_value
61
62   #
63   # a simple utility to count the number of continuous and discrete variables in a set of instances.
64   # unused variables are ignored, and counts include all active indices. returns a pair - num-discrete,
65   # num-continuous.
66   #
67
68   def compute_variable_counts(self):
69
70      num_continuous_vars = 0
71      num_discrete_vars = 0
72     
73      for stage in self._scenario_tree._stages:
74
75         if stage != self._scenario_tree._stages[-1]: # no blending over the final stage
76
77            for tree_node in stage._tree_nodes:
78
79               for (variable, index_template, variable_indices) in stage._variables:
80
81                  variable_name = variable.name
82
83                  variable_type = variable.domain
84
85                  for index in variable_indices:
86
87                     is_used = True # until proven otherwise                     
88                     for scenario in tree_node._scenarios:
89                        instance = self._instances[scenario._name]
90                        if getattr(instance,variable_name)[index].status == VarStatus.unused:
91                           is_used = False
92
93                     if is_used is True:                       
94
95                        # JPW TBD: ideally, we want an "is-discrete" check in the logic below, which COOPR doesn't currently support.
96                        if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
97                           num_discrete_vars = num_discrete_vars + 1
98                        else:
99                           num_continuous_vars = num_continuous_vars + 1
100
101      return (num_discrete_vars, num_continuous_vars)
102
103   #
104   # checkpoint the current PH state via pickle'ing. the input iteration count
105   # simply serves as a tag to create the output file name. everything with the
106   # exception of the _ph_plugin, _solver_manager, and _solver attributes are
107   # pickled. currently, plugins fail in the pickle process, which is fine as
108   # JPW doesn't think you want to pickle plugins (particularly the solver and
109   # solver manager) anyway. For example, you might want to change those later,
110   # after restoration - and the PH state is independent of how scenario
111   # sub-problems are solved.
112   #
113
114   def checkpoint(self, iteration_count):
115
116      checkpoint_filename = "checkpoint."+str(iteration_count)
117
118      tmp_ph_plugin = self._ph_plugin
119      tmp_solver_manager = self._solver_manager
120      tmp_solver = self._solver
121
122      self._ph_plugin = None
123      self._solver_manager = None
124      self._solver = None
125     
126      checkpoint_file = open(checkpoint_filename, "w")
127      pickle.dump(self,checkpoint_file)
128      checkpoint_file.close()
129
130      self._ph_plugin = tmp_ph_plugin
131      self._solver_manager = tmp_solver_manager
132      self._solver = tmp_solver
133
134      print "Checkpoint written to file="+checkpoint_filename
135   
136   #
137   # ditto above, but count the number of fixed discrete and continuous variables.
138   # important: once a variable (value) is fixed, it is flagged as unused in the
139   # course of presolve - because it is no longer referenced. this makes sense,
140   # of course; it's just something to watch for. this is an obvious assumption
141   # that we won't be fixing unused variables, which should not be an issue.
142   #
143
144   def compute_fixed_variable_counts(self):
145
146      num_fixed_continuous_vars = 0
147      num_fixed_discrete_vars = 0
148     
149      for stage in self._scenario_tree._stages:
150         
151         if stage != self._scenario_tree._stages[-1]: # no blending over the final stage
152           
153            for tree_node in stage._tree_nodes:
154
155               for (variable, index_template, variable_indices) in stage._variables:
156
157                  variable_name = variable.name
158
159                  variable_type = variable.domain
160
161                  for index in variable_indices:
162
163                     # implicit assumption is that if a variable value is fixed in one
164                     # scenario, it is fixed in all scenarios.
165
166                     is_fixed = False # until proven otherwise
167                     for scenario in tree_node._scenarios:
168                        instance = self._instances[scenario._name]
169                        var_value = getattr(instance,variable_name)[index]
170                        if var_value.fixed is True:
171                           is_fixed = True
172
173                     if is_fixed is True:
174
175                        # JPW TBD: ideally, we want an "is-discrete" check in the logic below, which COOPR doesn't currently support.
176                        if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
177                           num_fixed_discrete_vars = num_fixed_discrete_vars + 1
178                        else:
179                           num_fixed_continuous_vars = num_fixed_continuous_vars + 1                           
180
181      return (num_fixed_discrete_vars, num_fixed_continuous_vars)
182
183   """ Constructor
184       Arguments:
185          max_iterations        the maximum number of iterations to run PH (>= 0). defaults to 0.
186          rho                   the global rho value (> 0). defaults to 0.
187          rho_setter            an optional name of a python file used to set particular variable rho values.
188          solver                the solver type that PH uses to solve scenario sub-problems. defaults to "cplex".
189          solver_manager        the solver manager type that coordinates scenario sub-problem solves. defaults to "serial".
190          keep_solver_files     do I keep intermediate solver files around (for debugging)? defaults to False.
191          output_solver_log     do I dump the solver log (as it is being generated) to the screen? defaults to False.
192          output_solver_results do I output (for debugging) the detailed solver results, including solutions, for scenario solves? defaults to False.
193          verbose               does the PH object stream debug/status output? defaults to False.
194          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).
195          checkpoint_interval   how many iterations between writing a checkpoint file containing the entire PH state? defaults to 0, indicating never.
196
197   """
198   def __init__(self, *args, **kwds):
199
200      # PH configuration parameters
201      self._rho = 0.0 # a default, global values for rhos.
202      self._rho_setter = None # filename for the modeler to set rho on a per-variable or per-scenario basis.
203      self._max_iterations = 0
204
205      # PH reporting parameters
206      self._verbose = False # do I flood the screen with status output?
207      self._output_continuous_variable_stats = True # when in verbose mode, do I output weights/averages for continuous variables?
208      self._output_solver_results = False
209      self._output_times = False
210
211      # PH run-time variables
212      self._current_iteration = 0 # the 'k'
213      self._xbar = {} # current per-variable averages. maps (node_id, variable_name) -> value
214      self._initialized = False # am I ready to call "solve"? Set to True by the initialize() method.
215
216      # PH solver information / objects.
217      self._solver_type = "cplex"
218      self._solver_manager_type = "serial" # serial or pyro are the options currently available
219     
220      self._solver = None # will eventually be unnecessary once Bill eliminates the need for a solver in the solver manager constructor.
221      self._solver_manager = None 
222     
223      self._keep_solver_files = False
224      self._output_solver_log = False
225
226      # PH convergence computer/updater.
227      self._converger = None
228
229      # PH history
230      self._solutions = {}
231
232      # the checkpoint interval - expensive operation, but worth it for big models.
233      # 0 indicates don't checkpoint.
234      self._checkpoint_interval = 0
235
236      # all information related to the scenario tree (implicit and explicit).
237      self._model = None # not instantiated
238      self._model_instance = None # instantiated
239
240      self._scenario_tree = None
241     
242      self._scenario_data_directory = "" # this the prefix for all scenario data
243      self._instances = {} # maps scenario name to the corresponding model instance
244
245      # for various reasons (mainly hacks at this point), it's good to know whether we're minimizing or maximizing.
246      self._is_minimizing = None
247
248      # global handle to ph extension plugin
249      self._ph_plugin = ExtensionPoint(IPHExtension)
250
251      # PH timing statistics - relative to last invocation.
252      self._init_start_time = None # for initialization() method
253      self._init_end_time = None
254      self._solve_start_time = None # for solve() method
255      self._solve_end_time = None
256      self._cumulative_solve_time = None # seconds, over course of solve()
257      self._cumulative_xbar_time = None # seconds, over course of update_xbars()
258      self._cumulative_weight_time = None # seconds, over course of update_weights()
259
260      # do I disable warm-start for scenario sub-problem solves during PH iterations >= 1?
261      self._disable_warmstarts = False
262
263      # do I retain quadratic objective terms associated with binary variables? in general,
264      # there is no good reason to not linearize, but just in case, I introduced the option.
265      self._retain_quadratic_binary_terms = False
266
267      # PH default tolerances - for use in fixing and testing equality across scenarios,
268      # and other stuff.
269      self._integer_tolerance = 0.00001
270
271      # process the keyword options
272      for key in kwds.keys():
273         if key == "max_iterations":
274            self._max_iterations = kwds[key]
275         elif key == "rho":
276            self._rho = kwds[key]
277         elif key == "rho_setter":
278            self._rho_setter = kwds[key]           
279         elif key == "solver":
280            self._solver_type = kwds[key]
281         elif key == "solver_manager":
282            self._solver_manager_type = kwds[key]           
283         elif key == "keep_solver_files":
284            self._keep_solver_files = kwds[key]
285         elif key == "output_solver_results":
286            self._output_solver_results = kwds[key]
287         elif key == "output_solver_log":
288            self._output_solver_log = kwds[key]
289         elif key == "verbose":
290            self._verbose = kwds[key]
291         elif key == "output_times":
292            self._output_times = kwds[key]
293         elif key == "disable_warmstarts":
294            self._disable_warmstarts = kwds[key]
295         elif key == "retain_quadratic_binary_terms":
296            self._retain_quadratic_binary_terms = kwds[key]
297         elif key == "checkpoint_interval":
298            self._checkpoint_interval = kwds[key]           
299         else:
300            print "Unknown option=" + key + " specified in call to PH constructor"
301
302      # validate all "atomic" options (those that can be validated independently)
303      if self._max_iterations < 0:
304         raise ValueError, "Maximum number of PH iterations must be non-negative; value specified=" + `self._max_iterations`
305      if self._rho <= 0.0:
306         raise ValueError, "Value of rho paramter in PH must be non-zero positive; value specified=" + `self._rho`
307
308      # validate rho setter file if specified.
309      if self._rho_setter is not None:
310         if os.path.exists(self._rho_setter) is False:
311            raise ValueError, "The rho setter script file="+self._rho_setter+" does not exist"
312
313      # validate the checkpoint interval.
314      if self._checkpoint_interval < 0:
315         raise ValueError, "A negative checkpoint interval with value="+str(self._checkpoint_interval)+" was specified in call to PH constructor"
316
317      # construct the sub-problem solver.
318      self._solver = SolverFactory(self._solver_type)
319      if self._solver == None:
320         raise ValueError, "Unknown solver type=" + self._solver_type + " specified in call to PH constructor"
321      if self._keep_solver_files is True:
322         # TBD - this is a bit kludgy, as it requires PH to know/assume that
323         #       it is dealing with a system call solver. the solver factory should take keyword options as a fix.
324         self._solver.keepFiles = True
325
326      # construct the solver manager.
327      self._solver_manager = SolverManagerFactory(self._solver_manager_type)
328      if self._solver_manager is None:
329         raise ValueError, "Failed to create solver manager of type="+self._solver_manager_type+" specified in call to PH constructor"
330
331      # a set of all valid PH iteration indicies is generally useful for plug-ins, so create it here.
332      self._iteration_index_set = Set(name="PHIterations")
333      for i in range(0,self._max_iterations + 1):
334         self._iteration_index_set.add(i)     
335
336      # spit out parameterization if verbosity is enabled
337      if self._verbose is True:
338         print "PH solver configuration: "
339         print "   Max iterations=" + `self._max_iterations`
340         print "   Default global rho=" + `self._rho`
341         if self._rho_setter is not None:
342            print "   Rho initialization file=" + self._rho_setter
343         print "   Sub-problem solver type=" + `self._solver_type`
344         print "   Solver manager type=" + `self._solver_manager_type`
345         print "   Keep solver files? " + str(self._keep_solver_files)
346         print "   Output solver results? " + str(self._output_solver_results)
347         print "   Output solver log? " + str(self._output_solver_log)
348         print "   Output times? " + str(self._output_times)
349         print "   Checkpoint interval="+str(self._checkpoint_interval)
350
351   """ Initialize PH with model and scenario data, in preparation for solve().
352       Constructs and reads instances.
353   """
354   def initialize(self, scenario_data_directory_name=".", model=None, model_instance=None, scenario_tree=None, converger=None):
355
356      self._init_start_time = time.time()
357
358      if self._verbose is True:
359         print "Initializing PH"
360         print "   Scenario data directory=" + scenario_data_directory_name
361
362      if not os.path.exists(scenario_data_directory_name):
363         raise ValueError, "Scenario data directory=" + scenario_data_directory_name + " either does not exist or cannot be read"
364
365      self._scenario_data_directory_name = scenario_data_directory_name
366
367      # IMPT: The input model should be an *instance*, as it is very useful (critical!) to know
368      #       the dimensions of sets, be able to store suffixes on variable values, etc.
369      # TBD: There may be a way to see if a model is initialized - throw an exception if it is not!
370      if model is None:
371         raise ValueError, "A model must be supplied to the PH initialize() method"
372
373      if scenario_tree is None:
374         raise ValueError, "A scenario tree must be supplied to the PH initialize() method"
375
376      if converger is None:
377         raise ValueError, "A convergence computer must be supplied to the PH initialize() method"
378
379      self._model = model
380      self._model_instance = model_instance
381      self._scenario_tree = scenario_tree
382      self._converger = converger
383
384      model_objective = model._component[Objective]
385      self._is_minimizing = (model_objective[ model_objective.keys()[0] ].sense == minimize)
386
387      self._converger.reset()
388
389      # construct instances for each scenario
390      if self._verbose is True:
391         if self._scenario_tree._scenario_based_data == 1:
392            print "Scenario-based instance initialization enabled"
393         else:
394            print "Node-based instance initialization enabled"
395         
396      for scenario in self._scenario_tree._scenarios:
397
398         scenario_instance = None
399
400         if self._verbose is True:
401            print "Creating instance for scenario=" + scenario._name
402
403         try:
404            if self._scenario_tree._scenario_based_data == 1:
405               scenario_data_filename = self._scenario_data_directory_name + os.sep + scenario._name + ".dat"
406               if self._verbose is True:
407                  print "Data for scenario=" + scenario._name + " loads from file=" + scenario_data_filename
408               scenario_instance = (self._model).create(scenario_data_filename)
409            else:
410               scenario_instance = self._model.clone()
411               scenario_data = ModelData()
412               current_node = scenario._leaf_node
413               while current_node is not None:
414                  node_data_filename = self._scenario_data_directory_name + os.sep + current_node._name + ".dat"
415                  if self._verbose is True:
416                     print "Node data for scenario=" + scenario._name + " partially loading from file=" + node_data_filename
417                  scenario_data.add_data_file(node_data_filename)
418                  current_node = current_node._parent
419               scenario_data.read(model=scenario_instance)
420               scenario_instance._load_model_data(scenario_data)
421               scenario_instance.presolve()
422         except:
423            print "Encountered exception in model instance creation - traceback:"
424            traceback.print_exc()
425            raise RuntimeError, "Failed to create model instance for scenario=" + scenario._name
426
427         self._instances[scenario._name] = scenario_instance
428         self._instances[scenario._name].name = scenario._name
429
430      # TBD - Technically, we don't need to add these parameters to the models until after the iteration 0 solve (move later).
431
432      # create PH weight and xbar vectors, on a per-scenario basis, for each variable that is not in the
433      # final stage, i.e., for all variables that are being blended by PH. the parameters are created
434      # in the space of each scenario instance, so that they can be directly and automatically
435      # incorporated into the (appropriately modified) objective function.
436
437      if self._verbose is True:
438         print "Creating weight, average, and rho parameter vectors for scenario instances"
439
440      for (instance_name, instance) in self._instances.items():
441
442         # first, gather all unique variables referenced in any stage
443         # other than the last, independent of specific indices. this
444         # "gather" step is currently required because we're being lazy
445         # in terms of index management in the scenario tree - which
446         # should really be done in terms of sets of indices.
447         # NOTE: technically, the "instance variables" aren't really references
448         # to the variable in the instance - instead, the reference model. this
449         # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
450         instance_variables = {}
451         
452         for stage in self._scenario_tree._stages:
453           
454            if stage != self._scenario_tree._stages[-1]:
455               
456               for (reference_variable, index_template, reference_indices) in stage._variables:
457                 
458                  if reference_variable.name not in instance_variables.keys():
459                     
460                     instance_variables[reference_variable.name] = reference_variable
461
462         # for each blended variable, create a corresponding ph weight and average parameter in the instance.
463         # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor
464         # in the grand scheme of things.
465
466         for (variable_name, reference_variable) in instance_variables.items():
467
468            # PH WEIGHT
469
470            new_w_index = reference_variable._index # TBD - need to be careful with the shallow copy here
471            new_w_parameter_name = "PHWEIGHT_"+reference_variable.name
472            new_w_parameter = Param(new_w_index,name=new_w_parameter_name)
473            setattr(instance,new_w_parameter_name,new_w_parameter)
474
475            # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference
476            # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not
477            # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the
478            # issue in any case for now.
479            for index in new_w_index:
480               new_w_parameter[index] = 0.0
481
482            # PH AVG
483
484            new_avg_index = reference_variable._index # TBD - need to be careful with the shallow copy here
485            new_avg_parameter_name = "PHAVG_"+reference_variable.name
486            new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
487            setattr(instance,new_avg_parameter_name,new_avg_parameter)
488
489            for index in new_avg_index:
490               new_avg_parameter[index] = 0.0
491
492            # PH RHO
493
494            new_rho_index = reference_variable._index # TBD - need to be careful with the shallow copy here
495            new_rho_parameter_name = "PHRHO_"+reference_variable.name
496            new_rho_parameter = Param(new_rho_index,name=new_rho_parameter_name)
497            setattr(instance,new_rho_parameter_name,new_rho_parameter)
498
499            for index in new_rho_index:
500               new_rho_parameter[index] = self._rho
501
502            # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY.
503
504            # also controls whether weight updates proceed at any iteration.
505
506            new_blend_index = reference_variable._index # TBD - need to be careful with the shallow copy here
507            new_blend_parameter_name = "PHBLEND_"+reference_variable.name
508            new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary)
509            setattr(instance,new_blend_parameter_name,new_blend_parameter)
510
511            for index in new_blend_index:
512               new_blend_parameter[index] = 1
513           
514      # if specified, run the user script to initialize variable rhos at their whim.
515      if self._rho_setter is not None:
516
517         print "Executing user rho set script from filename=", self._rho_setter
518         execfile(self._rho_setter)
519
520      # create parameters to store variable statistics (of general utility) at each node in the scenario tree.
521
522      if self._verbose is True:
523         print "Creating variable statistic (min/avg/max) parameter vectors for scenario tree nodes"
524
525      for stage in self._scenario_tree._stages:
526         if stage != self._scenario_tree._stages[-1]:
527
528            # first, gather all unique variables referenced in this stage
529            # this "gather" step is currently required because we're being lazy
530            # in terms of index management in the scenario tree - which
531            # should really be done in terms of sets of indices.
532            stage_variables = {}
533            for (reference_variable, index_template, reference_index) in stage._variables:
534               if reference_variable.name not in stage_variables.keys():
535                  stage_variables[reference_variable.name] = reference_variable
536
537            # next, create min/avg/max parameters for each variable in the corresponding tree node.
538            # NOTE: the parameter names below could really be empty, as they are never referenced
539            #       explicitly.
540            for (variable_name, reference_variable) in stage_variables.items():
541               for tree_node in stage._tree_nodes:
542
543                  new_min_index = reference_variable._index # TBD - need to be careful with the shallow copy here (and below)
544                  new_min_parameter_name = "NODEMIN_"+reference_variable.name
545                  new_min_parameter = Param(new_min_index,name=new_min_parameter_name)
546                  for index in new_min_index:
547                     new_min_parameter[index] = 0.0
548                  tree_node._minimums[reference_variable.name] = new_min_parameter                                   
549
550                  new_avg_index = reference_variable._index
551                  new_avg_parameter_name = "NODEAVG_"+reference_variable.name
552                  new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
553                  for index in new_avg_index:
554                     new_avg_parameter[index] = 0.0
555                  tree_node._averages[reference_variable.name] = new_avg_parameter
556
557                  new_max_index = reference_variable._index
558                  new_max_parameter_name = "NODEMAX_"+reference_variable.name
559                  new_max_parameter = Param(new_max_index,name=new_max_parameter_name)
560                  for index in new_max_index:
561                     new_max_parameter[index] = 0.0
562                  tree_node._maximums[reference_variable.name] = new_max_parameter                                                     
563
564
565      # indicate that we're ready to run.
566      self._initialized = True
567
568      # cache the number of discrete and continuous variables in the master instance. this value
569      # is of general use, e.g., in the converger classes and in plugins.
570      (self._total_discrete_vars,self._total_continuous_vars) = self.compute_variable_counts()
571      if self._verbose is True:
572         print "Total number of discrete instance variables="+str(self._total_discrete_vars)
573         print "Total number of continuous instance variables="+str(self._total_continuous_vars)
574
575      # track the total number of fixed variables of each category at the end of each PH iteration.
576      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()     
577
578      if self._verbose is True:
579         print "PH successfully created model instances for all scenarios"
580
581      self._init_end_time = time.time()
582
583      if self._verbose is True:
584         print "PH is successfully initialized"
585         if self._output_times is True:
586            print "Initialization time=%8.2f seconds" % (self._init_end_time - self._init_start_time)
587
588      # let plugins know if they care.
589      if len(self._ph_plugin) == 1:
590         self._ph_plugin.service().post_ph_initialization(self)
591
592
593   """ Perform the non-weighted scenario solves and form the initial w and xbars.
594   """   
595   def iteration_0_solve(self):
596
597      if self._verbose is True:
598         print "------------------------------------------------"
599         print "Starting PH iteration 0 solves"
600
601      self._current_iteration = 0
602
603      solve_start_time = time.time()
604
605      # STEP 1: queue up the solves for all scenario sub-problems.
606      #         grab all the action handles for the subsequent barrier sync.
607
608      action_handles = []
609      action_handle_instance_map = {}
610
611      for scenario in self._scenario_tree._scenarios:
612         
613         instance = self._instances[scenario._name]
614         
615         if self._verbose is True:
616            print "Queuing solve for scenario=" + scenario._name
617           
618         # TBD - need to have the solver be able to solve a particular instance, with a specific objective
619
620         # there's nothing to warm-start from in iteration 0, so don't include the keyword.
621         # the reason you don't want to include it is that some solvers don't know how to handle
622         # the keyword at all (despite it being false). you might want to solve iteration 0 solves
623         # using some other solver.
624
625         new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
626         action_handle_instance_map[scenario._name] = new_action_handle
627
628         action_handles.append(new_action_handle)
629
630      # STEP 2: barrier sync for all scenario sub-problem solves.
631      if self._verbose is True:
632         print "Waiting for scenario sub-problem solves"
633      self._solver_manager.wait_all(action_handles)
634      if self._verbose is True:     
635         print "Scenario sub-problem solves completed"     
636
637      solve_end_time = time.time()
638      self._cumulative_solve_time += (solve_end_time - solve_start_time)
639
640      if self._output_times is True:
641         print "Aggregate sub-problem solve time=%8.2f" % (solve_end_time - solve_start_time)
642
643      # STEP 3: Load the results!
644      for scenario_name, action_handle in action_handle_instance_map.items():
645
646         if self._verbose is True:         
647            print "Processing results for scenario=", scenario_name, " completed successfully"         
648
649         instance = self._instances[scenario_name]
650         results = self._solver_manager.get_results(action_handle)
651
652         if len(results.solution) == 0:
653            results.write(num=1)
654            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
655
656         if self._output_solver_results is True:
657            print "Results for scenario=",scenario_name
658            results.write(num=1)           
659
660         instance.load(results)
661
662         if self._verbose is True:                 
663            print "Successfully loaded solution for scenario=",scenario_name
664
665
666      # TBD - save the solutions for history tracking
667
668      if self._verbose is True:
669         print "Successfully completed PH iteration 0 solves - solution statistics:"
670         print "         Scenario              Objective                  Value"
671         for scenario in self._scenario_tree._scenarios:
672            instance = self._instances[scenario._name]
673            for objective_name in instance._component[Objective]:
674               objective = instance._component[Objective][objective_name]
675               # TBD: I don't know how to deal with objective arrays, so I'll assume they aren't there
676               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
677         print "------------------------------------------------"
678
679   #
680   # recompute the averages, minimum, and maximum statistics for all variables to be blended by PH, i.e.,
681   # not appearing in the final stage. technically speaking, the min/max aren't required by PH, but they
682   # are used often enough to warrant their computation and it's basically free if you're computing the
683   # average.
684   #
685   def update_variable_statistics(self):
686
687      start_time = time.time()
688     
689      for stage in self._scenario_tree._stages:
690         
691         if stage != self._scenario_tree._stages[-1]: # no blending over the final stage
692           
693            for tree_node in stage._tree_nodes:
694               
695               for (variable, index_template, variable_indices) in stage._variables:
696                 
697                  variable_name = variable.name
698                     
699                  avg_parameter_name = "PHAVG_"+variable_name
700                 
701                  for index in variable_indices:
702                     min = float("inf")
703                     avg = 0.0
704                     max = float("-inf")
705                     node_probability = 0.0
706                     
707                     is_used = True # until proven otherwise                     
708                     for scenario in tree_node._scenarios:
709                       
710                        instance = self._instances[scenario._name]
711                       
712                        if getattr(instance,variable_name)[index].status == VarStatus.unused:
713                           is_used = False
714                        else:                       
715                           node_probability += scenario._probability
716                           var_value = getattr(instance, variable.name)[index].value
717                           if var_value < min:
718                              min = var_value
719                           avg += (scenario._probability * var_value)
720                           if var_value > max:
721                              max = var_value
722
723                     if is_used is True:
724                        tree_node._minimums[variable.name][index] = min
725                        tree_node._averages[variable.name][index] = avg / node_probability
726                        tree_node._maximums[variable.name][index] = max                       
727
728                        # distribute the newly computed average to the xbar variable in
729                        # each instance/scenario associated with this node. only do this
730                        # if the variable is used!
731                        for scenario in tree_node._scenarios:
732                           instance = self._instances[scenario._name]
733                           avg_parameter = getattr(instance, avg_parameter_name)
734                           avg_parameter[index] = avg / node_probability
735
736      end_time = time.time()
737      self._cumulative_xbar_time += (end_time - start_time)
738
739   def update_weights(self):
740
741      # because the weight updates rely on the xbars, and the xbars are node-based,
742      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
743      start_time = time.time()     
744
745      for stage in self._scenario_tree._stages:
746         
747         if stage != self._scenario_tree._stages[-1]: # no blending over the final stage, so no weights to worry about.
748           
749            for tree_node in stage._tree_nodes:
750
751               for (variable, index_template, variable_indices) in stage._variables:
752
753                  variable_name = variable.name
754                  blend_parameter_name = "PHBLEND_"+variable_name
755                  weight_parameter_name = "PHWEIGHT_"+variable_name
756                  rho_parameter_name = "PHRHO_"+variable_name
757
758                  for index in variable_indices:
759
760                     tree_node_average = tree_node._averages[variable.name][index]()
761
762                     for scenario in tree_node._scenarios:
763
764                        instance = self._instances[scenario._name]
765
766                        if getattr(instance,variable.name)[index].status != VarStatus.unused:
767
768                           # we are currently not updating weights if blending is disabled for a variable.
769                           # this is done on the premise that unless you are actively trying to move
770                           # the variable toward the mean, the weights will blow up and be huge by the
771                           # time that blending is activated.
772                           variable_blend_indicator = getattr(instance, blend_parameter_name)[index]()                             
773
774                           # get the weight and rho parameters for this variable/index combination.
775                           rho_value = getattr(instance, rho_parameter_name)[index]()
776                           current_variable_weight = getattr(instance, weight_parameter_name)[index]()
777                                               
778                           # if I'm maximizing, invert value prior to adding (hack to implement negatives).
779                           # probably fixed in Pyomo at this point - I just haven't checked in a long while.
780                           if self._is_minimizing is False:
781                              current_variable_weight = (-current_variable_weight)
782                           current_variable_value = getattr(instance,variable.name)[index]()
783                           new_variable_weight = current_variable_weight + variable_blend_indicator * rho_value * (current_variable_value - tree_node_average)
784                           # I have the correct updated value, so now invert if maximizing.
785                           if self._is_minimizing is False:
786                              new_variable_weight = (-new_variable_weight)
787                           getattr(instance, weight_parameter_name)[index].value = new_variable_weight
788
789      # we shouldn't have to re-simplify the expression, as we aren't adding any constant-variable terms - just modifying parameters.
790
791      end_time = time.time()
792      self._cumulative_weight_time += (end_time - start_time)
793
794   def form_iteration_k_objectives(self):
795
796      if self._verbose is True:
797         print "Forming PH-weighted objective"
798
799      # for each blended variable (i.e., those not appearing in the final stage),
800      # add the linear and quadratic penalty terms to the objective.
801      for instance_name, instance in self._instances.items():
802         
803         objective_name = instance._component[Objective].keys()[0]         
804         objective = instance._component[Objective][objective_name]         
805         objective_expression = objective._data[None].expr # TBD: we don't deal with indexed expressions (not even sure how to interpret them)
806         # the quadratic expression is really treated as just a list - eventually should be treated as a full expression.
807         quad_expression = 0.0
808         
809         for stage in self._scenario_tree._stages:
810
811            # skip the last stage, as no blending occurs
812            if stage != self._scenario_tree._stages[-1]:
813               # find the instance objective.
814               # TBD - for simiplicity, I'm assuming a single objective - we should select "active" objectives later.
815               #       also, can create a quadratic objective and leave the original alone.
816
817               for (reference_variable, index_template, variable_indices) in stage._variables:
818
819                  variable_name = reference_variable.name
820                  variable_type = reference_variable.domain
821
822                  w_parameter_name = "PHWEIGHT_"+variable_name
823                  w_parameter = instance._component[param._ParamBase][w_parameter_name]
824                 
825                  average_parameter_name = "PHAVG_"+variable_name
826                  average_parameter = instance._component[param._ParamBase][average_parameter_name]
827
828                  rho_parameter_name = "PHRHO_"+variable_name
829                  rho_parameter = instance._component[param._ParamBase][rho_parameter_name]
830
831                  blend_parameter_name = "PHBLEND_"+variable_name
832                  blend_parameter = instance._component[param._ParamBase][blend_parameter_name]                 
833
834                  instance_variable = instance._component[var._VarBase][variable_name]
835
836                  for index in variable_indices:
837
838                     if (instance_variable[index].status is not VarStatus.unused) and (instance_variable[index].fixed is False):
839
840                        # TBD - if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes.
841                        objective_expression += (w_parameter[index] * instance_variable[index])
842
843                        if isinstance(variable_type, BooleanSet) is True:
844
845                           if self._retain_quadratic_binary_terms is False:
846                              # this rather ugly form of the linearized quadratic expression term is required
847                              # due to a pyomo bug - the form (rho/2) * (x+y+z) chokes in presolve when distributing
848                              # over the sum.
849                              new_term = (blend_parameter[index] * rho_parameter[index] / 2.0 * instance_variable[index]) - \
850                                         (blend_parameter[index] * rho_parameter[index] * average_parameter[index] * instance_variable[index]) + \
851                                         (blend_parameter[index] * rho_parameter[index] / 2.0 * average_parameter[index] * average_parameter[index])                             
852                              if objective.sense is minimize:
853                                 objective_expression += new_term
854                              else:
855                                 objective_expression -= new_term                                 
856                           else:
857                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)
858
859                        else:
860
861                           quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)                           
862                   
863         # strictly speaking, this probably isn't necessary - parameter coefficients won't get
864         # pre-processed out of the expression tree. however, if the under-the-hood should change,
865         # we'll be covered.
866         objective_expression.simplify(instance)
867         instance._component[Objective][objective_name]._data[None].expr = objective_expression
868         # if we are linearizing everything, then nothing will appear in the quadratic expression -
869         # don't add the empty "0.0" expression to the objective.
870         if quad_expression != 0.0:
871           instance._component[Objective][objective_name]._quad_subexpr = quad_expression
872
873   def iteration_k_solve(self):
874
875     if self._verbose is True:
876        print "------------------------------------------------"       
877        print "Starting PH iteration " + str(self._current_iteration) + " solves"
878
879     # cache the objective values generated by PH for output at the end of this function.
880     ph_objective_values = {}
881
882     solve_start_time = time.time()
883
884     # STEP 1: queue up the solves for all scenario sub-problems.
885     #         grab all the action handles for the subsequent barrier sync.
886
887     action_handles = []
888     action_handle_instance_map = {}
889
890     for scenario in self._scenario_tree._scenarios:     
891
892        instance = self._instances[scenario._name]
893
894        if self._verbose is True:
895           print "Queuing solve for scenario=" + scenario._name
896
897        # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
898        # don't do this, you won't see any chance in the output files as you vary the problem parameters!
899        # ditto for instance fixing!
900        instance.presolve()
901
902        # TBD - need to have the solver be able to solve a particular instance, with a specific objective
903
904        # once past iteration 0, there is always a feasible solution from which to warm-start.
905        # however, you might want to disable warm-start when the solver is behaving badly (which does happen).
906        new_action_handle = None
907        if self._disable_warmstarts is False:
908           new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
909        else:
910           new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)           
911
912        action_handle_instance_map[scenario._name] = new_action_handle
913
914        action_handles.append(new_action_handle)
915
916     # STEP 2: barrier sync for all scenario sub-problem solves.
917     if self._verbose is True:           
918        print "Waiting for scenario sub-problem solves"
919     self._solver_manager.wait_all(action_handles)
920     if self._verbose is True:               
921        print "Scenario sub-problem solves completed"
922       
923     solve_end_time = time.time()
924     self._cumulative_solve_time += (solve_end_time - solve_start_time)
925
926     if self._output_times is True:
927        print "Aggregate sub-problem solve time=%8.2f" % (solve_end_time - solve_start_time)
928
929     # STEP 3: Load the results!
930     for scenario_name, action_handle in action_handle_instance_map.items():
931
932        instance = self._instances[scenario_name]
933        results = self._solver_manager.get_results(action_handle)
934
935        if len(results.solution) == 0:
936           raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
937
938        if self._output_solver_results is True:
939           print "Results for scenario=", scenario_name
940           results.write(num=1)           
941
942        instance.load(results)
943
944        if self._verbose is True:                 
945           print "Successfully loaded solution for scenario=", scenario_name
946
947        # we're assuming there is a single solution.
948        # the "value" attribute is a pre-defined feature of any solution - it is relative to whatever
949        # objective was selected during optimization, which of course should be the PH objective.
950        ph_objective_values[instance.name] = float(results.solution(0).value)
951
952     if self._verbose is True:
953        print "Successfully completed PH iteration " + str(self._current_iteration) + " solves - solution statistics:"
954        print "  Scenario             PH Objective             Cost Objective"
955        for scenario in self._scenario_tree._scenarios:
956           instance = self._instances[scenario._name]
957           for objective_name in instance._component[Objective]:
958              objective = instance._component[Objective][objective_name]
959              # TBD: I don't know how to deal with objective arrays, so I'll assume they aren't there
960              print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], 0.0)
961
962   def solve(self):
963
964      self._solve_start_time = time.time()
965      self._cumulative_solve_time = 0.0
966      self._cumulative_xbar_time = 0.0
967      self._cumulative_weight_time = 0.0
968
969      print "Starting PH"
970
971      if self._initialized == False:
972         raise RuntimeError, "PH is not initialized - cannot invoke solve() method"
973
974      print "Initiating PH iteration=" + `self._current_iteration`
975
976      self.iteration_0_solve()
977
978      # update variable statistics prior to any output.
979      self.update_variable_statistics()     
980
981      if self._verbose is True:
982         print "Variable values following scenario solves:"
983         self.pprint(False,False,True,False)
984
985      # let plugins know if they care.
986      if len(self._ph_plugin) == 1:
987         self._ph_plugin.service().post_iteration_0_solves(self)
988
989      # update the fixed variable statistics.
990      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
991
992      if self._verbose is True:
993         print "Number of discrete variables fixed=",self._total_fixed_discrete_vars," (total=",self._total_discrete_vars,")"
994         print "Number of continuous variables fixed=",self._total_fixed_continuous_vars," (total=",self._total_continuous_vars,")"
995
996      self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
997      print "Convergence metric=%12.4f" % self._converger.lastMetric()
998
999      self.update_weights()
1000
1001      if self._max_iterations > 0:
1002         self.form_iteration_k_objectives()
1003
1004      # let plugins know if they care.
1005      if len(self._ph_plugin) == 1:
1006         self._ph_plugin.service().post_iteration_0(self)
1007
1008      # checkpoint if it's time - which it always is after iteration 0,
1009      # if the interval is >= 1!
1010      if (self._checkpoint_interval > 0):
1011         self.checkpoint(0)
1012
1013      # there is an upper bound on the number of iterations to execute -
1014      # the actual bound depends on the converger supplied by the user.
1015      for i in range(1, self._max_iterations+1):
1016
1017         self._current_iteration = self._current_iteration + 1                 
1018
1019         print "Initiating PH iteration=" + `self._current_iteration`         
1020
1021         if self._verbose is True:
1022            print "Variable averages and weights prior to scenario solves:"
1023            self.pprint(True,True,False,False)
1024
1025         self.iteration_k_solve()
1026
1027         # update variable statistics prior to any output.
1028         self.update_variable_statistics()
1029         
1030         if self._verbose is True:
1031            print "Variable values following scenario solves:"
1032            self.pprint(False,False,True,False)
1033
1034         # we don't technically have to do this at the last iteration,
1035         # but with checkpointing and re-starts, you're never sure
1036         # when you're executing the last iteration.
1037         self.update_weights()
1038
1039         # let plugins know if they care.
1040         if len(self._ph_plugin) == 1:
1041            self._ph_plugin.service().post_iteration_k_solves(self)
1042
1043         # update the fixed variable statistics.
1044         (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
1045
1046         if self._verbose is True:
1047            print "Number of discrete variables fixed=",self._total_fixed_discrete_vars," (total=",self._total_discrete_vars,")"
1048            print "Number of continuous variables fixed=",self._total_fixed_continuous_vars," (total=",self._total_continuous_vars,")"
1049
1050         # let plugins know if they care.
1051         if len(self._ph_plugin) == 1:
1052            self._ph_plugin.service().post_iteration_k(self)
1053
1054         # at this point, all the real work of an iteration is complete.
1055
1056         # checkpoint if it's time.
1057         if (self._checkpoint_interval > 0) and (i % self._checkpoint_interval is 0):
1058            self.checkpoint(i)
1059
1060         # check for early termination.
1061         self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1062         print "Convergence metric=%12.4f" % self._converger.lastMetric()
1063
1064         if self._converger.isConverged(self) is True:
1065            if self._total_discrete_vars == 0:
1066               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)
1067            else:
1068               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)+" or all discrete variables are fixed"               
1069            break
1070
1071         # if we're terminating due to exceeding the maximum iteration count, print a message
1072         # indicating so - otherwise, you get a quiet, information-free output trace.
1073         if i == self._max_iterations:
1074            print "Halting PH - reached maximal iteration count="+str(self._max_iterations)
1075
1076      print "PH complete"
1077
1078      if self._verbose is True:
1079         print "Convergence history:"
1080         self._converger.pprint()
1081
1082      print "Final variable values:"
1083      self.pprint(False,False,True,True)         
1084
1085      print "Final costs:"
1086      self._scenario_tree.pprintCosts(self._instances)
1087
1088      self._solve_end_time = time.time()
1089
1090      if (self._verbose is True) and (self._output_times is True):
1091         print "Overall run-time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1092
1093      # let plugins know if they care.
1094      if len(self._ph_plugin) == 1:
1095         self._ph_plugin.service().post_ph_execution(self)                                 
1096
1097   #
1098   # prints a summary of all collected time statistics
1099   #
1100   def print_time_stats(self):
1101
1102      print "PH run-time statistics (user):"
1103
1104      print "Initialization time=  %8.2f seconds" % (self._init_end_time - self._init_start_time)
1105      print "Overall solve time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1106      print "Scenario solve time=  %8.2f seconds" % self._cumulative_solve_time
1107      print "Average update time=  %8.2f seconds" % self._cumulative_xbar_time
1108      print "Weight update time=   %8.2f seconds" % self._cumulative_weight_time
1109
1110   #
1111   # a utility to determine whether to output weight / average / etc. information for
1112   # a variable/node combination. when the printing is moved into a callback/plugin,
1113   # this routine will go there. for now, we don't dive down into the node resolution -
1114   # just the variable/stage.
1115   #
1116   def should_print(self, stage, variable, variable_indices):
1117
1118      if self._output_continuous_variable_stats is False:
1119
1120         variable_type = variable.domain         
1121
1122         if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):
1123
1124            return False
1125
1126      return True
1127     
1128   #
1129   # pretty-prints the state of the current variable averages, weights, and values.
1130   # inputs are booleans indicating which components should be output.
1131   #
1132   def pprint(self, output_averages, output_weights, output_values, output_fixed):
1133
1134      # TBD - write a utility routine to figure out the longest identifier width for indicies and other names,
1135      #       to make the output a bit more readable.
1136
1137      if self._initialized is False:
1138         raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"         
1139     
1140      # print tree nodes and associated variable/xbar/ph information in stage-order
1141      for stage in self._scenario_tree._stages:
1142
1143         # we don't blend in the last stage, so we don't current care about printing the associated information.
1144         if stage != self._scenario_tree._stages[-1]:
1145
1146            print "\tStage=" + stage._name
1147
1148            num_outputs_this_stage = 0 # tracks the number of outputs on a per-index basis.
1149
1150            for (variable, index_template, variable_indices) in stage._variables:
1151
1152               variable_name = variable.name
1153
1154               if self.should_print(stage, variable, variable_indices) is True:
1155
1156                  num_outputs_this_variable = 0 # track, so we don't output the variable names unless there is an entry to report.
1157
1158                  for index in variable_indices:               
1159
1160                     weight_parameter_name = "PHWEIGHT_"+variable_name
1161
1162                     num_outputs_this_index = 0 # track, so we don't output the variable index more than once.
1163
1164                     for tree_node in stage._tree_nodes:                 
1165
1166                        # determine if the variable/index pair is used across the set of scenarios (technically,
1167                        # it should be good enough to check one scenario). ditto for "fixed" status. fixed does
1168                        # imply unused (see note below), but we care about the fixed status when outputting
1169                        # final solutions.
1170
1171                        is_used = True # should be consistent across scenarios, so one "unused" flags as invalid.
1172                        is_fixed = False
1173
1174                        for scenario in tree_node._scenarios:
1175                           instance = self._instances[scenario._name]
1176                           variable_value = getattr(instance,variable_name)[index]
1177                           if variable_value.status == VarStatus.unused:
1178                              is_used = False
1179                           if variable_value.fixed is True:
1180                              is_fixed = True
1181 
1182                        # IMPT: this is far from obvious, but variables that are fixed will - because
1183                        #       presolve will identify them as constants and eliminate them from all
1184                        #       expressions - be flagged as "unused" and therefore not output.
1185
1186                        if ((output_fixed is True) and (is_fixed is True)) or (is_used is True):
1187
1188                              minimum_value = tree_node._minimums[variable_name][index]()
1189                              maximum_value = tree_node._maximums[variable_name][index]()
1190
1191                              num_outputs_this_stage = num_outputs_this_stage + 1                           
1192                              num_outputs_this_variable = num_outputs_this_variable + 1
1193                              num_outputs_this_index = num_outputs_this_index + 1
1194
1195                              if num_outputs_this_variable == 1:
1196                                 print "\t\tVariable=",variable_name
1197
1198                              if num_outputs_this_index == 1:
1199                                 print "\t\t\tIndex:", indexToString(index)                             
1200
1201                              print "\t\t\t\tTree Node=",tree_node._name,"\t\t (Scenarios: ",                             
1202                              for scenario in tree_node._scenarios:
1203                                 print scenario._name," ",
1204                                 if scenario == tree_node._scenarios[-1]:
1205                                    print ")"
1206                           
1207                              if output_values is True:
1208                                 average_value = tree_node._averages[variable_name][index]()
1209                                 print "\t\t\t\tValues: ",                       
1210                                 for scenario in tree_node._scenarios:
1211                                    instance = self._instances[scenario._name]
1212                                    this_value = getattr(instance,variable_name)[index].value
1213                                    print "%12.4f" % this_value,
1214                                    if scenario == tree_node._scenarios[-1]:
1215                                       print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1216                                       print "    Avg=%12.4f" % (average_value),
1217                                       print ""
1218                              if output_weights:
1219                                 print "\t\t\t\tWeights: ",
1220                                 for scenario in tree_node._scenarios:
1221                                    instance = self._instances[scenario._name]
1222                                    print "%12.4f" % getattr(instance,weight_parameter_name)[index].value,
1223                                    if scenario == tree_node._scenarios[-1]:
1224                                       print ""
1225                              if output_averages:
1226                                 print "\t\t\t\tAverage: %12.4f" % (tree_node._averages[variable_name][index].value)
1227
1228            if num_outputs_this_stage == 0:
1229               print "\t\tNo non-converged variables in stage"
1230
1231            # cost variables aren't blended, so go through the gory computation of min/max/avg.
1232            # TBD: possibly move these into the tree node anyway, as they are useful in general.
1233            # TBD: these should also be probability-weighted - fix with above-mentioned TBD fix.
1234            # we currently always print these.
1235            cost_variable_name = stage._cost_variable[0].name
1236            cost_variable_index = stage._cost_variable[1]
1237            if cost_variable_index is None:
1238               print "\t\tCost Variable=" + cost_variable_name
1239            else:
1240               print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)           
1241            for tree_node in stage._tree_nodes:
1242               print "\t\t\tTree Node=" + tree_node._name + "\t\t (Scenarios: ",
1243               for scenario in tree_node._scenarios:
1244                  print scenario._name," ",
1245                  if scenario == tree_node._scenarios[-1]:
1246                     print ")"
1247               maximum_value = 0.0
1248               minimum_value = 0.0
1249               sum_values = 0.0
1250               num_values = 0
1251               first_time = True
1252               print "\t\t\tValues: ",                       
1253               for scenario in tree_node._scenarios:
1254                   instance = self._instances[scenario._name]
1255                   this_value = getattr(instance,cost_variable_name)[cost_variable_index].value
1256                   print "%12.4f" % this_value,
1257                   num_values += 1
1258                   sum_values += this_value
1259                   if first_time is True:
1260                      first_time = False
1261                      maximum_value = this_value
1262                      minimum_value = this_value
1263                   else:
1264                      if this_value > maximum_value:
1265                         maximum_value = this_value
1266                      if this_value < minimum_value:
1267                         minimum_value = this_value
1268                   if scenario == tree_node._scenarios[-1]:
1269                      print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1270                      print "    Avg=%12.4f" % (sum_values/num_values),
1271                      print ""
1272           
Note: See TracBrowser for help on using the repository browser.