Changeset 2434

Timestamp:

Mar 12, 2010 3:43:03 AM (9 years ago)

Author:

wehart

Message:

Merged revisions 2391-2433 via svnmerge from
​https://software.sandia.gov/svn/public/coopr/coopr.pysp/trunk

........

r2393 | jwatson | 2010-02-22 17:18:12 -0700 (Mon, 22 Feb 2010) | 3 lines

Single-threading the ph solver server (for obvious reasons - don't want to spawn a zillion solves) and fixing issues in the solver server test client. Solves occur and results can be shipped back to a client, w/o passing of the instance.

........

r2396 | jwatson | 2010-02-23 15:18:15 -0700 (Tue, 23 Feb 2010) | 3 lines

Added initial cut at a PH solver manager - only works for iteration 0 at the moment.

........

r2398 | jwatson | 2010-02-24 11:34:47 -0700 (Wed, 24 Feb 2010) | 3 lines

Re-factoring the PH instance augmentation routines, moving them from ph.py to phutils.py so they can be used by the PH solver server.

........

r2401 | jwatson | 2010-02-24 15:08:49 -0700 (Wed, 24 Feb 2010) | 3 lines

PH solver manager and servers now handle rho and weight/average update requests correctly.

........

r2402 | jwatson | 2010-02-25 09:33:25 -0700 (Thu, 25 Feb 2010) | 1 line

Added --all-scenarios option to PH solver server

........

r2403 | jwatson | 2010-02-25 09:42:02 -0700 (Thu, 25 Feb 2010) | 3 lines

Added a utility method to compute the aggregate cost for a scenario, by looping over the stage variables and querying instance variable values.

........

r2404 | jwatson | 2010-02-25 09:58:33 -0700 (Thu, 25 Feb 2010) | 3 lines

Eliminated debug print statement from scenario tree class. Used new compute_scenario_cost method of scenario tree class to compute non-PH-weighted objective in PH output (per-iteration).

........

r2405 | jwatson | 2010-02-25 19:54:36 -0700 (Thu, 25 Feb 2010) | 3 lines

Restructuring of PySP to facilitate full implementation of PH solver servers (penalty objectives in particular).

........

r2406 | jwatson | 2010-02-26 20:54:23 -0700 (Fri, 26 Feb 2010) | 3 lines

More tweaks to the PH solver manager/servers.

........

r2410 | jwatson | 2010-03-02 14:53:30 -0700 (Tue, 02 Mar 2010) | 5 lines

Fixed the PH solver server to automatically deregister a delegator object if it already exists in the name server.

Fixed various bugs in the PH solver server relating to quadratic expression output and weight/rho updates. It now replicates CPLEX results across the server network.

........

r2411 | jwatson | 2010-03-04 16:23:41 -0700 (Thu, 04 Mar 2010) | 1 line

Initial commit of PySP lot sizing example

........

r2412 | jwatson | 2010-03-04 16:24:36 -0700 (Thu, 04 Mar 2010) | 1 line

Miscellaneous EF fixes - mainly relating to output messages

........

r2413 | jwatson | 2010-03-08 10:20:36 -0700 (Mon, 08 Mar 2010) | 1 line

PySP EF restructuring / refactoring

........

r2414 | jwatson | 2010-03-08 20:18:14 -0700 (Mon, 08 Mar 2010) | 6 lines

Various EF writer changes in PySP, including:
1) Added a --verbose option, and by default am disabling echo of the scenario tree.
2) Changed the default behavior, due to recent experience with lot sizing and wind farm examples, to enable garbage collection.
3) Added a --disable-gc option to disable garbage collection, set to False by default.

........

r2418 | jwatson | 2010-03-11 15:49:52 -0700 (Thu, 11 Mar 2010) | 3 lines

More restructuring and simplification of the PySP extensive writer code, which is now rather compact.

........

r2423 | wehart | 2010-03-11 16:01:10 -0700 (Thu, 11 Mar 2010) | 3 lines

Rework of unit tests to (a) import pyutilib.th as 'unittest' and
(b) employ test skipping.

........

r2427 | jwatson | 2010-03-11 20:06:30 -0700 (Thu, 11 Mar 2010) | 3 lines

Added logic to automatically disable and re-enable garbage collection at various sane points in the PySP extensive form writer, in order to (noticeably) improve performance.

........

r2428 | jwatson | 2010-03-11 20:25:39 -0700 (Thu, 11 Mar 2010) | 3 lines

Some "final" performance improvements to the PySP extensive form writer. The lion's share of the remaining time is consumed by instance construction.

........

Location:

coopr.pysp/stable/2.3

Files:

• . (modified) (1 prop)
• coopr/pysp/__init__.py (modified) (1 diff)
• coopr/pysp/ef.py (modified) (13 diffs)
• coopr/pysp/ef_writer_script.py (modified) (4 diffs)
• coopr/pysp/ph.py (modified) (10 diffs)
• coopr/pysp/phinit.py (modified) (2 diffs)
• coopr/pysp/phobjective.py (copied) (copied from coopr.pysp/trunk/coopr/pysp/phobjective.py)
• coopr/pysp/phserver.py (modified) (16 diffs)
• coopr/pysp/phsolvermanager.py (copied) (copied from coopr.pysp/trunk/coopr/pysp/phsolvermanager.py)
• coopr/pysp/phutils.py (modified) (1 diff)
• coopr/pysp/scenariotree.py (modified) (11 diffs)
• coopr/pysp/tests/unit/test_ph.py (modified) (7 diffs)
• examples/pysp/lotsizing (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing)
• examples/pysp/lotsizing/README.txt (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/README.txt)
• examples/pysp/lotsizing/generator (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/generator)
• examples/pysp/lotsizing/generator/generate_instance.py (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/generator/generate_instance.py)
• examples/pysp/lotsizing/models (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/models)
• examples/pysp/lotsizing/models/ReferenceModel.py (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/models/ReferenceModel.py)
• examples/pysp/lotsizing/testinstance (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/testinstance)
• examples/pysp/lotsizing/testinstance/Node_1_1.dat (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/testinstance/Node_1_1.dat)
• examples/pysp/lotsizing/testinstance/Node_2_1.dat (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/testinstance/Node_2_1.dat)
• examples/pysp/lotsizing/testinstance/Node_2_2.dat (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/testinstance/Node_2_2.dat)
• examples/pysp/lotsizing/testinstance/Node_2_3.dat (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/testinstance/Node_2_3.dat)
• examples/pysp/lotsizing/testinstance/ReferenceModel.dat (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/testinstance/ReferenceModel.dat)
• examples/pysp/lotsizing/testinstance/ScenarioStructure.dat (copied) (copied from coopr.pysp/trunk/examples/pysp/lotsizing/testinstance/ScenarioStructure.dat)
• scripts/ph_test_client (modified) (2 diffs)

coopr.pysp/stable/2.3/coopr/pysp/__init__.py

@@ -21 +21 @@
 from ef_writer_script import *
 from phinit import *
+from phsolvermanager import *
+from phobjective import *
 
 pyutilib.component.core.PluginGlobals.pop_env()
-
 
 try:

coopr.pysp/stable/2.3/coopr/pysp/ef.py

@@ -5 +5 @@
 import traceback
 import copy
+import gc
 
 from coopr.pysp.scenariotree import *
@@ -15 +16 @@
 from coopr.pyomo.base.var import _VarValue, _VarBase
 
+#
 # brain-dead utility for determing if there is a binary to write in the 
 # composite model - need to know this, because CPLEX doesn't like empty
 # binary blocks in the LP file.
+#
+
 def binaries_present(master_model, scenario_instances):
 
@@ -34 +38 @@
    return False
 
+#
 # brain-dead utility for determing if there is a binary to write in the 
 # composite model - need to know this, because CPLEX doesn't like empty
 # integer blocks in the LP file.
+#
+
 def integers_present(master_model, scenario_instances):
 
@@ -53 +60 @@
    return False
 
