source: coopr.pysp/stable/coopr/pysp/ef.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: 30.3 KB
Line 
1import pyutilib
2import sys
3import os
4import time
5import traceback
6import copy
7import gc
8
9from coopr.pysp.scenariotree import *
10from coopr.pysp.convergence import *
11from coopr.pysp.ph import *
12from coopr.pysp.phutils import *
13
14from coopr.pyomo.base import *
15from coopr.pyomo.io import *
16
17from coopr.pyomo.base.var import _VarValue, _VarBase
18
19from coopr.opt.results.solution import Solution
20
21#
22# brain-dead utility for determing if there is a binary to write in the
23# composite model - need to know this, because CPLEX doesn't like empty
24# binary blocks in the LP file.
25#
26
27def binaries_present(master_model, scenario_instances):
28
29   # check the master model first.
30   for var in master_model.active_components(Var).values():
31      if isinstance(var.domain, BooleanSet):
32         return True
33
34   # scan the scenario instances next.
35   for scenario_name in scenario_instances.keys():
36      scenario_instance = scenario_instances[scenario_name]
37      for var in scenario_instance.active_components(Var).values():
38         if isinstance(var.domain, BooleanSet):
39            return True
40
41   return False
42
43#
44# brain-dead utility for determing if there is a binary to write in the
45# composite model - need to know this, because CPLEX doesn't like empty
46# integer blocks in the LP file.
47#
48
49def integers_present(master_model, scenario_instances):
50
51   # check the master model first.
52   for var in master_model.active_components(Var).values():
53      if isinstance(var.domain, IntegerSet):
54         return True
55
56   # scan the scenario instances next.
57   for scenario_name in scenario_instances.keys():
58      scenario_instance = scenario_instances[scenario_name]
59      for var in scenario_instance.active_components(Var).values():
60         if isinstance(var.domain, IntegerSet):
61            return True
62
63   return False
64
65#
66# a routine to create the extensive form, given an input scenario tree and instances.
67# IMPT: unlike scenario instances, the extensive form instance is *not* self-contained.
68#       in particular, it has binding constraints that cross the binding instance and
69#       the scenario instances. it is up to the caller to keep track of which scenario
70#       instances are associated with the extensive form. this might be something we
71#       encapsulate at some later time.
72# NOTE: if cvar terms are generated, then the input scenario tree is modified accordingly,
73#       i.e., with the addition of the "eta" variable at the root node and the excess
74#       variables at (for lack of a better place - they are per-scenario, but are not
75#       blended) the second stage.
76#
77
78def create_ef_instance(scenario_tree, scenario_instances,
79                       verbose_output = False,
80                       generate_weighted_cvar = False, cvar_weight = None, risk_alpha = None):
81
82   #
83   # validate cvar options, if specified.
84   #
85   if generate_weighted_cvar is True:
86      if (cvar_weight is None) or (cvar_weight < 0.0):
87         raise RuntimeError, "Weight of CVaR term must be >= 0.0 - value supplied="+str(cvar_weight)
88      if (risk_alpha is None) or (risk_alpha <= 0.0) or (risk_alpha >= 1.0):
89         raise RuntimeError, "CVaR risk alpha must be between 0 and 1, exclusive - value supplied="+str(risk_alpha)
90
91      if verbose_output is True:
92         print "Writing CVaR weighted objective"
93         print "CVaR term weight="+str(cvar_weight)
94         print "CVaR alpha="+str(risk_alpha)
95         print ""
96
97      # update the scenario tree with cvar-specific variable names, so
98      # they will get cloned accordingly in the master instance.
99      first_stage = scenario_tree._stages[0]
100      second_stage = scenario_tree._stages[1]
101      root_node = first_stage._tree_nodes[0]
102
103      # NOTE: because we currently don't have access to the reference
104      #       instance in this method, temporarily (and largely orphaned)
105      #       variables are constructed to supply to the scenario tree.
106      #       this decision should be ultimately revisited.
107      cvar_eta_variable_name = "CVAR_ETA"
108      cvar_eta_variable = Var(name=cvar_eta_variable_name)
109      cvar_eta_variable.construct()               
110
111      first_stage.add_variable(cvar_eta_variable, "*", [None])
112
113      cvar_excess_variable_name = "CVAR_EXCESS"
114      cvar_excess_variable = Var(name=cvar_excess_variable_name)
115      cvar_excess_variable.construct()
116
117      second_stage.add_variable(cvar_excess_variable, "*", [None])
118
119      # create the eta and excess variable on a per-scenario basis,
120      # in addition to the constraint relating to the two.
121      for scenario_name, scenario_instance in scenario_instances.items():
122
123         cvar_excess_variable_name = "CVAR_EXCESS"
124         cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals)
125         cvar_excess_variable.construct()         
126         setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable)
127
128         cvar_eta_variable_name = "CVAR_ETA"
129         cvar_eta_variable = Var(name=cvar_eta_variable_name)
130         cvar_eta_variable.construct()         
131         setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable)
132
133         compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS"
134         compute_excess_constraint = Constraint(name=compute_excess_constraint_name)
135         compute_excess_expression = cvar_excess_variable
136         for node in scenario_tree._scenario_map[scenario_name]._node_list:
137            (cost_variable, cost_variable_idx) = node._stage._cost_variable
138            compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx]
139         compute_excess_expression += cvar_eta_variable
140         compute_excess_constraint.add(None, (0.0, compute_excess_expression, None))
141         compute_excess_constraint._model = scenario_instance
142         setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint)
143
144         # re-process the scenario instance due to the newly added constraints/variables associated
145         # with CVaR. a complete preprocess is overkill, of course - the correct approach would be
146         # to just preprocess those newly added variables and constraints.
147         scenario_instance.preprocess()
148
149   #
150   # create the new and empty binding instance.
151   #
152
153   binding_instance = Model()
154   binding_instance.name = "MASTER"
155
156   # walk the scenario tree - create variables representing the common values for all scenarios
157   # associated with that node, along with equality constraints to enforce non-anticipativity.
158   # also create expected cost variables for each node, to be computed via constraints/objectives
159   # defined in a subsequent pass. master variables are created for all nodes but those in the
160   # last stage. expected cost variables are, for no particularly good reason other than easy
161   # coding, created for nodes in all stages.
162   if verbose_output is True:
163      print "Creating variables for master binding instance"
164
165   for stage in scenario_tree._stages:
166
167      # first loop is to create master (blended) variables across all stages but the last.
168      for (stage_variable, index_template, stage_variable_indices) in stage._variables:
169
170         if verbose_output is True:
171            print "Creating master variable and blending constraints for decision variable="+stage_variable.name+", indices="+str(index_template)
172
173         for tree_node in stage._tree_nodes:
174
175            if stage != scenario_tree._stages[-1]:     
176
177               master_variable_name = stage_variable.name               
178
179               # because there may be a single stage variable and multiple indices, check
180               # for the existence of the variable at this node - if you don't, you'll
181               # inadvertently over-write what was there previously!
182               master_variable = None
183               try:
184                  master_variable = getattr(binding_instance, master_variable_name)
185               except:
186                  new_master_variable_index = stage_variable._index
187                  new_master_variable = None
188                  if (len(new_master_variable_index) is 1) and (None in new_master_variable_index):
189                     new_master_variable = Var(name=stage_variable.name)
190                  else:
191                     new_master_variable = Var(new_master_variable_index, name=stage_variable.name)
192                  new_master_variable.construct()
193                  new_master_variable._model = binding_instance
194                  setattr(binding_instance, master_variable_name, new_master_variable)
195
196                  master_variable = new_master_variable
197
198               # TBD: we should create an indexed variable here, and then add entries within the loop.
199               for index in stage_variable_indices:
200
201                  is_used = True # until proven otherwise                     
202                  for scenario in tree_node._scenarios:
203                     instance = scenario_instances[scenario._name]
204                     if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
205                        is_used = False
206
207                  is_fixed = False # until proven otherwise
208                  for scenario in tree_node._scenarios:
209                     instance = scenario_instances[scenario._name]
210                     if getattr(instance,stage_variable.name)[index].fixed is True:
211                        is_fixed = True
212
213                  if (is_used is True) and (is_fixed is False):
214                           
215                     for scenario in tree_node._scenarios:
216                        scenario_instance = scenario_instances[scenario._name]
217                        scenario_variable = getattr(scenario_instance, stage_variable.name)
218                        new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
219                        new_constraint = Constraint(name=new_constraint_name)
220                        new_expr = master_variable[index] - scenario_variable[index]
221                        new_constraint.add(None, (0.0, new_expr, 0.0))
222                        new_constraint._model = binding_instance
223                        setattr(binding_instance, new_constraint_name, new_constraint)
224
225      # the second loop is for creating the stage cost variable in each tree node.
226      for tree_node in stage._tree_nodes:                       
227
228         # create a variable to represent the expected cost at this node -
229         # the constraint to compute this comes later.
230         expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
231         expected_cost_variable = Var(name=expected_cost_variable_name)
232         expected_cost_variable._model = binding_instance
233         setattr(binding_instance, expected_cost_variable_name, expected_cost_variable)
234
235   if generate_weighted_cvar is True:
236
237      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
238      cvar_cost_variable = Var(name=cvar_cost_variable_name)
239      cvar_cost_variable.construct()           
240      setattr(binding_instance, cvar_cost_variable_name, cvar_cost_variable)
241
242   binding_instance.preprocess()
243
244   # ditto above for the (non-expected) cost variable.
245   for stage in scenario_tree._stages:
246
247      (cost_variable,cost_variable_index) = stage._cost_variable
248
249      if verbose_output is True:
250         print "Creating master variable and blending constraints for cost variable="+cost_variable.name+", index="+str(cost_variable_index)
251
252      for tree_node in stage._tree_nodes:
253
254         # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)         
255
256         # this is undoubtedly wasteful, in that a cost variable
257         # for each tree node is created with *all* indices.         
258         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
259         new_cost_variable_index = cost_variable._index
260         new_cost_variable = None
261         if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index):
262            new_cost_variable = Var(name=new_cost_variable_name)
263         else:
264            new_cost_variable = Var(new_cost_variable_index, name=new_cost_variable_name)
265         new_cost_variable.construct()
266         new_cost_variable._model = binding_instance
267         setattr(binding_instance, new_cost_variable_name, new_cost_variable)         
268
269         # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
270         new_cost_variable[cost_variable_index].var = new_cost_variable
271         if cost_variable_index is not None:
272            # if the variable index is None, the variable is derived from a VarValue, so the
273            # name gets updated automagically.
274            new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
275
276         for scenario in tree_node._scenarios:
277
278            scenario_instance = scenario_instances[scenario._name]
279            scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
280            new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
281            new_constraint = Constraint(name=new_constraint_name)
282            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
283            new_constraint.add(None, (0.0, new_expr, 0.0))
284            new_constraint._model = binding_instance
285            setattr(binding_instance, new_constraint_name, new_constraint)
286
287   # create the constraints for computing the master per-node cost variables,
288   # i.e., the current node cost and the expected cost of the child nodes.
289   # if the root, then the constraint is just the objective.
290   for stage in scenario_tree._stages:
291
292      (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
293
294      for tree_node in stage._tree_nodes:
295
296         node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
297         node_expected_cost_variable = getattr(binding_instance, node_expected_cost_variable_name)
298
299         node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
300         node_cost_variable = getattr(binding_instance, node_cost_variable_name)                       
301           
302         constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
303
304         for child_node in tree_node._children:
305
306            child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
307            child_node_expected_cost_variable = getattr(binding_instance, child_node_expected_cost_variable_name)
308            constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
309
310         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
311         new_constraint = Constraint(name=new_constraint_name)
312         new_constraint.add(None, (0.0, constraint_expr, 0.0))
313         new_constraint._model = binding_instance                     
314         setattr(binding_instance, new_constraint_name, new_constraint)
315
316         if tree_node._parent is None:
317
318            an_instance = scenario_instances[scenario_instances.keys()[0]]
319            an_objective = an_instance.active_components(Objective)
320            opt_sense = an_objective[an_objective.keys()[0]].sense
321
322            opt_expression = node_expected_cost_variable
323
324            if generate_weighted_cvar is True:
325               cvar_cost_variable_name = "CVAR_COST_" + tree_node._name
326               cvar_cost_variable = getattr(binding_instance, cvar_cost_variable_name)
327               if cvar_weight == 0.0:
328                  # if the cvar weight is 0, then we're only doing cvar - no mean.
329                  opt_expression = cvar_cost_variable                             
330               else:
331                  opt_expression += cvar_weight * cvar_cost_variable           
332
333            new_objective = Objective(name="MASTER", sense=opt_sense)
334            new_objective._data[None].expr = opt_expression
335            setattr(binding_instance, "MASTER", new_objective)
336
337   # CVaR requires the addition of a variable per scenario to represent the cost excess,
338   # and a constraint to compute the cost excess relative to eta. we also replicate (following
339   # what we do for node cost variables) an eta variable for each scenario instance, and
340   # require equality with the master eta variable via constraints.
341   if generate_weighted_cvar is True:
342     
343      # add the constraint to compute the master CVaR variable value. iterate
344      # over scenario instances to create the expected excess component first.
345      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
346      cvar_cost_variable = getattr(binding_instance, cvar_cost_variable_name)
347      cvar_eta_variable_name = "CVAR_ETA"
348      cvar_eta_variable = getattr(binding_instance, cvar_eta_variable_name)
349     
350      cvar_cost_expression = cvar_cost_variable - cvar_eta_variable
351     
352      for scenario_name, scenario_instance in scenario_instances.items():
353
354         scenario_probability = scenario_tree._scenario_map[scenario_name]._probability
355
356         scenario_excess_variable_name = "CVAR_EXCESS"
357         scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name)
358
359         cvar_cost_expression = cvar_cost_expression - (scenario_probability * scenario_excess_variable) / (1.0 - risk_alpha)
360
361      compute_cvar_cost_constraint_name = "COMPUTE_CVAR_COST"
362      compute_cvar_cost_constraint = Constraint(name=compute_cvar_cost_constraint_name)
363      compute_cvar_cost_constraint.add(None, (0.0, cvar_cost_expression, 0.0))
364      compute_cvar_cost_constraint._model = binding_instance
365      setattr(binding_instance, compute_cvar_cost_constraint_name, compute_cvar_cost_constraint)
366
367   # always preprocess the binding instance, so that it is ready for solution.
368   binding_instance.preprocess()
369
370   return binding_instance
371
372#
373# write the EF binding instance and all sub-instances. currently only outputs the CPLEX LP file format.
374#
375
376def write_ef(binding_instance, scenario_instances, output_filename):
377
378   # create the output file.
379   problem_writer = cpxlp.ProblemWriter_cpxlp()
380   output_file = open(output_filename,"w")
381
382   problem_writer._output_prefixes = True # we always want prefixes
383
384   ################################################################################################
385   #### WRITE THE MASTER OBJECTIVE ################################################################
386   ################################################################################################
387
388   # write the objective for the master binding instance.
389   problem_writer._output_objectives = True
390   problem_writer._output_constraints = False
391   problem_writer._output_variables = False
392
393   print >>output_file, "\\ Begin objective block for master"
394   problem_writer._print_model_LP(binding_instance, output_file)
395   print >>output_file, "\\ End objective block for master"
396   print >>output_file, ""
397
398   ################################################################################################
399   #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
400   ################################################################################################
401
402   print >>output_file, "s.t."
403   print >>output_file, ""
404   
405   problem_writer._output_objectives = False
406   problem_writer._output_constraints = True
407   problem_writer._output_variables = False
408
409   print >>output_file, "\\ Begin constraint block for master"
410   problem_writer._print_model_LP(binding_instance, output_file)
411   print >>output_file, "\\ End constraint block for master",
412   print >>output_file, ""
413
414   for scenario_name in scenario_instances.keys():
415      instance = scenario_instances[scenario_name]
416      print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
417      problem_writer._print_model_LP(instance, output_file)
418      print >>output_file, "\\ End constraint block for scenario",scenario_name
419      print >>output_file, ""
420
421   ################################################################################################
422   #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
423   ################################################################################################
424
425   # write the variables for the master binding instance, and then for each scenario.
426   print >>output_file, "bounds"
427   print >>output_file, ""
428   
429   problem_writer._output_objectives = False
430   problem_writer._output_constraints = False
431   problem_writer._output_variables = True
432
433   # first step: write variable bounds
434
435   problem_writer._output_continuous_variables = True
436   problem_writer._output_integer_variables = False
437   problem_writer._output_binary_variables = False
438
439   print >>output_file, "\\ Begin variable bounds block for master"
440   problem_writer._print_model_LP(binding_instance, output_file)
441   print >>output_file, "\\ End variable bounds block for master"
442   print >>output_file, ""
443   
444   for scenario_name in scenario_instances.keys():
445      instance = scenario_instances[scenario_name]
446      print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name
447      problem_writer._print_model_LP(instance, output_file)
448      print >>output_file, "\\ End variable bounds block for scenario",scenario_name
449      print >>output_file, ""
450
451   # second step: write integer indicators.
452
453   problem_writer._output_continuous_variables = False
454   problem_writer._output_integer_variables = True
455
456   if integers_present(binding_instance, scenario_instances) is True:
457
458      print >>output_file, "general"
459      print >>output_file, ""
460
461      print >>output_file, "\\ Begin discrete variable block for master"
462      problem_writer._print_model_LP(binding_instance, output_file)
463      print >>output_file, "\\ End discrete variable block for master"
464      print >>output_file, ""
465   
466      for scenario_name in scenario_instances.keys():
467         instance = scenario_instances[scenario_name]
468         print >>output_file, "\\ Begin discrete variable block for scenario",scenario_name
469         problem_writer._print_model_LP(instance, output_file)
470         print >>output_file, "\\ End discrete variable block for scenario",scenario_name
471         print >>output_file, ""
472
473   # third step: write binary indicators.
474
475   problem_writer._output_integer_variables = False
476   problem_writer._output_binary_variables = True
477
478   if binaries_present(binding_instance, scenario_instances) is True:
479
480      print >>output_file, "binary"
481      print >>output_file, ""
482
483      print >>output_file, "\\ Begin binary variable block for master"
484      problem_writer._print_model_LP(binding_instance, output_file)
485      print >>output_file, "\\ End binary variable block for master"
486      print >>output_file, ""
487   
488      for scenario_name in scenario_instances.keys():
489         instance = scenario_instances[scenario_name]
490         print >>output_file, "\\ Begin binary variable block for scenario",scenario_name
491         problem_writer._print_model_LP(instance, output_file)
492         print >>output_file, "\\ End integer binary block for scenario",scenario_name
493         print >>output_file, ""
494
495   # wrap up.
496   print >>output_file, "end"
497
498   # clean up.
499   output_file.close()
500   
501
502#
503# the main extensive-form writer routine - including read of scenarios/etc.
504# returns a triple consisting of the scenario tree, master binding instance, and scenario instance map
505#
506
507def write_ef_from_scratch(model_directory, instance_directory, output_filename,
508                          verbose_output, linearize,
509                          generate_weighted_cvar, cvar_weight, risk_alpha):
510
511   start_time = time.time()
512
513   scenario_data_directory_name = instance_directory
514
515   print "Loading scenario and instance data"
516
517   #
518   # create and populate the core model
519   #
520   master_scenario_model = None
521   master_scenario_instance = None
522
523   if verbose_output:
524      print "Constructing reference model and instance"
525   
526   try:
527     
528      reference_model_filename = model_directory+os.sep+"ReferenceModel.py"   
529      modelimport = pyutilib.misc.import_file(reference_model_filename)
530      if "model" not in dir(modelimport):
531         print ""
532         print "Exiting ef module: No 'model' object created in module "+reference_model_filename
533         sys.exit(0)
534      if modelimport.model is None:
535         print ""
536         print "Exiting ef module: 'model' object equals 'None' in module "+reference_model_filename
537         sys.exit(0)
538   
539      master_scenario_model = modelimport.model
540
541   except IOError:
542     
543      print "***ERROR: Failed to load scenario reference model from file="+reference_model_filename
544      return None, None, None
545
546   try:
547     
548      reference_scenario_filename = instance_directory+os.sep+"ReferenceModel.dat"
549      master_scenario_instance = master_scenario_model.create(reference_scenario_filename)
550     
551   except IOError:
552     
553      print "***ERROR: Failed to load scenario reference instance data from file="+reference_scenario_filename
554      return None, None, None           
555
556   #
557   # create and populate the scenario tree model
558   #
559
560   from coopr.pysp.util.scenariomodels import scenario_tree_model
561
562   if verbose_output:
563      print "Constructing scenario tree instance"
564
565   scenario_tree_instance = scenario_tree_model.create(instance_directory+os.sep+"ScenarioStructure.dat")
566
567   #
568   # construct the scenario tree
569   #
570   if verbose_output:
571      print "Constructing scenario tree object"
572   
573   scenario_tree = ScenarioTree(scenarioinstance=master_scenario_instance,
574                                scenariotreeinstance=scenario_tree_instance)
575
576   #
577   # print the input tree for validation/information purposes.
578   #
579   if verbose_output is True:
580      scenario_tree.pprint()
581
582   #
583   # validate the tree prior to doing anything serious
584   #
585   if scenario_tree.validate() is False:
586      print "***Scenario tree is invalid****"
587      sys.exit(1)
588   else:
589      if verbose_output is True:
590         print "Scenario tree is valid!"
591
592   #
593   # construct instances for each scenario
594   #
595
596   # the construction of instances takes little overhead in terms of
597   # memory potentially lost in the garbage-collection sense (mainly
598   # only that due to parsing and instance simplification/prep-processing).
599   # to speed things along, disable garbage collection if it enabled in
600   # the first place through the instance construction process.
601   # IDEA: If this becomes too much for truly large numbers of scenarios,
602   #       we could manually collect every time X instances have been created.
603
604   re_enable_gc = False
605   if gc.isenabled() is True:
606      re_enable_gc = True
607      gc.disable()
608
609   scenario_instances = {}
610   
611   if scenario_tree._scenario_based_data == 1:
612      if verbose_output is True:
613         print "Scenario-based instance initialization enabled"
614   else:
615      if verbose_output is True:
616         print "Node-based instance initialization enabled"
617         
618   for scenario in scenario_tree._scenarios:
619
620      scenario_instance = construct_scenario_instance(scenario_tree,
621                                                      instance_directory,
622                                                      scenario._name,
623                                                      master_scenario_model,
624                                                      verbose=verbose_output,
625                                                      preprocess=True,
626                                                      linearize=linearize)
627
628      scenario_instances[scenario._name] = scenario_instance
629      # name each instance with the scenario name, so the prefixes in the EF make sense.
630      scenario_instance.name = scenario._name
631
632   if re_enable_gc is True:
633      gc.enable()
634
635   print "Creating extensive form binding instance"
636
637   binding_instance = create_ef_instance(scenario_tree, scenario_instances,
638                                         verbose_output = verbose_output,
639                                         generate_weighted_cvar = generate_weighted_cvar,
640                                         cvar_weight = cvar_weight,
641                                         risk_alpha = risk_alpha)
642
643   print "Starting to write extensive form"
644
645   write_ef(binding_instance, scenario_instances, output_filename)   
646
647   print "Output file written to file=",output_filename
648
649   end_time = time.time()
650
651   print "Total execution time=%8.2f seconds" %(end_time - start_time)
652
653   return scenario_tree, binding_instance, scenario_instances
654
655#
656# does what it says, with the added functionality of returning the master binding instance.
657#
658
659def create_and_write_ef(scenario_tree, scenario_instances, output_filename):
660
661   start_time = time.time()
662
663   binding_instance = create_ef_instance(scenario_tree, scenario_instances)
664
665   print "Starting to write extensive form"
666
667   write_ef(binding_instance, scenario_instances, output_filename)
668
669   print "Output file written to file=",output_filename
670
671   end_time = time.time()
672
673   print "Total execution time=%8.2f seconds" %(end_time - start_time)
674
675   return binding_instance
676
677#
678# a utility to load an EF solution into the corresponding instances.
679#
680
681def load_ef_solution(ef_results, binding_instance, scenario_instances):
682
683   # type is coopr.opt.results.results.SolverResults
684   if len(ef_results.solution) == 0:
685      raise RuntimeError, "Method load_ef_solution invoked with results object containing no solutions!"
686   elif len(ef_results.solution) > 1:
687      raise RuntimeError, "Method load_ef_solution invoked with results object containing more than one solution!"
688   
689   # type is coopr.opt.results.solution.Solution
690   solution = ef_results.solution[0]
691
692   # shotgun the ef solution into individual solutions for the binding and scenario instances.
693   sub_solutions = {} # map between instance name and the corresponding Solution
694   sub_solutions[binding_instance.name] = Solution()
695   for scenario_name, scenario in scenario_instances.items():
696      sub_solutions[scenario_name] = Solution()
697
698   for key, attr_dictionary in solution.variable.items():
699      tokens = string.split(key, '_', 1)
700      instance_name = tokens[0]
701      variable_name = tokens[1]
702      subsolution_variable = sub_solutions[instance_name].variable[variable_name]
703      for attr_name in attr_dictionary.keys():
704         attr_value = attr_dictionary[attr_name]
705         setattr(subsolution_variable, attr_name, attr_value)
706
707   # load the sub-solutions into the appropriate instances.
708   for instance_name, sub_solution in sub_solutions.items():
709      if instance_name == binding_instance.name:
710         binding_instance.load(sub_solution)
711      else:
712         scenario_instances[instance_name].load(sub_solution)
Note: See TracBrowser for help on using the repository browser.