source: coopr.pysp/trunk/coopr/pysp/scenariotree.py @ 3261

Last change on this file since 3261 was 3261, checked in by jwatson, 10 years ago

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.

File size: 42.5 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
33
34      # don't bother copying bounds for variables, as the values stored
35      # here are computed elsewhere - and that entity is responsible for
36      # ensuring feasibility. this also leaves room for specifying infeasible
37      # or partial solutions.
38      new_variable_name = variable.name     
39      new_variable_index = getattr(scenario_instance_map[self._scenarios[0]._name], variable.name)._index
40      new_variable = None
41      if (len(new_variable_index) is 1) and (None in new_variable_index):
42         new_variable = Var(name=new_variable_name)
43      else:
44         new_variable = Var(new_variable_index, name=new_variable_name)
45      new_variable.construct()
46
47      variable_indices = self._variable_indices[variable.name]     
48
49      # by default, deactive all variable values - we're copying the
50      # full index set of a variable, which is necessarily wasteful and
51      # incorrect. then, do a second pass and activate those indicies
52      # that are actually associated with this tree node.
53      for index in new_variable_index:
54         new_variable[index].deactivate()
55
56      for index in variable_indices:
57         new_variable[index].activate()
58
59      self._solutions[new_variable_name] = new_variable
60
61   #
62   # initialize the _solutions attribute of a tree node, from scratch.
63   #
64
65   def _initialize_solution_map(self, scenario_instance_map):
66
67      # clear whatever was there before.
68      self._solutions = {}
69
70      # NOTE: Given the expense, this should really be optional - don't
71      #       construct unless we're actually going to use!
72      # for each variable referenced in the stage, clone the variable
73      # for purposes of storing solutions. we are being wasteful in
74      # terms copying indices that may not be referenced in the stage.
75      # this is something that we might revisit if space/performance
76      # is an issue (space is the most likely issue)
77      for variable, match_template in self._stage._variables:
78         self._update_solution_map(variable, match_template, scenario_instance_map)
79
80      self._solution_map_initialized = True
81
82
83   """ Constructor
84   """
85   def __init__(self, name, conditional_probability, stage, initialize_solution=False, scenario_instance_map={}):
86
87      self._name = name
88      self._stage = stage
89      self._parent = None
90      self._children = [] # a collection of ScenarioTreeNodes
91      self._conditional_probability = conditional_probability # conditional on parent
92      self._scenarios = [] # a collection of all Scenarios passing through this node in the tree
93      self._variable_indices = {} # a map from a variable name to the indices blended at this node.
94
95      # general use statistics for the variables at each node.
96      # each attribute is a map between the variable name and a
97      # parameter (over the same index set) encoding the corresponding
98      # statistic computed over all scenarios for that node. the
99      # parameters are named as the source variable name suffixed
100      # by one of: "NODEMIN", "NODEAVG", and "NODEMAX".
101      # NOTE: the averages are probability_weighted - the min/max
102      #       values are not.
103      # NOTE: the parameter names are basically irrelevant, and the
104      #       convention is assumed to be enforced by whoever populates
105      #       these parameters.
106      self._averages = {}
107      self._minimums = {}
108      self._maximums = {}
109
110      # a flag indicating whether the _solutions attribute has been properly initialized.
111      self._solution_map_initialized = False
112
113      # solution (variable) values for this node. assumed to be distinct
114      # from self._averages, as the latter are not necessarily feasible.
115      # objects in the map are actual pyomo Var instances; keys are
116      # variable names.
117      self._solutions = {}
118
119      if initialize_solution is True:
120         self._initialize_solution_map(scenario_instance_map)
121
122   #
123   # given a set of scenario instances, compute the set of indices being blended
124   # for each variable at this node.
125   #
126
127   def defineVariableIndexSets(self, scenario_instances):
128
129      # find a representative scenario instance belonging to this node. the
130      # first scenario is as good as any.
131      scenario_instance = scenario_instances[self._scenarios[0]._name]
132
133      for reference_variable, match_template in self._stage._variables:
134
135         # the stage variable simply references the variable object in the
136         # reference scenario instance - we need to grab the variable in the
137         # scenario instance, as the index set might be different.
138         variable_name = reference_variable.name
139
140         instance_variable = getattr(scenario_instance, variable_name)
141
142         match_indices = extractVariableIndices(instance_variable, match_template)
143
144         self._variable_indices[variable_name] = match_indices
145
146   #
147   # copies the parameter values values from the _averages attribute
148   # into the _solutions attribute - only for active variable values.
149   #
150
151   def snapshotSolutionFromAverages(self, scenario_instance_map):
152
153      if self._solution_map_initialized is False:
154         self._initialize_solution_map(scenario_instance_map)
155
156      for variable_name, variable in self._solutions.items():
157
158         # try and grab the corresponding averages parameter - if it
159         # doesn't exist, throw an exception.
160         average_parameter = None
161         try:
162            average_parameter = self._averages[variable_name]
163         except:
164            raise RuntimeError, "No averages parameter present on tree node="+self._name+" for variable="+variable_name
165
166         for index in variable._index:
167            if variable[index].active is True:
168               variable[index] = value(average_parameter[index])
169
170   #
171   # computes the solution values from the composite scenario instances at this tree node.
172   # the input scenario_instance_map is a map from scenario name to instance objects.
173   #
174
175   def snapshotSolutionFromInstances(self, scenario_instance_map):
176
177      if self._solution_map_initialized is False:
178         self._initialize_solution_map(scenario_instance_map)
179
180      for variable_name, variable in self._solutions.items():
181         for index in variable._index:
182            if variable[index].active is True:
183               node_probability = 0.0
184               avg = 0.0
185               num_scenarios_with_index = 0
186               for scenario in self._scenarios:
187                  scenario_instance = scenario_instance_map[scenario._name]
188                  node_probability += scenario._probability
189                  scenario_variable = getattr(scenario_instance, variable.name)
190                  if index in scenario_variable:
191                     num_scenarios_with_index = num_scenarios_with_index + 1
192                     var_value = value(getattr(scenario_instance, variable.name)[index])
193                     avg += (scenario._probability * var_value)
194               if num_scenarios_with_index > 0:
195                  variable[index] = avg / node_probability
196
197   #
198   # a utility to compute the cost of the current node plus the expected costs of child nodes.
199   #
200
201   def computeExpectedNodeCost(self, scenario_instance_map):
202
203      # IMPT: This implicitly assumes convergence across the scenarios - if not, garbage results.
204      instance = scenario_instance_map[self._scenarios[0]._name]
205      stage_cost_variable = instance.active_components(Var)[self._stage._cost_variable[0].name][self._stage._cost_variable[1]]
206      if stage_cost_variable.value is not None:
207         my_cost = stage_cost_variable.value
208      else:
209         # depending on the situation (and the model), the stage cost variables might not have values.
210         my_cost = 0.0
211      child_cost = 0.0
212      for child in self._children:
213         child_cost += (child._conditional_probability * child.computeExpectedNodeCost(scenario_instance_map))
214      return my_cost + child_cost
215
216
217class Stage(object):
218
219   """ Constructor
220   """
221   def __init__(self, *args, **kwds):
222
223      self._name = ""
224
225      # a collection of ScenarioTreeNode objects.
226      self._tree_nodes = []
227
228      # a collection of pairs consisting of (1) a reference to a Pyomo model Var object (in the reference scenario instance) and
229      # (2) the original match template string (for output purposes). the specific indices that match belong to the tree node.
230      self._variables = []
231
232      # a tuple consisting of (1) a reference to a pyomo model Var that computes the stage-specific cost and (2) the corresponding
233      # index. the index *is* the sole index in the cost variable, as the cost variable refers to a single variable index.
234      self._cost_variable = (None, None)
235
236   #
237   # add a new variable to the stage, which will include updating the solution maps for each associated ScenarioTreeNode.
238   #
239   def add_variable(self, variable, match_template):
240
241      self._variables.append((variable, match_template))
242
243      for tree_node in self._tree_nodes:
244         tree_node._update_solution_map(variable, match_template)
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.