source: coopr.pysp/stable/coopr/pysp/wwphextension.py @ 3286

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

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

........

r3201 | jwatson | 2010-10-29 13:18:17 -0600 (Fri, 29 Oct 2010) | 3 lines


Inverting order of .dat files in PySP when loading from a node representation - now root-to-leaf, instead of leaf-to-root. This allows for deeper-in-the-tree nodes to over-write parameter values defined higher in the tree, which is a more "expected" behavior than the converse. The real answer is to throw an exception if a parameter is re-defined, but we're not there yet.

........

r3217 | jwatson | 2010-11-05 11:29:42 -0600 (Fri, 05 Nov 2010) | 3 lines


Various updates to support heteogeneous index sets in PH for different nodes in the scenario tree - more work / testing remains.

........

r3218 | jwatson | 2010-11-05 12:35:51 -0600 (Fri, 05 Nov 2010) | 3 lines


More changes associated with generalizing the PySP index structures from per-stage to per-node.

........

r3220 | jwatson | 2010-11-05 20:28:29 -0600 (Fri, 05 Nov 2010) | 3 lines


Various fixes to the WW PH extension, bringing it in compliance to previous commit changes.

........

r3221 | jwatson | 2010-11-05 20:43:40 -0600 (Fri, 05 Nov 2010) | 1 line


Removing some older PySP test problems from the repository

........

r3222 | jwatson | 2010-11-05 20:55:15 -0600 (Fri, 05 Nov 2010) | 1 line


Moving PySP forestry example to local sandbox, to streamline distribution.

........

r3261 | jwatson | 2010-11-29 15:26:46 -0700 (Mon, 29 Nov 2010) | 9 lines


Adding two options to the runef and runph pysp scripts, to facilitate scenario downsampling - the case where you have a big tree, but you don't want to use it all.


The options are:
--scenario-tree-downsample-fraction=X
--scenario-tree-random-seed


The options are fairly self-explanatory - the only possible nuance is that the downsample fraction is the fraction of scenarios retained.

........

r3263 | jwatson | 2010-12-01 10:47:21 -0700 (Wed, 01 Dec 2010) | 3 lines


Corrected issue with cvar generation introduced a while back.

........

r3264 | jwatson | 2010-12-01 11:16:01 -0700 (Wed, 01 Dec 2010) | 3 lines


Adding PySP extensive form tests.

........

r3265 | wehart | 2010-12-01 13:19:56 -0700 (Wed, 01 Dec 2010) | 2 lines


Auxmenting the filter

........

r3266 | wehart | 2010-12-01 13:58:59 -0700 (Wed, 01 Dec 2010) | 2 lines


Adding further diagnostics to the filter.

........

r3267 | wehart | 2010-12-01 14:12:15 -0700 (Wed, 01 Dec 2010) | 2 lines


Another attempt to fix this filter...

........

r3271 | wehart | 2010-12-01 15:51:23 -0700 (Wed, 01 Dec 2010) | 4 lines


Simplifying filter.


The error was really a Python 2.5 portability issue in EF. :(

........

r3275 | jwatson | 2010-12-02 19:35:52 -0700 (Thu, 02 Dec 2010) | 3 lines


Fixing Python 2.5-related issue with string.translate in phutils.py - using None as the translation table is not allowed in Python 2.5.

........

  • 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.