Changeset 2917


Ignore:
Timestamp:
Aug 10, 2010 7:04:15 PM (9 years ago)
Author:
jwatson
Message:

Cleaning up the PySP confidence interval computation code, adding error checking, comments, etc.

Location:
coopr.pysp/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • coopr.pysp/trunk/coopr/pysp/phinit.py

    r2772 r2917  
    296296#
    297297# Create the reference model / instance and scenario tree instance for PH.
     298# IMPT: This method should be moved into a more generic module - it has nothing
     299#       to do with PH, and is used elsewhere (by routines that shouldn't have
     300#       to know about PH).
    298301#
    299302
  • coopr.pysp/trunk/scripts/computeconf

    r2838 r2917  
    2121def run(args=None):
    2222
    23 # JPW: args=None is OK, as the arguments are by default propagated through sys.argv
    24 #   if args == None:
    25 #      print "Error: testconf run called with no args."
    26 #      sys.exit(1)
    2723   try:
    28       conf_options_parser = construct_ph_options_parser("testconf [options]")
     24      conf_options_parser = construct_ph_options_parser("computeconf [options]")
    2925      conf_options_parser.add_option("--Fraction-of-Scenarios-for-Solve",
    3026                                     help="The fraction of scenarios that are allocated to finding a solution. Default is 0.5",
     
    4541                                     type="float",
    4642                                     default=0.10)
    47       conf_options_parser.add_option("--solve-hatx-with-ef-only",
    48                                      help="Perform hatx solve via EF rather than PH. Default is False",
     43      conf_options_parser.add_option("--solve-xhat-with-ef-only",
     44                                     help="Perform xhat solve via EF rather than PH. Default is False",
    4945                                     action="store_true",
    50                                      dest="solve_hatx_with_ef_only",
     46                                     dest="solve_xhat_with_ef_only",
    5147                                     default=False)
    5248      conf_options_parser.add_option("--random-seed",
     
    6965      return
    7066
     67   # seed the generator is a user-supplied seed is provided. otherwise,
     68   # python will seed from the current system time.
    7169   if options.random_seed > 0:
    7270      random.seed(options.random_seed)
    7371
    74    reference_model, reference_instance, full_scenario_tree, si = load_reference_and_scenario_models(options)
    75 
    76    # HERE: If we get an objective function with maximization, bail right away.
    77 
    78    print "Starting confidence interval calculations"
    79 #   print "full_scenario_tree pprint:"
    80 #   full_scenario_tree.pprint()
    81    scenariocount = len(full_scenario_tree._stages[-1]._tree_nodes)
     72   # create the reference instances and the scenario tree - no scenario instances yet.
     73   print "Loading reference model and scenario tree"
     74   reference_model, reference_instance, full_scenario_tree, scenario_tree_instance = load_reference_and_scenario_models(options)
     75   if (reference_model is None) or (reference_instance is None) or (full_scenario_tree is None) or (scenario_tree_instance is None):
     76      print "***ERROR: Failed to initialize reference model/instance pair and/or the scenario tree."
     77      sys.exit(1)
     78
     79   # verify that we're dealing with a minimization problem - the script might be easily made to work with
     80   # maximization problems - we just haven't done the work to do that.
     81   objective_name = reference_model.active_components()[Objective].keys()[0]
     82   objective = reference_model.active_components()[Objective][objective_name]
     83   if objective.sense != minimize:
     84      print "***ERROR: Confidence intervals are available only for minimization problems"
     85      sys.exit(1)
     86
     87   print "Starting confidence interval calculation..."
     88   
     89   scenario_count = len(full_scenario_tree._stages[-1]._tree_nodes)
    8290   if len(full_scenario_tree._stages) > 2:
    83       print "Confidence intervals are available only for two stage problems. Stage count=",len(full_scenario_tree._stages)
     91      print "***ERROR: Confidence intervals are available only for two stage stochastic programs;"+str(len(full_scenario_tree._stages))+" stages specified"
    8492      sys.exit(1)
    8593
    86    # random permutation of the indexes
    87    IndList = range(scenariocount)
    88    random.shuffle(IndList)
    89 #   print "After shuffle, IndList=",IndList
    90    
    91    num_scenarios_for_solution = int(options.fraction_for_solve * scenariocount)
     94   # randomly permute the indices to extract a subset to compute xhat.
     95   index_list = range(scenario_count)
     96   random.shuffle(index_list)
     97   if options.verbose is True:
     98      print "Random permutation of the scenario indices="+str(index_list)
     99   
     100   # figure out the scenario counts for both the xhat and confidence interval computations.
     101   num_scenarios_for_solution = int(options.fraction_for_solve * scenario_count)
    92102   n_g = options.n_g
    93    num_scenarios_per_sample = int((scenariocount - num_scenarios_for_solution) / n_g)  #n in Morton's slides
    94    wastedscenarios = scenariocount - num_scenarios_for_solution - n_g * num_scenarios_per_sample
    95    
    96    print str(scenariocount)+" scenarios, from which "+str(num_scenarios_for_solution)+" will be used to find a solution."
    97    print str(n_g)+" groups of "+str(num_scenarios_per_sample)+" will be used to compute the confidence interval."
    98    if wastedscenarios > 0:
    99       print "(so ",wastedscenarios," will not be used.)"
    100 
    101    # create ph for finding the solution
    102    hatxph = ph_for_bundle(0, num_scenarios_for_solution, reference_model, full_scenario_tree, reference_instance, si, IndList, options)
    103 
     103   num_scenarios_per_sample = int((scenario_count - num_scenarios_for_solution) / n_g) # 'n' in Morton's slides
     104   wasted_scenarios = scenario_count - num_scenarios_for_solution - n_g * num_scenarios_per_sample
     105   
     106   print "Problem contains "+str(scenario_count)+" scenarios, of which "+str(num_scenarios_for_solution)+" will be used to find a solution xhat."
     107   print "A total of "+str(n_g)+" groups of "+str(num_scenarios_per_sample)+" scenarios will be used to compute the confidence interval on xhat."
     108   if wasted_scenarios > 0:
     109      print "A total of "+str(wasted_scenarios)+" scenarios will not be used."
     110
     111   # create a ph object for finding the solution. we do this even if we're solving the extensive form
     112   # directly, mainly out of convenience - we're leveraging the code in ph_for_bundle to create the
     113   # scenario tree and scenario instances.
     114   xhatph = ph_for_bundle(0, num_scenarios_for_solution, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options)
    104115   xhatobj = None
    105116
    106    if options.solve_hatx_with_ef_only is True:
    107       # DLW: LOOK HERE!
    108       hatex_ef = create_ef_instance(hatxph._scenario_tree, hatxph._instances)
    109       ef_results = write_and_solve_ef(hatex_ef, hatxph._instances, options)
    110       load_ef_solution(ef_results, hatex_ef, hatxph._instances)
     117   if options.solve_xhat_with_ef_only is True:
     118
     119      hatex_ef = create_ef_instance(xhatph._scenario_tree, xhatph._instances)
     120      ef_results = write_and_solve_ef(hatex_ef, xhatph._instances, options)
     121      load_ef_solution(ef_results, hatex_ef, xhatph._instances)
    111122      # the following is just for printing things in an organized manner.
    112       hatxph._scenario_tree.snapshotSolutionFromInstances(hatxph._instances)
    113 #      hatxph._scenario_tree.pprintSolution()
    114 #      hatxph._scenario_tree.pprintCosts(hatxph._instances)
     123      xhatph._scenario_tree.snapshotSolutionFromInstances(xhatph._instances)
    115124      xhatobj = float(ef_results.solution(0).objective['f'].value)  ## DLW to JPW: we need the gap too
    116125   else:
    117       print "SOLVING HATX VIA PH!"
    118       hatxph.solve()
    119    print "so we now have the hat{x} variables in the PHAVG structure of hatxph"
     126      print "SOLVING XHAT VIA PH!"
     127      xhatph.solve()
     128      # TBD - grab xhatobj; complicated by the fact that PH may not have converged.
     129   print "Now have the hat{x} variables in the PHAVG structure of xhatph"
    120130
    121131   print "now form and solve the problems for each sample"
    122132   # in order to handle the case of scenarios that are not equally likely, we will split the expectations for Gsupk
    123133   # BUT we are going to assume that the groups themselves are equally likely and just scale by n_g and n_g-1 for Gbar and VarG
    124    G_supk_of_hatx = [] # really not always needed... http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
     134   G_supk_of_xhat = [] # really not always needed... http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
    125135   Gbar = 0
    126    sum_E_f_of_Hatx = 0
     136   sum_E_f_of_Xhat = 0
    127137   for k in range(1, n_g+1):
    128138      start = num_scenarios_for_solution + (k-1)*num_scenarios_per_sample
     
    132142      # NOTE: We'll never run this ph - code could be refactored
    133143
    134       # this computes hatx for the EF associated with sample k
    135       phforGk = ph_for_bundle(start, stop, reference_model, full_scenario_tree, reference_instance, si, IndList, options)
     144      # this computes xhat for the EF associated with sample k
     145      phforGk = ph_for_bundle(start, stop, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options)
    136146      efforGk = create_ef_instance(phforGk._scenario_tree, phforGk._instances)     
    137147      ef_results = write_and_solve_ef(efforGk, phforGk._instances, options)
     
    151161      E_f_of_xstar_bound = E_f_of_xstar - efforgK_gap # HERE - watch for signs on gap - is it absolute, or can it be negative?
    152162
    153       # to get fj(hatx) we need the obj value for each scenario, j, in the bundle evaluated at hat x
     163      # to get fj(xhat) we need the obj value for each scenario, j, in the bundle evaluated at hat x
    154164      # this is no small thing. A lot of code here is copied from ph.py
    155165      action_handles = []
    156166      scenario_action_handle_map = {} # maps scenario names to action handles
    157167      action_handle_scenario_map = {} # maps action handles to scenario names
    158       E_f_of_Hatx = 0
     168      E_f_of_Xhat = 0
    159169      # to support parallel, loop to launch then loop to get and use results
    160170      for scenario in  phforGk._scenario_tree._scenarios:
    161171         instance = phforGk._instances[scenario._name]
    162          # do the fixing at hatx
    163          fix_first_stage_vars_for_instance_from_PHAVG(hatxph, instance)
     172         # do the fixing at xhat
     173         fix_first_stage_vars_for_instance_from_PHAVG(xhatph, instance)
    164174
    165175         instance.preprocess()
     
    204214         objval = phforGk._scenario_tree.compute_scenario_cost(instance)
    205215         # JPW: THIS PROBABILITY IS BAD - SCENARIO OBJECT IS WRONG, FROM PREVIOUS LOOP.
    206          E_f_of_Hatx += objval * scenario._probability
    207 
    208       print "E_f_of_Hatx=",E_f_of_Hatx
     216         E_f_of_Xhat += objval * scenario._probability
     217
     218      print "E_f_of_Xhat=",E_f_of_Xhat
    209219      print "E_f_of_xstar_bound=",E_f_of_xstar_bound
    210220
    211       G_supk_of_hatx.append(E_f_of_Hatx - E_f_of_xstar_bound)   
    212       Gbar += G_supk_of_hatx[k-1]
    213       sum_E_f_of_Hatx += E_f_of_Hatx
     221      G_supk_of_xhat.append(E_f_of_Xhat - E_f_of_xstar_bound)   
     222      Gbar += G_supk_of_xhat[k-1]
     223      sum_E_f_of_Xhat += E_f_of_Xhat
    214224   Gbar = Gbar / n_g
    215225   # second pass for variance calculation (because we like storing the G_supk)
    216226   VarG = 0
    217227   for k in range(0, n_g-1):
    218       VarG = VarG + (G_supk_of_hatx[k] - Gbar) * (G_supk_of_hatx[k] - Gbar)
     228      VarG = VarG + (G_supk_of_xhat[k] - Gbar) * (G_supk_of_xhat[k] - Gbar)
    219229   VarG = VarG / (n_g - 1)    # sample var
    220230   print "Gbar=",Gbar
    221231   print "StdG=",math.sqrt(VarG)
    222    print "Average f(hatx)=",sum_E_f_of_Hatx / n_g
     232   print "Average f(xhat)=",sum_E_f_of_Xhat / n_g
    223233
    224234
    225235   if options.append_file is not None:
    226236      output_file = open(options.append_file, "a")
    227       print >>output_file, "\ninstancedirectory, ",options.instance_directory, ", seed, ",options.random_seed, ", N, ",scenariocount, ", hatn, ",num_scenarios_for_solution, ", n_g, ", options.n_g, ", EoffofHatx, ", sum_E_f_of_Hatx / n_g, ", gbar, ", Gbar, ", sg, ", math.sqrt(VarG), ", objforhatx, ", xhatobj, ", n,",num_scenarios_per_sample,
     237      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_E_f_of_Xhat / n_g, ", gbar, ", Gbar, ", sg, ", math.sqrt(VarG), ", objforxhat, ", xhatobj, ", n,",num_scenarios_per_sample,
    228238      output_file.close()
    229239   
    230240
    231241#==============================================   
    232 def ph_for_bundle(BundleStart, BundleStop, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, IndList, options):
     242def ph_for_bundle(BundleStart, BundleStop, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options):
    233243   
    234244   scenarios_to_bundle = []
    235245   for i in range(BundleStart, BundleStop+1):   # python has zero based indexes
    236       scenarios_to_bundle.append(full_scenario_tree._scenarios[IndList[i]]._name)
     246      scenarios_to_bundle.append(full_scenario_tree._scenarios[index_list[i]]._name)
    237247
    238248#   print "Scenarios to bundle:"+str(scenarios_to_bundle)
Note: See TracChangeset for help on using the changeset viewer.