+#
+# a routine to create the extensive form, given an input scenario tree and instances.
+# IMPT: unlike scenario instances, the extensive form instance is *not* self-contained.
+#       in particular, it has binding constraints that cross the binding instance and
+#       the scenario instances. it is up to the caller to keep track of which scenario
+#       instances are associated with the extensive form. this might be something we
+#       encapsulate at some later time.
+#
+
+def create_ef_instance(scenario_tree, scenario_instances,
+                       verbose_output = False,
+                       generate_weighted_cvar = False, cvar_weight = None, risk_alpha = None):
+
+   #
+   # validate cvar options, if specified.
+   #
+   if generate_weighted_cvar is True:
+      if (cvar_weight is None) or (cvar_weight <= 0.0):
+         raise RuntimeError, "Weight of CVaR term must be >= 0.0 - value supplied="+str(cvar_weight)
+      if (risk_alpha is None) or (risk_alpha <= 0.0) or (risk_alpha >= 1.0):
+         raise RuntimeError, "CVaR risk alpha must be between 0 and 1, exclusive - value supplied="+str(risk_alpha)
+
+      if verbose_output is True:
+         print "Writing CVaR weighted objective"
+         print "CVaR term weight="+str(cvar_weight)
+         print "CVaR alpha="+str(risk_alpha)
+         print ""
+
+   #
+   # create the new and empty binding instance.
+   #
+
+   binding_instance = Model()
+   binding_instance.name = "MASTER"
+
+   # walk the scenario tree - create variables representing the common values for all scenarios
+   # associated with that node. the constraints will be created later. also create expected-cost
+   # variables for each node, to be computed via constraints/objectives defined in a subsequent pass.
+   # master variables are created for all nodes but those in the last stage. expected cost variables
+   # are, for no particularly good reason other than easy coding, created for nodes in all stages.
+   if verbose_output is True:
+      print "Creating variables for master binding instance"
+
+   for stage in scenario_tree._stages:
+
+      # first loop is to create master (blended) variables across all stages but the last.
+      for (stage_variable, index_template, stage_variable_indices) in stage._variables:
+
+         if verbose_output is True:
+            print "Creating master variable and blending constraints for decision variable="+stage_variable.name+", indices="+str(index_template)
+
+         for tree_node in stage._tree_nodes:
+
+            if stage != scenario_tree._stages[-1]:      
+
+               master_variable_name = tree_node._name + "_" + stage_variable.name
+
+               # because there may be a single stage variable and multiple indices, check
+               # for the existence of the variable at this node - if you don't, you'll
+               # inadvertently over-write what was there previously!
+               master_variable = None
+               try:
+                  master_variable = getattr(binding_instance, master_variable_name)
+               except:
+                  new_master_variable_index = stage_variable._index
+                  new_master_variable = None
+                  if (len(new_master_variable_index) is 1) and (None in new_master_variable_index):
+                     new_master_variable = Var(name=stage_variable.name)
+                  else:
+                     new_master_variable = Var(new_master_variable_index, name=stage_variable.name)
+                  new_master_variable.construct()
+                  new_master_variable._model = binding_instance
+                  setattr(binding_instance, master_variable_name, new_master_variable)
+
+                  master_variable = new_master_variable
+
+               for index in stage_variable_indices:
+
+                  is_used = True # until proven otherwise                     
+                  for scenario in tree_node._scenarios:
+                     instance = scenario_instances[scenario._name]
+                     if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
+                        is_used = False
+
+                  is_fixed = False # until proven otherwise
+                  for scenario in tree_node._scenarios:
+                     instance = scenario_instances[scenario._name]
+                     if getattr(instance,stage_variable.name)[index].fixed is True:
+                        is_fixed = True
+
+                  if (is_used is True) and (is_fixed is False):
+                           
+                     # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
+                     # and because presolve/simplification is name-based, the names *have* to be different.
+                     master_variable[index].var = master_variable
+                     master_variable[index].name = tree_node._name + "_" + master_variable[index].name
+
+                     for scenario in tree_node._scenarios:
+
+                        scenario_instance = scenario_instances[scenario._name]
+                        scenario_variable = getattr(scenario_instance, stage_variable.name)
+                        new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
+                        new_constraint = Constraint(name=new_constraint_name)
+                        new_expr = master_variable[index] - scenario_variable[index]
+                        new_constraint.add(None, (0.0, new_expr, 0.0))
+                        new_constraint._model = binding_instance
+                        setattr(binding_instance, new_constraint_name, new_constraint)
+
+      # the second loop is for creating the stage cost variable in each tree node.
+      for tree_node in stage._tree_nodes:                        
+
+         # create a variable to represent the expected cost at this node -
+         # the constraint to compute this comes later.
+         expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
+         expected_cost_variable = Var(name=expected_cost_variable_name)
+         expected_cost_variable._model = binding_instance
+         setattr(binding_instance, expected_cost_variable_name, expected_cost_variable)
+
+   # if we're generating the weighted CVaR objective term, create the
+   # corresponding variable and the master CVaR eta variable.
+   if generate_weighted_cvar is True:
+
+      root_node = scenario_tree._stages[0]._tree_nodes[0]
+      
+      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
+      cvar_cost_variable = Var(name=cvar_cost_variable_name)
+      setattr(binding_instance, cvar_cost_variable_name, cvar_cost_variable)
+      cvar_cost_variable.construct()
+
+      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
+      cvar_eta_variable = Var(name=cvar_eta_variable_name)
+      setattr(binding_instance, cvar_eta_variable_name, cvar_eta_variable)      
+      cvar_eta_variable.construct()
+
+   binding_instance.preprocess()
+
+   # ditto above for the (non-expected) cost variable.
+   for stage in scenario_tree._stages:
+
+      (cost_variable,cost_variable_index) = stage._cost_variable
+
+      print "Creating master variable and blending constraints for cost variable="+cost_variable.name+", index="+str(cost_variable_index)
+
+      for tree_node in stage._tree_nodes:
+
+         # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)         
+
+         # this is undoubtedly wasteful, in that a cost variable
+         # for each tree node is created with *all* indices.          
+         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
+         new_cost_variable_index = cost_variable._index
+         new_cost_variable = None
+         if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index):
+            new_cost_variable = Var(name=new_cost_variable_name)
+         else:
+            new_cost_variable = Var(new_cost_variable_index, name=new_cost_variable_name)
+         new_cost_variable.construct()
+         new_cost_variable._model = binding_instance
+         setattr(binding_instance, new_cost_variable_name, new_cost_variable)         
+
+         # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
+         new_cost_variable[cost_variable_index].var = new_cost_variable
+         if cost_variable_index is not None:
+            # if the variable index is None, the variable is derived from a VarValue, so the
+            # name gets updated automagically.
+            new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
+
+         for scenario in tree_node._scenarios:
+
+            scenario_instance = scenario_instances[scenario._name]
+            scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
+            new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
+            new_constraint = Constraint(name=new_constraint_name)
+            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
+            new_constraint.add(None, (0.0, new_expr, 0.0))
+            new_constraint._model = binding_instance
+            setattr(binding_instance, new_constraint_name, new_constraint)
+
+   # create the constraints for computing the master per-node cost variables,
+   # i.e., the current node cost and the expected cost of the child nodes.
+   # if the root, then the constraint is just the objective.
+
+   for stage in scenario_tree._stages:
+
+      (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
+
+      for tree_node in stage._tree_nodes:
+
+         node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
+         node_expected_cost_variable = getattr(binding_instance, node_expected_cost_variable_name)
+
+         node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
+         node_cost_variable = getattr(binding_instance, node_cost_variable_name)                        
+            
+         constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
+
+         for child_node in tree_node._children:
+
+            child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
+            child_node_expected_cost_variable = getattr(binding_instance, child_node_expected_cost_variable_name)
+            constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
+
+         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
+         new_constraint = Constraint(name=new_constraint_name)
+         new_constraint.add(None, (0.0, constraint_expr, 0.0))
+         new_constraint._model = binding_instance                     
+         setattr(binding_instance, new_constraint_name, new_constraint)
+
+         if tree_node._parent is None:
+
+            an_instance = scenario_instances[scenario_instances.keys()[0]]
+            an_objective = an_instance.active_components(Objective)
+            opt_sense = an_objective[an_objective.keys()[0]].sense
+
+            opt_expression = node_expected_cost_variable
+
+            if generate_weighted_cvar is True:
+               cvar_cost_variable_name = "CVAR_COST_" + tree_node._name
+               cvar_cost_variable = getattr(binding_instance, cvar_cost_variable_name)
+               opt_expression += cvar_weight * cvar_cost_variable            
+
+            new_objective = Objective(name="MASTER", sense=opt_sense)
+            new_objective._data[None].expr = node_expected_cost_variable
+            setattr(binding_instance, "MASTER", new_objective)
+
+   # CVaR requires the addition of a variable per scenario to represent the cost excess,
+   # and a constraint to compute the cost excess relative to eta. we also replicate (following
+   # what we do for node cost variables) an eta variable for each scenario instance, and
+   # require equality with the master eta variable via constraints.
+   if generate_weighted_cvar is True:
+      
+      root_node = scenario_tree._stages[0]._tree_nodes[0]
+
+      master_cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
+      master_cvar_eta_variable = getattr(binding_instance, master_cvar_eta_variable_name)
+      
+      for scenario_name in scenario_instances.keys():
+         scenario_instance = scenario_instances[scenario_name]
+
+         # unique names are required because the presolve isn't
+         # aware of the "owning" models for variables.
+         cvar_excess_variable_name = "CVAR_EXCESS_"+scenario_name
+         cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals)
+         setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable)
+         cvar_excess_variable.construct()
+
+         cvar_eta_variable_name = "CVAR_ETA"
+         cvar_eta_variable = Var(name=cvar_eta_variable_name)
+         setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable)
+         cvar_eta_variable.construct()
+
+         compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS"
+         compute_excess_constraint = Constraint(name=compute_excess_constraint_name)
+         compute_excess_expression = cvar_excess_variable
+         for node in scenario_tree._scenario_map[scenario_name]._node_list:
+            (cost_variable, cost_variable_idx) = node._stage._cost_variable
+            compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx]
+         compute_excess_expression += cvar_eta_variable
+         compute_excess_constraint.add(None, (0.0, compute_excess_expression, None))
+         compute_excess_constraint._model = scenario_instance
+         setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint)
+
+         eta_equality_constraint_name = "MASTER_ETA_EQUALITY_WITH_" + scenario_instance.name
+         eta_equality_constraint = Constraint(name=eta_equality_constraint_name)
+         eta_equality_expr = master_cvar_eta_variable - cvar_eta_variable
+         eta_equality_constraint.add(None, (0.0, eta_equality_expr, 0.0))
+         eta_equality_constraint._model = binding_instance
+         setattr(binding_instance, eta_equality_constraint_name, eta_equality_constraint)
+
+      # add the constraint to compute the master CVaR variable value. iterate
+      # over scenario instances to create the expected excess component first.
+      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
+      cvar_cost_variable = getattr(binding_instance, cvar_cost_variable_name)
+      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
+      cvar_eta_variable = getattr(binding_instance, cvar_eta_variable_name)
+      
+      cvar_cost_expression = cvar_cost_variable - cvar_eta_variable
+      
+      for scenario_name in scenario_instances.keys():
+         scenario_instance = scenario_instances[scenario_name]
+         scenario_probability = scenario_tree._scenario_map[scenario_name]._probability
+
+         scenario_excess_variable_name = "CVAR_EXCESS_"+scenario_name
+         scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name)
+
+         cvar_cost_expression = cvar_cost_expression - (scenario_probability * scenario_excess_variable) / (1.0 - risk_alpha)
+
+      compute_cvar_cost_constraint_name = "COMPUTE_CVAR_COST"
+      compute_cvar_cost_constraint = Constraint(name=compute_cvar_cost_constraint_name)
+      compute_cvar_cost_constraint.add(None, (0.0, cvar_cost_expression, 0.0))
+      compute_cvar_cost_constraint._model = binding_instance
+      setattr(binding_instance, compute_cvar_cost_constraint_name, compute_cvar_cost_constraint)
+
+   # after mucking with instances, presolve to collect terms required prior to output.         
+   # IMPT: Do the scenario instances first, as the master depends on variables in the scenarios.
+   for scenario_name in scenario_instances.keys():
+      scenario_instance = scenario_instances[scenario_name]   
+      scenario_instance.preprocess()               
+
+   binding_instance.preprocess()
+
+   return binding_instance
+
+#
+# write the EF binding instance and all sub-instances. currently only outputs the CPLEX LP file format.
+#
+
+def write_ef(binding_instance, scenario_instances, output_filename):
+
+   # create the output file.
+   problem_writer = cpxlp.ProblemWriter_cpxlp()
+   output_file = open(output_filename,"w")
+
+   problem_writer._output_prefixes = True # we always want prefixes
+
+   ################################################################################################
+   #### WRITE THE MASTER OBJECTIVE ################################################################
+   ################################################################################################
+
+   # write the objective for the master binding instance.
+   problem_writer._output_objectives = True
+   problem_writer._output_constraints = False
+   problem_writer._output_variables = False
+
+   print >>output_file, "\\ Begin objective block for master"
+   problem_writer._print_model_LP(binding_instance, output_file)
+   print >>output_file, "\\ End objective block for master"
+   print >>output_file, ""
+
+   ################################################################################################
+   #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
+   ################################################################################################
+
+   print >>output_file, "s.t."
+   print >>output_file, ""
+   
+   problem_writer._output_objectives = False
+   problem_writer._output_constraints = True
+   problem_writer._output_variables = False
+
+   print >>output_file, "\\ Begin constraint block for master"
+   problem_writer._print_model_LP(binding_instance, output_file)
+   print >>output_file, "\\ End constraint block for master",
+   print >>output_file, ""
+
+   for scenario_name in scenario_instances.keys():
+      instance = scenario_instances[scenario_name]
+      print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
+      problem_writer._print_model_LP(instance, output_file)
+      print >>output_file, "\\ End constraint block for scenario",scenario_name
+      print >>output_file, ""
+
+   ################################################################################################
+   #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
+   ################################################################################################
+
+   # write the variables for the master binding instance, and then for each scenario.
+   print >>output_file, "bounds"
+   print >>output_file, ""
+   
+   problem_writer._output_objectives = False
+   problem_writer._output_constraints = False
+   problem_writer._output_variables = True
+
+   # first step: write variable bounds
+
+   problem_writer._output_continuous_variables = True
+   problem_writer._output_integer_variables = False
+   problem_writer._output_binary_variables = False
+
+   print >>output_file, "\\ Begin variable bounds block for master"
+   problem_writer._print_model_LP(binding_instance, output_file)
+   print >>output_file, "\\ End variable bounds block for master"
+   print >>output_file, ""
+   
+   for scenario_name in scenario_instances.keys():
+      instance = scenario_instances[scenario_name]
+      print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name 
+      problem_writer._print_model_LP(instance, output_file)
+      print >>output_file, "\\ End variable bounds block for scenario",scenario_name
+      print >>output_file, ""
+
+   # second step: write integer indicators.
+
+   problem_writer._output_continuous_variables = False
+   problem_writer._output_integer_variables = True
+
+   if integers_present(binding_instance, scenario_instances) is True:
+
+      print >>output_file, "integer"
+      print >>output_file, ""
+
+      print >>output_file, "\\ Begin integer variable block for master"
+      problem_writer._print_model_LP(binding_instance, output_file)
+      print >>output_file, "\\ End integer variable block for master"
+      print >>output_file, ""
+   
+      for scenario_name in scenario_instances.keys():
+         instance = scenario_instances[scenario_name]
+         print >>output_file, "\\ Begin integer variable block for scenario",scenario_name 
+         problem_writer._print_model_LP(instance, output_file)
+         print >>output_file, "\\ End integer variable block for scenario",scenario_name
+         print >>output_file, ""
+
+   # third step: write binary indicators.
+
+   problem_writer._output_integer_variables = False
+   problem_writer._output_binary_variables = True
+
+   if binaries_present(binding_instance, scenario_instances) is True:
+
+      print >>output_file, "binary"
+      print >>output_file, ""
+
+      print >>output_file, "\\ Begin binary variable block for master"
+      problem_writer._print_model_LP(binding_instance, output_file)
+      print >>output_file, "\\ End binary variable block for master"
+      print >>output_file, ""
+   
+      for scenario_name in scenario_instances.keys():
+         instance = scenario_instances[scenario_name]
+         print >>output_file, "\\ Begin binary variable block for scenario",scenario_name 
+         problem_writer._print_model_LP(instance, output_file)
+         print >>output_file, "\\ End integer binary block for scenario",scenario_name
+         print >>output_file, ""
+
+   # wrap up.
+   print >>output_file, "end"
+
+   # clean up.
+   output_file.close()
+   
+
+#
 # the main extensive-form writer routine - including read of scenarios/etc.
