Changeset 3047


Ignore:
Timestamp:
Sep 26, 2010 12:24:02 PM (11 years ago)
Author:
wehart
Message:

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

........

r3014 | jwatson | 2010-09-15 12:21:07 -0500 (Wed, 15 Sep 2010) | 5 lines


Fixing bug introduced in the computeconf script, introduced in my zealous attempt to eliminate redundant preprocess() calls - but you still need one!


Also simplified the names of a few of the command-line arguments, and made the extensive-form solve the default.

........

r3031 | jwatson | 2010-09-20 17:35:24 -0500 (Mon, 20 Sep 2010) | 3 lines


Adding t-table entries to PySP computeconf script, to automate the computation of (1-alpha) confidence intervals.

........

r3032 | jwatson | 2010-09-20 17:37:42 -0500 (Mon, 20 Sep 2010) | 3 lines


Adding error-handling logic (pretty exception trapping/output) to PySP computeconf script.

........

r3033 | jwatson | 2010-09-20 17:50:14 -0500 (Mon, 20 Sep 2010) | 3 lines


Fixed an off-by-one error in the computation of the gap variance estimator.

........

r3034 | jwatson | 2010-09-20 22:28:20 -0500 (Mon, 20 Sep 2010) | 3 lines


Added computation of confidence interval width to the PySP computeconf script, leveraging the built-in t-table - if entries for the given degrees of freedom are in the table. Other output cleanup as well.

........

r3035 | jwatson | 2010-09-21 10:49:05 -0500 (Tue, 21 Sep 2010) | 3 lines


Added boundary issue checks in computeconf script, in cases where there aren't enough samples to populate n_g batches to compute a confidence interval.

........

r3036 | jwatson | 2010-09-21 14:56:22 -0500 (Tue, 21 Sep 2010) | 3 lines


Forgot to add gbar (the base offset) to the confidence interval width computation!

........

r3038 | jwatson | 2010-09-22 16:12:01 -0500 (Wed, 22 Sep 2010) | 3 lines


Fairly extensive re-work and simplification of the PySP extensive form writer, focusing on a cleaner implementation of the cvar term generation.

........

r3039 | jwatson | 2010-09-23 15:06:39 -0500 (Thu, 23 Sep 2010) | 3 lines


Modified CI computation script to deal with some weirdness in CPLEX abruptly terminating when running Sarah's GEP problem with CVAR as the optimization objective. There is a feasible solution, but no indication of why CPLEX stopped. Reporting bug to IBM.

........

