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

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

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

........

r3120 | jwatson | 2010-10-18 21:11:21 -0600 (Mon, 18 Oct 2010) | 3 lines


Adding PySP options for linearizing expressions.

........

r3134 | jwatson | 2010-10-21 15:45:49 -0600 (Thu, 21 Oct 2010) | 3 lines


Bug fix associated with linear expressions.

........

r3137 | khunter | 2010-10-22 10:19:44 -0600 (Fri, 22 Oct 2010) | 2 lines


Add name as contributor.

........

r3138 | jwatson | 2010-10-22 14:11:43 -0600 (Fri, 22 Oct 2010) | 3 lines


Added "--simplify-expressions" option to runph, in order to eliminate the memory and time costs associated with simplifying expressions (e.g., in formulation of the PH objective) for which simpification is very unlikely to help. By default, expression simplification is disabled. This may cause issues with certain writers, e.g., NL - which is why I have retained the option.

........

r3139 | jwatson | 2010-10-22 15:20:49 -0600 (Fri, 22 Oct 2010) | 3 lines


When forming the PH linear terms for all variables and the proximal terms for binary variables, use value() to immediately extract the fixed value of the underlying parameters. This speeds up the associated expression tree construction time, which was non-trivial (5% of the total run-time, for example). Because we re-form the expressions each iteration, this is OK.

........

r3143 | jwatson | 2010-10-22 22:07:07 -0600 (Fri, 22 Oct 2010) | 3 lines


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

........

r3145 | jwatson | 2010-10-22 22:21:06 -0600 (Fri, 22 Oct 2010) | 3 lines


Update of PySP baseline output files, due to small numerical differences (4th significant digit) associated with some code mods I recently performed in the course of optimization.

........

r3146 | jwatson | 2010-10-23 13:27:38 -0600 (Sat, 23 Oct 2010) | 3 lines


Speed improvements in weight update.

........

r3147 | jwatson | 2010-10-23 13:38:52 -0600 (Sat, 23 Oct 2010) | 3 lines


More speed improvements to PH weight computation procedure.

........

r3152 | jwatson | 2010-10-24 13:20:03 -0600 (Sun, 24 Oct 2010) | 3 lines


Further code optimizations to PH weight and statistic update routines.

........

r3155 | jwatson | 2010-10-24 22:15:29 -0600 (Sun, 24 Oct 2010) | 3 lines


Added --linearize-expression and associated processing to runef PySP script.

........

r3156 | jwatson | 2010-10-24 22:19:44 -0600 (Sun, 24 Oct 2010) | 3 lines


Correcting shift from "integer" to "general" variable blocks when writing the LP file for an extensive form.

........

r3165 | jwatson | 2010-10-26 18:46:24 -0600 (Tue, 26 Oct 2010) | 3 lines


Update of PySP CHANGELOG for 2.4 release.

........

r3168 | jwatson | 2010-10-26 20:12:23 -0600 (Tue, 26 Oct 2010) | 3 lines


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.