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

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

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

........

r1956 | jwatson | 2009-12-02 17:56:53 -0700 (Wed, 02 Dec 2009) | 3 lines


Added --scenario-solver-options and --ef-solver-options options to the "runph" script.

........

r1957 | dlwoodr | 2009-12-03 14:17:35 -0700 (Thu, 03 Dec 2009) | 2 lines


Documentation updates for pysp

........

r1974 | wehart | 2009-12-06 17:20:56 -0700 (Sun, 06 Dec 2009) | 2 lines


Updating PyPI categories

........

r1978 | jwatson | 2009-12-10 21:29:33 -0700 (Thu, 10 Dec 2009) | 3 lines


Eliminated exception-handling logic when loading user-defined extension modules in PH. The range of exceptions is too large, and for debugging purposes, it is more useful to see the raw trace output.

........

r1979 | jwatson | 2009-12-10 22:23:17 -0700 (Thu, 10 Dec 2009) | 5 lines


Biggest enhancement: The efwriter command-line script now has the option to output a CVaR-weighted objective term. Not extensively tested, but behaves sane on a number of small test cases.


Improved exception handling and error diagnostics in both the runph and efwriter scripts.

........

r1985 | jwatson | 2009-12-12 10:45:17 -0700 (Sat, 12 Dec 2009) | 3 lines


Modified PH to only use warm-starts if a solver has the capability!

........

r1998 | jwatson | 2009-12-13 15:17:58 -0700 (Sun, 13 Dec 2009) | 3 lines


Changed references to _component to active_component.

........

r2026 | wehart | 2009-12-21 23:27:06 -0700 (Mon, 21 Dec 2009) | 2 lines


Attempting to update PH. I'm not sure if this works, since I don't know how to test PH.

........

r2029 | jwatson | 2009-12-22 09:52:21 -0700 (Tue, 22 Dec 2009) | 3 lines


Some fixes to the ef writer based on Bill's recent changes to _initialize_constraint.

........

r2035 | jwatson | 2009-12-22 21:10:32 -0700 (Tue, 22 Dec 2009) | 3 lines


Added --scenario-mipgap option to PH script. Added _mipgap attribute to PH object. This attribute is mirrored to the solver plugin at the initiation of each iteration, after any PH extensions have the opportunity to provide a new value to the attribute. It is currently made use of by the WW PH extension.

........

r2037 | dlwoodr | 2009-12-23 14:38:43 -0700 (Wed, 23 Dec 2009) | 2 lines


add this example from Fernando

........

r2038 | dlwoodr | 2009-12-23 14:46:56 -0700 (Wed, 23 Dec 2009) | 3 lines


finish the job: we can now duplicate Fernando's example

........

r2039 | jwatson | 2009-12-23 15:13:54 -0700 (Wed, 23 Dec 2009) | 3 lines


Missed fix with new constraint initialization syntax in PH linearization.

........

r2058 | jwatson | 2009-12-29 10:57:58 -0700 (Tue, 29 Dec 2009) | 3 lines


Missed some _initialize_constraint function calls within the PySP EF writer during the recent switch to the corresponding "add" method.

........

r2059 | jwatson | 2009-12-29 10:58:34 -0700 (Tue, 29 Dec 2009) | 3 lines


Enabling garbage collection by default in PH.

........

r2060 | jwatson | 2009-12-29 10:59:05 -0700 (Tue, 29 Dec 2009) | 3 lines


Elimnating some debug output.

........

r2061 | jwatson | 2009-12-29 11:07:47 -0700 (Tue, 29 Dec 2009) | 3 lines


Fixing some option documentation in PH.

........

r2062 | jwatson | 2009-12-29 12:00:37 -0700 (Tue, 29 Dec 2009) | 3 lines


Added ef-mipgap option to PH scripts.

........

  • Property svn:executable set to *
File size: 44.3 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.active_components(Var).keys():
103                  raise RuntimeError, "Unknown variable="+variable_name+" referenced in ww ph extension suffix file="+self._suffix_filename
104               variable = reference_instance.active_components(Var)[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.active_components(Var).keys():
135                  raise RuntimeError, "Unknown variable="+variable_string+" referenced in ww ph extension suffix file="+self._suffix_filename
136
137               variable = reference_instance.active_components(Var)[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._mipgap = 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.active_components(Var)[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._mipgap = 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.active_components(Var)[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._mipgap = 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.active_components(Var).keys():
751            raise RuntimeError, "Unknown variable="+variable_name+" referenced while slamming. "
752         variable = reference_instance.active_components(Var)[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.