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

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

Various fixes to the WW PH extension, bringing it in compliance to previous commit changes.
  • Property svn:executable set to *
File size: 49.7 KB
Line 
1#  _________________________________________________________________________
2#
3#  Coopr: A COmmon Optimization Python Repository
4#  Copyright (c) 2009 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 types
12from pyutilib.component.core import *
13from coopr.pysp import phextension
14from coopr.pysp.phutils import *
15from coopr.pyomo.base import *
16
17from coopr.pyomo.base.set_types import *
18
19from coopr.pyomo.base.var import VarStatus
20
21import string
22import os
23
24#=========================
25def slam_priority_descend_compare(a, b):
26   # used to sort the variable-suffix map for slamming priority
27   value_a = getattr(a, "SlammingPriority")
28   value_b = getattr(b, "SlammingPriority")
29   return cmp(value_b, value_a)
30
31#==================================================
32class wwphextension(SingletonPlugin):
33
34   implements (phextension.IPHExtension)
35
36   def __init__(self, *args, **kwds):
37
38      # TBD - migrate all of the self attributes defined on-the-fly
39      #       in the post-post-initialization routine here!
40
41      self._configuration_filename = None 
42      self._suffix_filename = None
43
44#=========================
45   def process_suffix_file(self, ph):
46
47      # for each suffix, maintain an inverse map from the suffix name to a list of
48      # objects with that suffix. an object could be a variable (_VarBase) or a
49      # variable value, depending on the resolution of the object for which the
50      # suffix is defined.
51      self._suffix_to_variable_map = {}
52
53      self.slam_list = []
54     
55      if self._suffix_filename is None:
56         return
57
58      if os.path.exists(self._suffix_filename) is False:
59         raise RuntimeError, "***WW PH Extension: The suffix file "+self._suffix_filename+" either does not exist or cannot be read"
60
61      print "WW PH Extension: Loading variable suffixes from file=", self._suffix_filename
62
63      reference_instance = ph._model_instance
64
65      suffix_file = open(self._suffix_filename,'r')
66      for line in suffix_file.readlines():
67         line = line.strip()
68         if len(line) > 0 and line[0] != '#':
69            pieces = line.split()
70            if len(pieces) != 3:
71               raise RuntimeError, "Illegal line=["+line+"] encountered in ww ph extension suffix file="+self._suffix_filename+"; format is variable, suffix, suffix-value."
72
73            variable_string = pieces[0]
74            suffix_name = pieces[1]
75            suffix_value = pieces[2]
76
77            # decide what type of suffix value we're dealing with.
78            is_int = False
79            is_bool = False
80            converted_value = None
81            try:
82               converted_value = bool(suffix_value)
83               is_bool = True
84            except valueerror:
85               pass
86            try:
87               converted_value = int(suffix_value)
88               is_int = True
89            except ValueError:
90               pass
91
92            if (is_int is False) and (is_bool is False):
93               raise RuntimeError, "WW PH Extension unable to deduce type of data referenced in ww ph extension suffix file="+self._suffix_filename+"; value="+suffix_value+" for "+variable_name
94
95            # determine if we're dealing with complete variables or indexed variables.
96            if isVariableNameIndexed(variable_string) is True:
97
98               variable_name, index_template = extractVariableNameAndIndex(variable_string)
99               
100               # verify that the root variable exists and grab it.
101               if variable_name not in reference_instance.active_components(Var).keys():
102                  raise RuntimeError, "Unknown variable="+variable_name+" referenced in ww ph extension suffix file="+self._suffix_filename
103               variable = reference_instance.active_components(Var)[variable_name]
104
105               # extract all "real", i.e., fully specified, indices matching the index template.
106               match_indices = extractVariableIndices(variable, index_template)
107
108               # there is a possibility that no indices match the input template.
109               # if so, let the user know about it.
110               if len(match_indices) == 0:
111                  raise RuntimeError, "No indices match template="+str(index_template)+" for variable="+variable_name               
112
113               # add the suffix to all variable values identified.
114               for index in match_indices:
115
116                  variable_value = variable[index]
117
118                  # if the suffix already exists, we don't care - we're stomping it.
119                  if hasattr(variable_value, suffix_name) is True:
120                     delattr(variable_value, suffix_name)
121 
122                  # set the suffix on the variable value.
123                  setattr(variable_value, suffix_name, converted_value)
124
125                  # place the variable value in the suffix->variable map, for easy searching elsewhere in this plugin.
126                  if suffix_name not in self._suffix_to_variable_map.keys():
127                     self._suffix_to_variable_map[suffix_name] = []
128                  self._suffix_to_variable_map[suffix_name].append(variable_value)                 
129
130            else:
131
132               # verify that the variable exists.
133               if variable_string not in reference_instance.active_components(Var).keys():
134                  raise RuntimeError, "Unknown variable="+variable_string+" referenced in ww ph extension suffix file="+self._suffix_filename
135
136               variable = reference_instance.active_components(Var)[variable_string]
137
138               # 9/14/2009 - now forcing the user to explicit specify the full
139               # match template (e.g., "foo[*,*]") instead of just the variable
140               # name (e.g., "foo") to represent the set of all indices.
141               
142               # if the variable is a singleton - that is, non-indexed - no brackets is fine.
143               # we'll just tag the var[None] variable value with the (suffix,value) pair.
144               if None not in variable._index:
145                  raise RuntimeError, "Illegal match template="+variable_string+" specified in ww ph extension suffix file="+self._suffix_filename                 
146
147               # if the suffix already exists, we don't care - we're stomping it.
148               if hasattr(variable, suffix_name) is True:
149                  delattr(variable, suffix_name)
150 
151               # set the suffix on the variable.
152               setattr(variable, suffix_name, converted_value)
153
154               # place the variable in the suffix->variable map, for easy searching elsewhere in this plugin.
155               if suffix_name not in self._suffix_to_variable_map.keys():
156                  self._suffix_to_variable_map[suffix_name] = []
157               self._suffix_to_variable_map[suffix_name].append(variable)
158
159#      print "pre-sort suffix->variable map:"
160#      for key,value_list in self._suffix_to_variable_map.items():
161#         print "key=",key,":",
162#         for value in value_list:
163#            print value.name,"",
164#         print ""
165
166      if "SlammingPriority" in self._suffix_to_variable_map:
167         self.slam_list = self._suffix_to_variable_map["SlammingPriority"]
168
169      self.slam_list.sort(slam_priority_descend_compare)
170
171#==================================================
172   def post_instance_creation(self, ph):
173       """ Called after PH initialization has created the scenario instances, but before any PH-related weights/variables/parameters/etc are defined!"""
174       # we don't muck with the instances.
175       pass
176
177#==================================================
178   def post_ph_initialization(self, ph):
179
180      # set up "global" record keeping.
181      self.cumulative_discrete_fixed_count = 0
182      self.cumulative_continuous_fixed_count = 0
183
184      # we always track convergence of continuous variables, but we may not want to fix/slam them.                                                           
185      self.fix_continuous_variables = False
186
187      # there are occasions where we might want to fix any values at the end of the
188      # run if the scenarios agree - even if the normal fixing criterion (e.g.,
189      # converged for N iterations) don't apply. one example is when the term-diff
190      # is 0, at which point you really do have a solution. currently only applies
191      # to discrete variables.
192      self.fix_discrete_variables_at_exit = False
193
194      # set up the mipgap parameters (zero means ignore)
195      # note: because we lag one iteration, the final will be less than requested
196      # initial and final refer to PH iteration 1 and PH iteration X, where
197      # X is the iteration at which the convergence metric hits 0.
198      self.Iteration0MipGap = 0.0
199      self.InitialMipGap = 0.10
200      self.FinalMipGap = 0.001
201
202      # "Read" the defaults for parameters that control fixing
203      # (these defaults can be overridden at the variable level)
204      # for all six of these, zero means don't do it.
205      self.Iter0FixIfConvergedAtLB = 0 # 1 or 0
206      self.Iter0FixIfConvergedAtUB = 0  # 1 or 0
207      self.Iter0FixIfConvergedAtNB = 0  # 1 or 0 (converged to a non-bound)
208      # TBD: check for range errors for all six of these
209      self.FixWhenItersConvergedAtLB = 10 
210      self.FixWhenItersConvergedAtUB = 10
211      self.FixWhenItersConvergedAtNB = 12  # converged to a non-bound
212      self.FixWhenItersConvergedContinuous = 0
213     
214      # "default" slamming parms:
215      # TBD: These should get ovverides from suffixes
216      # notice that for a particular var, all could be False
217      # TBD: the user might also want to give slamming priorities
218      self.CanSlamToLB = False
219      self.CanSlamToMin = False
220      self.CanSlamToAnywhere = True
221      self.CanSlamToMax = False
222      self.CanSlamToUB = False
223      self.PH_Iters_Between_Cycle_Slams = 1  # zero means "slam at will"
224      self.SlamAfterIter = len(ph._scenario_tree._stages[-1]._tree_nodes)
225
226      self.CycleLengthSlamThreshold = len(ph._scenario_tree._stages[-1]._tree_nodes) 
227      self.W_hash_history_len = max(100, self.CycleLengthSlamThreshold+1)
228
229      self.ReportPotentialCycles = 0 # do I report potential cycles, i.e., those too short to base slamming on?
230
231      # end of parms
232
233      self._last_slam_iter = -1    # dynamic
234
235      # constants for W vector hashing (low cost rand() is all we need)
236      # note: July 09, dlw is planning to eschew pre-computed random vector
237      # another note: the top level reset is OK, but really it should
238      #   done way down at the index level with a serial number (stored)
239      #   to avoid correlated hash errors
240      # the hash seed was deleted in 1.1 and we seed with the
241      self.W_hash_seed = 17  # we will reset for dynamic rand vector generation
242      self.W_hash_rand_val = self.W_hash_seed # current rand
243      self.W_hash_a = 1317       # a,b, and c are for a simple generator
244      self.W_hash_b = 27699
245      self.W_hash_c = 131072  # that period should be usually long enough for us!
246                              # (assuming fewer than c scenarios)
247
248      # set up tree storage for statistics that we care about tracking.
249      for stage in ph._scenario_tree._stages[:-1]:
250
251         # first, gather all unique variables referenced in self stage
252         # self "gather" step is currently required because we're being lazy
253         # in terms of index management in the scenario tree
254##         stage_variables = {}
255##         for (reference_variable, index_template, reference_index) in stage._variables:
256##            if reference_variable.name not in stage_variables.keys():
257##               stage_variables[reference_variable.name] = reference_variable
258##
259         for tree_node in stage._tree_nodes:
260               
261            tree_node._num_iters_converged = {}
262            tree_node._last_converged_val = {}
263            tree_node._w_hash = {}
264            tree_node._fixed_var_flag = {}
265
266            for (reference_variable, index_template) in stage._variables:
267
268               variable_name = reference_variable.name
269               variable_type = reference_variable.domain
270
271               variable_indices = tree_node._variable_indices[variable_name]
272
273               new_param_index = getattr(ph._instances[tree_node._scenarios[0]._name], variable_name)._index
274
275               # create the parameters for all indices, even though this might be a bit
276               # wasteful. in practice, this probably won't be the case, and determining
277               # the reduced index set would probably take longer that it is worth.
278
279               new_stat_index = new_param_index
280               new_stat_parameter_name = "NODESTAT_NUM_ITERS_CONVERGED_"+reference_variable.name
281               new_stat_parameter = None
282               # this bit of ugliness is due to Pyomo not correctly handling the Param construction
283               # case when the supplied index set consists strictly of None, i.e., the source variable
284               # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.                     
285               if (len(new_stat_index) is 1) and (None in new_stat_index):
286                  new_stat_parameter = Param(name=new_stat_parameter_name, mutable=True)
287               else:
288                  new_stat_parameter = Param(new_stat_index, name=new_stat_parameter_name, mutable=True)
289               for newindex in new_stat_index:
290                  new_stat_parameter[newindex] = 0
291               tree_node._num_iters_converged[reference_variable.name] = new_stat_parameter
292                     
293               # need to know to what we have most recently converged
294               new_conv_index = new_param_index
295               new_conv_parameter_name = "NODESTAT_LAST_CONVERGED_VAL_"+reference_variable.name
296               new_conv_parameter = None
297               if (len(new_conv_index) is 1) and (None in new_conv_index):
298                  new_conv_parameter = Param(name=new_conv_parameter_name, mutable=True)
299               else:
300                  new_conv_parameter = Param(new_conv_index, name=new_conv_parameter_name, mutable=True)
301               for newindex in new_conv_index:
302                  new_conv_parameter[newindex] = 0.5 # not an int, so harmless
303               tree_node._last_converged_val[reference_variable.name] = new_conv_parameter
304                     
305               # need to know to what has been fixed
306               new_fix_index = new_param_index
307               new_fix_parameter_name = "NODESTAT_FIXED_FLAG_VAL_"+reference_variable.name
308               new_fix_parameter = None
309               if (len(new_fix_index) is 1) and (None in new_fix_index):
310                  new_fix_parameter = Param(name=new_fix_parameter_name, mutable=True)
311               else:
312                  new_fix_parameter = Param(new_fix_index, name=new_fix_parameter_name, mutable=True)
313               for newindex in new_fix_index:
314                  new_fix_parameter[newindex] = False
315               tree_node._fixed_var_flag[reference_variable.name] = new_fix_parameter
316
317               # now make the w hash value storage array
318               new_hash_index = new_param_index
319               new_hash_parameter_name = "W_HASH_STORAGE_"+reference_variable.name
320               new_hash_parameter = None
321               if (len(new_hash_index) is 1) and (None in new_hash_index):
322                  new_hash_parameter = Param(ph._iteration_index_set, name=new_hash_parameter_name, mutable=True)
323               else:
324                  new_hash_parameter = Param(new_hash_index, ph._iteration_index_set, name=new_hash_parameter_name, mutable=True)
325               for newindex in new_hash_index:
326                  for i in range(0, ph._max_iterations+1):
327                     # the following if-then block is a complete hack, due to the
328                     # fact that we can't index by None if the Param is unary.
329                     if new_hash_parameter.dim() == 1:
330                        new_hash_parameter[i] = 0
331                     else:
332                        new_hash_parameter[newindex,i] = 0
333               tree_node._w_hash[reference_variable.name] = new_hash_parameter
334                 
335               # JPW has no idea why the following code block is here, or if it is necessary.
336               for index in variable_indices:
337
338                  # IMPT: it has bitten us before in this plug-in, so I'll state the obvious:
339                  #       variable values in the last stage of a stochastic program will *not*
340                  #       have a defined _stage attribute.
341                  # NOTE: the use of the reference variable here, in contrast to a instance-specific
342                  #       variable, is potentially dangerous.
343                  reference_variable[index]._stage = stage
344
345
346      # store the total variable counts for future reporting.
347      (self.total_discrete_vars,self.total_continuous_vars) = ph.compute_variable_counts()
348
349      if self._configuration_filename is not None:
350         if os.path.exists(self._configuration_filename) is True:
351            print "WW PH Extension: Loading user-specified configuration from file=" + self._configuration_filename
352            execfile(self._configuration_filename)
353         else:
354            raise RuntimeError, "***WW PH Extension: The configuration file "+self._configuration_filename+" either does not exist or cannot be read"
355      else:
356         print "WW PH Extension: No user-specified configuration file supplied - using defaults"
357
358      # process any suffix information, if it exists.
359      self.process_suffix_file(ph)         
360
361      # set up the mip gap for iteration 0.
362      if self.Iteration0MipGap > 0.0:
363         print "Setting mipgap to "+str(self.Iteration0MipGap)
364         ph._mipgap = self.Iteration0MipGap
365         
366#==================================================
367   def post_iteration_0_solves(self, ph):
368
369      for stage in ph._scenario_tree._stages[:-1]:
370         
371         for tree_node in stage._tree_nodes:
372
373            for (variable, index_template) in stage._variables:
374                 
375               variable_name = variable.name
376               variable_type = variable.domain
377
378               variable_indices = tree_node._variable_indices[variable_name]
379
380               for index in variable_indices:
381
382                  # determine if this index is used - otherwise, don't waste time
383                  # fixing and cycle checking. for one, the code will crash :-) with
384                  # none values during the cycle checking computation!
385
386                  is_used = True # until proven otherwise                     
387                  for scenario in tree_node._scenarios:
388                     instance = ph._instances[scenario._name]
389                     if getattr(instance,variable_name)[index].status == VarStatus.unused:
390                        is_used = False
391
392                  if is_used is True:
393
394                     # unlikely, but variables might be fixed even at this stage, depending on
395                     # what weirdness users do prior to the iteration 0 solves.
396                     instance_fixed_count = 0
397                     for scenario in tree_node._scenarios:
398                        instance = ph._instances[scenario._name]
399                        if getattr(instance,variable_name)[index].fixed is True:
400                           instance_fixed_count += 1
401                     if ((instance_fixed_count > 0) and (instance_fixed_count < len(tree_node._scenarios))):
402                        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
403
404                     if instance_fixed_count == 0:
405
406                        if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):                           
407                           node_min = self.Int_If_Close_Enough(ph, value(tree_node._minimums[variable_name][index]))
408                           node_max = self.Int_If_Close_Enough(ph, value(tree_node._maximums[variable_name][index]))
409
410                           # update convergence prior to checking for fixing.
411                           self._int_convergence_tracking(ph, tree_node, variable_name, index, node_min, node_max)
412                           attrvariable = ph._model_instance.active_components(Var)[variable_name][index]
413                           if hasattr(attrvariable, 'Iter0FixIfConvergedAtLB'):
414                              lb = getattr(attrvariable, 'Iter0FixIfConvergedAtLB')
415                           else:
416                              lb = self.Iter0FixIfConvergedAtLB
417                           if hasattr(attrvariable, 'Iter0FixIfConvergedatUB'):
418                              ub = getattr(attrvariable, 'Iter0FixIfConvergedAtUB')
419                           else:
420                              ub = self.Iter0FixIfConvergedAtUB
421                           if hasattr(attrvariable, 'Iter0FixIfConvergedAtNB'):
422                              nb = getattr(attrvariable, 'Iter0FixIfConvergedAtNB')
423                           else:
424                              nb = self.Iter0FixIfConvergedAtNB
425                           if self._should_fix_discrete_due_to_conv(tree_node, variable, index, lb, ub, nb):
426                              self._fix_var(ph, tree_node, variable, index, node_min)
427                           elif self.W_hash_history_len > 0:   # if not fixed, then hash - no slamming at iteration 0
428                              self._w_hash_acct(ph, tree_node, variable_name, index) # obviously not checking for cycles at iteration 0!
429
430                        else:
431                             
432                           node_min = value(tree_node._minimums[variable_name][index])
433                           node_max = value(tree_node._maximums[variable_name][index])
434                             
435                           self._continuous_convergence_tracking(ph, tree_node, variable_name, index, node_min, node_max)
436                             
437# jpw: not sure if we care about cycle detection in continuous variables?
438#                           if self.W_hash_history_len > 0: 
439#                              self._w_hash_acct(ph, tree_node, variable_name, index)
440                             
441
442#==================================================
443   def post_iteration_0(self, ph):
444 
445      self._met0 = ph._converger.lastMetric()
446
447      if (self.InitialMipGap > 0) and (self.FinalMipGap >= 0) and (self.InitialMipGap > self.FinalMipGap):
448         gap = self.InitialMipGap
449         print "Setting mipgap to "+str(gap)
450         ph._mipgap = gap
451
452#==================================================         
453
454   def pre_iteration_k_solves(self, ph):
455        """ Called immediately before the iteration k solves!"""
456        # we don't muck with the PH objectives
457        pass
458
459#==================================================
460   def post_iteration_k_solves(self, ph):
461
462      for stage in ph._scenario_tree._stages[:-1]:
463         
464         for tree_node in stage._tree_nodes:
465
466            for (variable, index_template) in stage._variables:
467
468               variable_name = variable.name
469               variable_type = variable.domain
470
471               variable_indices = tree_node._variable_indices[variable_name]
472
473               for index in variable_indices:
474
475                  # determine if this index is used - otherwise, don't waste time
476                  # fixing and cycle checking. for one, the code will crash :-) with
477                  # None values during the cycle checking computation!
478
479                  is_used = True # until proven otherwise                     
480                  for scenario in tree_node._scenarios:
481                     instance = ph._instances[scenario._name]
482                     if getattr(instance,variable_name)[index].status == VarStatus.unused:
483                        is_used = False
484
485                  if is_used is True:                       
486
487                     # determine if the variable is already fixed.
488                     instance_fixed_count = 0
489                     for scenario in tree_node._scenarios:
490                        instance = ph._instances[scenario._name]
491                        if getattr(instance,variable_name)[index].fixed is True:
492                           instance_fixed_count += 1
493                     if ((instance_fixed_count > 0) and (instance_fixed_count < len(tree_node._scenarios))):
494                        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
495
496                     if instance_fixed_count == 0:
497
498                        if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):                           
499                           node_min = self.Int_If_Close_Enough(ph, value(tree_node._minimums[variable_name][index]))
500                           node_max = self.Int_If_Close_Enough(ph, value(tree_node._maximums[variable_name][index]))
501
502                           # update convergence prior to checking for fixing.
503                           self._int_convergence_tracking(ph, tree_node, variable_name, index, node_min, node_max)
504
505                           # now check on permissions to converge to various placed (e.g., lb is lb permission)
506                           attrvariable = ph._model_instance.active_components(Var)[variable_name][index]
507                           if hasattr(attrvariable, 'FixWhenItersConvergedAtLB'):
508                              lb = getattr(attrvariable, 'FixWhenItersConvergedAtLB')
509                           else:
510                              lb = self.FixWhenItersConvergedAtLB
511                           if hasattr(attrvariable, 'FixWhenItersConvergedAtUB'):
512                              ub = getattr(attrvariable, 'FixWhenItersConvergedAtUB')
513                           else:
514                              ub = self.FixWhenItersConvergedAtUB
515                           if hasattr(attrvariable, 'FixWhenItersConvergedAtNB'):
516                              nb = getattr(attrvariable, 'FixWhenItersConvergedAtNB')
517                           else:
518                              nb = self.FixWhenItersConvergedAtNB
519                                 
520                           if self._should_fix_discrete_due_to_conv(tree_node, variable, index, lb, ub, nb):
521                                 
522                              self._fix_var(ph, tree_node, variable, index, node_min)
523                                 
524                           else: # if not fixed, then hash and slam as necessary.
525                                 
526                              if self.W_hash_history_len > 0:   
527                                 self._w_hash_acct(ph, tree_node, variable_name, index)
528                                 computed_cycle_length,msg = self.hash_hit_len(ph, tree_node, variable_name, index, False)
529                                 if (computed_cycle_length >= self.CycleLengthSlamThreshold) and ((ph._current_iteration - self._last_slam_iter) > self.PH_Iters_Between_Cycle_Slams):
530                                    # TBD: we may not want to slam immediately - it may disappear on it's own after a few iterations, depending on what other variables do.
531                                    # note: we are *not* slamming the offending variable, but a selected variable
532                                    if index is None:
533                                       print "Cycle issue detected with variable="+variable_name
534                                    else:
535                                       print "Cycle issue detected with variable="+variable_name+"["+str(index)+"]"
536                                    print msg
537                                    print "Cycle length exceeds iteration slamming threshold="+str(self.CycleLengthSlamThreshold)+"; choosing a variable to slam"
538                                    self._pick_one_and_slam_it(ph)
539                                 elif (computed_cycle_length > 1) and (computed_cycle_length < self.CycleLengthSlamThreshold):
540                                    # there was a (potential) cycle, but the slam threshold wasn't reached.
541                                    if self.ReportPotentialCycles is True:
542                                       if index is None:
543                                          print "Cycle issue detected with variable="+variable_name
544                                       else:                                         
545                                          print "Cycle issue detected with variable="+variable_name+"["+str(index)+"]"                                         
546                                       print msg
547                                       print "Taking no action to break cycle - length="+str(computed_cycle_length)+" does not exceed slam threshold="+str(self.CycleLengthSlamThreshold)
548                                 elif (computed_cycle_length >= self.CycleLengthSlamThreshold) and ((ph._current_iteration - self._last_slam_iter) > self.PH_Iters_Between_Cycle_Slams):
549                                    # we could have slammed, but we recently did and are taking a break to see if things settle down on their own.
550                                    if index is None:
551                                       print "Cycle issue detected with variable="+variable_name
552                                    else:                                         
553                                       print "Cycle issue detected with variable="+variable_name+"["+str(index)+"]"                                         
554                                    print msg
555                                    print "Taking no action to break cycle - length="+str(computed_cycle_length)+" does exceed slam threshold="+str(self.CycleLengthSlamThreshold)+ \
556                                          ", but another variable was slammed within the past "+str(self.PH_Iters_Between_Cycle_Slams)+" iterations" 
557                        else:
558
559                           # obviously don't round in the continuous case.
560                           node_min = value(tree_node._minimums[variable_name][index])
561                           node_max = value(tree_node._maximums[variable_name][index])
562
563                           # update convergence prior to checking for fixing.
564                           self._continuous_convergence_tracking(ph, tree_node, variable_name, index, node_min, node_max)
565
566                           if self._should_fix_continuous_due_to_conv(tree_node, variable, index):
567                              # fixing to max value for safety (could only be an issue due to tolerances).
568                              self._fix_var(ph, tree_node, variable, index, node_max)
569                              # note: we currently don't clam continuous variables!
570
571      # TBD: the 1 might need to be parameterized - TBD - the 1 should be the PH ITERATIONS BETWEEN CYCLE SLAMS
572      if (ph._current_iteration > self.SlamAfterIter) and \
573         ((ph._current_iteration - self._last_slam_iter) > self.PH_Iters_Between_Cycle_Slams) and \
574         (ph._converger.isImproving(self.PH_Iters_Between_Cycle_Slams)):
575         print "Slamming criteria are satisifed - accelerating convergence"
576         self._pick_one_and_slam_it(ph)
577         self._just_slammed_ = True
578      else:
579         self._just_slammed_ = False
580     
581#==================================================
582   def post_iteration_k(self, ph):
583
584      # note: we are lagging one iteration
585      # linear
586      if (self.InitialMipGap > 0 and self.FinalMipGap >= 0) and self.InitialMipGap > self.FinalMipGap:
587         m0 = self._met0
588         m = ph._converger.lastMetric()
589         mlast = ph._converger._convergence_threshold
590         g0 = self.InitialMipGap
591         glast = self.FinalMipGap
592         gap = ((m-m0)/(m0-mlast) + g0/(g0-glast))* (g0-glast)
593         if gap > g0:
594            print "***WARNING: Setting mipgap to thresholded maximal initial mapgap value="+str(g0)+"; unthresholded value="+str(gap)           
595            gap = g0
596         else:
597            print "Setting mipgap to "+str(gap)
598         ph._mipgap = gap
599
600
601#==================================================
602   def post_ph_execution(self, ph):
603
604      if self.fix_discrete_variables_at_exit is True:
605         print "WW PH extension: Fixing all discrete variables that are converged at termination"
606         self._fix_all_converged_discrete_variables(ph)
607
608#=========================
609   def Int_If_Close_Enough(self, ph, x):
610       # if x is close enough to the nearest integer, return the integer
611       # else return x
612       if abs(round(x)-x) <= ph._integer_tolerance:
613          return int(round(x))
614       else:
615          return x
616
617#=========================   
618   def _int_convergence_tracking(self, ph, tree_node, variable_name, index, node_min, node_max):
619       # keep track of cumulative iters of convergence to the same int
620       if (node_min == node_max) and (type(node_min) is types.IntType): 
621          if node_min == value(tree_node._last_converged_val[variable_name][index]):
622             tree_node._num_iters_converged[variable_name][index] = value(tree_node._num_iters_converged[variable_name][index]) + 1
623          else:
624             tree_node._num_iters_converged[variable_name][index] = 1
625             tree_node._last_converged_val[variable_name][index] = node_min
626       else:
627          tree_node._num_iters_converged[variable_name][index] = 0
628          tree_node._last_converged_val[variable_name][index] = 0.5 
629
630#=========================   
631   def _continuous_convergence_tracking(self, ph, tree_node, variable_name, index, node_min, node_max):
632       # keep track of cumulative iters of convergence to the same value within tolerance.
633       if abs(node_max - node_min) <= ph._integer_tolerance:
634          if abs(node_min - value(tree_node._last_converged_val[variable_name][index])) <= ph._integer_tolerance:
635             tree_node._num_iters_converged[variable_name][index] = value(tree_node._num_iters_converged[variable_name][index]) + 1
636          else:
637             tree_node._num_iters_converged[variable_name][index] = 1
638             tree_node._last_converged_val[variable_name][index] = node_min
639       else:
640          tree_node._num_iters_converged[variable_name][index] = 0
641          tree_node._last_converged_val[variable_name][index] = 0.2342343243223423 # TBD - avoid the magic constant!
642
643#=========================         
644   def _w_hash_acct(self, ph, tree_node, variable_name, index):
645      # do the w hash accounting work
646      # we hash on the variable ph weights, and not the values; the latter may not shift for some time, while the former should.
647      self.W_hash_rand_val = self.W_hash_seed
648      for scenario in tree_node._scenarios:
649         instance = ph._instances[scenario._name]
650         weight_parameter_name = "PHWEIGHT_"+variable_name
651         if index is None:
652            tree_node._w_hash[variable_name][ph._current_iteration] += value(getattr(instance,weight_parameter_name)[index]) * self.W_hash_rand_val
653         else:
654            tree_node._w_hash[variable_name][index,ph._current_iteration] += value(getattr(instance,weight_parameter_name)[index]) * self.W_hash_rand_val
655         self.W_hash_rand_val = (self.W_hash_b + self.W_hash_a * self.W_hash_rand_val) % self.W_hash_c
656
657#=========================
658   def dump_w_hash(self, ph, tree_node, stage):
659       # debug code
660      print "Stage=",stage._name," tree node=",tree_node._name
661      print "PH Iteration      Variable                          PH Weight Hash Value"
662      for (variable, index_template, variable_indices) in stage._variables:
663
664          variable_name = variable.name
665          variable_type = variable.domain
666
667          # TBD - should we cycle-detect on continuous vars?
668          if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
669             for index in variable_index:
670                if index is None:
671                   print "%4d        %50ls %20.5f" % (ph._current_iteration, tree_node._w_hash[variable_name][ph._current_iteration], value(tree_node._w_hash[variable_name][ph._current_iteration]))
672                else:
673                   print "%4d        %50ls %20.5f" % (ph._current_iteration, tree_node._w_hash[variable_name][index,ph._current_iteration], value(tree_node._w_hash[variable_name][index,ph._current_iteration]))
674                                                                                                                           
675#=========================
676   def hash_hit_len(self, ph, tree_node, variable_name, index, report_possible_cycles):
677      # return cycles back to closest hash hit for hashval or 0 if no hash hit
678
679      # if the values are converged, then don't report a cycle - often, the weights at convergence are 0s, and even
680      # if they aren't, they won't move if the values are uniform.
681      if (value(tree_node._num_iters_converged[variable_name][index]) == 0) and (value(tree_node._fixed_var_flag[variable_name][index]) is False):
682         current_hash_value = None
683         if index is None:
684            current_hash_value = value(tree_node._w_hash[variable_name][ph._current_iteration])
685         else:
686            current_hash_value = value(tree_node._w_hash[variable_name][index,ph._current_iteration])
687         # scan starting from the farthest point back in history to the closest - this is required to
688         # identify the longest possible cycles, which is what we want.
689         for i in range(max(ph._current_iteration - self.W_hash_history_len - 1, 1), ph._current_iteration - 1, 1):
690             this_hash_value = None
691             if index is None:
692                this_hash_value = value(tree_node._w_hash[variable_name][i])
693             else:
694                this_hash_value = value(tree_node._w_hash[variable_name][index,i])
695             if abs(this_hash_value - current_hash_value) <= ph._integer_tolerance:
696                if report_possible_cycles is True:
697                   if index is None:
698                      print "Possible cycle detected via PH weight hashing - variable="+variable_name+", node="+tree_node._name
699                   else:
700                      print "Possible cycle detected via PH weight hashing - variable="+variable_name+indexToString(index)+" node="+ tree_node._name
701                msg = "Current hash value="+str(current_hash_value)+" matched (within tolerance) hash value="+str(this_hash_value)+" found at PH iteration="+str(i)+"; cycle length="+str(ph._current_iteration - i)
702                return ph._current_iteration - i, msg
703      return 0, ""
704
705#=========================
706   def _fix_var(self, ph, tree_node, variable, index, fix_value):
707       # fix the variable, account for it and maybe output some trace information
708       # note: whether you fix at current values or not can severly impact feasibility later
709       # in the game. my original thought was below - but that didn't work out. if you
710       # make integers, well, integers, all appears to be well.
711       # IMPT: we are leaving the values for individual variables alone, at potentially
712       #       smallish and heterogeneous values. if we fix/round, we run the risk of
713       #       infeasibilities due to numerical issues. the computed value below is
714       #       strictly for output purposes. dlw note: as of aug 1 '09,
715       #       node_min and node_max should be
716       #       int if they should be (so to speak)
717
718       fixing_reported = False # to track whether you have already output the fix message for one scenario.
719
720       for scenario in tree_node._scenarios:
721
722          instance = ph._instances[scenario._name]
723                       
724          getattr(instance,variable.name)[index].fixed = True
725          getattr(instance,variable.name)[index].value = fix_value
726          tree_node._fixed_var_flag[variable.name][index] = True
727
728          variable_type = variable.domain         
729
730          if fixing_reported is False:
731             # pretty-print the index, string the trailing spaces from the strings.
732             if index is None:
733                print "Fixing variable="+variable.name+" at tree node="+tree_node._name+" to value="+str(fix_value)+"; converged for "+str(value(tree_node._num_iters_converged[variable.name][index]))+" iterations"
734             else:
735                print "Fixing variable="+variable.name+indexToString(index)+" at tree node="+tree_node._name+" to value="+str(fix_value)+"; converged for "+str(value(tree_node._num_iters_converged[variable.name][index]))+" iterations"               
736             fixing_reported = True
737             if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
738                self.cumulative_discrete_fixed_count = self.cumulative_discrete_fixed_count + 1
739             else:
740                self.cumulative_continuous_fixed_count = self.cumulative_continuous_fixed_count + 1                                                           
741
742#=========================
743   # the last 3 input arguments are the number of iterations the variable is required to
744   # be at the respective bound (or lack thereof) before fixing can be initiated.
745   def _should_fix_discrete_due_to_conv(self, tree_node, variable, index, lb_iters, ub_iters, nb_iters):
746      # return True if this should be fixed due to convergence
747      variable_name = variable.name
748
749      # jpw: i don't think this logic is correct - shouldn't "non-bound" be moved after the lb/ub checks - this doesn't check a bound!
750      # dlw reply: i meant it to mean "without regard to bound" so i have updated the document
751      if nb_iters > 0 and value(tree_node._num_iters_converged[variable_name][index]) >= nb_iters:
752            return True
753      else:
754         # there is a possibility that the variable doesn't have a bound specified, in which
755         # case we should obviously ignore the corresponding lb_iters/ub_iters/nb_iters - which
756         # should be none as well!
757         lb = None
758         ub = None
759         if variable[index].lb is not None:
760            lb = variable[index].lb()
761         if variable[index].ub is not None:
762            ub = variable[index].ub()
763         conval = value(tree_node._last_converged_val[variable_name][index])
764         # note: if they are converged node_max == node_min
765         if (lb is not None) and (lb_iters > 0) and (value(tree_node._num_iters_converged[variable_name][index]) >= lb_iters) and (conval == lb):
766            return True
767         elif (ub is not None) and (ub_iters > 0) and (value(tree_node._num_iters_converged[variable_name][index]) >= ub_iters) and (conval == ub):
768            return True
769      # if we are still here, nothing triggered fixing
770      return False
771
772#=========================
773   def _should_fix_continuous_due_to_conv(self, tree_node, variable, index):
774
775      if self.fix_continuous_variables is True:
776         if self.FixWhenItersConvergedContinuous > 0 and value(tree_node._num_iters_converged[variable.name][index]) >= self.FixWhenItersConvergedContinuous:
777               return True
778
779      # if we are still here, nothing triggered fixing
780      return False   
781
782#=========================
783   def _slam(self, ph, tree_node, variable, index):
784      # this function returns a boolean indicating if it slammed
785      # TBD in the distant future: also: slam it to somewhere it sort of wants to go
786      # e.g., the anywhere case could be to the mode
787      #   or if more than one dest is True, pick the one closest to the average
788      #   as of sept 09, it is written with the implicit assumption that only one
789      #   destination is True or that if not, then min/max trumps lb/ub and anywhere trumps all
790     
791      fix_value = False  # assume the worst
792      variable_type = variable.domain
793      variable_name = variable.name
794      if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
795         node_min = self.Int_If_Close_Enough(ph, value(tree_node._minimums[variable_name][index]))
796         node_max = self.Int_If_Close_Enough(ph, value(tree_node._maximums[variable_name][index]))
797         anywhere = round(value(tree_node._averages[variable.name][index]))
798      else:   
799         node_min = value(tree_node._minimums[variable_name][index])
800         node_max = value(tree_node._maximums[variable_name][index])
801         anywhere = value(tree_node._averages[variable.name][index])
802
803      slam_basis_string = ""
804      if self.CanSlamToLB is True:
805         fix_value = variable[index].lb()
806         slam_basis_string = "lower bound"
807      if self.CanSlamToMin is True:
808         fix_value = node_min
809         slam_basis_string = "node minimum"
810      if self.CanSlamToUB is True:
811         fix_value = variable[index].ub()
812         slam_basis_string = "upper bound"
813      if self.CanSlamToMax is True:
814         fix_value = node_max
815         slam_basis_string = "node maximum"
816      if self.CanSlamToAnywhere is True:
817         fix_value = anywhere
818         slam_basis_string = "node average (anywhere)"
819      if fix_value is False:
820         print "Warning: Not allowed to slam variable="+variable.name+str(index)+" at tree node="+tree_node._name
821         return False
822      else:
823         if index is None:
824            print "Slamming variable="+variable.name+" at tree node="+tree_node._name+" to value="+str(fix_value)+"; value="+slam_basis_string
825         else:
826            print "Slamming variable="+variable.name+indexToString(index)+" at tree node="+tree_node._name+" to value="+str(fix_value)+"; value="+slam_basis_string
827         self._fix_var(ph, tree_node, variable, index, fix_value)
828         return True
829
830#=========================
831   def _pick_one_and_slam_it(self, ph):
832
833      reference_instance = ph._model_instance
834
835      for vl in self.slam_list:
836         variable_string = str(vl)
837         full_index = None # JPW is not entirely sure of python scoping rules, so I'm declaring this outside of the if-then block.
838         variable_name = None
839         if isVariableNameIndexed(variable_string) is True:
840            pieces = variable_string.split('[')
841            variable_name = string.strip(pieces[0])
842            full_index = pieces[1].rstrip(']')
843            # the full_index is a string - tuplize it!
844            full_index = tupleizeIndexString(full_index)
845         else:
846            variable_name = variable_string
847            full_index = None
848
849         # verify that the root variable exists and grab it.
850         if variable_name not in reference_instance.active_components(Var).keys():
851            raise RuntimeError, "Unknown variable="+variable_name+" referenced while slamming. "
852         variable = reference_instance.active_components(Var)[variable_name]
853
854         didone = False;   # did we find at least one node to slam in?
855         # it is possible (even likely) that the slam list contains variable values that
856         # reside in the final stage - which, because of the initialization loops in
857         # the post_ph_initialization() method, will not have a _stage attribute defined.
858         # check for the presence of this attribute and skip if not present, as it
859         # doesn't make sense to slam variable values in the final stage anyway.
860         if hasattr(variable[full_index],'_stage') is True:
861            for tree_node in variable[full_index]._stage._tree_nodes:
862               # determine if the variable is already fixed (the trusting version...).
863               if value(tree_node._fixed_var_flag[variable_name][full_index]) is False:
864                  didone = self._slam(ph, tree_node, variable, full_index)
865            if didone:
866               self._last_slam_iter = ph._current_iteration
867               return
868         
869      print "Warning: Nothing free with a non-zero slam priority - no variable will be slammed"
870
871#==========================
872# a simple utility to fix any discrete variables to their common value, assuming they
873# are at a common value
874#==========================
875
876   def _fix_all_converged_discrete_variables(self, ph):
877
878      num_variables_fixed = 0
879
880      for stage in ph._scenario_tree._stages[:-1]: # no blending over the final stage
881         
882         for tree_node in stage._tree_nodes:
883
884            for (variable, index_template, variable_indices) in stage._variables:
885
886               variable_name = variable.name
887               variable_type = variable.domain
888
889               for index in variable_indices:
890
891                  # determine if this index is used - otherwise, don't waste time
892                  # fixing and cycle checking. for one, the code will crash :-) with
893                  # None values during the cycle checking computation!
894
895                  is_used = True # until proven otherwise                     
896                  for scenario in tree_node._scenarios:
897                     instance = ph._instances[scenario._name]
898                     if getattr(instance,variable_name)[index].status == VarStatus.unused:
899                        is_used = False
900
901                  if is_used is True:                       
902
903                     # determine if the variable is already fixed.
904                     instance_fixed_count = 0
905                     for scenario in tree_node._scenarios:
906                        instance = ph._instances[scenario._name]
907                        if getattr(instance,variable_name)[index].fixed is True:
908                           instance_fixed_count += 1
909                     if ((instance_fixed_count > 0) and (instance_fixed_count < len(tree_node._scenarios))):
910                        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
911
912                     if instance_fixed_count == 0:
913
914                        if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):                           
915                           node_min = self.Int_If_Close_Enough(ph, value(tree_node._minimums[variable_name][index]))
916                           node_max = self.Int_If_Close_Enough(ph, value(tree_node._maximums[variable_name][index]))
917
918                           if node_min == node_max:
919                              self._fix_var(ph, tree_node, variable, index, node_min)
920                              num_variables_fixed = num_variables_fixed + 1
921
922      print "Total number of variables fixed at PH termination due to convergence="+str(num_variables_fixed)
Note: See TracBrowser for help on using the repository browser.