-def write_ef_from_scratch(model_directory, instance_directory, output_filename, \
+#
+
+def write_ef_from_scratch(model_directory, instance_directory, output_filename,
+                          verbose_output,
                           generate_weighted_cvar, cvar_weight, risk_alpha):
 
@@ -61 +505 @@
    scenario_data_directory_name = instance_directory
 
-   print "Initializing extensive form writer"
-   print ""
-
-   ################################################################################################
-   #### INITIALIZATION ############################################################################
-   ################################################################################################
-
-   #
-   # validate cvar options, if specified.
-   #
-   if generate_weighted_cvar is True:
-      if cvar_weight <= 0.0:
-         raise RuntimeError, "Weight of CVaR term must be >= 0.0 - value supplied="+str(cvar_weight)
-      if (risk_alpha <= 0.0) or (risk_alpha >= 1.0):
-         raise RuntimeError, "CVaR risk alpha must be between 0 and 1, exclusive - value supplied="+str(risk_alpha)
-
-      print "Writing CVaR weighted objective"
-      print "CVaR term weight="+str(cvar_weight)
-      print "CVaR alpha="+str(risk_alpha)
+   if verbose_output is True:
+      print "Initializing extensive form writer"
       print ""
-   
+
    #
    # create and populate the core model
@@ -115 +542 @@
    except IOError:
       
-      print "***ERROR: Failed to load scenario reference instance data from file="+reference_instance_filename
+      print "***ERROR: Failed to load scenario reference instance data from file="+reference_scenario_filename
       return            
 
@@ -135 +562 @@
    # print the input tree for validation/information purposes.
    #
-   scenario_tree.pprint()
+   if verbose_output is True:
+      scenario_tree.pprint()
 
    #
@@ -145 +573 @@
       sys.exit(1)
    else:
-      print "Scenario tree is valid!"
+      if verbose_output is True:
+         print "Scenario tree is valid!"
    print ""
 
@@ -152 +581 @@
    #
 
-   instances = {}
+   # the construction of instances takes little overhead in terms of
+   # memory potentially lost in the garbage-collection sense (mainly
+   # only that due to parsing and instance simplification/prep-processing).
+   # to speed things along, disable garbage collection if it enabled in
+   # the first place through the instance construction process.
+   # IDEA: If this becomes too much for truly large numbers of scenarios,
+   #       we could manually collect every time X instances have been created.
+
+   re_enable_gc = False
+   if gc.isenabled() is True:
+      re_enable_gc = True
+      gc.disable()
+
+   scenario_instances = {}
    
    if scenario_tree._scenario_based_data == 1:
-      print "Scenario-based instance initialization enabled"
+      if verbose_output is True:
+         print "Scenario-based instance initialization enabled"
    else:
-      print "Node-based instance initialization enabled"
+      if verbose_output is True:
+         print "Node-based instance initialization enabled"
          
    for scenario in scenario_tree._scenarios:
@@ -168 +612 @@
          if scenario_tree._scenario_based_data == 1:
             scenario_data_filename = scenario_data_directory_name + os.sep + scenario._name + ".dat"
-#            print "Data for scenario=" + scenario._name + " loads from file=" + scenario_data_filename
             scenario_instance = master_scenario_model.create(scenario_data_filename)
          else:
@@ -176 +619 @@
             while current_node is not None:
                node_data_filename = scenario_data_directory_name + os.sep + current_node._name + ".dat"
-#               print "Node data for scenario=" + scenario._name + " partially loading from file=" + node_data_filename
                scenario_data.add(node_data_filename)
                current_node = current_node._parent
@@ -191 +633 @@
       
       scenario_instance.preprocess()
-      instances[scenario._name] = scenario_instance
-
-   print ""
-
-   ################################################################################################
-   #### CREATE THE MASTER / BINDING INSTANCE ######################################################
-   ################################################################################################
-
-   master_binding_instance = Model()
-   master_binding_instance.name = "MASTER"
-
-   # walk the scenario tree - create variables representing the common values for all scenarios
-   # associated with that node. the constraints will be created later. also create expected-cost
-   # variables for each node, to be computed via constraints/objectives defined in a subsequent pass.
-   # master variables are created for all nodes but those in the last stage. expected cost variables
-   # are, for no particularly good reason other than easy coding, created for nodes in all stages.
-   print "Creating variables for master binding instance"
-
-   for stage in scenario_tree._stages:
-
-      for (stage_variable, index_template, stage_variable_indices) in stage._variables:
-
-         print "Creating master variable and blending constraints for decision variable=", stage_variable, ", indices=",stage_variable_indices
-
-         for tree_node in stage._tree_nodes:
-
-            if stage != scenario_tree._stages[-1]:      
-
-               master_variable_name = tree_node._name + "_" + stage_variable.name
-
-               # because there may be a single stage variable and multiple indices, check
-               # for the existence of the variable at this node - if you don't, you'll
-               # inadvertently over-write what was there previously!
-               master_variable = None
-               try:
-                  master_variable = getattr(master_binding_instance, master_variable_name)
-               except:
-                  new_master_variable_index = stage_variable._index
-                  new_master_variable = None
-                  if (len(new_master_variable_index) is 1) and (None in new_master_variable_index):
-                     new_master_variable = Var(name=stage_variable.name)
-                  else:
-                     new_master_variable = Var(new_master_variable_index, name=stage_variable.name)
-                  new_master_variable.construct()
-                  new_master_variable._model = master_binding_instance
-                  setattr(master_binding_instance, master_variable_name, new_master_variable)
-
-                  # TBD - TECHNICALLY, WE NEED TO COPY BOUNDS - BUT WE REALLY DON'T, AS THEY ARE ON THE PER-INSTNACE VARS!
-
-                  master_variable = new_master_variable
-
-               for index in stage_variable_indices:
-
-                  is_used = True # until proven otherwise                     
-                  for scenario in tree_node._scenarios:
-                     instance = instances[scenario._name]
-                     if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
-                        is_used = False
-
-                  is_fixed = False # until proven otherwise
-                  for scenario in tree_node._scenarios:
-                     instance = instances[scenario._name]
-                     if getattr(instance,stage_variable.name)[index].fixed is True:
-                        is_fixed = True
-
-                  if (is_used is True) and (is_fixed is False):
-                           
-                     # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
-                     # and because presolve/simplification is name-based, the names *have* to be different.
-                     master_variable[index].var = master_variable
-                     master_variable[index].name = tree_node._name + "_" + master_variable[index].name
-
-                     for scenario in tree_node._scenarios:
-
-                        scenario_instance = instances[scenario._name]
-                        scenario_variable = getattr(scenario_instance, stage_variable.name)
-                        new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
-                        new_constraint = Constraint(name=new_constraint_name)
-                        new_expr = master_variable[index] - scenario_variable[index]
-                        new_constraint.add(None, (0.0, new_expr, 0.0))
-                        new_constraint._model = master_binding_instance
-                        setattr(master_binding_instance, new_constraint_name, new_constraint)
-
-            # create a variable to represent the expected cost at this node -
-            # the constraint to compute this comes later.
-            expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
-            expected_cost_variable = Var(name=expected_cost_variable_name)
-            expected_cost_variable._model = master_binding_instance
-            setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable)
-
-   # if we're generating the weighted CVaR objective term, create the corresponding variable and
-   # the master CVaR eta variable.
-   if generate_weighted_cvar is True:
-      root_node = scenario_tree._stages[0]._tree_nodes[0]
-      
-      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
-      cvar_cost_variable = Var(name=cvar_cost_variable_name)
-      setattr(master_binding_instance, cvar_cost_variable_name, cvar_cost_variable)
-      cvar_cost_variable.construct()
-
-      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
-      cvar_eta_variable = Var(name=cvar_eta_variable_name)
-      setattr(master_binding_instance, cvar_eta_variable_name, cvar_eta_variable)      
-      cvar_eta_variable.construct()
-
-   master_binding_instance.preprocess()
-
-   # ditto above for the (non-expected) cost variable.
-   for stage in scenario_tree._stages:
-
-      (cost_variable,cost_variable_index) = stage._cost_variable
-
-      print "Creating master variable and blending constraints for cost variable=", cost_variable, ", index=", cost_variable_index      
-
-      for tree_node in stage._tree_nodes:
-
-         # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)         
-
-         # this is undoubtedly wasteful, in that a cost variable
-         # for each tree node is created with *all* indices.          
-         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
-         new_cost_variable_index = cost_variable._index
-         new_cost_variable = None
-         if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index):
-            new_cost_variable = Var(name=new_cost_variable_name)
-         else:
-            new_cost_variable = Var(new_cost_variable_index, name=new_cost_variable_name)
-         new_cost_variable.construct()
-         new_cost_variable._model = master_binding_instance
-         setattr(master_binding_instance, new_cost_variable_name, new_cost_variable)         
-
-         # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
-         new_cost_variable[cost_variable_index].var = new_cost_variable
-         if cost_variable_index is not None:
-            # if the variable index is None, the variable is derived from a VarValue, so the
-            # name gets updated automagically.
-            new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
-
-         for scenario in tree_node._scenarios:
-
-            scenario_instance = instances[scenario._name]
-            scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
-            new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
-            new_constraint = Constraint(name=new_constraint_name)
-            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
-            new_constraint.add(None, (0.0, new_expr, 0.0))
-            new_constraint._model = master_binding_instance
-            setattr(master_binding_instance, new_constraint_name, new_constraint)
-
-   # create the constraints for computing the master per-node cost variables,
-   # i.e., the current node cost and the expected cost of the child nodes.
-   # if the root, then the constraint is just the objective.
-
-   for stage in scenario_tree._stages:
-
-      (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
-
-      for tree_node in stage._tree_nodes:
-
-         node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
-         node_expected_cost_variable = getattr(master_binding_instance, node_expected_cost_variable_name)
-
-         node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
-         node_cost_variable = getattr(master_binding_instance, node_cost_variable_name)                        
-            
-         constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
-
-         for child_node in tree_node._children:
-
-            child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
-            child_node_expected_cost_variable = getattr(master_binding_instance, child_node_expected_cost_variable_name)
-            constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
-
-         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
-         new_constraint = Constraint(name=new_constraint_name)
-         new_constraint.add(None, (0.0, constraint_expr, 0.0))
-         new_constraint._model = master_binding_instance                     
-         setattr(master_binding_instance, new_constraint_name, new_constraint)
-
-         if tree_node._parent is None:
-
-            an_instance = instances[instances.keys()[0]]
-            an_objective = an_instance.active_components(Objective)
-            opt_sense = an_objective[an_objective.keys()[0]].sense
-
-            opt_expression = node_expected_cost_variable
-
-            if generate_weighted_cvar is True:
-               cvar_cost_variable_name = "CVAR_COST_" + tree_node._name
-               cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name)
-               opt_expression += cvar_weight * cvar_cost_variable
-
-            new_objective = Objective(name="MASTER", sense=opt_sense)
-            new_objective._data[None].expr = opt_expression
-            setattr(master_binding_instance, "MASTER", new_objective)
-
-   # CVaR requires the addition of a variable per scenario to represent the cost excess,
-   # and a constraint to compute the cost excess relative to eta. we also replicate (following
-   # what we do for node cost variables) an eta variable for each scenario instance, and
-   # require equality with the master eta variable via constraints.
-   if generate_weighted_cvar is True:
-      
-      root_node = scenario_tree._stages[0]._tree_nodes[0]
-
-      master_cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
-      master_cvar_eta_variable = getattr(master_binding_instance, master_cvar_eta_variable_name)
-      
-      for scenario_name in instances.keys():
-         scenario_instance = instances[scenario_name]
-
-         # unique names are required because the presolve isn't
-         # aware of the "owning" models for variables.
-         cvar_excess_variable_name = "CVAR_EXCESS_"+scenario_name
-         cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals)
-         setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable)
-         cvar_excess_variable.construct()
-
-         cvar_eta_variable_name = "CVAR_ETA"
-         cvar_eta_variable = Var(name=cvar_eta_variable_name)
-         setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable)
-         cvar_eta_variable.construct()
-
-         compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS"
-         compute_excess_constraint = Constraint(name=compute_excess_constraint_name)
-         compute_excess_expression = cvar_excess_variable
-         for node in scenario_tree._scenario_map[scenario_name]._node_list:
-            (cost_variable, cost_variable_idx) = node._stage._cost_variable
-            compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx]
-         compute_excess_expression += cvar_eta_variable
-         compute_excess_constraint.add(None, (0.0, compute_excess_expression, None))
-         compute_excess_constraint._model = scenario_instance
-         setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint)
-
-         eta_equality_constraint_name = "MASTER_ETA_EQUALITY_WITH_" + scenario_instance.name
-         eta_equality_constraint = Constraint(name=eta_equality_constraint_name)
-         eta_equality_expr = master_cvar_eta_variable - cvar_eta_variable
-         eta_equality_constraint.add(None, (0.0, eta_equality_expr, 0.0))
-         eta_equality_constraint._model = master_binding_instance
-         setattr(master_binding_instance, eta_equality_constraint_name, eta_equality_constraint)
-
-      # add the constraint to compute the master CVaR variable value. iterate
-      # over scenario instances to create the expected excess component first.
-      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
-      cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name)
-      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
-      cvar_eta_variable = getattr(master_binding_instance, cvar_eta_variable_name)
-      
-      cvar_cost_expression = cvar_cost_variable - cvar_eta_variable
-      
-      for scenario_name in instances.keys():
-         scenario_instance = instances[scenario_name]
-         scenario_probability = scenario_tree._scenario_map[scenario_name]._probability
-
-         scenario_excess_variable_name = "CVAR_EXCESS_"+scenario_name
-         scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name)
-
-         cvar_cost_expression = cvar_cost_expression - (scenario_probability * scenario_excess_variable) / (1.0 - risk_alpha)
-
-      compute_cvar_cost_constraint_name = "COMPUTE_CVAR_COST"
-      compute_cvar_cost_constraint = Constraint(name=compute_cvar_cost_constraint_name)
-      compute_cvar_cost_constraint.add(None, (0.0, cvar_cost_expression, 0.0))
-      compute_cvar_cost_constraint._model = master_binding_instance
-      setattr(master_binding_instance, compute_cvar_cost_constraint_name, compute_cvar_cost_constraint)
-
-   # after mucking with instances, presolve to collect terms required prior to output.         
-   # IMPT: Do the scenario instances first, as the master depends on variables in the scenarios.
-   for scenario_name in instances.keys():
-      scenario_instance = instances[scenario_name]   
-      scenario_instance.preprocess()         
-
-   master_binding_instance.preprocess()
-
-   ################################################################################################
-   #### WRITE THE COMPOSITE MODEL #################################################################
-   ################################################################################################
+      scenario_instances[scenario._name] = scenario_instance
+
+   if re_enable_gc is True:
+      gc.enable()
+
+   print ""
+
+   binding_instance = create_ef_instance(scenario_tree, scenario_instances,
+                                         verbose_output = verbose_output,
+                                         generate_weighted_cvar = generate_weighted_cvar,
+                                         cvar_weight = cvar_weight,
+                                         risk_alpha = risk_alpha)   
 
    print ""
    print "Starting to write extensive form"
 
-   # create the output file.
-   problem_writer = cpxlp.ProblemWriter_cpxlp()
-   output_file = open(output_filename,"w")
-
-   problem_writer._output_prefixes = True # we always want prefixes
-
-   ################################################################################################
-   #### WRITE THE MASTER OBJECTIVE ################################################################
-   ################################################################################################
-
-   # write the objective for the master binding instance.
-   problem_writer._output_objectives = True
-   problem_writer._output_constraints = False
-   problem_writer._output_variables = False
-
-   print >>output_file, "\\ Begin objective block for master"
-   problem_writer._print_model_LP(master_binding_instance, output_file)
-   print >>output_file, "\\ End objective block for master"
-   print >>output_file, ""
-
-   ################################################################################################
-   #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
-   ################################################################################################
-
-   print >>output_file, "s.t."
-   print >>output_file, ""
-   
-   problem_writer._output_objectives = False
-   problem_writer._output_constraints = True
-   problem_writer._output_variables = False
-
-   print >>output_file, "\\ Begin constraint block for master"
-   problem_writer._print_model_LP(master_binding_instance, output_file)
-   print >>output_file, "\\ End constraint block for master",
-   print >>output_file, ""
-
-   for scenario_name in instances.keys():
-      instance = instances[scenario_name]
-      print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
-      problem_writer._print_model_LP(instance, output_file)
-      print >>output_file, "\\ End constraint block for scenario",scenario_name
-      print >>output_file, ""
-
-   ################################################################################################
-   #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
-   ################################################################################################
-
-   # write the variables for the master binding instance, and then for each scenario.
-   print >>output_file, "bounds"
-   print >>output_file, ""
-   
-   problem_writer._output_objectives = False
-   problem_writer._output_constraints = False
-   problem_writer._output_variables = True
-
-   # first step: write variable bounds
-
-   problem_writer._output_continuous_variables = True
-   problem_writer._output_integer_variables = False
-   problem_writer._output_binary_variables = False
-
-   print >>output_file, "\\ Begin variable bounds block for master"
-   problem_writer._print_model_LP(master_binding_instance, output_file)
-   print >>output_file, "\\ End variable bounds block for master"
-   print >>output_file, ""
-
-   for scenario_name in instances.keys():
-      instance = instances[scenario_name]
-      print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name 
-      problem_writer._print_model_LP(instance, output_file)
-      print >>output_file, "\\ End variable bounds block for scenario",scenario_name
-      print >>output_file, ""
-
-   # second step: write integer indicators.
-
-   problem_writer._output_continuous_variables = False
-   problem_writer._output_integer_variables = True
-
-   if integers_present(master_binding_instance, instances) is True:
-
-      print >>output_file, "integer"
-      print >>output_file, ""
-
-      print >>output_file, "\\ Begin integer variable block for master"
-      problem_writer._print_model_LP(master_binding_instance, output_file)
-      print >>output_file, "\\ End integer variable block for master"
-      print >>output_file, ""
-   
-      for scenario_name in instances.keys():
-         instance = instances[scenario_name]
-         print >>output_file, "\\ Begin integer variable block for scenario",scenario_name 
-         problem_writer._print_model_LP(instance, output_file)
-         print >>output_file, "\\ End integer variable block for scenario",scenario_name
-         print >>output_file, ""
-
-   # third step: write binary indicators.
-
-   problem_writer._output_integer_variables = False
-   problem_writer._output_binary_variables = True
-
-   if binaries_present(master_binding_instance, instances) is True:
-
-      print >>output_file, "binary"
-      print >>output_file, ""
-
-      print >>output_file, "\\ Begin binary variable block for master"
-      problem_writer._print_model_LP(master_binding_instance, output_file)
-      print >>output_file, "\\ End binary variable block for master"
-      print >>output_file, ""
-   
-      for scenario_name in instances.keys():
-         instance = instances[scenario_name]
-         print >>output_file, "\\ Begin binary variable block for scenario",scenario_name 
-         problem_writer._print_model_LP(instance, output_file)
-         print >>output_file, "\\ End integer binary block for scenario",scenario_name
-         print >>output_file, ""
-
-   # wrap up.
-   print >>output_file, "end"
-
-   # clean up.
-   output_file.close()
+   write_ef(binding_instance, scenario_instances, output_filename)   
 
    print ""
