Changeset 2938


Ignore:
Timestamp:
Aug 14, 2010 6:42:01 PM (9 years ago)
Author:
wehart
Message:

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

........

r2772 | wehart | 2010-07-05 15:29:34 -0600 (Mon, 05 Jul 2010) | 3 lines


Eliminating the direct use of cProfile, which is not
backwards compatible to Python 2.4

........

r2781 | jwatson | 2010-07-09 09:24:27 -0600 (Fri, 09 Jul 2010) | 3 lines


In the PySP farmer example, changing "MeanYield?" to "Yield" - MeanYield? is obviously incorrect (it's just Yield), as I discovered while writing the PySP journal article.

........

r2790 | khunter | 2010-07-13 16:16:46 -0600 (Tue, 13 Jul 2010) | 3 lines


NFC: Remove some whitespace, so it's clear what the *actual* changes
are next commit.

........

r2791 | khunter | 2010-07-13 16:18:38 -0600 (Tue, 13 Jul 2010) | 3 lines


Error message: Automatic cast to string in case, fixing possible
ugly traceback from an already-known error point.

........

r2792 | khunter | 2010-07-13 16:20:59 -0600 (Tue, 13 Jul 2010) | 3 lines


Fix bug where pprint dies if stage names are non-strings.
Basically, cast 'em to a string.

........

r2800 | khunter | 2010-07-14 15:00:10 -0600 (Wed, 14 Jul 2010) | 3 lines


NFC: remove EOL whitespace so as not to obfuscate actual work in
next commit

........

r2801 | khunter | 2010-07-14 15:02:15 -0600 (Wed, 14 Jul 2010) | 3 lines


Allow for the possibility of numeric names. (Fix pprint
assumption.) Similar to r2792

........

r2808 | khunter | 2010-07-15 17:11:18 -0600 (Thu, 15 Jul 2010) | 3 lines


No wording or pdf output changes, but reformat LaTeX so that it's
easier to see what changes in next commit.

........

r2809 | khunter | 2010-07-15 17:13:02 -0600 (Thu, 15 Jul 2010) | 3 lines


Minor grammar fixes, and resolution of a word-wrap issue for a file
name.

........

r2810 | khunter | 2010-07-15 17:16:21 -0600 (Thu, 15 Jul 2010) | 4 lines


Found an almost helpful error message while working on my model;
figure I can't be the only one, so add an error message that asks a
question about a perhaps noob stepping stone issue.

........

r2815 | khunter | 2010-07-19 17:35:13 -0600 (Mon, 19 Jul 2010) | 2 lines


NFC: remove EOL whitespace ...

........

r2816 | khunter | 2010-07-19 17:37:49 -0600 (Mon, 19 Jul 2010) | 2 lines


Remove unnecessary variable set and bound check.

........

r2817 | khunter | 2010-07-19 17:39:27 -0600 (Mon, 19 Jul 2010) | 5 lines


Ref #4093 "PySP missing import"


Fix for --linearize-nonbinary-penalty-terms where stage cost
variable is indexed.

........

r2818 | khunter | 2010-07-19 17:54:53 -0600 (Mon, 19 Jul 2010) | 3 lines


Make use of Python functions as first-class citizens to clean up
some uber long lines and add just smidgen of readability.

........

r2832 | jwatson | 2010-07-21 13:07:47 -0600 (Wed, 21 Jul 2010) | 3 lines


Adding some enhanced verbosity to the ef writer to aid debug.

........

r2833 | jwatson | 2010-07-21 13:54:51 -0600 (Wed, 21 Jul 2010) | 3 lines


Fixes to pretty-print and solution computation in the PySP scenario tree object to deal with nodes in the same stage that have heterogeneous indices for the same variable name. Comes up in EPA sensor placement formulation. Requires some real thought - runs for now.

........

r2834 | jwatson | 2010-07-21 15:03:15 -0600 (Wed, 21 Jul 2010) | 3 lines


Various updates to PySP computeconf while debugging/exercising the main script.

........

r2835 | jwatson | 2010-07-21 17:29:53 -0600 (Wed, 21 Jul 2010) | 3 lines


More computeconf fixes.

........

r2836 | jwatson | 2010-07-22 08:07:16 -0600 (Thu, 22 Jul 2010) | 3 lines


Eliminating a legacy debug output unconditionally generated when compressing a scenario tree.

........

r2837 | jwatson | 2010-07-22 09:12:02 -0600 (Thu, 22 Jul 2010) | 3 lines


Various touch-ups to compute conf PySP script - it looks good to go!

........

r2838 | jwatson | 2010-07-22 13:03:51 -0600 (Thu, 22 Jul 2010) | 3 lines


Fixed bug with random seed initialization.

........

r2842 | jwatson | 2010-07-23 15:36:23 -0600 (Fri, 23 Jul 2010) | 3 lines


Adding --keep-solver-files option to the PySP runef script, to aid debugging.

........

r2852 | jwatson | 2010-07-23 21:36:51 -0600 (Fri, 23 Jul 2010) | 3 lines


Enforced the PySP convention that underscore characters not be part of scenario names. A ValueError? is now thrown if such a scenario name is encountered. It is better to do this up-front, as otherwise it causes parse problems when reading solutions later on.

........

r2917 | jwatson | 2010-08-10 17:04:15 -0600 (Tue, 10 Aug 2010) | 3 lines


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

........

r2918 | jwatson | 2010-08-10 18:55:47 -0600 (Tue, 10 Aug 2010) | 3 lines


More cleanup on the PySP confidence interval computation script.

........

r2919 | jwatson | 2010-08-10 19:54:36 -0600 (Tue, 10 Aug 2010) | 3 lines


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

........

r2923 | jwatson | 2010-08-13 15:36:17 -0600 (Fri, 13 Aug 2010) | 3 lines


Removing a redundant preprocess() method invocation for scenario instances created via node-based initialization in Progressive Hedging.

........

r2931 | wehart | 2010-08-14 15:54:35 -0600 (Sat, 14 Aug 2010) | 2 lines


Categorizing slow PySP tests as 'nightly'

........

Location:
coopr.pysp/stable/2.4
Files:
26 edited

Legend:

Unmodified
Added
Removed
  • coopr.pysp/stable/2.4

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

    r2734 r2938  
    519519   master_scenario_model = None
    520520   master_scenario_instance = None
     521
     522   if verbose_output:
     523      print "Constructing reference model and instance"
    521524   
    522525   try:
     
    556559   from coopr.pysp.util.scenariomodels import scenario_tree_model
    557560
     561   if verbose_output:
     562      print "Constructing scenario tree instance"
     563
    558564   scenario_tree_instance = scenario_tree_model.create(instance_directory+os.sep+"ScenarioStructure.dat")
    559565
     
    561567   # construct the scenario tree
    562568   #
     569   if verbose_output:
     570      print "Constructing scenario tree object"
     571   
    563572   scenario_tree = ScenarioTree(scenarioinstance=master_scenario_instance,
    564573                                scenariotreeinstance=scenario_tree_instance)
  • coopr.pysp/stable/2.4/coopr/pysp/ef_writer_script.py

    r2734 r2938  
    1717import textwrap
    1818import traceback
    19 import cProfile
     19try:
     20    import cProfile as profile
     21except ImportError:
     22    import profile
    2023import pstats
    2124import gc
     
    9194                     default="serial")
    9295   parser.add_option("--solver-options",
    93                      help="Solver options for the extension form problem",
     96                     help="Solver options for the extension form problem.",
    9497                     action="append",
    9598                     dest="solver_options",
     
    97100                     default=[])
    98101   parser.add_option("--mipgap",
    99                      help="Specifies the mipgap for the EF solve",
     102                     help="Specifies the mipgap for the EF solve.",
    100103                     action="store",
    101104                     dest="mipgap",
     
    103106                     default=None)   
    104107   parser.add_option("--output-solver-log",
    105                      help="Output solver log during the extensive form solve",
     108                     help="Output solver log during the extensive form solve.",
    106109                     action="store_true",
    107110                     dest="output_solver_log",
     111                     default=False)
     112   parser.add_option("--keep-solver-files",
     113                     help="Retain temporary input and output files for solve.",
     114                     action="store_true",
     115                     dest="keep_solver_files",
    108116                     default=False)   
    109117   parser.add_option("--profile",
     
    157165         else:
    158166            ef_solver.mipgap = options.mipgap
     167      if options.keep_solver_files is True:
     168         ef_solver.keepFiles = True         
    159169
    160170      ef_solver_manager = SolverManagerFactory(options.solver_manager_type)
     
    209219        #
    210220        tfile = pyutilib.services.TempfileManager.create_tempfile(suffix=".profile")
    211         tmp = cProfile.runctx('run_ef_writer(options,args)',globals(),locals(),tfile)
     221        tmp = profile.runctx('run_ef_writer(options,args)',globals(),locals(),tfile)
    212222        p = pstats.Stats(tfile).strip_dirs()
    213223        p.sort_stats('time', 'cum')
  • coopr.pysp/stable/2.4/coopr/pysp/ph.py

    r2734 r2938  
    5353         variable_name = variable_value.name
    5454         variable_index = None
    55      
     55
    5656      new_rho_value = None
    5757      if isinstance(rho_expression, float):
     
    128128      self._solver_manager = None
    129129      self._solver = None
    130      
     130
    131131      checkpoint_file = open(checkpoint_filename, "w")
    132132      pickle.dump(self,checkpoint_file)
     
    138138
    139139      print "Checkpoint written to file="+checkpoint_filename
    140    
     140
    141141   #
    142142   # a simple utility to count the number of continuous and discrete variables in a set of instances.
     
    149149      num_continuous_vars = 0
    150150      num_discrete_vars = 0
    151      
     151
    152152      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage
    153153
     
    158158               for index in variable_indices:
    159159
    160                   is_used = True # until proven otherwise                     
     160                  is_used = True # until proven otherwise
    161161                  for scenario in tree_node._scenarios:
    162162                     instance = self._instances[scenario._name]
     
    164164                        is_used = False
    165165
    166                   if is_used is True:                       
     166                  if is_used is True:
    167167
    168168                     if isinstance(variable.domain, IntegerSet) or isinstance(variable.domain, BooleanSet):
     
    185185      num_fixed_continuous_vars = 0
    186186      num_fixed_discrete_vars = 0
    187      
     187
    188188      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage
    189          
     189
    190190         for tree_node in stage._tree_nodes:
    191191
     
    195195
    196196                  # implicit assumption is that if a variable value is fixed in one
    197                   # scenario, it is fixed in all scenarios. 
     197                  # scenario, it is fixed in all scenarios.
    198198
    199199                  is_fixed = False # until proven otherwise
     
    209209                        num_fixed_discrete_vars = num_fixed_discrete_vars + 1
    210210                     else:
    211                         num_fixed_continuous_vars = num_fixed_continuous_vars + 1                           
     211                        num_fixed_continuous_vars = num_fixed_continuous_vars + 1
    212212
    213213      return (num_fixed_discrete_vars, num_fixed_continuous_vars)
     
    216216   # we end up (necessarily) "littering" the scenario instances with extra constraints.
    217217   # these need to and should be cleaned up after PH, for purposes of post-PH manipulation,
    218    # e.g., writing the extensive form. 
     218   # e.g., writing the extensive form.
    219219   def _cleanup_scenario_instances(self):
    220220
     
    228228         instance.preprocess()
    229229
    230    # create PH weight and xbar vectors, on a per-scenario basis, for each variable that is not in the 
     230   # create PH weight and xbar vectors, on a per-scenario basis, for each variable that is not in the
    231231   # final stage, i.e., for all variables that are being blended by PH. the parameters are created
    232232   # in the space of each scenario instance, so that they can be directly and automatically
     
    248248      num_values = 0
    249249      first_time = True
    250      
     250
    251251      first_stage = self._scenario_tree._stages[0]
    252252      (cost_variable, cost_variable_index) = first_stage._cost_variable
     
    283283               weights_to_transmit.append(getattr(scenario_instance, weight_parameter_name))
    284284               average_parameter_name = "PHAVG_"+variable_name
    285                averages_to_transmit.append(getattr(scenario_instance, average_parameter_name))               
     285               averages_to_transmit.append(getattr(scenario_instance, average_parameter_name))
    286286
    287287         self._solver_manager.transmit_weights_and_averages(scenario_instance, weights_to_transmit, averages_to_transmit)
     
    312312   # a utility to transmit - across the PH solver manager - the current scenario
    313313   # tree node statistics to each of my problem instances. done prior to each
    314    # PH iteration k. 
     314   # PH iteration k.
    315315   #
    316316
     
    329329            tree_node_maximums[tree_node._name] = tree_node._maximums
    330330
    331          self._solver_manager.transmit_tree_node_statistics(scenario_instance, tree_node_minimums, tree_node_maximums)         
     331         self._solver_manager.transmit_tree_node_statistics(scenario_instance, tree_node_minimums, tree_node_maximums)
    332332
    333333   #
     
    340340
    341341         self._solver_manager.enable_ph_objective(scenario_instance)
    342          
     342
    343343
    344344   """ Constructor
     
    369369      self._report_solutions = False # do I report solutions after each PH iteration?
    370370      self._report_weights = False # do I report PH weights prior to each PH iteration?
    371       self._report_only_statistics = False # do I report only variable statistics when outputting solutions and weights? 
     371      self._report_only_statistics = False # do I report only variable statistics when outputting solutions and weights?
    372372      self._output_continuous_variable_stats = True # when in verbose mode, do I output weights/averages for continuous variables?
    373373      self._output_solver_results = False
     
    383383      self._solver_type = "cplex"
    384384      self._solver_manager_type = "serial" # serial or pyro are the options currently available
    385      
    386       self._solver = None 
    387       self._solver_manager = None 
    388      
     385
     386      self._solver = None
     387      self._solver_manager = None
     388
    389389      self._keep_solver_files = False
    390390      self._output_solver_log = False
     
    405405
    406406      self._scenario_tree = None
    407      
     407
    408408      self._scenario_data_directory = "" # this the prefix for all scenario data
    409409      self._instances = {} # maps scenario name to the corresponding model instance
     
    466466            self._rho_setter = kwds[key]
    467467         elif key == "bounds_setter":
    468             self._bounds_setter = kwds[key]                       
     468            self._bounds_setter = kwds[key]
    469469         elif key == "solver":
    470470            self._solver_type = kwds[key]
     
    474474            scenario_solver_options = kwds[key]
    475475         elif key == "scenario_mipgap":
    476             self._mipgap = kwds[key]           
     476            self._mipgap = kwds[key]
    477477         elif key == "keep_solver_files":
    478478            self._keep_solver_files = kwds[key]
     
    488488            self._report_weights = kwds[key]
    489489         elif key == "report_only_statistics":
    490             self._report_only_statistics = kwds[key]                                   
     490            self._report_only_statistics = kwds[key]
    491491         elif key == "output_times":
    492492            self._output_times = kwds[key]
     
    504504            self._checkpoint_interval = kwds[key]
    505505         elif key == "output_scenario_tree_solution":
    506             self._output_scenario_tree_solution = kwds[key]                       
     506            self._output_scenario_tree_solution = kwds[key]
    507507         else:
    508508            print "Unknown option=" + key + " specified in call to PH constructor"
     
    518518      # validate the linearization (number of pieces) and breakpoint distribution parameters.
    519519      if self._linearize_nonbinary_penalty_terms < 0:
    520          raise ValueError, "Value of linearization parameter for nonbinary penalty terms must be non-negative; value specified=" + `self._linearize_nonbinary_penalty_terms`         
     520         raise ValueError, "Value of linearization parameter for nonbinary penalty terms must be non-negative; value specified=" + `self._linearize_nonbinary_penalty_terms`
    521521      if self._breakpoint_strategy < 0:
    522522         raise ValueError, "Value of the breakpoint distribution strategy parameter must be non-negative; value specified=" + str(self._breakpoint_strategy)
    523523      if self._breakpoint_strategy > 3:
    524524         raise ValueError, "Unknown breakpoint distribution strategy specified - valid values are between 0 and 2, inclusive; value specified=" + str(self._breakpoint_strategy)
    525    
     525
    526526      # validate rho setter file if specified.
    527527      if self._rho_setter is not None:
     
    532532      if self._bounds_setter is not None:
    533533         if os.path.exists(self._bounds_setter) is False:
    534             raise ValueError, "The bounds setter script file="+self._bounds_setter+" does not exist"     
     534            raise ValueError, "The bounds setter script file="+self._bounds_setter+" does not exist"
    535535
    536536      # validate the checkpoint interval.
     
    540540      # construct the sub-problem solver.
    541541      if self._verbose is True:
    542          print "Constructing solver type="+self._solver_type         
     542         print "Constructing solver type="+self._solver_type
    543543      self._solver = SolverFactory(self._solver_type)
    544544      if self._solver == None:
     
    563563      self._iteration_index_set = Set(name="PHIterations")
    564564      for i in range(0,self._max_iterations + 1):
    565          self._iteration_index_set.add(i)     
     565         self._iteration_index_set.add(i)
    566566
    567567      # spit out parameterization if verbosity is enabled
     
    573573            print "   Rho initialization file=" + self._rho_setter
    574574         if self._bounds_setter is not None:
    575             print "   Variable bounds initialization file=" + self._bounds_setter           
     575            print "   Variable bounds initialization file=" + self._bounds_setter
    576576         print "   Sub-problem solver type=" + `self._solver_type`
    577577         print "   Solver manager type=" + `self._solver_manager_type`
     
    620620
    621621      # construct the instances for each scenario.
    622       #       
     622      #
    623623      # garbage collection noticeably slows down PH when dealing with
    624624      # large numbers of scenarios. disable prior to instance construction,
     
    627627      re_enable_gc = gc.isenabled()
    628628      gc.disable()
    629      
     629
    630630      if self._verbose is True:
    631631         if self._scenario_tree._scenario_based_data == 1:
     
    633633         else:
    634634            print "Node-based instance initialization enabled"
    635          
     635
    636636      for scenario in self._scenario_tree._scenarios:
    637637
     
    655655      # the scenario tree prior to creating PH-related parameters,
    656656      # variables, and the like.
    657       for plugin in self._ph_plugins:     
    658          plugin.post_instance_creation(self)         
     657      for plugin in self._ph_plugins:
     658         plugin.post_instance_creation(self)
    659659
    660660      # create ph-specific parameters (weights, xbar, etc.) for each instance.
     
    664664
    665665      self._create_ph_scenario_parameters()
    666            
     666
    667667      # if specified, run the user script to initialize variable rhos at their whim.
    668668      if self._rho_setter is not None:
     
    699699            for tree_node in stage._tree_nodes:
    700700
    701                new_min_index = reference_variable._index 
     701               new_min_index = reference_variable._index
    702702               new_min_parameter_name = "NODEMIN_"+reference_variable.name
    703703               # this bit of ugliness is due to Pyomo not correctly handling the Param construction
     
    711711               for index in new_min_index:
    712712                  new_min_parameter[index] = 0.0
    713                tree_node._minimums[reference_variable.name] = new_min_parameter                                   
    714 
    715                new_avg_index = reference_variable._index 
     713               tree_node._minimums[reference_variable.name] = new_min_parameter
     714
     715               new_avg_index = reference_variable._index
    716716               new_avg_parameter_name = "NODEAVG_"+reference_variable.name
    717717               new_avg_parameter = None
     
    721721                  new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
    722722               for index in new_avg_index:
    723                   new_avg_parameter[index] = 0.0                     
     723                  new_avg_parameter[index] = 0.0
    724724               tree_node._averages[reference_variable.name] = new_avg_parameter
    725725
    726                new_max_index = reference_variable._index 
     726               new_max_index = reference_variable._index
    727727               new_max_parameter_name = "NODEMAX_"+reference_variable.name
    728728               new_max_parameter = None
     
    732732                  new_max_parameter = Param(new_max_index,name=new_max_parameter_name)
    733733               for index in new_max_index:
    734                   new_max_parameter[index] = 0.0                     
     734                  new_max_parameter[index] = 0.0
    735735               tree_node._maximums[reference_variable.name] = new_max_parameter
    736736
    737737      # the objective functions are modified throughout the course of PH iterations.
    738738      # save the original, as a baseline to modify in subsequent iterations. reserve
    739       # the original objectives, for subsequent modification. 
     739      # the original objectives, for subsequent modification.
    740740      self._original_objective_expression = {}
    741741      for instance_name, instance in self._instances.items():
     
    778778
    779779   """ Perform the non-weighted scenario solves and form the initial w and xbars.
    780    """   
     780   """
    781781   def iteration_0_solve(self):
    782782
     
    800800
    801801      for scenario in self._scenario_tree._scenarios:
    802          
     802
    803803         instance = self._instances[scenario._name]
    804          
     804
    805805         if self._verbose is True:
    806806            print "Queuing solve for scenario=" + scenario._name
     
    812812         #       to the base scenario instances - and the variable values/etc. need to be collectged.
    813813         instance.preprocess()
    814            
     814
    815815         # there's nothing to warm-start from in iteration 0, so don't include the keyword in the solve call.
    816          # the reason you don't want to include it is that some solvers don't know how to handle the keyword 
     816         # the reason you don't want to include it is that some solvers don't know how to handle the keyword
    817817         # at all (despite it being false). you might want to solve iteration 0 solves using some other solver.
    818818
     
    834834
    835835         action_handle = self._solver_manager.wait_any()
    836          results = self._solver_manager.get_results(action_handle)         
     836         results = self._solver_manager.get_results(action_handle)
    837837         scenario_name = action_handle_scenario_map[action_handle]
    838          instance = self._instances[scenario_name]         
     838         instance = self._instances[scenario_name]
    839839
    840840         if self._verbose is True:
     
    855855            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
    856856
    857          if self._verbose is True:                 
     857         if self._verbose is True:
    858858            print "Successfully loaded solution for scenario="+scenario_name
    859859
    860860         num_results_so_far = num_results_so_far + 1
    861          
    862       if self._verbose is True:     
    863          print "Scenario sub-problem solves completed"     
     861
     862      if self._verbose is True:
     863         print "Scenario sub-problem solves completed"
    864864
    865865      solve_end_time = time.time()
     
    888888
    889889      start_time = time.time()
    890      
     890
    891891      # compute statistics over all stages, even the last. this is necessary in order to
    892892      # successfully snapshot a scenario tree solution from the average values.
    893893      for stage in self._scenario_tree._stages:
    894          
     894
    895895         for tree_node in stage._tree_nodes:
    896                
     896
    897897            for (variable, index_template, variable_indices) in stage._variables:
    898                  
     898
    899899               variable_name = variable.name
    900                      
     900
    901901               avg_parameter_name = "PHAVG_"+variable_name
    902                  
     902
    903903               for index in variable_indices:
    904904                  min = float("inf")
     
    906906                  max = float("-inf")
    907907                  node_probability = 0.0
    908                      
    909                   is_used = True # until proven otherwise                     
     908
     909                  is_used = True # until proven otherwise
    910910                  for scenario in tree_node._scenarios:
    911                        
     911
    912912                     instance = self._instances[scenario._name]
    913                        
     913
    914914                     if getattr(instance,variable_name)[index].status == VarStatus.unused:
    915915                        is_used = False
    916                      else:                       
     916                     else:
    917917                        node_probability += scenario._probability
    918918                        var_value = getattr(instance, variable.name)[index].value
     
    927927                     tree_node._minimums[variable.name][index] = min
    928928                     tree_node._averages[variable.name][index] = avg / node_probability
    929                      tree_node._maximums[variable.name][index] = max                       
     929                     tree_node._maximums[variable.name][index] = max
    930930
    931931                     # distribute the newly computed average to the xbar variable in
     
    947947      # because the weight updates rely on the xbars, and the xbars are node-based,
    948948      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
    949       start_time = time.time()     
     949      start_time = time.time()
    950950
    951951      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.
    952          
     952
    953953         for tree_node in stage._tree_nodes:
    954954
     
    974974                        # the variable toward the mean, the weights will blow up and be huge by the
    975975                        # time that blending is activated.
    976                         variable_blend_indicator = getattr(instance, blend_parameter_name)[index]()                             
     976                        variable_blend_indicator = getattr(instance, blend_parameter_name)[index]()
    977977
    978978                        # get the weight and rho parameters for this variable/index combination.
    979979                        rho_value = getattr(instance, rho_parameter_name)[index]()
    980980                        current_variable_weight = getattr(instance, weight_parameter_name)[index]()
    981                                                
     981
    982982                        # if I'm maximizing, invert value prior to adding (hack to implement negatives).
    983983                        # probably fixed in Pyomo at this point - I just haven't checked in a long while.
     
    10141014
    10151015      if self._verbose is True:
    1016          print "------------------------------------------------"       
     1016         print "------------------------------------------------"
    10171017         print "Starting PH iteration " + str(self._current_iteration) + " solves"
    10181018
     
    10331033
    10341034      # STEP 0: set up all global solver options.
    1035       self._solver.mipgap = self._mipgap     
    1036 
    1037       # STEP 1: queue up the solves for all scenario sub-problems and 
     1035      self._solver.mipgap = self._mipgap
     1036
     1037      # STEP 1: queue up the solves for all scenario sub-problems and
    10381038      #         grab all of the action handles for the subsequent barrier sync.
    10391039
     
    10421042      action_handle_scenario_map = {} # maps action handles to scenario names
    10431043
    1044       for scenario in self._scenario_tree._scenarios:     
     1044      for scenario in self._scenario_tree._scenarios:
    10451045
    10461046         instance = self._instances[scenario._name]
     
    10601060            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
    10611061         else:
    1062             new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)           
     1062            new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
    10631063
    10641064         scenario_action_handle_map[scenario._name] = new_action_handle
    1065          action_handle_scenario_map[new_action_handle] = scenario._name         
     1065         action_handle_scenario_map[new_action_handle] = scenario._name
    10661066
    10671067         action_handles.append(new_action_handle)
     
    10771077
    10781078         action_handle = self._solver_manager.wait_any()
    1079          results = self._solver_manager.get_results(action_handle)         
     1079         results = self._solver_manager.get_results(action_handle)
    10801080         scenario_name = action_handle_scenario_map[action_handle]
    1081          instance = self._instances[scenario_name]         
     1081         instance = self._instances[scenario_name]
    10821082
    10831083         if self._verbose is True:
     
    10981098            print "Time loading results into instance="+str(end_time-start_time)+" seconds"
    10991099
    1100          if self._verbose is True:                 
     1100         if self._verbose is True:
    11011101            print "Successfully loaded solution for scenario="+scenario_name
    11021102
     
    11071107
    11081108         num_results_so_far = num_results_so_far + 1
    1109          
    1110       if self._verbose is True:     
    1111          print "Scenario sub-problem solves completed"     
     1109
     1110      if self._verbose is True:
     1111         print "Scenario sub-problem solves completed"
    11121112
    11131113      solve_end_time = time.time()
     
    11501150
    11511151      # update variable statistics prior to any output.
    1152       self.update_variable_statistics()     
     1152      self.update_variable_statistics()
    11531153
    11541154      if (self._verbose is True) or (self._report_solutions is True):
     
    11571157
    11581158      # let plugins know if they care.
    1159       for plugin in self._ph_plugins:     
     1159      for plugin in self._ph_plugins:
    11601160         plugin.post_iteration_0_solves(self)
    11611161
    11621162      # update the fixed variable statistics.
    1163       (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
     1163      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
    11641164
    11651165      if self._verbose is True:
     
    11701170      self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
    11711171      first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()
    1172       print "Convergence metric=%12.4f  First stage cost avg=%12.4f  Max-Min=%8.2f" % (self._converger.lastMetric(), first_stage_avg, first_stage_max-first_stage_min)               
     1172      print "Convergence metric=%12.4f  First stage cost avg=%12.4f  Max-Min=%8.2f" % (self._converger.lastMetric(), first_stage_avg, first_stage_max-first_stage_min)
    11731173
    11741174      self.update_weights()
    11751175
    11761176      # let plugins know if they care.
    1177       for plugin in self._ph_plugins:     
     1177      for plugin in self._ph_plugins:
    11781178         plugin.post_iteration_0(self)
    11791179
     
    12001200      for i in range(1, self._max_iterations+1):
    12011201
    1202          self._current_iteration = self._current_iteration + 1                 
    1203 
    1204          print "Initiating PH iteration=" + `self._current_iteration`         
     1202         self._current_iteration = self._current_iteration + 1
     1203
     1204         print "Initiating PH iteration=" + `self._current_iteration`
    12051205
    12061206         if (self._verbose is True) or (self._report_weights is True):
     
    12221222         # update variable statistics prior to any output.
    12231223         self.update_variable_statistics()
    1224          
     1224
    12251225         if (self._verbose is True) or (self._report_solutions is True):
    12261226            print "Variable values following scenario solves:"
     
    12281228
    12291229         # we don't technically have to do this at the last iteration,
    1230          # but with checkpointing and re-starts, you're never sure 
     1230         # but with checkpointing and re-starts, you're never sure
    12311231         # when you're executing the last iteration.
    12321232         self.update_weights()
     
    12371237
    12381238         # update the fixed variable statistics.
    1239          (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
     1239         (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
    12401240
    12411241         if self._verbose is True:
     
    12551255         # check for early termination.
    12561256         self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
    1257          first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()         
    1258          print "Convergence metric=%12.4f  First stage cost avg=%12.4f  Max-Min=%8.2f" % (self._converger.lastMetric(), first_stage_avg, first_stage_max-first_stage_min)         
     1257         first_stage_min, first_stage_avg, first_stage_max = self._extract_first_stage_cost_statistics()
     1258         print "Convergence metric=%12.4f  First stage cost avg=%12.4f  Max-Min=%8.2f" % (self._converger.lastMetric(), first_stage_avg, first_stage_max-first_stage_min)
    12591259
    12601260         if self._converger.isConverged(self) is True:
     
    12621262               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)
    12631263            else:
    1264                print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)+" or all discrete variables are fixed"               
     1264               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)+" or all discrete variables are fixed"
    12651265            break
    12661266
     
    12801280      if self._verbose is True:
    12811281         print "Number of discrete variables fixed before final plugin calls="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
    1282          print "Number of continuous variables fixed before final plugin calls="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"           
     1282         print "Number of continuous variables fixed before final plugin calls="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
    12831283
    12841284      # let plugins know if they care. do this before
     
    12911291      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
    12921292
    1293       self._solve_end_time = time.time()     
     1293      self._solve_end_time = time.time()
    12941294
    12951295      print "PH complete"
     
    12991299
    13001300      print "Final number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
    1301       print "Final number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"         
     1301      print "Final number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
    13021302
    13031303      print "Final variable values:"
    1304       self.pprint(False, False, True, True, output_only_statistics=self._report_only_statistics)         
     1304      self.pprint(False, False, True, True, output_only_statistics=self._report_only_statistics)
    13051305
    13061306      print "Final costs:"
     
    13431343      if self._output_continuous_variable_stats is False:
    13441344
    1345          variable_type = variable.domain         
     1345         variable_type = variable.domain
    13461346
    13471347         if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):
     
    13501350
    13511351      return True
    1352      
     1352
    13531353   #
    13541354   # pretty-prints the state of the current variable averages, weights, and values.
     
    13581358
    13591359      if self._initialized is False:
    1360          raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"         
    1361      
     1360         raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"
     1361
    13621362      # print tree nodes and associated variable/xbar/ph information in stage-order
    13631363      # we don't blend in the last stage, so we don't current care about printing the associated information.
    13641364      for stage in self._scenario_tree._stages[:-1]:
    13651365
    1366          print "\tStage=" + stage._name
     1366         print "\tStage=" + str(stage._name)
    13671367
    13681368         num_outputs_this_stage = 0 # tracks the number of outputs on a per-index basis.
     
    13761376               num_outputs_this_variable = 0 # track, so we don't output the variable names unless there is an entry to report.
    13771377
    1378                for index in variable_indices:               
     1378               for index in variable_indices:
    13791379
    13801380                  weight_parameter_name = "PHWEIGHT_"+variable_name
     
    13821382                  num_outputs_this_index = 0 # track, so we don't output the variable index more than once.
    13831383
    1384                   for tree_node in stage._tree_nodes:                 
     1384                  for tree_node in stage._tree_nodes:
    13851385
    13861386                     # determine if the variable/index pair is used across the set of scenarios (technically,
     
    13991399                        if variable_value.fixed is True:
    14001400                           is_fixed = True
    1401  
     1401
    14021402                     # IMPT: this is far from obvious, but variables that are fixed will - because
    14031403                     #       presolve will identify them as constants and eliminate them from all
     
    14221422                              (fabs(maximum_value) > self._integer_tolerance):
    14231423
    1424                               num_outputs_this_stage = num_outputs_this_stage + 1                           
     1424                              num_outputs_this_stage = num_outputs_this_stage + 1
    14251425                              num_outputs_this_variable = num_outputs_this_variable + 1
    14261426                              num_outputs_this_index = num_outputs_this_index + 1
     
    14371437                                 print "\t\t\t\tTree Node="+tree_node._name,
    14381438                              if output_only_statistics is False:
    1439                                  print "\t\t (Scenarios: ",                             
     1439                                 print "\t\t (Scenarios: ",
    14401440                                 for scenario in tree_node._scenarios:
    14411441                                    print scenario._name," ",
    14421442                                    if scenario == tree_node._scenarios[-1]:
    14431443                                       print ")"
    1444                            
     1444
    14451445                              if output_values is True:
    14461446                                 average_value = tree_node._averages[variable_name][index]()
    14471447                                 if output_only_statistics is False:
    1448                                     print "\t\t\t\tValues: ",                       
     1448                                    print "\t\t\t\tValues: ",
    14491449                                 for scenario in tree_node._scenarios:
    14501450                                    instance = self._instances[scenario._name]
     
    14581458                                          # is to avoid updating our regression test baseline output.
    14591459                                          print "    Min=%12.4f" % (minimum_value),
    1460                                           print "    Avg=%12.4f" % (average_value),                                         
     1460                                          print "    Avg=%12.4f" % (average_value),
    14611461                                          print "    Max=%12.4f" % (maximum_value),
    14621462                                       else:
     
    14851485            print "\t\tCost Variable=" + cost_variable_name
    14861486         else:
    1487             print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)           
     1487            print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)
    14881488         for tree_node in stage._tree_nodes:
    14891489            print "\t\t\tTree Node=" + tree_node._name,
     
    15221522                   if output_only_statistics is True:
    15231523                      print "    Min=%12.4f" % (minimum_value),
    1524                       print "    Avg=%12.4f" % (sum_values/num_values),                                         
     1524                      print "    Avg=%12.4f" % (sum_values/num_values),
    15251525                      print "    Max=%12.4f" % (maximum_value),
    15261526                   else:
     
    15281528                      print "    Avg=%12.4f" % (sum_values/num_values),
    15291529                   print ""
    1530            
     1530
  • coopr.pysp/stable/2.4/coopr/pysp/phinit.py

    r2461 r2938  
    2323
    2424# for profiling
    25 import cProfile
     25try:
     26    import cProfile as profile
     27except ImportError:
     28    import profile
    2629import pstats
    2730
     
    293296#
    294297# 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).
    295301#
    296302
     
    411417   # validate the tree prior to doing anything serious
    412418   #
    413    print ""
    414419   if scenario_tree.validate() is False:
    415420      print "***ERROR: Scenario tree is invalid****"
     
    418423      if options.verbose is True:
    419424         print "Scenario tree is valid!"
    420    print ""
    421425
    422426   #
     
    630634        #
    631635        tfile = pyutilib.services.TempfileManager.create_tempfile(suffix=".profile")
    632         tmp = cProfile.runctx('exec_ph(options)',globals(),locals(),tfile)
     636        tmp = profile.runctx('exec_ph(options)',globals(),locals(),tfile)
    633637        p = pstats.Stats(tfile).strip_dirs()
    634638        p.sort_stats('time', 'cum')
  • coopr.pysp/stable/2.4/coopr/pysp/phobjective.py

    r2410 r2938  
    1414from math import fabs, log, exp
    1515from coopr.pyomo import *
     16from phutils import indexToString
    1617
    1718# IMPT: In general, the breakpoint computation codes can return a 2-list even if the lb equals
     
    2425
    2526def compute_uniform_breakpoints(lb, node_min, xavg, node_max, ub, num_breakpoints_per_side, tolerance):
    26    
     27
    2728   breakpoints = []
    2829
     
    3839      if (fabs(this_lb-lb) > tolerance) and (fabs(this_lb-xavg) > tolerance):
    3940         breakpoints.append(this_lb)
    40       current_x += left_step           
     41      current_x += left_step
    4142
    4243   # add the mean - it's always a breakpoint. unless!
     
    5152      this_lb = current_x
    5253      this_ub = current_x+right_step
    53       if (fabs(this_ub-xavg) > tolerance) and (fabs(this_ub-ub) > tolerance):         
     54      if (fabs(this_ub-xavg) > tolerance) and (fabs(this_ub-ub) > tolerance):
    5455         breakpoints.append(this_ub)
    5556      current_x += right_step
     
    6465#
    6566# routine to compute linearization breakpoints uniformly between the current node min/max bounds.
    66 #   
     67#
    6768
    6869def compute_uniform_between_nodestat_breakpoints(lb, node_min, xavg, node_max, ub, num_breakpoints, tolerance):
     
    7475
    7576   # add the node-min - the second breakpoint. but only if it is different than the lower bound and the mean.
    76    if (fabs(node_min-lb) > tolerance) and (fabs(node_min-xavg) > tolerance):     
     77   if (fabs(node_min-lb) > tolerance) and (fabs(node_min-xavg) > tolerance):
    7778      breakpoints.append(node_min)
    7879
     
    8485      if (fabs(this_lb-node_min) > tolerance) and (fabs(this_lb-node_max) > tolerance) and (fabs(this_lb-xavg) > tolerance):
    8586         breakpoints.append(this_lb)
    86       current_x += step           
     87      current_x += step
    8788
    8889   # add the node-max - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
    89    if (fabs(node_max-ub) > tolerance) and (fabs(node_max-xavg) > tolerance):           
    90       breakpoints.append(node_max)     
     90   if (fabs(node_max-ub) > tolerance) and (fabs(node_max-xavg) > tolerance):
     91      breakpoints.append(node_max)
    9192
    9293   # add the upper bound - the last breakpoint.
     
    104105#
    105106# routine to compute linearization breakpoints using "Woodruff" relaxation of the compute_uniform_between_nodestat_breakpoints.
    106 #   
     107#
    107108
    108109def compute_uniform_between_woodruff_breakpoints(lb, node_min, xavg, node_max, ub, num_breakpoints, tolerance):
     
    115116   # be either three "differences" from the mean, or else "halfway to the bound", whichever is closer to the mean.
    116117   left = max(xavg - 3.0 * (xavg - node_min), xavg - 0.5 * (xavg - lb))
    117    right = min(xavg + 3.0 * (node_max - xavg), xavg + 0.5 * (ub - xavg))     
     118   right = min(xavg + 3.0 * (node_max - xavg), xavg + 0.5 * (ub - xavg))
    118119
    119120   # add the left bound - the second breakpoint. but only if it is different than the lower bound and the mean.
    120    if (fabs(left-lb) > tolerance) and (fabs(left-xavg) > tolerance):     
     121   if (fabs(left-lb) > tolerance) and (fabs(left-xavg) > tolerance):
    121122      breakpoints.append(left)
    122123
     
    128129      if (fabs(this_lb-left) > tolerance) and (fabs(this_lb-right) > tolerance) and (fabs(this_lb-xavg) > tolerance):
    129130         breakpoints.append(this_lb)
    130       current_x += step           
     131      current_x += step
    131132
    132133   # add the right bound - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
    133    if (fabs(right-ub) > tolerance) and (fabs(right-xavg) > tolerance):           
    134       breakpoints.append(right)     
     134   if (fabs(right-ub) > tolerance) and (fabs(right-xavg) > tolerance):
     135      breakpoints.append(right)
    135136
    136137   # add the upper bound - the last breakpoint.
     
    146147#
    147148# routine to compute linearization breakpoints based on an exponential distribution from the mean in each direction.
    148 #   
     149#
    149150
    150151def compute_exponential_from_mean_breakpoints(lb, node_min, xavg, node_max, ub, num_breakpoints_per_side, tolerance):
     
    174175   for i in range(1,num_breakpoints_per_side+1):
    175176      current_x = xavg + current_offset
    176       if (fabs(current_x-xavg) > tolerance) and (fabs(current_x-ub) > tolerance):         
     177      if (fabs(current_x-xavg) > tolerance) and (fabs(current_x-ub) > tolerance):
    177178         breakpoints.append(current_x)
    178179      current_offset *= base
    179180
    180181   # add the upper bound - the last breakpoint.
    181    breakpoints.append(ub)         
    182 
    183    return breakpoints     
     182   breakpoints.append(ub)
     183
     184   return breakpoints
    184185
    185186#
     
    214215def form_standard_objective(instance_name, instance, original_objective_expression, scenario_tree):
    215216
    216    objective_name = instance.active_components(Objective).keys()[0]         
     217   objective_name = instance.active_components(Objective).keys()[0]
    217218   objective = instance.active_components(Objective)[objective_name]
    218219   # clone the objective, because we're going to augment (repeatedly) the original objective.
    219    objective_expression = original_objective_expression.clone() 
     220   objective_expression = original_objective_expression.clone()
    220221   objective_expression.simplify(instance)
    221222   instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression
     
    235236   new_instance_attributes = []
    236237
    237    objective_name = instance.active_components(Objective).keys()[0]         
     238   objective_name = instance.active_components(Objective).keys()[0]
    238239   objective = instance.active_components(Objective)[objective_name]
    239240   # clone the objective, because we're going to augment (repeatedly) the original objective.
    240    objective_expression = original_objective_expression.clone() 
     241   objective_expression = original_objective_expression.clone()
    241242   # the quadratic expression is really treated as just a list - eventually should be treated as a full expression.
    242243   quad_expression = 0.0
     
    244245   # for each blended variable (i.e., those not appearing in the final stage),
    245246   # add the linear and quadratic penalty terms to the objective.
    246    
     247
    247248   for stage in scenario_tree._stages[:-1]: # skip the last stage, as no blending occurs
    248249
     
    261262         w_parameter_name = "PHWEIGHT_"+variable_name
    262263         w_parameter = instance.active_components(Param)[w_parameter_name]
    263            
     264
    264265         average_parameter_name = "PHAVG_"+variable_name
    265266         average_parameter = instance.active_components(Param)[average_parameter_name]
     
    272273
    273274         node_min_parameter = variable_tree_node._minimums[variable_name]
    274          node_max_parameter = variable_tree_node._maximums[variable_name]               
     275         node_max_parameter = variable_tree_node._maximums[variable_name]
    275276
    276277         quad_penalty_term_variable = None
     
    285286            if (instance_variable[index].status is not VarStatus.unused) and (instance_variable[index].fixed is False):
    286287
    287                # add the linear (w-weighted) term is a consistent fashion, independent of variable type. 
     288               # add the linear (w-weighted) term is a consistent fashion, independent of variable type.
    288289               # if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes for reasons currently unknown.
    289290               objective_expression += (w_parameter[index] * instance_variable[index])
     
    303304                        new_term = (blend_parameter[index] * rho_parameter[index] / 2.0 * instance_variable[index]) - \
    304305                                   (blend_parameter[index] * rho_parameter[index] * average_parameter[index] * instance_variable[index]) + \
    305                                    (blend_parameter[index] * rho_parameter[index] / 2.0 * average_parameter[index] * average_parameter[index])                             
     306                                   (blend_parameter[index] * rho_parameter[index] / 2.0 * average_parameter[index] * average_parameter[index])
    306307                        if objective.sense is minimize:
    307308                           objective_expression += new_term
    308309                        else:
    309                            objective_expression -= new_term                                 
     310                           objective_expression -= new_term
    310311                     else:
    311312                        quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)
     
    321322                        else:
    322323                           objective_expression -= (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
    323                            
     324
    324325                        # the constraints are somewhat nastier.
    325326
     
    328329                        x = instance_variable[index]
    329330
    330                         lb = None
    331                         ub = None
    332 
    333                         if x.lb is None:
    334                            raise ValueError, "No lower bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
    335                         else:
    336                            lb = x.lb()
    337 
    338                         if x.ub is None:
    339                            raise ValueError, "No upper bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
    340                         else:
    341                            ub = x.ub()
     331                        if x.lb is None or x.ub is None:
     332                            var = variable_name + indexToString(index)
     333                            msg = "Missing bound for '%s'\n"                  \
     334                                  'Both lower and upper bounds required when' \
     335                                  ' piece-wise approximating quadratic '      \
     336                                  'penalty terms'
     337                            raise ValueError, msg % var
     338                        lb = x.lb()
     339                        ub = x.ub()
    342340
    343341                        node_min = node_min_parameter[index]()
     
    345343
    346344                        # compute the breakpoint sequence according to the specified strategy.
    347                         breakpoints = []
    348                         if breakpoint_strategy == 0:
    349                            breakpoints = compute_uniform_breakpoints(lb, node_min, xavg(), node_max, ub, \
    350                                                                      linearize_nonbinary_penalty_terms, \
    351                                                                      tolerance)
    352                         elif breakpoint_strategy == 1:
    353                            breakpoints = compute_uniform_between_nodestat_breakpoints(lb, node_min, xavg(), node_max, ub, \
    354                                                                                       linearize_nonbinary_penalty_terms, \
    355                                                                                       tolerance)
    356                         elif breakpoint_strategy == 2:
    357                            breakpoints = compute_uniform_between_woodruff_breakpoints(lb, node_min, xavg(), node_max, ub, \
    358                                                                                       linearize_nonbinary_penalty_terms, \
    359                                                                                       tolerance)
    360                         elif breakpoint_strategy == 3:
    361                            breakpoints = compute_exponential_from_mean_breakpoints(lb, node_min, xavg(), node_max, ub, \
    362                                                                                    linearize_nonbinary_penalty_terms, \
    363                                                                                    tolerance)
    364                         else:
    365                            raise ValueError, "A breakpoint distribution strategy="+str(breakpoint_strategy)+" is currently not supported within PH!"
     345                        try:
     346                            strategy = (
     347                              compute_uniform_breakpoints,
     348                              compute_uniform_between_nodestat_breakpoints,
     349                              compute_uniform_between_woodruff_breakpoints,
     350                              compute_exponential_from_mean_breakpoints,
     351                            )[ breakpoint_strategy ]
     352                            args = ( lb, node_min, xavg(), node_max, ub, \
     353                                linearize_nonbinary_penalty_terms, tolerance )
     354                            breakpoints = strategy( *args )
     355                        except ValueError, e:
     356                            msg = 'A breakpoint distribution strategy (%s) '  \
     357                                  'is currently not supported within PH!'
     358                            raise ValueError, msg % breakpoint_strategy
    366359
    367360                        for i in range(0,len(breakpoints)-1):
     
    381374                                                                                     tolerance)
    382375                           piece_constraint.add(None, piece_expression)
    383                            setattr(instance, piece_constraint_name, piece_constraint)                                   
     376                           setattr(instance, piece_constraint_name, piece_constraint)
    384377
    385378                     else:
    386379
    387                         quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)                           
    388              
     380                        quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)
     381
    389382   # strictly speaking, this probably isn't necessary - parameter coefficients won't get
    390383   # pre-processed out of the expression tree. however, if the under-the-hood should change,
  • coopr.pysp/stable/2.4/coopr/pysp/phserver.py

    r2410 r2938  
    3535
    3636# for profiling
    37 import cProfile
     37try:
     38    import cProfile as profile
     39except ImportError:
     40    import profile
    3841import pstats
    3942
     
    504507        #
    505508        tfile = pyutilib.services.TempfileManager.create_tempfile(suffix=".profile")
    506         tmp = cProfile.runctx('exec_ph(options)',globals(),locals(),tfile)
     509        tmp = profile.runctx('exec_ph(options)',globals(),locals(),tfile)
    507510        p = pstats.Stats(tfile).strip_dirs()
    508511        p.sort_stats('time', 'cum')
  • coopr.pysp/stable/2.4/coopr/pysp/phutils.py

    r2465 r2938  
    7777      return return_index[0]
    7878   else:
    79       return return_index     
     79      return return_index
    8080
    8181#
     
    104104
    105105   return_index = ()
    106    
     106
    107107   for index in indices:
    108108
     
    121121         transformed_index = int(index)
    122122      except ValueError:
    123          transformed_index = index                 
     123         transformed_index = index
    124124      return_index = return_index + (transformed_index,)
    125125
     
    140140   # ditto with the index template. one-dimensional
    141141   # indices in pyomo are not tuples, but anything
    142    # else is. 
     142   # else is.
    143143
    144144   if type(index) != tuple:
     
    192192      else:
    193193         print "Node-based instance initialization enabled"
    194    
     194
    195195   scenario = scenario_tree_instance.get_scenario(scenario_name)
    196196   scenario_instance = None
     
    217217         scenario_data.read(model=scenario_instance)
    218218         scenario_instance.load(scenario_data)
    219          scenario_instance.preprocess()
    220219   except:
    221220      print "Encountered exception in model instance creation - traceback:"
     
    234233def create_ph_parameters(instance, scenario_tree, default_rho, linearizing_penalty_terms):
    235234
    236    new_penalty_variable_names = []   
    237    
     235   new_penalty_variable_names = []
     236
    238237   # first, gather all unique variables referenced in any stage
    239238   # other than the last, independent of specific indices. this
     
    245244   # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
    246245   instance_variables = {}
    247    
     246
    248247   for stage in scenario_tree._stages[:-1]:
    249      
     248
    250249      for (reference_variable, index_template, reference_indices) in stage._variables:
    251            
     250
    252251         if reference_variable.name not in instance_variables.keys():
    253                
     252
    254253            instance_variables[reference_variable.name] = reference_variable
    255254
     
    283282      # PH AVG
    284283
    285       new_avg_index = reference_variable._index 
     284      new_avg_index = reference_variable._index
    286285      new_avg_parameter_name = "PHAVG_"+reference_variable.name
    287286      new_avg_parameter = None
     
    297296      # PH RHO
    298297
    299       new_rho_index = reference_variable._index 
     298      new_rho_index = reference_variable._index
    300299      new_rho_parameter_name = "PHRHO_"+reference_variable.name
    301300      new_rho_parameter = None
     
    318317         new_penalty_term_variable = None
    319318         if (len(new_avg_index) is 1) and (None in new_avg_index):
    320             new_penalty_term_variable = Var(name=new_penalty_term_variable_name, bounds=(0.0,None))                 
     319            new_penalty_term_variable = Var(name=new_penalty_term_variable_name, bounds=(0.0,None))
    321320         else:
    322321            new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None))
     
    325324         new_penalty_variable_names.append(new_penalty_term_variable_name)
    326325
    327       # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY. 
     326      # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY.
    328327
    329328      # also controls whether weight updates proceed at any iteration.
     
    341340         new_blend_parameter[index] = 1
    342341
    343    return new_penalty_variable_names     
     342   return new_penalty_variable_names
  • coopr.pysp/stable/2.4/coopr/pysp/scenariotree.py

    r2461 r2938  
    2323
    2424   #
    25    # initialize the _solutions attribute of a tree node. 
     25   # initialize the _solutions attribute of a tree node.
    2626   #
    2727
     
    3131      #       construct unless we're actually going to use!
    3232      # for each variable referenced in the stage, clone the variable
    33       # for purposes of storing solutions. we are being wasteful in 
     33      # for purposes of storing solutions. we are being wasteful in
    3434      # terms copying indices that may not be referenced in the stage.
    3535      # this is something that we might revisit if space/performance
     
    4040         # here are computed elsewhere - and that entity is responsible for
    4141         # ensuring feasibility. this also leaves room for specifying infeasible
    42          # or partial solutions.       
     42         # or partial solutions.
    4343         new_variable_index = variable._index
    4444         new_variable_name = variable.name
     
    6363
    6464      self._solutions_initialized = True
    65      
     65
    6666
    6767   """ Constructor
     
    8181      # statistic computed over all scenarios for that node. the
    8282      # parameters are named as the source variable name suffixed
    83       # by one of: "NODEMIN", "NODEAVG", and "NODEMAX". 
     83      # by one of: "NODEMIN", "NODEAVG", and "NODEMAX".
    8484      # NOTE: the averages are probability_weighted - the min/max
    8585      #       values are not.
     
    8787      #       convention is assumed to be enforced by whoever populates
    8888      #       these parameters.
    89       self._averages = {} 
     89      self._averages = {}
    9090      self._minimums = {}
    9191      self._maximums = {}
     
    111111
    112112      if self._solutions_initialized is False:
    113          self._initialize_solutions()     
     113         self._initialize_solutions()
    114114
    115115      for variable_name, variable in self._solutions.items():
     
    122122         except:
    123123            raise RuntimeError, "No averages parameter present on tree node="+self._name+" for variable="+variable_name
    124            
     124
    125125         for index in variable._index:
    126126            if variable[index].active is True:
     
    135135
    136136      if self._solutions_initialized is False:
    137          self._initialize_solutions()     
     137         self._initialize_solutions()
    138138
    139139      for variable_name, variable in self._solutions.items():
     
    142142               node_probability = 0.0
    143143               avg = 0.0
     144               num_scenarios_with_index = 0
    144145               for scenario in self._scenarios:
    145146                  scenario_instance = scenario_instance_map[scenario._name]
    146147                  node_probability += scenario._probability
    147                   var_value = getattr(scenario_instance, variable.name)[index].value
    148                   avg += (scenario._probability * var_value)
    149                variable[index].value = avg / node_probability
    150            
     148                  scenario_variable = getattr(scenario_instance, variable.name)
     149                  if index in scenario_variable:
     150                     num_scenarios_with_index = num_scenarios_with_index + 1
     151                     var_value = getattr(scenario_instance, variable.name)[index].value
     152                     avg += (scenario._probability * var_value)
     153               if num_scenarios_with_index > 0:
     154                  variable[index].value = avg / node_probability
     155
    151156   #
    152157   # a utility to compute the cost of the current node plus the expected costs of child nodes.
     
    224229               if variable_name not in self._reference_instance.active_components(Var):
    225230                  raise ValueError, "Variable=" + variable_name + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
    226                variable = self._reference_instance.active_components(Var)[variable_name]               
     231               variable = self._reference_instance.active_components(Var)[variable_name]
    227232
    228233               # extract all "real", i.e., fully specified, indices matching the index template.
    229                match_indices = extractVariableIndices(variable, index_template)               
     234               match_indices = extractVariableIndices(variable, index_template)
    230235
    231236               # there is a possibility that no indices match the input template.
     
    233238               if len(match_indices) == 0:
    234239                  raise RuntimeError, "No indices match template="+str(index_template)+" for variable="+variable_name+" ; encountered in scenario tree specification for model="+self._reference_instance.name
    235                  
     240
    236241               stage._variables.append((variable, index_template, match_indices))
    237242
     
    247252               # match template (e.g., "foo[*,*]") instead of just the variable
    248253               # name (e.g., "foo") to represent the set of all indices.
    249                
     254
    250255               # if the variable is a singleton - that is, non-indexed - no brackets is fine.
    251256               # we'll just tag the var[None] variable value with the (suffix,value) pair.
    252257               if None not in variable._index:
    253                   raise RuntimeError, "Variable="+variable_string+" is an indexed variable, and templates must specify an index match; encountered in scenario tree specification for model="+self._reference_instance.name                 
     258                  raise RuntimeError, "Variable="+variable_string+" is an indexed variable, and templates must specify an index match; encountered in scenario tree specification for model="+self._reference_instance.name
    254259
    255260               match_indices = []
     
    263268            raise ValueError, "Unknown stage=" + stage_id + " specified in scenario tree constructor (stage->cost variable map)"
    264269         stage = self._stage_map[stage_id]
    265          
     270
    266271         cost_variable_string = stage_cost_variable_ids[stage_id].value # de-reference is required to access the parameter value
    267272
     
    279284            if cost_variable_name not in self._reference_instance.active_components(Var):
    280285               raise ValueError, "Variable=" + cost_variable_name + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
    281             cost_variable = self._reference_instance.active_components(Var)[cost_variable_name]               
     286            cost_variable = self._reference_instance.active_components(Var)[cost_variable_name]
    282287
    283288            # extract all "real", i.e., fully specified, indices matching the index template.
     
    285290
    286291            # only one index can be supplied for a stage cost variable.
    287             if len(match_indices) != 1:
    288                raise RuntimeError, "Only one index can be specified for a stage cost variable - "+str(len(match_indices))+"match template="+index_template+" for variable="+cost_variable_name+" ; encountered in scenario tree specification for model="+self._reference_instance.name
     292            if len(match_indices) > 1:
     293                msg = 'Only one index can be specified for a stage cost '     \
     294                      'variable - %s match template "%s" for variable "%s" ;' \
     295                      ' encountered in scenario tree specification for model' \
     296                      ' "%s"'
     297                raise RuntimeError, msg % (
     298                  len(match_indices),
     299                  index_template,
     300                  cost_variable_name,
     301                  self._reference_instance.name
     302                )
     303            elif len(match_indices) == 0:
     304                msg = 'Stage cost index not found: %s[%s]\n'                  \
     305                      'Do you have an off-by-one miscalculation, or did you ' \
     306                      'forget to specify it in ReferenceModel.dat?'
     307                raise RuntimeError, msg % ( cost_variable_name, index_template )
    289308
    290309            cost_variable_index = match_indices[0]
     
    298317               raise ValueError, "Cost variable=" + cost_variable_name + " associated with stage=" + stage_id + " is not present in model=" + self._reference_instance.name
    299318            cost_variable = self._reference_instance.active_components(Var)[cost_variable_name]
    300            
     319
    301320         # store the validated info.
    302321         stage._cost_variable = (cost_variable, cost_variable_index)
     
    307326           scenariotreeinstance - the pyomo model specifying all scenario tree (text) data.
    308327           scenariobundlelist   - a list of scenario names to retain, i.e., cull the rest to create a reduced tree!
    309    """ 
     328   """
    310329   def __init__(self, *args, **kwds):
    311330      self._name = None # TBD - some arbitrary identifier
     
    340359            self._reference_instance = kwds[key]
    341360         elif key == "scenariotreeinstance":
    342             scenario_tree_instance = kwds[key]           
     361            scenario_tree_instance = kwds[key]
    343362         elif key == "scenariobundlelist":
    344             scenario_bundle_list = kwds[key]           
     363            scenario_bundle_list = kwds[key]
    345364         else:
    346365            print "Unknown option=" + key + " specified in call to ScenarioTree constructor"
     
    370389         raise ValueError, "An ordered set of stage IDs must be supplied in the ScenarioTree constructor"
    371390
    372       #     
     391      #
    373392      # construct the actual tree objects
    374393      #
     
    390409            raise ValueError, "Unknown stage=" + stage_name + " assigned to tree node=" + tree_node._name
    391410
    392          new_tree_node = ScenarioTreeNode(tree_node_name, 
    393                                           node_probability_map[tree_node_name].value, 
     411         new_tree_node = ScenarioTreeNode(tree_node_name,
     412                                          node_probability_map[tree_node_name].value,
    394413                                          self._stage_map[stage_name],
    395414                                          self._reference_instance)
     
    418437      # two-pass logic necessary when constructing scenarios.
    419438      for scenario_name in scenario_ids:
     439
     440         # IMPT: the name of the scenario is assumed to have no '_' (underscore) characters in the identifier.
     441         #       this is required when writing the extensive form, e.g., the scenario is used in the extensive
     442         #       form as a prefix on variable and constraint names. this convention simplifies parsing on the
     443         #       back end; if the underscore isn't used as a reserved character, then some other separator
     444         #       symbol would be required, or we would have to engage in some complex prefix matching with
     445         #       all possible scenario names.
     446         if string.find(scenario_name, "_") != -1:
     447            raise ValueError, "By convention, scenario names in PySP cannot contain underscore (_) characters; the scenario in violation="+scenario_name
     448         
    420449         new_scenario = Scenario()
    421450         new_scenario._name=scenario_name
     
    454483      # active scenario tree components and compress the tree.
    455484      if scenario_bundle_list is not None:
    456          print "Compressing scenario tree!"
    457485         self.compress(scenario_bundle_list)
    458486
     
    485513   #
    486514   # utility for compressing or culling a scenario tree based on
    487    # a provided list of scenarios (specified by name) to retain - 
     515   # a provided list of scenarios (specified by name) to retain -
    488516   # all non-referenced components are eliminated. this particular
    489517   # method compresses *in-place*, i.e., via direct modification
     
    502530         scenario.retain = True
    503531
    504          # chase all nodes comprising this scenario, 
     532         # chase all nodes comprising this scenario,
    505533         # marking them for retention.
    506534         for node in scenario._node_list:
     
    515543            pass
    516544         else:
    517             scenarios_to_delete.append(scenario)             
     545            scenarios_to_delete.append(scenario)
    518546            del self._scenario_map[scenario._name]
    519547
     
    631659            print "There are no scenarios associated with tree node=" + tree_node._name
    632660            return False
    633      
     661
    634662      return True
    635663
     
    653681      for tree_node in self._tree_nodes:
    654682
    655          tree_node.snapshotSolutionFromInstances(scenario_instance_map)         
     683         tree_node.snapshotSolutionFromInstances(scenario_instance_map)
    656684
    657685   #
     
    668696      else:
    669697         print "Model=" + "Unassigned"
    670       print "----------------------------------------------------"         
     698      print "----------------------------------------------------"
    671699      print "Tree Nodes:"
    672700      print ""
     
    675703         print "\tName=" + tree_node_name
    676704         if tree_node._stage is not None:
    677             print "\tStage=" + tree_node._stage._name
     705            print "\tStage=" + str(tree_node._stage._name)
    678706         else:
    679707            print "\t Stage=None"
     
    703731      for stage_name in sorted(self._stage_map.keys()):
    704732         stage = self._stage_map[stage_name]
    705          print "\tName=" + stage_name
     733         print "\tName=" + str(stage_name)
    706734         print "\tTree Nodes: "
    707735         for tree_node in sorted(stage._tree_nodes, cmp=lambda x,y: cmp(x._name, y._name)):
     
    713741            else:
    714742               print "\t\t",variable.name,":",index_template
    715          print "\tCost Variable: "           
     743         print "\tCost Variable: "
    716744         if stage._cost_variable[1] is None:
    717745            print "\t\t" + stage._cost_variable[0].name
    718746         else:
    719747            print "\t\t" + stage._cost_variable[0].name + indexToString(stage._cost_variable[1])
    720          print ""           
     748         print ""
    721749      print "----------------------------------------------------"
    722750      print "Scenarios:"
     
    732760         for tree_node in scenario._node_list:
    733761            print "\t\t" + tree_node._name
    734          print ""                       
     762         print ""
    735763      print "----------------------------------------------------"
    736764
     
    741769   def pprintSolution(self, epsilon=1.0e-5):
    742770
    743       print "----------------------------------------------------"         
     771      print "----------------------------------------------------"
    744772      print "Tree Nodes:"
    745773      print ""
     
    766794            else:
    767795               for index in indices:
    768                   if solution_variable[index].active is True:
     796                  if (solution_variable[index].active is True) and (index in solution_variable):
    769797                     value = solution_variable[index]()
    770                      if fabs(value) > epsilon:
     798                     if (value is not None) and (fabs(value) > epsilon):
    771799                        print "\t\t"+variable.name+indexToString(index)+"="+str(value)
    772          print ""           
     800         print ""
    773801
    774802   #
     
    787815         print "Model=" + "Unassigned"
    788816
    789       print "----------------------------------------------------"     
     817      print "----------------------------------------------------"
    790818      print "Tree Nodes:"
    791819      print ""
     
    845873         print "\tTotal scenario cost=%10.4f" % aggregate_cost
    846874         print ""
    847       print "----------------------------------------------------"     
     875      print "----------------------------------------------------"
  • coopr.pysp/stable/2.4/coopr/pysp/tests/unit/test_ph.py

    r2423 r2938  
    4747# Define a testing class, using the unittest.TestCase class.
    4848#
     49
     50# @unittest.category('nightly')
    4951class TestPH(unittest.TestCase):
    5052
     
    125127        self.cleanup()       
    126128
     129TestPH = unittest.category('nightly')(TestPH)
     130
     131
    127132if __name__ == "__main__":
    128133    unittest.main()
  • coopr.pysp/stable/2.4/doc/pysp/pyspbody.tex

    r2162 r2938  
    1616\section{Overview}
    1717
    18 The pysp package extends the pyomo modeling language to support multi-stage stochastic programs with
    19 enumerated scenarios. Pyomo and pysp are Python version 2.6 programs.  In order to
    20 specify a program, the user must provide a reference model and a scenario tree.
    21 
    22 Provided and the necessary paths have been communicated to the operating system, the command to execute the pysp package is
    23 of the form:
     18The pysp package extends the pyomo modeling language to support multi-stage
     19stochastic programs with enumerated scenarios. Pyomo and pysp are Python version
     202.6 programs. In order to specify a program, the user must provide a reference
     21model and a scenario tree.
     22
     23Provided and the necessary paths have been communicated to the operating system,
     24the command to execute the pysp package is of the form:
     25
    2426\begin{verbatim}
    2527runph
    2628\end{verbatim}
    27 It is possible,
    28 and generally necessary, to provide command line arguments. The simplest argument causes the program to output help text:
     29
     30It is possible, and generally necessary, to provide command line arguments. The
     31simplest argument causes the program to output help text:
     32
    2933\begin{verbatim}
    3034runph --help
    3135\end{verbatim}
    32 but notice that there are two dashes before the word ``help.'' Command line arguments are summarized in Section~\ref{cmdargsec}.
    33 
    34 The underlying algorithm in pysp
    35 is based on Progressive Hedging (PH) \cite{RockafellarWets}, which decomposes the problem into sub-problems, one for
    36 each scenario.
    37 The algorithm progressively computes {\em weights} corresponding to each variable to force convergence
    38 and also makes use of a {\em proximal} term that provides a penalty for the squared deviation from the
    39 mean solution from the last PH iteration.
     36
     37but notice that there are two dashes before the word ``help.'' Command line
     38arguments are summarized in Section~\ref{cmdargsec}.
     39
     40The underlying algorithm in pysp is based on Progressive Hedging (PH)
     41\cite{RockafellarWets}, which decomposes the problem into sub-problems, one for
     42each scenario. The algorithm progressively computes {\em weights} corresponding
     43to each variable to force convergence and also makes use of a {\em proximal}
     44term that provides a penalty for the squared deviation from the mean solution
     45from the last PH iteration.
    4046
    4147\subsection{Reference Model}
    4248
    43 The reference model describes the problem for a canonical scenario. It does not make use of, or describe, a scenario index or any information about uncertainty. Typically, it is just the model that would be used if there were only a single scenario. It is given as a pyomo file. Data from an arbitrary scenario is needed to instantiate.
    44 
    45 The objective function needs to be separated by stages. The term for each stage should be ``assigned'' (i.e., constrained to be equal to)
    46 a variable. These variables names are reported in ScenarioStructure.dat so that they can be used for reporting purposes.
     49The reference model describes the problem for a canonical scenario. It does not
     50make use of, or describe, a scenario index or any information about uncertainty.
     51Typically, it is just the model that would be used if there were only a single
     52scenario. It is given as a pyomo file. Data from an arbitrary scenario is needed
     53to instantiate.
     54
     55The objective function needs to be separated by stages. The term for each stage
     56should be ``assigned'' (i.e., constrained to be equal to) a variable. These
     57variable names are reported in ScenarioStructure.dat so that they can be used
     58for reporting purposes.
    4759
    4860\subsection{Scenario Tree}
    4961
    50 The scenario tree provides information about the time stages and the nature of the uncertainties. In order to specify a tree, we must
    51 indicate the time stages at which information becomes available. We also specify the nodes of a tree to indicate which variables
    52 are associated with which realization at each stage. The data for each scenario is provided in separate data files, one for each scenario.
     62The scenario tree provides information about the time stages and the nature of
     63the uncertainties. In order to specify a tree, we must indicate the time stages
     64at which information becomes available. We also specify the nodes of a tree to
     65indicate which variables are associated with which realization at each stage.
     66The data for each scenario is provided in separate data files, one for each
     67scenario.
    5368
    5469\section{File Structure}
     
    6176\end{itemize}
    6277
    63 In this list we use ``Sname'' as the generic scenario name. The file scenariostructure.dat gives the names of all the scenarios and for each scenario there is a data file with the same name and the suffix ``.dat'' that contains the full specification of data for the scenario.
     78In this list we use ``Sname'' as the generic scenario name. The file
     79\verb|ScenarioStructure.dat| gives the names of all the scenarios and for each
     80scenario there is a data file with the same name and the suffix ``.dat'' that
     81contains the full specification of data for the scenario.
    6482
    6583\subsection{ScenarioStructure.dat}
     
    6886
    6987\begin{itemize}
    70 \item set Scenarios: List of the names of the scenarios. These names will subsequently be used as indexes in this data file and these names will
    71 also be used as the root file names for the scenario data files (each of these will have a .dat extension) if the
    72 parameter ScenarioBasedData is set to True, which is the default.
    73 \item[]
    74 item set Stages: List of the names of the time stages, which must be given in time order. In the sequel we will use {\sc StageName} to represent a node name used
    75 as an index.
    76 \item[]
    77 \item set Nodes: List of the names of the nodes in the scenario tree. In the sequel we will use {\sc NodeName} to represent a node name used
    78 as an index.
    79 \item[]
    80 \item param NodeStage: A list of pairs of nodes and stages to indicate the stage for each node.
    81 \item[]
    82 \item param Parent: A list of node pairs to indicate the parent of each node that has a parent (the root node will not be listed).
    83 \item[]
    84 \item set Children[{\sc NodeName}]: For each node that has children, provide the list of children. No sets will be give for leaf nodes.
    85 \item[]
    86 \item param ConditionalProbability: For each node in the scenario tree, give the conditional probability. For the root node it must be given as 1 and for
    87 the children of any node with children, the conditional probabilities must sum to 1.
    88 \item[]
    89 \item param ScenarioLeafNode: A list of scenario and node pairs to indicate the leaf node for each scenario.
    90 \item[]
    91 \item set StageVariables[{\sc StageName}]: For each stage, list the pyomo model variables associated with that stage.
     88  \item set Scenarios: List of the names of the scenarios. These names will
     89  subsequently be used as indices in this data file and these names will also be
     90  used as the root file names for the scenario data files (each of these will
     91  have a .dat extension) if the parameter ScenarioBasedData is set to True,
     92  which is the default.
     93
     94  \item set Stages: List of the names of the time stages, which must be given in
     95  time order. In the sequel we will use {\sc StageName} to represent a node name
     96  used as an index.
     97
     98  \item set Nodes: List of the names of the nodes in the scenario tree. In the
     99  sequel we will use {\sc NodeName} to represent a node name used as an index.
     100
     101  \item param NodeStage: A list of pairs of nodes and stages to indicate the
     102  stage for each node.
     103
     104  \item param Parent: A list of node pairs to indicate the parent of each node
     105  that has a parent (the root node will not be listed).
     106
     107  \item set Children[{\sc NodeName}]: For each node that has children, provide
     108  the list of children. No sets will be give for leaf nodes.
     109
     110  \item param ConditionalProbability: For each node in the scenario tree, give
     111  the conditional probability. For the root node it must be given as 1 and for
     112  the children of any node with children, the conditional probabilities must sum
     113  to 1.
     114
     115  \item param ScenarioLeafNode: A list of scenario and node pairs to indicate
     116  the leaf node for each scenario.
     117
     118  \item set StageVariables[{\sc StageName}]: For each stage, list the pyomo
     119  model variables associated with that stage.
    92120\end{itemize}
    93121
    94 Data
    95 to instantiate these sets and parameters is provided by users in the file ScenarioStructure.dat, which can be given in AMPL \cite{ampl} format.
    96 
    97 The default behavior is one file per scenario and each file has the full data for the scenario. An
    98 alternative is to specify just the data that changes from the root node in one file per tree node.
    99 To select this option, add the following line to ScenarioStructure.dat:
     122Data to instantiate these sets and parameters is provided by users in the file
     123ScenarioStructure.dat, which can be given in AMPL \cite{ampl} format.
     124
     125The default behavior is one file per scenario and each file has the full data
     126for the scenario. An alternative is to specify just the data that changes from
     127the root node in one file per tree node. To select this option, add the
     128following line to ScenarioStructure.dat:
    100129
    101130\verb|param ScenarioBasedData := False ;|
    102131
    103 This will set it up to want a per-node file, something along the lines of what's in \verb|examples/pysp/farmer/NODEDATA|.
    104 
    105 Advanced users may be interested in seeing the file \verb|coopr/pysp/utils/scenariomodels.py|, which defines the python sets and parameters needed to describe stochastic elements. This file should not be edited.
     132This will set it up to want a per-node file, something along the lines of what's
     133in \verb|examples/pysp/farmer/NODEDATA|.
     134
     135Advanced users may be interested in seeing the file
     136\verb|coopr/pysp/utils/scenariomodels.py|, which defines the python sets and
     137parameters needed to describe stochastic elements. This file should not be
     138edited.
    106139
    107140\section{Command Line Arguments \label{cmdargsec}}
    108141
    109 The basic PH algorithm is controlled by parameters that are set as command line arguments. Note that options begin with a double dash.
     142The basic PH algorithm is controlled by parameters that are set as command line
     143arguments. Note that options begin with a double dash.
    110144\begin{itemize}
    111   \item \verb|-h|, \verb|--help|\\            Show help message and exit.
    112   \item \verb|--verbose|\\             Generate verbose output for both initialization and execution. Default is False.
    113   \item \verb|--report-solutions|\\     Always report PH solutions after each iteration. Enabled if --verbose is enabled. Default is False.
    114   \item \verb|--report-weights|\\    Always report PH weights prior to each iteration. Enabled if --verbose is enabled. Default is False.
    115   \item \verb|--model-directory|=MODEL\_DIRECTORY\\
    116                         The directory in which all model (reference and scenario) definitions are stored. I.e., the ``.py'' files. Default is ".".
    117   \item \verb|--instance-directory|=INSTANCE\_DIRECTORY\\
    118                         The directory in which all instance (reference and scenario) definitions are stored. I.e., the ``.dat'' files. Default is ".".
    119   \item \verb|--solver|=SOLVER\_TYPE\\  The type of solver used to solve scenario sub-problems. Default is cplex.
    120   \item \verb|--solver-manager|=SOLVER\_MANAGER\_TYPE\\
    121                         The type of solver manager used to coordinate scenario sub-problem solves. Default is serial. This option is changed in parallel applications
    122 as described in Section~\ref{parallelsec}.
    123   \item \verb|--max-iterations|=MAX\_ITERATIONS\\
    124                         The maximal number of PH iterations. Default is 100.
    125   \item \verb|--default-rho|=DEFAULT\_RHO\\
    126                         The default (global) rho for all blended variables. Default is 1.
    127   \item \verb|--rho-cfgfile|=RHO\_CFGFILE\\
    128                         The name of a configuration script to compute PH rho values. Default is None.
    129   \item \verb|--enable-termdiff|-convergence\\
    130                         Terminate PH based on the termdiff convergence metric. The convergcne metric is the unscaled sum of differences between variable values and the mean. Default is True.
    131   \item \verb|--enable-normalized|-termdiff-convergence\\
    132                         Terminate PH based on the normalized termdiff
    133                         convergence metric. Each term in the termdiff sum is normalized by the average value (NOTE: it is NOT normalized by the number of scenarios). Default is False.
    134   \item \verb|--termdiff-threshold|=TERMDIFF\_THRESHOLD\\
    135                         The convergence threshold used in the term-diff and
    136                         normalized term-diff convergence criteria. Default is
    137                         0.01, which is too low for most problems.
    138   \item \verb|--enable-free-discrete-count-convergence|\\
    139                         Terminate PH based on the free discrete variable count convergence metric. Default is False.
    140   \item \verb|--free-discrete-count-threshold|=FREE\_DISCRETE\_COUNT\_THRESHOLD\\
    141                         The convergence threshold used in the criterion based on when the free discrete variable count convergence criterion. Default is 20.
    142   \item \verb|--enable-ww-extensions|\\
    143                         Enable the Watson-Woodruff PH extensions plugin. Default is False.
    144   \item \verb|--ww-extension-cfgfile|=WW\_EXTENSION\_CFGFILE\\
    145                         The name of a configuration file for the Watson-Woodruff PH extensions plugin. Default is wwph.cfg.
    146   \item \verb|--ww-extension-suffixfile|=WW\_EXTENSION\_SUFFIXFILE\\
    147                         The name of a variable suffix file for the Watson-Woodruff PH extensions plugin. Default is wwph.suffixes.
    148    \item \verb|--user-defined-extension|=EXTENSIONFILE. Here, "EXTENSIONFILE" is the module name, which is in either the current directory (most likely) or somewhere on your PYTHONPATH. A simple example is "testphextension" plugin that simply prints a message to the screen for each callback. The file testphextension.py can be found in the sources directory and is shown in Section~\ref{ExtensionDetailsSec}. A test of this would be to specify "-user-defined-extension=testphextension", assuming testphextension.py is in your PYTHONPATH or current directory.
    149  
    150 Note that both PH extensions (WW PH and your own) can co-exist; however, the WW plugin will be invoked first.
    151  
    152 \item \verb|--scenario-solver-options| The options are specified just as in pyomo, e.g., \verb|--scenario-solver-options="mip_tolerances_mipgap=0.2"| to set the mipgap for all scenario sub-problem solves to 20\% for the CPLEX solver. The options are specified in a quote deliminted string that is passed to the sub-problem solver.
    153  Whatever options specified are persistent across all solves. 
    154 \item \verb|--ef-solver-options| The options are specified just as in pyomo, e.g., \verb|--scenario-solver-options="mip_tolerances_mipgap=0.2"| to set the mipgap for all scenario sub-problem solves to 20\% for the CPLEX solver. The options are specified in a quote deliminted string that is passed to the EF problem solver.
    155   \item \verb|--write-ef|\\            Upon termination, write the extensive form of the model - accounting for all fixed variables.
    156   \item \verb|--solve-ef|\\            Following write of the extensive form model, solve it.
    157   \item \verb|--ef-output-file|=EF\_OUTPUT\_FILE\\
    158                         The name of the extensive form output file (currently only LP format is supported), if writing of the extensive form is enabled. Default is efout.lp.
    159   \item \verb|--suppress-continuous-variable-output|\\
    160                         Eliminate PH-related output involving continuous variables. Default: no output.
    161   \item \verb|--keep-solver-files|\\   Retain temporary input and output files for scenario sub-problem solves.  Default: files not kept.
    162   \item \verb|--output-solver-logs|\\  Output solver logs during scenario sub-problem solves. Default: no output.
    163   \item \verb|--output-ef-solver-log|\\
    164                         Output solver log during the extensive form solve. Default: no output.
    165   \item \verb|--output-solver-results|\\
    166                         Output solutions obtained after each scenario sub-problem solve. Default: no output.
    167   \item \verb|--output-times|\\        Output timing statistics for various PH components. Default: no output.
    168   \item \verb|--disable-warmstarts|\\
    169                   Disable warm-start of scenario sub-problem solves in PH iterations >= 1.
    170                   Default=False (i.e., warm starts are the default).
    171   \item \verb|--drop-proximal-terms|\\
    172                   Eliminate proximal terms (i.e., the quadratic penalty terms) from the weighted PH objective.
    173                   Default=False (i.e., but default, the proximal terms are included).
    174   \item \verb|--retain-quadratic-binary-terms|\\
    175                   Do not linearize PH objective terms involving binary decision variables.
    176                   Default=False (i.e., the proximal term for binary variables is linearized by default; this can have some impact on the relaxations during the branch and bound solution process).
    177   \item \verb|--linearize-nonbinary-penalty-terms|=BPTS\\
    178                   Approximate the PH quadratic term for non-binary variables with a piece-wise linear function. The argument
    179                   BPTS gives the number of breakpoints in the linear approximation. The default=0. Reasonable non-zero
    180                   values are usually in the range of 3 to 7. Note that if a breakpoint would be very close to a variable
    181                   bound, then the break point is ommited. IMPORTANT: this option requires that all variables have bounds that
    182                   are either established in the reference model or by code specfied using the bounds-cfgfile command line option. See Section~\ref{LinearSec}
    183                   for more information about linearizing the proximal term.
    184   \item \verb|--breakpoint-strategy|=BREAKPOINT\_STRATEGY
    185                         Specify the strategy to distribute breakpoints on the
    186                         [lb, ub] interval of each variable when linearizing. 0
    187                         indicates uniform distribution. 1 indicates
    188                         breakpoints at the node min and max, uniformly in-
    189                         between. 2 indicates more aggressive concentration of
    190                         breakpoints near the observed node min/max.
    191   \item \verb|--bounds-cfgfile|=BOUNDS\_CFGFILE\\
    192                   The argument BOUNDS\_CFGFILE specifies the name of an executable pyomo file that sets bounds. The devault is that there is no file.
    193                   When specified, the code in this file is executed after the initialization of scenario data so the
    194                   bounds can be based on data from all scenarios. The config subdirectory of the farmer example contains a simple example of
    195                   such a file (boundsetter.cfg).
    196   \item \verb|--checkpoint-interval|\\
    197                   The number of iterations between writing of a checkpoint file. Default is 0, indicating never.
    198   \item \verb|--restore-from-checkpoint|\\
    199                   The name of the checkpoint file from which PH should be initialized. Default is not to restore from a checkpoint.
    200   \item \verb|--profile=PROFILE|\\     Enable profiling of Python code.  The value of this option is the number of functions that are summarized. The default is no profiling.
    201   \item \verb|--enable-gc|\\           Enable the python garbage collecter. The default is no garbage collection.
     145  \item \verb|-h|, \verb|--help|                                             \\
     146  Show help message and exit.
     147
     148  \item \verb|--verbose|                                                     \\
     149  Generate verbose output for both initialization and execution. Default is
     150  False.
     151
     152  \item \verb|--report-solutions|                                            \\
     153  Always report PH solutions after each iteration. Enabled if --verbose is
     154  enabled. Default is False.
     155
     156  \item \verb|--report-weights|                                              \\
     157  Always report PH weights prior to each iteration. Enabled if --verbose is
     158  enabled. Default is False.
     159
     160  \item \verb|--model-directory|=MODEL\_DIRECTORY                            \\
     161  The directory in which all model (reference and scenario) definitions are
     162  stored. I.e., the ``.py'' files. Default is ".".
     163
     164  \item \verb|--instance-directory|=INSTANCE\_DIRECTORY                      \\
     165  The directory in which all instance (reference and scenario) definitions are
     166  stored. I.e., the ``.dat'' files. Default is ".".
     167
     168  \item \verb|--solver|=SOLVER\_TYPE                                         \\
     169  The type of solver used to solve scenario sub-problems. Default is cplex.
     170
     171  \item \verb|--solver-manager|=SOLVER\_MANAGER\_TYPE                        \\
     172  The type of solver manager used to coordinate scenario sub-problem solves.
     173  Default is serial. This option is changed in parallel applications as
     174  described in Section~\ref{parallelsec}.
     175
     176  \item \verb|--max-iterations|=MAX\_ITERATIONS                              \\
     177  The maximal number of PH iterations. Default is 100.
     178
     179  \item \verb|--default-rho|=DEFAULT\_RHO                                    \\
     180  The default (global) rho for all blended variables. Default is 1.
     181
     182  \item \verb|--rho-cfgfile|=RHO\_CFGFILE                                    \\
     183  The name of a configuration script to compute PH rho values. Default is None.
     184
     185  \item \verb|--enable-termdiff|-convergence                                 \\
     186  Terminate PH based on the termdiff convergence metric. The convergcne metric
     187  is the unscaled sum of differences between variable values and the mean.
     188  Default is True.
     189
     190  \item \verb|--enable-normalized|-termdiff-convergence                      \\
     191  Terminate PH based on the normalized termdiff convergence metric. Each term in
     192  the termdiff sum is normalized by the average value (NOTE: it is NOT
     193  normalized by the number of scenarios). Default is False.
     194
     195  \item \verb|--termdiff-threshold|=TERMDIFF\_THRESHOLD                      \\
     196  The convergence threshold used in the term-diff and normalized term-diff
     197  convergence criteria. Default is 0.01, which is too low for most problems.
     198
     199  \item \verb|--enable-free-discrete-count-convergence|                      \\
     200  Terminate PH based on the free discrete variable count convergence metric.
     201  Default is False.
     202
     203  \item \verb|--free-discrete-count-threshold|=FREE\_DISCRETE\_COUNT\_THRESHOLD \\
     204  The convergence threshold used in the criterion based on when the free
     205  discrete variable count convergence criterion. Default is 20.
     206
     207  \item \verb|--enable-ww-extensions|                                                                                                                                          \\
     208  Enable the Watson-Woodruff PH extensions plugin. Default is False.
     209
     210  \item \verb|--ww-extension-cfgfile|=WW\_EXTENSION\_CFGFILE                 \\
     211  The name of a configuration file for the Watson-Woodruff PH extensions plugin.
     212  Default is wwph.cfg.
     213
     214  \item \verb|--ww-extension-suffixfile|=WW\_EXTENSION\_SUFFIXFILE           \\
     215  The name of a variable suffix file for the Watson-Woodruff PH extensions
     216  plugin. Default is wwph.suffixes.
     217
     218  \item \verb|--user-defined-extension|=EXTENSIONFILE                        \\
     219  Here, "EXTENSIONFILE" is the module name, which is in either the current
     220  directory (most likely) or somewhere on your PYTHONPATH. A simple example is
     221  "testphextension" plugin that simply prints a message to the screen for each
     222  callback. The file testphextension.py can be found in the sources directory
     223  and is shown in Section~\ref{ExtensionDetailsSec}. A test of this would be to
     224  specify "-user-defined-extension=testphextension", assuming testphextension.py
     225  is in your PYTHONPATH or current directory. Note that both PH extensions (WW
     226  PH and your own) can co-exist; however, the WW plugin will be invoked first.
     227
     228  \item \verb|--scenario-solver-options|                                     \\
     229  The options are specified just as in pyomo, e.g.,
     230  \verb|--scenario-solver-options="mip_tolerances_mipgap=0.2"| to set the mipgap
     231  for all scenario sub-problem solves to 20\% for the CPLEX solver. The options
     232  are specified in a quote deliminted string that is passed to the sub-problem
     233  solver. Whatever options specified are persistent across all solves.
     234
     235  \item \verb|--ef-solver-options|                                           \\
     236  The options are specified just as in pyomo, e.g.,
     237  \verb|--scenario-solver-options="mip_tolerances_mipgap=0.2"| to set the mipgap
     238  for all scenario sub-problem solves to 20\% for the CPLEX solver. The options
     239  are specified in a quote deliminted string that is passed to the EF problem
     240  solver.
     241
     242  \item \verb|--write-ef|                                                    \\
     243  Upon termination, write the extensive form of the model - accounting for all
     244  fixed variables.
     245
     246  \item \verb|--solve-ef|                                                    \\
     247  Following write of the extensive form model, solve it.
     248
     249  \item \verb|--ef-output-file|=EF\_OUTPUT\_FILE                             \\
     250  The name of the extensive form output file (currently only LP format is
     251  supported), if writing of the extensive form is enabled. Default is efout.lp.
     252
     253  \item \verb|--suppress-continuous-variable-output|                         \\
     254  Eliminate PH-related output involving continuous variables. Default: no
     255  output.
     256
     257  \item \verb|--keep-solver-files|                                           \\
     258  Retain temporary input and output files for scenario sub-problem solves.
     259  Default: files not kept.
     260
     261  \item \verb|--output-solver-logs|                                          \\
     262  Output solver logs during scenario sub-problem solves. Default: no output.
     263
     264  \item \verb|--output-ef-solver-log|                                        \\
     265  Output solver log during the extensive form solve. Default: no output.
     266
     267  \item \verb|--output-solver-results|                                       \\
     268  Output solutions obtained after each scenario sub-problem solve. Default: no
     269  output.
     270
     271  \item \verb|--output-times|                                                \\
     272  Output timing statistics for various PH components. Default: no output.
     273
     274  \item \verb|--disable-warmstarts|                                          \\
     275  Disable warm-start of scenario sub-problem solves in PH iterations >= 1.
     276  Default=False (i.e., warm starts are the default).
     277
     278  \item \verb|--drop-proximal-terms|                                         \\
     279  Eliminate proximal terms (i.e., the quadratic penalty terms) from the weighted
     280  PH objective. Default=False (i.e., but default, the proximal terms are
     281  included).
     282
     283  \item \verb|--retain-quadratic-binary-terms|                               \\
     284  Do not linearize PH objective terms involving binary decision variables.
     285  Default=False (i.e., the proximal term for binary variables is linearized by
     286  default; this can have some impact on the relaxations during the branch and
     287  bound solution process).
     288
     289  \item \verb|--linearize-nonbinary-penalty-terms|=BPTS                      \\
     290  Approximate the PH quadratic term for non-binary variables with a piece-wise
     291  linear function. The argument BPTS gives the number of breakpoints in the
     292  linear approximation. The default=0. Reasonable non-zero values are usually in
     293  the range of 3 to 7. Note that if a breakpoint would be very close to a
     294  variable bound, then the break point is ommited. IMPORTANT: this option
     295  requires that all variables have bounds that are either established in the
     296  reference model or by code specfied using the bounds-cfgfile command line
     297  option. See Section~\ref{LinearSec} for more information about linearizing the
     298  proximal term.
     299
     300  \item \verb|--breakpoint-strategy|=BREAKPOINT\_STRATEGY                    \\
     301  Specify the strategy to distribute breakpoints on the [lb, ub] interval of
     302  each variable when linearizing. 0 indicates uniform distribution. 1 indicates
     303  breakpoints at the node min and max, uniformly in- between. 2 indicates more
     304  aggressive concentration of breakpoints near the observed node min/max.
     305
     306  \item \verb|--bounds-cfgfile|=BOUNDS\_CFGFILE                              \\
     307  The argument BOUNDS\_CFGFILE specifies the name of an executable pyomo file
     308  that sets bounds. The devault is that there is no file. When specified, the
     309  code in this file is executed after the initialization of scenario data so the
     310  bounds can be based on data from all scenarios. The config subdirectory of the
     311  farmer example contains a simple example of such a file (boundsetter.cfg).
     312
     313  \item \verb|--checkpoint-interval|                                         \\
     314  The number of iterations between writing of a checkpoint file. Default is 0,
     315  indicating never.
     316
     317  \item \verb|--restore-from-checkpoint|                                     \\
     318  The name of the checkpoint file from which PH should be initialized. Default
     319  is not to restore from a checkpoint.
     320
     321  \item \verb|--profile=PROFILE|                                             \\
     322  Enable profiling of Python code. The value of this option is the number of
     323  functions that are summarized. The default is no profiling.
     324
     325  \item \verb|--enable-gc|                                                   \\
     326  Enable the python garbage collecter. The default is no garbage collection.
    202327\end{itemize}
    203328
    204329\section{Extensions via Callbacks \label{CallbackSec}}
    205330
    206 Basic PH can converge slowly, so it is usually advisable to extend it or modify it. In pysp, this is done via the pyomo plug-in
    207 mechanism. The basic PH implementation provides callbacks that enable access to the data structures used by the algorithm.
    208 In \S\ref{WWExtensionSec} we describe extensions that are provided with the release. In \S\ref{ExtensionDetailsSec},
    209 we provide information to power users who may wish to modify or replace the extensions.
     331Basic PH can converge slowly, so it is usually advisable to extend it or modify
     332it. In pysp, this is done via the pyomo plug-in mechanism. The basic PH
     333implementation provides callbacks that enable access to the data structures used
     334by the algorithm. In \S\ref{WWExtensionSec} we describe extensions that are
     335provided with the release. In \S\ref{ExtensionDetailsSec}, we provide
     336information to power users who may wish to modify or replace the extensions.
    210337
    211338\subsection{Watson and Woodruff Extensions \label{WWExtensionSec}}
    212339
    213 Watson and Woodruff describe innovations for accelerating PH \cite{phinnovate}, most of which are generalized and
    214 implemented in the
    215 file \verb|wwextension.py|, but users generally do not need to know this file name. To invoke the program
     340Watson and Woodruff describe innovations for accelerating PH \cite{phinnovate},
     341most of which are generalized and implemented in the file \verb|wwextension.py|,
     342but users generally do not need to know this file name. To invoke the program
    216343with these additional features, invoke the software with a command of the form:
     344
    217345\begin{verbatim}
    218346runph --enable-ww-extensions
    219347\end{verbatim}
    220 Many of the examples described in \S\ref{ExampleSec} use this plug-in. The main concept
    221 is that some integer variables should be fixed as the algorithm progresses for two reasons:
     348
     349Many of the examples described in \S\ref{ExampleSec} use this plug-in. The main
     350concept is that some integer variables should be fixed as the algorithm
     351progresses for two reasons:
    222352
    223353\begin{itemize}
    224 \item Convergence detection:
    225 A detailed analysis of PH algorithm behavior on a variety of problem indicates that individual decision variables
    226 frequently converge to specific,
    227 fixed values scenarios in early PH iterations. Further, despite interactions among the
    228 the variables, the value frequently does not change in subsequent PH
    229 iterations. Such variable ``fixing'' behaviors lead to a potentially powerful, albeit obvious, heuristic: once
    230 a particular variable has been the same in all scenarios for some number of iterations, fix it to that value.
    231 For problems where the constraints effectively limit $x$
    232 from both sides, these methods may result in PH encountering infeasible scenario sub-problems even though the
    233 problem is ultimately feasible.
    234 \item Cycle detection: When there are integer variables, cycling is sometimes
    235 encountered, consequently, cycle detection and avoidance mechanisms are required to force
    236 eventual convergence of the PH algorithm in the mixed-integer case. To detect cycles, we focus on
    237 repeated occurrences of the weights, implemented using a simple hashing scheme \cite{tabuhash} to
    238 minimize impact on run-time. Once a cycle in the weight vectors associated with any decision variable is detected,
    239 the value of that variable is fixed.
     354  \item Convergence detection: A detailed analysis of PH algorithm behavior on a
     355  variety of problem indicates that individual decision variables frequently
     356  converge to specific, fixed values scenarios in early PH iterations. Further,
     357  despite interactions among the the variables, the value frequently does not
     358  change in subsequent PH iterations. Such variable ``fixing'' behaviors lead to
     359  a potentially powerful, albeit obvious, heuristic: once a particular variable
     360  has been the same in all scenarios for some number of iterations, fix it to
     361  that value. For problems where the constraints effectively limit $x$ from both
     362  sides, these methods may result in PH encountering infeasible scenario
     363  sub-problems even though the problem is ultimately feasible.
     364
     365  \item Cycle detection: When there are integer variables, cycling is sometimes
     366  encountered, consequently, cycle detection and avoidance mechanisms are
     367  required to force eventual convergence of the PH algorithm in the
     368  mixed-integer case. To detect cycles, we focus on repeated occurrences of the
     369  weights, implemented using a simple hashing scheme \cite{tabuhash} to minimize
     370  impact on run-time. Once a cycle in the weight vectors associated with any
     371  decision variable is detected, the value of that variable is fixed.
    240372\end{itemize}
    241373
    242 Fixing variables aggressively can result in shorter solution times, but can also result in solutions that are not as good. Furthermore, for some problems, aggressive fixing can result in infeasible sub-problems even though the problem is
    243 ultimately feasible. Many of the parameters discussed in the next subsections control fixing of variables. This is
    244 discussed in a tutorial in section~\ref{WWTutorialSec}.
     374Fixing variables aggressively can result in shorter solution times, but can also
     375result in solutions that are not as good. Furthermore, for some problems,
     376aggressive fixing can result in infeasible sub-problems even though the problem
     377is ultimately feasible. Many of the parameters discussed in the next subsections
     378control fixing of variables. This is discussed in a tutorial in
     379section~\ref{WWTutorialSec}.
    245380
    246381\subsubsection{Variable Specific Parameters}
    247382
    248 The plug-in makes use of parameters to control behavior at the variable level. Global defaults (to override the defaults stated here) should be set using methods described in \S\ref{ParmSec}. Values for each variable should be set using methods described in \S\ref{SuffixSec}.
    249 Note that for variable fixing based on convergence detection, iteration zero is treated separately. The parameters
    250 are as follows:
     383The plug-in makes use of parameters to control behavior at the variable level.
     384Global defaults (to override the defaults stated here) should be set using
     385methods described in \S\ref{ParmSec}. Values for each variable should be set
     386using methods described in \S\ref{SuffixSec}. Note that for variable fixing
     387based on convergence detection, iteration zero is treated separately. The
     388parameters are as follows:
    251389
    252390\begin{itemize}
    253 \item fix\_continuous\_variables: True or False. If true, fixing applies to all variables. If false, then fixing applies only to discrete variables.
    254 \item Iter0FixIfConvergedAtLB: 1 (True) or 0 (False). If 1, then discrete variables that are at their lower bound in all scenarios after
    255 the iteration zero solves will be fixed at that bound.
    256 \item Iter0FixIfConvergedAtUB: 1 (True) or 0 (False). If 1, then discrete variables that are at their upper bound in all sce<narios after
    257 the iteration zero solves will be fixed at that bound.
    258 \item Iter0FixIfConvergedAtNB: = 1  1 (True) or 0 (False). If 1, then discrete variables that are at the same value in all scenarios after
    259 the iteration zero solves will be fixed at that value, without regard to whether it is a bound. If this is true, it takes precedence. A value of zero, on the other hand, implies that variables will not be fixed at at a non-bound.
    260 \item FixWhenItersConvergedAtLB:  The number of consecutive PH iterations that discrete variables must be their lower bound in all scenarios before they
    261 will be fixed at that bound. A value of zero implies that variables will not be fixed at the bound.
    262 \item FixWhenItersConvergedAtUB: The number of consecutive PH iterations that discrete variables must be their upper bound in all scenarios before they
    263 will be fixed at that bound. A value of zero implies that variables will not be fixed at the bound.
    264 \item FixWhenItersConvergedAtNB: The number of consecutive PH iterations that discrete variables must be at the same, consistent value in all scenarios before they
    265 will be fixed at that value, without regard to whether it is a bound. If this is true, it takes precedence. A value of zero, on the other hand, implies that variables will not be fixed at at a non-bound.
    266 \item FixWhenItersConvergedContinuous: The number of consecutive PH iterations that continuous variables must be at the same, consistent value in all scenarios before they
    267 will be fixed at that value. A value of zero implies that continuous variables will not be fixed.
    268 \item CanSlamToLB: True or False. If True, then slamming can be to the lower bound for any variable.
    269 \item CanSlamToMin: True or False. If True, then slamming can be to the minimum across scenarios for any variable.
    270 \item CanSlamToAnywhere: True or False. If True, then slamming can be to any value.
    271 \item CanSlamToMax: True or False. If True, then slamming can be to the maximum across scenarios for any variable.
    272 \item CanSlamToUB: True of False. If True, then slamming can be to the upper bound for any variable.
    273 \item DisableCycleDetection: True or False. If True, then cycle detection and the associated slamming are completely disabled. This cannot be changed to False on the fly because a value of True at startup causes creation of the cycle detection storage to be bypassed.
     391  \item fix\_continuous\_variables: True or False. If true, fixing applies to
     392  all variables. If false, then fixing applies only to discrete variables.
     393
     394  \item Iter0FixIfConvergedAtLB: 1 (True) or 0 (False). If 1, then discrete
     395  variables that are at their lower bound in all scenarios after the iteration
     396  zero solves will be fixed at that bound.
     397
     398  \item Iter0FixIfConvergedAtUB: 1 (True) or 0 (False). If 1, then discrete
     399  variables that are at their upper bound in all sce<narios after the iteration
     400  zero solves will be fixed at that bound.
     401
     402  \item Iter0FixIfConvergedAtNB: = 1 1 (True) or 0 (False). If 1, then discrete
     403  variables that are at the same value in all scenarios after the iteration zero
     404  solves will be fixed at that value, without regard to whether it is a bound.
     405  If this is true, it takes precedence. A value of zero, on the other hand,
     406  implies that variables will not be fixed at at a non-bound.
     407
     408  \item FixWhenItersConvergedAtLB: The number of consecutive PH iterations that
     409  discrete variables must be their lower bound in all scenarios before they will
     410  be fixed at that bound. A value of zero implies that variables will not be
     411  fixed at the bound.
     412
     413  \item FixWhenItersConvergedAtUB: The number of consecutive PH iterations that
     414  discrete variables must be their upper bound in all scenarios before they will
     415  be fixed at that bound. A value of zero implies that variables will not be
     416  fixed at the bound.
     417
     418  \item FixWhenItersConvergedAtNB: The number of consecutive PH iterations that
     419  discrete variables must be at the same, consistent value in all scenarios
     420  before they will be fixed at that value, without regard to whether it is a
     421  bound. If this is true, it takes precedence. A value of zero, on the other
     422  hand, implies that variables will not be fixed at at a non-bound.
     423
     424  \item FixWhenItersConvergedContinuous: The number of consecutive PH iterations
     425  that continuous variables must be at the same, consistent value in all
     426  scenarios before they will be fixed at that value. A value of zero implies
     427  that continuous variables will not be fixed.
     428
     429  \item CanSlamToLB: True or False. If True, then slamming can be to the lower
     430  bound for any variable.
     431
     432  \item CanSlamToMin: True or False. If True, then slamming can be to the
     433  minimum across scenarios for any variable.
     434
     435  \item CanSlamToAnywhere: True or False. If True, then slamming can be to any
     436  value.
     437
     438  \item CanSlamToMax: True or False. If True, then slamming can be to the
     439  maximum across scenarios for any variable.
     440
     441  \item CanSlamToUB: True of False. If True, then slamming can be to the upper
     442  bound for any variable.
     443
     444  \item DisableCycleDetection: True or False. If True, then cycle detection and
     445  the associated slamming are completely disabled. This cannot be changed to
     446  False on the fly because a value of True at startup causes creation of the
     447  cycle detection storage to be bypassed.
    274448\end{itemize}
    275449
    276 In the event that multiple slam targets are True, then Min and Max trump LB and UB while Anywhere trumps all.
     450In the event that multiple slam targets are True, then Min and Max trump LB and
     451UB while Anywhere trumps all.
    277452
    278453\subsection{General Parameters}
    279454
    280 The plug-in also makes use of the following parameters, which should be set using methods described in \S\ref{ParmSec}.
     455The plug-in also makes use of the following parameters, which should be set
     456using methods described in \S\ref{ParmSec}.
     457
    281458\begin{itemize}
    282 \item Iteration0Mipgap: Gives the mipgap to be sent to the solver for iteration zero solves. The default is zero, which causes nothing to be sent to the solver (i.e., the solver uses its default mipgap).
    283 \item InitalMipGap: Gives the mipgap to be sent to the solver for iteration one solves. The default is zero, which causes nothing to be sent to the solver (i.e., the solver uses its default mipgap). If not zero, then this gap will be used as the starting point for the mipgap to change on each PH iteration in proportion to the convergence termination criterion so that when the criterion for termination is met
    284 the mipgap will be at (or near) the parameter value of FinalMipGap. If the InitialMipGap is significantly higher than the Iteration0MipGap parameter, the PH algorithm may
    285 perform poorly. This is because the iteration k-1 solutions are used to warm start iteration k-1 solves. If the iteration 1 mipgap is much higher than the iteration 0 mipgap, the iteration zero
    286 solution, although not optimal for the iteration one objective, might be within the mipgap. Default: 0.10.
    287 \item FinalMipGap: The target for the mipgap when PH has converged. Default: 0.001.
    288 \item PH\_Iters\_Between\_Cycle\_Slams: controls the number of iterations to wait after a variable is slammed due to
    289 hash hits that suggest convergence. Zero indicates unlimited slams per cycle. Default: 1.
    290 \item SlamAfterIter: Iteration number after which one variable every other iteration will be slammed to force convergence. Default: the number of scenarios.
    291 \item hash\_hit\_len\_to\_slam: Ignore possible cycles for which the only evidence of a cycle is less than this. Also, ignore cycles if any variables have been fixed in the previous hash\_hit\_len\_to\_slam PH iterations. Default:
    292 the number of scenarios. This default is often not a good choice. For many problems with a lot of scenarios,
    293 a value like 10 or 20 might be a lot better if rapid convergence is desired.
    294 \item W\_hash\_history\_len: This obscure parameter controls how far back the code will look to see if there is a possible cycle. Default: max(100, number of scenarios).
     459  \item Iteration0Mipgap: Gives the mipgap to be sent to the solver for
     460  iteration zero solves. The default is zero, which causes nothing to be sent to
     461  the solver (i.e., the solver uses its default mipgap).
     462
     463  \item InitalMipGap: Gives the mipgap to be sent to the solver for iteration
     464  one solves. The default is zero, which causes nothing to be sent to the solver
     465  (i.e., the solver uses its default mipgap). If not zero, then this gap will be
     466  used as the starting point for the mipgap to change on each PH iteration in
     467  proportion to the convergence termination criterion so that when the criterion
     468  for termination is met the mipgap will be at (or near) the parameter value of
     469  FinalMipGap. If the InitialMipGap is significantly higher than the
     470  Iteration0MipGap parameter, the PH algorithm may perform poorly. This is
     471  because the iteration k-1 solutions are used to warm start iteration k-1
     472  solves. If the iteration 1 mipgap is much higher than the iteration 0 mipgap,
     473  the iteration zero solution, although not optimal for the iteration one
     474  objective, might be within the mipgap. Default: 0.10.
     475
     476  \item FinalMipGap: The target for the mipgap when PH has converged. Default:
     477  0.001.
     478
     479  \item PH\_Iters\_Between\_Cycle\_Slams: controls the number of iterations to
     480  wait after a variable is slammed due to hash hits that suggest convergence.
     481  Zero indicates unlimited slams per cycle. Default: 1.
     482
     483  \item SlamAfterIter: Iteration number after which one variable every other
     484  iteration will be slammed to force convergence. Default: the number of
     485  scenarios.
     486
     487  \item hash\_hit\_len\_to\_slam: Ignore possible cycles for which the only
     488  evidence of a cycle is less than this. Also, ignore cycles if any variables
     489  have been fixed in the previous hash\_hit\_len\_to\_slam PH iterations.
     490  Default: the number of scenarios. This default is often not a good choice. For
     491  many problems with a lot of scenarios, a value like 10 or 20 might be a lot
     492  better if rapid convergence is desired.
     493
     494  \item W\_hash\_history\_len: This obscure parameter controls how far back the
     495  code will look to see if there is a possible cycle. Default: max(100, number
     496  of scenarios).
    295497\end{itemize}
    296498
    297499\subsubsection{Setting Parameter Values \label{ParmSec}}
    298500
    299 The parameters of PH and of any callbacks can be changed using the file wwph.cfg, which is executed by the python interpreter after PH has initialized. This is
    300 a potentially powerful and/or dangerous procedure because it gives an opportunity
    301 to change anything using the full features of python and pyomo.
     501The parameters of PH and of any callbacks can be changed using the file
     502wwph.cfg, which is executed by the python interpreter after PH has initialized.
     503This is a potentially powerful and/or dangerous procedure because it gives an
     504opportunity to change anything using the full features of python and pyomo.
    302505
    303506\subsubsection{Setting Suffix Values \label{SuffixSec}}
    304507
    305 Suffixes are set using the data file named \verb|wwph.suffixes| using this syntax:
     508Suffixes are set using the data file named \verb|wwph.suffixes| using this
     509syntax:
    306510
    307511\begin{verbatim}
     
    309513\end{verbatim}
    310514
    311 where VARSPEC is replaced by a variable specification, SUFFIX is replaced by a suffix name and VALUE is replaced by the value of the suffix for
    312 that variable or those variables. Here is an example:
     515where VARSPEC is replaced by a variable specification, SUFFIX is replaced by a
     516suffix name and VALUE is replaced by the value of the suffix for that variable
     517or those variables. Here is an example:
    313518
    314519\begin{verbatim}
     
    322527\subsection{Callback Details \label{ExtensionDetailsSec}}
    323528
    324 Most users of pysp can skip this subsection.
    325 A callback class definition named iphextension is in the file iphextension.py and can be used to implement callbacks at a
    326 variety of points in PH. For example, the method post\_iteration\_0\_solves is called immediately after all iteration zero
    327 solves, but before averages and weights have been computed while the method post\_iteration\_0 is called after averages and
    328 weights based on iteration zero have been computed. The file iphextension is in the coopr/pysp directory and is not
    329 intended to be edited by users.
    330 
    331 The user defines a class derived from SingletonPlugin that implements iphextension. Its name is given to ph as an option. This class will be automatically instantiated by ph.
    332 It has access to data and methods in the PH class, which are defined in the file ph.py. An example of such a class is in the
    333 file named testphextension.py in the pysp example directory.
    334 
    335 A user defined extension file can be incorporated by using the command line option:
    336 \verb|--user-defined-extension=EXTENSIONFILE|. Here, "EXTENSIONFILE" is the module name, which is in either the current directory (most likely) or somewhere on your PYTHONPATH. A simple example is "testphextension" plugin that simply prints a message to the screen for each callback. The file testphextension.py can be found in the sources directory and given in Section~\ref{CallbackSec}. An easy test of this would be to specify "-user-defined-extension=testphextension" and you should
    337 note the the ``.py'' file extension is not included on the runph command line.
    338  
    339 Both your own extension and the WWPH extensions can co-exist; however, the WW plugin will be invoked first at each callback point if both are included.
     529Most users of pysp can skip this subsection. A callback class definition named
     530iphextension is in the file iphextension.py and can be used to implement
     531callbacks at a variety of points in PH. For example, the method
     532post\_iteration\_0\_solves is called immediately after all iteration zero
     533solves, but before averages and weights have been computed while the method
     534post\_iteration\_0 is called after averages and weights based on iteration zero
     535have been computed. The file iphextension is in the coopr/pysp directory and is
     536not intended to be edited by users.
     537
     538The user defines a class derived from SingletonPlugin that implements
     539iphextension. Its name is given to ph as an option. This class will be
     540automatically instantiated by ph. It has access to data and methods in the PH
     541class, which are defined in the file ph.py. An example of such a class is in the
     542file named testphextension.py in the pysp example directory.
     543
     544A user defined extension file can be incorporated by using the command line
     545option: \verb|--user-defined-extension=EXTENSIONFILE|. Here, "EXTENSIONFILE" is
     546the module name, which is in either the current directory (most likely) or
     547somewhere on your PYTHONPATH. A simple example is "testphextension" plugin that
     548simply prints a message to the screen for each callback. The file
     549testphextension.py can be found in the sources directory and given in
     550Section~\ref{CallbackSec}. An easy test of this would be to specify
     551"-user-defined-extension=testphextension" and you should note the the ``.py''
     552file extension is not included on the runph command line.
     553
     554Both your own extension and the WWPH extensions can co-exist; however, the WW
     555plugin will be invoked first at each callback point if both are included.
    340556
    341557Here are the callbacks:
    342558\begin{itemize}
    343 \item post\_ph\_initialization: Called after PH data structures have been intialized but before iteration zero solves.
    344 \item post\_iteration\_0\_solves: Called after iteration zero solutions and some statistics such as averages have been computed, but before weights are updated.
    345 \item post\_iteration\_0: Called after all processing for iteration zero is complete.
    346 \item post\_iteration\_k\_solves: Called after solutions some statistics such as averages have been computed, but before weights are updated for iterations after iteration zero.
    347 \item post\_iteration\_k: Called after all processing for each iteration after iteration 0 is complete.
    348 \item post\_ph\_execution: Called execution is complete.
     559  \item post\_ph\_initialization: Called after PH data structures have been
     560  intialized but before iteration zero solves.
     561
     562  \item post\_iteration\_0\_solves: Called after iteration zero solutions and
     563  some statistics such as averages have been computed, but before weights are
     564  updated.
     565
     566  \item post\_iteration\_0: Called after all processing for iteration zero is
     567  complete.
     568
     569  \item post\_iteration\_k\_solves: Called after solutions some statistics such
     570  as averages have been computed, but before weights are updated for iterations
     571  after iteration zero.
     572
     573  \item post\_iteration\_k: Called after all processing for each iteration after
     574  iteration 0 is complete.
     575
     576  \item post\_ph\_execution: Called execution is complete.
    349577\end{itemize}
    350578
    351 Users interested in writing their own extensions will probably want to refer to the source file ph.py to get a deeper understanding of when the callback occur.
     579Users interested in writing their own extensions will probably want to refer to
     580the source file ph.py to get a deeper understanding of when the callback occur.
    352581
    353582\section{Examples \label{ExampleSec}}
     
    357586\subsection{Farmer Example}
    358587
    359 This two-stage example is composed of models and data for the "Farmer" stochastic program, introduced in Section 1.1 of
    360 "Introduction to Stochastic Programming" by Birge and
    361 Louveaux \cite{spbook}.
     588This two-stage example is composed of models and data for the "Farmer"
     589stochastic program, introduced in Section 1.1 of "Introduction to Stochastic
     590Programming" by Birge and Louveaux \cite{spbook}.
    362591
    363592\begin{itemize}
    364 \item ReferenceModel.py: a single-scenario model for the SP
    365 \item ReferenceModel.dat: a single-scenario data file for the SP (any scenario will do - used to flush out variable and constraint index sets)
    366 \item ScenarioStructure.dat: data file defining the scenario tree.
    367 \item AboveAverageScenario.dat: one of the scenario data files.
    368 \item BelowAverageScenario.dat: one of the scenario data files.
    369 \item AverageScenario.dat: one of the scenario data files.
     593  \item ReferenceModel.py: a single-scenario model for the SP
     594
     595  \item ReferenceModel.dat: a single-scenario data file for the SP (any scenario
     596  will do - used to flush out variable and constraint index sets)
     597
     598  \item ScenarioStructure.dat: data file defining the scenario tree.
     599
     600  \item AboveAverageScenario.dat: one of the scenario data files.
     601
     602  \item BelowAverageScenario.dat: one of the scenario data files.
     603
     604  \item AverageScenario.dat: one of the scenario data files.
    370605\end{itemize}
    371606
    372 The command runph executes PH, assuming the ReferenceModel.* and ScenarioStructure.* files are present and correct.
     607The command runph executes PH, assuming the ReferenceModel.* and
     608ScenarioStructure.* files are present and correct. This example is probably in a
     609directory with a name something like:
     610
     611\begin{verbatim}
     612coopr\examples\pysp\farmer
     613\end{verbatim}
     614
     615The data is in a subdirectory called scenariodata and the model is in the models
     616subdirectory. To invoke PH for this problem, connect to this farmer directory
     617and use the command:
     618
     619\begin{verbatim}
     620runph --model-directory=models --instance-directory=scenariodata
     621\end{verbatim}
     622
     623\subsection{Sizes Example}
     624
     625This two-stage example is composed of models and data for the "Sizes" stochastic
     626program \cite{sizes,lokwood}, which consists of the following files:
     627
     628\begin{itemize}
     629  \item wwph.cfg: replace default algorithm parameter values for the Watson and
     630  Woodruff extensions.
     631
     632  \item wwph.suffixes: sets algorithm parameter values at the variables level
     633  for the Watson and Woodruff extensions.
     634
     635  \item ReferenceModel.py: a single-scenario model for the SP
     636
     637  \item ReferenceModel.dat: a single-scenario data file for the SP (any scenario
     638  will do - used to flush out variable and constraint index sets)
     639
     640  \item ScenarioStructure.dat: data file defining the scenario tree.
     641
     642  \item Scenario1.dat: one of the scenario data files.
     643
     644  \item Scenario2.dat: one of the scenario data files.
     645
     646  \item ...
     647\end{itemize}
     648
    373649This example is probably in a directory with a name something like:
    374650
    375651\begin{verbatim}
    376 coopr\examples\pysp\farmer
    377 \end{verbatim}
    378 
    379 The data is in a subdirectory called scenariodata and the model is in the models subdirectory. To
    380 invoke PH for this problem, connect to this farmer directory and use the command:
    381 
    382 \begin{verbatim}
    383 runph --model-directory=models --instance-directory=scenariodata
    384 \end{verbatim}
    385 
    386 \subsection{Sizes Example}
    387 
    388 This two-stage example is composed of models and data for the "Sizes" stochastic program \cite{sizes,lokwood}, which consists of the following files:
     652coopr\packages\coopr\examples\pysp\sizes
     653\end{verbatim}
     654
     655The data for a three scenario version is in a subdirectory called SIZES3 and a
     656ten scenario dataset is in SIZES10.
     657
     658To invoke PH for the 10 scenario problem, connect to the sizes directory and use
     659the command:
     660
     661\begin{verbatim}
     662runph --model-directory=models --instance-directory=SIZES10
     663\end{verbatim}
     664
     665
     666\subsection{Forestry Example}
     667
     668This four-stage example is composed of models and data for the ``forestry''
     669stochastic program \cite{}, which consists of the following files:
    389670
    390671\begin{itemize}
    391 \item wwph.cfg: replace default algorithm parameter values for the Watson and Woodruff extensions.
    392 \item wwph.suffixes: sets algorithm parameter values at the variables level for the Watson and Woodruff extensions.
    393 \item ReferenceModel.py: a single-scenario model for the SP
    394 \item ReferenceModel.dat: a single-scenario data file for the SP (any scenario will do - used to flush out variable and constraint index sets)
    395 \item ScenarioStructure.dat: data file defining the scenario tree.
    396 \item Scenario1.dat: one of the scenario data files.
    397 \item Scenario2.dat: one of the scenario data files.
    398 \item ...
     672  \item wwph.cfg: replace default algorithm parameter values for the Watson and
     673  Woodruff extensions.
     674
     675  \item wwph.suffixes: sets algorithm parameter values at the variables level
     676  for the Watson and Woodruff extensions.
     677
     678  \item ReferenceModel.py: a single-scenario model for the SP
     679
     680  \item ReferenceModel.dat: a single-scenario data file for the SP (any scenario
     681  will do - used to flush out variable and constraint index sets)
     682
     683  \item ScenarioStructure.dat: data file defining the scenario tree.
     684
     685  \item Scenario1.dat: one of the scenario data files.
     686
     687  \item Scenario2.dat: one of the scenario data files.
     688
     689  \item ...
    399690\end{itemize}
    400691
     
    402693
    403694\begin{verbatim}
    404 coopr\packages\coopr\examples\pysp\sizes
    405 \end{verbatim}
    406 
    407 The data for a three scenario version is in a subdirectory called SIZES3 and a ten scenario
    408 dataset is in SIZES10.
    409 
    410 To invoke PH for the 10 scenario problem, connect to the sizes directory and use the command:
    411 
    412 \begin{verbatim}
    413 runph --model-directory=models --instance-directory=SIZES10
    414 \end{verbatim}
    415 
    416 
    417 \subsection{Forestry Example}
    418 
    419 This four-stage example is composed of models and data for the ``forestry'' stochastic program \cite{}, which consists of the following files:
    420 
    421 \begin{itemize}
    422 \item wwph.cfg: replace default algorithm parameter values for the Watson and Woodruff extensions.
    423 \item wwph.suffixes: sets algorithm parameter values at the variables level for the Watson and Woodruff extensions.
    424 \item ReferenceModel.py: a single-scenario model for the SP
    425 \item ReferenceModel.dat: a single-scenario data file for the SP (any scenario will do - used to flush out variable and constraint index sets)
    426 \item ScenarioStructure.dat: data file defining the scenario tree.
    427 \item Scenario1.dat: one of the scenario data files.
    428 \item Scenario2.dat: one of the scenario data files.
    429 \item ...
    430 \end{itemize}
    431 
    432 This example is probably in a directory with a name something like:
    433 
    434 \begin{verbatim}
    435695coopr\packages\coopr\examples\pysp\forestry
    436696\end{verbatim}
    437697
    438 
    439 There are two families of instances: ``Chile'' and ``Davis,'' each with four stages and eighteen scenarios. This is also a small
    440 two-stage, four scenario instances in the subdirectory DAVIS2STAGE.
    441 
    442 This full 18 scenario problem instance takes too long without the Watson Woodruff extensions, but to
    443 invoke PH for this problem without them, connect to the forestry directory and use the command:
     698There are two families of instances: ``Chile'' and ``Davis,'' each with four
     699stages and eighteen scenarios. This is also a small two-stage, four scenario
     700instances in the subdirectory DAVIS2STAGE.
     701
     702This full 18 scenario problem instance takes too long without the Watson
     703Woodruff extensions, but to invoke PH for this problem without them, connect to
     704the forestry directory and use the command:
    444705
    445706\begin{verbatim}
    446707 runph --models-directory=models --instancs-directory=chile
    447708\end{verbatim}
     709
    448710and to run with the extensions, use
     711
    449712\begin{verbatim}
    450713runph --model-directory=models --instance-directory=davis \
     
    454717
    455718
    456 \section{Tutorial: Parameters for Watson and Woodruff PH Extensions \label{WWTutorialSec}}
    457 
    458 The parameters for the PH extensions in WWPHExtensions.py provide the user with considerable control over
    459 how and under what conditions variables are fixed during the PH algorithm execution. Often, some variables converge
    460 to consistent values during early iterations and can be fixed at these values without affecting quality very much and without
    461 affecting feasibility at all. Also, the algorithm may need to fix some variables during execution in order to break cycles.
    462 In both cases, guidance from the user concerning which classes of variables can and/or should be fixed under various circumstances can be very helpful.
    463 
    464 The overarching goal is to support industrial and government users who solve the same model repeatedly with different, but
    465 somewhat similar, data each time. In such settings, it behooves the modeler to consider tradeoffs between speed, solution quality
    466 and feasibility and create at least one good set of directives and parameters for the PH algorithm. In some cases, a user may want more than
    467 one set of directives and parameters depending on whether speed of execution or quality of solution are more important. Iteration zero is controlled separately because often the absence of the quadratic term results in faster solves for this
    468 iteration and fixing variables after the iteration has the maximum possible impact on speedup.
     719\section{Tutorial: Parameters for Watson and Woodruff PH Extensions
     720\label{WWTutorialSec}}
     721
     722The parameters for the PH extensions in WWPHExtensions.py provide the user with
     723considerable control over how and under what conditions variables are fixed
     724during the PH algorithm execution. Often, some variables converge to consistent
     725values during early iterations and can be fixed at these values without
     726affecting quality very much and without affecting feasibility at all. Also, the
     727algorithm may need to fix some variables during execution in order to break
     728cycles. In both cases, guidance from the user concerning which classes of
     729variables can and/or should be fixed under various circumstances can be very
     730helpful.
     731
     732The overarching goal is to support industrial and government users who solve the
     733same model repeatedly with different, but somewhat similar, data each time. In
     734such settings, it behooves the modeler to consider tradeoffs between speed,
     735solution quality and feasibility and create at least one good set of directives
     736and parameters for the PH algorithm. In some cases, a user may want more than
     737one set of directives and parameters depending on whether speed of execution or
     738quality of solution are more important. Iteration zero is controlled separately
     739because often the absence of the quadratic term results in faster solves for
     740this iteration and fixing variables after the iteration has the maximum possible
     741impact on speedup.
    469742
    470743In order to discuss these issues, we consider an example.
     
    472745\subsection{Sizes Example}
    473746
    474 The Sizes example is simple and small. In fact, the instances that we will consider are so small that there is really
    475 no need to use the PH algorithm since the extensive form can be solved directly in a few minutes. However, it serves as a
    476 vehicle for introducing the concepts.
    477 
    478 This description and formulation comes from the original paper by Jorjani, Scott and Woodruff \cite{sizes}.
    479 
    480 If demand constraints for the current period are
    481 based on firm orders but future demands are based on forecasts or
    482 conjecture, then they should not be treated in the same fashion. We
    483 must recognize that future demand constraints are stochastic. That is
    484 to say that they should be modeled as random variables.
    485 It may not be reasonable or useful to consider the entire demand probability
    486 distribution functions. It may not be reasonable because there may not be
    487 sufficient data to estimate an entire distribution. It may not be useful
    488 because the essence of the stochastics may be captured by specifying a
    489 small number of representative {\em scenarios}. We assume that scenarios
    490 are specified by giving a full set of random variable realizations and a
    491 corresponding probability.
    492 We index the scenario set, ${\cal L}$, by $\ell$ and
    493 refer to the probability of occurrence of $\ell$ (or, more accurately, a
    494 realization ``near'' scenario $\ell$) as $\Prl$.
    495 We refer to solution systems that satisfy constraints with probability one as
    496 {\em admissible}.
    497 
    498 In addition to modeling stochastics, we would like to model {\em
    499 recourse} as well. That is, we would like to model the ability of
    500 decision makers to make use of new information (e.g., orders) at the
    501 start of each planning period.  We allow our solution vectors to
    502 depend on the scenario that is realized, but we do not want to assume
    503 prescience. We refer to a system of solution vectors as {\em
    504 implementable} if for all decision times $t$, the solution vector
    505 elements corresponding to period $1,\ldots,t$ are constant with
    506 respect to information that becomes available only after stage $t$.
    507 We refer to the set of implementable solutions as ${\cal N}_{\cal L}$.  It is
    508 possible to require implementable solutions by adding {\em
    509 non-anticipatitivity constraints}, but we will instead make use of
    510 solution procedures that implicitly guarantee implementable solutions.
    511 
    512 In the stochastic, multi-period formulation that follows the objective
    513 is to minimize expected costs.  We invoke the network equivalence
    514 given earlier to drop the explicit requirement that $\Vector{x}$ and
    515 $\Vector{y}$ be integers. Variables and data are subscripted with a
    516 period index $t$ that takes on values up to $T$. To model the idea
    517 that sleeves produced in one period can be used as-is or cut in
    518 subsequent periods, we use $x_{ijt}$ to indicate that sleeves of
    519 length index $i$ are to be used without cutting in period $t$ if $i=j$
    520 and with cutting otherwise.  The $\Vector{y}$ vector gives production
    521 quantities for each length in each period without regard to the period
    522 in which they will be used (and perhaps cut). The formulation is
    523 essentially an extension of ILP except that a capacity constraint
    524 must be added in the multiple period formulation. Holding costs could
    525 be added, but an additional subscript becomes necessary without the
    526 benefit of any additional insight. As an aside, note that the addition
    527 of holding costs would add a large number of continuous variables, but
    528 no new integers so the impact on computational performance would not
    529 be catastrophic.
     747The Sizes example is simple and small. In fact, the instances that we will
     748consider are so small that there is really no need to use the PH algorithm since
     749the extensive form can be solved directly in a few minutes. However, it serves
     750as a vehicle for introducing the concepts.
     751
     752This description and formulation comes from the original paper by Jorjani, Scott
     753and Woodruff \cite{sizes}.
     754
     755If demand constraints for the current period are based on firm orders but future
     756demands are based on forecasts or conjecture, then they should not be treated in
     757the same fashion. We must recognize that future demand constraints are
     758stochastic. That is to say that they should be modeled as random variables. It
     759may not be reasonable or useful to consider the entire demand probability
     760distribution functions. It may not be reasonable because there may not be
     761sufficient data to estimate an entire distribution. It may not be useful because
     762the essence of the stochastics may be captured by specifying a small number of
     763representative {\em scenarios}. We assume that scenarios are specified by giving
     764a full set of random variable realizations and a corresponding probability. We
     765index the scenario set, ${\cal L}$, by $\ell$ and refer to the probability of
     766occurrence of $\ell$ (or, more accurately, a realization ``near'' scenario
     767$\ell$) as $\Prl$. We refer to solution systems that satisfy constraints with
     768probability one as {\em admissible}.
     769
     770In addition to modeling stochastics, we would like to model {\em recourse} as
     771well. That is, we would like to model the ability of decision makers to make use
     772of new information (e.g., orders) at the start of each planning period. We allow
     773our solution vectors to depend on the scenario that is realized, but we do not
     774want to assume prescience. We refer to a system of solution vectors as {\em
     775implementable} if for all decision times $t$, the solution vector elements
     776corresponding to period $1,\ldots,t$ are constant with respect to information
     777that becomes available only after stage $t$. We refer to the set of
     778implementable solutions as ${\cal N}_{\cal L}$. It is possible to require
     779implementable solutions by adding {\em non-anticipatitivity constraints}, but we
     780will instead make use of solution procedures that implicitly guarantee
     781implementable solutions.
     782
     783In the stochastic, multi-period formulation that follows the objective is to
     784minimize expected costs. We invoke the network equivalence given earlier to drop
     785the explicit requirement that $\Vector{x}$ and $\Vector{y}$ be integers.
     786Variables and data are subscripted with a period index $t$ that takes on values
     787up to $T$. To model the idea that sleeves produced in one period can be used
     788as-is or cut in subsequent periods, we use $x_{ijt}$ to indicate that sleeves of
     789length index $i$ are to be used without cutting in period $t$ if $i=j$ and with
     790cutting otherwise. The $\Vector{y}$ vector gives production quantities for each
     791length in each period without regard to the period in which they will be used
     792(and perhaps cut). The formulation is essentially an extension of ILP except
     793that a capacity constraint must be added in the multiple period formulation.
     794Holding costs could be added, but an additional subscript becomes necessary
     795without the benefit of any additional insight. As an aside, note that the
     796addition of holding costs would add a large number of continuous variables, but
     797no new integers so the impact on computational performance would not be
     798catastrophic.
    530799
    531800\[
     
    536805subject to
    537806\begin{eqnarray}
    538 \sum_{j \geq i}x_{ijt\ell} & \geq & d_{it\ell} 
     807\sum_{j \geq i}x_{ijt\ell} & \geq & d_{it\ell}
    539808          \;\;\;\; \ell \in {\cal L}, \; i=1,\ldots,N, \; t=1,\ldots,T \label{Dconstr} \\
    540 \sum_{t'\leq t}\left[\sum_{j \leq i}x_{ijt'\ell} - y_{it'\ell}\right] & \leq & 0 
     809\sum_{t'\leq t}\left[\sum_{j \leq i}x_{ijt'\ell} - y_{it'\ell}\right] & \leq & 0
    541810          \;\;\;\; \ell \in {\cal L}, \; i=1,\ldots,N, \; t=1,\ldots,T \label{Pconstr} \\
    542 y_{it\ell} - Mz_{it\ell} & \leq & 0 
     811y_{it\ell} - Mz_{it\ell} & \leq & 0
    543812          \;\;\;\; \ell \in {\cal L}, \;i=1,\ldots,N, \; t=1,\ldots,T \label{Mconstr} \\
    544813\sum_{i=1}^{N}y_{it\ell} & \leq & c_{t\ell}
    545814          \;\;\;\; \ell \in {\cal L}, \; t=1,\ldots,T \label{Cconstr} \\
    546 z_{it\ell} & \in & \{0,1\} 
     815z_{it\ell} & \in & \{0,1\}
    547816          \;\;\;\; \ell \in {\cal L}, \; i=1,\ldots,N, \; t=1,\ldots,T \label{intconstr} \\
    548 \Vector{x}, \Vector{y}, \Vector{z} & \in & {\cal N}_{\cal L} 
     817\Vector{x}, \Vector{y}, \Vector{z} & \in & {\cal N}_{\cal L}
    549818\end{eqnarray}
    550819
    551 Bear in mind that solution vector elements corresponding
    552 to periods two through $T$ are not actually intended for use, they
    553 are computed just to see the effect that period 1 (the current period)
    554 decision would have on future optimal behavior. At the start of period 2
    555 -- or at any other time -- the decision maker would run the model again
    556 with updated demands for the current period and new scenario estimates.
     820Bear in mind that solution vector elements corresponding to periods two through
     821$T$ are not actually intended for use, they are computed just to see the effect
     822that period 1 (the current period) decision would have on future optimal
     823behavior. At the start of period 2 -- or at any other time -- the decision maker
     824would run the model again with updated demands for the current period and new
     825scenario estimates.
    557826
    558827\subsection{ReferenceModel.py}
     
    614883   for i in range(1,model.NumSizes()+1):
    615884      for j in range(1, i+1):
    616          ans.add((i,j))   
     885         ans.add((i,j))
    617886   return ans
    618887
     
    673942
    674943def enforce_production_second_stage_rule(i, model):
    675    # The production capacity per time stage serves as a simple upper bound for "M".   
     944   # The production capacity per time stage serves as a simple upper bound for "M".
    676945   return (None, \
    677946           model.NumProducedSecondStage[i] - model.Capacity * model.ProduceSizeSecondStage[i], \
     
    692961   return (None, \
    693962           sum(model.NumProducedSecondStage[i] for i in model.ProductSizes) - model.Capacity, \
    694            0.0)   
     963           0.0)
    695964
    696965model.EnforceCapacityLimitFirstStage = Constraint(rule=enforce_capacity_first_stage_rule)
     
    7391008             sum(model.UnitReductionCost * model.NumUnitsCutSecondStage[i,j] \
    7401009                 for (i,j) in model.NumUnitsCutDomain if i != j)
    741    return (model.SecondStageCost - production_costs - cut_costs) == 0.0   
     1010   return (model.SecondStageCost - production_costs - cut_costs) == 0.0
    7421011
    7431012model.ComputeSecondStageCost = Constraint(rule=second_stage_cost_rule)
     
    7561025\subsection{ReferenceModel.dat}
    7571026
    758 This file establishes the data for a representative instance. The main thing to notice is that the indexes for
    759 sizes happen to be given as integers, which is not required: they could have been strings.
     1027This file establishes the data for a representative instance. The main thing to
     1028notice is that the indexes for sizes happen to be given as integers, which is
     1029not required: they could have been strings.
    7601030
    7611031{\small
     
    7891059set Stages := FirstStage SecondStage ;
    7901060
    791 set Nodes := RootNode 
     1061set Nodes := RootNode
    7921062             Scenario1Node
    7931063             Scenario2Node
     
    8011071             Scenario10Node ;
    8021072
    803 param NodeStage := RootNode             FirstStage 
     1073param NodeStage := RootNode             FirstStage
    8041074                   Scenario1Node        SecondStage
    8051075                   Scenario2Node        SecondStage
     
    8361106                                Scenario10Node    0.10 ;
    8371107
    838 set Scenarios := Scenario1 
    839                  Scenario2 
    840                  Scenario3 
    841                  Scenario4 
    842                  Scenario5 
    843                  Scenario6 
    844                  Scenario7 
    845                  Scenario8 
    846                  Scenario9 
    847                  Scenario10 ; 
    848 
    849 param ScenarioLeafNode := Scenario1  Scenario1Node 
    850                           Scenario2  Scenario2Node 
    851                           Scenario3  Scenario3Node 
    852                           Scenario4  Scenario4Node 
    853                           Scenario5  Scenario5Node 
    854                           Scenario6  Scenario6Node 
    855                           Scenario7  Scenario7Node 
    856                           Scenario8  Scenario8Node 
    857                           Scenario9  Scenario9Node 
     1108set Scenarios := Scenario1
     1109                 Scenario2
     1110                 Scenario3
     1111                 Scenario4
     1112                 Scenario5
     1113                 Scenario6
     1114                 Scenario7
     1115                 Scenario8
     1116                 Scenario9
     1117                 Scenario10 ;
     1118
     1119param ScenarioLeafNode := Scenario1  Scenario1Node
     1120                          Scenario2  Scenario2Node
     1121                          Scenario3  Scenario3Node
     1122                          Scenario4  Scenario4Node
     1123                          Scenario5  Scenario5Node
     1124                          Scenario6  Scenario6Node
     1125                          Scenario7  Scenario7Node
     1126                          Scenario8  Scenario8Node
     1127                          Scenario9  Scenario9Node
    8581128                          Scenario10 Scenario10Node ;
    8591129
     
    8671137param StageCostVariable := FirstStage FirstStageCost
    8681138                           SecondStage SecondStageCost ;
    869                                    
     1139
    8701140\end{verbatim}
    8711141}
     
    8731143\subsection{wwph.cfg}
    8741144
    875 This file provides default values. Even though most of them will be overridden at the variable level, it makes sense
    876 to provide instance specific defaults. For this problem, it does not seem helpful or prudent to do any fixing of continuous
    877 variables since they only used to compute the objective function terms. The \verb|ProduceSizeFirstStage| variables are binary
    878 and fixing their values would seem to determine the value of the other values. The other first stage variables are general integers. The \verb|ProduceSizeFirstStage| can safely be fixed because provided that the largest size is not fixed at zero, there is little risk of infeasibility since larger sizes can be cut (at a cost) to meet demand for smaller sizes. Consequently, it is safe to let the algorithm slam those
    879 variables as needed. Slamming the other variables is riskier because they could get slammed to values that don't make
    880 sense given the values of the \verb|ProduceSizeFirstStage| variables.
    881 
    882 Fixing variables that have converged will speed the algorithm, perhaps resulting in a less desirable
    883 objective function value. It would make sense to tend to fix the \verb|ProduceSizeFirstStage| before fixing the others to avoid
    884 inconsistencies and because the The \verb|ProduceSizeFirstStage| variables have a strong tendency to imply values for the other
    885 variables.
    886 
    887 Here is a sensible and internally consistent, if a bit conservative, wwph.cfg file:
     1145This file provides default values. Even though most of them will be overridden
     1146at the variable level, it makes sense to provide instance specific defaults. For
     1147this problem, it does not seem helpful or prudent to do any fixing of continuous
     1148variables since they only used to compute the objective function terms. The
     1149\verb|ProduceSizeFirstStage| variables are binary and fixing their values would
     1150seem to determine the value of the other values. The other first stage variables
     1151are general integers. The \verb|ProduceSizeFirstStage| can safely be fixed
     1152because provided that the largest size is not fixed at zero, there is little
     1153risk of infeasibility since larger sizes can be cut (at a cost) to meet demand
     1154for smaller sizes. Consequently, it is safe to let the algorithm slam those
     1155variables as needed. Slamming the other variables is riskier because they could
     1156get slammed to values that don't make sense given the values of the
     1157\verb|ProduceSizeFirstStage| variables.
     1158
     1159Fixing variables that have converged will speed the algorithm, perhaps resulting
     1160in a less desirable objective function value. It would make sense to tend to fix
     1161the \verb|ProduceSizeFirstStage| before fixing the others to avoid
     1162inconsistencies and because the The \verb|ProduceSizeFirstStage| variables have
     1163a strong tendency to imply values for the other variables.
     1164
     1165Here is a sensible and internally consistent, if a bit conservative, wwph.cfg
     1166file:
    8881167
    8891168{\small
     
    9071186self.FixWhenItersConvergedAtNB = 0  # converged to a non-bound
    9081187self.FixWhenItersConvergedContinuous = 0
    909      
    910 # "default" slamming parms 
     1188
     1189# "default" slamming parms
    9111190# True and False are the options (case sensitive)
    9121191self.CanSlamToLB = False
     
    9271206}
    9281207
    929 There are a number of things to notice about the contents of this file. Since it is executed
    930 as python code, the syntax matters. Users should changes values of numbers of change True to False, but
    931 everything else should not be edited with the exception of comments, which is any text after a sharp sign
    932 (sometimes called a pound sign). Changes to this file should be tested incrementally because errors are trapped by
    933 the python interpreter and may be difficult for non-python programmers to decipher.
    934 
    935 This particular example file allows variables to be fixed if all scenarios have agreed on the
    936 upper bound for five iterations in a row. Since the \verb|ProduceSizeFirstStage| variables are the only discrete
    937 variables in the first stage with an upper bound, they are the only variables affected by this.
    938 
    939 This example allows slamming, but only to the max across scenarios. This is different than the upper bound, even for
    940 binary variables, because a variable could be selected for slamming for which all scenarios have agreed on the same value
    941 (which could be the lower lower bound). Data providing an override for this default as well as control over the selection priority for variables to slam are provided in the wwph.suffixes file.
     1208There are a number of things to notice about the contents of this file. Since it
     1209is executed as python code, the syntax matters. Users should changes values of
     1210numbers of change True to False, but everything else should not be edited with
     1211the exception of comments, which is any text after a sharp sign (sometimes
     1212called a pound sign). Changes to this file should be tested incrementally
     1213because errors are trapped by the python interpreter and may be difficult for
     1214non-python programmers to decipher.
     1215
     1216This particular example file allows variables to be fixed if all scenarios have
     1217agreed on the upper bound for five iterations in a row. Since the
     1218\verb|ProduceSizeFirstStage| variables are the only discrete variables in the
     1219first stage with an upper bound, they are the only variables affected by this.
     1220
     1221This example allows slamming, but only to the max across scenarios. This is
     1222different than the upper bound, even for binary variables, because a variable
     1223could be selected for slamming for which all scenarios have agreed on the same
     1224value (which could be the lower lower bound). Data providing an override for
     1225this default as well as control over the selection priority for variables to
     1226slam are provided in the wwph.suffixes file.
    9421227
    9431228\subsection{wwph.suffixes}
     
    9471232# Optional suffixes to help control PH
    9481233
    949 # Text triples specifying (variable/variable-index, suffix-name, suffix-value) tuples.
     1234# Text triples specifying (variable/variable-index, suffix-name, suffix-value)
     1235# tuples.
    9501236
    9511237# override the defaults given in wwph.cfg as needed
     
    9601246ProduceSizeFirstStage[1] FixWhenItersConvergedAtLB 8
    9611247
    962 # if the production quantities have been fixed a long time, fix them 
     1248# if the production quantities have been fixed a long time, fix them
    9631249NumProducedFirstStage[*] FixWhenItersConvergedAtNB 20
    9641250
     
    9931279}
    9941280
    995 The first thing to notice is that this is a data file and not a file of Python comments.
    996 Consequently, simpler messages are provided if there are errors, but the file can be verbose
    997 and a computer program or spreadsheet might be required to produce this data file. The next thing to notice
    998 is that indexes can be specified either using valid index values, or wildcards.
    999 
    1000 This file overrides the defaults to allow some fixing after iteration zero. Fixing the smallest
    1001 size (index 1) to zero cannot result in infeasibility because larger sizes can be cut to
    1002 satisfy demand for that size. Along the same lines, aggressively fixing the
    1003 production indicator for larger sizes to 1 is also
    1004 safe and perhaps not sub-optimal if all scenarios ``want'' those sizes anyway.
     1281The first thing to notice is that this is a data file and not a file of Python
     1282comments. Consequently, simpler messages are provided if there are errors, but
     1283the file can be verbose and a computer program or spreadsheet might be required
     1284to produce this data file. The next thing to notice is that indexes can be
     1285specified either using valid index values, or wildcards.
     1286
     1287This file overrides the defaults to allow some fixing after iteration zero.
     1288Fixing the smallest size (index 1) to zero cannot result in infeasibility
     1289because larger sizes can be cut to satisfy demand for that size. Along the same
     1290lines, aggressively fixing the production indicator for larger sizes to 1 is
     1291also safe and perhaps not sub-optimal if all scenarios ``want'' those sizes
     1292anyway.
    10051293
    10061294\subsection{rhosetter.cfg}
     
    10181306   self.setRhoAllScenarios(model_instance.NumProducedFirstStage[i], model_instance.UnitProductionCosts[i] * MyRhoFactor)
    10191307   for j in model_instance.ProductSizes:
    1020       if j <= i: 
     1308      if j <= i:
    10211309         self.setRhoAllScenarios(model_instance.NumUnitsCutFirstStage[i,j], model_instance.UnitReductionCost * MyRhoFactor)
    10221310
     
    10391327\end{verbatim}
    10401328
    1041 Since this instance is so small by modern standards, the enhancement are not needed and may even increase the total solution
    1042 time. It is provided to illustrate the features of the extensions, not to illustrate why you might need them. Much larger instances are are required for that.
     1329Since this instance is so small by modern standards, the enhancement are not
     1330needed and may even increase the total solution time. It is provided to
     1331illustrate the features of the extensions, not to illustrate why you might need
     1332them. Much larger instances are are required for that.
    10431333
    10441334\section{Linearizing the Proximal Term \label{LinearSec}}
     
    10461336\subsection{Introduction}
    10471337
    1048 For a decision variable $x$ the proximal term added to the objective function for each scenario subproblem at PH iteration $k$ is
     1338For a decision variable $x$ the proximal term added to the objective function
     1339for each scenario subproblem at PH iteration $k$ is
    10491340$$
    10501341\left(x - \bar{x}^{(k-1)}\right)^{2}
    10511342$$
    1052 where $\bar{x}^{(k-1)}$ is the average over the scenarios from the last iteration. Expanding the square reveals that the only
    1053 quadratic term is $x^{2}$. For binary variables, this is equal to $x$, although this is not the case when a relaxation is solved.
    1054 For binary variables, the default behaviour is to replace the quadratic term with $x$, but an option allows the quadratic to be
    1055 retained because this can effect the subproblem solution time due to the use of the quadratic term when the relaxations are solved
    1056 as part of the branch and bound process.
    1057 
    1058 For non-binary variables an option exists to replace the $x^2$ term with a piece-wise approximation and the number of breakpoints in the
    1059 approximation is under user control. This can have a significant effect on CPU time required to solve sub-problems because the presence of
    1060 the quadratic term increases the solution times significantly. However, linearization results in a time-quality tradeoff because increasing
    1061 the number of breakpoints increases the fidelity but each breakpoint introduces another (unseen) binary variable so solution times are generally
    1062 increased.
    1063 
    1064 A few strategies for placing the breakpoints are supported as command a line options:
    1065 \verb|--breakpoint-strategy|=BREAKPOINT\_STRATEGY. A value of 0
    1066 indicates uniform distribution. 1 indicates
    1067 breakpoints at the min and max for the node in the scenario tree, uniformly in-
    1068 between. A value of 2 indicates more aggressive concentration of
    1069 breakpoints near the observed node min/max.
    1070 
    1071 Upper and lower bounds for variables must be specified when the linearization option is chosen. This can be done by specifying bounds in the
    1072 reference model or by using the bounds-cfgfile command line option. It is, of course, best to use meaningful bounds provided
    1073 in the reference model; however, the modeller must be careful not use estimated bounds that are too tight since that could preclude
    1074 an optimal (or even a feasible) solution to the overall stochastic problem even though it might not cut off any optimal solutions for the
    1075 particular scenario. The use of a bounds cfgfile is an advanced topic, but enables the modeller to use bounds that cannot create
    1076 infeasibilities.
     1343where $\bar{x}^{(k-1)}$ is the average over the scenarios from the last
     1344iteration. Expanding the square reveals that the only quadratic term is $x^{2}$.
     1345For binary variables, this is equal to $x$, although this is not the case when a
     1346relaxation is solved. For binary variables, the default behaviour is to replace
     1347the quadratic term with $x$, but an option allows the quadratic to be retained
     1348because this can effect the subproblem solution time due to the use of the
     1349quadratic term when the relaxations are solved as part of the branch and bound
     1350process.
     1351
     1352For non-binary variables an option exists to replace the $x^2$ term with a
     1353piece-wise approximation and the number of breakpoints in the approximation is
     1354under user control. This can have a significant effect on CPU time required to
     1355solve sub-problems because the presence of the quadratic term increases the
     1356solution times significantly. However, linearization results in a time-quality
     1357tradeoff because increasing the number of breakpoints increases the fidelity but
     1358each breakpoint introduces another (unseen) binary variable so solution times
     1359are generally increased.
     1360
     1361A few strategies for placing the breakpoints are supported as command a line
     1362options: \verb|--breakpoint-strategy|=BREAKPOINT\_STRATEGY. A value of 0
     1363indicates uniform distribution. 1 indicates breakpoints at the min and max for
     1364the node in the scenario tree, uniformly in- between. A value of 2 indicates
     1365more aggressive concentration of breakpoints near the observed node min/max.
     1366
     1367Upper and lower bounds for variables must be specified when the linearization
     1368option is chosen. This can be done by specifying bounds in the reference model
     1369or by using the bounds-cfgfile command line option. It is, of course, best to
     1370use meaningful bounds provided in the reference model; however, the modeller
     1371must be careful not use estimated bounds that are too tight since that could
     1372preclude an optimal (or even a feasible) solution to the overall stochastic
     1373problem even though it might not cut off any optimal solutions for the
     1374particular scenario. The use of a bounds cfgfile is an advanced topic, but
     1375enables the modeller to use bounds that cannot create infeasibilities.
    10771376
    10781377\subsection{Bounds-cfgfile}
    10791378
    1080 Using a bounds cfgfile is an advanced topic. The modeler is writing python/pyomo code that will be executed by the
    1081 ph.py file that is the core of the PH algorithm. The first executable line in a bounds file is typically:
     1379Using a bounds cfgfile is an advanced topic. The modeler is writing python/pyomo
     1380code that will be executed by the ph.py file that is the core of the PH
     1381algorithm. The first executable line in a bounds file is typically:
    10821382
    10831383\begin{verbatim}
     
    10851385\end{verbatim}
    10861386
    1087 This establishes a local variable called ``model\_instance'' that can be used to refer to the model (of course, a different name can be
    1088 used, such as MODINS). The object ``self'' on the right hand side of this assignment refers to the core PH object. The model, in turn contains the parameters and
    1089 variables defined in the ReferenceModel.py file that can be accessed by name. For example, with the Farmer model, the cfg file
    1090 sets the uppter bound on DevotedAcreage to be value of the paramter TOTAL\_ACREAGE at intialization (since this is Python, the
    1091 parentheses after TOTAL\_ACREAGE cause the value in TOTAL\_ACREAGE to be assigned rather than the name of the parameter object:
    1092 
    1093 \begin{verbatim}
    1094 # only need to set upper bounds on first-stage variables, i.e., those being blended.
    1095 
     1387This establishes a local variable called ``model\_instance'' that can be used to
     1388refer to the model (of course, a different name can be used, such as MODINS).
     1389The object ``self'' on the right hand side of this assignment refers to the core
     1390PH object. The model, in turn contains the parameters and variables defined in
     1391the ReferenceModel.py file that can be accessed by name. For example, with the
     1392Farmer model, the cfg file sets the uppter bound on DevotedAcreage to be value
     1393of the paramter TOTAL\_ACREAGE at intialization (since this is Python, the
     1394parentheses after TOTAL\_ACREAGE cause the value in TOTAL\_ACREAGE to be
     1395assigned rather than the name of the parameter object:
     1396
     1397\begin{verbatim}
     1398# only need to set upper bounds on first-stage variables, i.e., those being
     1399# blended.
    10961400model_instance = self._model_instance
    10971401
    1098 # the values supplied to the method 
     1402# the values supplied to the method
    10991403upper_bound = float(model_instance.TOTAL_ACREAGE())
    11001404
    11011405for index in model_instance.CROPS:
    1102    # the lower and upper bounds are expected to be floats, and trigger an exception if not.
     1406   # the lower and upper bounds are expected to be floats, and trigger an
     1407   # exception if not.
    11031408   self.setVariableBoundsAllScenarios("DevotedAcreage", index, 0.0, upper_bound)
    11041409\end{verbatim}
    11051410
    1106 The same thing could be accomplished by setting the upper bound in the model file, but it does serve as a simple example of a bounds cfg file.
     1411The same thing could be accomplished by setting the upper bound in the model
     1412file, but it does serve as a simple example of a bounds cfg file.
    11071413
    11081414%\subsection{Sizes Example}
     
    11111417\section{Solving sub-problems in Parallel \label{parallelsec}}
    11121418
    1113 Pyomo makes use of the pyro facility to enable sub-problems to easily be assigned to be solved in parallel. This capability is suppored by pysp. We will refer to a single master computer
    1114 and multiple slave computers in the discussion here, but actually, the master computing processes can (and for synchronous parallel, should) be on a
    1115 processor that also runs a slave process.
     1419Pyomo makes use of the pyro facility to enable sub-problems to easily be
     1420assigned to be solved in parallel. This capability is suppored by pysp. We will
     1421refer to a single master computer and multiple slave computers in the discussion
     1422here, but actually, the master computing processes can (and for synchronous
     1423parallel, should) be on a processor that also runs a slave process.
    11161424
    11171425Here are the commands in order:
     1426
    11181427\begin{enumerate}
    1119 \item On the master: \verb|coopr-ns|
    1120 \item On the master: \verb|dispatch_srvr|
    1121 \item On each slave: \verb|pyro_mip_server|
    1122 \item On the master: \verb|runph ... --solver-manager=pyro ...|
     1428  \item On the master: \verb|coopr-ns|
     1429  \item On the master: \verb|dispatch_srvr|
     1430  \item On each slave: \verb|pyro_mip_server|
     1431  \item On the master: \verb|runph ... --solver-manager=pyro ...|
    11231432\end{enumerate}
    1124 Note that the command \verb|coopr-ns| and the argument \verb|solver-manger| have a dash in the middle, while
    1125 the commands \verb|dispatch_srvr| and \verb|pyro_mip_server| have underscores. The first three commands launch processes
    1126 that have no internal mechanism for termination; i.e., they will be terminated only if they crash or if they are
    1127 killed by an external process. It is common to launch these processes with output redirection, such as \verb|coopr-ns >& cooprns.log|. The
    1128 \verb|runph| command is a normal runph command with the usual arguments with the additional specification that subproblem solves should
    1129 be directed to the pyro solver manager.
     1433
     1434Note that the command \verb|coopr-ns| and the argument \verb|solver-manger| have
     1435a dash in the middle, while the commands \verb|dispatch_srvr| and
     1436\verb|pyro_mip_server| have underscores. The first three commands launch
     1437processes that have no internal mechanism for termination; i.e., they will be
     1438terminated only if they crash or if they are killed by an external process. It
     1439is common to launch these processes with output redirection, such as
     1440\verb|coopr-ns >& cooprns.log|. The \verb|runph| command is a normal runph
     1441command with the usual arguments with the additional specification that
     1442subproblem solves should be directed to the pyro solver manager.
  • coopr.pysp/stable/2.4/examples/pysp/farmer/asyncphdriver.py

    r1646 r2938  
    1313
    1414# for profiling
    15 import cProfile
     15try:
     16    import cProfile as profile
     17except ImportError:
     18    import profile
    1619import pstats
    1720
     
    110113      traceback.print_exc()     
    111114else:
    112    cProfile.run('run_ph()','profile.stats')
     115   profile.run('run_ph()','profile.stats')
    113116   p=pstats.Stats('profile.stats')
    114117   p.sort_stats('time')
  • coopr.pysp/stable/2.4/examples/pysp/farmer/farmer_lp.dat

    r1739 r2938  
    1717param PlantingCostPerAcre := WHEAT 150 CORN 230 SUGAR_BEETS 260 ;
    1818
    19 param MeanYield := WHEAT 2.5 CORN 3 SUGAR_BEETS 20 ;
     19param Yield := WHEAT 2.5 CORN 3 SUGAR_BEETS 20 ;
    2020
    2121
  • coopr.pysp/stable/2.4/examples/pysp/farmer/farmer_lp.py

    r2734 r2938  
    3434model.PlantingCostPerAcre = Param(model.CROPS, within=PositiveReals)
    3535
    36 model.MeanYield = Param(model.CROPS, within=NonNegativeReals)
     36model.Yield = Param(model.CROPS, within=NonNegativeReals)
    3737
    3838#
     
    6060
    6161def cattle_feed_rule(i, model):
    62     return model.CattleFeedRequirement[i] <= (model.MeanYield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i]
     62    return model.CattleFeedRequirement[i] <= (model.Yield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i]
    6363
    6464model.EnforceCattleFeedRequirement = Constraint(model.CROPS, rule=cattle_feed_rule)
    6565
    6666def limit_amount_sold_rule(i, model):
    67     return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.MeanYield[i] * model.DevotedAcreage[i]) <= 0.0
     67    return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.Yield[i] * model.DevotedAcreage[i]) <= 0.0
    6868
    6969model.LimitAmountSold = Constraint(model.CROPS, rule=limit_amount_sold_rule)
  • coopr.pysp/stable/2.4/examples/pysp/farmer/maxmodels/ReferenceModel.py

    r2734 r2938  
    3434model.PlantingCostPerAcre = Param(model.CROPS, within=PositiveReals)
    3535
    36 model.MeanYield = Param(model.CROPS, within=NonNegativeReals)
     36model.Yield = Param(model.CROPS, within=NonNegativeReals)
    3737
    3838#
     
    6363
    6464def cattle_feed_rule(i, model):
    65     return model.CattleFeedRequirement[i] <= (model.MeanYield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i]   
     65    return model.CattleFeedRequirement[i] <= (model.Yield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i]   
    6666
    6767model.EnforceCattleFeedRequirement = Constraint(model.CROPS, rule=cattle_feed_rule)
    6868
    6969def limit_amount_sold_rule(i, model):
    70     return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.MeanYield[i] * model.DevotedAcreage[i]) <= 0.0
     70    return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.Yield[i] * model.DevotedAcreage[i]) <= 0.0
    7171
    7272model.LimitAmountSold = Constraint(model.CROPS, rule=limit_amount_sold_rule)
  • coopr.pysp/stable/2.4/examples/pysp/farmer/models/ReferenceModel.py

    r2734 r2938  
    3434model.PlantingCostPerAcre = Param(model.CROPS, within=PositiveReals)
    3535
    36 model.MeanYield = Param(model.CROPS, within=NonNegativeReals)
     36model.Yield = Param(model.CROPS, within=NonNegativeReals)
    3737
    3838#
     
    6060
    6161def cattle_feed_rule(i, model):
    62     return model.CattleFeedRequirement[i] <= (model.MeanYield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i]   
     62    return model.CattleFeedRequirement[i] <= (model.Yield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i]   
    6363
    6464model.EnforceCattleFeedRequirement = Constraint(model.CROPS, rule=cattle_feed_rule)
    6565
    6666def limit_amount_sold_rule(i, model):
    67     return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.MeanYield[i] * model.DevotedAcreage[i]) <= 0.0
     67    return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.Yield[i] * model.DevotedAcreage[i]) <= 0.0
    6868
    6969model.LimitAmountSold = Constraint(model.CROPS, rule=limit_amount_sold_rule)
  • coopr.pysp/stable/2.4/examples/pysp/farmer/nodedata/AboveAverageNode.dat

    r1446 r2938  
    11# above mean scenario
    22
    3 param MeanYield := WHEAT 3.0 CORN 3.6 SUGAR_BEETS 24 ;
     3param Yield := WHEAT 3.0 CORN 3.6 SUGAR_BEETS 24 ;
    44
    55
  • coopr.pysp/stable/2.4/examples/pysp/farmer/nodedata/AverageNode.dat

    r1446 r2938  
    11# "mean" scenario
    22
    3 param MeanYield := WHEAT 2.5 CORN 3 SUGAR_BEETS 20 ;
     3param Yield := WHEAT 2.5 CORN 3 SUGAR_BEETS 20 ;
    44
    55
  • coopr.pysp/stable/2.4/examples/pysp/farmer/nodedata/BelowAverageNode.dat

    r1446 r2938  
    11# below-mean scenario
    22
    3 param MeanYield := WHEAT 2.0 CORN 2.4 SUGAR_BEETS 16 ;
     3param Yield := WHEAT 2.0 CORN 2.4 SUGAR_BEETS 16 ;
    44
    55
  • coopr.pysp/stable/2.4/examples/pysp/farmer/nodedata/ReferenceModel.dat

    r1446 r2938  
    1717param PlantingCostPerAcre := WHEAT 150 CORN 230 SUGAR_BEETS 260 ;
    1818
    19 param MeanYield := WHEAT 3.0 CORN 3.6 SUGAR_BEETS 24 ;
     19param Yield := WHEAT 3.0 CORN 3.6 SUGAR_BEETS 24 ;
  • coopr.pysp/stable/2.4/examples/pysp/farmer/scenariodata/AboveAverageScenario.dat

    r1446 r2938  
    1919param PlantingCostPerAcre := WHEAT 150 CORN 230 SUGAR_BEETS 260 ;
    2020
    21 param MeanYield := WHEAT 3.0 CORN 3.6 SUGAR_BEETS 24 ;
     21param Yield := WHEAT 3.0 CORN 3.6 SUGAR_BEETS 24 ;
    2222
    2323
  • coopr.pysp/stable/2.4/examples/pysp/farmer/scenariodata/AverageScenario.dat

    r1446 r2938  
    1919param PlantingCostPerAcre := WHEAT 150 CORN 230 SUGAR_BEETS 260 ;
    2020
    21 param MeanYield := WHEAT 2.5 CORN 3 SUGAR_BEETS 20 ;
     21param Yield := WHEAT 2.5 CORN 3 SUGAR_BEETS 20 ;
    2222
    2323
  • coopr.pysp/stable/2.4/examples/pysp/farmer/scenariodata/BelowAverageScenario.dat

    r1446 r2938  
    1919param PlantingCostPerAcre := WHEAT 150 CORN 230 SUGAR_BEETS 260 ;
    2020
    21 param MeanYield := WHEAT 2.0 CORN 2.4 SUGAR_BEETS 16 ;
     21param Yield := WHEAT 2.0 CORN 2.4 SUGAR_BEETS 16 ;
    2222
    2323
  • coopr.pysp/stable/2.4/examples/pysp/farmer/scenariodata/ReferenceModel.dat

    r1446 r2938  
    1717param PlantingCostPerAcre := WHEAT 150 CORN 230 SUGAR_BEETS 260 ;
    1818
    19 param MeanYield := WHEAT 3.0 CORN 3.6 SUGAR_BEETS 24 ;
     19param Yield := WHEAT 3.0 CORN 3.6 SUGAR_BEETS 24 ;
  • coopr.pysp/stable/2.4/examples/pysp/sizes/SIZES10/phdriver.py

    r1706 r2938  
    1313
    1414# for profiling
    15 import cProfile
     15try:
     16    import cProfile as profile
     17except ImportError:
     18    import profile
    1619import pstats
    1720
     
    117120      traceback.print_exc()     
    118121else:
    119    cProfile.run('run_ph()','profile.stats')
     122   profile.run('run_ph()','profile.stats')
    120123   p=pstats.Stats('profile.stats')
    121124   p.sort_stats('time')
  • coopr.pysp/stable/2.4/scripts/computeconf

    r2471 r2938  
    1313import os
    1414import random
     15import math
    1516from coopr.pysp.scenariotree import *
    1617from coopr.pysp.phinit import *
     
    2021def run(args=None):
    2122
    22 # JPW: args=None is OK, as the arguments are by default propagated through sys.argv
    23 #   if args == None:
    24 #      print "Error: testconf run called with no args."
    25 #      sys.exit(1)
    2623   try:
    27       conf_options_parser = construct_ph_options_parser("testconf [options]")
     24      conf_options_parser = construct_ph_options_parser("computeconf [options]")
    2825      conf_options_parser.add_option("--Fraction-of-Scenarios-for-Solve",
    2926                                     help="The fraction of scenarios that are allocated to finding a solution. Default is 0.5",
     
    4441                                     type="float",
    4542                                     default=0.10)
    46       conf_options_parser.add_option("--solve-hatx-with-ef-only",
    47                                      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",
    4845                                     action="store_true",
    49                                      dest="solve_hatx_with_ef_only",
    50                                      default=False)     
     46                                     dest="solve_xhat_with_ef_only",
     47                                     default=False)
     48      conf_options_parser.add_option("--random-seed",
     49                                     help="Seed the random number generator used to select samples. Defaults to 0, indicating time seed will be used.",
     50                                     action="store",
     51                                     dest="random_seed",
     52                                     type="int",
     53                                     default=0)
     54      conf_options_parser.add_option("--append-file",
     55                                     help="File to which summary run information is appended, for output tracking purposes.",
     56                                     action="store",
     57                                     dest="append_file",
     58                                     type="string",
     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)           
     77
    5178      (options, args) = conf_options_parser.parse_args(args=args)
    5279   except SystemExit:
     
    5582      return
    5683
    57    # HEY!!!!! Look at this!!!!
    58    random.seed(17)
    59 
    60    reference_model, reference_instance, full_scenario_tree, si = load_reference_and_scenario_models(options)
    61 
    62    # HERE: If we get an objective function with maximization, bail right away.
    63 
    64    print "Confidence Interval Calculations"
    65    print "full_scenario_tree pprint:"
    66    full_scenario_tree.pprint()
    67    scenariocount = len(full_scenario_tree._stages[-1]._tree_nodes)
     84   # seed the generator is a user-supplied seed is provided. otherwise,
     85   # python will seed from the current system time.
     86   if options.random_seed > 0:
     87      random.seed(options.random_seed)
     88
     89   # create the reference instances and the scenario tree - no scenario instances yet.
     90   print "Loading reference model and scenario tree"
     91   reference_model, reference_instance, full_scenario_tree, scenario_tree_instance = load_reference_and_scenario_models(options)
     92   if (reference_model is None) or (reference_instance is None) or (full_scenario_tree is None) or (scenario_tree_instance is None):
     93      print "***ERROR: Failed to initialize reference model/instance pair and/or the scenario tree."
     94      sys.exit(1)
     95
     96   # verify that we're dealing with a minimization problem - the script might be easily made to work with
     97   # maximization problems - we just haven't done the work to do that.
     98   objective_name = reference_model.active_components()[Objective].keys()[0]
     99   objective = reference_model.active_components()[Objective][objective_name]
     100   if objective.sense != minimize:
     101      print "***ERROR: Confidence intervals are available only for minimization problems"
     102      sys.exit(1)
     103
     104   print "Starting confidence interval calculation..."
     105   
     106   scenario_count = len(full_scenario_tree._stages[-1]._tree_nodes)
    68107   if len(full_scenario_tree._stages) > 2:
    69       print "Confidence intervals are available only for two stage problems. Stage count=",len(full_scenario_tree._stages)
     108      print "***ERROR: Confidence intervals are available only for two stage stochastic programs;"+str(len(full_scenario_tree._stages))+" stages specified"
    70109      sys.exit(1)
    71110
    72    # random permutation of the indexes
    73    IndList = range(scenariocount)
    74    random.shuffle(IndList)
    75    print "After shuffle, IndList=",IndList
    76    
    77    num_scenarios_for_solution = int(options.fraction_for_solve * scenariocount)
     111   # randomly permute the indices to extract a subset to compute xhat.
     112   index_list = range(scenario_count)
     113   random.shuffle(index_list)
     114   if options.verbose is True:
     115      print "Random permutation of the scenario indices="+str(index_list)
     116   
     117   # figure out the scenario counts for both the xhat and confidence interval computations.
     118   num_scenarios_for_solution = int(options.fraction_for_solve * scenario_count)
    78119   n_g = options.n_g
    79    num_scenarios_per_sample = int((scenariocount - num_scenarios_for_solution) / n_g)  #n in Morton's slides
    80    wastedscenarios = scenariocount - num_scenarios_for_solution - n_g * num_scenarios_per_sample
    81    
    82    print scenariocount, "Scenarios,from which ",num_scenarios_for_solution, " are to be used to find a solution and "
    83    print n_g," groups of ",num_scenarios_per_sample," are to be used for the confidence interval."
    84    if wastedscenarios > 0:
    85       print "(so ",wastedscenarios," will not be used.)"
    86 
    87    # create ph for finding the solution
    88    hatxph = ph_for_bundle(0, num_scenarios_for_solution, reference_model, full_scenario_tree, reference_instance, si, IndList, options)
    89 
    90    if options.solve_hatx_with_ef_only is True:
    91       # DLW: LOOK HERE!
    92       hatex_ef = create_ef_instance(hatxph._scenario_tree, hatxph._instances)
    93       ef_results = write_and_solve_ef(hatex_ef, hatxph._instances, options)
    94       load_ef_solution(ef_results, hatex_ef, hatxph._instances)
    95       hatxph._scenario_tree.snapshotSolutionFromInstances(hatxph._instances)
    96       hatxph._scenario_tree.pprintSolution()
     120   num_scenarios_per_sample = int((scenario_count - num_scenarios_for_solution) / n_g) # 'n' in Morton's slides
     121   wasted_scenarios = scenario_count - num_scenarios_for_solution - n_g * num_scenarios_per_sample
     122   
     123   print "Problem contains "+str(scenario_count)+" scenarios, of which "+str(num_scenarios_for_solution)+" will be used to find a solution xhat."
     124   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."
     125   if wasted_scenarios > 0:
     126      print "A total of "+str(wasted_scenarios)+" scenarios will not be used."
     127
     128   # create a ph object for finding the solution. we do this even if we're solving the extensive form
     129   # directly, mainly out of convenience - we're leveraging the code in ph_for_bundle to create the
     130   # scenario tree and scenario instances.
     131   print "Loading scenario instances and initializing scenario tree for xhat scenario bundle."
     132   xhat_ph = ph_for_bundle(0, num_scenarios_for_solution, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options)
     133   xhat_obj = None
     134
     135   if options.solve_xhat_with_ef_only is True:
     136      print "Creating the xhat extensive form."
     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)                                   
     142      print "Solving the xhat extensive form."
     143      ef_results = write_and_solve_ef(hatex_ef, xhat_ph._instances, options)
     144      load_ef_solution(ef_results, hatex_ef, xhat_ph._instances)
     145      # IMPT: the following method populates the _solution variables on the scenario tree
     146      #       nodes by forming an average of the corresponding variable values for all
     147      #       instances particpating in that node. if you don't do this, the scenario tree
     148      #       doesn't have the solution - and we need this below for variable fixing.
     149      xhat_ph._scenario_tree.snapshotSolutionFromInstances(xhat_ph._instances)
     150      xhat_obj = float(ef_results.solution(0).objective['f'].value)  ## DLW to JPW: we need the gap too
     151      print "Extensive form objective value given xhat="+str(xhat_obj)
    97152   else:
    98       print "SOLVING HATX VIA PH!"
    99       hatxph.solve()
    100    print "so we now have the hat{x} variables in the PHAVG structure of hatxph"
    101 
    102    print "now form and solve the problems for each sample"
     153      print "Solving for xhat via Progressive Hedging."     
     154      xhat_ph.solve()
     155      # TBD - grab xhat_obj; complicated by the fact that PH may not have converged.
     156      # TBD - also note sure if PH calls snapshotSolutionFromAverages.
     157   print "The computation of xhat is complete - starting to compute confidence interval via sub-sampling."
     158
    103159   # in order to handle the case of scenarios that are not equally likely, we will split the expectations for Gsupk
    104160   # 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
    105    G_supk_of_hatx = [] # really not always needed... http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
    106    Gbar = 0
     161   g_supk_of_xhat = [] # really not always needed... http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
     162   g_bar = 0
     163   sum_xstar_obj_given_xhat = 0
     164   
    107165   for k in range(1, n_g+1):
    108       start = num_scenarios_for_solution + (k-1)*num_scenarios_per_sample
    109       stop = start + num_scenarios_per_sample - 1
    110       print "Sample k=", k
    111       # NOTE: We'll never run this ph - code could be refactored
    112       phforGk = ph_for_bundle(start, stop, reference_model, full_scenario_tree, reference_instance, si, IndList, options)
    113       efforGk = create_ef_instance(phforGk._scenario_tree, phforGk._instances)     
    114       ef_results = write_and_solve_ef(efforGk, phforGk._instances, options)
    115       load_ef_solution(ef_results, efforGk, phforGk._instances)
    116       print "HEY JPW!!! we need to add or subtract the gap"
     166     
     167      start_index = num_scenarios_for_solution + (k-1)*num_scenarios_per_sample
     168      stop_index = start_index + num_scenarios_per_sample - 1
     169      print "Computing statistics for sample k="+str(k)+"."
     170      if options.verbose is True:
     171         print "Bundle start index="+str(start_index)+", stop index="+str(stop_index)+"."
     172     
     173      # compute this xstar solution for the EF associated with sample k.
     174      print "Loading scenario instances and initializing scenario tree for xstar scenario bundle."     
     175      gk_ph = ph_for_bundle(start_index, stop_index, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options)
     176      print "Creating the xstar extensive form."     
     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)
     182      print "Solving the xstar extensive form."     
     183      ef_results = write_and_solve_ef(gk_ef, gk_ph._instances, options)
     184      load_ef_solution(ef_results, gk_ef, gk_ph._instances)
     185
     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.
     188      gk_ph._scenario_tree.snapshotSolutionFromInstances(gk_ph._instances)
     189     
     190      # extract the objective function value corresponding to the xstar solution, along with any gap information.
    117191      # JPW doesn't like the 'f' below - seems bad hard-coding a literal.
    118       E_f_of_xstar = float(ef_results.solution(0).objective['f'].value)  ## DLW to JPW: we need the gap too
    119       print "E_f_of_xstar=",E_f_of_xstar
    120       efforgK_gap = ef_results.solution(0).gap # assuming this is the absolute gap
    121       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?
    122 
    123       # to get fj(hatx) we need the obj value for each scenario, j, in the bundle evaluated at hat x
    124       # this is no small thing. A lot of code here is copied from ph.py
    125       action_handles = []
    126       scenario_action_handle_map = {} # maps scenario names to action handles
    127       action_handle_scenario_map = {} # maps action handles to scenario names
    128       E_f_of_Hatx = 0
    129       # to support parallel, loop to launch then loop to get and use results
    130       for scenario in  phforGk._scenario_tree._scenarios:
    131          instance = phforGk._instances[scenario._name]
    132          # do the fixing at hatx
    133          fix_first_stage_vars_for_instance_from_PHAVG(hatxph, instance)
    134 
     192      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.
     193      xstar_obj_gap = ef_results.solution(0).gap # assuming this is the absolute gap
     194      xstar_obj_bound = xstar_obj - xstar_obj_gap # HERE - watch for signs on gap - is it absolute, or can it be negative?     
     195      print "Extensive form objective value given xstar="+str(xstar_obj)+"; gap="+str(xstar_obj_gap)
     196
     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."
     201      for scenario in  gk_ph._scenario_tree._scenarios:
     202         instance = gk_ph._instances[scenario._name]
     203         fix_first_stage_vars(xhat_ph._scenario_tree, instance)
    135204         instance.preprocess()
    136 
    137          new_action_handle = phforGk._solver_manager.queue(instance, opt=phforGk._solver, tee=phforGk._output_solver_log)
    138          scenario_action_handle_map[scenario._name] = new_action_handle
    139          action_handle_scenario_map[new_action_handle] = scenario._name
    140 
    141          action_handles.append(new_action_handle)
    142 
    143       # loop to get and use results
    144       num_results_so_far = 0
    145       while (num_results_so_far < len(phforGk._scenario_tree._scenarios)):
    146 
    147          action_handle = phforGk._solver_manager.wait_any()
    148          results = phforGk._solver_manager.get_results(action_handle)         
    149          scenario_name = action_handle_scenario_map[action_handle]
    150          instance = phforGk._instances[scenario_name]         
    151 
    152          if phforGk._verbose is True:
    153             print "Sampling Results obtained for scenario="+scenario_name
    154 
    155          if len(results.solution) == 0:
    156             results.write(num=1)
    157             raise RuntimeError, "Sampling Solve failed for scenario="+scenario_name+"; no solutions generated"
    158 
    159          if phforGk._output_solver_results is True:
    160             print "Sampling Results for scenario=",scenario_name
    161             results.write(num=1)
    162 
    163          start_time = time.time()
    164          instance.load(results)
    165          end_time = time.time()
    166          if phforGk._output_times is True:
    167             print "Time loading results into instance="+str(end_time-start_time)+" seconds"
    168 
    169          if phforGk._verbose is True:                 
    170             print "Successfully loaded solution for scenario="+scenario_name
    171 
    172          num_results_so_far = num_results_so_far + 1
    173 
    174 
    175          objval = phforGk._scenario_tree.compute_scenario_cost(instance)
    176          E_f_of_Hatx += objval * scenario._probability
    177 
    178       G_supk_of_hatx.append(E_f_of_Hatx - E_f_of_xstar_bound)   
    179       Gbar += G_supk_of_hatx[k-1]
    180    Gbar = Gbar / n_g
    181    # second pass for variance calculation (because we like storing the G_supk)
    182    VarG = 0
     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)
     211      print "Sample extensive form objective value given xhat="+str(xstar_obj_given_xhat)
     212
     213      g_supk_of_xhat.append(xstar_obj_given_xhat - xstar_obj_bound)   
     214      g_bar += g_supk_of_xhat[k-1]
     215      sum_xstar_obj_given_xhat += xstar_obj_given_xhat
     216     
     217   g_bar = g_bar / n_g
     218   # second pass for variance calculation (because we like storing the g_supk)
     219   g_var = 0
    183220   for k in range(0, n_g-1):
    184       VarG = VarG + (G_supk_of_hatx[k] - Gbar) * (G_supk_of_hatx[k] - Gbar)
    185    VarG = VarG / (n_g - 1)    # sample var
    186    print "Gbar=",Gbar
    187    print "VarG=",VarG     
    188 
    189 #==============================================   
    190 def ph_for_bundle(BundleStart, BundleStop, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, IndList, options):
     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
     223   print "g_bar=",g_bar
     224   print "g_stddev=",math.sqrt(g_var)
     225   print "Average f(xhat)=",sum_xstar_obj_given_xhat / n_g
     226
     227   if options.append_file is not None:
     228      output_file = open(options.append_file, "a")
     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,
     230      output_file.close()
     231
     232#
     233# routine to create a down-sampled (bundled) scenario tree and the associated PH object.
     234#
     235def ph_for_bundle(bundle_start, bundle_stop, reference_model, full_scenario_tree, reference_instance, scenario_tree_instance, index_list, options):
    191236   
    192237   scenarios_to_bundle = []
    193    for i in range(BundleStart, BundleStop+1):   # python has zero based indexes
    194       scenarios_to_bundle.append(full_scenario_tree._scenarios[IndList[i]]._name)
    195 
    196    print "Scenarios to bundle:"+str(scenarios_to_bundle)
     238   for i in range(bundle_start, bundle_stop+1):
     239      scenarios_to_bundle.append(full_scenario_tree._scenarios[index_list[i]]._name)
     240
     241   if options.verbose is True:
     242      print "Creating PH object for scenario bundle="+str(scenarios_to_bundle)
    197243
    198244   scenario_tree_for_soln = ScenarioTree(scenarioinstance=reference_instance,
    199                             scenariotreeinstance=scenario_tree_instance,
    200                             scenariobundlelist=scenarios_to_bundle)
    201 
     245                                         scenariotreeinstance=scenario_tree_instance,
     246                                         scenariobundlelist=scenarios_to_bundle)
    202247   if scenario_tree_for_soln.validate() is False:
    203       print "***ERROR: Bundle Scenario tree is invalid****"
    204       sys.exit(0)
    205    else:
    206       print "Scenario tree for solution is valid!"
     248      print "***ERROR: Bundled scenario tree is invalid!!!"
     249      sys.exit(1)
    207250
    208251   ph = create_ph_from_scratch(options, reference_model, reference_instance, scenario_tree_for_soln)
    209252   return ph
    210253
    211 #==============================================   
    212 def fix_first_stage_vars_for_instance_from_PHAVG(ph, instance):
    213 # use the values from PHAVG of the ph object to fix the first stage vars in the instance
    214    stage = ph._scenario_tree._stages[0]   # root node only
    215    # use the first instance from PH because all should have the same PHAVG
    216    phinstance = ph._instances[ph._scenario_tree._scenarios[0]._name]
     254#
     255# fixes the first stage variables in the input instance to the values
     256# of the solution indicated in the input scenario tree.
     257#
     258def fix_first_stage_vars(scenario_tree, instance):
     259
     260   stage = scenario_tree._stages[0]
     261   root_node = stage._tree_nodes[0] # there should be only one root node!
    217262
    218263   for (variable, index_template, variable_indices) in stage._variables:
    219 
    220        variable_name = variable.name
    221        variable_type = variable.domain
    222        # HERE: we want to not take anything from PHAVG, but rather from the scenario tree itself - the
    223        #       only difficulty is finding the root node object.
    224        avg_parameter_name = "PHAVG_"+variable_name
    225        avg_parameter = getattr(phinstance, avg_parameter_name)
    226 
    227        for index in variable_indices:
    228           if getattr(phinstance,variable_name)[index].status != VarStatus.unused:
    229              fix_value = avg_parameter[index]()
    230              if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
    231                 fix_value = int(round(fix_value))
    232              getattr(instance,variable.name)[index].fixed = True
    233              getattr(instance,variable.name)[index].value = fix_value
    234    print "DLW says: first stage vars are fixed; maybe we need to delete any constraints with only first stage vars due to precision issues"
     264     
     265      variable_name = variable.name
     266      variable_type = variable.domain
     267
     268      scenario_tree_var = root_node._solutions[variable_name]
     269
     270      for index in variable_indices:
     271         if getattr(instance, variable_name)[index].status != VarStatus.unused:
     272            fix_value = scenario_tree_var[index]()
     273            if isinstance(variable_type, IntegerSet) or isinstance(variable_type, BooleanSet):
     274               fix_value = int(round(fix_value))
     275            getattr(instance,variable.name)[index].fixed = True
     276            getattr(instance,variable.name)[index].value = fix_value
     277
     278#   print "DLW says: first stage vars are fixed; maybe we need to delete any constraints with only first stage vars due to precision issues"
    235279
    236280#==============================================   
    237281def write_and_solve_ef(master_instance, scenario_instances, options):
    238282
    239    ef_file_name = "dlwef.lp"
    240    write_ef(master_instance, scenario_instances, os.path.expanduser(ef_file_name))
    241 
    242    print ""
    243    print "Solving extensive form written to file="+os.path.expanduser(ef_file_name)
    244    print ""
     283   ef_file_name = os.path.expanduser("dlwef.lp")
     284   write_ef(master_instance, scenario_instances, ef_file_name)
    245285
    246286   ef_solver = SolverFactory(options.solver_type)
     
    260300      raise ValueError, "Failed to create solver manager of type="+options.solver_type+" for use in extensive form solve"
    261301
    262    print "Queuing extensive form solve"
    263    ef_action_handle = ef_solver_manager.queue(os.path.expanduser(ef_file_name), opt=ef_solver, warmstart=False, tee=options.output_ef_solver_log)
    264    print "Waiting for extensive form solve"
     302   ef_action_handle = ef_solver_manager.queue(ef_file_name, opt=ef_solver, warmstart=False, tee=options.output_ef_solver_log)
    265303   ef_results = ef_solver_manager.wait_for(ef_action_handle)
    266 #   print "Extensive form solve results:"
    267 #   ef_results.write(num=1)
     304   os.remove(ef_file_name)
     305
    268306   return ef_results
    269307
    270 #=================================================================================================
    271 # JPW: Any executable code can be placed at the end of a file, and it will be executed.
    272 print "STARTING TO RUN"
     308
     309#
     310# the main script routine starts here - and is obviously pretty simple!
     311#
     312
    273313run()
    274 print "DONE"
Note: See TracChangeset for help on using the changeset viewer.