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

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

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

........

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


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

........

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


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


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

........

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

........

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


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

........

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


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

........

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


Update the WW Comp Mgt Sci Reference

........

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


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

........

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


Fixes due to mutability change in coopr.pyomo.

........

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


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

........

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


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

........

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


Forgot some debug output in the last commit.

........

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


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

........

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


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

........

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


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

........

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


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

........

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


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

........

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


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

........

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


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

........

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


NFC: remove EOL whitespace

........

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


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

........

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


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

........

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


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

........

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


NFC: EOL whitespace removal

........

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


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

........

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


Minor touch-ups relating to traceback facilities.

........

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