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

Last change on this file since 3104 was 3092, checked in by jwatson, 11 years ago

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

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