@@ -605 +663 @@
    print ""
 
-def write_ef(scenario_tree, instances, output_filename):
+def create_and_write_ef(scenario_tree, scenario_instances, output_filename):
 
    start_time = time.time()
 
-   ################################################################################################
-   #### CREATE THE MASTER / BINDING INSTANCE ######################################################
-   ################################################################################################
-
-   master_binding_instance = Model()
-   master_binding_instance.name = "MASTER"
-
-   # walk the scenario tree - create variables representing the common values for all scenarios
-   # associated with that node. the constraints will be created later. also create expected-cost
-   # variables for each node, to be computed via constraints/objectives defined in a subsequent pass.
-   # master variables are created for all nodes but those in the last stage. expected cost variables
-   # are, for no particularly good reason other than easy coding, created for nodes in all stages.
-   print "Creating variables for master binding instance"
-
-   for stage in scenario_tree._stages:
-
-      for (stage_variable, index_template, stage_variable_indices) in stage._variables:
-
-         print "Creating master variable and blending constraints for decision variable=", stage_variable, ", indices=", index_template
-
-         for tree_node in stage._tree_nodes:
-
-            if stage != scenario_tree._stages[-1]:      
-
-               master_variable_name = tree_node._name + "_" + stage_variable.name
-
-               # because there may be a single stage variable and multiple indices, check
-               # for the existence of the variable at this node - if you don't, you'll
-               # inadvertently over-write what was there previously!
-               master_variable = None
-               try:
-                  master_variable = getattr(master_binding_instance, master_variable_name)
-               except:
-                  new_master_variable_index = stage_variable._index
-                  new_master_variable = None
-                  if (len(new_master_variable_index) is 1) and (None in new_master_variable_index):
-                     new_master_variable = Var(name=stage_variable.name)
-                  else:
-                     new_master_variable = Var(new_master_variable_index, name=stage_variable.name)
-                  new_master_variable.construct()
-                  new_master_variable._model = master_binding_instance
-                  setattr(master_binding_instance, master_variable_name, new_master_variable)
-
-                  # TBD - TECHNICALLY, WE NEED TO COPY BOUNDS - BUT WE REALLY DON'T, AS THEY ARE ON THE PER-INSTNACE VARS!
-
-                  master_variable = new_master_variable
-
-               for index in stage_variable_indices:
-
-                  is_used = True # until proven otherwise                     
-                  for scenario in tree_node._scenarios:
-                     instance = instances[scenario._name]
-                     if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
-                        is_used = False
-
-                  is_fixed = False # until proven otherwise
-                  for scenario in tree_node._scenarios:
-                     instance = instances[scenario._name]
-                     if getattr(instance,stage_variable.name)[index].fixed is True:
-                        is_fixed = True
-
-                  if (is_used is True) and (is_fixed is False):
-                           
-                     # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
-                     # and because presolve/simplification is name-based, the names *have* to be different.
-                     master_variable[index].var = master_variable
-                     master_variable[index].name = tree_node._name + "_" + master_variable[index].name
-
-                     for scenario in tree_node._scenarios:
-
-                        scenario_instance = instances[scenario._name]
-                        scenario_variable = getattr(scenario_instance, stage_variable.name)
-                        new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
-                        new_constraint = Constraint(name=new_constraint_name)
-                        new_expr = master_variable[index] - scenario_variable[index]
-                        new_constraint.add(None, (0.0, new_expr, 0.0))
-                        new_constraint._model = master_binding_instance
-                        setattr(master_binding_instance, new_constraint_name, new_constraint)
-
-            # create a variable to represent the expected cost at this node -
-            # the constraint to compute this comes later.
-            expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
-            expected_cost_variable = Var(name=expected_cost_variable_name)
-            expected_cost_variable._model = master_binding_instance
-            setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable)
-
-   master_binding_instance.preprocess()
-
-   # ditto above for the (non-expected) cost variable.
-   for stage in scenario_tree._stages:
-
-      (cost_variable,cost_variable_index) = stage._cost_variable
-
-      print "Creating master variable and blending constraints for cost variable=", cost_variable, ", index=", cost_variable_index      
-
-      for tree_node in stage._tree_nodes:
-
-         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
-
-         # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)
-
-         # this is undoubtedly wasteful, in that a cost variable
-         # for each tree node is created with *all* indices. 
-         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
-         new_cost_variable_index = cost_variable._index
-         new_cost_variable = None
-         if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index):
-            new_cost_variable = Var(name=new_cost_variable_name)
-         else:
-            new_cost_variable = Var(new_cost_variable_index, new_cost_variable_name)
-         new_cost_variable.construct()
-         new_cost_variable._model = master_binding_instance
-         setattr(master_binding_instance, new_cost_variable_name, new_cost_variable)                  
-
-         # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
-         new_cost_variable[cost_variable_index].var = new_cost_variable
-         if cost_variable_index is not None:
-            # if the variable index is None, the variable is derived from a VarValue, so the
-            # name gets updated automagically.
-            new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
-
-         for scenario in tree_node._scenarios:
-
-            scenario_instance = instances[scenario._name]
-            scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
-            new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
-            new_constraint = Constraint(name=new_constraint_name)
-            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
-            new_constraint.add(None, (0.0, new_expr, 0.0))
-            new_constraint._model = master_binding_instance
-            setattr(master_binding_instance, new_constraint_name, new_constraint)
-
-   # create the constraints for computing the master per-node cost variables,
-   # i.e., the current node cost and the expected cost of the child nodes.
-   # if the root, then the constraint is just the objective.
-
-   for stage in scenario_tree._stages:
-
-      (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
-
-      for tree_node in stage._tree_nodes:
-
-         node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
-         node_expected_cost_variable = getattr(master_binding_instance, node_expected_cost_variable_name)
-
-         node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
-         node_cost_variable = getattr(master_binding_instance, node_cost_variable_name)                        
-            
-         constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
-
-         for child_node in tree_node._children:
-
-            child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
-            child_node_expected_cost_variable = getattr(master_binding_instance, child_node_expected_cost_variable_name)
-            constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
-
-         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
-         new_constraint = Constraint(name=new_constraint_name)
-         new_constraint.add(None, (0.0, constraint_expr, 0.0))
-         new_constraint._model = master_binding_instance                     
-         setattr(master_binding_instance, new_constraint_name, new_constraint)
-
-         if tree_node._parent is None:
-
-            an_instance = instances[instances.keys()[0]]
-            an_objective = an_instance.active_components(Objective)
-            opt_sense = an_objective[an_objective.keys()[0]].sense
-
-            new_objective = Objective(name="MASTER", sense=opt_sense)
-            new_objective._data[None].expr = node_expected_cost_variable
-            setattr(master_binding_instance, "MASTER", new_objective)
-
-   master_binding_instance.preprocess()
-
-   ################################################################################################
-   #### WRITE THE COMPOSITE MODEL #################################################################
-   ################################################################################################
+   binding_instance = create_ef_instance(scenario_tree, scenario_instances)
 
    print ""
    print "Starting to write extensive form"
 
-   # create the output file.
-   problem_writer = cpxlp.ProblemWriter_cpxlp()
-   output_file = open(output_filename,"w")
-
-   problem_writer._output_prefixes = True # we always want prefixes
-
-   ################################################################################################
-   #### WRITE THE MASTER OBJECTIVE ################################################################
-   ################################################################################################
-
-   # write the objective for the master binding instance.
-   problem_writer._output_objectives = True
-   problem_writer._output_constraints = False
-   problem_writer._output_variables = False
-
-   print >>output_file, "\\ Begin objective block for master"
-   problem_writer._print_model_LP(master_binding_instance, output_file)
-   print >>output_file, "\\ End objective block for master"
-   print >>output_file, ""
-
-   ################################################################################################
-   #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
-   ################################################################################################
-
-   print >>output_file, "s.t."
-   print >>output_file, ""
-   
-   problem_writer._output_objectives = False
-   problem_writer._output_constraints = True
-   problem_writer._output_variables = False
-
-   print >>output_file, "\\ Begin constraint block for master"
-   problem_writer._print_model_LP(master_binding_instance, output_file)
-   print >>output_file, "\\ End constraint block for master",
-   print >>output_file, ""
-
-   for scenario_name in instances.keys():
-      instance = instances[scenario_name]
-      print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
-      problem_writer._print_model_LP(instance, output_file)
-      print >>output_file, "\\ End constraint block for scenario",scenario_name
-      print >>output_file, ""
-
-   ################################################################################################
-   #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
-   ################################################################################################
-
-   # write the variables for the master binding instance, and then for each scenario.
-   print >>output_file, "bounds"
-   print >>output_file, ""
-   
-   problem_writer._output_objectives = False
-   problem_writer._output_constraints = False
-   problem_writer._output_variables = True
-
-   # first step: write variable bounds
-
-   problem_writer._output_continuous_variables = True
-   problem_writer._output_integer_variables = False
-   problem_writer._output_binary_variables = False
-
-   print >>output_file, "\\ Begin variable bounds block for master"
-   problem_writer._print_model_LP(master_binding_instance, output_file)
-   print >>output_file, "\\ End variable bounds block for master"
-   print >>output_file, ""
-   
-   for scenario_name in instances.keys():
-      instance = instances[scenario_name]
-      print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name 
-      problem_writer._print_model_LP(instance, output_file)
-      print >>output_file, "\\ End variable bounds block for scenario",scenario_name
-      print >>output_file, ""
-
-   # second step: write integer indicators.
-
-   problem_writer._output_continuous_variables = False
-   problem_writer._output_integer_variables = True
-
-   if integers_present(master_binding_instance, instances) is True:
-
-      print >>output_file, "integer"
-      print >>output_file, ""
-
-      print >>output_file, "\\ Begin integer variable block for master"
-      problem_writer._print_model_LP(master_binding_instance, output_file)
-      print >>output_file, "\\ End integer variable block for master"
-      print >>output_file, ""
-   
-      for scenario_name in instances.keys():
-         instance = instances[scenario_name]
-         print >>output_file, "\\ Begin integer variable block for scenario",scenario_name 
-         problem_writer._print_model_LP(instance, output_file)
-         print >>output_file, "\\ End integer variable block for scenario",scenario_name
-         print >>output_file, ""
-
-   # third step: write binary indicators.
-
-   problem_writer._output_integer_variables = False
-   problem_writer._output_binary_variables = True
-
-   if binaries_present(master_binding_instance, instances) is True:
-
-      print >>output_file, "binary"
-      print >>output_file, ""
-
-      print >>output_file, "\\ Begin binary variable block for master"
-      problem_writer._print_model_LP(master_binding_instance, output_file)
-      print >>output_file, "\\ End binary variable block for master"
-      print >>output_file, ""
-   
-      for scenario_name in instances.keys():
-         instance = instances[scenario_name]
-         print >>output_file, "\\ Begin binary variable block for scenario",scenario_name 
-         problem_writer._print_model_LP(instance, output_file)
-         print >>output_file, "\\ End integer binary block for scenario",scenario_name
-         print >>output_file, ""
-
-   # wrap up.
-   print >>output_file, "end"
-
-   # clean up.
-   output_file.close()
+   write_ef(binding_instance, scenario_instances, output_filename)
 
    print ""

coopr.pysp/stable/2.3/coopr/pysp/ef_writer_script.py

@@ -24 +24 @@
 
 #
-# Setup command-line options
+# utility method to construct an option parser for ef writer arguments
 #
 
-parser = OptionParser()
-parser.add_option("--model-directory",
-                  help="The directory in which all model (reference and scenario) definitions are stored. Default is \".\".",
-                  action="store",
-                  dest="model_directory",
-                  type="string",
-                  default=".")
-parser.add_option("--instance-directory",
-                  help="The directory in which all instance (reference and scenario) definitions are stored. Default is \".\".",
-                  action="store",
-                  dest="instance_directory",
-                  type="string",
-                  default=".")
-parser.add_option("--generate-weighted-cvar",
-                  help="Add a weighted CVaR term to the primary objective",
-                  action="store_true",
-                  dest="generate_weighted_cvar",
-                  default=False)
-parser.add_option("--cvar-weight",
-                  help="The weight associated with the CVaR term in the risk-weighted objective formulation.",
-                  action="store",
-                  dest="cvar_weight",
-                  type="float",
-                  default=1.0)
-parser.add_option("--risk-alpha",
-                  help="The probability threshold associated with cvar (or any future) risk-oriented performance metrics.",
-                  action="store",
-                  dest="risk_alpha",
-                  type="float",
-                  default=0.95)
-parser.add_option("--output-file",
-                  help="Specify the name of the extensive form output file",
-                  action="store",
-                  dest="output_file",
-                  type="string",
-                  default="efout.lp")
-parser.add_option("--profile",
-                  help="Enable profiling of Python code.  The value of this option is the number of functions that are summarized.",
-                  action="store",
-                  dest="profile",
-                  default=0)
-parser.usage="efwriter [options]"
+def construct_ef_writer_options_parser(usage_string):
 
+   parser = OptionParser()
+   parser.add_option("--verbose",
+                     help="Generate verbose output, beyond the usual status output. Default is False.",
+                     action="store_true",
+                     dest="verbose",
+                     default=False)
+   parser.add_option("--model-directory",
+                     help="The directory in which all model (reference and scenario) definitions are stored. Default is \".\".",
+                     action="store",
+                     dest="model_directory",
+                     type="string",
+                     default=".")
+   parser.add_option("--instance-directory",
+                     help="The directory in which all instance (reference and scenario) definitions are stored. Default is \".\".",
+                     action="store",
+                     dest="instance_directory",
+                     type="string",
+                     default=".")
+   parser.add_option("--generate-weighted-cvar",
+                     help="Add a weighted CVaR term to the primary objective",
+                     action="store_true",
+                     dest="generate_weighted_cvar",
+                     default=False)
+   parser.add_option("--cvar-weight",
+                     help="The weight associated with the CVaR term in the risk-weighted objective formulation.",
+                     action="store",
+                     dest="cvar_weight",
+                     type="float",
+                     default=1.0)
+   parser.add_option("--risk-alpha",
+                     help="The probability threshold associated with cvar (or any future) risk-oriented performance metrics.",
+                     action="store",
+                     dest="risk_alpha",
+                     type="float",
+                     default=0.95)
+   parser.add_option("--output-file",
+                     help="Specify the name of the extensive form output file",
+                     action="store",
+                     dest="output_file",
+                     type="string",
+                     default="efout.lp")
+   parser.add_option("--profile",
+                     help="Enable profiling of Python code.  The value of this option is the number of functions that are summarized.",
+                     action="store",
+                     dest="profile",
+                     default=0)
+   parser.add_option("--disable-gc",
+                     help="Disable the python garbage collecter. Default is False.",
+                     action="store_true",
+                     dest="disable_gc",
+                     default=False)
+   parser.usage=usage_string
+
+   return parser
+   
 def run_ef_writer(options, args):
 
@@ -84 +98 @@
       risk_alpha = options.risk_alpha
 
-   write_ef_from_scratch(os.path.expanduser(options.model_directory), os.path.expanduser(options.instance_directory), os.path.expanduser(options.output_file), \
+   write_ef_from_scratch(os.path.expanduser(options.model_directory),
+                         os.path.expanduser(options.instance_directory),
+                         os.path.expanduser(options.output_file),
+                         options.verbose,
                          generate_weighted_cvar, cvar_weight, risk_alpha)
-
-   return
 
 def run(args=None):
@@ -96 +111 @@
     #
 
-    # for a one-pass execution, garbage collection doesn't make
-    # much sense - so I'm disabling it. Because: It drops the run-time
-    # by a factor of 3-4 on bigger instances.
-    gc.disable()       
-
     #
     # Parse command-line options.
     #
     try:
-       (options, args) = parser.parse_args(args=args)
+       options_parser = construct_ef_writer_options_parser("efwriter [options]")
+       (options, args) = options_parser.parse_args(args=args)
     except SystemExit:
        # the parser throws a system exit if "-h" is specified - catch
        # it to exit gracefully.
        return
+
+    if options.disable_gc is True:
+       gc.disable()
+    else:
+       gc.enable()
 
     if options.profile > 0:
@@ -140 +156 @@
 
     gc.enable()
+    
     return ans
 

coopr.pysp/stable/2.3/coopr/pysp/ph.py

@@ -27 +27 @@
 from scenariotree import *
 from phutils import *
+from phobjective import *
 
 from pyutilib.component.core import ExtensionPoint
@@ -33 +34 @@
 
 class ProgressiveHedging(object):
-
-   #
-   # routine to compute linearization breakpoints uniformly between the bounds and the mean.
-   #
-
-   # IMPT: In general, the breakpoint computation codes can return a 2-list even if the lb equals
-   #       the ub. This case happens quite often in real models (although typically lb=xvag=ub).
-   #       See the code for constructing the pieces on how this case is handled in the linearization.
-   
-   def compute_uniform_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints_per_side):
-      
-      breakpoints = []
-
-      # add the lower bound - the first breakpoint.
-      breakpoints.append(lb)
-
-      # determine the breakpoints to the left of the mean.
-      left_step = (xavg - lb) / num_breakpoints_per_side
-      current_x = lb
-      for i in range(1,num_breakpoints_per_side+1):
-         this_lb = current_x
-         this_ub = current_x+left_step
-         if (fabs(this_lb-lb) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
-            breakpoints.append(this_lb)
-         current_x += left_step            
-
-      # add the mean - it's always a breakpoint. unless!
-      # the lb or ub and the avg are the same.
-      if (fabs(lb-xavg) > self._integer_tolerance) and (fabs(ub-xavg) > self._integer_tolerance):
-         breakpoints.append(xavg)
-
-      # determine the breakpoints to the right of the mean.
-      right_step = (ub - xavg) / num_breakpoints_per_side
-      current_x = xavg
-      for i in range(1,num_breakpoints_per_side+1):
-         this_lb = current_x
-         this_ub = current_x+right_step
-         if (fabs(this_ub-xavg) > self._integer_tolerance) and (fabs(this_ub-ub) > self._integer_tolerance):         
-            breakpoints.append(this_ub)
-         current_x += right_step
-
-      # add the upper bound - the last breakpoint.
-      # the upper bound should always be different than the lower bound (I say with some
-      # hesitation - it really depends on what plugins are doing to modify the bounds dynamically).
-      breakpoints.append(ub)
-
-      return breakpoints
-
-   #
-   # routine to compute linearization breakpoints uniformly between the current node min/max bounds.
-   #   
-
-   def compute_uniform_between_nodestat_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints):
-
-      breakpoints = []
-
-      # add the lower bound - the first breakpoint.
-      breakpoints.append(lb)
-
-      # add the node-min - the second breakpoint. but only if it is different than the lower bound and the mean.
-      if (fabs(node_min-lb) > self._integer_tolerance) and (fabs(node_min-xavg) > self._integer_tolerance):      
-         breakpoints.append(node_min)
-
-      step = (node_max - node_min) / num_breakpoints
-      current_x = node_min
-      for i in range(1,num_breakpoints+1):
-         this_lb = current_x
-         this_ub = current_x+step
-         if (fabs(this_lb-node_min) > self._integer_tolerance) and (fabs(this_lb-node_max) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
-            breakpoints.append(this_lb)
-         current_x += step            
-
-      # add the node-max - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
-      if (fabs(node_max-ub) > self._integer_tolerance) and (fabs(node_max-xavg) > self._integer_tolerance):            
-         breakpoints.append(node_max)      
-
-      # add the upper bound - the last breakpoint.
-      breakpoints.append(ub)
-
-      # add the mean - it's always a breakpoint. unless! -
-      # it happens to be equal to (within tolerance) the lower or upper bounds.
-      # sort to insert it in the correct place.
-      if (fabs(xavg - lb) > self._integer_tolerance) and (fabs(xavg - ub) > self._integer_tolerance):
-         breakpoints.append(xavg)
-      breakpoints.sort()
-
-      return breakpoints
-
-   #
-   # routine to compute linearization breakpoints using "Woodruff" relaxation of the compute_uniform_between_nodestat_breakpoints.
-   #   
-
-   def compute_uniform_between_woodruff_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints):
-
-      breakpoints = []
-
-      # add the lower bound - the first breakpoint.
-      breakpoints.append(lb)
-
-      # be either three "differences" from the mean, or else "halfway to the bound", whichever is closer to the mean.
-      left = max(xavg - 3.0 * (xavg - node_min), xavg - 0.5 * (xavg - lb))
-      right = min(xavg + 3.0 * (node_max - xavg), xavg + 0.5 * (ub - xavg))      
-
-      # add the left bound - the second breakpoint. but only if it is different than the lower bound and the mean.
-      if (fabs(left-lb) > self._integer_tolerance) and (fabs(left-xavg) > self._integer_tolerance):      
-         breakpoints.append(left)
-
-      step = (right - left) / num_breakpoints
-      current_x = left
-      for i in range(1,num_breakpoints+1):
-         this_lb = current_x
-         this_ub = current_x+step
-         if (fabs(this_lb-left) > self._integer_tolerance) and (fabs(this_lb-right) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
-            breakpoints.append(this_lb)
-         current_x += step            
-
-      # add the right bound - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
-      if (fabs(right-ub) > self._integer_tolerance) and (fabs(right-xavg) > self._integer_tolerance):            
-         breakpoints.append(right)      
-
-      # add the upper bound - the last breakpoint.
-      breakpoints.append(ub)
-
-      # add the mean - it's always a breakpoint.
-      # sort to insert it in the correct place.
-      breakpoints.append(xavg)
-      breakpoints.sort()
-
-      return breakpoints
-
-   #
-   # routine to compute linearization breakpoints based on an exponential distribution from the mean in each direction.
-   #   
-
-   def compute_exponential_from_mean_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints_per_side):
-
-      breakpoints = []
-
-      # add the lower bound - the first breakpoint.
-      breakpoints.append(lb)
-
-      # determine the breakpoints to the left of the mean.
-      left_delta = xavg - lb
-      base = exp(log(left_delta) / num_breakpoints_per_side)
-      current_offset = base
-      for i in range(1,num_breakpoints_per_side+1):
-         current_x = xavg - current_offset
-         if (fabs(current_x-lb) > self._integer_tolerance) and (fabs(current_x-xavg) > self._integer_tolerance):
-            breakpoints.append(current_x)
-         current_offset *= base
-
-      # add the mean - it's always a breakpoint.
-      breakpoints.append(xavg)
-
-      # determine the breakpoints to the right of the mean.
-      right_delta = ub - xavg
-      base = exp(log(right_delta) / num_breakpoints_per_side)
-      current_offset = base
-      for i in range(1,num_breakpoints_per_side+1):
-         current_x = xavg + current_offset
-         if (fabs(current_x-xavg) > self._integer_tolerance) and (fabs(current_x-ub) > self._integer_tolerance):         
-            breakpoints.append(current_x)
-         current_offset *= base
-
-      # add the upper bound - the last breakpoint.
-      breakpoints.append(ub)         
-
-      return breakpoints      
 
    #
@@ -379 +212 @@
       return (num_fixed_discrete_vars, num_fixed_continuous_vars)
 
-   # a utility to create piece-wise linear constraint expressions for a given variable, for
-   # use in constructing the augmented (penalized) PH objective. lb and ub are the bounds
-   # on this piece, variable is the actual instance variable, and average is the instance
-   # parameter specifying the average of this variable across instances sharing passing
-   # through a common tree node. lb and ub are floats.
-   # IMPT: There are cases where lb=ub, in which case the slope is 0 and the intercept
-   #       is simply the penalty at the lower(or upper) bound.
-   def _create_piecewise_constraint_expression(self, lb, ub, instance_variable, variable_average, quad_variable):
-
-      penalty_at_lb = (lb - variable_average()) * (lb - variable_average())
-      penalty_at_ub = (ub - variable_average()) * (ub - variable_average())
-      slope = None
-      if fabs(ub-lb) > self._integer_tolerance:
-         slope = (penalty_at_ub - penalty_at_lb) / (ub - lb)
-      else:
-         slope = 0.0
-      intercept = penalty_at_lb - slope * lb
-      expression = (0.0, quad_variable - slope * instance_variable - intercept, None)
-
-      return expression
-
    # when the quadratic penalty terms are approximated via piecewise linear segments,
    # we end up (necessarily) "littering" the scenario instances with extra constraints.
@@ -424 +236 @@
       for (instance_name, instance) in self._instances.items():
 
-         # first, gather all unique variables referenced in any stage
-         # other than the last, independent of specific indices. this
-         # "gather" step is currently required because we're being lazy
-         # in terms of index management in the scenario tree - which
-         # should really be done in terms of sets of indices.
-         # NOTE: technically, the "instance variables" aren't really references
-         # to the variable in the instance - instead, the reference model. this
-         # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
-         instance_variables = {}
-         
-         for stage in self._scenario_tree._stages[:-1]:
-            
-            for (reference_variable, index_template, reference_indices) in stage._variables:
-                  
-               if reference_variable.name not in instance_variables.keys():
-                     
-                  instance_variables[reference_variable.name] = reference_variable
-
-         # for each blended variable, create a corresponding ph weight and average parameter in the instance.
-         # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor
-         # in the grand scheme of things.
-
-         for (variable_name, reference_variable) in instance_variables.items():
-
-            # PH WEIGHT
-
-            new_w_index = reference_variable._index
-            new_w_parameter_name = "PHWEIGHT_"+reference_variable.name
-            # this bit of ugliness is due to Pyomo not correctly handling the Param construction
-            # case when the supplied index set consists strictly of None, i.e., the source variable
-            # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
-            new_w_parameter = None
-            if (len(new_w_index) is 1) and (None in new_w_index):
-               new_w_parameter = Param(name=new_w_parameter_name)
-            else:
-               new_w_parameter = Param(new_w_index,name=new_w_parameter_name)
-            setattr(instance,new_w_parameter_name,new_w_parameter)
-
-            # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference
-            # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not
-            # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the
-            # issue in any case for now. NOTE: I'm not sure this is still the case in Pyomo.
-            for index in new_w_index:
-               new_w_parameter[index] = 0.0
-
-            # PH AVG
-
-            new_avg_index = reference_variable._index 
-            new_avg_parameter_name = "PHAVG_"+reference_variable.name
-            new_avg_parameter = None
-            if (len(new_avg_index) is 1) and (None in new_avg_index):
-               new_avg_parameter = Param(name=new_avg_parameter_name)
-            else:
-               new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
-            setattr(instance,new_avg_parameter_name,new_avg_parameter)
-
-            for index in new_avg_index:
-               new_avg_parameter[index] = 0.0
-
-            # PH RHO
-
-            new_rho_index = reference_variable._index 
-            new_rho_parameter_name = "PHRHO_"+reference_variable.name
-            new_rho_parameter = None
-            if (len(new_avg_index) is 1) and (None in new_avg_index):
-               new_rho_parameter = Param(name=new_rho_parameter_name)
-            else:
-               new_rho_parameter = Param(new_rho_index,name=new_rho_parameter_name)
-            setattr(instance,new_rho_parameter_name,new_rho_parameter)
-
-            for index in new_rho_index:
-               new_rho_parameter[index] = self._rho
-
-            # PH LINEARIZED PENALTY TERM
-
-            if self._linearize_nonbinary_penalty_terms > 0:
-
-               new_penalty_term_variable_index = reference_variable._index
-               new_penalty_term_variable_name = "PHQUADPENALTY_"+reference_variable.name
-               # this is a quadratic penalty term - the lower bound is 0!
-               new_penalty_term_variable = None
-               if (len(new_avg_index) is 1) and (None in new_avg_index):
-                  new_penalty_term_variable = Var(name=new_penalty_term_variable_name, bounds=(0.0,None))                  
-               else:
-                  new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None))
-               new_penalty_term_variable.construct()
-               setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable)
-               self._instance_augmented_attributes[instance_name].append(new_penalty_term_variable_name)
-
-            # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY. 
-
-            # also controls whether weight updates proceed at any iteration.
-
-            new_blend_index = reference_variable._index
-            new_blend_parameter_name = "PHBLEND_"+reference_variable.name
-            new_blend_parameter = None
-            if (len(new_avg_index) is 1) and (None in new_avg_index):
-               new_blend_parameter = Param(name=new_blend_parameter_name, within=Binary)
-            else:
-               new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary)
-            setattr(instance,new_blend_parameter_name,new_blend_parameter)
-
-            for index in new_blend_index:
-               new_blend_parameter[index] = 1
+         new_penalty_variable_names = create_ph_parameters(instance, self._scenario_tree, self._rho, self._linearize_nonbinary_penalty_terms)
+         self._instance_augmented_attributes[instance_name].extend(new_penalty_variable_names)
 
    # a simple utility to extract the first-stage cost statistics, e.g., min, average, and max.