Location:
coopr.pysp/stable
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • coopr.pysp/stable

  • coopr.pysp/stable/coopr/pysp/ef.py

    r3006 r3047  
    6969#       instances are associated with the extensive form. this might be something we
    7070#       encapsulate at some later time.
     71# NOTE: if cvar terms are generated, then the input scenario tree is modified accordingly,
     72#       i.e., with the addition of the "eta" variable at the root node and the excess
     73#       variables at (for lack of a better place - they are per-scenario, but are not
     74#       blended) the second stage.
    7175#
    7276
     
    8993         print "CVaR alpha="+str(risk_alpha)
    9094         print ""
     95
     96      # update the scenario tree with cvar-specific variable names, so
     97      # they will get cloned accordingly in the master instance.
     98      first_stage = scenario_tree._stages[0]
     99      second_stage = scenario_tree._stages[1]
     100      root_node = first_stage._tree_nodes[0]
     101
     102      # NOTE: because we currently don't have access to the reference
     103      #       instance in this method, temporarily (and largely orphaned)
     104      #       variables are constructed to supply to the scenario tree.
     105      #       this decision should be ultimately revisited.
     106      cvar_eta_variable_name = "CVAR_ETA"
     107      cvar_eta_variable = Var(name=cvar_eta_variable_name)
     108      cvar_eta_variable.construct()               
     109
     110      first_stage.add_variable(cvar_eta_variable, "*", [None])
     111
     112      cvar_excess_variable_name = "CVAR_EXCESS"
     113      cvar_excess_variable = Var(name=cvar_excess_variable_name)
     114      cvar_excess_variable.construct()
     115
     116      second_stage.add_variable(cvar_excess_variable, "*", [None])
     117
     118      # create the eta and excess variable on a per-scenario basis,
     119      # in addition to the constraint relating to the two.
     120      for scenario_name, scenario_instance in scenario_instances.items():
     121
     122         cvar_excess_variable_name = "CVAR_EXCESS"
     123         cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals)
     124         cvar_excess_variable.construct()         
     125         setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable)
     126
     127         cvar_eta_variable_name = "CVAR_ETA"
     128         cvar_eta_variable = Var(name=cvar_eta_variable_name)
     129         cvar_eta_variable.construct()         
     130         setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable)
     131
     132         compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS"
     133         compute_excess_constraint = Constraint(name=compute_excess_constraint_name)
     134         compute_excess_expression = cvar_excess_variable
     135         for node in scenario_tree._scenario_map[scenario_name]._node_list:
     136            (cost_variable, cost_variable_idx) = node._stage._cost_variable
     137            compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx]
     138         compute_excess_expression += cvar_eta_variable
     139         compute_excess_constraint.add(None, (0.0, compute_excess_expression, None))
     140         compute_excess_constraint._model = scenario_instance
     141         setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint)
     142
     143         # re-process the scenario instance due to the newly added constraints/variables associated
     144         # with CVaR. a complete preprocess is overkill, of course - the correct approach would be
     145         # to just preprocess those newly added variables and constraints.
     146         scenario_instance.preprocess()
    91147
    92148   #
     
    118174            if stage != scenario_tree._stages[-1]:     
    119175
    120                master_variable_name = tree_node._name + "_" + stage_variable.name
     176               master_variable_name = stage_variable.name               
    121177
    122178               # because there may be a single stage variable and multiple indices, check
     
    156212                  if (is_used is True) and (is_fixed is False):
    157213                           
    158                      # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
    159                      # and because presolve/simplification is name-based, the names *have* to be different.
    160                      master_variable[index].var = master_variable
    161                      master_variable[index].name = tree_node._name + "_" + master_variable[index].name
    162 
    163214                     for scenario in tree_node._scenarios:
    164215                        scenario_instance = scenario_instances[scenario._name]
     
    181232         setattr(binding_instance, expected_cost_variable_name, expected_cost_variable)
    182233
    183    # if we're generating the weighted CVaR objective term, create the
    184    # corresponding variable and the master CVaR eta variable.
    185234   if generate_weighted_cvar is True:
    186235
    187       root_node = scenario_tree._stages[0]._tree_nodes[0]
    188      
    189236      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
    190237      cvar_cost_variable = Var(name=cvar_cost_variable_name)
     238      cvar_cost_variable.construct()           
    191239      setattr(binding_instance, cvar_cost_variable_name, cvar_cost_variable)
    192       cvar_cost_variable.construct()
    193 
    194       cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
    195       cvar_eta_variable = Var(name=cvar_eta_variable_name)
    196       setattr(binding_instance, cvar_eta_variable_name, cvar_eta_variable)     
    197       cvar_eta_variable.construct()
    198240
    199241   binding_instance.preprocess()
     
    245287   # i.e., the current node cost and the expected cost of the child nodes.
    246288   # if the root, then the constraint is just the objective.
    247 
    248289   for stage in scenario_tree._stages:
    249290
     
    284325               cvar_cost_variable = getattr(binding_instance, cvar_cost_variable_name)
    285326               if cvar_weight == 0.0:
     327                  # if the cvar weight is 0, then we're only doing cvar - no mean.
    286328                  opt_expression = cvar_cost_variable                             
    287329               else:
     
    298340   if generate_weighted_cvar is True:
    299341     
    300       root_node = scenario_tree._stages[0]._tree_nodes[0]
    301 
    302       master_cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
    303       master_cvar_eta_variable = getattr(binding_instance, master_cvar_eta_variable_name)
    304      
    305       for scenario_name in scenario_instances.keys():
    306          scenario_instance = scenario_instances[scenario_name]
    307 
    308          # unique names are required because the presolve isn't
    309          # aware of the "owning" models for variables.
    310          cvar_excess_variable_name = "CVAR_EXCESS_"+scenario_name
    311          cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals)
    312          setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable)
    313          cvar_excess_variable.construct()
    314 
    315          cvar_eta_variable_name = "CVAR_ETA"
    316          cvar_eta_variable = Var(name=cvar_eta_variable_name)
    317          setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable)
    318          cvar_eta_variable.construct()
    319 
    320          compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS"
    321          compute_excess_constraint = Constraint(name=compute_excess_constraint_name)
    322          compute_excess_expression = cvar_excess_variable
    323          for node in scenario_tree._scenario_map[scenario_name]._node_list:
    324             (cost_variable, cost_variable_idx) = node._stage._cost_variable
    325             compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx]
    326          compute_excess_expression += cvar_eta_variable
    327          compute_excess_constraint.add(None, (0.0, compute_excess_expression, None))
    328          compute_excess_constraint._model = scenario_instance
    329          setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint)
    330 
    331          eta_equality_constraint_name = "MASTER_ETA_EQUALITY_WITH_" + scenario_instance.name
    332          eta_equality_constraint = Constraint(name=eta_equality_constraint_name)
    333          eta_equality_expr = master_cvar_eta_variable - cvar_eta_variable
    334          eta_equality_constraint.add(None, (0.0, eta_equality_expr, 0.0))
    335          eta_equality_constraint._model = binding_instance
    336          setattr(binding_instance, eta_equality_constraint_name, eta_equality_constraint)
    337 
    338          # re-process the scenario instance due to the newly added constraints/variables associated
    339          # with CVaR. a complete preprocess is overkill, of course - the correct approach would be
    340          # to just preprocess those newly added variables and constraints.
    341          scenario_instance.preprocess()
    342 
    343342      # add the constraint to compute the master CVaR variable value. iterate
    344343      # over scenario instances to create the expected excess component first.
    345344      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
    346345      cvar_cost_variable = getattr(binding_instance, cvar_cost_variable_name)
    347       cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
     346      cvar_eta_variable_name = "CVAR_ETA"
    348347      cvar_eta_variable = getattr(binding_instance, cvar_eta_variable_name)
    349348     
    350349      cvar_cost_expression = cvar_cost_variable - cvar_eta_variable
    351350     
    352       for scenario_name in scenario_instances.keys():
    353          scenario_instance = scenario_instances[scenario_name]
     351      for scenario_name, scenario_instance in scenario_instances.items():
     352
    354353         scenario_probability = scenario_tree._scenario_map[scenario_name]._probability
    355354
    356          scenario_excess_variable_name = "CVAR_EXCESS_"+scenario_name
     355         scenario_excess_variable_name = "CVAR_EXCESS"
    357356         scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name)
    358357
  • coopr.pysp/stable/coopr/pysp/phutils.py

    r3006 r3047  
    180180#
    181181# construct a scenario instance!
     182# IMPT: the returned instance is guaranteed to be preprocessed().
    182183#
    183184def construct_scenario_instance(scenario_tree_instance,
     
    217218         scenario_data.read(model=scenario_instance)
    218219         scenario_instance.load(scenario_data)
     220         scenario_instance.preprocess()         
    219221   except:
    220222      print "Encountered exception in model instance creation - traceback:"
  • coopr.pysp/stable/coopr/pysp/scenariotree.py

    r2852 r3047  
    2323
    2424   #
    25    # initialize the _solutions attribute of a tree node.
    26    #
    27 
    28    def _initialize_solutions(self):
     25   # update the _solutions attribute of a tree node, given a specific
     26   # variable/match-template/variable-index triple.
     27   #
     28   def _update_solution_map(self, variable, match_template, variable_indices):
     29
     30      # don't bother copying bounds for variables, as the values stored
     31      # here are computed elsewhere - and that entity is responsible for
     32      # ensuring feasibility. this also leaves room for specifying infeasible
     33      # or partial solutions.
     34      new_variable_index = variable._index
     35      new_variable_name = variable.name
     36      new_variable = None
     37      if (len(new_variable_index) is 1) and (None in new_variable_index):
     38         new_variable = Var(name=new_variable_name)
     39      else:
     40         new_variable = Var(new_variable_index, name=new_variable_name)
     41      new_variable.construct()
     42
     43      # by default, deactive all variable values - we're copying the
     44      # full index set of a variable, which is necessarily wasteful and
     45      # incorrect. then, do a second pass and activate those indicies
     46      # that are actually associated with this tree node.
     47      for index in new_variable_index:
     48         new_variable[index].deactivate()
     49
     50      for index in variable_indices:
     51         new_variable[index].activate()
     52
     53      self._solutions[new_variable_name] = new_variable     
     54
     55   #
     56   # initialize the _solutions attribute of a tree node, from scratch.
     57   #
     58
     59   def _initialize_solution_map(self):
     60
     61      # clear whatever was there before.
     62      self._solutions = {}
    2963
    3064      # NOTE: Given the expense, this should really be optional - don't
     
    3670      # is an issue (space is the most likely issue)
    3771      for variable, match_template, variable_indices in self._stage._variables:
    38 
    39          # don't bother copying bounds for variables, as the values stored
    40          # here are computed elsewhere - and that entity is responsible for
    41          # ensuring feasibility. this also leaves room for specifying infeasible
    42          # or partial solutions.
    43          new_variable_index = variable._index
    44          new_variable_name = variable.name
    45          new_variable = None
    46          if (len(new_variable_index) is 1) and (None in new_variable_index):
    47             new_variable = Var(name=new_variable_name)
    48          else:
    49             new_variable = Var(new_variable_index, name=new_variable_name)
    50          new_variable.construct()
    51 
    52          # by default, deactive all variable values - we're copying the
    53          # full index set of a variable, which is necessarily wasteful and
    54          # incorrect. then, do a second pass and activate those indicies
    55          # that are actually associated with this tree node.
    56          for index in new_variable_index:
    57             new_variable[index].deactivate()
    58 
    59          for index in variable_indices:
    60             new_variable[index].activate()
    61 
    62          self._solutions[new_variable_name] = new_variable
    63 
    64       self._solutions_initialized = True
     72         self._update_solution_map(variable, match_template, variable_indices)
     73
     74      self._solution_map_initialized = True
    6575
    6676
     
    92102
    93103      # a flag indicating whether the _solutions attribute has been properly initialized.
    94       self._solutions_initialized = False
     104      self._solution_map_initialized = False
    95105
    96106      # solution (variable) values for this node. assumed to be distinct
     
    101111
    102112      if initialize_solution is True:
    103          self._initialize_solutions()
     113         self._initialize_solution_map()
    104114
    105115   #
     
    110120   def snapshotSolutionFromAverages(self):
    111121
    112       if self._solutions_initialized is False:
    113          self._initialize_solutions()
     122      if self._solution_map_initialized is False:
     123         self._initialize_solution_map()
    114124
    115125      for variable_name, variable in self._solutions.items():
     
    134144   def snapshotSolutionFromInstances(self, scenario_instance_map):
    135145
    136       if self._solutions_initialized is False:
    137          self._initialize_solutions()
     146      if self._solution_map_initialized is False:
     147         self._initialize_solution_map()
    138148
    139149      for variable_name, variable in self._solutions.items():
     
    174184   """
    175185   def __init__(self, *args, **kwds):
     186
    176187      self._name = ""
    177       self._tree_nodes = []      # a collection of ScenarioTreeNodes
    178       # a collection of pairs consisting of (1) references to pyomo model Vars, (2) the original match template string (for output purposes),
    179       # and (3) a *list* of the corresponding indices.
    180       # the variables are references to those objects belonging to the instance in the parent ScenarioTree.
     188     
     189      # a collection of ScenarioTreeNode objects.
     190      self._tree_nodes = []     
     191     
     192      # a collection of triples consisting of (1) a reference to a Pyomo model Var object, (2) the original match
     193      # template string (for output purposes), and (3) a *list* of the corresponding indices. the variables are
     194      # references to those objects belonging to the Pyomo reference scenario instance associated with the parent
     195      # ScenarioTree of this Stage.
    181196      # NOTE: if the variable index is none, it is assumed that the entire variable is blended.
    182197      self._variables = []
    183       # a tuple consisting of (1) a reference to a pyomo model Var that computes the stage-specific cost and (2) the corresponding index.
    184       # the index *is* the sole index in the cost variable, as the cost variable refers to a single variable index.
     198     
     199      # a tuple consisting of (1) a reference to a pyomo model Var that computes the stage-specific cost and (2) the corresponding
     200      # index. the index *is* the sole index in the cost variable, as the cost variable refers to a single variable index.
    185201      self._cost_variable = (None, None)
     202
     203   #
     204   # add a new variable to the stage, which will include updating the solution maps for each associated ScenarioTreeNode.
     205   #
     206   def add_variable(self, variable, match_template, indices):
     207
     208      self._variables.append((variable, match_template, indices))
     209
     210      for tree_node in self._tree_nodes:
     211         tree_node._update_solution_map(variable, match_template, indices)
    186212
    187213class Scenario(object):
     
    190216   """
    191217   def __init__(self, *args, **kwds):
     218     
    192219      self._name = None
    193220      self._leaf_node = None  # allows for construction of node list
     
    328355   """
    329356   def __init__(self, *args, **kwds):
    330       self._name = None # TBD - some arbitrary identifier
    331       self._reference_instance = None # TBD - the reference (deterministic) base model
     357     
     358      self._name = None # some arbitrary identifier
     359      self._reference_instance = None # the reference (deterministic) base model
    332360
    333361      # the core objects defining the scenario tree.
     
    341369      self._scenario_map = {}
    342370
    343       # mapping of stages to sets of variables which belong in the corresponding stage.
    344       self._stage_variables = {}
    345 
    346371      # a boolean indicating how data for scenario instances is specified.
    347372      # possibly belongs elsewhere, e.g., in the PH algorithm.
    348373      self._scenario_based_data = None
    349 
    350       # every stage has a cost variable - this is a variable/index pair.
    351       self._cost_variable = None
    352374
    353375      scenario_tree_instance = None
  • coopr.pysp/stable/scripts/computeconf

    r3006 r3047  
    2020from coopr.pysp.ef import *
    2121
     22# this is a hack, in order to pick up the UndefinedData class. this is needed currently, as
     23# CPLEX is periodically barfing on cvar formulations, yielding an undefined gap. technically,
     24# the gap is defined and the solution is feasible, but a correct fix to the CPLEX plugin
     25# would yield a complete failure to solve cvar problems. see related hacks below, searching
     26# for CVARHACK.
     27from coopr.opt.results.container import *
     28
     29# to avoid the pain of user lookup of parameter in t-tables, we provide decent coverage automatically.
     30# feel free to add more values!!!! maps degrees-of-freedom to (alpha,t-statistic) pairs.
     31
     32t_table_values = {
     33
     341 :  {0.25  : 1.000  ,0.1   : 3.078  ,0.05  : 6.314  ,0.025 : 12.706 ,0.01  : 31.821 ,0.005 : 63.657 ,0.001 : 318.309},
     352 :  {0.25  : 0.816  ,0.1   : 1.886  ,0.05  : 2.920  ,0.025 : 4.303  ,0.01  : 6.965  ,0.005 : 9.925  ,0.001 : 22.327},
     363 :  {0.25  : 0.765  ,0.1   : 1.638  ,0.05  : 2.353  ,0.025 : 3.182  ,0.01  : 4.541  ,0.005 : 5.841  ,0.001 : 10.215},
     374 :  {0.25  : 0.741 , 0.1   : 1.533 , 0.05  : 2.132 , 0.025 : 2.776 , 0.01  : 3.747 , 0.005 : 4.604 , 0.001 : 7.173},
     385 :  {0.25  : 0.727 , 0.1   : 1.476 , 0.05  : 2.015 , 0.025 : 2.571 , 0.01  : 3.365 , 0.005 : 4.032 , 0.001 : 5.893},
     396 :  {0.25  : 0.718 , 0.1   : 1.440 , 0.05  : 1.943 , 0.025 : 2.447 , 0.01  : 3.143 , 0.005 : 3.707 , 0.001 : 5.208},
     407 :  {0.25  : 0.711 , 0.1   : 1.415 , 0.05  : 1.895 , 0.025 : 2.365 , 0.01  : 2.998 , 0.005 : 3.499 , 0.001 : 4.785},
     418 :  {0.25  : 0.706 , 0.1   : 1.397 , 0.05  : 1.860 , 0.025 : 2.306 , 0.01  : 2.896 , 0.005 : 3.355 , 0.001 : 4.501},
     429 :  {0.25  : 0.703 , 0.1   : 1.383 , 0.05  : 1.833 , 0.025 : 2.262 , 0.01  : 2.821 , 0.005 : 3.250 , 0.001 : 4.297},
     4310 : {0.25  : 0.700 , 0.1   : 1.372 , 0.05  : 1.812 , 0.025 : 2.228 , 0.01  : 2.764 , 0.005 : 3.169 , 0.001 : 4.144},
     4411 : {0.25  : 0.697 , 0.1   : 1.363 , 0.05  : 1.796 , 0.025 : 2.201 , 0.01  : 2.718 , 0.005 : 3.106 , 0.001 : 4.025},
     4512 : {0.25  : 0.695 , 0.1   : 1.356 , 0.05  : 1.782 , 0.025 : 2.179 , 0.01  : 2.681 , 0.005 : 3.055 , 0.001 : 3.930},
     4613 : {0.25  : 0.694 , 0.1   : 1.350 , 0.05  : 1.771 , 0.025 : 2.160 , 0.01  : 2.650 , 0.005 : 3.012 , 0.001 : 3.852},
     4714 : {0.25  : 0.692 , 0.1   : 1.345 , 0.05  : 1.761 , 0.025 : 2.145 , 0.01  : 2.624 , 0.005 : 2.977 , 0.001 : 3.787},
     4815 : {0.25  : 0.691 , 0.1   : 1.341 , 0.05  : 1.753 , 0.025 : 2.131 , 0.01  : 2.602 , 0.005 : 2.947 , 0.001 : 3.733},
     4916 : {0.25  : 0.690 , 0.1   : 1.337 , 0.05  : 1.746 , 0.025 : 2.120 , 0.01  : 2.583 , 0.005 : 2.921 , 0.001 : 3.686},
     5017 : {0.25  : 0.689 , 0.1   : 1.333 , 0.05  : 1.740 , 0.025 : 2.110 , 0.01  : 2.567 , 0.005 : 2.898 , 0.001 : 3.646},
     5118 : {0.25  : 0.688 , 0.1   : 1.330 , 0.05  : 1.734 , 0.025 : 2.101 , 0.01  : 2.552 , 0.005 : 2.878 , 0.001 : 3.610},
     5219 : {0.25  : 0.688 , 0.1   : 1.328 , 0.05  : 1.729 , 0.025 : 2.093 , 0.01  : 2.539 , 0.005 : 2.861 , 0.001 : 3.579},
     5320 : {0.25  : 0.687 , 0.1   : 1.325 , 0.05  : 1.725 , 0.025 : 2.086 , 0.01  : 2.528 , 0.005 : 2.845 , 0.001 : 3.552},
     5421 : {0.25  : 0.686 , 0.1   : 1.323 , 0.05  : 1.721 , 0.025 : 2.080 , 0.01  : 2.518 , 0.005 : 2.831 , 0.001 : 3.527},
     5522 : {0.25  : 0.686 , 0.1   : 1.321 , 0.05  : 1.717 , 0.025 : 2.074 , 0.01  : 2.508 , 0.005 : 2.819 , 0.001 : 3.505},
     5623 : {0.25  : 0.685 , 0.1   : 1.319 , 0.05  : 1.714 , 0.025 : 2.069 , 0.01  : 2.500 , 0.005 : 2.807 , 0.001 : 3.485},
     5724 : {0.25  : 0.685 , 0.1   : 1.318 , 0.05  : 1.711 , 0.025 : 2.064 , 0.01  : 2.492 , 0.005 : 2.797 , 0.001 : 3.467},
     5825 : {0.25  : 0.684 , 0.1   : 1.316 , 0.05  : 1.708 , 0.025 : 2.060 , 0.01  : 2.485 , 0.005 : 2.787 , 0.001 : 3.450},
     5926 : {0.25  : 0.684 , 0.1   : 1.315 , 0.05  : 1.706 , 0.025 : 2.056 , 0.01  : 2.479 , 0.005 : 2.779 , 0.001 : 3.435},
     6027 : {0.25  : 0.684 , 0.1   : 1.314 , 0.05  : 1.703 , 0.025 : 2.052 , 0.01  : 2.473 , 0.005 : 2.771 , 0.001 : 3.421},
     6128 : {0.25  : 0.683 , 0.1   : 1.313 , 0.05  : 1.701 , 0.025 : 2.048 , 0.01  : 2.467 , 0.005 : 2.763 , 0.001 : 3.408},
     6229 : {0.25  : 0.683 , 0.1   : 1.311 , 0.05  : 1.699 , 0.025 : 2.045 , 0.01  : 2.462 , 0.005 : 2.756 , 0.001 : 3.396},
     6330 : {0.25  : 0.683 , 0.1   : 1.310 , 0.05  : 1.697 , 0.025 : 2.042 , 0.01  : 2.457 , 0.005 : 2.750 , 0.001 : 3.385},
     6431 : {0.25  : 0.682 , 0.1   : 1.309 , 0.05  : 1.696 , 0.025 : 2.040 , 0.01  : 2.453 , 0.005 : 2.744 , 0.001 : 3.375},
     6532 : {0.25  : 0.682 , 0.1   : 1.309 , 0.05  : 1.694 , 0.025 : 2.037 , 0.01  : 2.449 , 0.005 : 2.738 , 0.001 : 3.365},
     6633 : {0.25  : 0.682 , 0.1   : 1.308 , 0.05  : 1.692 , 0.025 : 2.035 , 0.01  : 2.445 , 0.005 : 2.733 , 0.001 : 3.356},
     6734 : {0.25  : 0.682 , 0.1   : 1.307 , 0.05  : 1.691 , 0.025 : 2.032 , 0.01  : 2.441 , 0.005 : 2.728 , 0.001 : 3.348},
     6835 : {0.25  : 0.682 , 0.1   : 1.306 , 0.05  : 1.690 , 0.025 : 2.030 , 0.01  : 2.438 , 0.005 : 2.724 , 0.001 : 3.340},
     6936 : {0.25  : 0.681 , 0.1   : 1.306 , 0.05  : 1.688 , 0.025 : 2.028 , 0.01  : 2.434 , 0.005 : 2.719 , 0.001 : 3.333},
     7037 : {0.25  : 0.681 , 0.1   : 1.305 , 0.05  : 1.687 , 0.025 : 2.026 , 0.01  : 2.431 , 0.005 : 2.715 , 0.001 : 3.326},
     7138 : {0.25  : 0.681 , 0.1   : 1.304 , 0.05  : 1.686 , 0.025 : 2.024 , 0.01  : 2.429 , 0.005 : 2.712 , 0.001 : 3.319},
     7239 : {0.25  : 0.681 , 0.1   : 1.304 , 0.05  : 1.685 , 0.025 : 2.023 , 0.01  : 2.426 , 0.005 : 2.708 , 0.001 : 3.313},
     7340 : {0.25  : 0.681 , 0.1   : 1.303 , 0.05  : 1.684 , 0.025 : 2.021 , 0.01  : 2.423 , 0.005 : 2.704 , 0.001 : 3.307},
     7441 : {0.25  : 0.681 , 0.1   : 1.303 , 0.05  : 1.683 , 0.025 : 2.020 , 0.01  : 2.421 , 0.005 : 2.701 , 0.001 : 3.301},
     7542 : {0.25  : 0.680 , 0.1   : 1.302 , 0.05  : 1.682 , 0.025 : 2.018 , 0.01  : 2.418 , 0.005 : 2.698 , 0.001 : 3.296},
     7643 : {0.25  : 0.680 , 0.1   : 1.302 , 0.05  : 1.681 , 0.025 : 2.017 , 0.01  : 2.416 , 0.005 : 2.695 , 0.001 : 3.291},
     7744 : {0.25  : 0.680 , 0.1   : 1.301 , 0.05  : 1.680 , 0.025 : 2.015 , 0.01  : 2.414 , 0.005 : 2.692 , 0.001 : 3.286},
     7845 : {0.25  : 0.680 , 0.1   : 1.301 , 0.05  : 1.679 , 0.025 : 2.014 , 0.01  : 2.412 , 0.005 : 2.690 , 0.001 : 3.281},
     7946 : {0.25  : 0.680 , 0.1   : 1.300 , 0.05  : 1.679 , 0.025 : 2.013 , 0.01  : 2.410 , 0.005 : 2.687 , 0.001 : 3.277},
     8047 : {0.25  : 0.680 , 0.1   : 1.300 , 0.05  : 1.678 , 0.025 : 2.012 , 0.01  : 2.408 , 0.005 : 2.685 , 0.001 : 3.273},
     8148 : {0.25  : 0.680 , 0.1   : 1.299 , 0.05  : 1.677 , 0.025 : 2.011 , 0.01  : 2.407 , 0.005 : 2.682 , 0.001 : 3.269},
     8249 : {0.25  : 0.680 , 0.1   : 1.299 , 0.05  : 1.677 , 0.025 : 2.010 , 0.01  : 2.405 , 0.005 : 2.680 , 0.001 : 3.265},
     8350 : {0.25  : 0.679 , 0.1   : 1.299 , 0.05  : 1.676 , 0.025 : 2.009 , 0.01  : 2.403 , 0.005 : 2.678 , 0.001 : 3.261}
     84
     85}
     86
    2287def run(args=None):
    2388
    2489   try:
    2590      conf_options_parser = construct_ph_options_parser("computeconf [options]")
    26       conf_options_parser.add_option("--Fraction-of-Scenarios-for-Solve",
     91      conf_options_parser.add_option("--fraction-scenarios-for-solve",
    2792                                     help="The fraction of scenarios that are allocated to finding a solution. Default is 0.5",
    2893                                     action="store",
     
    3095                                     type="float",
    3196                                     default=0.5)
    32       conf_options_parser.add_option("--Number-of-Samples-for-Confidence-Interval",
     97      conf_options_parser.add_option("--number-samples-for-confidence-interval",
    3398                                     help="The number of samples of scenarios that are allocated to the confidence inteval (n_g). Default is 10",
    3499                                     action="store",
     
    36101                                     type="int",
    37102                                     default=10)
    38       conf_options_parser.add_option("--Confidence-Interval-Alpha",
    39                                      help="The alpha level for the confidence interval. Default is 0.10",
    40                                      action="store",
    41                                      dest="Conf_Alpha",
     103      conf_options_parser.add_option("--confidence-interval-alpha",
     104                                     help="The alpha level for the confidence interval. Default is 0.05",
     105                                     action="store",
     106                                     dest="confidence_interval_alpha",
    42107                                     type="float",
    43                                      default=0.10)
    44       conf_options_parser.add_option("--solve-xhat-with-ef-only",
    45                                      help="Perform xhat solve via EF rather than PH. Default is False",
     108                                     default=0.05)
     109      conf_options_parser.add_option("--solve-xhat-with-ph",
     110                                     help="Perform xhat solve via PH rather than an EF solve. Default is False",
    46111                                     action="store_true",
    47                                      dest="solve_xhat_with_ef_only",
     112                                     dest="solve_xhat_with_ph",
    48113                                     default=False)
    49114      conf_options_parser.add_option("--random-seed",
     
    97162   reference_model, reference_instance, full_scenario_tree, scenario_tree_instance = load_reference_and_scenario_models(options)
    98163   if (reference_model is None) or (reference_instance is None) or (full_scenario_tree is None) or (scenario_tree_instance is None):
    99       print "***ERROR: Failed to initialize reference model/instance pair and/or the scenario tree."
    100       sys.exit(1)
     164      raise RuntimeError, "***ERROR: Failed to initialize reference model/instance pair and/or the scenario tree."
    101165
    102166   # verify that we're dealing with a minimization problem - the script might be easily made to work with
     
    105169   objective = reference_model.active_components()[Objective][objective_name]
    106170   if objective.sense != minimize:
    107       print "***ERROR: Confidence intervals are available only for minimization problems"
    108       sys.exit(1)
     171      raise RuntimeError, "***ERROR: Confidence intervals are available only for minimization problems"
    109172
    110173   print "Starting confidence interval calculation..."
     
    112175   scenario_count = len(full_scenario_tree._stages[-1]._tree_nodes)
    113176   if len(full_scenario_tree._stages) > 2:
    114       print "***ERROR: Confidence intervals are available only for two stage stochastic programs;"+str(len(full_scenario_tree._stages))+" stages specified"
    115       sys.exit(1)
     177      raise RuntimeError, "***ERROR: Confidence intervals are available only for two stage stochastic programs;"+str(len(full_scenario_tree._stages))+" stages specified"
    116178
    117179   # randomly permute the indices to extract a subset to compute xhat.
     
    126188   num_scenarios_per_sample = int((scenario_count - num_scenarios_for_solution) / n_g) # 'n' in Morton's slides
    127189   wasted_scenarios = scenario_count - num_scenarios_for_solution - n_g * num_scenarios_per_sample
     190
     191   if num_scenarios_per_sample is 0:
     192      raise RuntimeError, "Computed number of scenarios per sample group equals 0 - "+str(scenario_count - num_scenarios_for_solution)+" scenarios cannot be divided into "+str(n_g)+" groups!"
    128193   
    129194   print "Problem contains "+str(scenario_count)+" scenarios, of which "+str(num_scenarios_for_solution)+" will be used to find a solution xhat."
     
    135200   # directly, mainly out of convenience - we're leveraging the code in ph_for_bundle to create the
    136201   # scenario tree and scenario instances.
     202   print ""
    137203   print "Loading scenario instances and initializing scenario tree for xhat scenario bundle."
    138204   xhat_ph = ph_for_bundle(0, num_scenarios_for_solution, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options)
    139205   xhat_obj = None
    140206
    141    if options.solve_xhat_with_ef_only is True:
     207   if options.solve_xhat_with_ph is False:
    142208      print "Creating the xhat extensive form."
    143209      if options.verbose is True:
     
    185251      start_index = num_scenarios_for_solution + (k-1)*num_scenarios_per_sample
    186252      stop_index = start_index + num_scenarios_per_sample - 1
     253      print ""
    187254      print "Computing statistics for sample k="+str(k)+"."
    188255      if options.verbose is True:
     
    210277      xstar_obj = float(ef_results.solution(0).objective['f'].value)  ## DLW to JPW: we need the gap too, and to add/subtract is as necessary.
    211278      xstar_obj_gap = ef_results.solution(0).gap # assuming this is the absolute gap
    212       xstar_obj_bound = xstar_obj - xstar_obj_gap # HERE - watch for signs on gap - is it absolute, or can it be negative?     
    213       print "Extensive form objective value given xstar="+str(xstar_obj)+"; gap="+str(xstar_obj_gap)
     279      # CVARHACK: is CPLEX barfed, keep trucking and bury our head in the sand.
     280      if type(xstar_obj_gap) is UndefinedData:
     281         xstar_obj_bound = xstar_obj
     282      else:
     283         xstar_obj_bound = xstar_obj - xstar_obj_gap
     284      # TBD: ADD VERBOSE OUTPUT HERE
    214285
    215286      # to get f(xhat) for this sample, fix the first-stage variables and re-solve the extensive form.
     
    236307   # second pass for variance calculation (because we like storing the g_supk)
    237308   g_var = 0
    238    for k in range(0, n_g-1):
     309   for k in range(0, n_g):
    239310      g_var = g_var + (g_supk_of_xhat[k] - g_bar) * (g_supk_of_xhat[k] - g_bar)
    240311   g_var = g_var / (n_g - 1)    # sample var
     312   print ""
     313   print "Raw results:"
    241314   print "g_bar=",g_bar
    242315   print "g_stddev=",math.sqrt(g_var)
    243316   print "Average f(xhat)=",sum_xstar_obj_given_xhat / n_g
    244317
     318   if n_g in t_table_values:
     319      print ""
     320      print "Results summary:"
     321      t_table_entries = t_table_values[n_g]
     322      for key in sorted(t_table_entries.keys()):
     323         print "Confidence interval width for alpha="+str(key)+" is "+str(g_bar + (t_table_entries[key] * math.sqrt(g_var) / math.sqrt(n_g)))
     324   else:
     325      print "No built-in t-table entries for "+str(n_g)+" degrees of freedom - cannot calculate confidence interval width"
     326
     327   if options.write_xhat_solution is True:
     328      print ""           
     329      print "xhat solution:"
     330      scenario_tree = xhat_ph._scenario_tree
     331      first_stage = scenario_tree._stages[0]
     332      root_node = first_stage._tree_nodes[0]
     333      for key, value in root_node._solutions.items():
     334         for idx in value:
     335            if value[idx]() != 0.0:
     336               print key, idx, value[idx]()     
     337
    245338   if options.append_file is not None:
    246339      output_file = open(options.append_file, "a")
    247340      print >>output_file, "\ninstancedirectory, ",options.instance_directory, ", seed, ",options.random_seed, ", N, ",scenario_count, ", hatn, ",num_scenarios_for_solution, ", n_g, ", options.n_g, ", Eoffofxhat, ", sum_xstar_obj_given_xhat / n_g, ", gbar, ", g_bar, ", sg, ", math.sqrt(g_var), ", objforxhat, ", xhat_obj, ", n,",num_scenarios_per_sample,
     341
     342      if n_g in t_table_values:
     343         t_table_entries = t_table_values[n_g]
     344         for key in sorted(t_table_entries.keys()):
     345            print >>output_file, " , alpha="+str(key)+" , "+str(g_bar + (t_table_entries[key] * math.sqrt(g_var) / math.sqrt(n_g))),
    248346
    249347      if options.write_xhat_solution is True:
     
    257355                  print >>output_file, key, idx, value[idx](),
    258356      output_file.close()
     357      print ""                 
     358      print "Results summary appended to file="+options.append_file
    259359
    260360#
     
    274374                                         scenariobundlelist=scenarios_to_bundle)
    275375   if scenario_tree_for_soln.validate() is False:
    276       print "***ERROR: Bundled scenario tree is invalid!!!"
    277       sys.exit(1)
     376      raise RuntimeError, "***ERROR: Bundled scenario tree is invalid!!!"
    278377
    279378   ph = create_ph_from_scratch(options, reference_model, reference_instance, scenario_tree_for_soln)
     
    325424      else:
    326425         ef_solver.mipgap = options.ef_mipgap
     426   if options.keep_solver_files is True:
     427      ef_solver.keepFiles = True                   
    327428
    328429   ef_solver_manager = SolverManagerFactory(options.solver_manager_type)
     
    342443#
    343444
    344 run()
     445try:
     446   run()
     447except ValueError, str:
     448    print "VALUE ERROR:"
     449    print str
     450except IOError, str:
     451    print "IO ERROR:"
     452    print str
     453except pyutilib.common.ApplicationError, str:
     454    print "APPLICATION ERROR:"
     455    print str
     456except RuntimeError, str:
     457    print "RUN-TIME ERROR:"   
     458    print str       
     459except:
     460   print "Encountered unhandled exception"
     461   traceback.print_exc()   
Note: See TracChangeset for help on using the changeset viewer.