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

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

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

........

r3048 | jwatson | 2010-09-27 09:33:53 -0500 (Mon, 27 Sep 2010) | 3 lines


Modifying PySP CSV solution writer to eliminate leading and trailing parantheses from index names, change embedded commas with colons, and eliminate embedded spaces.

........

r3055 | khunter | 2010-09-30 14:11:13 -0500 (Thu, 30 Sep 2010) | 8 lines


Reorganization of the runph command line options into logical
groups. I might've messed up any intentional ordering within the
groups as I also alphabetized the options, so please to correct as
necessary.


Added -m and -i because Joe commented more than once on the
verbosity of the run(ef|ph) command lines.

........

r3056 | dlwoodr | 2010-09-30 16:41:52 -0500 (Thu, 30 Sep 2010) | 1 line

........

r3057 | dlwoodr | 2010-09-30 16:43:26 -0500 (Thu, 30 Sep 2010) | 1 line


The last check-in simply fixed the missing references problem. This document still needs to be brought up to date with the software (some command line options are not yet documented)

........

r3059 | dlwoodr | 2010-09-30 17:25:25 -0500 (Thu, 30 Sep 2010) | 1 line


Added most of the new command line options and made a few minor edits. Changed the version given as the title to 1.11 (from 1.1).

........

r3060 | dlwoodr | 2010-09-30 17:31:50 -0500 (Thu, 30 Sep 2010) | 1 line


Update the WW Comp Mgt Sci Reference

........

r3065 | jwatson | 2010-10-01 14:46:57 -0500 (Fri, 01 Oct 2010) | 3 lines


Reworked some of the options groups in the runph script, and created some new aliases.

........

r3072 | wehart | 2010-10-04 11:14:32 -0500 (Mon, 04 Oct 2010) | 2 lines


Fixes due to mutability change in coopr.pyomo.

........

r3073 | jwatson | 2010-10-04 20:35:24 -0500 (Mon, 04 Oct 2010) | 3 lines


Complete set of changes to make PySP compatible with latest immutable parameters change.

........

r3074 | jwatson | 2010-10-04 20:45:14 -0500 (Mon, 04 Oct 2010) | 3 lines


Added routines to the PySP scenario tree object to return the stage for a constraint and variable.

........

r3075 | jwatson | 2010-10-04 21:16:58 -0500 (Mon, 04 Oct 2010) | 3 lines


Forgot some debug output in the last commit.

........

r3078 | jwatson | 2010-10-05 09:46:33 -0500 (Tue, 05 Oct 2010) | 3 lines


Fixing long-present bug in network flow model, exposed by Bill's recent commit.

........

r3084 | jwatson | 2010-10-12 22:12:17 -0500 (Tue, 12 Oct 2010) | 3 lines


Fixed bug in constraint stage determination code in PySP - forgot to ignore constant terms, which of course don't have variables!

........

r3088 | jwatson | 2010-10-13 12:52:13 -0500 (Wed, 13 Oct 2010) | 3 lines


Added --traceback option to runph, in order to dump the full stack trace when an exception is thrown.

........

r3089 | jwatson | 2010-10-13 12:57:45 -0500 (Wed, 13 Oct 2010) | 3 lines


Added --traceback option to runef, mirroring recent modification to pysp.

........

r3092 | jwatson | 2010-10-13 14:13:23 -0500 (Wed, 13 Oct 2010) | 3 lines


Update to PySP to suppress canonical expression representations when using ASL. Otherwise, things choke badly.

........

r3093 | jwatson | 2010-10-13 14:38:37 -0500 (Wed, 13 Oct 2010) | 3 lines


More PySP fixes associated with ipopt integration. Nothing wrong with ipopt per se - it was just doing things we didn't expect, but were perfectly reasonable.

........

r3096 | jwatson | 2010-10-13 21:32:35 -0500 (Wed, 13 Oct 2010) | 3 lines


Enhancements to PySP to deal with reporting issues caused by the lack of objective function values in SOL files (in general, and those produced by ipopt in particular).

........

r3100 | khunter | 2010-10-14 15:54:49 -0500 (Thu, 14 Oct 2010) | 2 lines


NFC: remove EOL whitespace

........

r3101 | khunter | 2010-10-14 16:06:40 -0500 (Thu, 14 Oct 2010) | 3 lines


Reorganization of the runef command line options into logical
groups. Also alphabetized the options within the groups.

........

r3102 | khunter | 2010-10-14 16:07:47 -0500 (Thu, 14 Oct 2010) | 2 lines


Add -m and -i because there über common options.

........

r3104 | khunter | 2010-10-15 00:08:51 -0500 (Fri, 15 Oct 2010) | 3 lines


Fix an erroneous exit in the no traceback specified case, if not
exception has occurred.

........

r3105 | khunter | 2010-10-15 01:41:47 -0500 (Fri, 15 Oct 2010) | 2 lines


NFC: EOL whitespace removal

........

r3106 | khunter | 2010-10-15 01:43:00 -0500 (Fri, 15 Oct 2010) | 3 lines


Drive-by mild clean-up of an init function while I researching
another issue.

........

r3107 | jwatson | 2010-10-15 09:27:03 -0500 (Fri, 15 Oct 2010) | 3 lines


Minor touch-ups relating to traceback facilities.

........

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