source: trunk/examples/pysp/forestry/testphextension.py @ 1475

Revision 1475, 13.5 KB checked in by jwatson, 5 years ago (diff)

Various fixes to fix-lag extension in pysp examples.

  • Property svn:executable set to *
Line 
1# ver 0.0106
2#  _________________________________________________________________________
3#
4#  Coopr: A COmmon Optimization Python Repository
5#  Copyright (c) 2009 Sandia Corporation.
6#  Self software is distributed under the BSD License.
7#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8#  the U.S. Government retains certain rights in self software.
9#  For more information, see the Coopr README.txt file.
10#  _________________________________________________________________________
11
12from pyutilib.plugin import *
13from coopr.pysp import phextension
14from coopr.pyomo.base import *
15
16from coopr.pyomo.base.set_types import *
17
18from coopr.pyomo.base.var import VarStatus
19
20import string
21
22class TestPHExtension(SingletonPlugin):
23
24   implements (phextension.IPHExtension)
25
26   def post_ph_initialization(self, ph):
27
28      print "****Called after PH initialization!"
29
30      # create a set containing all PH iteration indices - 0 through ph._max_iterations
31      # the index set should (I think) be stored persisently somewhere - so I put it into
32      # the ph class for lack of a better place.
33      ph._iteration_index_set = Set(name="PHIterations")
34      for i in range(0,ph._max_iterations + 1):
35         ph._iteration_index_set.add(i)
36
37      # set up "global" record keeping
38      self.cumulative_fixed_by_consensus_count = 0
39
40      # constants for W vector hashing (low cost rand() is all we need)
41      # note: July 09, dlw is planning to eschew pre-computed random vector
42      self.W_hash_history_len = 10
43      self.W_hash_seed = 17  # we will reset for dynamic rand vector generation
44      self.W_hash_rand_val = 0; # we will use self to avoid a function call
45      self.W_hash_a = 1317       # a,b, and c are for a simple generator
46      self.W_hash_b = 27699
47      self.W_hash_c = 131072  # that period should be usually long enough for us!
48                              # (assuming fewer than c scenarios)
49
50      # set up tree storage
51      for stage in ph._scenario_tree._stages:
52         
53         if stage != ph._scenario_tree._stages[-1]:
54
55            # first, gather all unique variables referenced in self stage
56            # self "gather" step is currently required because we're being lazy
57            # in terms of index management in the scenario tree
58            stage_variables = {}
59            for (reference_variable, reference_index) in stage._variables:
60               if reference_variable.name not in stage_variables.keys():
61                  stage_variables[reference_variable.name] = reference_variable
62
63            for tree_node in stage._tree_nodes:
64               tree_node._num_iters_converged = {}
65               tree_node._W_hash = {}
66
67            # next, create parameters for each variable in the corresponding tree node.
68            # NOTE: the parameter names below could really be empty, as they are never referenced
69            #       explicitly.
70            # TBD: can we make them only for the booleans and integers? Should we?
71            for (variable_name, reference_variable) in stage_variables.items():
72               for tree_node in stage._tree_nodes:
73                  new_stat_index = reference_variable._index # TBD - need to be careful with the shallow copy
74                  new_stat_parameter_name = "NODESTAT_NUM_ITERS_CONVERGED_"+reference_variable.name
75                  new_stat_parameter = Param(new_stat_index,name=new_stat_parameter_name)
76                  for index in new_stat_index:
77                     new_stat_parameter[index] = 0
78                  tree_node._num_iters_converged[reference_variable.name] = new_stat_parameter
79                  # now make the w hash value storage array
80                  # TBD: it needs to be two dimensional
81                  new_hash_index = reference_variable._index # TBD - need to be careful with the shallow copy
82                  new_hash_parameter_name = "W_HASH_STORAGE_"+reference_variable.name
83                  new_hash_parameter = Param(new_hash_index, ph._iteration_index_set, name=new_hash_parameter_name)
84                  for index in new_hash_index:    # do a need a new copy?
85                     for i in range(1, self.W_hash_history_len+1):
86                        new_hash_parameter[index,i] = 0
87                  tree_node._W_hash[reference_variable.name] = new_hash_parameter
88
89   def post_iteration_0_solves(self, ph):
90
91      print "ITERATION 0 SOLVES ARE DONE - TOTAL CPU TIME="+str(ph._cumulative_solve_time)
92
93      for stage in ph._scenario_tree._stages:
94         
95         if stage != ph._scenario_tree._stages[-1]: # no blending over the final stage
96           
97            for tree_node in stage._tree_nodes:
98
99               for (variable, idx) in stage._variables:
100                 
101                  variable_name = variable.name
102                  variable_index = variable._index
103                  # The index, if not None, specifies a single index.
104                  if idx != None:
105                     variable_index = []
106                     variable_index.append(idx)
107
108                  variable_type = variable.domain
109
110#                  print "var NAME=",variable_name," TYPE=",type(variable_type)
111                  if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
112
113                     for index in variable_index:
114
115                        # THE LB/UB REQUIRE A VARIABLE INDEX; THEY ARE ATTRIBUTES OF A VARIABLE VALUE.
116                        variable_lb = variable[index].lb
117                        variable_ub = variable[index].ub                                             
118
119                        node_min = tree_node._minimums[variable_name][index]()
120                        node_max = tree_node._maximums[variable_name][index]()
121
122                        should_fix = False
123
124                        # keep some records
125
126                        # keep track of cummulative iters of convernce
127                        if (node_min == node_max):
128                           ++tree_node._num_iters_converged[variable_name][index]
129                           
130                        # hash the W vectors to enable a test for cycling
131                        if self.W_hash_history_len > 0:
132                           self.W_hash_rand_val = self.W_hash_seed  # reset seed
133                           for scenario in tree_node._scenarios:
134                              # TBD: python has one-based indexes (self is iter 0)
135                              instance = ph._instances[scenario._name]
136                              weight_parameter_name = "PHWEIGHT_"+variable_name
137                              tree_node._W_hash[variable_name][index,1] += getattr(instance, weight_parameter_name)[index].value * self.W_hash_rand_val
138                              self.W_hash_rand_val = (self.W_hash_b + self.W_hash_a * self.W_hash_rand_val) % self.W_hash_c
139                           
140                        # DLW TBD: too simple
141                        # for one thing, test a suffix to see if fixing is allowed (similar to slamming)
142                        #   and if so, to what (lb, ub, any, none)
143                        if tree_node._num_iters_converged[variable_name][index]() > 0:
144                           should_fix = True
145                           fix_value = node_min
146                        # DLW TBD end too simple   
147
148                        if should_fix is True:
149                           for scenario in tree_node._scenarios:
150                           
151                              instance = ph._instances[scenario._name]
152                       
153                              if getattr(instance,variable.name)[index].status != VarStatus.unused:
154                                 print "debug dlw: fixing ", variable.name, " ", index
155                                 getattr(instance,variable.name)[index].value = fix_value
156                                 getattr(instance,variable.name)[index].fixed = True
157                                 ++self.cumulative_fixed_by_consensus_count
158                         
159
160   def post_iteration_0(self, ph):
161
162      print "ITERATION 0 IS DONE - TOTAL CPU TIME="+str(ph._cumulative_solve_time)
163
164   def post_iteration_k_solves(self, ph):
165
166      print "TEST PH EXTENSION CALLBACK FOR ITERATION K SOLVES!"
167
168      for stage in ph._scenario_tree._stages:
169         
170         if stage != ph._scenario_tree._stages[-1]: # no blending over the final stage
171           
172            for tree_node in stage._tree_nodes:
173
174               for (variable, idx) in stage._variables:
175
176                  variable_name = variable.name
177                  variable_index = variable._index
178                  # The index, if not None, specifies a single index.
179                  if idx != None:
180                     variable_index = []
181                     variable_index.append(idx)
182
183                  variable_type = variable.domain
184
185                  # JPW TBD: I doubt this is ultimately sufficient - we really care whether the domain
186                  #          is 'discrete', but Pyomo doesn't have that notion abstracted at this point.
187                  # JPW TBD: We don't fix continuous variables at the moment, but we could with a
188                  #          threshold value for equality.
189                  if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
190
191                     for index in variable_index:
192
193                        # NOTE: variable statistics are updated immediately following the iteration k solves
194                        #       in the main PH code, and just prior to invoking this callback.
195                        node_min = tree_node._minimums[variable_name][index]()
196                        node_max = tree_node._maximums[variable_name][index]()
197
198                        # keep track of cumulative iterations converged.
199                        # JPW TBD: Should be at the same value as well (secondary check).
200                        # JPW TBD: Parameterize the tolerance on comparison in the case of discrete values.
201                        if abs(node_max - node_min) <= 0.00001:
202                           tree_node._num_iters_converged[variable_name][index].value = tree_node._num_iters_converged[variable_name][index].value + 1
203
204                        # determine if the variable is already fixed.
205                        instance_fixed_count = 0
206                        for scenario in tree_node._scenarios:
207                           instance = ph._instances[scenario._name]
208                           if getattr(instance,variable_name)[index].fixed is True:
209                              ++instance_fixed_count
210                        if ((instance_fixed_count > 0) and (instance_fixed_count < len(tree_node._scenarios))):
211                           raise RuntimeError, "Variable="+variable_name+str(index)+" is fixed in "+str(instance_fixed_count)+" scenarios, but less than the number of scenarios at tree node="+tree_node._name
212
213                        if instance_fixed_count == 0:
214
215                           # until proven otherwise.
216                           should_fix = False
217
218                           # DLW TBD: too simple
219                           # for one thing, test a suffix to see if fixing is allowed (similar to slamming)
220                           #   and if so, to what (lb, ub, any, none)
221                           if tree_node._num_iters_converged[variable_name][index]() >= 1:
222                             
223                              should_fix = True
224                              # whether you fix at current values or not can severly impact feasibility later
225                              # in the game. my original thought was below - but that didn't work out. if you
226                              # make integers, well, integers, all appears to be well.
227                              # IMPT: we are leaving the values for individual variables alone, at potentially
228                              #       smallish and heterogeneous values. if we fix/round, we run the risk of
229                              #       infeasibilities due to numerical issues. the computed value below is
230                              #       strictly for output purposes.
231                              fix_value = int(round(node_min))
232
233                           if should_fix is True:
234
235                              fixing_reported = False                           
236
237                              for scenario in tree_node._scenarios:
238
239                                 instance = ph._instances[scenario._name]
240                       
241                                 if getattr(instance,variable.name)[index].status != VarStatus.unused:
242
243                                    getattr(instance,variable.name)[index].fixed = True
244                                    getattr(instance,variable.name)[index].value = fix_value
245
246                                    # JPW TBD: make the output switchable (add verbose flag)
247                                    if fixing_reported is False:
248                                       # pretty-print the index, string the trailing spaces from the strings.
249                                       # JPW TBD: promote the indexToString function to a ph utility, and incorporate below.
250                                       print "Fixing variable=",variable.name,index," at tree node=",tree_node._name,", stage=",stage._name," to value=",fix_value,"; converged for ", tree_node._num_iters_converged[variable_name][index](), " iterations"
251                                       fixing_reported = True
252                                       self.cumulative_fixed_by_consensus_count = self.cumulative_fixed_by_consensus_count + 1
253
254
255
256
257      print "Total number of variables fixed by consensus=",self.cumulative_fixed_by_consensus_count
258
259   def post_iteration_k(self, ph):
260       print" Callback for iteration k"
261
262
Note: See TracBrowser for help on using the repository browser.