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

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

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

........

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


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

........

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


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

........

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


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

........

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


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

........

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


Removing some older PySP test problems from the repository

........

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


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

........

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


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


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


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

........

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


Corrected issue with cvar generation introduced a while back.

........

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


Adding PySP extensive form tests.

........

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


Auxmenting the filter

........

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


Adding further diagnostics to the filter.

........

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


Another attempt to fix this filter...

........

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


Simplifying filter.


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

........

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


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

........

File size: 13.8 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 string
12import traceback
13import os
14
15from coopr.pyomo import *
16
17from coopr.pyomo.expr.linear_repn import linearize_model_expressions
18
19#
20# a simple utility function to pretty-print an index tuple into a [x,y] form string.
21#
22
23def indexToString(index):
24
25   # if the input type is a string or an int, then this isn't a tuple!
26   if isinstance(index, (str, int)):
27      return "["+str(index)+"]"
28
29   result = "["
30   for i in range(0,len(index)):
31      result += str(index[i])
32      if i != len(index) - 1:
33         result += ","
34   result += "]"
35   return result
36
37#
38# a simple utility to determine if a variable name contains an index specification.
39# in other words, is the reference to a complete variable (e.g., "foo") - which may
40# or may not be indexed - or a specific index or set of indices (e.g., "foo[1]" or
41# or "foo[1,*]".
42#
43
44def isVariableNameIndexed(variable_name):
45
46   left_bracket_count = variable_name.count('[')
47   right_bracket_count = variable_name.count(']')
48
49   if (left_bracket_count == 1) and (right_bracket_count == 1):
50      return True
51   elif (left_bracket_count == 1) or (right_bracket_count == 1):
52      raise ValueError, "Illegally formed variable name="+variable_name+"; if indexed, variable names must contain matching left and right brackets"
53   else:
54      return False
55
56#
57# takes a string indexed of the form "('foo', 'bar')" and returns a proper tuple ('foo','bar')
58#
59
60def tupleizeIndexString(index_string):
61
62   index_string=index_string.lstrip('(')
63   index_string=index_string.rstrip(')')
64   pieces = index_string.split(',')
65   return_index = ()
66   for piece in pieces:
67      piece = string.strip(piece)
68      piece = piece.lstrip('\'')
69      piece = piece.rstrip('\'')
70      transformed_component = None
71      try:
72         transformed_component = int(piece)
73      except ValueError:
74         transformed_component = piece
75      return_index = return_index + (transformed_component,)
76
77   # IMPT: if the tuple is a singleton, return the element itself.
78   if len(return_index) == 1:
79      return return_index[0]
80   else:
81      return return_index
82
83#
84# related to above, extract the index from the variable name.
85# will throw an exception if the variable name isn't indexed.
86# the returned variable name is a string, while the returned
87# index is a tuple. integer values are converted to integers
88# if the conversion works!
89#
90
91def extractVariableNameAndIndex(variable_name):
92
93   if isVariableNameIndexed(variable_name) is False:
94      raise ValueError, "Non-indexed variable name passed to function extractVariableNameAndIndex()"
95
96   pieces = variable_name.split('[')
97   name = string.strip(pieces[0])
98   full_index = pieces[1].rstrip(']')
99
100   # even nested tuples in pyomo are "flattened" into
101   # one-dimensional tuples. to accomplish flattening
102   # replace all parens in the string with commas and
103   # proceed with the split.
104   full_index=string.translate(full_index, string.maketrans("",""), "()")
105   indices = full_index.split(',')
106
107   return_index = ()
108
109   for index in indices:
110
111      # unlikely, but strip white-space from the string.
112      index=string.strip(index)
113
114      # if the tuple contains nested tuples, then the nested
115      # tuples have single quotes - "'" characters - around
116      # strings. remove these, as otherwise you have an
117      # illegal index.
118      index=string.translate(index, string.maketrans("",""), "\'")
119
120      # if the index is an integer, make it one!
121      transformed_index = None
122      try:
123         transformed_index = int(index)
124      except ValueError:
125         transformed_index = index
126      return_index = return_index + (transformed_index,)
127
128   # IMPT: if the tuple is a singleton, return the element itself.
129   if len(return_index) == 1:
130      return name, return_index[0]
131   else:
132      return name, return_index
133
134#
135# determine if the input index is an instance of the template,
136# which may or may not contain wildcards.
137#
138
139def indexMatchesTemplate(index, index_template):
140
141   # if the input index is not a tuple, make it one.
142   # ditto with the index template. one-dimensional
143   # indices in pyomo are not tuples, but anything
144   # else is.
145
146   if type(index) != tuple:
147      index = (index,)
148   if type(index_template) != tuple:
149      index_template = (index_template,)
150
151   if len(index) != len(index_template):
152      return False
153
154   for i in range(0,len(index_template)):
155      if index_template[i] == '*':
156         # anything matches
157         pass
158      else:
159         if index_template[i] != index[i]:
160            return False
161
162   return True
163
164#
165# given a variable (the real object, not the name) and an index,
166# "shotgun" the index and see which variable indices match the
167# input index. the cardinality could be > 1 if slices are
168# specified, e.g., [*,1].
169#
170
171def extractVariableIndices(variable, index_template):
172
173   result = []
174
175   for index in variable._index:
176
177      if indexMatchesTemplate(index, index_template) is True:
178         result.append(index)
179
180   return result
181
182#
183# construct a scenario instance!
184#
185def construct_scenario_instance(scenario_tree_instance,
186                                scenario_data_directory_name,
187                                scenario_name,
188                                reference_model,
189                                verbose,
190                                preprocess = True,
191                                linearize = True,
192                                simplify = True):
193
194   if verbose is True:
195      if scenario_tree_instance._scenario_based_data == 1:
196         print "Scenario-based instance initialization enabled"
197      else:
198         print "Node-based instance initialization enabled"
199
200   scenario = scenario_tree_instance.get_scenario(scenario_name)
201   scenario_instance = None
202
203   if verbose is True:
204      print "Creating instance for scenario=" + scenario._name
205
206   try:
207      if scenario_tree_instance._scenario_based_data == 1:
208         scenario_data_filename = scenario_data_directory_name + os.sep + scenario._name + ".dat"
209         if verbose is True:
210            print "Data for scenario=" + scenario._name + " loads from file=" + scenario_data_filename
211         scenario_instance = reference_model.create(scenario_data_filename, simplify=simplify)
212      else:
213         scenario_instance = reference_model.clone()
214
215         data_files = []
216         current_node = scenario._leaf_node
217         while current_node is not None:
218            node_data_filename = scenario_data_directory_name + os.sep + current_node._name + ".dat"
219            data_files.append(node_data_filename)
220            current_node = current_node._parent         
221
222         # to make sure we read from root node to leaf node
223         data_files.reverse() 
224
225         scenario_data = ModelData()
226         for data_file in data_files:
227            if verbose is True:
228               print "Node data for scenario=" + scenario._name + " partially loading from file=" + data_file
229            scenario_data.add(data_file)
230
231         scenario_data.read(model=scenario_instance)
232         scenario_instance.load(scenario_data, simplify=simplify)
233         if preprocess is True:
234            scenario_instance.preprocess()         
235   except:
236      print "Encountered exception in model instance creation - traceback:"
237      traceback.print_exc()
238      raise RuntimeError, "Failed to create model instance for scenario=" + scenario._name
239
240   if linearize is True:
241      # IMPT: The model *must* be preprocessed in order for linearization to work. This is because
242      #       linearization relies on the canonical expression representation extraction routine,
243      #       which in turn relies on variables being identified/categorized (e.g., into "Used").
244      if preprocess is False:
245         scenario_instance.preprocess()
246      linearize_model_expressions(scenario_instance)
247
248   return scenario_instance
249
250
251# creates all PH parameters for a problem instance, given a scenario tree
252# (to identify the relevant variables), a default rho (simply for initialization),
253# and a boolean indicating if quadratic penalty terms are to be linearized.
254# returns a list of any created variables, specifically when linearizing -
255# this is useful to clean-up reporting.
256
257def create_ph_parameters(instance, scenario_tree, default_rho, linearizing_penalty_terms):
258
259   new_penalty_variable_names = []
260
261   # first, gather all unique variables referenced in any stage
262   # other than the last, independent of specific indices. this
263   # "gather" step is currently required because we're being lazy
264   # in terms of index management in the scenario tree - which
265   # should really be done in terms of sets of indices.
266   # NOTE: technically, the "instance variables" aren't really references
267   # to the variable in the instance - instead, the reference model. this
268   # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
269   instance_variables = {}
270
271   for stage in scenario_tree._stages[:-1]:
272
273      for (reference_variable, index_template) in stage._variables:
274
275         if reference_variable.name not in instance_variables.keys():
276
277            instance_variables[reference_variable.name] = reference_variable
278
279   # for each blended variable, create a corresponding ph weight and average parameter in the instance.
280   # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor
281   # in the grand scheme of things.
282
283   for (variable_name, reference_variable) in instance_variables.items():
284
285      # PH WEIGHT
286
287      new_w_index = reference_variable._index
288      new_w_parameter_name = "PHWEIGHT_"+reference_variable.name
289      # this bit of ugliness is due to Pyomo not correctly handling the Param construction
290      # case when the supplied index set consists strictly of None, i.e., the source variable
291      # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
292      new_w_parameter = None
293      if (len(new_w_index) is 1) and (None in new_w_index):
294         new_w_parameter = Param(name=new_w_parameter_name, mutable=True, nochecking=True)
295      else:
296         new_w_parameter = Param(new_w_index, name=new_w_parameter_name, mutable=True, nochecking=True)
297      setattr(instance,new_w_parameter_name,new_w_parameter)
298
299      # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference
300      # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not
301      # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the
302      # issue in any case for now. NOTE: I'm not sure this is still the case in Pyomo.
303      for index in new_w_index:
304         new_w_parameter[index] = 0.0
305
306      # PH AVG
307
308      new_avg_index = reference_variable._index
309      new_avg_parameter_name = "PHAVG_"+reference_variable.name
310      new_avg_parameter = None
311      if (len(new_avg_index) is 1) and (None in new_avg_index):
312         new_avg_parameter = Param(name=new_avg_parameter_name, mutable=True, nochecking=True)
313      else:
314         new_avg_parameter = Param(new_avg_index, name=new_avg_parameter_name, mutable=True, nochecking=True)
315      setattr(instance,new_avg_parameter_name,new_avg_parameter)
316
317      for index in new_avg_index:
318         new_avg_parameter[index] = 0.0
319
320      # PH RHO
321
322      new_rho_index = reference_variable._index
323      new_rho_parameter_name = "PHRHO_"+reference_variable.name
324      new_rho_parameter = None
325      if (len(new_avg_index) is 1) and (None in new_avg_index):
326         new_rho_parameter = Param(name=new_rho_parameter_name, mutable=True, nochecking=True)
327      else:
328         new_rho_parameter = Param(new_rho_index, name=new_rho_parameter_name, mutable=True, nochecking=True)
329      setattr(instance,new_rho_parameter_name,new_rho_parameter)
330
331      for index in new_rho_index:
332         new_rho_parameter[index] = default_rho
333
334      # PH LINEARIZED PENALTY TERM
335
336      if linearizing_penalty_terms > 0:
337
338         new_penalty_term_variable_index = reference_variable._index
339         new_penalty_term_variable_name = "PHQUADPENALTY_"+reference_variable.name
340         # this is a quadratic penalty term - the lower bound is 0!
341         new_penalty_term_variable = None
342         if (len(new_avg_index) is 1) and (None in new_avg_index):
343            new_penalty_term_variable = Var(name=new_penalty_term_variable_name, bounds=(0.0,None))
344         else:
345            new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None))
346         new_penalty_term_variable.construct()
347         setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable)
348         new_penalty_variable_names.append(new_penalty_term_variable_name)
349
350      # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY.
351
352      # also controls whether weight updates proceed at any iteration.
353
354      new_blend_index = reference_variable._index
355      new_blend_parameter_name = "PHBLEND_"+reference_variable.name
356      new_blend_parameter = None
357      if (len(new_avg_index) is 1) and (None in new_avg_index):
358         new_blend_parameter = Param(name=new_blend_parameter_name, within=Binary, mutable=True, nochecking=True)
359      else:
360         new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary, mutable=True, nochecking=True)
361      setattr(instance,new_blend_parameter_name,new_blend_parameter)
362
363      for index in new_blend_index:
364         new_blend_parameter[index] = 1
365
366   return new_penalty_variable_names
Note: See TracBrowser for help on using the repository browser.