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

Last change on this file since 5756 was 5756, checked in by wehart, 7 years ago

Misc documentation change

  • Property svn:executable set to *
File size: 66.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
20from pyutilib.misc import flatten_tuple
21
22import math
23import string
24import os
25
26#=========================
27def slam_priority_descend_compare(a, b):
28    # used to sort the variable-suffix map for slamming priority
29    value_a = a.getattrvalue("SlammingPriority")
30    value_b = b.getattrvalue("SlammingPriority")
31    return cmp(value_b, value_a)
32
33#==================================================
34class wwphextension(SingletonPlugin):
35
36    implements (phextension.IPHExtension)
37
38    def __init__(self, *args, **kwds):
39
40        # TBD - migrate all of the self attributes defined on-the-fly
41        #       in the post-post-initialization routine here!
42
43        self._configuration_filename = None
44        self._suffix_filename = None
45        self._annotation_filename = None
46
47#=========================
48    def process_suffix_file(self, ph):
49
50        # for each suffix, maintain an inverse map from the suffix name to a list of
51        # objects with that suffix. an object could be a variable (Var) or a
52        # variable value, depending on the resolution of the object for which the
53        # suffix is defined.
54        self._suffix_to_variable_map = {}
55
56        self.slam_list = []
57
58        if self._suffix_filename is None:
59            return
60
61        if os.path.exists(self._suffix_filename) is False:
62            raise RuntimeError, "***WW PH Extension: The suffix file "+self._suffix_filename+" either does not exist or cannot be read"
63
64        print "WW PH Extension: Loading variable suffixes from file="+self._suffix_filename
65
66        reference_instance = ph._model_instance
67
68        suffix_file = open(self._suffix_filename,'r')
69        for line in suffix_file.readlines():
70            line = line.strip()
71            if len(line) > 0 and line[0] != '#':
72                pieces = line.split()
73                if len(pieces) != 3:
74                    raise RuntimeError, "Illegal line=["+line+"] encountered in ww ph extension suffix file="+self._suffix_filename+"; format is variable, suffix, suffix-value."
75
76                variable_string = pieces[0]
77                suffix_name = pieces[1]
78                suffix_value = pieces[2]
79
80                # decide what type of suffix value we're dealing with.
81                is_int = False
82                is_bool = False
83                converted_value = None
84                try:
85                    converted_value = bool(suffix_value)
86                    is_bool = True
87                except valueerror:
88                    pass
89                try:
90                    converted_value = int(suffix_value)
91                    is_int = True
92                except ValueError:
93                    pass
94
95                if (is_int is False) and (is_bool is False):
96                    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
97
98                # determine if we're dealing with complete variables or indexed variables.
99                if isVariableNameIndexed(variable_string) is True:
100
101                    variable_name, index_template = extractVariableNameAndIndex(variable_string)
102
103                    # verify that the root variable exists and grab it.
104                    if variable_name not in reference_instance.active_components(Var).keys():
105                        raise RuntimeError, "Unknown variable="+variable_name+" referenced in ww ph extension suffix file="+self._suffix_filename
106                    variable = reference_instance.active_components(Var)[variable_name]
107
108                    # if the suffix isn't already an attribute of the variable, add it.
109                    if variable.attribute_defined(suffix_name) is False:
110                        variable.declare_attribute(suffix_name)
111
112                    # extract all "real", i.e., fully specified, indices matching the index template.
113                    match_indices = extractVariableIndices(variable, index_template)
114
115                    # there is a possibility that no indices match the input template.
116                    # if so, let the user know about it.
117                    if len(match_indices) == 0:
118                        raise RuntimeError, "No indices match template="+str(index_template)+" for variable="+variable_name
119
120                    # add the suffix to all variable values identified.
121                    for index in match_indices:
122
123                        variable_value = variable[index]
124
125                        # set the suffix on the variable value.
126                        variable_value.setattrvalue(suffix_name, converted_value)
127
128                        # place the variable value in the suffix->variable map, for easy searching elsewhere in this plugin.
129                        if suffix_name not in self._suffix_to_variable_map.keys():
130                            self._suffix_to_variable_map[suffix_name] = []
131                        self._suffix_to_variable_map[suffix_name].append(variable_value)
132
133                else:
134
135                    # verify that the variable exists.
136                    if variable_string not in reference_instance.active_components(Var).keys():
137                        raise RuntimeError, "Unknown variable="+variable_string+" referenced in ww ph extension suffix file="+self._suffix_filename
138
139                    variable = reference_instance.active_components(Var)[variable_string]
140
141                    # if the suffix isn't already an attribute of the variable, add it.
142                    if variable.attribute_defined(suffix_name) is False:
143                        variable.declare_attribute(suffix_name)
144
145                    # 9/14/2009 - now forcing the user to explicit specify the full
146                    # match template (e.g., "foo[*,*]") instead of just the variable
147                    # name (e.g., "foo") to represent the set of all indices.
148
149                    # if the variable is a singleton - that is, non-indexed - no brackets is fine.
150                    # we'll just tag the var[None] variable value with the (suffix,value) pair.
151                    if None not in variable._index:
152                        raise RuntimeError, "Illegal match template="+variable_string+" specified in ww ph extension suffix file="+self._suffix_filename
153
154                    # set the suffix on all variable values.
155                    for index in variable:
156
157                        variable_value = variable[index]
158
159                        # set the suffix on the variable value.
160                        variable_value.setattrvalue(suffix_name, converted_value)
161
162                        # place the variable in the suffix->variable map, for easy searching elsewhere in this plugin.
163                        if suffix_name not in self._suffix_to_variable_map.keys():
164                            self._suffix_to_variable_map[suffix_name] = []
165                        self._suffix_to_variable_map[suffix_name].append(variable_value)
166
167#      print "pre-sort suffix->variable map:"
168#      for key,value_list in self._suffix_to_variable_map.items():
169#         print "key=",key,":",
170#         for value in value_list:
171#            print value.name,"",
172#         print ""
173
174        if "SlammingPriority" in self._suffix_to_variable_map:
175            self.slam_list = self._suffix_to_variable_map["SlammingPriority"]
176
177        self.slam_list.sort(slam_priority_descend_compare)
178
179#=========================
180    def process_annotation_file(self, ph):
181        # note: these suffixes can have a string value
182        # very similar to suffixes, except for the following:
183        # annotation names are from a restricted list, and
184        # annotations can have various types that might depend on the name of the annotation
185        # not all type checking will be done here, but some might be
186
187        # note: "variable_value" is a misnomer: this thing is the variable
188
189        self._annotation_to_variable_map = {}
190        # obj effect families are special
191        self._obj_family_to_variable_map = {}
192        self._obj_family_normalized_rho = {}
193
194        if self._annotation_filename is None:
195            return
196
197        AnnotationTypes = {}
198        AnnotationTypes['going_price'] = None # Real
199        AnnotationTypes['obj_effect_family_name'] = None # String
200        AnnotationTypes['obj_effect_family_factor'] = None # Real
201        AnnotationTypes['decision_hierarchy_level'] = None # non-negative int
202        AnnotationTypes['feasibility_direction'] = ['down', 'up', 'either', 'None']
203        AnnotationTypes['relax_int'] = None # int
204        AnnotationTypes['reasonable_int'] = None # int
205        AnnotationTypes['low_int'] = None # int
206
207        if os.path.exists(self._annotation_filename) is False:
208            raise RuntimeError, "***WW PH Extension: The annotation file "+self._annotation_filename+" either does not exist or cannot be read"
209
210        print "WW PH Extension: Loading variable annotations from file="+self._annotation_filename
211
212        reference_instance = ph._model_instance
213
214        annotation_file = open(self._annotation_filename,'r')
215        for line in annotation_file.readlines():
216            line = line.strip()
217            if len(line) > 0 and line[0] != '#':
218                pieces = line.split()
219                if len(pieces) != 3:
220                    raise RuntimeError, "Illegal line=["+line+"] encountered in ww ph extension annotation file="+self._annotation_filename+"; format is variable, annotation, annotation-value."
221
222                variable_string = pieces[0]
223                annotation_name = pieces[1]
224                annotation_value = pieces[2]
225
226                # check for some input errors
227                if annotation_name not in AnnotationTypes.keys():
228                    print "Error encountered."
229                    print "Here are the annotations that can be given (they are case sensitive):"
230                    for i in AnnotationTypes:
231                        print i
232                    raise RuntimeError, "Line=["+line+"] in ww ph extension annotation file="+self._annotation_filename+"; contains unknown annotation: "+annotation_name
233
234                # check for some more input errors
235                if AnnotationTypes[annotation_name] is not None:
236                    if annotation_value not in AnnotationTypes[annotation_name]:
237                        raise RuntimeError, "Line=["+line+"] in ww ph extension annotation file="+self._annotation_filename+"; contains unknown annotation value="+annotation_value+" for: "+annotation_name
238
239                # if this is a new obj effect family, then we need new maps
240                if annotation_name == 'obj_effect_family_name':
241                    if annotation_value not in self._obj_family_to_variable_map.keys():
242                        self._obj_family_to_variable_map[annotation_value] = []
243                        self._obj_family_normalized_rho[annotation_value] = 0.0
244               
245
246                # determine if we're dealing with complete variables or indexed variables.
247                if isVariableNameIndexed(variable_string) is True:
248
249                    variable_name, index_template = extractVariableNameAndIndex(variable_string)
250
251                    # verify that the root variable exists and grab it.
252                    if variable_name not in reference_instance.active_components(Var).keys():
253                        raise RuntimeError, "Unknown variable="+variable_name+" referenced in ww ph extension annotation file="+self._annotation_filename
254                    variable = reference_instance.active_components(Var)[variable_name]
255
256                    # if the annotation isn't already an attribute of the variable, add it.
257                    if variable.attribute_defined(annotation_name) is False:
258                        variable.declare_attribute(annotation_name)
259
260                    # extract all "real", i.e., fully specified, indices matching the index template.
261                    match_indices = extractVariableIndices(variable, index_template)
262
263                    # there is a possibility that no indices match the input template.
264                    # if so, let the user know about it.
265                    if len(match_indices) == 0:
266                        raise RuntimeError, "No indices match template="+str(index_template)+" for variable="+variable_name
267
268                    # add the annotation to all variable values identified.
269                    for index in match_indices:
270
271                        variable_value = variable[index]
272
273                        # set the annotation on the variable value.
274                        variable_value.setattrvalue(annotation_name, annotation_value)
275
276                        # place the variable value in the annotation->variable map, for easy searching elsewhere in this plugin.
277                        if annotation_name not in self._annotation_to_variable_map.keys():
278                            self._annotation_to_variable_map[annotation_name] = []
279                        self._annotation_to_variable_map[annotation_name].append(variable_value)
280                        # if this is an obj effect family annotation, then we need a special map
281                        if annotation_name == 'obj_effect_family_name':
282                            self._obj_family_to_variable_map[annotation_value].append(variable_value)
283                else:
284
285                    # verify that the variable exists.
286                    if variable_string not in reference_instance.active_components(Var).keys():
287                        raise RuntimeError, "Unknown variable="+variable_string+" referenced in ww ph extension annotation file="+self._annotation_filename
288
289                    variable = reference_instance.active_components(Var)[variable_string]
290
291                    # if the annotation isn't already an attribute of the variable, add it.
292                    if variable.attribute_defined(annotation_name) is False:
293                        variable.declare_attribute(annotation_name)
294
295                    # 9/14/2009 - now forcing the user to explicit specify the full
296                    # match template (e.g., "foo[*,*]") instead of just the variable
297                    # name (e.g., "foo") to represent the set of all indices.
298
299                    # if the variable is a singleton - that is, non-indexed - no brackets is fine.
300                    # we'll just tag the var[None] variable value with the (annotation,value) pair.
301                    if None not in variable._index:
302                        raise RuntimeError, "Illegal match template="+variable_string+" specified in ww ph extension annotation file="+self._annotation_filename
303
304                    # set the annotation on all variable values.
305                    for index in variable:
306
307                        variable_value = variable[index]
308
309                        # set the annotation on the variable value.
310                        variable_value.setattrvalue(annotation_name, annotation_value)
311
312                        # place the variable in the annotation->variable map, for easy searching elsewhere in this plugin.
313                        if annotation_name not in self._annotation_to_variable_map.keys():
314                            self._annotation_to_variable_map[annotation_name] = []
315                        self._annotation_to_variable_map[annotation_name].append(variable_value)
316                        # if this is an obj effect family annotation, then we need a special map
317                        if annotation_name == 'obj_effect_family_name':
318                            self._obj_family_to_variable_map[annotation_value].append(variable_value)
319
320#==================================================
321    def post_instance_creation(self, ph):
322        """ Called after PH initialization has created the scenario instances, but before any PH-related weights/variables/parameters/etc are defined!"""
323        # we don't muck with the instances.
324        pass
325
326#==================================================
327    def post_ph_initialization(self, ph):
328
329        # set up "global" record keeping.
330        self.cumulative_discrete_fixed_count = 0
331        self.cumulative_continuous_fixed_count = 0
332
333        # we always track convergence of continuous variables, but we may not want to fix/slam them.
334        self.fix_continuous_variables = False
335
336        # there are occasions where we might want to fix any values at the end of the
337        # run if the scenarios agree - even if the normal fixing criterion (e.g.,
338        # converged for N iterations) don't apply. one example is when the term-diff
339        # is 0, at which point you really do have a solution. currently only applies
340        # to discrete variables.
341        self.fix_converged_discrete_variables_at_exit = False
342
343        # set up the mipgap parameters (zero means ignore)
344        # note: because we lag one iteration, the final will be less than requested
345        # initial and final refer to PH iteration 1 and PH iteration X, where
346        # X is the iteration at which the convergence metric hits 0.
347        self.Iteration0MipGap = 0.0
348        self.InitialMipGap = 0.10
349        self.FinalMipGap = 0.001
350
351        # "Read" the defaults for parameters that control fixing
352        # (these defaults can be overridden at the variable level)
353        # for all six of these, zero means don't do it.
354        self.Iter0FixIfConvergedAtLB = 0 # 1 or 0
355        self.Iter0FixIfConvergedAtUB = 0  # 1 or 0
356        self.Iter0FixIfConvergedAtNB = 0  # 1 or 0 (converged to a non-bound)
357        # TBD: check for range errors for all six of these
358        self.FixWhenItersConvergedAtLB = 10
359        self.FixWhenItersConvergedAtUB = 10
360        self.FixWhenItersConvergedAtNB = 12  # converged to a non-bound
361        self.FixWhenItersConvergedContinuous = 0
362
363        # "default" slamming parms:
364        # TBD: These should get ovverides from suffixes
365        # notice that for a particular var, all could be False
366        self.CanSlamToLB = False
367        self.CanSlamToMin = False
368        self.CanSlamToAnywhere = True
369        self.CanSlamToMax = False
370        self.CanSlamToUB = False
371        self.PH_Iters_Between_Cycle_Slams = 1  # zero means "slam at will"
372        self.SlamAfterIter = len(ph._scenario_tree._stages[-1]._tree_nodes)
373
374        # default params associated with fixing due to weight vector oscillation.
375        self.CheckWeightOscillationAfterIter = 0
376        self.FixIfWeightOscillationCycleLessThan = 10       
377
378        # flags enabling various rho computation schemes.
379        self.ComputeRhosWithSEP = False
380
381        self.CycleLengthSlamThreshold = len(ph._scenario_tree._stages[-1]._tree_nodes)
382        self.W_hash_history_len = max(100, self.CycleLengthSlamThreshold+1)
383
384        self.ReportPotentialCycles = 0 # do I report potential cycles, i.e., those too short to base slamming on?
385
386        # end of parms
387
388        self._last_slam_iter = -1    # dynamic
389
390        # constants for W vector hashing (low cost rand() is all we need)
391        # note: July 09, dlw is planning to eschew pre-computed random vector
392        # another note: the top level reset is OK, but really it should
393        #   done way down at the index level with a serial number (stored)
394        #   to avoid correlated hash errors
395        # the hash seed was deleted in 1.1 and we seed with the
396        self.W_hash_seed = 17  # we will reset for dynamic rand vector generation
397        self.W_hash_rand_val = self.W_hash_seed # current rand
398        self.W_hash_a = 1317       # a,b, and c are for a simple generator
399        self.W_hash_b = 27699
400        self.W_hash_c = 131072  # that period should be usually long enough for us!
401                                # (assuming fewer than c scenarios)
402
403        # set up tree storage for statistics that we care about tracking.
404        for stage in ph._scenario_tree._stages[:-1]:
405
406            for tree_node in stage._tree_nodes:
407
408                # we're adding a lot of statistics / tracking data to each tree node.
409                # these are all maps from variable name to a parameter that tracks the corresponding information.
410                tree_node._num_iters_converged = {}
411                tree_node._last_converged_val = {}
412                tree_node._w_hash = {}
413                tree_node._w_sign_vector = {} # sign vector for weights at the last PH iteration
414                tree_node._w_last_sign_flip_iter = {} # the number of iterations since the last flip in the sign (TBD - of any scenario in the vector)?               
415                tree_node._fixed_var_flag = {}
416
417                for variable_name, (reference_variable, index_template) in stage._variables.iteritems():
418
419                    if not reference_variable.attribute_defined("my_stage"):
420                        reference_variable.declare_attribute("my_stage")
421
422                    # creates the minimalistic Param objects for the tree nodes, in terms of the number of indices.
423                    variable_indices = tree_node._variable_indices[variable_name]
424                    parameter_index_set = Set(initialize=variable_indices)
425                    parameter_index_set.construct()
426
427                    new_stat_parameter_name = "NODESTAT_NUM_ITERS_CONVERGED_"+reference_variable.name
428                    new_stat_parameter = None
429                    # this bit of ugliness is due to Pyomo not correctly handling the Param construction
430                    # case when the supplied index set consists strictly of None, i.e., the source variable
431                    # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
432                    if (len(parameter_index_set) is 1) and (None in parameter_index_set):
433                        new_stat_parameter = Param(name=new_stat_parameter_name, default=0.0)
434                    else:
435                        new_stat_parameter = Param(parameter_index_set, name=new_stat_parameter_name, default=0.0)
436                    for newindex in parameter_index_set:
437                        new_stat_parameter[newindex] = 0
438                    tree_node._num_iters_converged[reference_variable.name] = new_stat_parameter
439
440                    # need to know to what value we have most recently converged
441                    new_conv_parameter_name = "NODESTAT_LAST_CONVERGED_VAL_"+reference_variable.name
442                    new_conv_parameter = None
443                    if (len(parameter_index_set) is 1) and (None in parameter_index_set):
444                        new_conv_parameter = Param(name=new_conv_parameter_name, default=0.0)
445                    else:
446                        new_conv_parameter = Param(parameter_index_set, name=new_conv_parameter_name, default=0.0)
447                    for newindex in parameter_index_set:
448                        new_conv_parameter[newindex] = 0.5 # not an int, so harmless
449                    tree_node._last_converged_val[reference_variable.name] = new_conv_parameter
450
451                    # need to know to what has been fixed
452                    new_fix_parameter_name = "NODESTAT_FIXED_FLAG_VAL_"+reference_variable.name
453                    new_fix_parameter = None
454                    if (len(parameter_index_set) is 1) and (None in parameter_index_set):
455                        new_fix_parameter = Param(name=new_fix_parameter_name, default=0.0)
456                    else:
457                        new_fix_parameter = Param(parameter_index_set, name=new_fix_parameter_name, default=0.0)
458                    for newindex in parameter_index_set:
459                        new_fix_parameter[newindex] = False
460                    tree_node._fixed_var_flag[reference_variable.name] = new_fix_parameter
461
462                    # now make the w hash value storage array
463                    new_hash_parameter_name = "W_HASH_STORAGE_"+reference_variable.name
464                    new_hash_parameter = None
465                    if (len(parameter_index_set) is 1) and (None in parameter_index_set):
466                        new_hash_parameter = Param(ph._iteration_index_set, name=new_hash_parameter_name, default=0.0)
467                        for new_index in ph._iteration_index_set:
468                            new_hash_parameter[new_index] = 0                       
469                    else:
470                        new_hash_parameter = Param(parameter_index_set * ph._iteration_index_set, name=new_hash_parameter_name, default=0.0)
471                        for new_index in new_hash_parameter:
472                            new_hash_parameter[new_index] = 0
473
474                    tree_node._w_hash[reference_variable.name] = new_hash_parameter
475
476                    # the following are easier - we aren't maintaining history.
477
478                    tree_node._w_sign_vector[reference_variable.name] = {}
479                    tree_node._w_last_sign_flip_iter[reference_variable.name] = {}
480                    for index in parameter_index_set:
481                       tree_node._w_sign_vector[reference_variable.name][index] = []
482                       tree_node._w_last_sign_flip_iter[reference_variable.name][index] = 0                   
483
484                    # JPW has no idea why the following code block is here, or if it is necessary.
485                    for index in variable_indices:
486
487                        # IMPT: it has bitten us before in this plug-in, so I'll state the obvious:
488                        #       variable values in the last stage of a stochastic program will *not*
489                        #       have a defined _stage attribute.
490                        # NOTE: the use of the reference variable here, in contrast to a instance-specific
491                        #       variable, is potentially dangerous.
492                        reference_variable[index].setattrvalue("my_stage",stage)
493
494        if self._configuration_filename is not None:
495            if os.path.exists(self._configuration_filename) is True:
496                print "WW PH Extension: Loading user-specified configuration from file=" + self._configuration_filename
497                try:
498                    execfile(self._configuration_filename)
499                except:
500                    raise RuntimeError, "Failed to load WW PH configuration file="+self._configuration_filename
501            else:
502                raise RuntimeError, "***WW PH Extension: The configuration file "+self._configuration_filename+" either does not exist or cannot be read"
503        else:
504            print "WW PH Extension: No user-specified configuration file supplied - using defaults"
505
506        # process any suffix and/or annotation data, if they exists.
507        self.process_suffix_file(ph)
508        self.process_annotation_file(ph)
509
510        # set up the mip gap for iteration 0.
511        if self.Iteration0MipGap > 0.0:
512            print "Setting mipgap to "+str(self.Iteration0MipGap)
513            ph._mipgap = self.Iteration0MipGap
514
515#==================================================
516    def post_iteration_0_solves(self, ph):
517
518        for stage in ph._scenario_tree._stages[:-1]:
519
520            for tree_node in stage._tree_nodes:
521
522                for variable_name, (variable, index_template) in stage._variables.iteritems():
523
524                    variable_type = variable.domain
525
526                    variable_indices = tree_node._variable_indices[variable_name]
527
528                    for index in variable_indices:
529
530                        # determine if this index is used - otherwise, don't waste time
531                        # fixing and cycle checking. for one, the code will crash :-) with
532                        # none values during the cycle checking computation!
533
534                        is_used = True # until proven otherwise
535                        node_average = value(tree_node._averages[variable_name][index])
536                        deviation_from_average = 0.0
537                        for scenario in tree_node._scenarios:
538                            instance = ph._instances[scenario._name]
539                            var_value = getattr(instance,variable_name)[index]
540                            if var_value.status == VarStatus.unused:
541                                is_used = False
542                            # IMPT: This is wrong - we really need the absolute tree node probability for this to work in the multi-stage case.
543                            scenario_probability = ph._scenario_tree._scenario_map[scenario._name]._probability
544                            deviation_from_average += (scenario_probability * math.fabs(var_value.value - node_average))
545                        deviation_from_average /= tree_node._conditional_probability
546
547                        if is_used is True:
548
549                            # unlikely, but variables might be fixed even at this stage, depending on
550                            # what weirdness users do prior to the iteration 0 solves.
551                            instance_fixed_count = 0
552                            for scenario in tree_node._scenarios:
553                                instance = ph._instances[scenario._name]
554                                if getattr(instance,variable_name)[index].fixed is True:
555                                    instance_fixed_count += 1
556                            if ((instance_fixed_count > 0) and (instance_fixed_count < len(tree_node._scenarios))):
557                                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
558
559                            if instance_fixed_count == 0:
560
561                                if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
562                                    node_min = self.Int_If_Close_Enough(ph, value(tree_node._minimums[variable_name][index]))
563                                    node_max = self.Int_If_Close_Enough(ph, value(tree_node._maximums[variable_name][index]))
564
565                                    # update convergence prior to checking for fixing.
566                                    self._int_convergence_tracking(ph, tree_node, variable_name, index, node_min, node_max)
567                                    attrvariable = ph._model_instance.active_components(Var)[variable_name]
568                                    attrvarvalue = attrvariable[index]
569
570                                    if (attrvariable.attribute_defined('Iter0FixIfConvergedAtLB') is True) and (attrvarvalue.getattrvalue('Iter0FixIfConvergedAtLB') is not None):
571                                        lb = attrvarvalue.getattrvalue('Iter0FixIfConvergedAtLB')
572                                    else:
573                                        lb = self.Iter0FixIfConvergedAtLB
574
575                                    if (attrvariable.attribute_defined('Iter0FixIfConvergedAtUB') is True) and (attrvarvalue.getattrvalue('Iter0FixIfConvergedAtUB') is not None):
576                                        ub = attrvarvalue.getattrvalue('Iter0FixIfConvergedAtUB')
577                                    else:
578                                        ub = self.Iter0FixIfConvergedAtUB
579
580                                    if (attrvariable.attribute_defined('Iter0FixIfConvergedAtNB') is True) and (attrvarvalue.getattrvalue('Iter0FixIfConvergedAtNB') is not None):
581                                        nb = attrvarvalue.getattrvalue('Iter0FixIfConvergedAtNB')
582                                    else:
583                                        nb = self.Iter0FixIfConvergedAtNB
584
585                                    if self._should_fix_discrete_due_to_conv(tree_node, variable, index, lb, ub, nb):
586                                        self._fix_var(ph, tree_node, variable, index, node_min)
587                                    elif self.W_hash_history_len > 0:   # if not fixed, then hash - no slamming at iteration 0
588                                        self._w_history_accounting(ph, tree_node, variable_name, index) # obviously not checking for cycles at iteration 0!
589
590                                else:
591
592                                    node_min = value(tree_node._minimums[variable_name][index])
593                                    node_max = value(tree_node._maximums[variable_name][index])
594
595                                    self._continuous_convergence_tracking(ph, tree_node, variable_name, index, node_min, node_max)
596
597                            if self.ComputeRhosWithSEP is True:
598
599                                numerator = 1.0
600
601                                node_min = self.Int_If_Close_Enough(ph, value(tree_node._minimums[variable_name][index]))
602                                node_max = self.Int_If_Close_Enough(ph, value(tree_node._maximums[variable_name][index]))
603
604                                if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
605                                    denominator = max(node_max - node_min + 1, 1)
606                                else:
607                                    denominator = max(deviation_from_average, 1)
608
609                                # CostForRho are the costs to be used as the numerator in the rho computation below.
610                                reference_variable = getattr(ph._model_instance, variable_name)
611                                cost_for_rho = reference_variable[index].getattrvalue("CostForRho")
612                                if cost_for_rho is not None:
613                                    ph.setRhoAllScenarios(variable[index], cost_for_rho * numerator / denominator)
614
615# jpw: not sure if we care about cycle detection in continuous variables?
616#                           if self.W_hash_history_len > 0:
617#                              self._w_history_accounting(ph, tree_node, variable_name, index)
618
619
620#==================================================
621    def post_iteration_0(self, ph):
622
623        self._met0 = ph._converger.lastMetric()
624
625        if (self.InitialMipGap > 0) and (self.FinalMipGap >= 0) and (self.InitialMipGap > self.FinalMipGap):
626            gap = self.InitialMipGap
627            print "Setting mipgap to "+str(gap)
628            ph._mipgap = gap
629
630#==================================================
631
632    def pre_iteration_k_solves(self, ph):
633        """ Called immediately before the iteration k solves!"""
634        # we don't muck with the PH objectives
635        pass
636
637#==================================================
638    def post_iteration_k_solves(self, ph):
639
640        for stage in ph._scenario_tree._stages[:-1]:
641
642            for tree_node in stage._tree_nodes:
643
644                for variable_name, (variable, index_template) in stage._variables.iteritems():
645
646                    variable_type = variable.domain
647
648                    variable_indices = tree_node._variable_indices[variable_name]
649
650                    for index in variable_indices:
651
652                        # determine if this index is used - otherwise, don't waste time
653                        # fixing and cycle checking. for one, the code will crash :-) with
654                        # None values during the cycle checking computation!
655
656                        is_used = True # until proven otherwise
657                        for scenario in tree_node._scenarios:
658                            instance = ph._instances[scenario._name]
659                            if getattr(instance,variable_name)[index].status == VarStatus.unused:
660                                is_used = False
661
662                        if is_used is True:
663
664                            # determine if the variable is already fixed.
665                            instance_fixed_count = 0
666                            for scenario in tree_node._scenarios:
667                                instance = ph._instances[scenario._name]
668                                if getattr(instance,variable_name)[index].fixed is True:
669                                    instance_fixed_count += 1
670                            if ((instance_fixed_count > 0) and (instance_fixed_count < len(tree_node._scenarios))):
671                                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
672
673                            if instance_fixed_count == 0:
674
675                                if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
676                                    node_min = self.Int_If_Close_Enough(ph, value(tree_node._minimums[variable_name][index]))
677                                    node_max = self.Int_If_Close_Enough(ph, value(tree_node._maximums[variable_name][index]))
678
679                                    # update convergence prior to checking for fixing.
680                                    self._int_convergence_tracking(ph, tree_node, variable_name, index, node_min, node_max)
681
682                                    # now check on permissions to converge to various placed (e.g., lb is lb permission)
683                                    attrvariable = ph._model_instance.active_components(Var)[variable_name]
684                                    attrvarvalue = attrvariable[index]
685
686                                    if (attrvariable.attribute_defined('FixWhenItersConvergedAtLB') is True) and (attrvarvalue.getattrvalue('FixWhenItersConvergedAtLB') is not None):
687                                        lb = attrvarvalue.getattrvalue('FixWhenItersConvergedAtLB')
688                                    else:
689                                        lb = self.FixWhenItersConvergedAtLB
690
691                                    if (attrvariable.attribute_defined('FixWhenItersConvergedAtUB') is True) and (attrvarvalue.getattrvalue('FixWhenItersConvergedAtUB') is not None):
692                                        ub = attrvarvalue.getattrvalue('FixWhenItersConvergedAtUB')
693                                    else:
694                                        ub = self.FixWhenItersConvergedAtUB
695
696                                    if (attrvariable.attribute_defined('FixWhenItersConvergedAtNB') is True) and (attrvarvalue.getattrvalue('FixWhenItersConvergedAtNB') is not None):
697                                        nb = attrvarvalue.getattrvalue('FixWhenItersConvergedAtNB')
698                                    else:
699                                        nb = self.FixWhenItersConvergedAtNB
700
701                                    if self._should_fix_discrete_due_to_conv(tree_node, variable, index, lb, ub, nb):
702
703                                        self._fix_var(ph, tree_node, variable, index, node_min)
704
705                                    else: # if not fixed, then hash and slam as necessary.
706
707                                        if self.W_hash_history_len > 0:
708                                            self._w_history_accounting(ph, tree_node, variable_name, index)
709                                            computed_cycle_length,msg = self.hash_hit_len(ph, tree_node, variable_name, index, False)
710                                            if (computed_cycle_length >= self.CycleLengthSlamThreshold) and ((ph._current_iteration - self._last_slam_iter) > self.PH_Iters_Between_Cycle_Slams):
711                                                # 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.
712                                                # note: we are *not* slamming the offending variable, but a selected variable
713                                                if index is None:
714                                                    print "Cycle issue detected with variable="+variable_name
715                                                else:
716                                                    print "Cycle issue detected with variable="+variable_name+"["+str(index)+"]"
717                                                print msg
718                                                print "Cycle length exceeds iteration slamming threshold="+str(self.CycleLengthSlamThreshold)+"; choosing a variable to slam"
719                                                self._pick_one_and_slam_it(ph)
720                                            elif (computed_cycle_length > 1) and (computed_cycle_length < self.CycleLengthSlamThreshold):
721                                                # there was a (potential) cycle, but the slam threshold wasn't reached.
722                                                if self.ReportPotentialCycles is True:
723                                                    if index is None:
724                                                        print "Cycle issue detected with variable="+variable_name
725                                                    else:
726                                                        print "Cycle issue detected with variable="+variable_name+"["+str(index)+"]"
727                                                    print msg
728                                                    print "Taking no action to break cycle - length="+str(computed_cycle_length)+" does not exceed slam threshold="+str(self.CycleLengthSlamThreshold)
729                                            elif (computed_cycle_length >= self.CycleLengthSlamThreshold) and ((ph._current_iteration - self._last_slam_iter) > self.PH_Iters_Between_Cycle_Slams):
730                                                # we could have slammed, but we recently did and are taking a break to see if things settle down on their own.
731                                                if index is None:
732                                                    print "Cycle issue detected with variable="+variable_name
733                                                else:
734                                                    print "Cycle issue detected with variable="+variable_name+"["+str(index)+"]"
735                                                print msg
736                                                print "Taking no action to break cycle - length="+str(computed_cycle_length)+" does exceed slam threshold="+str(self.CycleLengthSlamThreshold)+ \
737                                                      ", but another variable was slammed within the past "+str(self.PH_Iters_Between_Cycle_Slams)+" iterations"
738                                else:
739
740                                    # obviously don't round in the continuous case.
741                                    node_min = value(tree_node._minimums[variable_name][index])
742                                    node_max = value(tree_node._maximums[variable_name][index])
743
744                                    # update w statistics for whatever nefarious purposes are enabled.
745                                    if self.W_hash_history_len > 0:   
746                                        self._w_history_accounting(ph, tree_node, variable_name, index)                                                               
747
748                                    # update convergence prior to checking for fixing.
749                                    self._continuous_convergence_tracking(ph, tree_node, variable_name, index, node_min, node_max)
750
751                                    if self._should_fix_continuous_due_to_conv(tree_node, variable, index):
752                                        # fixing to max value for safety (could only be an issue due to tolerances).
753                                        self._fix_var(ph, tree_node, variable, index, node_max)
754                                        # note: we currently don't slam continuous variables!
755
756        # TBD: the 1 might need to be parameterized - TBD - the 1 should be the PH ITERATIONS BETWEEN CYCLE SLAMS
757        if (ph._current_iteration > self.SlamAfterIter) and \
758           ((ph._current_iteration - self._last_slam_iter) > self.PH_Iters_Between_Cycle_Slams) and \
759           (ph._converger.isImproving(self.PH_Iters_Between_Cycle_Slams)):
760            print "Slamming criteria are satisifed - accelerating convergence"
761            self._pick_one_and_slam_it(ph)
762            self._just_slammed_ = True
763        else:
764            self._just_slammed_ = False
765
766        ### THIS IS EXPERIMENTAL - CODE MAY BELONG SOMEWHERE ELSE ###
767        for stage in ph._scenario_tree._stages[:-1]:
768         
769           for tree_node in stage._tree_nodes:
770
771              for variable_name, (variable, index_template) in stage._variables.iteritems():
772
773                 variable_type = variable.domain
774
775                 variable_indices = tree_node._variable_indices[variable_name]
776
777                 for index in variable_indices:
778
779                    last_flip_iter = tree_node._w_last_sign_flip_iter[variable_name][index]
780                    flip_duration = ph._current_iteration - last_flip_iter
781
782                    if (self.CheckWeightOscillationAfterIter > 0) and (ph._current_iteration >= self.CheckWeightOscillationAfterIter):
783                       if (last_flip_iter is 0) or (flip_duration >= self.FixIfWeightOscillationCycleLessThan): 
784                          pass
785                       else:
786                          if self._slam(ph, tree_node, variable, index) is True:
787                             tree_node._w_last_sign_flip_iter[variable_name][index] = 0 
788                             return           
789
790#==================================================
791    def post_iteration_k(self, ph):
792
793        # note: we are lagging one iteration
794        # linear
795        if (self.InitialMipGap > 0 and self.FinalMipGap >= 0) and self.InitialMipGap > self.FinalMipGap:
796            m0 = self._met0
797            m = ph._converger.lastMetric()
798            mlast = ph._converger._convergence_threshold
799            g0 = self.InitialMipGap
800            glast = self.FinalMipGap
801            gap = ((m-m0)/(m0-mlast) + g0/(g0-glast))* (g0-glast)
802            if gap > g0:
803                print "***WARNING: Setting mipgap to thresholded maximal initial mapgap value="+str(g0)+"; unthresholded value="+str(gap)
804                gap = g0
805            else:
806                print "Setting mipgap to "+str(gap)
807            ph._mipgap = gap
808
809
810#==================================================
811    def post_ph_execution(self, ph):
812
813        if self.fix_converged_discrete_variables_at_exit is True:
814            print "WW PH extension: Fixing all discrete variables that are converged at termination"
815            self._fix_all_converged_discrete_variables(ph)
816
817#=========================
818    def Int_If_Close_Enough(self, ph, x):
819        # if x is close enough to the nearest integer, return the integer
820        # else return x
821        if abs(round(x)-x) <= ph._integer_tolerance:
822            return int(round(x))
823        else:
824            return x
825
826#=========================
827    def _int_convergence_tracking(self, ph, tree_node, variable_name, index, node_min, node_max):
828        # keep track of cumulative iters of convergence to the same int
829        if (node_min == node_max) and (type(node_min) is types.IntType):
830            if node_min == value(tree_node._last_converged_val[variable_name][index]):
831                tree_node._num_iters_converged[variable_name][index] = value(tree_node._num_iters_converged[variable_name][index]) + 1
832            else:
833                tree_node._num_iters_converged[variable_name][index] = 1
834                tree_node._last_converged_val[variable_name][index] = node_min
835        else:
836            tree_node._num_iters_converged[variable_name][index] = 0
837            tree_node._last_converged_val[variable_name][index] = 0.5
838
839#=========================
840    def _continuous_convergence_tracking(self, ph, tree_node, variable_name, index, node_min, node_max):
841        # keep track of cumulative iters of convergence to the same value within tolerance.
842        if abs(node_max - node_min) <= ph._integer_tolerance:
843            if abs(node_min - value(tree_node._last_converged_val[variable_name][index])) <= ph._integer_tolerance:
844                tree_node._num_iters_converged[variable_name][index] = value(tree_node._num_iters_converged[variable_name][index]) + 1
845            else:
846                tree_node._num_iters_converged[variable_name][index] = 1
847                tree_node._last_converged_val[variable_name][index] = node_min
848        else:
849            tree_node._num_iters_converged[variable_name][index] = 0
850            tree_node._last_converged_val[variable_name][index] = 0.2342343243223423 # TBD - avoid the magic constant!
851
852#=========================
853    def _w_history_accounting(self, ph, tree_node, variable_name, index):
854        # do the w hash accounting work
855        # we hash on the variable ph weights, and not the values; the latter may not shift for some time, while the former should.
856        self.W_hash_rand_val = self.W_hash_seed
857
858        new_sign_vector = []
859        old_sign_vector = tree_node._w_sign_vector[variable_name][index]
860
861        for scenario in tree_node._scenarios:
862           instance = ph._instances[scenario._name]
863           weight_parameter_name = "PHWEIGHT_"+variable_name
864           if index is None:
865              weight_value = value(getattr(instance,weight_parameter_name)[index])
866              weight_sign = True
867              if weight_value < 0.0:
868                 weight_sign = False
869              tree_node._w_hash[variable_name][ph._current_iteration] = value(tree_node._w_hash[variable_name][ph._current_iteration]) + weight_value * self.W_hash_rand_val
870              new_sign_vector.append(weight_sign)
871           else:
872              weight_value = value(getattr(instance,weight_parameter_name)[index])
873              weight_sign = True
874              if weight_value < 0.0:
875                 weight_sign = False
876              # the issue of "who" should flatten tuples is up in the air. for efficiency purposes (so we don't try to do it all the time), we will do it as the caller, and not expect the callee to waste (most of) their time.
877              tree_node._w_hash[variable_name][index,ph._current_iteration] = value(tree_node._w_hash[variable_name][flatten_tuple((index,ph._current_iteration))]) + weight_value * self.W_hash_rand_val
878              new_sign_vector.append(weight_sign)
879           self.W_hash_rand_val = (self.W_hash_b + self.W_hash_a * self.W_hash_rand_val) % self.W_hash_c
880
881        num_flips = 0
882        for i in xrange(0,len(old_sign_vector)):
883           if new_sign_vector[i] != old_sign_vector[i]:
884              num_flips += 1
885
886        tree_node._w_sign_vector[variable_name][index] = new_sign_vector
887
888        if num_flips >= 1:
889           tree_node._w_last_sign_flip_iter[variable_name][index] = ph._current_iteration       
890
891#=========================
892    def dump_w_hash(self, ph, tree_node, stage):
893        # debug code
894        print "Stage=",stage._name," tree node=",tree_node._name
895        print "PH Iteration      Variable                          PH Weight Hash Value"
896        for variable_name, (variable, index_template) in stage._variables.iteritems():
897
898            variable_type = variable.domain
899
900            variable_indices = tree_node._variable_indices[variable_name]
901
902            # TBD - should we cycle-detect on continuous vars?
903            if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
904                for index in variable_index:
905                    if index is None:
906                        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]))
907                    else:
908                        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]))
909
910#=========================
911    def hash_hit_len(self, ph, tree_node, variable_name, index, report_possible_cycles):
912        # return cycles back to closest hash hit for hashval or 0 if no hash hit
913
914        # if the values are converged, then don't report a cycle - often, the weights at convergence are 0s, and even
915        # if they aren't, they won't move if the values are uniform.
916        if (value(tree_node._num_iters_converged[variable_name][index]) == 0) and (value(tree_node._fixed_var_flag[variable_name][index]) is False):
917            current_hash_value = None
918            if index is None:
919                current_hash_value = value(tree_node._w_hash[variable_name][ph._current_iteration])
920            else:
921                current_hash_value = value(tree_node._w_hash[variable_name][flatten_tuple((index,ph._current_iteration))])
922            # scan starting from the farthest point back in history to the closest - this is required to
923            # identify the longest possible cycles, which is what we want.
924            for i in range(max(ph._current_iteration - self.W_hash_history_len - 1, 1), ph._current_iteration - 1, 1):
925                this_hash_value = None
926                if index is None:
927                    this_hash_value = value(tree_node._w_hash[variable_name][i])
928                else:
929                    this_hash_value = value(tree_node._w_hash[variable_name][flatten_tuple((index,i))])
930                if abs(this_hash_value - current_hash_value) <= ph._integer_tolerance:
931                    if report_possible_cycles is True:
932                        if index is None:
933                            print "Possible cycle detected via PH weight hashing - variable="+variable_name+", node="+tree_node._name
934                        else:
935                            print "Possible cycle detected via PH weight hashing - variable="+variable_name+indexToString(index)+" node="+ tree_node._name
936                    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)
937                    return ph._current_iteration - i, msg
938        return 0, ""
939
940#=========================
941    def _fix_var(self, ph, tree_node, variable, index, fix_value):
942        # fix the variable, account for it and maybe output some trace information
943        # note: whether you fix at current values or not can severly impact feasibility later
944        # in the game. my original thought was below - but that didn't work out. if you
945        # make integers, well, integers, all appears to be well.
946        # IMPT: we are leaving the values for individual variables alone, at potentially
947        #       smallish and heterogeneous values. if we fix/round, we run the risk of
948        #       infeasibilities due to numerical issues. the computed value below is
949        #       strictly for output purposes. dlw note: as of aug 1 '09,
950        #       node_min and node_max should be
951        #       int if they should be (so to speak)
952
953        fixing_reported = False # to track whether you have already output the fix message for one scenario.
954
955        for scenario in tree_node._scenarios:
956
957            instance = ph._instances[scenario._name]
958
959            getattr(instance,variable.name)[index].fixed = True
960            getattr(instance,variable.name)[index].value = fix_value
961            tree_node._fixed_var_flag[variable.name][index] = True
962
963            variable_type = variable.domain
964
965            if fixing_reported is False:
966                # pretty-print the index, string the trailing spaces from the strings.
967                if index is None:
968                    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"
969                else:
970                    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"
971                fixing_reported = True
972                if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
973                    self.cumulative_discrete_fixed_count = self.cumulative_discrete_fixed_count + 1
974                else:
975                    self.cumulative_continuous_fixed_count = self.cumulative_continuous_fixed_count + 1
976
977            # record the specific variable fixed for the PH client.
978            ph._instance_variables_fixed_detail[scenario._name].append((variable.name, index))
979
980#=========================
981    # the last 3 input arguments are the number of iterations the variable is required to
982    # be at the respective bound (or lack thereof) before fixing can be initiated.
983    def _should_fix_discrete_due_to_conv(self, tree_node, variable, index, lb_iters, ub_iters, nb_iters):
984        # return True if this should be fixed due to convergence
985        variable_name = variable.name
986
987        # 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!
988        # dlw reply: i meant it to mean "without regard to bound" so i have updated the document
989        if nb_iters > 0 and value(tree_node._num_iters_converged[variable_name][index]) >= nb_iters:
990            return True
991        else:
992            # there is a possibility that the variable doesn't have a bound specified, in which
993            # case we should obviously ignore the corresponding lb_iters/ub_iters/nb_iters - which
994            # should be none as well!
995            lb = None
996            ub = None
997            if variable[index].lb is not None:
998                lb = variable[index].lb()
999            if variable[index].ub is not None:
1000                ub = variable[index].ub()
1001            conval = value(tree_node._last_converged_val[variable_name][index])
1002            # note: if they are converged node_max == node_min
1003            if (lb is not None) and (lb_iters > 0) and (value(tree_node._num_iters_converged[variable_name][index]) >= lb_iters) and (conval == lb):
1004                return True
1005            elif (ub is not None) and (ub_iters > 0) and (value(tree_node._num_iters_converged[variable_name][index]) >= ub_iters) and (conval == ub):
1006                return True
1007        # if we are still here, nothing triggered fixing
1008        return False
1009
1010#=========================
1011    def _should_fix_continuous_due_to_conv(self, tree_node, variable, index):
1012
1013        if self.fix_continuous_variables is True:
1014            if self.FixWhenItersConvergedContinuous > 0 and value(tree_node._num_iters_converged[variable.name][index]) >= self.FixWhenItersConvergedContinuous:
1015                return True
1016
1017        # if we are still here, nothing triggered fixing
1018        return False
1019
1020#=========================
1021    def _slam(self, ph, tree_node, variable, index):
1022        # this function returns a boolean indicating if it slammed
1023        # TBD in the distant future: also: slam it to somewhere it sort of wants to go
1024        # e.g., the anywhere case could be to the mode
1025        #   or if more than one dest is True, pick the one closest to the average
1026        #   as of sept 09, it is written with the implicit assumption that only one
1027        #   destination is True or that if not, then min/max trumps lb/ub and anywhere trumps all
1028
1029        fix_value = False  # assume the worst
1030        variable_type = variable.domain
1031        variable_name = variable.name
1032        if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
1033            node_min = self.Int_If_Close_Enough(ph, value(tree_node._minimums[variable_name][index]))
1034            node_max = self.Int_If_Close_Enough(ph, value(tree_node._maximums[variable_name][index]))
1035            anywhere = round(value(tree_node._averages[variable.name][index]))
1036        else:
1037            node_min = value(tree_node._minimums[variable_name][index])
1038            node_max = value(tree_node._maximums[variable_name][index])
1039            anywhere = value(tree_node._averages[variable.name][index])
1040
1041        slam_basis_string = ""
1042        if self.CanSlamToLB is True:
1043            fix_value = variable[index].lb()
1044            slam_basis_string = "lower bound"
1045        if self.CanSlamToMin is True:
1046            fix_value = node_min
1047            slam_basis_string = "node minimum"
1048        if self.CanSlamToUB is True:
1049            fix_value = variable[index].ub()
1050            slam_basis_string = "upper bound"
1051        if self.CanSlamToMax is True:
1052            fix_value = node_max
1053            slam_basis_string = "node maximum"
1054        if self.CanSlamToAnywhere is True:
1055            fix_value = anywhere
1056            slam_basis_string = "node average (anywhere)"
1057        if fix_value is False:
1058            print "Warning: Not allowed to slam variable="+variable.name+str(index)+" at tree node="+tree_node._name
1059            return False
1060        else:
1061            if index is None:
1062                print "Slamming variable="+variable.name+" at tree node="+tree_node._name+" to value="+str(fix_value)+"; value="+slam_basis_string
1063            else:
1064                print "Slamming variable="+variable.name+indexToString(index)+" at tree node="+tree_node._name+" to value="+str(fix_value)+"; value="+slam_basis_string
1065            self._fix_var(ph, tree_node, variable, index, fix_value)
1066            return True
1067
1068#=========================
1069    def _pick_one_and_slam_it(self, ph):
1070
1071        reference_instance = ph._model_instance
1072
1073        for vl in self.slam_list:
1074            variable_string = str(vl)
1075            full_index = None # JPW is not entirely sure of python scoping rules, so I'm declaring this outside of the if-then block.
1076            variable_name = None
1077            if isVariableNameIndexed(variable_string) is True:
1078                pieces = variable_string.split('[')
1079                variable_name = string.strip(pieces[0])
1080                full_index = pieces[1].rstrip(']')
1081                # the full_index is a string - tuplize it!
1082                full_index = tupleizeIndexString(full_index)
1083            else:
1084                variable_name = variable_string
1085                full_index = None
1086
1087            # verify that the root variable exists and grab it.
1088            if variable_name not in reference_instance.active_components(Var).keys():
1089                raise RuntimeError, "Unknown variable="+variable_name+" referenced while slamming. "
1090            variable = reference_instance.active_components(Var)[variable_name]
1091
1092            didone = False;   # did we find at least one node to slam in?
1093            # it is possible (even likely) that the slam list contains variable values that
1094            # reside in the final stage - which, because of the initialization loops in
1095            # the post_ph_initialization() method, will not have a _stage attribute defined.
1096            # check for the presence of this attribute and skip if not present, as it
1097            # doesn't make sense to slam variable values in the final stage anyway.
1098            if variable.attribute_defined("my_stage") is True:
1099                variable_stage = variable[full_index].getattrvalue("my_stage")
1100
1101                if variable_stage is not None: # None => the variable isn't used.
1102
1103                   for tree_node in variable_stage._tree_nodes:
1104                   # determine if the variable is already fixed (the trusting version...).
1105                       if value(tree_node._fixed_var_flag[variable_name][full_index]) is False:
1106                           didone = self._slam(ph, tree_node, variable, full_index)
1107                   if didone:
1108                       self._last_slam_iter = ph._current_iteration
1109                       return
1110
1111        print "Warning: Nothing free with a non-zero slam priority - no variable will be slammed"
1112        # DLW says: Hey, look at this: if we were to start deleting from the slam list this would be wrong."
1113        if len(self.slam_list) == 0:
1114            print "   (No Slamming Priorities were specified in a suffix file.)"
1115
1116#==========================
1117# a simple utility to fix any discrete variables to their common value, assuming they
1118# are at a common value
1119#==========================
1120
1121    def _fix_all_converged_discrete_variables(self, ph):
1122
1123        num_variables_fixed = 0
1124
1125        for stage in ph._scenario_tree._stages[:-1]: # no blending over the final stage
1126
1127            for tree_node in stage._tree_nodes:
1128
1129                for variable_name, (variable, index_template) in stage._variables.iteritems():
1130
1131                    variable_type = variable.domain
1132
1133                    variable_indices = tree_node._variable_indices[variable_name]
1134
1135                    for index in variable_indices:
1136
1137                        # determine if this index is used - otherwise, don't waste time
1138                        # fixing and cycle checking. for one, the code will crash :-) with
1139                        # None values during the cycle checking computation!
1140
1141                        is_used = True # until proven otherwise
1142                        for scenario in tree_node._scenarios:
1143                            instance = ph._instances[scenario._name]
1144                            if getattr(instance,variable_name)[index].status == VarStatus.unused:
1145                                is_used = False
1146
1147                        if is_used is True:
1148
1149                            # determine if the variable is already fixed.
1150                            instance_fixed_count = 0
1151                            for scenario in tree_node._scenarios:
1152                                instance = ph._instances[scenario._name]
1153                                if getattr(instance,variable_name)[index].fixed is True:
1154                                    instance_fixed_count += 1
1155                            if ((instance_fixed_count > 0) and (instance_fixed_count < len(tree_node._scenarios))):
1156                                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
1157
1158                            if instance_fixed_count == 0:
1159
1160                                if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
1161                                    node_min = self.Int_If_Close_Enough(ph, value(tree_node._minimums[variable_name][index]))
1162                                    node_max = self.Int_If_Close_Enough(ph, value(tree_node._maximums[variable_name][index]))
1163
1164                                    if node_min == node_max:
1165                                        self._fix_var(ph, tree_node, variable, index, node_min)
1166                                        num_variables_fixed = num_variables_fixed + 1
1167
1168        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.