@@ -554 +264 @@
                minimum_value = this_value
       return minimum_value, sum_values/num_values, maximum_value
-      
+
+   # a utility to transmit - across the PH solver manager - the current weights
+   # and averages for each of my problem instances. used to set up iteration K solves.
+   def _transmit_weights_and_averages(self):
+
+      for scenario_name, scenario_instance in self._instances.items():
+
+         weights_to_transmit = []
+         averages_to_transmit = []
+
+         for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.
+
+            for (variable, index_template, variable_indices) in stage._variables:
+
+               variable_name = variable.name
+               weight_parameter_name = "PHWEIGHT_"+variable_name
+               weights_to_transmit.append(getattr(scenario_instance, weight_parameter_name))
+               average_parameter_name = "PHAVG_"+variable_name
+               averages_to_transmit.append(getattr(scenario_instance, average_parameter_name))               
+
+         self._solver_manager.transmit_weights_and_averages(scenario_instance, weights_to_transmit, averages_to_transmit)
+
+   #
+   # a utility to transmit - across the PH solver manager - the current rho values
+   # for each of my problem instances. mainly after PH iteration 0 is complete,
+   # in preparation for the iteration K solves.
+   #
+
+   def _transmit_rhos(self):
+
+      for scenario_name, scenario_instance in self._instances.items():
+
+         rhos_to_transmit = []
+
+         for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no rhos to worry about.
+
+            for (variable, index_template, variable_indices) in stage._variables:
+
+               variable_name = variable.name
+               rho_parameter_name = "PHRHO_"+variable_name
+               rhos_to_transmit.append(getattr(scenario_instance, rho_parameter_name))
+
+         self._solver_manager.transmit_rhos(scenario_instance, rhos_to_transmit)
+
+   #
+   # a utility to transmit - across the PH solver manager - the current scenario
+   # tree node statistics to each of my problem instances. done prior to each
+   # PH iteration k. 
+   #
+
+   def _transmit_tree_node_statistics(self):
+
+      for scenario_name, scenario_instance in self._instances.items():
+
+         tree_node_minimums = {}
+         tree_node_maximums = {}
+
+         scenario = self._scenario_tree._scenario_map[scenario_name]
+
+         for tree_node in scenario._node_list:
+
+            tree_node_minimums[tree_node._name] = tree_node._minimums
+            tree_node_maximums[tree_node._name] = tree_node._maximums
+
+         self._solver_manager.transmit_tree_node_statistics(scenario_instance, tree_node_minimums, tree_node_maximums)         
+
+   #
+   # a utility to enable - across the PH solver manager - weighted penalty objectives.
+   #
+
+   def _enable_ph_objectives(self):
+
+      for scenario_name, scenario_instance in self._instances.items():
+
+         self._solver_manager.enable_ph_objective(scenario_instance)
+         
 
    """ Constructor
@@ -596 +381 @@
       self._solver_manager_type = "serial" # serial or pyro are the options currently available
       
-      self._solver = None # will eventually be unnecessary once Bill eliminates the need for a solver in the solver manager constructor.
+      self._solver = None 
       self._solver_manager = None 
       
@@ -1187 +972 @@
    def form_iteration_k_objectives(self):
 
-      # for each blended variable (i.e., those not appearing in the final stage),
-      # add the linear and quadratic penalty terms to the objective.
       for instance_name, instance in self._instances.items():
-         
-         objective_name = instance.active_components(Objective).keys()[0]         
-         objective = instance.active_components(Objective)[objective_name]
-         # clone the objective, because we're going to augment (repeatedly) the original objective.
-         objective_expression = self._original_objective_expression[instance_name].clone() 
-         # the quadratic expression is really treated as just a list - eventually should be treated as a full expression.
-         quad_expression = 0.0
-         
-         for stage in self._scenario_tree._stages[:-1]: # skip the last stage, as no blending occurs
-
-            variable_tree_node = None
-            for node in stage._tree_nodes:
-               for scenario in node._scenarios:
-                  if scenario._name == instance_name:
-                     variable_tree_node = node
-                     break
-
-            for (reference_variable, index_template, variable_indices) in stage._variables:
-
-               variable_name = reference_variable.name
-               variable_type = reference_variable.domain
-
-               w_parameter_name = "PHWEIGHT_"+variable_name
-               w_parameter = instance.active_components(Param)[w_parameter_name]
-                  
-               average_parameter_name = "PHAVG_"+variable_name
-               average_parameter = instance.active_components(Param)[average_parameter_name]
-
-               rho_parameter_name = "PHRHO_"+variable_name
-               rho_parameter = instance.active_components(Param)[rho_parameter_name]
-
-               blend_parameter_name = "PHBLEND_"+variable_name
-               blend_parameter = instance.active_components(Param)[blend_parameter_name]
-
-               node_min_parameter = variable_tree_node._minimums[variable_name]
-               node_max_parameter = variable_tree_node._maximums[variable_name]               
-
-               quad_penalty_term_variable = None
-               if self._linearize_nonbinary_penalty_terms > 0:
-                  quad_penalty_term_variable_name = "PHQUADPENALTY_"+variable_name
-                  quad_penalty_term_variable = instance.active_components(Var)[quad_penalty_term_variable_name]
-
-               instance_variable = instance.active_components(Var)[variable_name]
-
-               for index in variable_indices:
-
-                  if (instance_variable[index].status is not VarStatus.unused) and (instance_variable[index].fixed is False):
-
-                     # add the linear (w-weighted) term is a consistent fashion, independent of variable type. 
-                     # if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes for reasons currently unknown.
-                     objective_expression += (w_parameter[index] * instance_variable[index])
-
-                     # there are some cases in which a user may want to avoid the proximal term completely.
-                     # it's of course only a good idea when there are at least bounds (both lb and ub) on
-                     # the variables to be blended.
-                     if self._drop_proximal_terms is False:
-
-                        # deal with binaries
-                        if isinstance(variable_type, BooleanSet) is True:
-
-                           if self._retain_quadratic_binary_terms is False:
-                              # this rather ugly form of the linearized quadratic expression term is required
-                              # due to a pyomo bug - the form (rho/2) * (x+y+z) chokes in presolve when distributing
-                              # over the sum.
-                              new_term = (blend_parameter[index] * rho_parameter[index] / 2.0 * instance_variable[index]) - \
-                                         (blend_parameter[index] * rho_parameter[index] * average_parameter[index] * instance_variable[index]) + \
-                                         (blend_parameter[index] * rho_parameter[index] / 2.0 * average_parameter[index] * average_parameter[index])                              
-                              if objective.sense is minimize:
-                                 objective_expression += new_term
-                              else:
-                                 objective_expression -= new_term                                 
-                           else:
-                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)
-
-                        # deal with everything else
-                        else:
-
-                           if self._linearize_nonbinary_penalty_terms > 0:
-
-                              # the variables are easy - just a single entry.
-                              if objective.sense is minimize:
-                                 objective_expression += (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
-                              else:
-                                 objective_expression -= (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
-                                 
-                              # the constraints are somewhat nastier.
-
-                              # TBD - DEFINE CONSTRAINTS ON-THE-FLY?? (INDIVIDUALLY NAMED FOR NOW - CREATE INDEX SETS!) - OR A LEAST AN INDEX SET PER "PIECE"
-                              xavg = average_parameter[index]
-                              x = instance_variable[index]
-
-                              lb = None
-                              ub = None
-
-                              if x.lb is None:
-                                 raise ValueError, "No lower bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
-                              else:
-                                 lb = x.lb()
-
-                              if x.ub is None:
-                                 raise ValueError, "No upper bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
-                              else:
-                                 ub = x.ub()
-
-                              node_min = node_min_parameter[index]()
-                              node_max = node_max_parameter[index]()
-
-                              # compute the breakpoint sequence according to the specified strategy.
-                              breakpoints = []
-                              if self._breakpoint_strategy == 0:
-                                 breakpoints = self.compute_uniform_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
-                              elif self._breakpoint_strategy == 1:
-                                 breakpoints = self.compute_uniform_between_nodestat_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
-                              elif self._breakpoint_strategy == 2:
-                                 breakpoints = self.compute_uniform_between_woodruff_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
-                              elif self._breakpoint_strategy == 3:
-                                 breakpoints = self.compute_exponential_from_mean_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)                                                                                                   
-                              else:
-                                 raise ValueError, "A breakpoint distribution strategy="+str(self._breakpoint_strategy)+" is currently not supported within PH!"
-
-                              for i in range(0,len(breakpoints)-1):
-
-                                 this_lb = breakpoints[i]
-                                 this_ub = breakpoints[i+1]
-
-                                 piece_constraint_name = "QUAD_PENALTY_PIECE_"+str(i)+"_"+variable_name+str(index)
-                                 if hasattr(instance, piece_constraint_name) is False:
-                                    # this is the first time the constraint is being added - add it to the list of PH-specific constraints for this instance.
-                                    self._instance_augmented_attributes[instance_name].append(piece_constraint_name)                                                                           
-                                 piece_constraint = Constraint(name=piece_constraint_name)
-                                 piece_constraint.model = instance
-                                 piece_expression = self._create_piecewise_constraint_expression(this_lb, this_ub, x, xavg, quad_penalty_term_variable[index])
-                                 piece_constraint.add(None, piece_expression)
-                                 setattr(instance, piece_constraint_name, piece_constraint)                                    
-
-                           else:
-
-                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)                           
-                   
-         # strictly speaking, this probably isn't necessary - parameter coefficients won't get
-         # pre-processed out of the expression tree. however, if the under-the-hood should change,
-         # we'll be covered.
-         objective_expression.simplify(instance)
-         instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression
-         # if we are linearizing everything, then nothing will appear in the quadratic expression -
-         # don't add the empty "0.0" expression to the objective. otherwise, the output file won't
-         # be properly generated.
-         if quad_expression != 0.0:
-           instance.active_components(Objective)[objective_name]._quad_subexpr = quad_expression
+
+         new_attrs = form_ph_objective(instance_name, \
+                                       instance, \
+                                       self._original_objective_expression[instance_name], \
+                                       self._scenario_tree, \
+                                       self._linearize_nonbinary_penalty_terms, \
+                                       self._drop_proximal_terms, \
+                                       self._retain_quadratic_binary_terms, \
+                                       self._breakpoint_strategy, \
+                                       self._integer_tolerance)
+         self._instance_augmented_attributes[instance_name].extend(new_attrs)
 
    def iteration_k_solve(self):
@@ -1350 +995 @@
 
       solve_start_time = time.time()
+
+      # STEP -1: if using a PH solver manager, propagate current weights/averages to the appropriate solver servers.
+      #          ditto the tree node statistics, which are necessary if linearizing (so an optimization could be
+      #          performed here).
+      # NOTE: We aren't currently propagating rhos, as they generally don't change - we need to
+      #       have a flag, though, indicating whether the rhos have changed, so they can be
+      #       transmitted if needed.
+      if self._solver_manager_type == "ph":
+         self._transmit_weights_and_averages()
+         self._transmit_tree_node_statistics()
 
       # STEP 0: set up all global solver options.
@@ -1443 +1098 @@
             for objective_name in instance.active_components(Objective):
                objective = instance.active_components(Objective)[objective_name]
-               print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], 0.0)
+               print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], self._scenario_tree.compute_scenario_cost(instance))
 
    def solve(self):
@@ -1489 +1144 @@
       for plugin in self._ph_plugins:      
          plugin.post_iteration_0(self)
+
+      # if using a PH solver server, trasnsmit the rhos prior to the iteration
+      # k solve sequence. for now, we are assuming that the rhos don't change
+      # on a per-iteration basis, but that assumption can be easily relaxed.
+      # it is important to do this after the plugins have a chance to do their
+      # computation.
+      if self._solver_manager_type == "ph":
+         self._transmit_rhos()
+         self._enable_ph_objectives()
 
       # checkpoint if it's time - which it always is after iteration 0,

coopr.pysp/stable/2.3/coopr/pysp/phinit.py

@@ -531 +531 @@
       print "Writing EF for remainder problem"
       print ""
-      write_ef(ph._scenario_tree, ph._instances, os.path.expanduser(options.ef_output_file))
+      create_and_write_ef(ph._scenario_tree, ph._instances, os.path.expanduser(options.ef_output_file))
 
    #
@@ -643 +643 @@
     
     return ans
-

coopr.pysp/stable/2.3/coopr/pysp/phserver.py

@@ -6 +6 @@
 #  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
 #  the U.S. Government retains certain rights in this software.
-#  For more information, see the Coopr README.txt file.
+#  Fr more information, see the Coopr README.txt file.
 #  _________________________________________________________________________
 
@@ -23 +23 @@
 from coopr.opt.base import SolverFactory
 from coopr.pysp.scenariotree import *
+from phutils import *
+from phobjective import *
 
 # Pyro
 import Pyro.core
 import Pyro.naming
+from Pyro.errors import PyroError,NamingError
 
 # garbage collection control.
@@ -35 +38 @@
 import pstats
 
+import pickle
+
+# disable multi-threading. for a solver server, this
+# would truly be a bad idea, as you (1) likely have limited
+# solver licenses and (2) likely want to dedicate all compute
+# resources on this node to a single solve.
+Pyro.config.PYRO_MULTITHREADED=0
 
 #
@@ -41 +51 @@
 class PHSolverServer(object):
 
-   def __init__(self, scenario_instances, solver):
-
-      self._scenario_instances = scenario_instances
+   def __init__(self, scenario_instances, solver, scenario_tree):
+
+      # the obvious stuff!
+      self._instances = scenario_instances # a map from scenario name to pyomo model instance
       self._solver = solver
+      self._scenario_tree = scenario_tree
+
+      # the objective functions are modified throughout the course of PH iterations.
+      # save the original, as a baseline to modify in subsequent iterations. reserve
+      # the original objectives, for subsequent modification. 
+      self._original_objective_expression = {}
+      for instance_name, instance in self._instances.items():
+         objective_name = instance.active_components(Objective).keys()[0]
+         expr = instance.active_components(Objective)[objective_name]._data[None].expr
+         if isinstance(expr, Expression) is False:
+            expr = _IdentityExpression(expr)
+         self._original_objective_expression[instance_name] = expr
+
+      # the server is in one of two modes - solve the baseline instances, or
+      # those augmented with the PH penalty terms. default is standard.
+      # NOTE: This is currently a global flag for all scenarios handled
+      #       by this server - easy enough to extend if we want.
+      self._solving_standard_objective = True
+
+   def enable_standard_objective(self):
+
+#      print "RECEIVED REQUEST TO ENABLE STANDARD OBJECTIVE FOR ALL SCENARIOS"
+
+      self._solving_standard_objective = True
+
+   def enable_ph_objective(self):
+
+#      print "RECEIVED REQUEST TO ENABLE PH OBJECTIVE FOR ALL SCENARIOS"
+
+      self._solving_standard_objective = False
 
    def solve(self, scenario_name):
 
-      print "RECEIVED SOLVE REQUEST!"
-      print "SCENARIO TO SOLVE=",scenario_name
-      if scenario_name not in self._scenario_instances:
+#      print "RECEIVED REQUEST TO SOLVE SCENARIO INSTANCE="+scenario_name
+#      if self._solving_standard_objective is True:
+#         print "OBJECTIVE=STANDARD"
+#      else:
+#         print "OBJECTIVE=PH"
+         
+      if scenario_name not in self._instances:
          print "***ERROR: Requested instance to solve not in PH solver server instance collection!"
          return None
-      scenario_instance = self._scenario_instances[scenario_name]
-      print "SCENARIO INSTANCE=",scenario_instance
-      print "SUPPORTS MULTITHREADING=",Pyro.util.supports_multithreading()
-      print "SOLVING"
-      self._solver.solve(scenario_instance)
-      print "DONE WITH SOLVE"
-      return scenario_instance
-#      return result
-
+      scenario_instance = self._instances[scenario_name]
+
+      # form the desired objective, depending on the solve mode.
+      if self._solving_standard_objective is True:
+         form_standard_objective(scenario_name, scenario_instance, self._original_objective_expression[scenario_name], self._scenario_tree)         
+      else:
+         # TBD - command-line drive the various options (integer tolerance, breakpoint strategies, etc.)
+         form_ph_objective(scenario_name, scenario_instance, \
+                           self._original_objective_expression[scenario_name], self._scenario_tree, \
+                           False, False, False, 0, 0.00001)
+
+      # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
+      # don't do this, you won't see any chance in the output files as you vary the problem parameters!
+      # ditto for instance fixing!
+      scenario_instance.preprocess()
+
+      results = self._solver.solve(scenario_instance)
+
+      print "Successfully solved scenario instance="+scenario_name
+#      print "RESULTS:"
+#      results.write(num=1)
+      encoded_results = pickle.dumps(results)
+      
+      return encoded_results
+
+   def update_weights_and_averages(self, scenario_name, new_weights, new_averages):
+
+#      print "RECEIVED REQUEST TO UPDATE WEIGHTS AND AVERAGES FOR SCENARIO=",scenario_name
+
+      if scenario_name not in self._instances:
+         print "***ERROR: Received request to update weights for instance not in PH solver server instance collection!"
+         return None
+      scenario_instance = self._instances[scenario_name]
+      
+      for weight_update in new_weights:
+         
+         weight_index = weight_update._index
+
+         target_weight_parameter = getattr(scenario_instance, weight_update.name)
+
+         for index in weight_index:
+            target_weight_parameter[index] = weight_update[index]()
+
+      for average_update in new_averages:
+         
+         average_index = average_update._index
+
+         target_average_parameter = getattr(scenario_instance, average_update.name)
+
+         for index in average_index:
+            target_average_parameter[index] = average_update[index]()
+
+   def update_rhos(self, scenario_name, new_rhos):
+
+#      print "RECEIVED REQUEST TO UPDATE RHOS FOR SCENARIO=",scenario_name
+
+      if scenario_name not in self._instances:
+         print "***ERROR: Received request to update weights for instance not in PH solver server instance collection!"
+         return None
+      scenario_instance = self._instances[scenario_name]      
+
+      for rho_update in new_rhos:
+
+         rho_index = rho_update._index
+
+         target_rho_parameter = getattr(scenario_instance, rho_update.name)
+
+         for index in rho_index: 
+            target_rho_parameter[index] = rho_update[index]() # the value operator is crucial!
+
+   def update_tree_node_statistics(self, scenario_name, new_node_minimums, new_node_maximums):
+
+#      print "RECEIVED REQUEST TO UPDATE TREE NODE STATISTICS SCENARIO=",scenario_name
+
+      for tree_node_name, tree_node_minimums in new_node_minimums.items():
+
+         tree_node = self._scenario_tree._tree_node_map[tree_node_name]
+         tree_node._minimums = tree_node_minimums
+
+      for tree_node_name, tree_node_maximums in new_node_maximums.items():
+
+         tree_node = self._scenario_tree._tree_node_map[tree_node_name]
+         tree_node._maximums = tree_node_maximums         
+      
+            
 #
 # utility method to construct an option parser for ph arguments, to be
@@ -80 +201 @@
                      dest="scenarios",
                      default=[])
+   parser.add_option("--all-scenarios",
+                     help="Indicate that the server is responsible for solving all scenarios",
+                     action="store_true",
+                     dest="all_scenarios",
+                     default=False)
    parser.add_option("--report-solutions",
                      help="Always report PH solutions after each iteration. Enabled if --verbose is enabled. Default is False.",
@@ -151 +277 @@
                      type="int",
                      default=0)
-   parser.add_option("--enable-gc",
-                     help="Enable the python garbage collecter. Default is True.",
-                     action="store_true",
-                     dest="enable_gc",
-                     default=True)
+   parser.add_option("--disable-gc",
+                     help="Disable the python garbage collecter. Default is False.",
+                     action="store_true",
+                     dest="disable_gc",
+                     default=False)
 
    parser.usage=usage_string
@@ -182 +308 @@
          print ""
          print "***ERROR: Exiting test driver: No 'model' object created in module "+reference_model_filename
-         return None, None, None
+         return None, None, None, None
 
       if model_import.model is None:
          print ""
          print "***ERROR: Exiting test driver: 'model' object equals 'None' in module "+reference_model_filename
-         return None, None, None
+         return None, None, None, None
  
       reference_model = model_import.model
    except IOError:
       print "***ERROR: Failed to load scenario reference model from file="+reference_model_filename
-      return None, None, None
+      return None, None, None, None
 
    try:
@@ -201 +327 @@
    except IOError:
       print "***ERROR: Failed to load scenario reference instance data from file="+reference_instance_filename
-      return None, None, None
+      return None, None, None, None
 
    #
@@ -231 +357 @@
 #
 
-def run_server(options, scenario_instances, solver):
+def run_server(options, scenario_instances, solver, scenario_tree):
 
    start_time = time.time()
 
-   print "RUNNING SERVER NOW!"
-   Pyro.core.initServer(banner=1)   
-
-   print "NUMBER OF SCENARIO INSTANCES="+str(len(scenario_instances))
-
-   print "ATTEMPTING TO LOCATE NAME SERVER"
+   Pyro.core.initServer(banner=0)   
+
    locator = Pyro.naming.NameServerLocator()
-   ns = locator.getNS()
-   print "FOUND NAME SERVER"
+   ns = None
+   try:
+      ns = locator.getNS()
+   except Pyro.errors.NamingError:
+      raise RuntimeError, "PH solver server failed to locate Pyro name server"      
 
    solver_daemon = Pyro.core.Daemon()
    solver_daemon.useNameServer(ns)
 
-   daemon_object = PHSolverServer(scenario_instances, solver)
+   daemon_object = PHSolverServer(scenario_instances, solver, scenario_tree)
    delegator_object = Pyro.core.ObjBase()
    delegator_object.delegateTo(daemon_object)
@@ -255 +380 @@
    # each scenario processed by this daemon.
    for scenario_name, scenario_instance in scenario_instances.items():
-      solver_daemon.connect(delegator_object, scenario_name)
-
-   while True:
-      solver_daemon.handleRequests()
-      print "HANDLED A REQUEST!"
-
-   solver_daemon.disconnect()
-
-   print "DONE RUNNING SERVER!"
+      # first, see if it's already registered - if no, cull the old entry.
+      # NOTE: This is a hack, as the cleanup / disconnect code doesn't
+      #       seem to be invoked when the daemon is killed.
+      try:
+         ns.unregister(scenario_name)
+      except NamingError:
+         pass
+
+      try:
+         solver_daemon.connect(delegator_object, scenario_name)
+      except Pyro.errors.NamingError:
+         raise RuntimeError, "Entry in name server already exists for scenario="+scenario_name
+
+   try:
+      solver_daemon.requestLoop()
+   except KeyboardInterrupt, SystemExit:
+      pass
+
+   solver_daemon.disconnect(delegator_object)
+   solver_daemon.shutdown()
 
    end_time = time.time()
@@ -273 +409 @@
 def exec_server(options):
 
-   # we need the base model (not so much the reference instance, but that -
-   # can't hurt too much until proven otherwise, that is) to construct
+   # we need the base model (not so much the reference instance, but that 
+   # can't hurt too much - until proven otherwise, that is) to construct
    # the scenarios that this server is responsible for.
-   print "LOADING REFERENCE AND SCENARIO MODELS/INSTANCES"
    reference_model, reference_instance, scenario_tree, scenario_tree_instance = load_reference_and_scenario_models(options)
    if reference_model is None or reference_instance is None or scenario_tree is None:
@@ -284 +419 @@
    # construct the set of scenario instances based on the command-line.
    scenario_instances = {}
-   if len(options.scenarios) == 0:
-      print "***Unable to launch PH solver server - no scenario specified!"
+   if (options.all_scenarios is True) and (len(options.scenarios) > 0):
+      print "***WARNING: Both individual scenarios and all-scenarios were specified on the command-line; proceeding using all scenarios"
+
+   if (len(options.scenarios) == 0) and (options.all_scenarios is False):
+      print "***Unable to launch PH solver server - no scenario(s) specified!"
       return
 
-   for scenario_name in options.scenarios:
+   scenarios_to_construct = []
+   if options.all_scenarios is True:
+      scenarios_to_construct.extend(scenario_tree._scenario_map.keys())
+   else:
+      scenarios_to_construct.extend(options.scenarios)
+
+   for scenario_name in scenarios_to_construct:
       print "Creating instance for scenario="+scenario_name
       if scenario_tree.contains_scenario(scenario_name) is False:
@@ -294 +438 @@
          return
 
+      # create the baseline scenario instance
       scenario_instance = construct_scenario_instance(scenario_tree,
                                                       options.instance_directory,
@@ -301 +446 @@
       if scenario_instance is None:
          print "***Unable to launch PH solver server - failed to create instance for scenario="+scenario_name
+         return
+
+      # augment the instance with PH-specific parameters (weights, rhos, etc).
+      # TBD: The default rho of 1.0 is kind of bogus. Need to somehow propagate
+      #      this value and the linearization parameter as a command-line argument.
+      new_penalty_variable_names = create_ph_parameters(scenario_instance, scenario_tree, 1.0, False)
+
       scenario_instances[scenario_name] = scenario_instance
 
@@ -318 +470 @@
       solver._report_timing = True
 
-   print "SOLVING!"
-   results = solver.solve(scenario_instance)
-   print "TYPE OF RESULTS=",type(results)
-   results.write(num=1)
-#   print "RESULTS=",results
-   print "DONE SOLVING!"
-
    # spawn the daemon.
-   run_server(options, scenario_instances, solver)
+   run_server(options, scenario_instances, solver, scenario_tree)
 
 def run(args=None):
@@ -349 +494 @@
     # much sense - so it is disabled by default. Because: It drops
     # the run-time by a factor of 3-4 on bigger instances.
-    if options.enable_gc is False:
+    if options.disable_gc is True:
        gc.disable()
     else:

coopr.pysp/stable/2.3/coopr/pysp/phutils.py

@@ -224 +224 @@
 
    return scenario_instance
+
+
+# creates all PH parameters for a problem instance, given a scenario tree
+# (to identify the relevant variables), a default rho (simply for initialization),
+# and a boolean indicating if quadratic penalty terms are to be linearized.
+# returns a list of any created variables, specifically when linearizing -
+# this is useful to clean-up reporting.
+
+def create_ph_parameters(instance, scenario_tree, default_rho, linearizing_penalty_terms):
+
+   new_penalty_variable_names = []   
+   
+   # first, gather all unique variables referenced in any stage
+   # other than the last, independent of specific indices. this
+   # "gather" step is currently required because we're being lazy
+   # in terms of index management in the scenario tree - which
+   # should really be done in terms of sets of indices.
+   # NOTE: technically, the "instance variables" aren't really references
+   # to the variable in the instance - instead, the reference model. this
+   # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
+   instance_variables = {}
+   
+   for stage in scenario_tree._stages[:-1]:
+      
+      for (reference_variable, index_template, reference_indices) in stage._variables:
+            
+         if reference_variable.name not in instance_variables.keys():
+               
+            instance_variables[reference_variable.name] = reference_variable
+
+   # for each blended variable, create a corresponding ph weight and average parameter in the instance.
+   # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor
+   # in the grand scheme of things.
+
+   for (variable_name, reference_variable) in instance_variables.items():
+
+      # PH WEIGHT
+
+      new_w_index = reference_variable._index
+      new_w_parameter_name = "PHWEIGHT_"+reference_variable.name
+      # this bit of ugliness is due to Pyomo not correctly handling the Param construction
+      # case when the supplied index set consists strictly of None, i.e., the source variable
+      # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
+      new_w_parameter = None
+      if (len(new_w_index) is 1) and (None in new_w_index):
+         new_w_parameter = Param(name=new_w_parameter_name)
+      else:
+         new_w_parameter = Param(new_w_index,name=new_w_parameter_name)
+      setattr(instance,new_w_parameter_name,new_w_parameter)
+
+      # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference
+      # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not
+      # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the
+      # issue in any case for now. NOTE: I'm not sure this is still the case in Pyomo.
+      for index in new_w_index:
+         new_w_parameter[index] = 0.0
+
+      # PH AVG
+
+      new_avg_index = reference_variable._index 
+      new_avg_parameter_name = "PHAVG_"+reference_variable.name
+      new_avg_parameter = None
+      if (len(new_avg_index) is 1) and (None in new_avg_index):
+         new_avg_parameter = Param(name=new_avg_parameter_name)
+      else:
+         new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
+      setattr(instance,new_avg_parameter_name,new_avg_parameter)
+
+      for index in new_avg_index:
+         new_avg_parameter[index] = 0.0
+
+      # PH RHO
+
+      new_rho_index = reference_variable._index 
+      new_rho_parameter_name = "PHRHO_"+reference_variable.name
+      new_rho_parameter = None
+      if (len(new_avg_index) is 1) and (None in new_avg_index):
+         new_rho_parameter = Param(name=new_rho_parameter_name)
+      else:
+         new_rho_parameter = Param(new_rho_index,name=new_rho_parameter_name)
+      setattr(instance,new_rho_parameter_name,new_rho_parameter)
+
+      for index in new_rho_index:
+         new_rho_parameter[index] = default_rho
+
+      # PH LINEARIZED PENALTY TERM
+
+      if linearizing_penalty_terms > 0:
+
+         new_penalty_term_variable_index = reference_variable._index
+         new_penalty_term_variable_name = "PHQUADPENALTY_"+reference_variable.name
+         # this is a quadratic penalty term - the lower bound is 0!
+         new_penalty_term_variable = None
+         if (len(new_avg_index) is 1) and (None in new_avg_index):
+            new_penalty_term_variable = Var(name=new_penalty_term_variable_name, bounds=(0.0,None))                  
+         else:
+            new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None))
+         new_penalty_term_variable.construct()
+         setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable)
+         new_penalty_variable_names.append(new_penalty_term_variable_name)
+
+      # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY. 
+
+      # also controls whether weight updates proceed at any iteration.
+
+      new_blend_index = reference_variable._index
+      new_blend_parameter_name = "PHBLEND_"+reference_variable.name
+      new_blend_parameter = None
+      if (len(new_avg_index) is 1) and (None in new_avg_index):
+         new_blend_parameter = Param(name=new_blend_parameter_name, within=Binary)
+      else:
+         new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary)
+      setattr(instance,new_blend_parameter_name,new_blend_parameter)
+
+      for index in new_blend_index:
+         new_blend_parameter[index] = 1
+
+   return new_penalty_variable_names      

coopr.pysp/stable/2.3/coopr/pysp/scenariotree.py

@@ -11 +11 @@
 import sys
 import types
-from coopr.pyomo import *
 import copy
 import os.path
 import traceback
 
+from coopr.pyomo import *
 from phutils import *
 
@@ -23 +23 @@
 
    """
