Changeset 2919


Ignore:
Timestamp:
Aug 10, 2010 9:54:36 PM (9 years ago)
Author:
jwatson
Message:

Simplified the PySP confidence interval code, specifically using the extensive form library to leverage common routines.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • coopr.pysp/trunk/scripts/computeconf

    r2918 r2919  
    5757                                     dest="append_file",
    5858                                     type="string",
    59                                      default=None)     
     59                                     default=None)
     60      conf_options_parser.add_option("--generate-weighted-cvar",
     61                                     help="Add a weighted CVaR term to the primary objective",
     62                                     action="store_true",
     63                                     dest="generate_weighted_cvar",
     64                                     default=False)
     65      conf_options_parser.add_option("--cvar-weight",
     66                                     help="The weight associated with the CVaR term in the risk-weighted objective formulation. Default is 1.0. If the weight is 0, then *only* a non-weighted CVaR cost will appear in the EF objective - the expected cost component will be dropped.",
     67                                     action="store",
     68                                     dest="cvar_weight",
     69                                     type="float",
     70                                     default=1.0)
     71      conf_options_parser.add_option("--risk-alpha",
     72                                     help="The probability threshold associated with cvar (or any future) risk-oriented performance metrics. Default is 0.95.",
     73                                     action="store",
     74                                     dest="risk_alpha",
     75                                     type="float",
     76                                     default=0.95)           
    6077
    6178      (options, args) = conf_options_parser.parse_args(args=args)
     
    118135   if options.solve_xhat_with_ef_only is True:
    119136      print "Creating the xhat extensive form."
    120       hatex_ef = create_ef_instance(xhat_ph._scenario_tree, xhat_ph._instances)
     137      hatex_ef = create_ef_instance(xhat_ph._scenario_tree,
     138                                    xhat_ph._instances,
     139                                    generate_weighted_cvar=options.generate_weighted_cvar,
     140                                    cvar_weight=options.cvar_weight,
     141                                    risk_alpha=options.risk_alpha)                                   
    121142      print "Solving the xhat extensive form."
    122143      ef_results = write_and_solve_ef(hatex_ef, xhat_ph._instances, options)
     
    140161   g_supk_of_xhat = [] # really not always needed... http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
    141162   g_bar = 0
    142    sum_E_f_of_xhat = 0
     163   sum_xstar_obj_given_xhat = 0
    143164   
    144165   for k in range(1, n_g+1):
     166     
    145167      start_index = num_scenarios_for_solution + (k-1)*num_scenarios_per_sample
    146168      stop_index = start_index + num_scenarios_per_sample - 1
     
    153175      gk_ph = ph_for_bundle(start_index, stop_index, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options)
    154176      print "Creating the xstar extensive form."     
    155       gk_ef = create_ef_instance(gk_ph._scenario_tree, gk_ph._instances)
     177      gk_ef = create_ef_instance(gk_ph._scenario_tree,
     178                                 gk_ph._instances,
     179                                 generate_weighted_cvar=options.generate_weighted_cvar,
     180                                 cvar_weight=options.cvar_weight,
     181                                 risk_alpha=options.risk_alpha)
    156182      print "Solving the xstar extensive form."     
    157183      ef_results = write_and_solve_ef(gk_ef, gk_ph._instances, options)
    158184      load_ef_solution(ef_results, gk_ef, gk_ph._instances)
    159185
    160       # as in the computation of xhat, the following is required to form a solution
    161       # to the extensive form in the scenario tree itself.
     186      # as in the computation of xhat, the following is required to form a
     187      # solution to the extensive form in the scenario tree itself.
    162188      gk_ph._scenario_tree.snapshotSolutionFromInstances(gk_ph._instances)
    163189     
     
    169195      print "Extensive form objective value given xstar="+str(xstar_obj)+"; gap="+str(xstar_obj_gap)
    170196
    171       # to get fj(xhat) we need the obj value for each scenario, j, in the bundle evaluated at hat x
    172       # this is no small thing. A lot of code here is copied from ph.py
    173       action_handles = []
    174       scenario_action_handle_map = {} # maps scenario names to action handles
    175       action_handle_scenario_map = {} # maps action handles to scenario names
    176       xstar_obj_given_xhat = 0
    177       # to support parallel, loop to launch then loop to get and use results
     197      # to get f(xhat) for this sample, fix the first-stage variables and re-solve the extensive form.
     198      # note that the fixing yields side-effects on the original gk_ef, but that is fine as it isn't
     199      # used after this point.
     200      print "Solving the extensive form given the xhat solution."
    178201      for scenario in  gk_ph._scenario_tree._scenarios:
    179202         instance = gk_ph._instances[scenario._name]
    180          # do the fixing at xhat
    181203         fix_first_stage_vars(xhat_ph._scenario_tree, instance)
    182 
    183204         instance.preprocess()
    184 
    185          new_action_handle = gk_ph._solver_manager.queue(instance, opt=gk_ph._solver, tee=gk_ph._output_solver_log)
    186          scenario_action_handle_map[scenario._name] = new_action_handle
    187          action_handle_scenario_map[new_action_handle] = scenario._name
    188 
    189          action_handles.append(new_action_handle)
    190 
    191       # loop to get and use results
    192       num_results_so_far = 0
    193       while (num_results_so_far < len(gk_ph._scenario_tree._scenarios)):
    194 
    195          action_handle = gk_ph._solver_manager.wait_any()
    196          results = gk_ph._solver_manager.get_results(action_handle)         
    197          scenario_name = action_handle_scenario_map[action_handle]
    198          instance = gk_ph._instances[scenario_name]
    199 
    200          if gk_ph._verbose is True:
    201             print "Sampling Results obtained for scenario="+scenario_name
    202 
    203          if len(results.solution) == 0:
    204             results.write(num=1)
    205             raise RuntimeError, "Sampling Solve failed for scenario="+scenario_name+"; no solutions generated"
    206 
    207          if gk_ph._output_solver_results is True:
    208             print "Sampling Results for scenario=",scenario_name
    209             results.write(num=1)
    210 
    211          start_time = time.time()
    212          instance.load(results)
    213          end_time = time.time()
    214          if gk_ph._output_times is True:
    215             print "Time loading results into instance="+str(end_time-start_time)+" seconds"
    216 
    217          if gk_ph._verbose is True:                 
    218             print "Successfully loaded solution for scenario="+scenario_name
    219 
    220          num_results_so_far = num_results_so_far + 1
    221 
    222          objval = gk_ph._scenario_tree.compute_scenario_cost(instance)
    223          # JPW: THIS PROBABILITY IS BAD - SCENARIO OBJECT IS WRONG, FROM PREVIOUS LOOP.
    224          xstar_obj_given_xhat += objval * scenario._probability
    225 
     205      gk_ef.preprocess() # to account for the fixed variables in the scenario instances.
     206      ef_results = write_and_solve_ef(gk_ef, gk_ph._instances, options)
     207      load_ef_solution(ef_results, gk_ef, gk_ph._instances)
     208
     209      # we don't need the solution - just the objective value.
     210      xstar_obj_given_xhat = float(ef_results.solution(0).objective['f'].value)
    226211      print "Sample extensive form objective value given xhat="+str(xstar_obj_given_xhat)
    227       print "xstar_obj_bound=",xstar_obj_bound
    228212
    229213      g_supk_of_xhat.append(xstar_obj_given_xhat - xstar_obj_bound)   
    230214      g_bar += g_supk_of_xhat[k-1]
    231       sum_E_f_of_xhat += E_f_of_xhat
     215      sum_xstar_obj_given_xhat += xstar_obj_given_xhat
    232216     
    233217   g_bar = g_bar / n_g
    234218   # second pass for variance calculation (because we like storing the g_supk)
    235    VarG = 0
     219   g_var = 0
    236220   for k in range(0, n_g-1):
    237       VarG = VarG + (g_supk_of_xhat[k] - g_bar) * (g_supk_of_xhat[k] - g_bar)
    238    VarG = VarG / (n_g - 1)    # sample var
     221      g_var = g_var + (g_supk_of_xhat[k] - g_bar) * (g_supk_of_xhat[k] - g_bar)
     222   g_var = g_var / (n_g - 1)    # sample var
    239223   print "g_bar=",g_bar
    240    print "StdG=",math.sqrt(VarG)
    241    print "Average f(xhat)=",sum_E_f_of_xhat / n_g
    242 
     224   print "g_stddev=",math.sqrt(g_var)
     225   print "Average f(xhat)=",sum_xstar_obj_given_xhat / n_g
    243226
    244227   if options.append_file is not None:
    245228      output_file = open(options.append_file, "a")
    246       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, ", g_bar, ", sg, ", math.sqrt(VarG), ", objforxhat, ", xhat_obj, ", n,",num_scenarios_per_sample,
     229      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,
    247230      output_file.close()
    248231
Note: See TracChangeset for help on using the changeset viewer.