source: coopr.pysp/stable/2.1/coopr/pysp/asynchph.py @ 2068

Last change on this file since 2068 was 2068, checked in by wehart, 10 years ago

Merged revisions 1952-2067 via svnmerge from
https://software.sandia.gov/svn/public/coopr/coopr.pysp/trunk

........

r1956 | jwatson | 2009-12-02 17:56:53 -0700 (Wed, 02 Dec 2009) | 3 lines


Added --scenario-solver-options and --ef-solver-options options to the "runph" script.

........

r1957 | dlwoodr | 2009-12-03 14:17:35 -0700 (Thu, 03 Dec 2009) | 2 lines


Documentation updates for pysp

........

r1974 | wehart | 2009-12-06 17:20:56 -0700 (Sun, 06 Dec 2009) | 2 lines


Updating PyPI categories

........

r1978 | jwatson | 2009-12-10 21:29:33 -0700 (Thu, 10 Dec 2009) | 3 lines


Eliminated exception-handling logic when loading user-defined extension modules in PH. The range of exceptions is too large, and for debugging purposes, it is more useful to see the raw trace output.

........

r1979 | jwatson | 2009-12-10 22:23:17 -0700 (Thu, 10 Dec 2009) | 5 lines


Biggest enhancement: The efwriter command-line script now has the option to output a CVaR-weighted objective term. Not extensively tested, but behaves sane on a number of small test cases.


Improved exception handling and error diagnostics in both the runph and efwriter scripts.

........

r1985 | jwatson | 2009-12-12 10:45:17 -0700 (Sat, 12 Dec 2009) | 3 lines


Modified PH to only use warm-starts if a solver has the capability!

........

r1998 | jwatson | 2009-12-13 15:17:58 -0700 (Sun, 13 Dec 2009) | 3 lines


Changed references to _component to active_component.

........

r2026 | wehart | 2009-12-21 23:27:06 -0700 (Mon, 21 Dec 2009) | 2 lines


Attempting to update PH. I'm not sure if this works, since I don't know how to test PH.

........

r2029 | jwatson | 2009-12-22 09:52:21 -0700 (Tue, 22 Dec 2009) | 3 lines


Some fixes to the ef writer based on Bill's recent changes to _initialize_constraint.

........

r2035 | jwatson | 2009-12-22 21:10:32 -0700 (Tue, 22 Dec 2009) | 3 lines


Added --scenario-mipgap option to PH script. Added _mipgap attribute to PH object. This attribute is mirrored to the solver plugin at the initiation of each iteration, after any PH extensions have the opportunity to provide a new value to the attribute. It is currently made use of by the WW PH extension.

........

r2037 | dlwoodr | 2009-12-23 14:38:43 -0700 (Wed, 23 Dec 2009) | 2 lines


add this example from Fernando

........

r2038 | dlwoodr | 2009-12-23 14:46:56 -0700 (Wed, 23 Dec 2009) | 3 lines


finish the job: we can now duplicate Fernando's example

........

r2039 | jwatson | 2009-12-23 15:13:54 -0700 (Wed, 23 Dec 2009) | 3 lines


Missed fix with new constraint initialization syntax in PH linearization.

........

r2058 | jwatson | 2009-12-29 10:57:58 -0700 (Tue, 29 Dec 2009) | 3 lines


Missed some _initialize_constraint function calls within the PySP EF writer during the recent switch to the corresponding "add" method.

........

r2059 | jwatson | 2009-12-29 10:58:34 -0700 (Tue, 29 Dec 2009) | 3 lines


Enabling garbage collection by default in PH.

........

r2060 | jwatson | 2009-12-29 10:59:05 -0700 (Tue, 29 Dec 2009) | 3 lines


Elimnating some debug output.

........

r2061 | jwatson | 2009-12-29 11:07:47 -0700 (Tue, 29 Dec 2009) | 3 lines


Fixing some option documentation in PH.

........

r2062 | jwatson | 2009-12-29 12:00:37 -0700 (Tue, 29 Dec 2009) | 3 lines


Added ef-mipgap option to PH scripts.

........