-   def __init__(self, *args, **kwds):
-
-      self._name = ""
-      self._stage = None
+   def __init__(self, name, conditional_probability, stage, reference_instance):
+
+      self._name = name
+      self._stage = stage
       self._parent = None
       self._children = [] # a collection of ScenarioTreeNodes
-      self._conditional_probability = None # conditional on parent
+      self._conditional_probability = conditional_probability # conditional on parent
       self._scenarios = [] # a collection of all Scenarios passing through this node in the tree
 
@@ -44 +44 @@
       self._maximums = {}
 
+      # solution (variable) values for this node. assumed to be distinct
+      # from self._averages, as the latter are not necessarily feasible.
+      # objects in the map are actual variables. 
+      self._solution = {}
+
+      # for each variable referenced in the stage, clone the variable
+      # for purposes of storing solutions. we are being wasteful in 
+      # terms copying indices that may not be referenced in the stage.
+      # this is something that we might revisit if space/performance
+      # is an issue (space is the most likely issue)
+      for variable, match_template, variable_indices in self._stage._variables:
+
+         # don't bother copying bounds for variables, as the values stored
+         # here are computed elsewhere - and that entity is responsible for
+         # ensuring feasibility. this also leaves room for specifying infeasible
+         # or partial solutions.       
+
+         new_variable_index = variable._index
+         new_variable_name = variable._name
+         new_variable = None
+         if (len(new_variable_index) is 1) and (None in new_variable_index):
+            new_variable = Var(name=new_variable_name)
+         else:
+            new_variable = Var(new_variable_index, name=new_variable_name)
+
+         self._solution[new_variable_name] = new_variable
+
    #
    # a utility to compute the cost of the current node plus the expected costs of child nodes.
    #
