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

Last change on this file since 3261 was 3073, checked in by jwatson, 10 years ago

Complete set of changes to make PySP compatible with latest immutable parameters change.

File size: 49.4 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.component.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.component.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(scenario_data)
319               scenario_instance.preprocess()
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, mutable=True)
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, mutable=True)
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, mutable=True)
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, mutable=True)
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, mutable=True)
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, mutable=True)
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 = value(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 = value(getattr(instance, weight_parameter_name)[index])
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] = 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 = value(tree_node._averages[variable.name][index])
666
667                     if getattr(instance,variable.name)[index].status != VarStatus.unused:
668
669                        current_variable_weight = value(getattr(instance, weight_parameter_name)[index])
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] = 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.preprocess()
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.preprocess()
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 = value(tree_node._minimums[variable_name][index])
989                              maximum_value = 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 = 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" % value(getattr(instance,weight_parameter_name)[index])
1023                                    if scenario == tree_node._scenarios[-1]:
1024                                       print ""
1025                              if output_averages:
1026                                 print "\t\t\t\tAverage: %12.4f" % (value(tree_node._averages[variable_name][index]))
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.