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

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

Modifying more PH parameters to be of the "nochecking" variety. Also slightly improved the performance of the update_variable_statistics PH function.

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