Changeset 2918


Ignore:
Timestamp:
Aug 10, 2010 8:55:47 PM (9 years ago)
Author:
jwatson
Message:

More cleanup on the PySP confidence interval computation script.

Location:
coopr.pysp/trunk
Files:
2 edited

Legend:

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

    r2917 r2918  
    417417   # validate the tree prior to doing anything serious
    418418   #
    419    print ""
    420419   if scenario_tree.validate() is False:
    421420      print "***ERROR: Scenario tree is invalid****"
     
    424423      if options.verbose is True:
    425424         print "Scenario tree is valid!"
    426    print ""
    427425
    428426   #
  • coopr.pysp/trunk/scripts/computeconf

    r2917 r2918  
    112112   # directly, mainly out of convenience - we're leveraging the code in ph_for_bundle to create the
    113113   # 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)
    115    xhatobj = None
     114   print "Loading scenario instances and initializing scenario tree for xhat scenario bundle."
     115   xhat_ph = ph_for_bundle(0, num_scenarios_for_solution, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options)
     116   xhat_obj = None
    116117
    117118   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)
    122       # the following is just for printing things in an organized manner.
    123       xhatph._scenario_tree.snapshotSolutionFromInstances(xhatph._instances)
    124       xhatobj = float(ef_results.solution(0).objective['f'].value)  ## DLW to JPW: we need the gap too
     119      print "Creating the xhat extensive form."
     120      hatex_ef = create_ef_instance(xhat_ph._scenario_tree, xhat_ph._instances)
     121      print "Solving the xhat extensive form."
     122      ef_results = write_and_solve_ef(hatex_ef, xhat_ph._instances, options)
     123      load_ef_solution(ef_results, hatex_ef, xhat_ph._instances)
     124      # IMPT: the following method populates the _solution variables on the scenario tree
     125      #       nodes by forming an average of the corresponding variable values for all
     126      #       instances particpating in that node. if you don't do this, the scenario tree
     127      #       doesn't have the solution - and we need this below for variable fixing.
     128      xhat_ph._scenario_tree.snapshotSolutionFromInstances(xhat_ph._instances)
     129      xhat_obj = float(ef_results.solution(0).objective['f'].value)  ## DLW to JPW: we need the gap too
     130      print "Extensive form objective value given xhat="+str(xhat_obj)
    125131   else:
    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"
    130 
    131    print "now form and solve the problems for each sample"
     132      print "Solving for xhat via Progressive Hedging."     
     133      xhat_ph.solve()
     134      # TBD - grab xhat_obj; complicated by the fact that PH may not have converged.
     135      # TBD - also note sure if PH calls snapshotSolutionFromAverages.
     136   print "The computation of xhat is complete - starting to compute confidence interval via sub-sampling."
     137
    132138   # in order to handle the case of scenarios that are not equally likely, we will split the expectations for Gsupk
    133139   # 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
    134    G_supk_of_xhat = [] # really not always needed... http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
    135    Gbar = 0
    136    sum_E_f_of_Xhat = 0
     140   g_supk_of_xhat = [] # really not always needed... http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
     141   g_bar = 0
     142   sum_E_f_of_xhat = 0
     143   
    137144   for k in range(1, n_g+1):
    138       start = num_scenarios_for_solution + (k-1)*num_scenarios_per_sample
    139       stop = start + num_scenarios_per_sample - 1
    140       print "Computing statistics for sample k="+str(k)
     145      start_index = num_scenarios_for_solution + (k-1)*num_scenarios_per_sample
     146      stop_index = start_index + num_scenarios_per_sample - 1
     147      print "Computing statistics for sample k="+str(k)+"."
     148      if options.verbose is True:
     149         print "Bundle start index="+str(start_index)+", stop index="+str(stop_index)+"."
    141150     
    142       # NOTE: We'll never run this ph - code could be refactored
    143 
    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)
    146       efforGk = create_ef_instance(phforGk._scenario_tree, phforGk._instances)     
    147       ef_results = write_and_solve_ef(efforGk, phforGk._instances, options)
    148       load_ef_solution(ef_results, efforGk, phforGk._instances)
    149 
    150       # the following is just for printing things in an organized manner.
    151       phforGk._scenario_tree.snapshotSolutionFromInstances(phforGk._instances)
    152 #      phforGk._scenario_tree.pprintSolution()
    153 #      phforGk._scenario_tree.pprintCosts(phforGk._instances)           
     151      # compute this xstar solution for the EF associated with sample k.
     152      print "Loading scenario instances and initializing scenario tree for xstar scenario bundle."     
     153      gk_ph = ph_for_bundle(start_index, stop_index, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options)
     154      print "Creating the xstar extensive form."     
     155      gk_ef = create_ef_instance(gk_ph._scenario_tree, gk_ph._instances)
     156      print "Solving the xstar extensive form."     
     157      ef_results = write_and_solve_ef(gk_ef, gk_ph._instances, options)
     158      load_ef_solution(ef_results, gk_ef, gk_ph._instances)
     159
     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.
     162      gk_ph._scenario_tree.snapshotSolutionFromInstances(gk_ph._instances)
    154163     
    155 #      print "HEY JPW!!! we need to add or subtract the gap"
    156 #      # JPW doesn't like the 'f' below - seems bad hard-coding a literal.
    157       E_f_of_xstar = float(ef_results.solution(0).objective['f'].value)  ## DLW to JPW: we need the gap too
    158       print "E_f_of_xstar=",E_f_of_xstar
    159       efforgK_gap = ef_results.solution(0).gap # assuming this is the absolute gap
    160       print "GAP=",efforgK_gap
    161       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?
     164      # extract the objective function value corresponding to the xstar solution, along with any gap information.
     165      # JPW doesn't like the 'f' below - seems bad hard-coding a literal.
     166      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.
     167      xstar_obj_gap = ef_results.solution(0).gap # assuming this is the absolute gap
     168      xstar_obj_bound = xstar_obj - xstar_obj_gap # HERE - watch for signs on gap - is it absolute, or can it be negative?     
     169      print "Extensive form objective value given xstar="+str(xstar_obj)+"; gap="+str(xstar_obj_gap)
    162170
    163171      # to get fj(xhat) we need the obj value for each scenario, j, in the bundle evaluated at hat x
     
    166174      scenario_action_handle_map = {} # maps scenario names to action handles
    167175      action_handle_scenario_map = {} # maps action handles to scenario names
    168       E_f_of_Xhat = 0
     176      xstar_obj_given_xhat = 0
    169177      # to support parallel, loop to launch then loop to get and use results
    170       for scenario in  phforGk._scenario_tree._scenarios:
    171          instance = phforGk._instances[scenario._name]
     178      for scenario in  gk_ph._scenario_tree._scenarios:
     179         instance = gk_ph._instances[scenario._name]
    172180         # do the fixing at xhat
    173          fix_first_stage_vars_for_instance_from_PHAVG(xhatph, instance)
     181         fix_first_stage_vars(xhat_ph._scenario_tree, instance)
    174182
    175183         instance.preprocess()
    176184
    177          new_action_handle = phforGk._solver_manager.queue(instance, opt=phforGk._solver, tee=phforGk._output_solver_log)
     185         new_action_handle = gk_ph._solver_manager.queue(instance, opt=gk_ph._solver, tee=gk_ph._output_solver_log)
    178186         scenario_action_handle_map[scenario._name] = new_action_handle
    179187         action_handle_scenario_map[new_action_handle] = scenario._name
     
    183191      # loop to get and use results
    184192      num_results_so_far = 0
    185       while (num_results_so_far < len(phforGk._scenario_tree._scenarios)):
    186 
    187          action_handle = phforGk._solver_manager.wait_any()
    188          results = phforGk._solver_manager.get_results(action_handle)         
     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)         
    189197         scenario_name = action_handle_scenario_map[action_handle]
    190          instance = phforGk._instances[scenario_name]
    191 
    192          if phforGk._verbose is True:
     198         instance = gk_ph._instances[scenario_name]
     199
     200         if gk_ph._verbose is True:
    193201            print "Sampling Results obtained for scenario="+scenario_name
    194202
     
    197205            raise RuntimeError, "Sampling Solve failed for scenario="+scenario_name+"; no solutions generated"
    198206
    199          if phforGk._output_solver_results is True:
     207         if gk_ph._output_solver_results is True:
    200208            print "Sampling Results for scenario=",scenario_name
    201209            results.write(num=1)
     
    204212         instance.load(results)
    205213         end_time = time.time()
    206          if phforGk._output_times is True:
     214         if gk_ph._output_times is True:
    207215            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
    208216
    209          if phforGk._verbose is True:                 
     217         if gk_ph._verbose is True:                 
    210218            print "Successfully loaded solution for scenario="+scenario_name
    211219
    212220         num_results_so_far = num_results_so_far + 1
    213221
    214          objval = phforGk._scenario_tree.compute_scenario_cost(instance)
     222         objval = gk_ph._scenario_tree.compute_scenario_cost(instance)
    215223         # JPW: THIS PROBABILITY IS BAD - SCENARIO OBJECT IS WRONG, FROM PREVIOUS LOOP.
    216          E_f_of_Xhat += objval * scenario._probability
    217 
    218       print "E_f_of_Xhat=",E_f_of_Xhat
    219       print "E_f_of_xstar_bound=",E_f_of_xstar_bound
    220 
    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
    224    Gbar = Gbar / n_g
    225    # second pass for variance calculation (because we like storing the G_supk)
     224         xstar_obj_given_xhat += objval * scenario._probability
     225
     226      print "Sample extensive form objective value given xhat="+str(xstar_obj_given_xhat)
     227      print "xstar_obj_bound=",xstar_obj_bound
     228
     229      g_supk_of_xhat.append(xstar_obj_given_xhat - xstar_obj_bound)   
     230      g_bar += g_supk_of_xhat[k-1]
     231      sum_E_f_of_xhat += E_f_of_xhat
     232     
     233   g_bar = g_bar / n_g
     234   # second pass for variance calculation (because we like storing the g_supk)
    226235   VarG = 0
    227236   for k in range(0, n_g-1):
    228       VarG = VarG + (G_supk_of_xhat[k] - Gbar) * (G_supk_of_xhat[k] - Gbar)
     237      VarG = VarG + (g_supk_of_xhat[k] - g_bar) * (g_supk_of_xhat[k] - g_bar)
    229238   VarG = VarG / (n_g - 1)    # sample var
    230    print "Gbar=",Gbar
     239   print "g_bar=",g_bar
    231240   print "StdG=",math.sqrt(VarG)
    232    print "Average f(xhat)=",sum_E_f_of_Xhat / n_g
     241   print "Average f(xhat)=",sum_E_f_of_xhat / n_g
    233242
    234243
    235244   if options.append_file is not None:
    236245      output_file = open(options.append_file, "a")
    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,
     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,
    238247      output_file.close()
    239    
    240 
    241 #==============================================   
    242 def ph_for_bundle(BundleStart, BundleStop, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options):
     248
     249#
     250# routine to create a down-sampled (bundled) scenario tree and the associated PH object.
     251#
     252def ph_for_bundle(bundle_start, bundle_stop, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options):
    243253   
    244254   scenarios_to_bundle = []
    245    for i in range(BundleStart, BundleStop+1):   # python has zero based indexes
     255   for i in range(bundle_start, bundle_stop+1):
    246256      scenarios_to_bundle.append(full_scenario_tree._scenarios[index_list[i]]._name)
    247257
    248 #   print "Scenarios to bundle:"+str(scenarios_to_bundle)
     258   if options.verbose is True:
     259      print "Creating PH object for scenario bundle="+str(scenarios_to_bundle)
    249260
    250261   scenario_tree_for_soln = ScenarioTree(scenarioinstance=reference_instance,
    251                             scenariotreeinstance=scenario_tree_instance,
    252                             scenariobundlelist=scenarios_to_bundle)
    253 
     262                                         scenariotreeinstance=scenario_tree_instance,
     263                                         scenariobundlelist=scenarios_to_bundle)
    254264   if scenario_tree_for_soln.validate() is False:
    255       print "***ERROR: Bundle Scenario tree is invalid****"
    256       sys.exit(0)
     265      print "***ERROR: Bundled scenario tree is invalid!!!"
     266      sys.exit(1)
     267
    257268   ph = create_ph_from_scratch(options, reference_model, reference_instance, scenario_tree_for_soln)
    258269   return ph
    259270
    260 #==============================================   
    261 def fix_first_stage_vars_for_instance_from_PHAVG(ph, instance):
    262 # TBD: CORRECT FOR PH VERSUS EF SOURCE
    263 # use the values from PHAVG of the ph object to fix the first stage vars in the instance
    264 #   stage = ph._scenario_tree._stages[0]   # root node only
    265    # use the first instance from PH because all should have the same PHAVG
    266 #   phinstance = ph._instances[ph._scenario_tree._scenarios[0]._name]
    267 #   phinstance.pprint()
    268 
    269 #   for (variable, index_template, variable_indices) in stage._variables:
    270 
    271 #       variable_name = variable.name
    272 #       variable_type = variable.domain
    273 #       print "FIXING VARIABLE=",variable_name
    274        # HERE: we want to not take anything from PHAVG, but rather from the scenario tree itself - the
    275        #       only difficulty is finding the root node object.
    276 #       avg_parameter_name = "PHAVG_"+variable_name
    277 #       avg_parameter = getattr(phinstance, avg_parameter_name)
    278 
    279 #       for index in variable_indices:
    280 #          print "INDEX=",index
    281 #          if getattr(phinstance,variable_name)[index].status != VarStatus.unused:
    282 #             fix_value = avg_parameter[index]()
    283 #             print "FIX VALUE=",fix_value
    284 #             if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
    285 #                fix_value = int(round(fix_value))
    286 #             getattr(instance,variable.name)[index].fixed = True
    287 #             getattr(instance,variable.name)[index].value = fix_value
    288 
    289    phinstance = ph._instances[ph._scenario_tree._scenarios[0]._name]
    290 
    291    scenario_tree = ph._scenario_tree
    292    stage = ph._scenario_tree._stages[0]
    293    root_node = stage._tree_nodes[0] # there should be only one
     271#
     272# fixes the first stage variables in the input instance to the values
     273# of the solution indicated in the input scenario tree.
     274#
     275def fix_first_stage_vars(scenario_tree, instance):
     276
     277   stage = scenario_tree._stages[0]
     278   root_node = stage._tree_nodes[0] # there should be only one root node!
    294279
    295280   for (variable, index_template, variable_indices) in stage._variables:
     
    301286
    302287      for index in variable_indices:
    303          if getattr(phinstance,variable_name)[index].status != VarStatus.unused:
     288         if getattr(instance, variable_name)[index].status != VarStatus.unused:
    304289            fix_value = scenario_tree_var[index]()
    305290            if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
     
    313298def write_and_solve_ef(master_instance, scenario_instances, options):
    314299
    315    ef_file_name = "dlwef.lp"
    316    write_ef(master_instance, scenario_instances, os.path.expanduser(ef_file_name))
    317 
    318    print ""
    319    print "Solving extensive form written to file="+os.path.expanduser(ef_file_name)
    320    print ""
     300   ef_file_name = os.path.expanduser("dlwef.lp")
     301   write_ef(master_instance, scenario_instances, ef_file_name)
    321302
    322303   ef_solver = SolverFactory(options.solver_type)
     
    336317      raise ValueError, "Failed to create solver manager of type="+options.solver_type+" for use in extensive form solve"
    337318
    338    print "Queuing extensive form solve"
    339    ef_action_handle = ef_solver_manager.queue(os.path.expanduser(ef_file_name), opt=ef_solver, warmstart=False, tee=options.output_ef_solver_log)
    340    print "Waiting for extensive form solve"
     319   ef_action_handle = ef_solver_manager.queue(ef_file_name, opt=ef_solver, warmstart=False, tee=options.output_ef_solver_log)
    341320   ef_results = ef_solver_manager.wait_for(ef_action_handle)
    342 #   print "Extensive form solve results:"
    343 #   ef_results.write(num=1)
     321   os.remove(ef_file_name)
     322
    344323   return ef_results
    345324
    346 #=================================================================================================
    347 # JPW: Any executable code can be placed at the end of a file, and it will be executed.
     325
     326#
     327# the main script routine starts here - and is obviously pretty simple!
     328#
     329
    348330run()
Note: See TracChangeset for help on using the changeset viewer.