source: coopr.pysp/stable/coopr/pysp/ef.py @ 3286

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

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

........

r3201 | jwatson | 2010-10-29 13:18:17 -0600 (Fri, 29 Oct 2010) | 3 lines


Inverting order of .dat files in PySP when loading from a node representation - now root-to-leaf, instead of leaf-to-root. This allows for deeper-in-the-tree nodes to over-write parameter values defined higher in the tree, which is a more "expected" behavior than the converse. The real answer is to throw an exception if a parameter is re-defined, but we're not there yet.

........

r3217 | jwatson | 2010-11-05 11:29:42 -0600 (Fri, 05 Nov 2010) | 3 lines


Various updates to support heteogeneous index sets in PH for different nodes in the scenario tree - more work / testing remains.

........

r3218 | jwatson | 2010-11-05 12:35:51 -0600 (Fri, 05 Nov 2010) | 3 lines


More changes associated with generalizing the PySP index structures from per-stage to per-node.

........

r3220 | jwatson | 2010-11-05 20:28:29 -0600 (Fri, 05 Nov 2010) | 3 lines


Various fixes to the WW PH extension, bringing it in compliance to previous commit changes.

........

r3221 | jwatson | 2010-11-05 20:43:40 -0600 (Fri, 05 Nov 2010) | 1 line


Removing some older PySP test problems from the repository

........

r3222 | jwatson | 2010-11-05 20:55:15 -0600 (Fri, 05 Nov 2010) | 1 line


Moving PySP forestry example to local sandbox, to streamline distribution.

........

r3261 | jwatson | 2010-11-29 15:26:46 -0700 (Mon, 29 Nov 2010) | 9 lines


Adding two options to the runef and runph pysp scripts, to facilitate scenario downsampling - the case where you have a big tree, but you don't want to use it all.


The options are:
--scenario-tree-downsample-fraction=X
--scenario-tree-random-seed


The options are fairly self-explanatory - the only possible nuance is that the downsample fraction is the fraction of scenarios retained.

........

r3263 | jwatson | 2010-12-01 10:47:21 -0700 (Wed, 01 Dec 2010) | 3 lines


Corrected issue with cvar generation introduced a while back.

........

r3264 | jwatson | 2010-12-01 11:16:01 -0700 (Wed, 01 Dec 2010) | 3 lines


Adding PySP extensive form tests.

........

r3265 | wehart | 2010-12-01 13:19:56 -0700 (Wed, 01 Dec 2010) | 2 lines


Auxmenting the filter

........

r3266 | wehart | 2010-12-01 13:58:59 -0700 (Wed, 01 Dec 2010) | 2 lines


Adding further diagnostics to the filter.

........

r3267 | wehart | 2010-12-01 14:12:15 -0700 (Wed, 01 Dec 2010) | 2 lines


Another attempt to fix this filter...

........

r3271 | wehart | 2010-12-01 15:51:23 -0700 (Wed, 01 Dec 2010) | 4 lines


Simplifying filter.


The error was really a Python 2.5 portability issue in EF. :(

........

r3275 | jwatson | 2010-12-02 19:35:52 -0700 (Thu, 02 Dec 2010) | 3 lines


Fixing Python 2.5-related issue with string.translate in phutils.py - using None as the translation table is not allowed in Python 2.5.

........

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