source: coopr.pysp/stable/2.1/coopr/pysp/wwphextension.py @ 1952

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

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

........

r1884 | wehart | 2009-11-14 00:11:40 -0700 (Sat, 14 Nov 2009) | 2 lines


Update to current revision number.

........

r1890 | jwatson | 2009-11-14 17:08:09 -0700 (Sat, 14 Nov 2009) | 3 lines


Added option to PH to fix all discrete variables that are converged upon termination. Functionality is within the WW PH extension. Not for typical use, but useful for debugging.

........

r1895 | jwatson | 2009-11-16 10:44:48 -0700 (Mon, 16 Nov 2009) | 3 lines


Changed PySP package manager information from Bill to JPW.

........

r1911 | jwatson | 2009-11-19 10:24:07 -0700 (Thu, 19 Nov 2009) | 3 lines


Restructured PH code to allow for general breakpoint distribution strategies. Added option to runph (--breakpoint-strategy) to specify particular strategies.

........

r1912 | jwatson | 2009-11-19 13:24:04 -0700 (Thu, 19 Nov 2009) | 3 lines


Added PH breakpoint distribution strategy that uniformly concentrates pieces between the observed node-min and node-max values, as opposed to the full lb/ub range.

........

r1913 | jwatson | 2009-11-19 13:46:13 -0700 (Thu, 19 Nov 2009) | 3 lines


Added Woodruff PH breakpoint distribution strategy.

........

r1914 | jwatson | 2009-11-19 15:22:56 -0700 (Thu, 19 Nov 2009) | 3 lines


Added exponentially biased breakpoint distribution technique to PH.

........

r1922 | jwatson | 2009-11-20 15:50:27 -0700 (Fri, 20 Nov 2009) | 3 lines


Fix to make PH compliant with recent solution re-structuring.

........

r1925 | jwatson | 2009-11-20 22:42:49 -0700 (Fri, 20 Nov 2009) | 3 lines


Fixes to make PySP consistent with the variable lower/upper bound changes described in the previous commit. Fixed a few bugs in the farmer example now that bounds specified through parameters are working correctly.

........

r1928 | jwatson | 2009-11-21 11:33:10 -0700 (Sat, 21 Nov 2009) | 3 lines


Fixed an issue with construction of piecewise linear breakpoints in PH - if the average was at a bound, divide-by-zero was triggering.

........

r1929 | jwatson | 2009-11-21 11:45:34 -0700 (Sat, 21 Nov 2009) | 3 lines


Fixed another numerical issue in breakpoint computation in PH.

........

r1930 | jwatson | 2009-11-21 12:53:31 -0700 (Sat, 21 Nov 2009) | 3 lines


Adding 288 scenario forestry instance, for a serious challenge!

........

r1931 | jwatson | 2009-11-21 13:01:18 -0700 (Sat, 21 Nov 2009) | 3 lines


Added a feasible/new chile 48 instance (the old, hand-constructed instance seemed to be infeasible), and the EF for the 288 scenario instance.

........

r1932 | jwatson | 2009-11-21 13:37:45 -0700 (Sat, 21 Nov 2009) | 3 lines


Bug fix to PH when LB/UB aren't specified and you are linearizing.

........

r1933 | jwatson | 2009-11-21 14:24:35 -0700 (Sat, 21 Nov 2009) | 3 lines


Updating documentation on forestry PySP example.

........

r1934 | jwatson | 2009-11-22 15:38:31 -0700 (Sun, 22 Nov 2009) | 3 lines


Fix to profiler option in PH.

........

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