File size: 49.3 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 pyutilib.plugin.core
13import types
14from coopr.pyomo import *
15import copy
16import os.path
17import traceback
18import copy
19from coopr.opt import SolverResults,SolverStatus
20from coopr.opt.base import SolverFactory
21from coopr.opt.parallel import SolverManagerFactory
22import time
23import types
24
25from scenariotree import *
26from phutils import *
27
28from pyutilib.plugin.core import ExtensionPoint
29
30from coopr.pysp.phextension import IPHExtension
31
32class AsyncProgressiveHedging(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                     # determine if this index is used - otherwise, don't waste time
88                     # fixing and cycle checking. for one, the code will crash :-) with
89                     # None values during the cycle checking computation!
90
91                     is_used = True # until proven otherwise                     
92                     for scenario in tree_node._scenarios:
93                        instance = self._instances[scenario._name]
94                        if getattr(instance,variable_name)[index].status == VarStatus.unused:
95                           is_used = False
96
97                     if is_used is True:                       
98
99                        # JPW TBD: ideally, we want an "is-discrete" check in the logic below, which COOPR doesn't currently support.
100                        if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
101                           num_discrete_vars = num_discrete_vars + 1
102                        else:
103                           num_continuous_vars = num_continuous_vars + 1
104
105      return (num_discrete_vars, num_continuous_vars)
106
107   """ Constructor
108       Arguments:
109          max_iterations        the maximum number of iterations to run PH (>= 0). defaults to 0.
110          rho                   the global rho value (> 0). defaults to 0.
111          rho_setter            an optional name of a python file used to set particular variable rho values.
112          solver                the solver type that PH uses to solve scenario sub-problems. defaults to "cplex".
113          solver_manager        the solver manager type that coordinates scenario sub-problem solves. defaults to "serial".
114          keep_solver_files     do I keep intermediate solver files around (for debugging)? defaults to False.
115          output_solver_log     do I dump the solver log (as it is being generated) to the screen? defaults to False.
116          output_solver_results do I output (for debugging) the detailed solver results, including solutions, for scenario solves? defaults to False.
117          verbose               does the PH object stream debug/status output? defaults to False.
118          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).
119
120   """
121   def __init__(self, *args, **kwds):
122
123      # PH configuration parameters
124      # TBD - shift rho to a parameter of the instances (probably keep it here, and in the model)
125      self._rho = 0.0 # TBD morph to per-variable later, and possibly per-variable, per-scenario.
126      self._rho_setter = None # filename for the modeler to set rho.
127      self._max_iterations = 0
128
129      # PH reporting parameters
130      self._verbose = False # do I flood the screen with status output?
131      self._output_continuous_variable_stats = True # when in verbose mode, do I output weights/averages for continuous variables?
132      self._output_solver_results = False
133      self._output_times = False
134
135      # PH run-time variables
136      self._current_iteration = 0 # the 'k'
137      self._xbar = {} # current per-variable averages. maps (node_id, variable_name) -> value
138      self._initialized = False # am I ready to call "solve"? Set to True by the initialize() method.
139
140      # PH solver information / objects.
141      self._solver_type = "cplex"
142      self._solver_manager_type = "pyro" # only pyro makes sense for async currently.
143     
144      self._solver = None # will eventually be unnecessary once Bill eliminates the need for a solver in the solver manager constructor.
145      self._solver_manager = None 
146     
147      self._keep_solver_files = False
148      self._output_solver_log = False
149
150      # PH convergence computer/updater.
151      self._converger = None
152
153      # PH history
154      self._solutions = {}
155
156      # all information related to the scenario tree (implicit and explicit).
157      self._model = None # not instantiated
158      self._model_instance = None # instantiated
159
160      self._scenario_tree = None
161     
162      self._scenario_data_directory = "" # this the prefix for all scenario data
163      self._instances = {} # maps scenario name to the corresponding model instance
164
165      # for various reasons (mainly hacks at this point), it's good to know whether we're minimizing or maximizing.
166      self._is_minimizing = None
167
168      # global handle to ph extension plugin
169      self._ph_plugin = ExtensionPoint(IPHExtension)
170
171      # PH timing statistics - relative to last invocation.
172      self._init_start_time = None # for initialization() method
173      self._init_end_time = None
174      self._solve_start_time = None # for solve() method
175      self._solve_end_time = None
176      self._cumulative_solve_time = None # seconds, over course of solve()
177      self._cumulative_xbar_time = None # seconds, over course of update_xbars()
178      self._cumulative_weight_time = None # seconds, over course of update_weights()
179
180      # PH default tolerances - for use in fixing and testing equality across scenarios,
181      # and other stuff.
182      self._integer_tolerance = 0.00001
183
184      # process the keyword options
185      for key in kwds.keys():
186         if key == "max_iterations":
187            self._max_iterations = kwds[key]
188         elif key == "rho":
189            self._rho = kwds[key]
190         elif key == "rho_setter":
191            self._rho_setter = kwds[key]           
192         elif key == "solver":
193            self._solver_type = kwds[key]
194         elif key == "solver_manager":
195            self._solver_manager_type = kwds[key]           
196         elif key == "keep_solver_files":
197            self._keep_solver_files = kwds[key]
198         elif key == "output_solver_results":
199            self._output_solver_results = kwds[key]
200         elif key == "output_solver_log":
201            self._output_solver_log = kwds[key]
202         elif key == "verbose":
203            self._verbose = kwds[key]
204         elif key == "output_times":
205            self._output_times = kwds[key]           
206         else:
207            print "Unknown option=" + key + " specified in call to PH constructor"
208
209      # validate all "atomic" options (those that can be validated independently)
210      if self._max_iterations < 0:
211         raise ValueError, "Maximum number of PH iterations must be non-negative; value specified=" + `self._max_iterations`
212      if self._rho <= 0.0:
213         raise ValueError, "Value of rho paramter in PH must be non-zero positive; value specified=" + `self._rho`
214
215      # validate rho setter file if specified.
216      if self._rho_setter is not None:
217         if os.path.exists(self._rho_setter) is False:
218            raise ValueError, "The rho setter script file="+self._rho_setter+" does not exist"
219
220      # construct the sub-problem solver.
221      self._solver = SolverFactory(self._solver_type)
222      if self._solver == None:
223         raise ValueError, "Unknown solver type=" + self._solver_type + " specified in call to PH constructor"
224      if self._keep_solver_files is True:
225         # TBD - this is a bit kludgy, as it requires PH to know/assume that
226         #       it is dealing with a system call solver. the solver factory should take keyword options as a fix.
227         self._solver.keepFiles = True
228
229      # construct the solver manager.
230      if self._solver_manager_type != "pyro":
231         raise ValueError, "Only the pyro solver manager makes sense for asynch PH!"
232      print "SOLVER MANAGER TYPE=",self._solver_manager_type
233      self._solver_manager = SolverManagerFactory(self._solver_manager_type)
234      if self._solver_manager is None:
235         raise ValueError, "Failed to create solver manager of type="+self._solver_manager_type+" specified in call to PH constructor"
236
237      # a set of all valid PH iteration indicies is generally useful for plug-ins, so create it here.
238      self._iteration_index_set = Set(name="PHIterations")
239      for i in range(0,self._max_iterations + 1):
240         self._iteration_index_set.add(i)     
241
242      # spit out parameterization if verbosity is enabled
243      if self._verbose == True:
244         print "PH solver configuration: "
245         print "   Max iterations=" + `self._max_iterations`
246         print "   Rho=" + `self._rho`
247         print "   Sub-problem solver type=" + `self._solver_type`
248
249   """ Initialize PH with model and scenario data, in preparation for solve().
250       Constructs and reads instances.
251   """
252   def initialize(self, scenario_data_directory_name=".", model=None, model_instance=None, scenario_tree=None, converger=None):
253
254      self._init_start_time = time.time()
255
256      if self._verbose == True:
257         print "Initializing PH"
258         print "   Scenario data directory=" + scenario_data_directory_name
259
260      if not os.path.exists(scenario_data_directory_name):
261         raise ValueError, "Scenario data directory=" + scenario_data_directory_name + " either does not exist or cannot be read"
262
263      self._scenario_data_directory_name = scenario_data_directory_name
264
265      # IMPT: The input model should be an *instance*, as it is very useful (critical!) to know
266      #       the dimensions of sets, be able to store suffixes on variable values, etc.
267      # TBD: There may be a way to see if a model is initialized - throw an exception if it is not!
268      if model is None:
269         raise ValueError, "A model must be supplied to the PH initialize() method"
270
271      if scenario_tree is None:
272         raise ValueError, "A scenario tree must be supplied to the PH initialize() method"
273
274      if converger is None:
275         raise ValueError, "A convergence computer must be supplied to the PH initialize() method"
276
277      self._model = model
278      self._model_instance = model_instance
279      self._scenario_tree = scenario_tree
280      self._converger = converger
281
282      model_objective = model.active_components(Objective)
283      self._is_minimizing = (model_objective[ model_objective.keys()[0] ].sense == minimize)
284
285      self._converger.reset()
286
287      # construct instances for each scenario
288      if self._verbose is True:
289         if self._scenario_tree._scenario_based_data == 1:
290            print "Scenario-based instance initialization enabled"
291         else:
292            print "Node-based instance initialization enabled"
293         
294      for scenario in self._scenario_tree._scenarios:
295
296         scenario_instance = None
297
298         if self._verbose is True:
299            print "Creating instance for scenario=" + scenario._name
300
301         try:
302            if self._scenario_tree._scenario_based_data == 1:
303               scenario_data_filename = self._scenario_data_directory_name + os.sep + scenario._name + ".dat"
304               if self._verbose is True:
305                  print "Data for scenario=" + scenario._name + " loads from file=" + scenario_data_filename
306               scenario_instance = (self._model).create(scenario_data_filename)
307            else:
308               scenario_instance = self._model.clone()
309               scenario_data = ModelData()
310               current_node = scenario._leaf_node
311               while current_node is not None:
312                  node_data_filename = self._scenario_data_directory_name + os.sep + current_node._name + ".dat"
313                  if self._verbose is True:
314                     print "Node data for scenario=" + scenario._name + " partially loading from file=" + node_data_filename
315                  scenario_data.add_data_file(node_data_filename)
316                  current_node = current_node._parent
317               scenario_data.read(model=scenario_instance)
318               scenario_instance._load_model_data(scenario_data)
319               scenario_instance.presolve()
320         except:
321            print "Encountered exception in model instance creation - traceback:"
322            traceback.print_exc()
323            raise RuntimeError, "Failed to create model instance for scenario=" + scenario._name
324
325         self._instances[scenario._name] = scenario_instance
326         self._instances[scenario._name].name = scenario._name
327
328      if self._verbose is True:
329         print "PH successfully created model instances for all scenarios"
330
331      # TBD - Technically, we don't need to add these parameters to the models until after the iteration 0 solve (move later).
332
333      # create PH weight and xbar vectors, on a per-scenario basis, for each variable that is not in the
334      # final stage, i.e., for all variables that are being blended by PH. the parameters are created
335      # in the space of each scenario instance, so that they can be directly and automatically
336      # incorporated into the (appropriately modified) objective function.
337
338      if self._verbose is True:
339         print "Creating weight, average, and rho parameter vectors for scenario instances"
340
341      for (instance_name, instance) in self._instances.items():
342
343         # first, gather all unique variables referenced in any stage
344         # other than the last, independent of specific indices. this
345         # "gather" step is currently required because we're being lazy
346         # in terms of index management in the scenario tree - which
347         # should really be done in terms of sets of indices.
348         # NOTE: technically, the "instance variables" aren't really references
349         # to the variable in the instance - instead, the reference model. this
350         # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
351         instance_variables = {}
352         
353         for stage in self._scenario_tree._stages:
354           
355            if stage != self._scenario_tree._stages[-1]:
356               
357               for (reference_variable, index_template, reference_indices) in stage._variables:
358                 
359                  if reference_variable.name not in instance_variables.keys():
360                     
361                     instance_variables[reference_variable.name] = reference_variable
362
363         # for each blended variable, create a corresponding ph weight and average parameter in the instance.
364         # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor
365         # in the grand scheme of things.
366
367         for (variable_name, reference_variable) in instance_variables.items():
368
369            new_w_index = reference_variable._index # TBD - need to be careful with the shallow copy here
370            new_w_parameter_name = "PHWEIGHT_"+reference_variable.name
371            new_w_parameter = Param(new_w_index,name=new_w_parameter_name)
372            setattr(instance,new_w_parameter_name,new_w_parameter)
373
374            # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference
375            # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not
376            # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the
377            # issue in any case for now.
378            for index in new_w_index:
379               new_w_parameter[index] = 0.0
380
381            new_avg_index = reference_variable._index # TBD - need to be careful with the shallow copy here
382            new_avg_parameter_name = "PHAVG_"+reference_variable.name
383            new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
384            setattr(instance,new_avg_parameter_name,new_avg_parameter)
385
386            for index in new_avg_index:
387               new_avg_parameter[index] = 0.0
388
389            new_rho_index = reference_variable._index # TBD - need to be careful with the shallow copy here
390            new_rho_parameter_name = "PHRHO_"+reference_variable.name
391            new_rho_parameter = Param(new_rho_index,name=new_rho_parameter_name)
392            setattr(instance,new_rho_parameter_name,new_rho_parameter)
393
394            for index in new_rho_index:
395               new_rho_parameter[index] = self._rho
396
397      # if specified, run the user script to initialize variable rhos at their whim.
398      if self._rho_setter is not None:
399
400         print "Executing user rho set script from filename=", self._rho_setter
401         execfile(self._rho_setter)
402
403      # create parameters to store variable statistics (of general utility) at each node in the scenario tree.
404
405      if self._verbose is True:
406         print "Creating variable statistic (min/avg/max) parameter vectors for scenario tree nodes"
407
408      for stage in self._scenario_tree._stages:
409         if stage != self._scenario_tree._stages[-1]:
410
411            # first, gather all unique variables referenced in this stage
412            # this "gather" step is currently required because we're being lazy
413            # in terms of index management in the scenario tree - which
414            # should really be done in terms of sets of indices.
415            stage_variables = {}
416            for (reference_variable, index_template, reference_index) in stage._variables:
417               if reference_variable.name not in stage_variables.keys():
418                  stage_variables[reference_variable.name] = reference_variable
419
420            # next, create min/avg/max parameters for each variable in the corresponding tree node.
421            # NOTE: the parameter names below could really be empty, as they are never referenced
422            #       explicitly.
423            for (variable_name, reference_variable) in stage_variables.items():
424               for tree_node in stage._tree_nodes:
425
426                  new_min_index = reference_variable._index # TBD - need to be careful with the shallow copy here (and below)
427                  new_min_parameter_name = "NODEMIN_"+reference_variable.name
428                  new_min_parameter = Param(new_min_index,name=new_min_parameter_name)
429                  for index in new_min_index:
430                     new_min_parameter[index] = 0.0
431                  tree_node._minimums[reference_variable.name] = new_min_parameter                                   
432
433                  new_avg_index = reference_variable._index
434                  new_avg_parameter_name = "NODEAVG_"+reference_variable.name
435                  new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
436                  for index in new_avg_index:
437                     new_avg_parameter[index] = 0.0
438                  tree_node._averages[reference_variable.name] = new_avg_parameter
439
440                  new_max_index = reference_variable._index
441                  new_max_parameter_name = "NODEMAX_"+reference_variable.name
442                  new_max_parameter = Param(new_max_index,name=new_max_parameter_name)
443                  for index in new_max_index:
444                     new_max_parameter[index] = 0.0
445                  tree_node._maximums[reference_variable.name] = new_max_parameter                                                     
446
447
448      # indicate that we're ready to run.
449      self._initialized = True
450
451      self._init_end_time = time.time()
452
453      if self._verbose is True:
454         print "PH is successfully initialized"
455         if self._output_times is True:
456            print "Initialization time=" + str(self._init_end_time - self._init_start_time) + " seconds"
457
458      # let plugins know if they care.
459      if len(self._ph_plugin) == 1:
460         self._ph_plugin.service().post_ph_initialization(self)
461
462   """ Perform the non-weighted scenario solves and form the initial w and xbars.
463   """   
464   def iteration_0_solve(self):
465
466      if self._verbose == True:
467         print "------------------------------------------------"
468         print "Starting PH iteration 0 solves"
469
470      self._current_iteration = 0
471
472      solve_start_time = time.time()                     
473
474      # STEP 1: queue up the solves for all scenario sub-problems.
475      #         grab all the action handles for the subsequent barrier sync.
476
477      action_handles = []
478      action_handle_instance_map = {} 
479
480      for scenario in self._scenario_tree._scenarios:
481         
482         instance = self._instances[scenario._name]
483         
484         if self._verbose == True:
485            print "Queuing solve for scenario=" + scenario._name
486           
487         # TBD - need to have the solver be able to solve a particular instance, with a specific objective
488
489         # there's nothing to warm-start from in iteration 0.
490         new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=False, tee=self._output_solver_log)
491
492         action_handle_instance_map[scenario._name] = new_action_handle
493
494         action_handles.append(new_action_handle)
495
496      # STEP 2: barrier sync for all scenario sub-problem solves.
497      print "Waiting for scenario sub-problem solves"
498      self._solver_manager.wait_all(action_handles)
499
500      solve_end_time = time.time()
501      self._cumulative_solve_time += (solve_end_time - solve_start_time)
502
503      # STEP 3: Load the results!
504      for scenario_name, action_handle in action_handle_instance_map.items():
505
506         instance = self._instances[scenario_name]
507         results = self._solver_manager.get_results(action_handle)
508
509         if len(results.solution) == 0:
510            raise RuntimeError, "Solve failed for scenario="+scenario._name+"; no solutions generated"
511
512         if self._verbose is True:         
513            print "Solve completed successfully"
514
515#         if self._output_times is True:
516#            print "Solve time=%8.2f" % (solve_end_time - solve_start_time)
517
518#         if self._output_solver_results == True:
519         print "Results:"
520         results.write(num=1)           
521
522         instance.load(results)
523
524         if self._verbose == True:                 
525            print "Successfully loaded solution"
526
527      # TBD - save the solutions for history tracking
528
529      if self._verbose == True:
530         print "Successfully completed PH iteration 0 solves"
531         print "         Scenario              Objective                  Value"
532         for scenario in self._scenario_tree._scenarios:
533            instance = self._instances[scenario._name]
534            for objective_name in instance.active_components(Objective):
535               objective = instance.active_components(Objective)objective_name]
536               # TBD: I don't know how to deal with objective arrays, so I'll assume they aren't there
537               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
538         print "------------------------------------------------"
539         print ""
540
541   #
542   # recompute the averages, minimum, and maximum statistics for all variables to be blended by PH, i.e.,
543   # not appearing in the final stage. technically speaking, the min/max aren't required by PH, but they
544   # are used often enough to warrant their computation and it's basically free if you're computing the
545   # average.
546   #
547   def update_variable_statistics(self):
548
549      start_time = time.time()
550     
551      for stage in self._scenario_tree._stages:
552         
553         if stage != self._scenario_tree._stages[-1]: # no blending over the final stage
554           
555            for tree_node in stage._tree_nodes:
556               
557               for (variable, index_template, variable_indices) in stage._variables:
558                 
559                  variable_name = variable.name
560                     
561                  avg_parameter_name = "PHAVG_"+variable_name
562                 
563                  for index in variable_indices:
564                     min = float("inf")
565                     avg = 0.0
566                     max = float("-inf")
567                     node_probability = 0.0
568                     
569                     is_used = True # until proven otherwise                     
570                     for scenario in tree_node._scenarios:
571                       
572                        instance = self._instances[scenario._name]
573                       
574                        if getattr(instance,variable_name)[index].status == VarStatus.unused:
575                           is_used = False
576                        else:
577                           node_probability += scenario._probability
578                           var_value = getattr(instance, variable.name)[index].value
579                           if var_value < min:
580                              min = var_value
581                           avg += (scenario._probability * var_value)
582                           if var_value > max:
583                              max = var_value
584
585                     if is_used is True:
586                        tree_node._minimums[variable.name][index] = min
587                        tree_node._averages[variable.name][index] = avg / node_probability
588                        tree_node._maximums[variable.name][index] = max                       
589
590                        # distribute the newly computed average to the xbar variable in
591                        # each instance/scenario associated with this node. only do this
592                        # if the variable is used!
593                        for scenario in tree_node._scenarios:
594                           instance = self._instances[scenario._name]
595                           avg_parameter = getattr(instance, avg_parameter_name)
596                           avg_parameter[index] = avg / node_probability
597
598      end_time = time.time()
599      self._cumulative_xbar_time += (end_time - start_time)
600
601   def update_weights(self):
602
603      # because the weight updates rely on the xbars, and the xbars are node-based,
604      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
605      start_time = time.time()     
606
607      for stage in self._scenario_tree._stages:
608         
609         if stage != self._scenario_tree._stages[-1]: # no blending over the final stage, so no weights to worry about.
610           
611            for tree_node in stage._tree_nodes:
612
613               for (variable, index_template, variable_indices) in stage._variables:
614
615                  variable_name = variable.name
616                     
617                  weight_parameter_name = "PHWEIGHT_"+variable_name
618
619                  for index in variable_indices:
620
621                     tree_node_average = tree_node._averages[variable.name][index]()
622
623                     for scenario in tree_node._scenarios:
624
625                        instance = self._instances[scenario._name]
626                       
627                        if getattr(instance,variable.name)[index].status != VarStatus.unused:
628
629                           current_variable_weight = getattr(instance, weight_parameter_name)[index].value
630                           # if I'm maximizing, invert value prior to adding (hack to implement negatives)
631                           if self._is_minimizing is False:
632                              current_variable_weight = (-current_variable_weight)
633                           current_variable_value = getattr(instance,variable.name)[index].value
634                           new_variable_weight = current_variable_weight + self._rho * (current_variable_value - tree_node_average)
635                           # I have the correct updated value, so now invert if maximizing.
636                           if self._is_minimizing is False:
637                              new_variable_weight = (-new_variable_weight)
638                           getattr(instance, weight_parameter_name)[index].value = new_variable_weight
639
640      # we shouldn't have to re-simplify the expression, as we aren't adding any constant-variable terms - just modifying parameters.
641
642      end_time = time.time()
643      self._cumulative_weight_time += (end_time - start_time)
644
645   def update_scenario_weights(self, instance):
646
647      # because the weight updates rely on the xbars, and the xbars are node-based,
648      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
649      start_time = time.time()     
650
651      for stage in self._scenario_tree._stages:
652         
653         if stage != self._scenario_tree._stages[-1]: # no blending over the final stage, so no weights to worry about.
654           
655            for tree_node in stage._tree_nodes:
656
657               for (variable, index_template, variable_indices) in stage._variables:
658
659                  variable_name = variable.name
660                     
661                  weight_parameter_name = "PHWEIGHT_"+variable_name
662
663                  for index in variable_indices:
664
665                     tree_node_average = tree_node._averages[variable.name][index]()
666
667                     if getattr(instance,variable.name)[index].status != VarStatus.unused:
668
669                        current_variable_weight = getattr(instance, weight_parameter_name)[index].value
670                        # if I'm maximizing, invert value prior to adding (hack to implement negatives)
671                        if self._is_minimizing is False:
672                           current_variable_weight = (-current_variable_weight)
673                        current_variable_value = getattr(instance,variable.name)[index].value
674                        new_variable_weight = current_variable_weight + self._rho * (current_variable_value - tree_node_average)
675                        # I have the correct updated value, so now invert if maximizing.
676                        if self._is_minimizing is False:
677                           new_variable_weight = (-new_variable_weight)
678                        getattr(instance, weight_parameter_name)[index].value = new_variable_weight
679
680      # we shouldn't have to re-simplify the expression, as we aren't adding any constant-variable terms - just modifying parameters.
681
682      end_time = time.time()
683      self._cumulative_weight_time += (end_time - start_time)     
684
685   def form_iteration_k_objectives(self):
686
687      if self._verbose == True:
688         print "Forming PH-weighted objective"
689
690      # for each blended variable (i.e., those not appearing in the final stage),
691      # add the linear and quadratic penalty terms to the objective.
692      for instance_name, instance in self._instances.items():
693         
694         objective_name = instance.active_components(Objective).keys()[0]         
695         objective = instance.active_components(Objective)[objective_name]         
696         objective_expression = objective._data[None].expr # TBD: we don't deal with indexed expressions (not even sure how to interpret them)
697         # the quadratic expression is really treated as just a list - eventually should be treated as a full expression.
698         quad_expression = 0.0
699         
700         for stage in self._scenario_tree._stages:
701
702            # skip the last stage, as no blending occurs
703            if stage != self._scenario_tree._stages[-1]:
704               # find the instance objective.
705               # TBD - for simiplicity, I'm assuming a single objective - we should select "active" objectives later.
706               #       also, can create a quadratic objective and leave the original alone.
707
708               for (reference_variable, index_template, variable_indices) in stage._variables:
709
710                  variable_name = reference_variable.name
711
712                  w_parameter_name = "PHWEIGHT_"+variable_name
713                  w_parameter = instance.active_components(Param)[w_parameter_name]
714                 
715                  average_parameter_name = "PHAVG_"+variable_name
716                  average_parameter = instance.active_components(Param)[average_parameter_name]
717
718                  rho_parameter_name = "PHRHO_"+variable_name
719                  rho_parameter = instance.active_components(Param)[rho_parameter_name]
720
721                  for index in variable_indices:
722
723                     instance_variable = instance.active_components(Var)[variable_name][index]
724                     
725                     if (instance_variable.status is not VarStatus.unused) and (instance_variable.fixed is False):
726                       
727                        # TBD - if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes.
728                        objective_expression += (w_parameter[index] * instance_variable)
729                        quad_expression += (rho_parameter[index] * (instance_variable - average_parameter[index]) ** 2)
730
731         # strictly speaking, this probably isn't necessary - parameter coefficients won't get
732         # pre-processed out of the expression tree. however, if the under-the-hood should change,
733         # we'll be covered.
734         objective_expression.simplify(instance)
735         instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression
736         instance.active_components(Objective)[objective_name]._quad_subexpr = quad_expression
737
738   def iteration_k_plus_solves(self):
739
740     if self._verbose == True:
741        print "Starting PH iteration k+ solves"
742
743     # things progress at different rates - keep track of what's going on.
744     scenario_ks = {}
745     for scenario in self._scenario_tree._scenarios:
746        scenario_ks[scenario._name]=0
747
748     # keep track of action handles mapping to scenarios.
749     action_handle_instance_map = {}       
750
751     # STEP 1: queue up the solves for all scenario sub-problems.
752
753     for scenario in self._scenario_tree._scenarios:     
754
755        instance = self._instances[scenario._name]
756
757        if self._verbose == True:
758           print "Queuing solve for scenario=" + scenario._name
759
760        # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
761        # don't do this, you won't see any chance in the output files as you vary the problem parameters!
762        # ditto for instance fixing!
763        instance.presolve()
764
765        # once past iteration 0, there is always a feasible solution from which to warm-start.
766        new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
767
768        action_handle_instance_map[new_action_handle] = scenario._name
769
770     # STEP 2: wait for the first action handle to return, process it, and keep chugging.
771
772     while(True):
773
774        this_action_handle = self._solver_manager.wait_any()
775        scenario_name = action_handle_instance_map[this_action_handle]
776
777        scenario_ks[scenario_name] += 1
778
779        print "SOLVE FOR SCENARIO=",scenario_name," COMPLETED - NEW SOLVE COUNT=", scenario_ks[scenario_name]       
780
781        instance = self._instances[scenario_name]
782        results = self._solver_manager.get_results(this_action_handle)
783       
784        if len(results.solution) == 0:
785           raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
786
787        if self._verbose is True:         
788           print "Solve completed successfully"
789
790        if self._output_solver_results == True:
791           print "Results:"
792           results.write(num=1)           
793
794        instance.load(results)
795
796        if self._verbose == True:                 
797           print "Successfully loaded solution"
798
799        # we're assuming there is a single solution.
800        # the "value" attribute is a pre-defined feature of any solution - it is relative to whatever
801        # objective was selected during optimization, which of course should be the PH objective.
802        ph_objective_value = float(results.solution(0).value)
803
804        if self._verbose == True:
805           for objective_name in instance.active_components(Objective):
806              objective = instance.active_components(Objective)[objective_name]
807              print "%20s       %18.4f     %14.4f" % (scenario_name, ph_objective_value, 0.0)
808
809        # update variable statistics prior to any output.
810        self.update_variable_statistics()
811        self.update_scenario_weights(instance)
812
813        print "Variable averages and weights prior to scenario solves:"
814        self.pprint(True,True,False)       
815
816        # see if we've converged.
817        all_good = True
818        for scenario in self._scenario_tree._scenarios:     
819           instance = self._instances[scenario._name]
820           if scenario_ks[scenario._name] < self._max_iterations:
821              all_good = False
822              break
823        if all_good is True:
824           return
825
826        if scenario_ks[scenario_name] < self._max_iterations:
827
828           # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
829           # don't do this, you won't see any chance in the output files as you vary the problem parameters!
830           # ditto for instance fixing!
831           instance.presolve()
832
833           # once past iteration 0, there is always a feasible solution from which to warm-start.
834           new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
835
836           action_handle_instance_map[new_action_handle] = scenario_name
837
838           print "Queued up solve k=",str(scenario_ks[scenario_name]+1)," for scenario=",scenario_name
839
840        print "Variable values following scenario solve:"
841        self.pprint(False,False,True)           
842
843   def solve(self):
844
845      self._solve_start_time = time.time()
846      self._cumulative_solve_time = 0.0
847      self._cumulative_xbar_time = 0.0
848      self._cumulative_weight_time = 0.0
849
850      print "Starting PH"
851
852      if self._initialized == False:
853         raise RuntimeError, "PH is not initialized - cannot invoke solve() method"
854
855      print "Initiating PH iteration=" + `self._current_iteration`         
856
857      self.iteration_0_solve()
858
859      # update variable statistics prior to any output.
860      self.update_variable_statistics()
861
862      if self._verbose == True:
863         print "Variable values following scenario solves:"
864         self.pprint(False,False,True)
865
866      # let plugins know if they care.
867      if len(self._ph_plugin) == 1:
868         self._ph_plugin.service().post_iteration_0_solves(self)     
869
870      self._converger.update(self._current_iteration, self._scenario_tree, self._instances)
871      print "Convergence metric=%12.4f" % self._converger.lastMetric()
872
873      self.update_weights()
874
875      if self._max_iterations > 0:
876         self.form_iteration_k_objectives()
877
878      # let plugins know if they care.
879      if len(self._ph_plugin) == 1:
880         self._ph_plugin.service().post_iteration_0(self)
881
882      # iteration 0 is done, plugins are good to go. enter the async solve loop.
883      self.iteration_k_plus_solves()
884     
885      print "PH complete"
886
887      print "Final variable values:"
888      self.pprint(False,False,True)         
889
890      print "Final costs:"
891      self._scenario_tree.pprintCosts(self._instances)
892
893      self._solve_end_time = time.time()
894
895      if (self._verbose is True) and (self._output_times is True):
896         print "Overall solve time=" + str(self._solve_end_time - self._solve_start_time) + " seconds"
897
898      # let plugins know if they care.
899      if len(self._ph_plugin) == 1:
900         self._ph_plugin.service().post_ph_execution(self)                                 
901
902   #
903   # prints a summary of all collected time statistics
904   #
905   def print_time_stats(self):
906
907      print "PH run-time statistics (user):"
908
909      print "Initialization time=  %8.2f seconds" % (self._init_end_time - self._init_start_time)
910      print "Overall solve time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
911      print "Scenario solve time=  %8.2f seconds" % self._cumulative_solve_time
912      print "Average update time=  %8.2f seconds" % self._cumulative_xbar_time
913      print "Weight update time=   %8.2f seconds" % self._cumulative_weight_time
914
915   #
916   # a utility to determine whether to output weight / average / etc. information for
917   # a variable/node combination. when the printing is moved into a callback/plugin,
918   # this routine will go there. for now, we don't dive down into the node resolution -
919   # just the variable/stage.
920   #
921   def should_print(self, stage, variable, variable_indices):
922
923      if self._output_continuous_variable_stats is False:
924
925         variable_type = variable.domain         
926
927         if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):
928
929            return False
930
931      return True
932     
933   #
934   # pretty-prints the state of the current variable averages, weights, and values.
935   # inputs are booleans indicating which components should be output.
936   #
937   def pprint(self, output_averages, output_weights, output_values):
938
939      # TBD - write a utility routine to figure out the longest identifier width for indicies and other names,
940      #       to make the output a bit more readable.
941
942      if self._initialized is False:
943         raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"         
944     
945      # print tree nodes and associated variable/xbar/ph information in stage-order
946      for stage in self._scenario_tree._stages:
947
948         # we don't blend in the last stage, so we don't current care about printing the associated information.
949         if stage != self._scenario_tree._stages[-1]:
950
951            print "\tStage=" + stage._name
952
953            num_outputs_this_stage = 0 # tracks the number of outputs on a per-index basis.
954
955            for (variable, index_template, variable_indices) in stage._variables:
956
957               variable_name = variable.name
958
959               if self.should_print(stage, variable, variable_indices) is True:
960
961                  num_outputs_this_variable = 0 # track, so we don't output the variable names unless there is an entry to report.
962
963                  for index in variable_indices:               
964
965                     weight_parameter_name = "PHWEIGHT_"+variable_name
966
967                     num_outputs_this_index = 0 # track, so we don't output the variable index more than once.
968
969                     for tree_node in stage._tree_nodes:                 
970
971                        # determine if the variable/index pair is used across the set of scenarios (technically,
972                        # it should be good enough to check one scenario). will be deleted once we have
973                        # implemented a "cull" method to get rid of unused variables and variable indices.
974
975                        is_used = True # should be consistent across scenarios, so one "unused" flags as invalid.
976
977                        for scenario in tree_node._scenarios:
978                           instance = self._instances[scenario._name]
979                           if getattr(instance,variable_name)[index].status == VarStatus.unused:
980                              is_used = False
981 
982                        # IMPT: this is far from obvious, but variables that are fixed will - because
983                        #       presolve will identify them as constants and eliminate them from all
984                        #       expressions - be flagged as "unused" and therefore not output.
985   
986                        if is_used is True:
987
988                              minimum_value = tree_node._minimums[variable_name][index]()
989                              maximum_value = tree_node._maximums[variable_name][index]()
990
991                              num_outputs_this_stage = num_outputs_this_stage + 1                           
992                              num_outputs_this_variable = num_outputs_this_variable + 1
993                              num_outputs_this_index = num_outputs_this_index + 1
994
995                              if num_outputs_this_variable == 1:
996                                 print "\t\tVariable=",variable_name
997
998                              if num_outputs_this_index == 1:
999                                 print "\t\t\tIndex:", indexToString(index)                             
1000
1001                              print "\t\t\t\tTree Node=",tree_node._name,"\t\t (Scenarios: ",                             
1002                              for scenario in tree_node._scenarios:
1003                                 print scenario._name," ",
1004                                 if scenario == tree_node._scenarios[-1]:
1005                                    print ")"
1006                           
1007                              if output_values is True:
1008                                 average_value = tree_node._averages[variable_name][index]()
1009                                 print "\t\t\t\tValues: ",                       
1010                                 for scenario in tree_node._scenarios:
1011                                    instance = self._instances[scenario._name]
1012                                    this_value = getattr(instance,variable_name)[index].value
1013                                    print "%12.4f" % this_value,
1014                                    if scenario == tree_node._scenarios[-1]:
1015                                       print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1016                                       print "    Avg=%12.4f" % (average_value),
1017                                       print ""
1018                              if output_weights:
1019                                 print "\t\t\t\tWeights: ",
1020                                 for scenario in tree_node._scenarios:
1021                                    instance = self._instances[scenario._name]
1022                                    print "%12.4f" % getattr(instance,weight_parameter_name)[index].value,
1023                                    if scenario == tree_node._scenarios[-1]:
1024                                       print ""
1025                              if output_averages:
1026                                 print "\t\t\t\tAverage: %12.4f" % (tree_node._averages[variable_name][index].value)
1027
1028            if num_outputs_this_stage == 0:
1029               print "\t\tNo non-converged variables in stage"
1030
1031            # cost variables aren't blended, so go through the gory computation of min/max/avg.
1032            # TBD: possibly move these into the tree node anyway, as they are useful in general.
1033            # TBD: these should also be probability-weighted - fix with above-mentioned TBD fix.
1034            # we currently always print these.
1035            cost_variable_name = stage._cost_variable[0].name
1036            cost_variable_index = stage._cost_variable[1]
1037            if cost_variable_index is None:
1038               print "\t\tCost Variable=" + cost_variable_name
1039            else:
1040               print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)           
1041            for tree_node in stage._tree_nodes:
1042               print "\t\t\tTree Node=" + tree_node._name + "\t\t (Scenarios: ",
1043               for scenario in tree_node._scenarios:
1044                  print scenario._name," ",
1045                  if scenario == tree_node._scenarios[-1]:
1046                     print ")"
1047               maximum_value = 0.0
1048               minimum_value = 0.0
1049               sum_values = 0.0
1050               num_values = 0
1051               first_time = True
1052               print "\t\t\tValues: ",                       
1053               for scenario in tree_node._scenarios:
1054                   instance = self._instances[scenario._name]
1055                   this_value = getattr(instance,cost_variable_name)[cost_variable_index].value
1056                   print "%12.4f" % this_value,
1057                   num_values += 1
1058                   sum_values += this_value
1059                   if first_time is True:
1060                      first_time = False
1061                      maximum_value = this_value
1062                      minimum_value = this_value
1063                   else:
1064                      if this_value > maximum_value:
1065                         maximum_value = this_value
1066                      if this_value < minimum_value:
1067                         minimum_value = this_value
1068                   if scenario == tree_node._scenarios[-1]:
1069                      print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1070                      print "    Avg=%12.4f" % (sum_values/num_values),
1071                      print ""
1072           
Note: See TracBrowser for help on using the repository browser.