+
    def computeExpectedNodeCost(self, scenario_instance_map):
 
@@ -90 +118 @@
    """ Constructor
        Arguments:
-           scenarioinstance             - the reference scenario instance.
+           scenarioinstance     - the reference (deterministic) scenario instance.
            scenariotreeinstance - the pyomo model specifying all scenario tree (text) data.
            scenariobundlelist   - a list of scenario names to retain, i.e., cull the rest to create a reduced tree!
@@ -172 +200 @@
       # can't do a single pass because the objects may not exist.
       for tree_node_name in node_ids:
-         new_tree_node = ScenarioTreeNode()
-         new_tree_node._name = tree_node_name
+
+         if tree_node_name not in node_stage_ids:
+            raise ValueError, "No stage is assigned to tree node=" + tree_node._name
+
+         stage_name = node_stage_ids[tree_node_name].value
+         if stage_name not in self._stage_map.keys():
+            raise ValueError, "Unknown stage=" + stage_name + " assigned to tree node=" + tree_node._name
+
+         new_tree_node = ScenarioTreeNode(tree_node_name, 
+                                          node_probability_map[tree_node_name].value, 
+                                          self._stage_map[stage_name],
+                                          self._reference_instance)
 
          self._tree_nodes.append(new_tree_node)
          self._tree_node_map[tree_node_name] = new_tree_node
-
-         new_tree_node._conditional_probability = node_probability_map[tree_node_name].value
-
-         if tree_node_name not in node_stage_ids:
-            raise ValueError, "No stage is assigned to tree node=" + tree_node._name
-         else:
-            stage_name = node_stage_ids[new_tree_node._name].value
-            if stage_name not in self._stage_map.keys():
-               raise ValueError, "Unknown stage=" + stage_name + " assigned to tree node=" + tree_node._name
-            else:
-               new_tree_node._stage = self._stage_map[stage_name]
-               self._stage_map[stage_name]._tree_nodes.append(new_tree_node)
+         self._stage_map[stage_name]._tree_nodes.append(new_tree_node)
 
       # link up the tree nodes objects based on the child id sets.
@@ -347 +374 @@
    # is the indicated scenario in the tree?
    #
+
    def contains_scenario(self, name):
       return name in self._scenario_map.keys()
@@ -353 +381 @@
    # get the scenario object from the tree.
    #
+
    def get_scenario(self, name):
       return self._scenario_map[name]
+
+   #
+   # compute the scenario cost for the input instance, i.e.,
+   # the sum of all stage cost variables.
+   #
+
+   def compute_scenario_cost(self, instance):
+      aggregate_cost = 0.0
+      for stage in self._stages:
+         instance_cost_variable = instance.active_components(Var)[stage._cost_variable[0].name][stage._cost_variable[1]]()
+         aggregate_cost += instance_cost_variable
+      return aggregate_cost
 
    #
@@ -363 +404 @@
    # of the scenario tree structure.
    #
+
    def compress(self, scenario_bundle_list):
 
@@ -450 +492 @@
    # returns the root node of the scenario tree
    #
+
    def findRootNode(self):
 
@@ -461 +504 @@
    # the maximal length of identifiers in various categories.
    #
+
    def computeIdentifierMaxLengths(self):
 
@@ -471 +515 @@
    # a utility function to (partially, at the moment) validate a scenario tree
    #
+
    def validate(self):
 

coopr.pysp/stable/2.3/coopr/pysp/tests/unit/test_ph.py

@@ -18 +18 @@
 # Import the testing packages
 #
-import unittest
 import pyutilib.misc
-import pyutilib.th
+import pyutilib.th as unittest
 import pyutilib.subprocess
 import coopr.pysp
@@ -46 +45 @@
 
 #
-# Define a testing class, using the pyutilib.th.TestCase class.  This is
-# an extension of the unittest.TestCase class that adds additional testing
-# functions.
+# Define a testing class, using the unittest.TestCase class.
 #
-class TestPH(pyutilib.th.TestCase):
+class TestPH(unittest.TestCase):
 
     def cleanup(self):
@@ -61 +58 @@
            del sys.modules["ReferenceModel"]
 
+    @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available")
     def test_quadratic_farmer(self):
-
-        if cplex_available is False:
-           return
-
         pyutilib.misc.setup_redirect(current_directory+"farmer_quadratic.out")
         farmer_examples_dir = pysp_examples_dir + "farmer"
@@ -77 +71 @@
         self.failUnlessFileEqualsBaseline(current_directory+"farmer_quadratic.out",current_directory+"farmer_quadratic.baseline", filter=filter_time)
        
+    @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available")
     def test_linearized_farmer(self):
-
-        solver_string = ""
         if cplex_available:
             solver_string="cplex"
-        #elif glpk_available:
-            #solver_string="glpk"
-        else:
-           return
-
         pyutilib.misc.setup_redirect(current_directory+"farmer_linearized.out")
         farmer_examples_dir = pysp_examples_dir + "farmer"
@@ -98 +86 @@
         self.failUnlessFileEqualsBaseline(current_directory+"farmer_linearized.out",current_directory+"farmer_linearized.baseline", filter=filter_time)
 
+    @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available")
     def test_linearized_farmer_nodedata(self):
-
-        solver_string = ""
         if cplex_available:
             solver_string="cplex"
-        #elif glpk_available:
-            #solver_string="glpk"
-        else:
-           return
-
         pyutilib.misc.setup_redirect(current_directory+"farmer_linearized_nodedata.out")
         farmer_examples_dir = pysp_examples_dir + "farmer"
@@ -119 +101 @@
         self.failUnlessFileEqualsBaseline(current_directory+"farmer_linearized_nodedata.out",current_directory+"farmer_linearized_nodedata.baseline", filter=filter_time)        
 
+    @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available")
     def test_quadratic_sizes3(self):
 
         # Ignore this for now
         return
-
-        if cplex_available is False:
-           return
        
         pyutilib.misc.setup_redirect(current_directory+"sizes3_quadratic.out")
@@ -146 +126 @@
 
 if __name__ == "__main__":
-   unittest.main()
+    unittest.main()

coopr.pysp/stable/2.3/scripts/ph_test_client

@@ -4 +4 @@
 import Pyro.naming
 
-from coopr.opt.base import SolverFactory
-from coopr.opt import *
+import pickle
+
+import pyutilib.misc
+import pyutilib.component.core
+from coopr.opt.results import SolverResults
 
 import sys
@@ -38 +41 @@
 
 print "CALLING SOLVE METHOD"
-result = obj.solve(scenario_name)
+encoded_result = obj.solve(scenario_name)
+result = pickle.loads(encoded_result)
+print "SUCCESSFULLY OBTAINED RESULT!"
 
-print "RESULT=",result
+print "RESULTS:"
+result.write(num=1)
 
 print "DONE"