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

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

Some prep work for pickle'ing a PH object. No functional change quite yet, but preliminary suggests indicate it actually works!

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