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

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

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

File size: 13.7 KB
Line 
1#  _________________________________________________________________________
2#
3#  Coopr: A COmmon Optimization Python Repository
4#  Copyright (c) 2008 Sandia Corporation.
5#  This software is distributed under the BSD License.
6#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
7#  the U.S. Government retains certain rights in this software.
8#  For more information, see the Coopr README.txt file.
9#  _________________________________________________________________________
10
11import 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, None, "()")
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, None, "\'")
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.