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

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

Completing implementation of optional processing of expression simplification in PySP.

File size: 13.6 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         scenario_data = ModelData()
215         current_node = scenario._leaf_node
216         while current_node is not None:
217            node_data_filename = scenario_data_directory_name + os.sep + current_node._name + ".dat"
218            if verbose is True:
219               print "Node data for scenario=" + scenario._name + " partially loading from file=" + node_data_filename
220            scenario_data.add(node_data_filename)
221            current_node = current_node._parent
222         scenario_data.read(model=scenario_instance)
223         scenario_instance.load(scenario_data, simplify=simplify)
224         if preprocess is True:
225            scenario_instance.preprocess()         
226   except:
227      print "Encountered exception in model instance creation - traceback:"
228      traceback.print_exc()
229      raise RuntimeError, "Failed to create model instance for scenario=" + scenario._name
230
231   if linearize is True:
232      # IMPT: The model *must* be preprocessed in order for linearization to work. This is because
233      #       linearization relies on the canonical expression representation extraction routine,
234      #       which in turn relies on variables being identified/categorized (e.g., into "Used").
235      if preprocess is False:
236         scenario_instance.preprocess()
237      linearize_model_expressions(scenario_instance)
238
239   return scenario_instance
240
241
242# creates all PH parameters for a problem instance, given a scenario tree
243# (to identify the relevant variables), a default rho (simply for initialization),
244# and a boolean indicating if quadratic penalty terms are to be linearized.
245# returns a list of any created variables, specifically when linearizing -
246# this is useful to clean-up reporting.
247
248def create_ph_parameters(instance, scenario_tree, default_rho, linearizing_penalty_terms):
249
250   new_penalty_variable_names = []
251
252   # first, gather all unique variables referenced in any stage
253   # other than the last, independent of specific indices. this
254   # "gather" step is currently required because we're being lazy
255   # in terms of index management in the scenario tree - which
256   # should really be done in terms of sets of indices.
257   # NOTE: technically, the "instance variables" aren't really references
258   # to the variable in the instance - instead, the reference model. this
259   # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
260   instance_variables = {}
261
262   for stage in scenario_tree._stages[:-1]:
263
264      for (reference_variable, index_template, reference_indices) in stage._variables:
265
266         if reference_variable.name not in instance_variables.keys():
267
268            instance_variables[reference_variable.name] = reference_variable
269
270   # for each blended variable, create a corresponding ph weight and average parameter in the instance.
271   # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor
272   # in the grand scheme of things.
273
274   for (variable_name, reference_variable) in instance_variables.items():
275
276      # PH WEIGHT
277
278      new_w_index = reference_variable._index
279      new_w_parameter_name = "PHWEIGHT_"+reference_variable.name
280      # this bit of ugliness is due to Pyomo not correctly handling the Param construction
281      # case when the supplied index set consists strictly of None, i.e., the source variable
282      # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
283      new_w_parameter = None
284      if (len(new_w_index) is 1) and (None in new_w_index):
285         new_w_parameter = Param(name=new_w_parameter_name, mutable=True, nochecking=True)
286      else:
287         new_w_parameter = Param(new_w_index, name=new_w_parameter_name, mutable=True, nochecking=True)
288      setattr(instance,new_w_parameter_name,new_w_parameter)
289
290      # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference
291      # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not
292      # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the
293      # issue in any case for now. NOTE: I'm not sure this is still the case in Pyomo.
294      for index in new_w_index:
295         new_w_parameter[index] = 0.0
296
297      # PH AVG
298
299      new_avg_index = reference_variable._index
300      new_avg_parameter_name = "PHAVG_"+reference_variable.name
301      new_avg_parameter = None
302      if (len(new_avg_index) is 1) and (None in new_avg_index):
303         new_avg_parameter = Param(name=new_avg_parameter_name, mutable=True, nochecking=True)
304      else:
305         new_avg_parameter = Param(new_avg_index, name=new_avg_parameter_name, mutable=True, nochecking=True)
306      setattr(instance,new_avg_parameter_name,new_avg_parameter)
307
308      for index in new_avg_index:
309         new_avg_parameter[index] = 0.0
310
311      # PH RHO
312
313      new_rho_index = reference_variable._index
314      new_rho_parameter_name = "PHRHO_"+reference_variable.name
315      new_rho_parameter = None
316      if (len(new_avg_index) is 1) and (None in new_avg_index):
317         new_rho_parameter = Param(name=new_rho_parameter_name, mutable=True, nochecking=True)
318      else:
319         new_rho_parameter = Param(new_rho_index, name=new_rho_parameter_name, mutable=True, nochecking=True)
320      setattr(instance,new_rho_parameter_name,new_rho_parameter)
321
322      for index in new_rho_index:
323         new_rho_parameter[index] = default_rho
324
325      # PH LINEARIZED PENALTY TERM
326
327      if linearizing_penalty_terms > 0:
328
329         new_penalty_term_variable_index = reference_variable._index
330         new_penalty_term_variable_name = "PHQUADPENALTY_"+reference_variable.name
331         # this is a quadratic penalty term - the lower bound is 0!
332         new_penalty_term_variable = None
333         if (len(new_avg_index) is 1) and (None in new_avg_index):
334            new_penalty_term_variable = Var(name=new_penalty_term_variable_name, bounds=(0.0,None))
335         else:
336            new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None))
337         new_penalty_term_variable.construct()
338         setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable)
339         new_penalty_variable_names.append(new_penalty_term_variable_name)
340
341      # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY.
342
343      # also controls whether weight updates proceed at any iteration.
344
345      new_blend_index = reference_variable._index
346      new_blend_parameter_name = "PHBLEND_"+reference_variable.name
347      new_blend_parameter = None
348      if (len(new_avg_index) is 1) and (None in new_avg_index):
349         new_blend_parameter = Param(name=new_blend_parameter_name, within=Binary, mutable=True, nochecking=True)
350      else:
351         new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary, mutable=True, nochecking=True)
352      setattr(instance,new_blend_parameter_name,new_blend_parameter)
353
354      for index in new_blend_index:
355         new_blend_parameter[index] = 1
356
357   return new_penalty_variable_names
Note: See TracBrowser for help on using the repository browser.