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

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

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

........

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


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

........

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


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

........

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


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

........

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


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

........

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


Removing some older PySP test problems from the repository

........

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


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

........

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


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


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


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

........

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


Corrected issue with cvar generation introduced a while back.

........

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


Adding PySP extensive form tests.

........

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


Auxmenting the filter

........

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


Adding further diagnostics to the filter.

........

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


Another attempt to fix this filter...

........

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


Simplifying filter.


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

........

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


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

........

File size: 42.7 KB
Line 
1#  _________________________________________________________________________
2#
3#  Coopr: A COmmon Optimization Python Repository
4#  Copyright (c) 2008 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 sys
12import types
13import copy
14import os.path
15import traceback
16
17from math import fabs
18
19import random
20
21from coopr.pyomo import *
22from phutils import *
23
24class ScenarioTreeNode(object):
25
26   #
27   # update the _solutions attribute of a tree node, given a specific
28   # variable/match-template/variable-index triple.
29   #
30   def _update_solution_map(self, variable, match_template, scenario_instance_map):
31
32      # don't bother copying bounds for variables, as the values stored
33      # here are computed elsewhere - and that entity is responsible for
34      # ensuring feasibility. this also leaves room for specifying infeasible
35      # or partial solutions.
36      new_variable_name = variable.name     
37      new_variable_index = getattr(scenario_instance_map[self._scenarios[0]._name], variable.name)._index
38      new_variable = None
39      if (len(new_variable_index) is 1) and (None in new_variable_index):
40         new_variable = Var(name=new_variable_name)
41      else:
42         new_variable = Var(new_variable_index, name=new_variable_name)
43      new_variable.construct()
44
45      variable_indices = self._variable_indices[variable.name]     
46
47      # by default, deactive all variable values - we're copying the
48      # full index set of a variable, which is necessarily wasteful and
49      # incorrect. then, do a second pass and activate those indicies
50      # that are actually associated with this tree node.
51      for index in new_variable_index:
52         new_variable[index].deactivate()
53
54      for index in variable_indices:
55         new_variable[index].activate()
56
57      self._solutions[new_variable_name] = new_variable
58
59   #
60   # initialize the _solutions attribute of a tree node, from scratch.
61   #
62
63   def _initialize_solution_map(self, scenario_instance_map):
64
65      # clear whatever was there before.
66      self._solutions = {}
67
68      # NOTE: Given the expense, this should really be optional - don't
69      #       construct unless we're actually going to use!
70      # for each variable referenced in the stage, clone the variable
71      # for purposes of storing solutions. we are being wasteful in
72      # terms copying indices that may not be referenced in the stage.
73      # this is something that we might revisit if space/performance
74      # is an issue (space is the most likely issue)
75      for variable, match_template in self._stage._variables:
76         self._update_solution_map(variable, match_template, scenario_instance_map)
77
78      self._solution_map_initialized = True
79
80
81   """ Constructor
82   """
83   def __init__(self, name, conditional_probability, stage, initialize_solution=False, scenario_instance_map={}):
84
85      self._name = name
86      self._stage = stage
87      self._parent = None
88      self._children = [] # a collection of ScenarioTreeNodes
89      self._conditional_probability = conditional_probability # conditional on parent
90      self._scenarios = [] # a collection of all Scenarios passing through this node in the tree
91      self._variable_indices = {} # a map from a variable name to the indices blended at this node.
92
93      # general use statistics for the variables at each node.
94      # each attribute is a map between the variable name and a
95      # parameter (over the same index set) encoding the corresponding
96      # statistic computed over all scenarios for that node. the
97      # parameters are named as the source variable name suffixed
98      # by one of: "NODEMIN", "NODEAVG", and "NODEMAX".
99      # NOTE: the averages are probability_weighted - the min/max
100      #       values are not.
101      # NOTE: the parameter names are basically irrelevant, and the
102      #       convention is assumed to be enforced by whoever populates
103      #       these parameters.
104      self._averages = {}
105      self._minimums = {}
106      self._maximums = {}
107
108      # a flag indicating whether the _solutions attribute has been properly initialized.
109      self._solution_map_initialized = False
110
111      # solution (variable) values for this node. assumed to be distinct
112      # from self._averages, as the latter are not necessarily feasible.
113      # objects in the map are actual pyomo Var instances; keys are
114      # variable names.
115      self._solutions = {}
116
117      if initialize_solution is True:
118         self._initialize_solution_map(scenario_instance_map)
119
120   #
121   # given a set of scenario instances, compute the set of indices being blended
122   # for each variable at this node.
123   #
124
125   def defineVariableIndexSets(self, scenario_instances):
126
127      # find a representative scenario instance belonging to this node. the
128      # first scenario is as good as any.
129      scenario_instance = scenario_instances[self._scenarios[0]._name]
130
131      for reference_variable, match_template in self._stage._variables:
132
133         # the stage variable simply references the variable object in the
134         # reference scenario instance - we need to grab the variable in the
135         # scenario instance, as the index set might be different.
136         variable_name = reference_variable.name
137
138         instance_variable = getattr(scenario_instance, variable_name)
139
140         match_indices = extractVariableIndices(instance_variable, match_template)
141
142         self._variable_indices[variable_name] = match_indices
143
144   #
145   # copies the parameter values values from the _averages attribute
146   # into the _solutions attribute - only for active variable values.
147   #
148
149   def snapshotSolutionFromAverages(self, scenario_instance_map):
150
151      if self._solution_map_initialized is False:
152         self._initialize_solution_map(scenario_instance_map)
153
154      for variable_name, variable in self._solutions.items():
155
156         # try and grab the corresponding averages parameter - if it
157         # doesn't exist, throw an exception.
158         average_parameter = None
159         try:
160            average_parameter = self._averages[variable_name]
161         except:
162            raise RuntimeError, "No averages parameter present on tree node="+self._name+" for variable="+variable_name
163
164         for index in variable._index:
165            if variable[index].active is True:
166               variable[index] = value(average_parameter[index])
167
168   #
169   # computes the solution values from the composite scenario instances at this tree node.
170   # the input scenario_instance_map is a map from scenario name to instance objects.
171   #
172
173   def snapshotSolutionFromInstances(self, scenario_instance_map):
174
175      if self._solution_map_initialized is False:
176         self._initialize_solution_map(scenario_instance_map)
177
178      for variable_name, variable in self._solutions.items():
179         for index in variable._index:
180            if variable[index].active is True:
181               node_probability = 0.0
182               avg = 0.0
183               num_scenarios_with_index = 0
184               for scenario in self._scenarios:
185                  scenario_instance = scenario_instance_map[scenario._name]
186                  node_probability += scenario._probability
187                  scenario_variable = getattr(scenario_instance, variable.name)
188                  if index in scenario_variable:
189                     num_scenarios_with_index = num_scenarios_with_index + 1
190                     var_value = value(getattr(scenario_instance, variable.name)[index])
191                     avg += (scenario._probability * var_value)
192               if num_scenarios_with_index > 0:
193                  variable[index] = avg / node_probability
194
195   #
196   # a utility to compute the cost of the current node plus the expected costs of child nodes.
197   #
198
199   def computeExpectedNodeCost(self, scenario_instance_map):
200
201      # IMPT: This implicitly assumes convergence across the scenarios - if not, garbage results.
202      instance = scenario_instance_map[self._scenarios[0]._name]
203      stage_cost_variable = instance.active_components(Var)[self._stage._cost_variable[0].name][self._stage._cost_variable[1]]
204      if stage_cost_variable.value is not None:
205         my_cost = stage_cost_variable.value
206      else:
207         # depending on the situation (and the model), the stage cost variables might not have values.
208         my_cost = 0.0
209      child_cost = 0.0
210      for child in self._children:
211         child_cost += (child._conditional_probability * child.computeExpectedNodeCost(scenario_instance_map))
212      return my_cost + child_cost
213
214
215class Stage(object):
216
217   """ Constructor
218   """
219   def __init__(self, *args, **kwds):
220
221      self._name = ""
222
223      # a collection of ScenarioTreeNode objects.
224      self._tree_nodes = []
225
226      # a collection of pairs consisting of (1) a reference to a Pyomo model Var object (in the reference scenario instance) and
227      # (2) the original match template string (for output purposes). the specific indices that match belong to the tree node.
228      self._variables = []
229
230      # a tuple consisting of (1) a reference to a pyomo model Var that computes the stage-specific cost and (2) the corresponding
231      # index. the index *is* the sole index in the cost variable, as the cost variable refers to a single variable index.
232      self._cost_variable = (None, None)
233
234   #
235   # add a new variable to the stage, which will include updating the solution maps for each associated ScenarioTreeNode.
236   #
237   def add_variable(self, variable, match_template, scenario_instance_map):
238
239      self._variables.append((variable, match_template))
240      match_indices = extractVariableIndices(variable, match_template)
241
242      for tree_node in self._tree_nodes:
243         tree_node._variable_indices[variable.name] = match_indices
244         tree_node._update_solution_map(variable, match_template, scenario_instance_map)
245
246class Scenario(object):
247
248   """ Constructor
249   """
250   def __init__(self, *args, **kwds):
251
252      self._name = None
253      self._leaf_node = None  # allows for construction of node list
254      self._node_list = []    # sequence from parent to leaf of ScenarioTreeNodes
255      self._probability = 0.0 # the unconditional probability for this scenario, computed from the node list
256
257class ScenarioTree(object):
258
259   # a utility to construct the stage objects for this scenario tree.
260   # operates strictly by side effects, initializing the self
261   # _stages and _stage_map attributes.
262   def _construct_stages(self, stage_ids, stage_variable_ids, stage_cost_variable_ids):
263
264      # construct the stage objects, which will leave them
265      # largely uninitialized - no variable information, in particular.
266      for stage_name in stage_ids:
267         new_stage = Stage()
268         new_stage._name = stage_name
269         self._stages.append(new_stage)
270         self._stage_map[stage_name] = new_stage
271
272      # initialize the variable collections (blended and cost) for each stage.
273      # these are references, not copies.
274      for stage_id in stage_variable_ids.keys():
275
276         if stage_id not in self._stage_map.keys():
277            raise ValueError, "Unknown stage=" + stage_id + " specified in scenario tree constructor (stage->variable map)"
278
279         stage = self._stage_map[stage_id]
280         variable_ids = stage_variable_ids[stage_id]
281
282         for variable_string in variable_ids:
283
284            if isVariableNameIndexed(variable_string) is True:
285
286               variable_name, index_template = extractVariableNameAndIndex(variable_string)
287
288               # validate that the variable exists and extract the reference.
289               if variable_name not in self._reference_instance.active_components(Var):
290                  raise ValueError, "Variable=" + variable_name + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
291               variable = self._reference_instance.active_components(Var)[variable_name]
292
293               stage._variables.append((variable, index_template))
294
295            else:
296
297               # verify that the variable exists.
298               if variable_string not in self._reference_instance.active_components(Var).keys():
299                  raise RuntimeError, "Unknown variable=" + variable_string + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
300
301               variable = self._reference_instance.active_components(Var)[variable_string]
302
303               # 9/14/2009 - now forcing the user to explicit specify the full
304               # match template (e.g., "foo[*,*]") instead of just the variable
305               # name (e.g., "foo") to represent the set of all indices.
306
307               # if the variable is a singleton - that is, non-indexed - no brackets is fine.
308               # we'll just tag the var[None] variable value with the (suffix,value) pair.
309               if None not in variable._index:
310                  raise RuntimeError, "Variable="+variable_string+" is an indexed variable, and templates must specify an index match; encountered in scenario tree specification for model="+self._reference_instance.name
311
312               match_indices = []
313               match_indices.append(None)
314
315               stage._variables.append((variable, ""))
316
317      for stage_id in stage_cost_variable_ids.keys():
318
319         if stage_id not in self._stage_map.keys():
320            raise ValueError, "Unknown stage=" + stage_id + " specified in scenario tree constructor (stage->cost variable map)"
321         stage = self._stage_map[stage_id]
322
323         cost_variable_string = value(stage_cost_variable_ids[stage_id]) # de-reference is required to access the parameter value
324
325         # to be extracted from the string.
326         cost_variable_name = None
327         cost_variable = None
328         cost_variable_index = None
329
330         # do the extraction.
331         if isVariableNameIndexed(cost_variable_string) is True:
332
333            cost_variable_name, index_template = extractVariableNameAndIndex(cost_variable_string)
334
335            # validate that the variable exists and extract the reference.
336            if cost_variable_name not in self._reference_instance.active_components(Var):
337               raise ValueError, "Variable=" + cost_variable_name + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
338            cost_variable = self._reference_instance.active_components(Var)[cost_variable_name]
339
340            # extract all "real", i.e., fully specified, indices matching the index template.
341            match_indices = extractVariableIndices(cost_variable, index_template)
342
343            # only one index can be supplied for a stage cost variable.
344            if len(match_indices) > 1:
345                msg = 'Only one index can be specified for a stage cost '     \
346                      'variable - %s match template "%s" for variable "%s" ;' \
347                      ' encountered in scenario tree specification for model' \
348                      ' "%s"'
349                raise RuntimeError, msg % (
350                  len(match_indices),
351                  index_template,
352                  cost_variable_name,
353                  self._reference_instance.name
354                )
355            elif len(match_indices) == 0:
356                msg = 'Stage cost index not found: %s[%s]\n'                  \
357                      'Do you have an off-by-one miscalculation, or did you ' \
358                      'forget to specify it in ReferenceModel.dat?'
359                raise RuntimeError, msg % ( cost_variable_name, index_template )
360
361            cost_variable_index = match_indices[0]
362
363         else:
364
365            cost_variable_name = cost_variable_string
366
367            # validate that the variable exists and extract the reference
368            if cost_variable_name not in self._reference_instance.active_components(Var):
369               raise ValueError, "Cost variable=" + cost_variable_name + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
370            cost_variable = self._reference_instance.active_components(Var)[cost_variable_name]
371
372         # store the validated info.
373         stage._cost_variable = (cost_variable, cost_variable_index)
374
375   """ Constructor
376       Arguments:
377           scenarioinstance     - the reference (deterministic) scenario instance.
378           scenariotreeinstance - the pyomo model specifying all scenario tree (text) data.
379           scenariobundlelist   - a list of scenario names to retain, i.e., cull the rest to create a reduced tree!
380   """
381   def __init__(self, *args, **kwds):
382
383      self._name = None # some arbitrary identifier
384      self._reference_instance = kwds.pop( 'scenarioinstance', None ) # the reference (deterministic) base model
385
386      # the core objects defining the scenario tree.
387      self._tree_nodes = [] # collection of ScenarioTreeNodes
388      self._stages = [] # collection of Stages - assumed to be in time-order. the set (provided by the user) itself *must* be ordered.
389      self._scenarios = [] # collection of Scenarios
390
391      # dictionaries for the above.
392      self._tree_node_map = {}
393      self._stage_map = {}
394      self._scenario_map = {}
395
396      # a boolean indicating how data for scenario instances is specified.
397      # possibly belongs elsewhere, e.g., in the PH algorithm.
398      self._scenario_based_data = None
399
400      scenario_tree_instance = kwds.pop( 'scenariotreeinstance', None )
401      scenario_bundle_list = kwds.pop( 'scenariobundlelist', None )
402
403      # process the keyword options
404      for key in kwds.keys():
405         print >>sys.stderr, "Unknown option '%s' specified in call to ScenarioTree constructor" % key
406
407      if self._reference_instance is None:
408         raise ValueError, "A reference scenario instance must be supplied in the ScenarioTree constructor"
409
410      if scenario_tree_instance is None:
411         raise ValueError, "A scenario tree instance must be supplied in the ScenarioTree constructor"
412
413      node_ids = scenario_tree_instance.Nodes
414      node_child_ids = scenario_tree_instance.Children
415      node_stage_ids = scenario_tree_instance.NodeStage
416      node_probability_map = scenario_tree_instance.ConditionalProbability
417      stage_ids = scenario_tree_instance.Stages
418      stage_variable_ids = scenario_tree_instance.StageVariables
419      stage_cost_variable_ids = scenario_tree_instance.StageCostVariable
420      scenario_ids = scenario_tree_instance.Scenarios
421      scenario_leaf_ids = scenario_tree_instance.ScenarioLeafNode
422      scenario_based_data = scenario_tree_instance.ScenarioBasedData
423
424      # save the method for instance data storage.
425      self._scenario_based_data = scenario_based_data()
426
427      # the input stages must be ordered, for both output purposes and knowledge of the final stage.
428      if stage_ids.ordered is False:
429         raise ValueError, "An ordered set of stage IDs must be supplied in the ScenarioTree constructor"
430
431      #
432      # construct the actual tree objects
433      #
434
435      # construct the stage objects w/o any linkages first; link them up
436      # with tree nodes after these have been fully constructed.
437      self._construct_stages(stage_ids, stage_variable_ids, stage_cost_variable_ids)
438
439      # construct the tree node objects themselves in a first pass,
440      # and then link them up in a second pass to form the tree.
441      # can't do a single pass because the objects may not exist.
442      for tree_node_name in node_ids:
443
444         if tree_node_name not in node_stage_ids:
445            raise ValueError, "No stage is assigned to tree node=" + tree_node._name
446
447         stage_name = value(node_stage_ids[tree_node_name])
448         if stage_name not in self._stage_map.keys():
449            raise ValueError, "Unknown stage=" + stage_name + " assigned to tree node=" + tree_node._name
450
451         new_tree_node = ScenarioTreeNode(tree_node_name,
452                                          value(node_probability_map[tree_node_name]),
453                                          self._stage_map[stage_name],
454                                          self._reference_instance)
455
456         self._tree_nodes.append(new_tree_node)
457         self._tree_node_map[tree_node_name] = new_tree_node
458         self._stage_map[stage_name]._tree_nodes.append(new_tree_node)
459
460      # link up the tree nodes objects based on the child id sets.
461      for this_node in self._tree_nodes:
462         this_node._children = []
463         if this_node._name in node_child_ids: # otherwise, you're at a leaf and all is well.
464            child_ids = node_child_ids[this_node._name]
465            for child_id in child_ids:
466               if child_id in self._tree_node_map.keys():
467                  child_node = self._tree_node_map[child_id]
468                  this_node._children.append(self._tree_node_map[child_id])
469                  if child_node._parent is None:
470                     child_node._parent = this_node
471                  else:
472                     raise ValueError, "Multiple parents specified for tree node="+child_id+"; existing parent node="+child_node._parent._name+"; conflicting parent node="+this_node._name
473               else:
474                  raise ValueError, "Unknown child tree node=" + child_id + " specified for tree node=" + this_node._name
475
476      # at this point, the scenario tree nodes and the stages are set - no
477      # two-pass logic necessary when constructing scenarios.
478      for scenario_name in scenario_ids:
479
480         # IMPT: the name of the scenario is assumed to have no '_' (underscore) characters in the identifier.
481         #       this is required when writing the extensive form, e.g., the scenario is used in the extensive
482         #       form as a prefix on variable and constraint names. this convention simplifies parsing on the
483         #       back end; if the underscore isn't used as a reserved character, then some other separator
484         #       symbol would be required, or we would have to engage in some complex prefix matching with
485         #       all possible scenario names.
486         if string.find(scenario_name, "_") != -1:
487            raise ValueError, "By convention, scenario names in PySP cannot contain underscore (_) characters; the scenario in violation="+scenario_name
488
489         new_scenario = Scenario()
490         new_scenario._name=scenario_name
491
492         if scenario_name not in scenario_leaf_ids.keys():
493            raise ValueError, "No leaf tree node specified for scenario=" + scenario_name
494         else:
495            scenario_leaf_node_name = value(scenario_leaf_ids[scenario_name])
496            if scenario_leaf_node_name not in self._tree_node_map.keys():
497               raise ValueError, "Uknown tree node=" + scenario_leaf_node_name + " specified as leaf of scenario=" + scenario_name
498            else:
499               new_scenario._leaf_node = self._tree_node_map[scenario_leaf_node_name]
500
501         current_node = new_scenario._leaf_node
502         probability = 1.0
503         while current_node is not None:
504            new_scenario._node_list.append(current_node)
505            current_node._scenarios.append(new_scenario) # links the scenarios to the nodes to enforce necessary non-anticipativity
506            probability *= current_node._conditional_probability
507            current_node = current_node._parent
508         new_scenario._node_list.reverse()
509         new_scenario._probability = probability
510
511         self._scenarios.append(new_scenario)
512         self._scenario_map[scenario_name] = new_scenario
513
514      # for output purposes, it is useful to known the maximal length of identifiers
515      # in the scenario tree for any particular category. I'm building these up
516      # incrementally, as they are needed. 0 indicates unassigned.
517      self._max_scenario_id_length = 0
518
519      # does the actual traversal to populate the members.
520      self.computeIdentifierMaxLengths()
521
522      # if a sub-bundle of scenarios has been specified, mark the
523      # active scenario tree components and compress the tree.
524      if scenario_bundle_list is not None:
525         self.compress(scenario_bundle_list)
526
527   #
528   # given a set of scenario instances, compute the set of variable indices being blended at each node.
529   # this can't be done until the scenario instances are available, as different scenarios can have
530   # different index sets.
531   #
532
533   def defineVariableIndexSets(self, scenario_instances):
534
535      for tree_node in self._tree_nodes:
536
537         tree_node.defineVariableIndexSets(scenario_instances)
538
539   #
540   # is the indicated scenario in the tree?
541   #
542
543   def contains_scenario(self, name):
544      return name in self._scenario_map.keys()
545
546   #
547   # get the scenario object from the tree.
548   #
549
550   def get_scenario(self, name):
551      return self._scenario_map[name]
552
553   #
554   # compute the scenario cost for the input instance, i.e.,
555   # the sum of all stage cost variables.
556   #
557
558   def compute_scenario_cost(self, instance):
559      aggregate_cost = 0.0
560      for stage in self._stages:
561            instance_cost_variable = instance.active_components(Var)[stage._cost_variable[0].name][stage._cost_variable[1]]
562            if instance_cost_variable.value is not None:
563                # depending on the situation (and the model), the stage cost variables might not have values.
564               aggregate_cost += instance_cost_variable
565      return aggregate_cost
566
567   #
568   # utility for compressing or culling a scenario tree based on
569   # a provided list of scenarios (specified by name) to retain -
570   # all non-referenced components are eliminated. this particular
571   # method compresses *in-place*, i.e., via direct modification
572   # of the scenario tree structure.
573   #
574
575   def compress(self, scenario_bundle_list):
576
577      # scan for and mark all referenced scenarios and
578      # tree nodes in the bundle list - all stages will
579      # obviously remain.
580      for scenario_name in scenario_bundle_list:
581         if scenario_name not in self._scenario_map:
582            raise ValueError, "Scenario="+scenario_name+" selected for bundling not present in scenario tree"
583         scenario = self._scenario_map[scenario_name]
584         scenario.retain = True
585
586         # chase all nodes comprising this scenario,
587         # marking them for retention.
588         for node in scenario._node_list:
589            node.retain = True
590
591      # scan for any non-retained scenarios and tree nodes.
592      scenarios_to_delete = []
593      tree_nodes_to_delete = []
594      for scenario in self._scenarios:
595         if hasattr(scenario, "retain") is True:
596            delattr(scenario, "retain")
597            pass
598         else:
599            scenarios_to_delete.append(scenario)
600            del self._scenario_map[scenario._name]
601
602      for tree_node in self._tree_nodes:
603         if hasattr(tree_node, "retain") is True:
604            delattr(tree_node, "retain")
605            pass
606         else:
607            tree_nodes_to_delete.append(tree_node)
608            del self._tree_node_map[tree_node._name]
609
610      # JPW does not claim the following routines are
611      # the most efficient. rather, they get the job
612      # done while avoiding serious issues with
613      # attempting to remove elements from a list that
614      # you are iterating over.
615
616      # delete all references to unmarked scenarios
617      # and child tree nodes in the scenario tree node
618      # structures.
619      for tree_node in self._tree_nodes:
620         for scenario in scenarios_to_delete:
621            if scenario in tree_node._scenarios:
622               tree_node._scenarios.remove(scenario)
623         for node_to_delete in tree_nodes_to_delete:
624            if node_to_delete in tree_node._children:
625               tree_node._children.remove(node_to_delete)
626
627      # delete all references to unmarked tree nodes
628      # in the scenario tree stage structures.
629      for stage in self._stages:
630         for tree_node in tree_nodes_to_delete:
631            if tree_node in stage._tree_nodes:
632               stage._tree_nodes.remove(tree_node)
633
634      # delete all unreferenced entries from the core scenario
635      # tree data structures.
636      for scenario in scenarios_to_delete:
637         self._scenarios.remove(scenario)
638      for tree_node in tree_nodes_to_delete:
639         self._tree_nodes.remove(tree_node)
640
641
642      # re-normalize the conditional probabilities of the
643      # children at each tree node.
644      for tree_node in self._tree_nodes:
645         sum_child_probabilities = 0.0
646         for child_node in tree_node._children:
647            sum_child_probabilities += child_node._conditional_probability
648         for child_node in tree_node._children:
649            child_node._conditional_probability = child_node._conditional_probability / sum_child_probabilities
650
651      # re-compute the absolute scenario probabilities based
652      # on the re-normalized conditional node probabilities.
653      for scenario in self._scenarios:
654         probability = 1.0
655         for tree_node in scenario._node_list:
656            probability = probability * tree_node._conditional_probability
657         scenario._probability = probability
658
659   #
660   # utility for automatically selecting a proportion of scenarios from the
661   # tree to retain, eliminating the rest.
662   #
663
664   def downsample(self, fraction_to_retain, random_seed, verbose=False):
665
666      random.seed(random_seed)
667
668      random_sequence=range(0,len(self._scenarios))
669      random.shuffle(random_sequence)
670
671      number_to_retain = max(int(round(float(len(random_sequence)*fraction_to_retain))), 1)
672
673      scenario_bundle_list = []
674      for i in range(0,number_to_retain):
675         scenario_bundle_list.append(self._scenarios[random_sequence[i]]._name)
676
677      if verbose is True:
678         print "Downsampling scenario tree - retained scenarios: "+str(scenario_bundle_list)
679
680      self.compress(scenario_bundle_list)
681
682
683   #
684   # returns the root node of the scenario tree
685   #
686
687   def findRootNode(self):
688
689      for tree_node in self._tree_nodes:
690         if tree_node._parent is None:
691            return tree_node
692      return None
693
694   #
695   # a utility function to compute, based on the current scenario tree content,
696   # the maximal length of identifiers in various categories.
697   #
698
699   def computeIdentifierMaxLengths(self):
700
701      self._max_scenario_id_length = 0
702      for scenario in self._scenarios:
703         if len(scenario._name) > self._max_scenario_id_length:
704            self._max_scenario_id_length = len(scenario._name)
705
706   #
707   # a utility function to (partially, at the moment) validate a scenario tree
708   #
709
710   def validate(self):
711
712      # for any node, the sum of conditional probabilities of the children should sum to 1.
713      for tree_node in self._tree_nodes:
714         sum_probabilities = 0.0
715         if len(tree_node._children) > 0:
716            for child in tree_node._children:
717               sum_probabilities += child._conditional_probability
718            if abs(1.0 - sum_probabilities) > 0.000001:
719               print "The child conditional probabilities for tree node=" + tree_node._name + " sum to " + `sum_probabilities`
720               return False
721
722      # ensure that there is only one root node in the tree
723      num_roots = 0
724      root_ids = []
725      for tree_node in self._tree_nodes:
726         if tree_node._parent is None:
727            num_roots += 1
728            root_ids.append(tree_node._name)
729
730      if num_roots != 1:
731         print "Illegal set of root nodes detected: " + `root_ids`
732         return False
733
734      # there must be at least one scenario passing through each tree node.
735      for tree_node in self._tree_nodes:
736         if len(tree_node._scenarios) == 0:
737            print "There are no scenarios associated with tree node=" + tree_node._name
738            return False
739
740      return True
741
742   #
743   # copies the parameter values stored in any tree node _averages attribute
744   # into any tree node _solutions attribute - only for active variable values.
745   #
746
747   def snapshotSolutionFromAverages(self, scenario_instance_map):
748
749      for tree_node in self._tree_nodes:
750
751         tree_node.snapshotSolutionFromAverages(scenario_instance_map)
752
753   #
754   # computes the variable values at each tree node from the input scenario instances.
755   #
756
757   def snapshotSolutionFromInstances(self, scenario_instance_map):
758
759      for tree_node in self._tree_nodes:
760
761         tree_node.snapshotSolutionFromInstances(scenario_instance_map)
762
763   #
764   # a utility to determine the stage to which the input variable belongs.
765   # this is horribly inefficient, lacking an inverse map. fortunately,
766   # it isn't really called that often (yet). stage membership is determined
767   # by comparing the input variable name with the reference instance
768   # variable name (which is what the scenario tree refers to) and the
769   # associated indices.
770   #
771
772   def variableStage(self, variable, index):
773
774      for stage in self._stages:
775         for (stage_var, match_template) in stage._variables:
776            if (variable.name == stage_var.name) and (index in match_indices):
777               return stage
778
779         # IMPT: this is a temporary hack - the real fix is to force users to
780         # have every variable assigned to some stage in the model, either
781         # automatically or explicitly.
782         if (variable.name == stage._cost_variable[0].name):
783            return stage
784
785      raise RuntimeError, "The variable="+str(variable.name)+", index="+str(index)+" does not belong to any stage in the scenario tree"
786
787   #
788   # a utility to determine the stage to which the input constraint "belongs".
789   # a constraint belongs to the latest stage in which referenced variables
790   # in the constraint appears in that stage.
791   # input is a constraint is of type "Constraint", and an index of that
792   # constraint - which might be None in the case of singleton constraints.
793   # currently doesn't deal with SOS constraints, for no real good reason.
794   # returns an instance of a Stage object.
795   # IMPT: this method works on the canonical representation ("repn" attribute)
796   #       of a constraint. this implies that pre-processing of the instance
797   #       has been performed.
798   # NOTE: there is still the issue of whether the contained variables really
799   #       belong to the same model, but that is a different issue we won't
800   #       address right now (e.g., what does it mean for a constraint in an
801   #       extensive form binding instance to belong to a stage?).
802   #
803
804   def constraintStage(self, constraint, index):
805
806      largest_stage_index = -1
807      largest_stage = None
808
809      canonical_repn = constraint[index].repn
810      for degree, terms in canonical_repn.items():
811         if (degree != -1) and (degree != 0): # ignore constant terms and the variable definitions themselves.
812            for var_key, coefficient in terms.items():
813               var_value = canonical_repn[-1][var_key.keys()[0]]
814               var_stage = self.variableStage(var_value.var, var_value.index)
815               var_stage_index = self._stages.index(var_stage)
816
817               if var_stage_index > largest_stage_index:
818                  largest_stage_index = var_stage_index
819                  largest_stage = var_stage
820
821      return largest_stage
822
823   #
824   # a utility function to pretty-print the static/non-cost information associated with a scenario tree
825   #
826
827   def pprint(self):
828
829      print "Scenario Tree Detail"
830
831      print "----------------------------------------------------"
832      if self._reference_instance is not None:
833         print "Model=" + self._reference_instance.name
834      else:
835         print "Model=" + "Unassigned"
836      print "----------------------------------------------------"
837      print "Tree Nodes:"
838      print ""
839      for tree_node_name in sorted(self._tree_node_map.keys()):
840         tree_node = self._tree_node_map[tree_node_name]
841         print "\tName=" + tree_node_name
842         if tree_node._stage is not None:
843            print "\tStage=" + str(tree_node._stage._name)
844         else:
845            print "\t Stage=None"
846         if tree_node._parent is not None:
847            print "\tParent=" + tree_node._parent._name
848         else:
849            print "\tParent=" + "None"
850         if tree_node._conditional_probability is not None:
851            print "\tConditional probability=%4.4f" % tree_node._conditional_probability
852         else:
853            print "\tConditional probability=" + "***Undefined***"
854         print "\tChildren:"
855         if len(tree_node._children) > 0:
856            for child_node in sorted(tree_node._children, cmp=lambda x,y: cmp(x._name, y._name)):
857               print "\t\t" + child_node._name
858         else:
859            print "\t\tNone"
860         print "\tScenarios:"
861         if len(tree_node._scenarios) == 0:
862            print "\t\tNone"
863         else:
864            for scenario in sorted(tree_node._scenarios, cmp=lambda x,y: cmp(x._name, y._name)):
865               print "\t\t" + scenario._name
866         print ""
867      print "----------------------------------------------------"
868      print "Stages:"
869      for stage_name in sorted(self._stage_map.keys()):
870         stage = self._stage_map[stage_name]
871         print "\tName=" + str(stage_name)
872         print "\tTree Nodes: "
873         for tree_node in sorted(stage._tree_nodes, cmp=lambda x,y: cmp(x._name, y._name)):
874            print "\t\t" + tree_node._name
875         print "\tVariables: "
876         for (variable, index_template) in stage._variables:
877            print "\t\t",variable.name,":",index_template
878         print "\tCost Variable: "
879         if stage._cost_variable[1] is None:
880            print "\t\t" + stage._cost_variable[0].name
881         else:
882            print "\t\t" + stage._cost_variable[0].name + indexToString(stage._cost_variable[1])
883         print ""
884      print "----------------------------------------------------"
885      print "Scenarios:"
886      for scenario_name in sorted(self._scenario_map.keys()):
887         scenario = self._scenario_map[scenario_name]
888         print "\tName=" + scenario_name
889         print "\tProbability=%4.4f" % scenario._probability
890         if scenario._leaf_node is None:
891            print "\tLeaf node=None"
892         else:
893            print "\tLeaf node=" + scenario._leaf_node._name
894         print "\tTree node sequence:"
895         for tree_node in scenario._node_list:
896            print "\t\t" + tree_node._name
897         print ""
898      print "----------------------------------------------------"
899
900   #
901   # a utility function to pretty-print the solution associated with a scenario tree
902   #
903
904   def pprintSolution(self, epsilon=1.0e-5):
905
906      print "----------------------------------------------------"
907      print "Tree Nodes:"
908      print ""
909      for tree_node_name in sorted(self._tree_node_map.keys()):
910         tree_node = self._tree_node_map[tree_node_name]
911         print "\tName=" + tree_node_name
912         if tree_node._stage is not None:
913            print "\tStage=" + tree_node._stage._name
914         else:
915            print "\t Stage=None"
916         if tree_node._parent is not None:
917            print "\tParent=" + tree_node._parent._name
918         else:
919            print "\tParent=" + "None"
920         print "\tVariables: "
921         for (variable, index_template) in tree_node._stage._variables:
922            indices = tree_node._variable_indices[variable.name]
923            solution_variable = tree_node._solutions[variable.name]
924            if (len(indices) == 1) and (indices[0] == None):
925               # if this is a singleton variable, then it should necessarily be active -
926               # otherwise, it wouldn't be referenced in the stage!!!
927               value = solution_variable[None].value
928               if fabs(value) > epsilon:
929                  print "\t\t"+variable.name+"="+str(value)
930            else:
931               for index in indices:
932                  if (solution_variable[index].active is True) and (index in solution_variable):
933                     value = solution_variable[index].value
934                     if (value is not None) and (fabs(value) > epsilon):
935                        print "\t\t"+variable.name+indexToString(index)+"="+str(value)
936         print ""
937
938   #
939   # a utility function to pretty-print the cost information associated with a scenario tree
940   #
941
942   def pprintCosts(self, scenario_instance_map):
943
944      print "Scenario Tree Costs"
945      print "***WARNING***: Assumes full (or nearly so) convergence of scenario solutions at each node in the scenario tree - computed costs are invalid otherwise"
946
947      print "----------------------------------------------------"
948      if self._reference_instance is not None:
949         print "Model=" + self._reference_instance.name
950      else:
951         print "Model=" + "Unassigned"
952
953      print "----------------------------------------------------"
954      print "Tree Nodes:"
955      print ""
956      for tree_node_name in sorted(self._tree_node_map.keys()):
957         tree_node = self._tree_node_map[tree_node_name]
958         print "\tName=" + tree_node_name
959         if tree_node._stage is not None:
960            print "\tStage=" + tree_node._stage._name
961         else:
962            print "\t Stage=None"
963         if tree_node._parent is not None:
964            print "\tParent=" + tree_node._parent._name
965         else:
966            print "\tParent=" + "None"
967         if tree_node._conditional_probability is not None:
968            print "\tConditional probability=%4.4f" % tree_node._conditional_probability
969         else:
970            print "\tConditional probability=" + "***Undefined***"
971         print "\tChildren:"
972         if len(tree_node._children) > 0:
973            for child_node in sorted(tree_node._children, cmp=lambda x,y: cmp(x._name, y._name)):
974               print "\t\t" + child_node._name
975         else:
976            print "\t\tNone"
977         print "\tScenarios:"
978         if len(tree_node._scenarios) == 0:
979            print "\t\tNone"
980         else:
981            for scenario in sorted(tree_node._scenarios, cmp=lambda x,y: cmp(x._name, y._name)):
982               print "\t\t" + scenario._name
983         print "\tExpected node cost=%10.4f" % tree_node.computeExpectedNodeCost(scenario_instance_map)
984         print ""
985
986      print "----------------------------------------------------"
987      print "Scenarios:"
988      print ""
989      for scenario_name, scenario in sorted(self._scenario_map.iteritems()):
990         instance = scenario_instance_map[scenario_name]
991
992         print "\tName=" + scenario_name
993         print "\tProbability=%4.4f" % scenario._probability
994
995         if scenario._leaf_node is None:
996            print "\tLeaf Node=None"
997         else:
998            print "\tLeaf Node=" + scenario._leaf_node._name
999
1000         print "\tTree node sequence:"
1001         for tree_node in scenario._node_list:
1002            print "\t\t" + tree_node._name
1003
1004         aggregate_cost = 0.0
1005         for stage in self._stages:
1006            instance_cost_variable = instance.active_components(Var)[stage._cost_variable[0].name][stage._cost_variable[1]]
1007            if instance_cost_variable.value is not None:
1008               print "\tStage=%20s     Cost=%10.4f" % (stage._name, instance_cost_variable.value)
1009               cost = instance_cost_variable.value
1010            else:
1011               print "\tStage=%20s     Cost=%10s" % (stage._name, "Not Rprted.")
1012               cost = 0.0
1013            aggregate_cost += cost
1014         print "\tTotal scenario cost=%10.4f" % aggregate_cost
1015         print ""
1016      print "----------------------------------------------------"
1017
Note: See TracBrowser for help on using the repository browser.