Changeset 2434


Ignore:
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:
11 edited
15 copied

Legend:

Unmodified
Added
Removed
  • coopr.pysp/stable/2.3

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

    r2317 r2434  
    2121from ef_writer_script import *
    2222from phinit import *
     23from phsolvermanager import *
     24from phobjective import *
    2325
    2426pyutilib.component.core.PluginGlobals.pop_env()
    25 
    2627
    2728try:
  • coopr.pysp/stable/2.3/coopr/pysp/ef.py

    r2391 r2434  
    55import traceback
    66import copy
     7import gc
    78
    89from coopr.pysp.scenariotree import *
     
    1516from coopr.pyomo.base.var import _VarValue, _VarBase
    1617
     18#
    1719# brain-dead utility for determing if there is a binary to write in the
    1820# composite model - need to know this, because CPLEX doesn't like empty
    1921# binary blocks in the LP file.
     22#
     23
    2024def binaries_present(master_model, scenario_instances):
    2125
     
    3438   return False
    3539
     40#
    3641# brain-dead utility for determing if there is a binary to write in the
    3742# composite model - need to know this, because CPLEX doesn't like empty
    3843# integer blocks in the LP file.
     44#
     45
    3946def integers_present(master_model, scenario_instances):
    4047
     
    5360   return False
    5461
     62#
     63# a routine to create the extensive form, given an input scenario tree and instances.
     64# IMPT: unlike scenario instances, the extensive form instance is *not* self-contained.
     65#       in particular, it has binding constraints that cross the binding instance and
     66#       the scenario instances. it is up to the caller to keep track of which scenario
     67#       instances are associated with the extensive form. this might be something we
     68#       encapsulate at some later time.
     69#
     70
     71def create_ef_instance(scenario_tree, scenario_instances,
     72                       verbose_output = False,
     73                       generate_weighted_cvar = False, cvar_weight = None, risk_alpha = None):
     74
     75   #
     76   # validate cvar options, if specified.
     77   #
     78   if generate_weighted_cvar is True:
     79      if (cvar_weight is None) or (cvar_weight <= 0.0):
     80         raise RuntimeError, "Weight of CVaR term must be >= 0.0 - value supplied="+str(cvar_weight)
     81      if (risk_alpha is None) or (risk_alpha <= 0.0) or (risk_alpha >= 1.0):
     82         raise RuntimeError, "CVaR risk alpha must be between 0 and 1, exclusive - value supplied="+str(risk_alpha)
     83
     84      if verbose_output is True:
     85         print "Writing CVaR weighted objective"
     86         print "CVaR term weight="+str(cvar_weight)
     87         print "CVaR alpha="+str(risk_alpha)
     88         print ""
     89
     90   #
     91   # create the new and empty binding instance.
     92   #
     93
     94   binding_instance = Model()
     95   binding_instance.name = "MASTER"
     96
     97   # walk the scenario tree - create variables representing the common values for all scenarios
     98   # associated with that node. the constraints will be created later. also create expected-cost
     99   # variables for each node, to be computed via constraints/objectives defined in a subsequent pass.
     100   # master variables are created for all nodes but those in the last stage. expected cost variables
     101   # are, for no particularly good reason other than easy coding, created for nodes in all stages.
     102   if verbose_output is True:
     103      print "Creating variables for master binding instance"
     104
     105   for stage in scenario_tree._stages:
     106
     107      # first loop is to create master (blended) variables across all stages but the last.
     108      for (stage_variable, index_template, stage_variable_indices) in stage._variables:
     109
     110         if verbose_output is True:
     111            print "Creating master variable and blending constraints for decision variable="+stage_variable.name+", indices="+str(index_template)
     112
     113         for tree_node in stage._tree_nodes:
     114
     115            if stage != scenario_tree._stages[-1]:     
     116
     117               master_variable_name = tree_node._name + "_" + stage_variable.name
     118
     119               # because there may be a single stage variable and multiple indices, check
     120               # for the existence of the variable at this node - if you don't, you'll
     121               # inadvertently over-write what was there previously!
     122               master_variable = None
     123               try:
     124                  master_variable = getattr(binding_instance, master_variable_name)
     125               except:
     126                  new_master_variable_index = stage_variable._index
     127                  new_master_variable = None
     128                  if (len(new_master_variable_index) is 1) and (None in new_master_variable_index):
     129                     new_master_variable = Var(name=stage_variable.name)
     130                  else:
     131                     new_master_variable = Var(new_master_variable_index, name=stage_variable.name)
     132                  new_master_variable.construct()
     133                  new_master_variable._model = binding_instance
     134                  setattr(binding_instance, master_variable_name, new_master_variable)
     135
     136                  master_variable = new_master_variable
     137
     138               for index in stage_variable_indices:
     139
     140                  is_used = True # until proven otherwise                     
     141                  for scenario in tree_node._scenarios:
     142                     instance = scenario_instances[scenario._name]
     143                     if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
     144                        is_used = False
     145
     146                  is_fixed = False # until proven otherwise
     147                  for scenario in tree_node._scenarios:
     148                     instance = scenario_instances[scenario._name]
     149                     if getattr(instance,stage_variable.name)[index].fixed is True:
     150                        is_fixed = True
     151
     152                  if (is_used is True) and (is_fixed is False):
     153                           
     154                     # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
     155                     # and because presolve/simplification is name-based, the names *have* to be different.
     156                     master_variable[index].var = master_variable
     157                     master_variable[index].name = tree_node._name + "_" + master_variable[index].name
     158
     159                     for scenario in tree_node._scenarios:
     160
     161                        scenario_instance = scenario_instances[scenario._name]
     162                        scenario_variable = getattr(scenario_instance, stage_variable.name)
     163                        new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
     164                        new_constraint = Constraint(name=new_constraint_name)
     165                        new_expr = master_variable[index] - scenario_variable[index]
     166                        new_constraint.add(None, (0.0, new_expr, 0.0))
     167                        new_constraint._model = binding_instance
     168                        setattr(binding_instance, new_constraint_name, new_constraint)
     169
     170      # the second loop is for creating the stage cost variable in each tree node.
     171      for tree_node in stage._tree_nodes:                       
     172
     173         # create a variable to represent the expected cost at this node -
     174         # the constraint to compute this comes later.
     175         expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
     176         expected_cost_variable = Var(name=expected_cost_variable_name)
     177         expected_cost_variable._model = binding_instance
     178         setattr(binding_instance, expected_cost_variable_name, expected_cost_variable)
     179
     180   # if we're generating the weighted CVaR objective term, create the
     181   # corresponding variable and the master CVaR eta variable.
     182   if generate_weighted_cvar is True:
     183
     184      root_node = scenario_tree._stages[0]._tree_nodes[0]
     185     
     186      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
     187      cvar_cost_variable = Var(name=cvar_cost_variable_name)
     188      setattr(binding_instance, cvar_cost_variable_name, cvar_cost_variable)
     189      cvar_cost_variable.construct()
     190
     191      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
     192      cvar_eta_variable = Var(name=cvar_eta_variable_name)
     193      setattr(binding_instance, cvar_eta_variable_name, cvar_eta_variable)     
     194      cvar_eta_variable.construct()
     195
     196   binding_instance.preprocess()
     197
     198   # ditto above for the (non-expected) cost variable.
     199   for stage in scenario_tree._stages:
     200
     201      (cost_variable,cost_variable_index) = stage._cost_variable
     202
     203      print "Creating master variable and blending constraints for cost variable="+cost_variable.name+", index="+str(cost_variable_index)
     204
     205      for tree_node in stage._tree_nodes:
     206
     207         # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)         
     208
     209         # this is undoubtedly wasteful, in that a cost variable
     210         # for each tree node is created with *all* indices.         
     211         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
     212         new_cost_variable_index = cost_variable._index
     213         new_cost_variable = None
     214         if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index):
     215            new_cost_variable = Var(name=new_cost_variable_name)
     216         else:
     217            new_cost_variable = Var(new_cost_variable_index, name=new_cost_variable_name)
     218         new_cost_variable.construct()
     219         new_cost_variable._model = binding_instance
     220         setattr(binding_instance, new_cost_variable_name, new_cost_variable)         
     221
     222         # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
     223         new_cost_variable[cost_variable_index].var = new_cost_variable
     224         if cost_variable_index is not None:
     225            # if the variable index is None, the variable is derived from a VarValue, so the
     226            # name gets updated automagically.
     227            new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
     228
     229         for scenario in tree_node._scenarios:
     230
     231            scenario_instance = scenario_instances[scenario._name]
     232            scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
     233            new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
     234            new_constraint = Constraint(name=new_constraint_name)
     235            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
     236            new_constraint.add(None, (0.0, new_expr, 0.0))
     237            new_constraint._model = binding_instance
     238            setattr(binding_instance, new_constraint_name, new_constraint)
     239
     240   # create the constraints for computing the master per-node cost variables,
     241   # i.e., the current node cost and the expected cost of the child nodes.
     242   # if the root, then the constraint is just the objective.
     243
     244   for stage in scenario_tree._stages:
     245
     246      (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
     247
     248      for tree_node in stage._tree_nodes:
     249
     250         node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
     251         node_expected_cost_variable = getattr(binding_instance, node_expected_cost_variable_name)
     252
     253         node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
     254         node_cost_variable = getattr(binding_instance, node_cost_variable_name)                       
     255           
     256         constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
     257
     258         for child_node in tree_node._children:
     259
     260            child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
     261            child_node_expected_cost_variable = getattr(binding_instance, child_node_expected_cost_variable_name)
     262            constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
     263
     264         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
     265         new_constraint = Constraint(name=new_constraint_name)
     266         new_constraint.add(None, (0.0, constraint_expr, 0.0))
     267         new_constraint._model = binding_instance                     
     268         setattr(binding_instance, new_constraint_name, new_constraint)
     269
     270         if tree_node._parent is None:
     271
     272            an_instance = scenario_instances[scenario_instances.keys()[0]]
     273            an_objective = an_instance.active_components(Objective)
     274            opt_sense = an_objective[an_objective.keys()[0]].sense
     275
     276            opt_expression = node_expected_cost_variable
     277
     278            if generate_weighted_cvar is True:
     279               cvar_cost_variable_name = "CVAR_COST_" + tree_node._name
     280               cvar_cost_variable = getattr(binding_instance, cvar_cost_variable_name)
     281               opt_expression += cvar_weight * cvar_cost_variable           
     282
     283            new_objective = Objective(name="MASTER", sense=opt_sense)
     284            new_objective._data[None].expr = node_expected_cost_variable
     285            setattr(binding_instance, "MASTER", new_objective)
     286
     287   # CVaR requires the addition of a variable per scenario to represent the cost excess,
     288   # and a constraint to compute the cost excess relative to eta. we also replicate (following
     289   # what we do for node cost variables) an eta variable for each scenario instance, and
     290   # require equality with the master eta variable via constraints.
     291   if generate_weighted_cvar is True:
     292     
     293      root_node = scenario_tree._stages[0]._tree_nodes[0]
     294
     295      master_cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
     296      master_cvar_eta_variable = getattr(binding_instance, master_cvar_eta_variable_name)
     297     
     298      for scenario_name in scenario_instances.keys():
     299         scenario_instance = scenario_instances[scenario_name]
     300
     301         # unique names are required because the presolve isn't
     302         # aware of the "owning" models for variables.
     303         cvar_excess_variable_name = "CVAR_EXCESS_"+scenario_name
     304         cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals)
     305         setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable)
     306         cvar_excess_variable.construct()
     307
     308         cvar_eta_variable_name = "CVAR_ETA"
     309         cvar_eta_variable = Var(name=cvar_eta_variable_name)
     310         setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable)
     311         cvar_eta_variable.construct()
     312
     313         compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS"
     314         compute_excess_constraint = Constraint(name=compute_excess_constraint_name)
     315         compute_excess_expression = cvar_excess_variable
     316         for node in scenario_tree._scenario_map[scenario_name]._node_list:
     317            (cost_variable, cost_variable_idx) = node._stage._cost_variable
     318            compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx]
     319         compute_excess_expression += cvar_eta_variable
     320         compute_excess_constraint.add(None, (0.0, compute_excess_expression, None))
     321         compute_excess_constraint._model = scenario_instance
     322         setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint)
     323
     324         eta_equality_constraint_name = "MASTER_ETA_EQUALITY_WITH_" + scenario_instance.name
     325         eta_equality_constraint = Constraint(name=eta_equality_constraint_name)
     326         eta_equality_expr = master_cvar_eta_variable - cvar_eta_variable
     327         eta_equality_constraint.add(None, (0.0, eta_equality_expr, 0.0))
     328         eta_equality_constraint._model = binding_instance
     329         setattr(binding_instance, eta_equality_constraint_name, eta_equality_constraint)
     330
     331      # add the constraint to compute the master CVaR variable value. iterate
     332      # over scenario instances to create the expected excess component first.
     333      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
     334      cvar_cost_variable = getattr(binding_instance, cvar_cost_variable_name)
     335      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
     336      cvar_eta_variable = getattr(binding_instance, cvar_eta_variable_name)
     337     
     338      cvar_cost_expression = cvar_cost_variable - cvar_eta_variable
     339     
     340      for scenario_name in scenario_instances.keys():
     341         scenario_instance = scenario_instances[scenario_name]
     342         scenario_probability = scenario_tree._scenario_map[scenario_name]._probability
     343
     344         scenario_excess_variable_name = "CVAR_EXCESS_"+scenario_name
     345         scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name)
     346
     347         cvar_cost_expression = cvar_cost_expression - (scenario_probability * scenario_excess_variable) / (1.0 - risk_alpha)
     348
     349      compute_cvar_cost_constraint_name = "COMPUTE_CVAR_COST"
     350      compute_cvar_cost_constraint = Constraint(name=compute_cvar_cost_constraint_name)
     351      compute_cvar_cost_constraint.add(None, (0.0, cvar_cost_expression, 0.0))
     352      compute_cvar_cost_constraint._model = binding_instance
     353      setattr(binding_instance, compute_cvar_cost_constraint_name, compute_cvar_cost_constraint)
     354
     355   # after mucking with instances, presolve to collect terms required prior to output.         
     356   # IMPT: Do the scenario instances first, as the master depends on variables in the scenarios.
     357   for scenario_name in scenario_instances.keys():
     358      scenario_instance = scenario_instances[scenario_name]   
     359      scenario_instance.preprocess()               
     360
     361   binding_instance.preprocess()
     362
     363   return binding_instance
     364
     365#
     366# write the EF binding instance and all sub-instances. currently only outputs the CPLEX LP file format.
     367#
     368
     369def write_ef(binding_instance, scenario_instances, output_filename):
     370
     371   # create the output file.
     372   problem_writer = cpxlp.ProblemWriter_cpxlp()
     373   output_file = open(output_filename,"w")
     374
     375   problem_writer._output_prefixes = True # we always want prefixes
     376
     377   ################################################################################################
     378   #### WRITE THE MASTER OBJECTIVE ################################################################
     379   ################################################################################################
     380
     381   # write the objective for the master binding instance.
     382   problem_writer._output_objectives = True
     383   problem_writer._output_constraints = False
     384   problem_writer._output_variables = False
     385
     386   print >>output_file, "\\ Begin objective block for master"
     387   problem_writer._print_model_LP(binding_instance, output_file)
     388   print >>output_file, "\\ End objective block for master"
     389   print >>output_file, ""
     390
     391   ################################################################################################
     392   #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
     393   ################################################################################################
     394
     395   print >>output_file, "s.t."
     396   print >>output_file, ""
     397   
     398   problem_writer._output_objectives = False
     399   problem_writer._output_constraints = True
     400   problem_writer._output_variables = False
     401
     402   print >>output_file, "\\ Begin constraint block for master"
     403   problem_writer._print_model_LP(binding_instance, output_file)
     404   print >>output_file, "\\ End constraint block for master",
     405   print >>output_file, ""
     406
     407   for scenario_name in scenario_instances.keys():
     408      instance = scenario_instances[scenario_name]
     409      print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
     410      problem_writer._print_model_LP(instance, output_file)
     411      print >>output_file, "\\ End constraint block for scenario",scenario_name
     412      print >>output_file, ""
     413
     414   ################################################################################################
     415   #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
     416   ################################################################################################
     417
     418   # write the variables for the master binding instance, and then for each scenario.
     419   print >>output_file, "bounds"
     420   print >>output_file, ""
     421   
     422   problem_writer._output_objectives = False
     423   problem_writer._output_constraints = False
     424   problem_writer._output_variables = True
     425
     426   # first step: write variable bounds
     427
     428   problem_writer._output_continuous_variables = True
     429   problem_writer._output_integer_variables = False
     430   problem_writer._output_binary_variables = False
     431
     432   print >>output_file, "\\ Begin variable bounds block for master"
     433   problem_writer._print_model_LP(binding_instance, output_file)
     434   print >>output_file, "\\ End variable bounds block for master"
     435   print >>output_file, ""
     436   
     437   for scenario_name in scenario_instances.keys():
     438      instance = scenario_instances[scenario_name]
     439      print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name
     440      problem_writer._print_model_LP(instance, output_file)
     441      print >>output_file, "\\ End variable bounds block for scenario",scenario_name
     442      print >>output_file, ""
     443
     444   # second step: write integer indicators.
     445
     446   problem_writer._output_continuous_variables = False
     447   problem_writer._output_integer_variables = True
     448
     449   if integers_present(binding_instance, scenario_instances) is True:
     450
     451      print >>output_file, "integer"
     452      print >>output_file, ""
     453
     454      print >>output_file, "\\ Begin integer variable block for master"
     455      problem_writer._print_model_LP(binding_instance, output_file)
     456      print >>output_file, "\\ End integer variable block for master"
     457      print >>output_file, ""
     458   
     459      for scenario_name in scenario_instances.keys():
     460         instance = scenario_instances[scenario_name]
     461         print >>output_file, "\\ Begin integer variable block for scenario",scenario_name
     462         problem_writer._print_model_LP(instance, output_file)
     463         print >>output_file, "\\ End integer variable block for scenario",scenario_name
     464         print >>output_file, ""
     465
     466   # third step: write binary indicators.
     467
     468   problem_writer._output_integer_variables = False
     469   problem_writer._output_binary_variables = True
     470
     471   if binaries_present(binding_instance, scenario_instances) is True:
     472
     473      print >>output_file, "binary"
     474      print >>output_file, ""
     475
     476      print >>output_file, "\\ Begin binary variable block for master"
     477      problem_writer._print_model_LP(binding_instance, output_file)
     478      print >>output_file, "\\ End binary variable block for master"
     479      print >>output_file, ""
     480   
     481      for scenario_name in scenario_instances.keys():
     482         instance = scenario_instances[scenario_name]
     483         print >>output_file, "\\ Begin binary variable block for scenario",scenario_name
     484         problem_writer._print_model_LP(instance, output_file)
     485         print >>output_file, "\\ End integer binary block for scenario",scenario_name
     486         print >>output_file, ""
     487
     488   # wrap up.
     489   print >>output_file, "end"
     490
     491   # clean up.
     492   output_file.close()
     493   
     494
     495#
    55496# the main extensive-form writer routine - including read of scenarios/etc.
    56 def write_ef_from_scratch(model_directory, instance_directory, output_filename, \
     497#
     498
     499def write_ef_from_scratch(model_directory, instance_directory, output_filename,
     500                          verbose_output,
    57501                          generate_weighted_cvar, cvar_weight, risk_alpha):
    58502
     
    61505   scenario_data_directory_name = instance_directory
    62506
    63    print "Initializing extensive form writer"
    64    print ""
    65 
    66    ################################################################################################
    67    #### INITIALIZATION ############################################################################
    68    ################################################################################################
    69 
    70    #
    71    # validate cvar options, if specified.
    72    #
    73    if generate_weighted_cvar is True:
    74       if cvar_weight <= 0.0:
    75          raise RuntimeError, "Weight of CVaR term must be >= 0.0 - value supplied="+str(cvar_weight)
    76       if (risk_alpha <= 0.0) or (risk_alpha >= 1.0):
    77          raise RuntimeError, "CVaR risk alpha must be between 0 and 1, exclusive - value supplied="+str(risk_alpha)
    78 
    79       print "Writing CVaR weighted objective"
    80       print "CVaR term weight="+str(cvar_weight)
    81       print "CVaR alpha="+str(risk_alpha)
     507   if verbose_output is True:
     508      print "Initializing extensive form writer"
    82509      print ""
    83    
     510
    84511   #
    85512   # create and populate the core model
     
    115542   except IOError:
    116543     
    117       print "***ERROR: Failed to load scenario reference instance data from file="+reference_instance_filename
     544      print "***ERROR: Failed to load scenario reference instance data from file="+reference_scenario_filename
    118545      return           
    119546
     
    135562   # print the input tree for validation/information purposes.
    136563   #
    137    scenario_tree.pprint()
     564   if verbose_output is True:
     565      scenario_tree.pprint()
    138566
    139567   #
     
    145573      sys.exit(1)
    146574   else:
    147       print "Scenario tree is valid!"
     575      if verbose_output is True:
     576         print "Scenario tree is valid!"
    148577   print ""
    149578
     
    152581   #
    153582
    154    instances = {}
     583   # the construction of instances takes little overhead in terms of
     584   # memory potentially lost in the garbage-collection sense (mainly
     585   # only that due to parsing and instance simplification/prep-processing).
     586   # to speed things along, disable garbage collection if it enabled in
     587   # the first place through the instance construction process.
     588   # IDEA: If this becomes too much for truly large numbers of scenarios,
     589   #       we could manually collect every time X instances have been created.
     590
     591   re_enable_gc = False
     592   if gc.isenabled() is True:
     593      re_enable_gc = True
     594      gc.disable()
     595
     596   scenario_instances = {}
    155597   
    156598   if scenario_tree._scenario_based_data == 1:
    157       print "Scenario-based instance initialization enabled"
     599      if verbose_output is True:
     600         print "Scenario-based instance initialization enabled"
    158601   else:
    159       print "Node-based instance initialization enabled"
     602      if verbose_output is True:
     603         print "Node-based instance initialization enabled"
    160604         
    161605   for scenario in scenario_tree._scenarios:
     
    168612         if scenario_tree._scenario_based_data == 1:
    169613            scenario_data_filename = scenario_data_directory_name + os.sep + scenario._name + ".dat"
    170 #            print "Data for scenario=" + scenario._name + " loads from file=" + scenario_data_filename
    171614            scenario_instance = master_scenario_model.create(scenario_data_filename)
    172615         else:
     
    176619            while current_node is not None:
    177620               node_data_filename = scenario_data_directory_name + os.sep + current_node._name + ".dat"
    178 #               print "Node data for scenario=" + scenario._name + " partially loading from file=" + node_data_filename
    179621               scenario_data.add(node_data_filename)
    180622               current_node = current_node._parent
     
    191633     
    192634      scenario_instance.preprocess()
    193       instances[scenario._name] = scenario_instance
    194 
    195    print ""
    196 
    197    ################################################################################################
    198    #### CREATE THE MASTER / BINDING INSTANCE ######################################################
    199    ################################################################################################
    200 
    201    master_binding_instance = Model()
    202    master_binding_instance.name = "MASTER"
    203 
    204    # walk the scenario tree - create variables representing the common values for all scenarios
    205    # associated with that node. the constraints will be created later. also create expected-cost
    206    # variables for each node, to be computed via constraints/objectives defined in a subsequent pass.
    207    # master variables are created for all nodes but those in the last stage. expected cost variables
    208    # are, for no particularly good reason other than easy coding, created for nodes in all stages.
    209    print "Creating variables for master binding instance"
    210 
    211    for stage in scenario_tree._stages:
    212 
    213       for (stage_variable, index_template, stage_variable_indices) in stage._variables:
    214 
    215          print "Creating master variable and blending constraints for decision variable=", stage_variable, ", indices=",stage_variable_indices
    216 
    217          for tree_node in stage._tree_nodes:
    218 
    219             if stage != scenario_tree._stages[-1]:     
    220 
    221                master_variable_name = tree_node._name + "_" + stage_variable.name
    222 
    223                # because there may be a single stage variable and multiple indices, check
    224                # for the existence of the variable at this node - if you don't, you'll
    225                # inadvertently over-write what was there previously!
    226                master_variable = None
    227                try:
    228                   master_variable = getattr(master_binding_instance, master_variable_name)
    229                except:
    230                   new_master_variable_index = stage_variable._index
    231                   new_master_variable = None
    232                   if (len(new_master_variable_index) is 1) and (None in new_master_variable_index):
    233                      new_master_variable = Var(name=stage_variable.name)
    234                   else:
    235                      new_master_variable = Var(new_master_variable_index, name=stage_variable.name)
    236                   new_master_variable.construct()
    237                   new_master_variable._model = master_binding_instance
    238                   setattr(master_binding_instance, master_variable_name, new_master_variable)
    239 
    240                   # TBD - TECHNICALLY, WE NEED TO COPY BOUNDS - BUT WE REALLY DON'T, AS THEY ARE ON THE PER-INSTNACE VARS!
    241 
    242                   master_variable = new_master_variable
    243 
    244                for index in stage_variable_indices:
    245 
    246                   is_used = True # until proven otherwise                     
    247                   for scenario in tree_node._scenarios:
    248                      instance = instances[scenario._name]
    249                      if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
    250                         is_used = False
    251 
    252                   is_fixed = False # until proven otherwise
    253                   for scenario in tree_node._scenarios:
    254                      instance = instances[scenario._name]
    255                      if getattr(instance,stage_variable.name)[index].fixed is True:
    256                         is_fixed = True
    257 
    258                   if (is_used is True) and (is_fixed is False):
    259                            
    260                      # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
    261                      # and because presolve/simplification is name-based, the names *have* to be different.
    262                      master_variable[index].var = master_variable
    263                      master_variable[index].name = tree_node._name + "_" + master_variable[index].name
    264 
    265                      for scenario in tree_node._scenarios:
    266 
    267                         scenario_instance = instances[scenario._name]
    268                         scenario_variable = getattr(scenario_instance, stage_variable.name)
    269                         new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
    270                         new_constraint = Constraint(name=new_constraint_name)
    271                         new_expr = master_variable[index] - scenario_variable[index]
    272                         new_constraint.add(None, (0.0, new_expr, 0.0))
    273                         new_constraint._model = master_binding_instance
    274                         setattr(master_binding_instance, new_constraint_name, new_constraint)
    275 
    276             # create a variable to represent the expected cost at this node -
    277             # the constraint to compute this comes later.
    278             expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
    279             expected_cost_variable = Var(name=expected_cost_variable_name)
    280             expected_cost_variable._model = master_binding_instance
    281             setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable)
    282 
    283    # if we're generating the weighted CVaR objective term, create the corresponding variable and
    284    # the master CVaR eta variable.
    285    if generate_weighted_cvar is True:
    286       root_node = scenario_tree._stages[0]._tree_nodes[0]
    287      
    288       cvar_cost_variable_name = "CVAR_COST_" + root_node._name
    289       cvar_cost_variable = Var(name=cvar_cost_variable_name)
    290       setattr(master_binding_instance, cvar_cost_variable_name, cvar_cost_variable)
    291       cvar_cost_variable.construct()
    292 
    293       cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
    294       cvar_eta_variable = Var(name=cvar_eta_variable_name)
    295       setattr(master_binding_instance, cvar_eta_variable_name, cvar_eta_variable)     
    296       cvar_eta_variable.construct()
    297 
    298    master_binding_instance.preprocess()
    299 
    300    # ditto above for the (non-expected) cost variable.
    301    for stage in scenario_tree._stages:
    302 
    303       (cost_variable,cost_variable_index) = stage._cost_variable
    304 
    305       print "Creating master variable and blending constraints for cost variable=", cost_variable, ", index=", cost_variable_index     
    306 
    307       for tree_node in stage._tree_nodes:
    308 
    309          # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)         
    310 
    311          # this is undoubtedly wasteful, in that a cost variable
    312          # for each tree node is created with *all* indices.         
    313          new_cost_variable_name = tree_node._name + "_" + cost_variable.name
    314          new_cost_variable_index = cost_variable._index
    315          new_cost_variable = None
    316          if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index):
    317             new_cost_variable = Var(name=new_cost_variable_name)
    318          else:
    319             new_cost_variable = Var(new_cost_variable_index, name=new_cost_variable_name)
    320          new_cost_variable.construct()
    321          new_cost_variable._model = master_binding_instance
    322          setattr(master_binding_instance, new_cost_variable_name, new_cost_variable)         
    323 
    324          # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
    325          new_cost_variable[cost_variable_index].var = new_cost_variable
    326          if cost_variable_index is not None:
    327             # if the variable index is None, the variable is derived from a VarValue, so the
    328             # name gets updated automagically.
    329             new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
    330 
    331          for scenario in tree_node._scenarios:
    332 
    333             scenario_instance = instances[scenario._name]
    334             scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
    335             new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
    336             new_constraint = Constraint(name=new_constraint_name)
    337             new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
    338             new_constraint.add(None, (0.0, new_expr, 0.0))
    339             new_constraint._model = master_binding_instance
    340             setattr(master_binding_instance, new_constraint_name, new_constraint)
    341 
    342    # create the constraints for computing the master per-node cost variables,
    343    # i.e., the current node cost and the expected cost of the child nodes.
    344    # if the root, then the constraint is just the objective.
    345 
    346    for stage in scenario_tree._stages:
    347 
    348       (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
    349 
    350       for tree_node in stage._tree_nodes:
    351 
    352          node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
    353          node_expected_cost_variable = getattr(master_binding_instance, node_expected_cost_variable_name)
    354 
    355          node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
    356          node_cost_variable = getattr(master_binding_instance, node_cost_variable_name)                       
    357            
    358          constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
    359 
    360          for child_node in tree_node._children:
    361 
    362             child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
    363             child_node_expected_cost_variable = getattr(master_binding_instance, child_node_expected_cost_variable_name)
    364             constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
    365 
    366          new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
    367          new_constraint = Constraint(name=new_constraint_name)
    368          new_constraint.add(None, (0.0, constraint_expr, 0.0))
    369          new_constraint._model = master_binding_instance                     
    370          setattr(master_binding_instance, new_constraint_name, new_constraint)
    371 
    372          if tree_node._parent is None:
    373 
    374             an_instance = instances[instances.keys()[0]]
    375             an_objective = an_instance.active_components(Objective)
    376             opt_sense = an_objective[an_objective.keys()[0]].sense
    377 
    378             opt_expression = node_expected_cost_variable
    379 
    380             if generate_weighted_cvar is True:
    381                cvar_cost_variable_name = "CVAR_COST_" + tree_node._name
    382                cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name)
    383                opt_expression += cvar_weight * cvar_cost_variable
    384 
    385             new_objective = Objective(name="MASTER", sense=opt_sense)
    386             new_objective._data[None].expr = opt_expression
    387             setattr(master_binding_instance, "MASTER", new_objective)
    388 
    389    # CVaR requires the addition of a variable per scenario to represent the cost excess,
    390    # and a constraint to compute the cost excess relative to eta. we also replicate (following
    391    # what we do for node cost variables) an eta variable for each scenario instance, and
    392    # require equality with the master eta variable via constraints.
    393    if generate_weighted_cvar is True:
    394      
    395       root_node = scenario_tree._stages[0]._tree_nodes[0]
    396 
    397       master_cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
    398       master_cvar_eta_variable = getattr(master_binding_instance, master_cvar_eta_variable_name)
    399      
    400       for scenario_name in instances.keys():
    401          scenario_instance = instances[scenario_name]
    402 
    403          # unique names are required because the presolve isn't
    404          # aware of the "owning" models for variables.
    405          cvar_excess_variable_name = "CVAR_EXCESS_"+scenario_name
    406          cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals)
    407          setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable)
    408          cvar_excess_variable.construct()
    409 
    410          cvar_eta_variable_name = "CVAR_ETA"
    411          cvar_eta_variable = Var(name=cvar_eta_variable_name)
    412          setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable)
    413          cvar_eta_variable.construct()
    414 
    415          compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS"
    416          compute_excess_constraint = Constraint(name=compute_excess_constraint_name)
    417          compute_excess_expression = cvar_excess_variable
    418          for node in scenario_tree._scenario_map[scenario_name]._node_list:
    419             (cost_variable, cost_variable_idx) = node._stage._cost_variable
    420             compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx]
    421          compute_excess_expression += cvar_eta_variable
    422          compute_excess_constraint.add(None, (0.0, compute_excess_expression, None))
    423          compute_excess_constraint._model = scenario_instance
    424          setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint)
    425 
    426          eta_equality_constraint_name = "MASTER_ETA_EQUALITY_WITH_" + scenario_instance.name
    427          eta_equality_constraint = Constraint(name=eta_equality_constraint_name)
    428          eta_equality_expr = master_cvar_eta_variable - cvar_eta_variable
    429          eta_equality_constraint.add(None, (0.0, eta_equality_expr, 0.0))
    430          eta_equality_constraint._model = master_binding_instance
    431          setattr(master_binding_instance, eta_equality_constraint_name, eta_equality_constraint)
    432 
    433       # add the constraint to compute the master CVaR variable value. iterate
    434       # over scenario instances to create the expected excess component first.
    435       cvar_cost_variable_name = "CVAR_COST_" + root_node._name
    436       cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name)
    437       cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
    438       cvar_eta_variable = getattr(master_binding_instance, cvar_eta_variable_name)
    439      
    440       cvar_cost_expression = cvar_cost_variable - cvar_eta_variable
    441      
    442       for scenario_name in instances.keys():
    443          scenario_instance = instances[scenario_name]
    444          scenario_probability = scenario_tree._scenario_map[scenario_name]._probability
    445 
    446          scenario_excess_variable_name = "CVAR_EXCESS_"+scenario_name
    447          scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name)
    448 
    449          cvar_cost_expression = cvar_cost_expression - (scenario_probability * scenario_excess_variable) / (1.0 - risk_alpha)
    450 
    451       compute_cvar_cost_constraint_name = "COMPUTE_CVAR_COST"
    452       compute_cvar_cost_constraint = Constraint(name=compute_cvar_cost_constraint_name)
    453       compute_cvar_cost_constraint.add(None, (0.0, cvar_cost_expression, 0.0))
    454       compute_cvar_cost_constraint._model = master_binding_instance
    455       setattr(master_binding_instance, compute_cvar_cost_constraint_name, compute_cvar_cost_constraint)
    456 
    457    # after mucking with instances, presolve to collect terms required prior to output.         
    458    # IMPT: Do the scenario instances first, as the master depends on variables in the scenarios.
    459    for scenario_name in instances.keys():
    460       scenario_instance = instances[scenario_name]   
    461       scenario_instance.preprocess()         
    462 
    463    master_binding_instance.preprocess()
    464 
    465    ################################################################################################
    466    #### WRITE THE COMPOSITE MODEL #################################################################
    467    ################################################################################################
     635      scenario_instances[scenario._name] = scenario_instance
     636
     637   if re_enable_gc is True:
     638      gc.enable()
     639
     640   print ""
     641
     642   binding_instance = create_ef_instance(scenario_tree, scenario_instances,
     643                                         verbose_output = verbose_output,
     644                                         generate_weighted_cvar = generate_weighted_cvar,
     645                                         cvar_weight = cvar_weight,
     646                                         risk_alpha = risk_alpha)   
    468647
    469648   print ""
    470649   print "Starting to write extensive form"
    471650
    472    # create the output file.
    473    problem_writer = cpxlp.ProblemWriter_cpxlp()
    474    output_file = open(output_filename,"w")
    475 
    476    problem_writer._output_prefixes = True # we always want prefixes
    477 
    478    ################################################################################################
    479    #### WRITE THE MASTER OBJECTIVE ################################################################
    480    ################################################################################################
    481 
    482    # write the objective for the master binding instance.
    483    problem_writer._output_objectives = True
    484    problem_writer._output_constraints = False
    485    problem_writer._output_variables = False
    486 
    487    print >>output_file, "\\ Begin objective block for master"
    488    problem_writer._print_model_LP(master_binding_instance, output_file)
    489    print >>output_file, "\\ End objective block for master"
    490    print >>output_file, ""
    491 
    492    ################################################################################################
    493    #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
    494    ################################################################################################
    495 
    496    print >>output_file, "s.t."
    497    print >>output_file, ""
    498    
    499    problem_writer._output_objectives = False
    500    problem_writer._output_constraints = True
    501    problem_writer._output_variables = False
    502 
    503    print >>output_file, "\\ Begin constraint block for master"
    504    problem_writer._print_model_LP(master_binding_instance, output_file)
    505    print >>output_file, "\\ End constraint block for master",
    506    print >>output_file, ""
    507 
    508    for scenario_name in instances.keys():
    509       instance = instances[scenario_name]
    510       print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
    511       problem_writer._print_model_LP(instance, output_file)
    512       print >>output_file, "\\ End constraint block for scenario",scenario_name
    513       print >>output_file, ""
    514 
    515    ################################################################################################
    516    #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
    517    ################################################################################################
    518 
    519    # write the variables for the master binding instance, and then for each scenario.
    520    print >>output_file, "bounds"
    521    print >>output_file, ""
    522    
    523    problem_writer._output_objectives = False
    524    problem_writer._output_constraints = False
    525    problem_writer._output_variables = True
    526 
    527    # first step: write variable bounds
    528 
    529    problem_writer._output_continuous_variables = True
    530    problem_writer._output_integer_variables = False
    531    problem_writer._output_binary_variables = False
    532 
    533    print >>output_file, "\\ Begin variable bounds block for master"
    534    problem_writer._print_model_LP(master_binding_instance, output_file)
    535    print >>output_file, "\\ End variable bounds block for master"
    536    print >>output_file, ""
    537 
    538    for scenario_name in instances.keys():
    539       instance = instances[scenario_name]
    540       print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name
    541       problem_writer._print_model_LP(instance, output_file)
    542       print >>output_file, "\\ End variable bounds block for scenario",scenario_name
    543       print >>output_file, ""
    544 
    545    # second step: write integer indicators.
    546 
    547    problem_writer._output_continuous_variables = False
    548    problem_writer._output_integer_variables = True
    549 
    550    if integers_present(master_binding_instance, instances) is True:
    551 
    552       print >>output_file, "integer"
    553       print >>output_file, ""
    554 
    555       print >>output_file, "\\ Begin integer variable block for master"
    556       problem_writer._print_model_LP(master_binding_instance, output_file)
    557       print >>output_file, "\\ End integer variable block for master"
    558       print >>output_file, ""
    559    
    560       for scenario_name in instances.keys():
    561          instance = instances[scenario_name]
    562          print >>output_file, "\\ Begin integer variable block for scenario",scenario_name
    563          problem_writer._print_model_LP(instance, output_file)
    564          print >>output_file, "\\ End integer variable block for scenario",scenario_name
    565          print >>output_file, ""
    566 
    567    # third step: write binary indicators.
    568 
    569    problem_writer._output_integer_variables = False
    570    problem_writer._output_binary_variables = True
    571 
    572    if binaries_present(master_binding_instance, instances) is True:
    573 
    574       print >>output_file, "binary"
    575       print >>output_file, ""
    576 
    577       print >>output_file, "\\ Begin binary variable block for master"
    578       problem_writer._print_model_LP(master_binding_instance, output_file)
    579       print >>output_file, "\\ End binary variable block for master"
    580       print >>output_file, ""
    581    
    582       for scenario_name in instances.keys():
    583          instance = instances[scenario_name]
    584          print >>output_file, "\\ Begin binary variable block for scenario",scenario_name
    585          problem_writer._print_model_LP(instance, output_file)
    586          print >>output_file, "\\ End integer binary block for scenario",scenario_name
    587          print >>output_file, ""
    588 
    589    # wrap up.
    590    print >>output_file, "end"
    591 
    592    # clean up.
    593    output_file.close()
     651   write_ef(binding_instance, scenario_instances, output_filename)   
    594652
    595653   print ""
     
    605663   print ""
    606664
    607 def write_ef(scenario_tree, instances, output_filename):
     665def create_and_write_ef(scenario_tree, scenario_instances, output_filename):
    608666
    609667   start_time = time.time()
    610668
    611    ################################################################################################
    612    #### CREATE THE MASTER / BINDING INSTANCE ######################################################
    613    ################################################################################################
    614 
    615    master_binding_instance = Model()
    616    master_binding_instance.name = "MASTER"
    617 
    618    # walk the scenario tree - create variables representing the common values for all scenarios
    619    # associated with that node. the constraints will be created later. also create expected-cost
    620    # variables for each node, to be computed via constraints/objectives defined in a subsequent pass.
    621    # master variables are created for all nodes but those in the last stage. expected cost variables
    622    # are, for no particularly good reason other than easy coding, created for nodes in all stages.
    623    print "Creating variables for master binding instance"
    624 
    625    for stage in scenario_tree._stages:
    626 
    627       for (stage_variable, index_template, stage_variable_indices) in stage._variables:
    628 
    629          print "Creating master variable and blending constraints for decision variable=", stage_variable, ", indices=", index_template
    630 
    631          for tree_node in stage._tree_nodes:
    632 
    633             if stage != scenario_tree._stages[-1]:     
    634 
    635                master_variable_name = tree_node._name + "_" + stage_variable.name
    636 
    637                # because there may be a single stage variable and multiple indices, check
    638                # for the existence of the variable at this node - if you don't, you'll
    639                # inadvertently over-write what was there previously!
    640                master_variable = None
    641                try:
    642                   master_variable = getattr(master_binding_instance, master_variable_name)
    643                except:
    644                   new_master_variable_index = stage_variable._index
    645                   new_master_variable = None
    646                   if (len(new_master_variable_index) is 1) and (None in new_master_variable_index):
    647                      new_master_variable = Var(name=stage_variable.name)
    648                   else:
    649                      new_master_variable = Var(new_master_variable_index, name=stage_variable.name)
    650                   new_master_variable.construct()
    651                   new_master_variable._model = master_binding_instance
    652                   setattr(master_binding_instance, master_variable_name, new_master_variable)
    653 
    654                   # TBD - TECHNICALLY, WE NEED TO COPY BOUNDS - BUT WE REALLY DON'T, AS THEY ARE ON THE PER-INSTNACE VARS!
    655 
    656                   master_variable = new_master_variable
    657 
    658                for index in stage_variable_indices:
    659 
    660                   is_used = True # until proven otherwise                     
    661                   for scenario in tree_node._scenarios:
    662                      instance = instances[scenario._name]
    663                      if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
    664                         is_used = False
    665 
    666                   is_fixed = False # until proven otherwise
    667                   for scenario in tree_node._scenarios:
    668                      instance = instances[scenario._name]
    669                      if getattr(instance,stage_variable.name)[index].fixed is True:
    670                         is_fixed = True
    671 
    672                   if (is_used is True) and (is_fixed is False):
    673                            
    674                      # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
    675                      # and because presolve/simplification is name-based, the names *have* to be different.
    676                      master_variable[index].var = master_variable
    677                      master_variable[index].name = tree_node._name + "_" + master_variable[index].name
    678 
    679                      for scenario in tree_node._scenarios:
    680 
    681                         scenario_instance = instances[scenario._name]
    682                         scenario_variable = getattr(scenario_instance, stage_variable.name)
    683                         new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
    684                         new_constraint = Constraint(name=new_constraint_name)
    685                         new_expr = master_variable[index] - scenario_variable[index]
    686                         new_constraint.add(None, (0.0, new_expr, 0.0))
    687                         new_constraint._model = master_binding_instance
    688                         setattr(master_binding_instance, new_constraint_name, new_constraint)
    689 
    690             # create a variable to represent the expected cost at this node -
    691             # the constraint to compute this comes later.
    692             expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
    693             expected_cost_variable = Var(name=expected_cost_variable_name)
    694             expected_cost_variable._model = master_binding_instance
    695             setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable)
    696 
    697    master_binding_instance.preprocess()
    698 
    699    # ditto above for the (non-expected) cost variable.
    700    for stage in scenario_tree._stages:
    701 
    702       (cost_variable,cost_variable_index) = stage._cost_variable
    703 
    704       print "Creating master variable and blending constraints for cost variable=", cost_variable, ", index=", cost_variable_index     
    705 
    706       for tree_node in stage._tree_nodes:
    707 
    708          new_cost_variable_name = tree_node._name + "_" + cost_variable.name
    709 
    710          # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)
    711 
    712          # this is undoubtedly wasteful, in that a cost variable
    713          # for each tree node is created with *all* indices.
    714          new_cost_variable_name = tree_node._name + "_" + cost_variable.name
    715          new_cost_variable_index = cost_variable._index
    716          new_cost_variable = None
    717          if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index):
    718             new_cost_variable = Var(name=new_cost_variable_name)
    719          else:
    720             new_cost_variable = Var(new_cost_variable_index, new_cost_variable_name)
    721          new_cost_variable.construct()
    722          new_cost_variable._model = master_binding_instance
    723          setattr(master_binding_instance, new_cost_variable_name, new_cost_variable)                 
    724 
    725          # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
    726          new_cost_variable[cost_variable_index].var = new_cost_variable
    727          if cost_variable_index is not None:
    728             # if the variable index is None, the variable is derived from a VarValue, so the
    729             # name gets updated automagically.
    730             new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
    731 
    732          for scenario in tree_node._scenarios:
    733 
    734             scenario_instance = instances[scenario._name]
    735             scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
    736             new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
    737             new_constraint = Constraint(name=new_constraint_name)
    738             new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
    739             new_constraint.add(None, (0.0, new_expr, 0.0))
    740             new_constraint._model = master_binding_instance
    741             setattr(master_binding_instance, new_constraint_name, new_constraint)
    742 
    743    # create the constraints for computing the master per-node cost variables,
    744    # i.e., the current node cost and the expected cost of the child nodes.
    745    # if the root, then the constraint is just the objective.
    746 
    747    for stage in scenario_tree._stages:
    748 
    749       (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
    750 
    751       for tree_node in stage._tree_nodes:
    752 
    753          node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
    754          node_expected_cost_variable = getattr(master_binding_instance, node_expected_cost_variable_name)
    755 
    756          node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
    757          node_cost_variable = getattr(master_binding_instance, node_cost_variable_name)                       
    758            
    759          constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
    760 
    761          for child_node in tree_node._children:
    762 
    763             child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
    764             child_node_expected_cost_variable = getattr(master_binding_instance, child_node_expected_cost_variable_name)
    765             constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
    766 
    767          new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
    768          new_constraint = Constraint(name=new_constraint_name)
    769          new_constraint.add(None, (0.0, constraint_expr, 0.0))
    770          new_constraint._model = master_binding_instance                     
    771          setattr(master_binding_instance, new_constraint_name, new_constraint)
    772 
    773          if tree_node._parent is None:
    774 
    775             an_instance = instances[instances.keys()[0]]
    776             an_objective = an_instance.active_components(Objective)
    777             opt_sense = an_objective[an_objective.keys()[0]].sense
    778 
    779             new_objective = Objective(name="MASTER", sense=opt_sense)
    780             new_objective._data[None].expr = node_expected_cost_variable
    781             setattr(master_binding_instance, "MASTER", new_objective)
    782 
    783    master_binding_instance.preprocess()
    784 
    785    ################################################################################################
    786    #### WRITE THE COMPOSITE MODEL #################################################################
    787    ################################################################################################
     669   binding_instance = create_ef_instance(scenario_tree, scenario_instances)
    788670
    789671   print ""
    790672   print "Starting to write extensive form"
    791673
    792    # create the output file.
    793    problem_writer = cpxlp.ProblemWriter_cpxlp()
    794    output_file = open(output_filename,"w")
    795 
    796    problem_writer._output_prefixes = True # we always want prefixes
    797 
    798    ################################################################################################
    799    #### WRITE THE MASTER OBJECTIVE ################################################################
    800    ################################################################################################
    801 
    802    # write the objective for the master binding instance.
    803    problem_writer._output_objectives = True
    804    problem_writer._output_constraints = False
    805    problem_writer._output_variables = False
    806 
    807    print >>output_file, "\\ Begin objective block for master"
    808    problem_writer._print_model_LP(master_binding_instance, output_file)
    809    print >>output_file, "\\ End objective block for master"
    810    print >>output_file, ""
    811 
    812    ################################################################################################
    813    #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
    814    ################################################################################################
    815 
    816    print >>output_file, "s.t."
    817    print >>output_file, ""
    818    
    819    problem_writer._output_objectives = False
    820    problem_writer._output_constraints = True
    821    problem_writer._output_variables = False
    822 
    823    print >>output_file, "\\ Begin constraint block for master"
    824    problem_writer._print_model_LP(master_binding_instance, output_file)
    825    print >>output_file, "\\ End constraint block for master",
    826    print >>output_file, ""
    827 
    828    for scenario_name in instances.keys():
    829       instance = instances[scenario_name]
    830       print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
    831       problem_writer._print_model_LP(instance, output_file)
    832       print >>output_file, "\\ End constraint block for scenario",scenario_name
    833       print >>output_file, ""
    834 
    835    ################################################################################################
    836    #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
    837    ################################################################################################
    838 
    839    # write the variables for the master binding instance, and then for each scenario.
    840    print >>output_file, "bounds"
    841    print >>output_file, ""
    842    
    843    problem_writer._output_objectives = False
    844    problem_writer._output_constraints = False
    845    problem_writer._output_variables = True
    846 
    847    # first step: write variable bounds
    848 
    849    problem_writer._output_continuous_variables = True
    850    problem_writer._output_integer_variables = False
    851    problem_writer._output_binary_variables = False
    852 
    853    print >>output_file, "\\ Begin variable bounds block for master"
    854    problem_writer._print_model_LP(master_binding_instance, output_file)
    855    print >>output_file, "\\ End variable bounds block for master"
    856    print >>output_file, ""
    857    
    858    for scenario_name in instances.keys():
    859       instance = instances[scenario_name]
    860       print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name
    861       problem_writer._print_model_LP(instance, output_file)
    862       print >>output_file, "\\ End variable bounds block for scenario",scenario_name
    863       print >>output_file, ""
    864 
    865    # second step: write integer indicators.
    866 
    867    problem_writer._output_continuous_variables = False
    868    problem_writer._output_integer_variables = True
    869 
    870    if integers_present(master_binding_instance, instances) is True:
    871 
    872       print >>output_file, "integer"
    873       print >>output_file, ""
    874 
    875       print >>output_file, "\\ Begin integer variable block for master"
    876       problem_writer._print_model_LP(master_binding_instance, output_file)
    877       print >>output_file, "\\ End integer variable block for master"
    878       print >>output_file, ""
    879    
    880       for scenario_name in instances.keys():
    881          instance = instances[scenario_name]
    882          print >>output_file, "\\ Begin integer variable block for scenario",scenario_name
    883          problem_writer._print_model_LP(instance, output_file)
    884          print >>output_file, "\\ End integer variable block for scenario",scenario_name
    885          print >>output_file, ""
    886 
    887    # third step: write binary indicators.
    888 
    889    problem_writer._output_integer_variables = False
    890    problem_writer._output_binary_variables = True
    891 
    892    if binaries_present(master_binding_instance, instances) is True:
    893 
    894       print >>output_file, "binary"
    895       print >>output_file, ""
    896 
    897       print >>output_file, "\\ Begin binary variable block for master"
    898       problem_writer._print_model_LP(master_binding_instance, output_file)
    899       print >>output_file, "\\ End binary variable block for master"
    900       print >>output_file, ""
    901    
    902       for scenario_name in instances.keys():
    903          instance = instances[scenario_name]
    904          print >>output_file, "\\ Begin binary variable block for scenario",scenario_name
    905          problem_writer._print_model_LP(instance, output_file)
    906          print >>output_file, "\\ End integer binary block for scenario",scenario_name
    907          print >>output_file, ""
    908 
    909    # wrap up.
    910    print >>output_file, "end"
    911 
    912    # clean up.
    913    output_file.close()
     674   write_ef(binding_instance, scenario_instances, output_filename)
    914675
    915676   print ""
  • coopr.pysp/stable/2.3/coopr/pysp/ef_writer_script.py

    r2317 r2434  
    2424
    2525#
    26 # Setup command-line options
     26# utility method to construct an option parser for ef writer arguments
    2727#
    2828
    29 parser = OptionParser()
    30 parser.add_option("--model-directory",
    31                   help="The directory in which all model (reference and scenario) definitions are stored. Default is \".\".",
    32                   action="store",
    33                   dest="model_directory",
    34                   type="string",
    35                   default=".")
    36 parser.add_option("--instance-directory",
    37                   help="The directory in which all instance (reference and scenario) definitions are stored. Default is \".\".",
    38                   action="store",
    39                   dest="instance_directory",
    40                   type="string",
    41                   default=".")
    42 parser.add_option("--generate-weighted-cvar",
    43                   help="Add a weighted CVaR term to the primary objective",
    44                   action="store_true",
    45                   dest="generate_weighted_cvar",
    46                   default=False)
    47 parser.add_option("--cvar-weight",
    48                   help="The weight associated with the CVaR term in the risk-weighted objective formulation.",
    49                   action="store",
    50                   dest="cvar_weight",
    51                   type="float",
    52                   default=1.0)
    53 parser.add_option("--risk-alpha",
    54                   help="The probability threshold associated with cvar (or any future) risk-oriented performance metrics.",
    55                   action="store",
    56                   dest="risk_alpha",
    57                   type="float",
    58                   default=0.95)
    59 parser.add_option("--output-file",
    60                   help="Specify the name of the extensive form output file",
    61                   action="store",
    62                   dest="output_file",
    63                   type="string",
    64                   default="efout.lp")
    65 parser.add_option("--profile",
    66                   help="Enable profiling of Python code.  The value of this option is the number of functions that are summarized.",
    67                   action="store",
    68                   dest="profile",
    69                   default=0)
    70 parser.usage="efwriter [options]"
     29def construct_ef_writer_options_parser(usage_string):
    7130
     31   parser = OptionParser()
     32   parser.add_option("--verbose",
     33                     help="Generate verbose output, beyond the usual status output. Default is False.",
     34                     action="store_true",
     35                     dest="verbose",
     36                     default=False)
     37   parser.add_option("--model-directory",
     38                     help="The directory in which all model (reference and scenario) definitions are stored. Default is \".\".",
     39                     action="store",
     40                     dest="model_directory",
     41                     type="string",
     42                     default=".")
     43   parser.add_option("--instance-directory",
     44                     help="The directory in which all instance (reference and scenario) definitions are stored. Default is \".\".",
     45                     action="store",
     46                     dest="instance_directory",
     47                     type="string",
     48                     default=".")
     49   parser.add_option("--generate-weighted-cvar",
     50                     help="Add a weighted CVaR term to the primary objective",
     51                     action="store_true",
     52                     dest="generate_weighted_cvar",
     53                     default=False)
     54   parser.add_option("--cvar-weight",
     55                     help="The weight associated with the CVaR term in the risk-weighted objective formulation.",
     56                     action="store",
     57                     dest="cvar_weight",
     58                     type="float",
     59                     default=1.0)
     60   parser.add_option("--risk-alpha",
     61                     help="The probability threshold associated with cvar (or any future) risk-oriented performance metrics.",
     62                     action="store",
     63                     dest="risk_alpha",
     64                     type="float",
     65                     default=0.95)
     66   parser.add_option("--output-file",
     67                     help="Specify the name of the extensive form output file",
     68                     action="store",
     69                     dest="output_file",
     70                     type="string",
     71                     default="efout.lp")
     72   parser.add_option("--profile",
     73                     help="Enable profiling of Python code.  The value of this option is the number of functions that are summarized.",
     74                     action="store",
     75                     dest="profile",
     76                     default=0)
     77   parser.add_option("--disable-gc",
     78                     help="Disable the python garbage collecter. Default is False.",
     79                     action="store_true",
     80                     dest="disable_gc",
     81                     default=False)
     82   parser.usage=usage_string
     83
     84   return parser
     85   
    7286def run_ef_writer(options, args):
    7387
     
    8498      risk_alpha = options.risk_alpha
    8599
    86    write_ef_from_scratch(os.path.expanduser(options.model_directory), os.path.expanduser(options.instance_directory), os.path.expanduser(options.output_file), \
     100   write_ef_from_scratch(os.path.expanduser(options.model_directory),
     101                         os.path.expanduser(options.instance_directory),
     102                         os.path.expanduser(options.output_file),
     103                         options.verbose,
    87104                         generate_weighted_cvar, cvar_weight, risk_alpha)
    88 
    89    return
    90105
    91106def run(args=None):
     
    96111    #
    97112
    98     # for a one-pass execution, garbage collection doesn't make
    99     # much sense - so I'm disabling it. Because: It drops the run-time
    100     # by a factor of 3-4 on bigger instances.
    101     gc.disable()       
    102 
    103113    #
    104114    # Parse command-line options.
    105115    #
    106116    try:
    107        (options, args) = parser.parse_args(args=args)
     117       options_parser = construct_ef_writer_options_parser("efwriter [options]")
     118       (options, args) = options_parser.parse_args(args=args)
    108119    except SystemExit:
    109120       # the parser throws a system exit if "-h" is specified - catch
    110121       # it to exit gracefully.
    111122       return
     123
     124    if options.disable_gc is True:
     125       gc.disable()
     126    else:
     127       gc.enable()
    112128
    113129    if options.profile > 0:
     
    140156
    141157    gc.enable()
     158   
    142159    return ans
    143160
  • coopr.pysp/stable/2.3/coopr/pysp/ph.py

    r2391 r2434  
    2727from scenariotree import *
    2828from phutils import *
     29from phobjective import *
    2930
    3031from pyutilib.component.core import ExtensionPoint
     
    3334
    3435class ProgressiveHedging(object):
    35 
    36    #
    37    # routine to compute linearization breakpoints uniformly between the bounds and the mean.
    38    #
    39 
    40    # IMPT: In general, the breakpoint computation codes can return a 2-list even if the lb equals
    41    #       the ub. This case happens quite often in real models (although typically lb=xvag=ub).
    42    #       See the code for constructing the pieces on how this case is handled in the linearization.
    43    
    44    def compute_uniform_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints_per_side):
    45      
    46       breakpoints = []
    47 
    48       # add the lower bound - the first breakpoint.
    49       breakpoints.append(lb)
    50 
    51       # determine the breakpoints to the left of the mean.
    52       left_step = (xavg - lb) / num_breakpoints_per_side
    53       current_x = lb
    54       for i in range(1,num_breakpoints_per_side+1):
    55          this_lb = current_x
    56          this_ub = current_x+left_step
    57          if (fabs(this_lb-lb) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
    58             breakpoints.append(this_lb)
    59          current_x += left_step           
    60 
    61       # add the mean - it's always a breakpoint. unless!
    62       # the lb or ub and the avg are the same.
    63       if (fabs(lb-xavg) > self._integer_tolerance) and (fabs(ub-xavg) > self._integer_tolerance):
    64          breakpoints.append(xavg)
    65 
    66       # determine the breakpoints to the right of the mean.
    67       right_step = (ub - xavg) / num_breakpoints_per_side
    68       current_x = xavg
    69       for i in range(1,num_breakpoints_per_side+1):
    70          this_lb = current_x
    71          this_ub = current_x+right_step
    72          if (fabs(this_ub-xavg) > self._integer_tolerance) and (fabs(this_ub-ub) > self._integer_tolerance):         
    73             breakpoints.append(this_ub)
    74          current_x += right_step
    75 
    76       # add the upper bound - the last breakpoint.
    77       # the upper bound should always be different than the lower bound (I say with some
    78       # hesitation - it really depends on what plugins are doing to modify the bounds dynamically).
    79       breakpoints.append(ub)
    80 
    81       return breakpoints
    82 
    83    #
    84    # routine to compute linearization breakpoints uniformly between the current node min/max bounds.
    85    #   
    86 
    87    def compute_uniform_between_nodestat_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints):
    88 
    89       breakpoints = []
    90 
    91       # add the lower bound - the first breakpoint.
    92       breakpoints.append(lb)
    93 
    94       # add the node-min - the second breakpoint. but only if it is different than the lower bound and the mean.
    95       if (fabs(node_min-lb) > self._integer_tolerance) and (fabs(node_min-xavg) > self._integer_tolerance):     
    96          breakpoints.append(node_min)
    97 
    98       step = (node_max - node_min) / num_breakpoints
    99       current_x = node_min
    100       for i in range(1,num_breakpoints+1):
    101          this_lb = current_x
    102          this_ub = current_x+step
    103          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):
    104             breakpoints.append(this_lb)
    105          current_x += step           
    106 
    107       # add the node-max - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
    108       if (fabs(node_max-ub) > self._integer_tolerance) and (fabs(node_max-xavg) > self._integer_tolerance):           
    109          breakpoints.append(node_max)     
    110 
    111       # add the upper bound - the last breakpoint.
    112       breakpoints.append(ub)
    113 
    114       # add the mean - it's always a breakpoint. unless! -
    115       # it happens to be equal to (within tolerance) the lower or upper bounds.
    116       # sort to insert it in the correct place.
    117       if (fabs(xavg - lb) > self._integer_tolerance) and (fabs(xavg - ub) > self._integer_tolerance):
    118          breakpoints.append(xavg)
    119       breakpoints.sort()
    120 
    121       return breakpoints
    122 
    123    #
    124    # routine to compute linearization breakpoints using "Woodruff" relaxation of the compute_uniform_between_nodestat_breakpoints.
    125    #   
    126 
    127    def compute_uniform_between_woodruff_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints):
    128 
    129       breakpoints = []
    130 
    131       # add the lower bound - the first breakpoint.
    132       breakpoints.append(lb)
    133 
    134       # be either three "differences" from the mean, or else "halfway to the bound", whichever is closer to the mean.
    135       left = max(xavg - 3.0 * (xavg - node_min), xavg - 0.5 * (xavg - lb))
    136       right = min(xavg + 3.0 * (node_max - xavg), xavg + 0.5 * (ub - xavg))     
    137 
    138       # add the left bound - the second breakpoint. but only if it is different than the lower bound and the mean.
    139       if (fabs(left-lb) > self._integer_tolerance) and (fabs(left-xavg) > self._integer_tolerance):     
    140          breakpoints.append(left)
    141 
    142       step = (right - left) / num_breakpoints
    143       current_x = left
    144       for i in range(1,num_breakpoints+1):
    145          this_lb = current_x
    146          this_ub = current_x+step
    147          if (fabs(this_lb-left) > self._integer_tolerance) and (fabs(this_lb-right) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
    148             breakpoints.append(this_lb)
    149          current_x += step           
    150 
    151       # add the right bound - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
    152       if (fabs(right-ub) > self._integer_tolerance) and (fabs(right-xavg) > self._integer_tolerance):           
    153          breakpoints.append(right)     
    154 
    155       # add the upper bound - the last breakpoint.
    156       breakpoints.append(ub)
    157 
    158       # add the mean - it's always a breakpoint.
    159       # sort to insert it in the correct place.
    160       breakpoints.append(xavg)
    161       breakpoints.sort()
    162 
    163       return breakpoints
    164 
    165    #
    166    # routine to compute linearization breakpoints based on an exponential distribution from the mean in each direction.
    167    #   
    168 
    169    def compute_exponential_from_mean_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints_per_side):
    170 
    171       breakpoints = []
    172 
    173       # add the lower bound - the first breakpoint.
    174       breakpoints.append(lb)
    175 
    176       # determine the breakpoints to the left of the mean.
    177       left_delta = xavg - lb
    178       base = exp(log(left_delta) / num_breakpoints_per_side)
    179       current_offset = base
    180       for i in range(1,num_breakpoints_per_side+1):
    181          current_x = xavg - current_offset
    182          if (fabs(current_x-lb) > self._integer_tolerance) and (fabs(current_x-xavg) > self._integer_tolerance):
    183             breakpoints.append(current_x)
    184          current_offset *= base
    185 
    186       # add the mean - it's always a breakpoint.
    187       breakpoints.append(xavg)
    188 
    189       # determine the breakpoints to the right of the mean.
    190       right_delta = ub - xavg
    191       base = exp(log(right_delta) / num_breakpoints_per_side)
    192       current_offset = base
    193       for i in range(1,num_breakpoints_per_side+1):
    194          current_x = xavg + current_offset
    195          if (fabs(current_x-xavg) > self._integer_tolerance) and (fabs(current_x-ub) > self._integer_tolerance):         
    196             breakpoints.append(current_x)
    197          current_offset *= base
    198 
    199       # add the upper bound - the last breakpoint.
    200       breakpoints.append(ub)         
    201 
    202       return breakpoints     
    20336
    20437   #
     
    379212      return (num_fixed_discrete_vars, num_fixed_continuous_vars)
    380213
    381    # a utility to create piece-wise linear constraint expressions for a given variable, for
    382    # use in constructing the augmented (penalized) PH objective. lb and ub are the bounds
    383    # on this piece, variable is the actual instance variable, and average is the instance
    384    # parameter specifying the average of this variable across instances sharing passing
    385    # through a common tree node. lb and ub are floats.
    386    # IMPT: There are cases where lb=ub, in which case the slope is 0 and the intercept
    387    #       is simply the penalty at the lower(or upper) bound.
    388    def _create_piecewise_constraint_expression(self, lb, ub, instance_variable, variable_average, quad_variable):
    389 
    390       penalty_at_lb = (lb - variable_average()) * (lb - variable_average())
    391       penalty_at_ub = (ub - variable_average()) * (ub - variable_average())
    392       slope = None
    393       if fabs(ub-lb) > self._integer_tolerance:
    394          slope = (penalty_at_ub - penalty_at_lb) / (ub - lb)
    395       else:
    396          slope = 0.0
    397       intercept = penalty_at_lb - slope * lb
    398       expression = (0.0, quad_variable - slope * instance_variable - intercept, None)
    399 
    400       return expression
    401 
    402214   # when the quadratic penalty terms are approximated via piecewise linear segments,
    403215   # we end up (necessarily) "littering" the scenario instances with extra constraints.
     
    424236      for (instance_name, instance) in self._instances.items():
    425237
    426          # first, gather all unique variables referenced in any stage
    427          # other than the last, independent of specific indices. this
    428          # "gather" step is currently required because we're being lazy
    429          # in terms of index management in the scenario tree - which
    430          # should really be done in terms of sets of indices.
    431          # NOTE: technically, the "instance variables" aren't really references
    432          # to the variable in the instance - instead, the reference model. this
    433          # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
    434          instance_variables = {}
    435          
    436          for stage in self._scenario_tree._stages[:-1]:
    437            
    438             for (reference_variable, index_template, reference_indices) in stage._variables:
    439                  
    440                if reference_variable.name not in instance_variables.keys():
    441                      
    442                   instance_variables[reference_variable.name] = reference_variable
    443 
    444          # for each blended variable, create a corresponding ph weight and average parameter in the instance.
    445          # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor
    446          # in the grand scheme of things.
    447 
    448          for (variable_name, reference_variable) in instance_variables.items():
    449 
    450             # PH WEIGHT
    451 
    452             new_w_index = reference_variable._index
    453             new_w_parameter_name = "PHWEIGHT_"+reference_variable.name
    454             # this bit of ugliness is due to Pyomo not correctly handling the Param construction
    455             # case when the supplied index set consists strictly of None, i.e., the source variable
    456             # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
    457             new_w_parameter = None
    458             if (len(new_w_index) is 1) and (None in new_w_index):
    459                new_w_parameter = Param(name=new_w_parameter_name)
    460             else:
    461                new_w_parameter = Param(new_w_index,name=new_w_parameter_name)
    462             setattr(instance,new_w_parameter_name,new_w_parameter)
    463 
    464             # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference
    465             # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not
    466             # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the
    467             # issue in any case for now. NOTE: I'm not sure this is still the case in Pyomo.
    468             for index in new_w_index:
    469                new_w_parameter[index] = 0.0
    470 
    471             # PH AVG
    472 
    473             new_avg_index = reference_variable._index
    474             new_avg_parameter_name = "PHAVG_"+reference_variable.name
    475             new_avg_parameter = None
    476             if (len(new_avg_index) is 1) and (None in new_avg_index):
    477                new_avg_parameter = Param(name=new_avg_parameter_name)
    478             else:
    479                new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
    480             setattr(instance,new_avg_parameter_name,new_avg_parameter)
    481 
    482             for index in new_avg_index:
    483                new_avg_parameter[index] = 0.0
    484 
    485             # PH RHO
    486 
    487             new_rho_index = reference_variable._index
    488             new_rho_parameter_name = "PHRHO_"+reference_variable.name
    489             new_rho_parameter = None
    490             if (len(new_avg_index) is 1) and (None in new_avg_index):
    491                new_rho_parameter = Param(name=new_rho_parameter_name)
    492             else:
    493                new_rho_parameter = Param(new_rho_index,name=new_rho_parameter_name)
    494             setattr(instance,new_rho_parameter_name,new_rho_parameter)
    495 
    496             for index in new_rho_index:
    497                new_rho_parameter[index] = self._rho
    498 
    499             # PH LINEARIZED PENALTY TERM
    500 
    501             if self._linearize_nonbinary_penalty_terms > 0:
    502 
    503                new_penalty_term_variable_index = reference_variable._index
    504                new_penalty_term_variable_name = "PHQUADPENALTY_"+reference_variable.name
    505                # this is a quadratic penalty term - the lower bound is 0!
    506                new_penalty_term_variable = None
    507                if (len(new_avg_index) is 1) and (None in new_avg_index):
    508                   new_penalty_term_variable = Var(name=new_penalty_term_variable_name, bounds=(0.0,None))                 
    509                else:
    510                   new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None))
    511                new_penalty_term_variable.construct()
    512                setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable)
    513                self._instance_augmented_attributes[instance_name].append(new_penalty_term_variable_name)
    514 
    515             # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY.
    516 
    517             # also controls whether weight updates proceed at any iteration.
    518 
    519             new_blend_index = reference_variable._index
    520             new_blend_parameter_name = "PHBLEND_"+reference_variable.name
    521             new_blend_parameter = None
    522             if (len(new_avg_index) is 1) and (None in new_avg_index):
    523                new_blend_parameter = Param(name=new_blend_parameter_name, within=Binary)
    524             else:
    525                new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary)
    526             setattr(instance,new_blend_parameter_name,new_blend_parameter)
    527 
    528             for index in new_blend_index:
    529                new_blend_parameter[index] = 1
     238         new_penalty_variable_names = create_ph_parameters(instance, self._scenario_tree, self._rho, self._linearize_nonbinary_penalty_terms)
     239         self._instance_augmented_attributes[instance_name].extend(new_penalty_variable_names)
    530240
    531241   # a simple utility to extract the first-stage cost statistics, e.g., min, average, and max.
     
    554264               minimum_value = this_value
    555265      return minimum_value, sum_values/num_values, maximum_value
    556      
     266
     267   # a utility to transmit - across the PH solver manager - the current weights
     268   # and averages for each of my problem instances. used to set up iteration K solves.
     269   def _transmit_weights_and_averages(self):
     270
     271      for scenario_name, scenario_instance in self._instances.items():
     272
     273         weights_to_transmit = []
     274         averages_to_transmit = []
     275
     276         for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.
     277
     278            for (variable, index_template, variable_indices) in stage._variables:
     279
     280               variable_name = variable.name
     281               weight_parameter_name = "PHWEIGHT_"+variable_name
     282               weights_to_transmit.append(getattr(scenario_instance, weight_parameter_name))
     283               average_parameter_name = "PHAVG_"+variable_name
     284               averages_to_transmit.append(getattr(scenario_instance, average_parameter_name))               
     285
     286         self._solver_manager.transmit_weights_and_averages(scenario_instance, weights_to_transmit, averages_to_transmit)
     287
     288   #
     289   # a utility to transmit - across the PH solver manager - the current rho values
     290   # for each of my problem instances. mainly after PH iteration 0 is complete,
     291   # in preparation for the iteration K solves.
     292   #
     293
     294   def _transmit_rhos(self):
     295
     296      for scenario_name, scenario_instance in self._instances.items():
     297
     298         rhos_to_transmit = []
     299
     300         for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no rhos to worry about.
     301
     302            for (variable, index_template, variable_indices) in stage._variables:
     303
     304               variable_name = variable.name
     305               rho_parameter_name = "PHRHO_"+variable_name
     306               rhos_to_transmit.append(getattr(scenario_instance, rho_parameter_name))
     307
     308         self._solver_manager.transmit_rhos(scenario_instance, rhos_to_transmit)
     309
     310   #
     311   # a utility to transmit - across the PH solver manager - the current scenario
     312   # tree node statistics to each of my problem instances. done prior to each
     313   # PH iteration k.
     314   #
     315
     316   def _transmit_tree_node_statistics(self):
     317
     318      for scenario_name, scenario_instance in self._instances.items():
     319
     320         tree_node_minimums = {}
     321         tree_node_maximums = {}
     322
     323         scenario = self._scenario_tree._scenario_map[scenario_name]
     324
     325         for tree_node in scenario._node_list:
     326
     327            tree_node_minimums[tree_node._name] = tree_node._minimums
     328            tree_node_maximums[tree_node._name] = tree_node._maximums
     329
     330         self._solver_manager.transmit_tree_node_statistics(scenario_instance, tree_node_minimums, tree_node_maximums)         
     331
     332   #
     333   # a utility to enable - across the PH solver manager - weighted penalty objectives.
     334   #
     335
     336   def _enable_ph_objectives(self):
     337
     338      for scenario_name, scenario_instance in self._instances.items():
     339
     340         self._solver_manager.enable_ph_objective(scenario_instance)
     341         
    557342
    558343   """ Constructor
     
    596381      self._solver_manager_type = "serial" # serial or pyro are the options currently available
    597382     
    598       self._solver = None # will eventually be unnecessary once Bill eliminates the need for a solver in the solver manager constructor.
     383      self._solver = None
    599384      self._solver_manager = None
    600385     
     
    1187972   def form_iteration_k_objectives(self):
    1188973
    1189       # for each blended variable (i.e., those not appearing in the final stage),
    1190       # add the linear and quadratic penalty terms to the objective.
    1191974      for instance_name, instance in self._instances.items():
    1192          
    1193          objective_name = instance.active_components(Objective).keys()[0]         
    1194          objective = instance.active_components(Objective)[objective_name]
    1195          # clone the objective, because we're going to augment (repeatedly) the original objective.
    1196          objective_expression = self._original_objective_expression[instance_name].clone()
    1197          # the quadratic expression is really treated as just a list - eventually should be treated as a full expression.
    1198          quad_expression = 0.0
    1199          
    1200          for stage in self._scenario_tree._stages[:-1]: # skip the last stage, as no blending occurs
    1201 
    1202             variable_tree_node = None
    1203             for node in stage._tree_nodes:
    1204                for scenario in node._scenarios:
    1205                   if scenario._name == instance_name:
    1206                      variable_tree_node = node
    1207                      break
    1208 
    1209             for (reference_variable, index_template, variable_indices) in stage._variables:
    1210 
    1211                variable_name = reference_variable.name
    1212                variable_type = reference_variable.domain
    1213 
    1214                w_parameter_name = "PHWEIGHT_"+variable_name
    1215                w_parameter = instance.active_components(Param)[w_parameter_name]
    1216                  
    1217                average_parameter_name = "PHAVG_"+variable_name
    1218                average_parameter = instance.active_components(Param)[average_parameter_name]
    1219 
    1220                rho_parameter_name = "PHRHO_"+variable_name
    1221                rho_parameter = instance.active_components(Param)[rho_parameter_name]
    1222 
    1223                blend_parameter_name = "PHBLEND_"+variable_name
    1224                blend_parameter = instance.active_components(Param)[blend_parameter_name]
    1225 
    1226                node_min_parameter = variable_tree_node._minimums[variable_name]
    1227                node_max_parameter = variable_tree_node._maximums[variable_name]               
    1228 
    1229                quad_penalty_term_variable = None
    1230                if self._linearize_nonbinary_penalty_terms > 0:
    1231                   quad_penalty_term_variable_name = "PHQUADPENALTY_"+variable_name
    1232                   quad_penalty_term_variable = instance.active_components(Var)[quad_penalty_term_variable_name]
    1233 
    1234                instance_variable = instance.active_components(Var)[variable_name]
    1235 
    1236                for index in variable_indices:
    1237 
    1238                   if (instance_variable[index].status is not VarStatus.unused) and (instance_variable[index].fixed is False):
    1239 
    1240                      # add the linear (w-weighted) term is a consistent fashion, independent of variable type.
    1241                      # if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes for reasons currently unknown.
    1242                      objective_expression += (w_parameter[index] * instance_variable[index])
    1243 
    1244                      # there are some cases in which a user may want to avoid the proximal term completely.
    1245                      # it's of course only a good idea when there are at least bounds (both lb and ub) on
    1246                      # the variables to be blended.
    1247                      if self._drop_proximal_terms is False:
    1248 
    1249                         # deal with binaries
    1250                         if isinstance(variable_type, BooleanSet) is True:
    1251 
    1252                            if self._retain_quadratic_binary_terms is False:
    1253                               # this rather ugly form of the linearized quadratic expression term is required
    1254                               # due to a pyomo bug - the form (rho/2) * (x+y+z) chokes in presolve when distributing
    1255                               # over the sum.
    1256                               new_term = (blend_parameter[index] * rho_parameter[index] / 2.0 * instance_variable[index]) - \
    1257                                          (blend_parameter[index] * rho_parameter[index] * average_parameter[index] * instance_variable[index]) + \
    1258                                          (blend_parameter[index] * rho_parameter[index] / 2.0 * average_parameter[index] * average_parameter[index])                             
    1259                               if objective.sense is minimize:
    1260                                  objective_expression += new_term
    1261                               else:
    1262                                  objective_expression -= new_term                                 
    1263                            else:
    1264                               quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)
    1265 
    1266                         # deal with everything else
    1267                         else:
    1268 
    1269                            if self._linearize_nonbinary_penalty_terms > 0:
    1270 
    1271                               # the variables are easy - just a single entry.
    1272                               if objective.sense is minimize:
    1273                                  objective_expression += (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
    1274                               else:
    1275                                  objective_expression -= (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
    1276                                  
    1277                               # the constraints are somewhat nastier.
    1278 
    1279                               # TBD - DEFINE CONSTRAINTS ON-THE-FLY?? (INDIVIDUALLY NAMED FOR NOW - CREATE INDEX SETS!) - OR A LEAST AN INDEX SET PER "PIECE"
    1280                               xavg = average_parameter[index]
    1281                               x = instance_variable[index]
    1282 
    1283                               lb = None
    1284                               ub = None
    1285 
    1286                               if x.lb is None:
    1287                                  raise ValueError, "No lower bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
    1288                               else:
    1289                                  lb = x.lb()
    1290 
    1291                               if x.ub is None:
    1292                                  raise ValueError, "No upper bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
    1293                               else:
    1294                                  ub = x.ub()
    1295 
    1296                               node_min = node_min_parameter[index]()
    1297                               node_max = node_max_parameter[index]()
    1298 
    1299                               # compute the breakpoint sequence according to the specified strategy.
    1300                               breakpoints = []
    1301                               if self._breakpoint_strategy == 0:
    1302                                  breakpoints = self.compute_uniform_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
    1303                               elif self._breakpoint_strategy == 1:
    1304                                  breakpoints = self.compute_uniform_between_nodestat_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
    1305                               elif self._breakpoint_strategy == 2:
    1306                                  breakpoints = self.compute_uniform_between_woodruff_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
    1307                               elif self._breakpoint_strategy == 3:
    1308                                  breakpoints = self.compute_exponential_from_mean_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)                                                                                                   
    1309                               else:
    1310                                  raise ValueError, "A breakpoint distribution strategy="+str(self._breakpoint_strategy)+" is currently not supported within PH!"
    1311 
    1312                               for i in range(0,len(breakpoints)-1):
    1313 
    1314                                  this_lb = breakpoints[i]
    1315                                  this_ub = breakpoints[i+1]
    1316 
    1317                                  piece_constraint_name = "QUAD_PENALTY_PIECE_"+str(i)+"_"+variable_name+str(index)
    1318                                  if hasattr(instance, piece_constraint_name) is False:
    1319                                     # this is the first time the constraint is being added - add it to the list of PH-specific constraints for this instance.
    1320                                     self._instance_augmented_attributes[instance_name].append(piece_constraint_name)                                                                           
    1321                                  piece_constraint = Constraint(name=piece_constraint_name)
    1322                                  piece_constraint.model = instance
    1323                                  piece_expression = self._create_piecewise_constraint_expression(this_lb, this_ub, x, xavg, quad_penalty_term_variable[index])
    1324                                  piece_constraint.add(None, piece_expression)
    1325                                  setattr(instance, piece_constraint_name, piece_constraint)                                   
    1326 
    1327                            else:
    1328 
    1329                               quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)                           
    1330                    
    1331          # strictly speaking, this probably isn't necessary - parameter coefficients won't get
    1332          # pre-processed out of the expression tree. however, if the under-the-hood should change,
    1333          # we'll be covered.
    1334          objective_expression.simplify(instance)
    1335          instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression
    1336          # if we are linearizing everything, then nothing will appear in the quadratic expression -
    1337          # don't add the empty "0.0" expression to the objective. otherwise, the output file won't
    1338          # be properly generated.
    1339          if quad_expression != 0.0:
    1340            instance.active_components(Objective)[objective_name]._quad_subexpr = quad_expression
     975
     976         new_attrs = form_ph_objective(instance_name, \
     977                                       instance, \
     978                                       self._original_objective_expression[instance_name], \
     979                                       self._scenario_tree, \
     980                                       self._linearize_nonbinary_penalty_terms, \
     981                                       self._drop_proximal_terms, \
     982                                       self._retain_quadratic_binary_terms, \
     983                                       self._breakpoint_strategy, \
     984                                       self._integer_tolerance)
     985         self._instance_augmented_attributes[instance_name].extend(new_attrs)
    1341986
    1342987   def iteration_k_solve(self):
     
    1350995
    1351996      solve_start_time = time.time()
     997
     998      # STEP -1: if using a PH solver manager, propagate current weights/averages to the appropriate solver servers.
     999      #          ditto the tree node statistics, which are necessary if linearizing (so an optimization could be
     1000      #          performed here).
     1001      # NOTE: We aren't currently propagating rhos, as they generally don't change - we need to
     1002      #       have a flag, though, indicating whether the rhos have changed, so they can be
     1003      #       transmitted if needed.
     1004      if self._solver_manager_type == "ph":
     1005         self._transmit_weights_and_averages()
     1006         self._transmit_tree_node_statistics()
    13521007
    13531008      # STEP 0: set up all global solver options.
     
    14431098            for objective_name in instance.active_components(Objective):
    14441099               objective = instance.active_components(Objective)[objective_name]
    1445                print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], 0.0)
     1100               print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], self._scenario_tree.compute_scenario_cost(instance))
    14461101
    14471102   def solve(self):
     
    14891144      for plugin in self._ph_plugins:     
    14901145         plugin.post_iteration_0(self)
     1146
     1147      # if using a PH solver server, trasnsmit the rhos prior to the iteration
     1148      # k solve sequence. for now, we are assuming that the rhos don't change
     1149      # on a per-iteration basis, but that assumption can be easily relaxed.
     1150      # it is important to do this after the plugins have a chance to do their
     1151      # computation.
     1152      if self._solver_manager_type == "ph":
     1153         self._transmit_rhos()
     1154         self._enable_ph_objectives()
    14911155
    14921156      # checkpoint if it's time - which it always is after iteration 0,
  • coopr.pysp/stable/2.3/coopr/pysp/phinit.py

    r2361 r2434  
    531531      print "Writing EF for remainder problem"
    532532      print ""
    533       write_ef(ph._scenario_tree, ph._instances, os.path.expanduser(options.ef_output_file))
     533      create_and_write_ef(ph._scenario_tree, ph._instances, os.path.expanduser(options.ef_output_file))
    534534
    535535   #
     
    643643   
    644644    return ans
    645 
  • coopr.pysp/stable/2.3/coopr/pysp/phserver.py

    r2361 r2434  
    66#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
    77#  the U.S. Government retains certain rights in this software.
    8 #  For more information, see the Coopr README.txt file.
     8#  Fr more information, see the Coopr README.txt file.
    99#  _________________________________________________________________________
    1010
     
    2323from coopr.opt.base import SolverFactory
    2424from coopr.pysp.scenariotree import *
     25from phutils import *
     26from phobjective import *
    2527
    2628# Pyro
    2729import Pyro.core
    2830import Pyro.naming
     31from Pyro.errors import PyroError,NamingError
    2932
    3033# garbage collection control.
     
    3538import pstats
    3639
     40import pickle
     41
     42# disable multi-threading. for a solver server, this
     43# would truly be a bad idea, as you (1) likely have limited
     44# solver licenses and (2) likely want to dedicate all compute
     45# resources on this node to a single solve.
     46Pyro.config.PYRO_MULTITHREADED=0
    3747
    3848#
     
    4151class PHSolverServer(object):
    4252
    43    def __init__(self, scenario_instances, solver):
    44 
    45       self._scenario_instances = scenario_instances
     53   def __init__(self, scenario_instances, solver, scenario_tree):
     54
     55      # the obvious stuff!
     56      self._instances = scenario_instances # a map from scenario name to pyomo model instance
    4657      self._solver = solver
     58      self._scenario_tree = scenario_tree
     59
     60      # the objective functions are modified throughout the course of PH iterations.
     61      # save the original, as a baseline to modify in subsequent iterations. reserve
     62      # the original objectives, for subsequent modification.
     63      self._original_objective_expression = {}
     64      for instance_name, instance in self._instances.items():
     65         objective_name = instance.active_components(Objective).keys()[0]
     66         expr = instance.active_components(Objective)[objective_name]._data[None].expr
     67         if isinstance(expr, Expression) is False:
     68            expr = _IdentityExpression(expr)
     69         self._original_objective_expression[instance_name] = expr
     70
     71      # the server is in one of two modes - solve the baseline instances, or
     72      # those augmented with the PH penalty terms. default is standard.
     73      # NOTE: This is currently a global flag for all scenarios handled
     74      #       by this server - easy enough to extend if we want.
     75      self._solving_standard_objective = True
     76
     77   def enable_standard_objective(self):
     78
     79#      print "RECEIVED REQUEST TO ENABLE STANDARD OBJECTIVE FOR ALL SCENARIOS"
     80
     81      self._solving_standard_objective = True
     82
     83   def enable_ph_objective(self):
     84
     85#      print "RECEIVED REQUEST TO ENABLE PH OBJECTIVE FOR ALL SCENARIOS"
     86
     87      self._solving_standard_objective = False
    4788
    4889   def solve(self, scenario_name):
    4990
    50       print "RECEIVED SOLVE REQUEST!"
    51       print "SCENARIO TO SOLVE=",scenario_name
    52       if scenario_name not in self._scenario_instances:
     91#      print "RECEIVED REQUEST TO SOLVE SCENARIO INSTANCE="+scenario_name
     92#      if self._solving_standard_objective is True:
     93#         print "OBJECTIVE=STANDARD"
     94#      else:
     95#         print "OBJECTIVE=PH"
     96         
     97      if scenario_name not in self._instances:
    5398         print "***ERROR: Requested instance to solve not in PH solver server instance collection!"
    5499         return None
    55       scenario_instance = self._scenario_instances[scenario_name]
    56       print "SCENARIO INSTANCE=",scenario_instance
    57       print "SUPPORTS MULTITHREADING=",Pyro.util.supports_multithreading()
    58       print "SOLVING"
    59       self._solver.solve(scenario_instance)
    60       print "DONE WITH SOLVE"
    61       return scenario_instance
    62 #      return result
    63 
     100      scenario_instance = self._instances[scenario_name]
     101
     102      # form the desired objective, depending on the solve mode.
     103      if self._solving_standard_objective is True:
     104         form_standard_objective(scenario_name, scenario_instance, self._original_objective_expression[scenario_name], self._scenario_tree)         
     105      else:
     106         # TBD - command-line drive the various options (integer tolerance, breakpoint strategies, etc.)
     107         form_ph_objective(scenario_name, scenario_instance, \
     108                           self._original_objective_expression[scenario_name], self._scenario_tree, \
     109                           False, False, False, 0, 0.00001)
     110
     111      # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
     112      # don't do this, you won't see any chance in the output files as you vary the problem parameters!
     113      # ditto for instance fixing!
     114      scenario_instance.preprocess()
     115
     116      results = self._solver.solve(scenario_instance)
     117
     118      print "Successfully solved scenario instance="+scenario_name
     119#      print "RESULTS:"
     120#      results.write(num=1)
     121      encoded_results = pickle.dumps(results)
     122     
     123      return encoded_results
     124
     125   def update_weights_and_averages(self, scenario_name, new_weights, new_averages):
     126
     127#      print "RECEIVED REQUEST TO UPDATE WEIGHTS AND AVERAGES FOR SCENARIO=",scenario_name
     128
     129      if scenario_name not in self._instances:
     130         print "***ERROR: Received request to update weights for instance not in PH solver server instance collection!"
     131         return None
     132      scenario_instance = self._instances[scenario_name]
     133     
     134      for weight_update in new_weights:
     135         
     136         weight_index = weight_update._index
     137
     138         target_weight_parameter = getattr(scenario_instance, weight_update.name)
     139
     140         for index in weight_index:
     141            target_weight_parameter[index] = weight_update[index]()
     142
     143      for average_update in new_averages:
     144         
     145         average_index = average_update._index
     146
     147         target_average_parameter = getattr(scenario_instance, average_update.name)
     148
     149         for index in average_index:
     150            target_average_parameter[index] = average_update[index]()
     151
     152   def update_rhos(self, scenario_name, new_rhos):
     153
     154#      print "RECEIVED REQUEST TO UPDATE RHOS FOR SCENARIO=",scenario_name
     155
     156      if scenario_name not in self._instances:
     157         print "***ERROR: Received request to update weights for instance not in PH solver server instance collection!"
     158         return None
     159      scenario_instance = self._instances[scenario_name]     
     160
     161      for rho_update in new_rhos:
     162
     163         rho_index = rho_update._index
     164
     165         target_rho_parameter = getattr(scenario_instance, rho_update.name)
     166
     167         for index in rho_index:
     168            target_rho_parameter[index] = rho_update[index]() # the value operator is crucial!
     169
     170   def update_tree_node_statistics(self, scenario_name, new_node_minimums, new_node_maximums):
     171
     172#      print "RECEIVED REQUEST TO UPDATE TREE NODE STATISTICS SCENARIO=",scenario_name
     173
     174      for tree_node_name, tree_node_minimums in new_node_minimums.items():
     175
     176         tree_node = self._scenario_tree._tree_node_map[tree_node_name]
     177         tree_node._minimums = tree_node_minimums
     178
     179      for tree_node_name, tree_node_maximums in new_node_maximums.items():
     180
     181         tree_node = self._scenario_tree._tree_node_map[tree_node_name]
     182         tree_node._maximums = tree_node_maximums         
     183     
     184           
    64185#
    65186# utility method to construct an option parser for ph arguments, to be
     
    80201                     dest="scenarios",
    81202                     default=[])
     203   parser.add_option("--all-scenarios",
     204                     help="Indicate that the server is responsible for solving all scenarios",
     205                     action="store_true",
     206                     dest="all_scenarios",
     207                     default=False)
    82208   parser.add_option("--report-solutions",
    83209                     help="Always report PH solutions after each iteration. Enabled if --verbose is enabled. Default is False.",
     
    151277                     type="int",
    152278                     default=0)
    153    parser.add_option("--enable-gc",
    154                      help="Enable the python garbage collecter. Default is True.",
    155                      action="store_true",
    156                      dest="enable_gc",
    157                      default=True)
     279   parser.add_option("--disable-gc",
     280                     help="Disable the python garbage collecter. Default is False.",
     281                     action="store_true",
     282                     dest="disable_gc",
     283                     default=False)
    158284
    159285   parser.usage=usage_string
     
    182308         print ""
    183309         print "***ERROR: Exiting test driver: No 'model' object created in module "+reference_model_filename
    184          return None, None, None
     310         return None, None, None, None
    185311
    186312      if model_import.model is None:
    187313         print ""
    188314         print "***ERROR: Exiting test driver: 'model' object equals 'None' in module "+reference_model_filename
    189          return None, None, None
     315         return None, None, None, None
    190316 
    191317      reference_model = model_import.model
    192318   except IOError:
    193319      print "***ERROR: Failed to load scenario reference model from file="+reference_model_filename
    194       return None, None, None
     320      return None, None, None, None
    195321
    196322   try:
     
    201327   except IOError:
    202328      print "***ERROR: Failed to load scenario reference instance data from file="+reference_instance_filename
    203       return None, None, None
     329      return None, None, None, None
    204330
    205331   #
     
    231357#
    232358
    233 def run_server(options, scenario_instances, solver):
     359def run_server(options, scenario_instances, solver, scenario_tree):
    234360
    235361   start_time = time.time()
    236362
    237    print "RUNNING SERVER NOW!"
    238    Pyro.core.initServer(banner=1)   
    239 
    240    print "NUMBER OF SCENARIO INSTANCES="+str(len(scenario_instances))
    241 
    242    print "ATTEMPTING TO LOCATE NAME SERVER"
     363   Pyro.core.initServer(banner=0)   
     364
    243365   locator = Pyro.naming.NameServerLocator()
    244    ns = locator.getNS()
    245    print "FOUND NAME SERVER"
     366   ns = None
     367   try:
     368      ns = locator.getNS()
     369   except Pyro.errors.NamingError:
     370      raise RuntimeError, "PH solver server failed to locate Pyro name server"     
    246371
    247372   solver_daemon = Pyro.core.Daemon()
    248373   solver_daemon.useNameServer(ns)
    249374
    250    daemon_object = PHSolverServer(scenario_instances, solver)
     375   daemon_object = PHSolverServer(scenario_instances, solver, scenario_tree)
    251376   delegator_object = Pyro.core.ObjBase()
    252377   delegator_object.delegateTo(daemon_object)
     
    255380   # each scenario processed by this daemon.
    256381   for scenario_name, scenario_instance in scenario_instances.items():
    257       solver_daemon.connect(delegator_object, scenario_name)
    258 
    259    while True:
    260       solver_daemon.handleRequests()
    261       print "HANDLED A REQUEST!"
    262 
    263    solver_daemon.disconnect()
    264 
    265    print "DONE RUNNING SERVER!"
     382      # first, see if it's already registered - if no, cull the old entry.
     383      # NOTE: This is a hack, as the cleanup / disconnect code doesn't
     384      #       seem to be invoked when the daemon is killed.
     385      try:
     386         ns.unregister(scenario_name)
     387      except NamingError:
     388         pass
     389
     390      try:
     391         solver_daemon.connect(delegator_object, scenario_name)
     392      except Pyro.errors.NamingError:
     393         raise RuntimeError, "Entry in name server already exists for scenario="+scenario_name
     394
     395   try:
     396      solver_daemon.requestLoop()
     397   except KeyboardInterrupt, SystemExit:
     398      pass
     399
     400   solver_daemon.disconnect(delegator_object)
     401   solver_daemon.shutdown()
    266402
    267403   end_time = time.time()
     
    273409def exec_server(options):
    274410
    275    # we need the base model (not so much the reference instance, but that -
    276    # can't hurt too much until proven otherwise, that is) to construct
     411   # we need the base model (not so much the reference instance, but that
     412   # can't hurt too much - until proven otherwise, that is) to construct
    277413   # the scenarios that this server is responsible for.
    278    print "LOADING REFERENCE AND SCENARIO MODELS/INSTANCES"
    279414   reference_model, reference_instance, scenario_tree, scenario_tree_instance = load_reference_and_scenario_models(options)
    280415   if reference_model is None or reference_instance is None or scenario_tree is None:
     
    284419   # construct the set of scenario instances based on the command-line.
    285420   scenario_instances = {}
    286    if len(options.scenarios) == 0:
    287       print "***Unable to launch PH solver server - no scenario specified!"
     421   if (options.all_scenarios is True) and (len(options.scenarios) > 0):
     422      print "***WARNING: Both individual scenarios and all-scenarios were specified on the command-line; proceeding using all scenarios"
     423
     424   if (len(options.scenarios) == 0) and (options.all_scenarios is False):
     425      print "***Unable to launch PH solver server - no scenario(s) specified!"
    288426      return
    289427
    290    for scenario_name in options.scenarios:
     428   scenarios_to_construct = []
     429   if options.all_scenarios is True:
     430      scenarios_to_construct.extend(scenario_tree._scenario_map.keys())
     431   else:
     432      scenarios_to_construct.extend(options.scenarios)
     433
     434   for scenario_name in scenarios_to_construct:
    291435      print "Creating instance for scenario="+scenario_name
    292436      if scenario_tree.contains_scenario(scenario_name) is False:
     
    294438         return
    295439
     440      # create the baseline scenario instance
    296441      scenario_instance = construct_scenario_instance(scenario_tree,
    297442                                                      options.instance_directory,
     
    301446      if scenario_instance is None:
    302447         print "***Unable to launch PH solver server - failed to create instance for scenario="+scenario_name
     448         return
     449
     450      # augment the instance with PH-specific parameters (weights, rhos, etc).
     451      # TBD: The default rho of 1.0 is kind of bogus. Need to somehow propagate
     452      #      this value and the linearization parameter as a command-line argument.
     453      new_penalty_variable_names = create_ph_parameters(scenario_instance, scenario_tree, 1.0, False)
     454
    303455      scenario_instances[scenario_name] = scenario_instance
    304456
     
    318470      solver._report_timing = True
    319471
    320    print "SOLVING!"
    321    results = solver.solve(scenario_instance)
    322    print "TYPE OF RESULTS=",type(results)
    323    results.write(num=1)
    324 #   print "RESULTS=",results
    325    print "DONE SOLVING!"
    326 
    327472   # spawn the daemon.
    328    run_server(options, scenario_instances, solver)
     473   run_server(options, scenario_instances, solver, scenario_tree)
    329474
    330475def run(args=None):
     
    349494    # much sense - so it is disabled by default. Because: It drops
    350495    # the run-time by a factor of 3-4 on bigger instances.
    351     if options.enable_gc is False:
     496    if options.disable_gc is True:
    352497       gc.disable()
    353498    else:
  • coopr.pysp/stable/2.3/coopr/pysp/phutils.py

    r2391 r2434  
    224224
    225225   return scenario_instance
     226
     227
     228# creates all PH parameters for a problem instance, given a scenario tree
     229# (to identify the relevant variables), a default rho (simply for initialization),
     230# and a boolean indicating if quadratic penalty terms are to be linearized.
     231# returns a list of any created variables, specifically when linearizing -
     232# this is useful to clean-up reporting.
     233
     234def create_ph_parameters(instance, scenario_tree, default_rho, linearizing_penalty_terms):
     235
     236   new_penalty_variable_names = []   
     237   
     238   # first, gather all unique variables referenced in any stage
     239   # other than the last, independent of specific indices. this
     240   # "gather" step is currently required because we're being lazy
     241   # in terms of index management in the scenario tree - which
     242   # should really be done in terms of sets of indices.
     243   # NOTE: technically, the "instance variables" aren't really references
     244   # to the variable in the instance - instead, the reference model. this
     245   # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
     246   instance_variables = {}
     247   
     248   for stage in scenario_tree._stages[:-1]:
     249     
     250      for (reference_variable, index_template, reference_indices) in stage._variables:
     251           
     252         if reference_variable.name not in instance_variables.keys():
     253               
     254            instance_variables[reference_variable.name] = reference_variable
     255
     256   # for each blended variable, create a corresponding ph weight and average parameter in the instance.
     257   # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor
     258   # in the grand scheme of things.
     259
     260   for (variable_name, reference_variable) in instance_variables.items():
     261
     262      # PH WEIGHT
     263
     264      new_w_index = reference_variable._index
     265      new_w_parameter_name = "PHWEIGHT_"+reference_variable.name
     266      # this bit of ugliness is due to Pyomo not correctly handling the Param construction
     267      # case when the supplied index set consists strictly of None, i.e., the source variable
     268      # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed.
     269      new_w_parameter = None
     270      if (len(new_w_index) is 1) and (None in new_w_index):
     271         new_w_parameter = Param(name=new_w_parameter_name)
     272      else:
     273         new_w_parameter = Param(new_w_index,name=new_w_parameter_name)
     274      setattr(instance,new_w_parameter_name,new_w_parameter)
     275
     276      # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference
     277      # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not
     278      # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the
     279      # issue in any case for now. NOTE: I'm not sure this is still the case in Pyomo.
     280      for index in new_w_index:
     281         new_w_parameter[index] = 0.0
     282
     283      # PH AVG
     284
     285      new_avg_index = reference_variable._index
     286      new_avg_parameter_name = "PHAVG_"+reference_variable.name
     287      new_avg_parameter = None
     288      if (len(new_avg_index) is 1) and (None in new_avg_index):
     289         new_avg_parameter = Param(name=new_avg_parameter_name)
     290      else:
     291         new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
     292      setattr(instance,new_avg_parameter_name,new_avg_parameter)
     293
     294      for index in new_avg_index:
     295         new_avg_parameter[index] = 0.0
     296
     297      # PH RHO
     298
     299      new_rho_index = reference_variable._index
     300      new_rho_parameter_name = "PHRHO_"+reference_variable.name
     301      new_rho_parameter = None
     302      if (len(new_avg_index) is 1) and (None in new_avg_index):
     303         new_rho_parameter = Param(name=new_rho_parameter_name)
     304      else:
     305         new_rho_parameter = Param(new_rho_index,name=new_rho_parameter_name)
     306      setattr(instance,new_rho_parameter_name,new_rho_parameter)
     307
     308      for index in new_rho_index:
     309         new_rho_parameter[index] = default_rho
     310
     311      # PH LINEARIZED PENALTY TERM
     312
     313      if linearizing_penalty_terms > 0:
     314
     315         new_penalty_term_variable_index = reference_variable._index
     316         new_penalty_term_variable_name = "PHQUADPENALTY_"+reference_variable.name
     317         # this is a quadratic penalty term - the lower bound is 0!
     318         new_penalty_term_variable = None
     319         if (len(new_avg_index) is 1) and (None in new_avg_index):
     320            new_penalty_term_variable = Var(name=new_penalty_term_variable_name, bounds=(0.0,None))                 
     321         else:
     322            new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None))
     323         new_penalty_term_variable.construct()
     324         setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable)
     325         new_penalty_variable_names.append(new_penalty_term_variable_name)
     326
     327      # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY.
     328
     329      # also controls whether weight updates proceed at any iteration.
     330
     331      new_blend_index = reference_variable._index
     332      new_blend_parameter_name = "PHBLEND_"+reference_variable.name
     333      new_blend_parameter = None
     334      if (len(new_avg_index) is 1) and (None in new_avg_index):
     335         new_blend_parameter = Param(name=new_blend_parameter_name, within=Binary)
     336      else:
     337         new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary)
     338      setattr(instance,new_blend_parameter_name,new_blend_parameter)
     339
     340      for index in new_blend_index:
     341         new_blend_parameter[index] = 1
     342
     343   return new_penalty_variable_names     
  • coopr.pysp/stable/2.3/coopr/pysp/scenariotree.py

    r2391 r2434  
    1111import sys
    1212import types
    13 from coopr.pyomo import *
    1413import copy
    1514import os.path
    1615import traceback
    1716
     17from coopr.pyomo import *
    1818from phutils import *
    1919
     
    2323
    2424   """
    25    def __init__(self, *args, **kwds):
    26 
    27       self._name = ""
    28       self._stage = None
     25   def __init__(self, name, conditional_probability, stage, reference_instance):
     26
     27      self._name = name
     28      self._stage = stage
    2929      self._parent = None
    3030      self._children = [] # a collection of ScenarioTreeNodes
    31       self._conditional_probability = None # conditional on parent
     31      self._conditional_probability = conditional_probability # conditional on parent
    3232      self._scenarios = [] # a collection of all Scenarios passing through this node in the tree
    3333
     
    4444      self._maximums = {}
    4545
     46      # solution (variable) values for this node. assumed to be distinct
     47      # from self._averages, as the latter are not necessarily feasible.
     48      # objects in the map are actual variables.
     49      self._solution = {}
     50
     51      # for each variable referenced in the stage, clone the variable
     52      # for purposes of storing solutions. we are being wasteful in
     53      # terms copying indices that may not be referenced in the stage.
     54      # this is something that we might revisit if space/performance
     55      # is an issue (space is the most likely issue)
     56      for variable, match_template, variable_indices in self._stage._variables:
     57
     58         # don't bother copying bounds for variables, as the values stored
     59         # here are computed elsewhere - and that entity is responsible for
     60         # ensuring feasibility. this also leaves room for specifying infeasible
     61         # or partial solutions.       
     62
     63         new_variable_index = variable._index
     64         new_variable_name = variable._name
     65         new_variable = None
     66         if (len(new_variable_index) is 1) and (None in new_variable_index):
     67            new_variable = Var(name=new_variable_name)
     68         else:
     69            new_variable = Var(new_variable_index, name=new_variable_name)
     70
     71         self._solution[new_variable_name] = new_variable
     72
    4673   #
    4774   # a utility to compute the cost of the current node plus the expected costs of child nodes.
    4875   #
     76
    4977   def computeExpectedNodeCost(self, scenario_instance_map):
    5078
     
    90118   """ Constructor
    91119       Arguments:
    92            scenarioinstance             - the reference scenario instance.
     120           scenarioinstance     - the reference (deterministic) scenario instance.
    93121           scenariotreeinstance - the pyomo model specifying all scenario tree (text) data.
    94122           scenariobundlelist   - a list of scenario names to retain, i.e., cull the rest to create a reduced tree!
     
    172200      # can't do a single pass because the objects may not exist.
    173201      for tree_node_name in node_ids:
    174          new_tree_node = ScenarioTreeNode()
    175          new_tree_node._name = tree_node_name
     202
     203         if tree_node_name not in node_stage_ids:
     204            raise ValueError, "No stage is assigned to tree node=" + tree_node._name
     205
     206         stage_name = node_stage_ids[tree_node_name].value
     207         if stage_name not in self._stage_map.keys():
     208            raise ValueError, "Unknown stage=" + stage_name + " assigned to tree node=" + tree_node._name
     209
     210         new_tree_node = ScenarioTreeNode(tree_node_name,
     211                                          node_probability_map[tree_node_name].value,
     212                                          self._stage_map[stage_name],
     213                                          self._reference_instance)
    176214
    177215         self._tree_nodes.append(new_tree_node)
    178216         self._tree_node_map[tree_node_name] = new_tree_node
    179 
    180          new_tree_node._conditional_probability = node_probability_map[tree_node_name].value
    181 
    182          if tree_node_name not in node_stage_ids:
    183             raise ValueError, "No stage is assigned to tree node=" + tree_node._name
    184          else:
    185             stage_name = node_stage_ids[new_tree_node._name].value
    186             if stage_name not in self._stage_map.keys():
    187                raise ValueError, "Unknown stage=" + stage_name + " assigned to tree node=" + tree_node._name
    188             else:
    189                new_tree_node._stage = self._stage_map[stage_name]
    190                self._stage_map[stage_name]._tree_nodes.append(new_tree_node)
     217         self._stage_map[stage_name]._tree_nodes.append(new_tree_node)
    191218
    192219      # link up the tree nodes objects based on the child id sets.
     
    347374   # is the indicated scenario in the tree?
    348375   #
     376
    349377   def contains_scenario(self, name):
    350378      return name in self._scenario_map.keys()
     
    353381   # get the scenario object from the tree.
    354382   #
     383
    355384   def get_scenario(self, name):
    356385      return self._scenario_map[name]
     386
     387   #
     388   # compute the scenario cost for the input instance, i.e.,
     389   # the sum of all stage cost variables.
     390   #
     391
     392   def compute_scenario_cost(self, instance):
     393      aggregate_cost = 0.0
     394      for stage in self._stages:
     395         instance_cost_variable = instance.active_components(Var)[stage._cost_variable[0].name][stage._cost_variable[1]]()
     396         aggregate_cost += instance_cost_variable
     397      return aggregate_cost
    357398
    358399   #
     
    363404   # of the scenario tree structure.
    364405   #
     406
    365407   def compress(self, scenario_bundle_list):
    366408
     
    450492   # returns the root node of the scenario tree
    451493   #
     494
    452495   def findRootNode(self):
    453496
     
    461504   # the maximal length of identifiers in various categories.
    462505   #
     506
    463507   def computeIdentifierMaxLengths(self):
    464508
     
    471515   # a utility function to (partially, at the moment) validate a scenario tree
    472516   #
     517
    473518   def validate(self):
    474519
  • coopr.pysp/stable/2.3/coopr/pysp/tests/unit/test_ph.py

    r2391 r2434  
    1818# Import the testing packages
    1919#
    20 import unittest
    2120import pyutilib.misc
    22 import pyutilib.th
     21import pyutilib.th as unittest
    2322import pyutilib.subprocess
    2423import coopr.pysp
     
    4645
    4746#
    48 # Define a testing class, using the pyutilib.th.TestCase class.  This is
    49 # an extension of the unittest.TestCase class that adds additional testing
    50 # functions.
     47# Define a testing class, using the unittest.TestCase class.
    5148#
    52 class TestPH(pyutilib.th.TestCase):
     49class TestPH(unittest.TestCase):
    5350
    5451    def cleanup(self):
     
    6158           del sys.modules["ReferenceModel"]
    6259
     60    @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available")
    6361    def test_quadratic_farmer(self):
    64 
    65         if cplex_available is False:
    66            return
    67 
    6862        pyutilib.misc.setup_redirect(current_directory+"farmer_quadratic.out")
    6963        farmer_examples_dir = pysp_examples_dir + "farmer"
     
    7771        self.failUnlessFileEqualsBaseline(current_directory+"farmer_quadratic.out",current_directory+"farmer_quadratic.baseline", filter=filter_time)
    7872       
     73    @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available")
    7974    def test_linearized_farmer(self):
    80 
    81         solver_string = ""
    8275        if cplex_available:
    8376            solver_string="cplex"
    84         #elif glpk_available:
    85             #solver_string="glpk"
    86         else:
    87            return
    88 
    8977        pyutilib.misc.setup_redirect(current_directory+"farmer_linearized.out")
    9078        farmer_examples_dir = pysp_examples_dir + "farmer"
     
    9886        self.failUnlessFileEqualsBaseline(current_directory+"farmer_linearized.out",current_directory+"farmer_linearized.baseline", filter=filter_time)
    9987
     88    @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available")
    10089    def test_linearized_farmer_nodedata(self):
    101 
    102         solver_string = ""
    10390        if cplex_available:
    10491            solver_string="cplex"
    105         #elif glpk_available:
    106             #solver_string="glpk"
    107         else:
    108            return
    109 
    11092        pyutilib.misc.setup_redirect(current_directory+"farmer_linearized_nodedata.out")
    11193        farmer_examples_dir = pysp_examples_dir + "farmer"
     
    119101        self.failUnlessFileEqualsBaseline(current_directory+"farmer_linearized_nodedata.out",current_directory+"farmer_linearized_nodedata.baseline", filter=filter_time)       
    120102
     103    @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available")
    121104    def test_quadratic_sizes3(self):
    122105
    123106        # Ignore this for now
    124107        return
    125 
    126         if cplex_available is False:
    127            return
    128108       
    129109        pyutilib.misc.setup_redirect(current_directory+"sizes3_quadratic.out")
     
    146126
    147127if __name__ == "__main__":
    148    unittest.main()
     128    unittest.main()
  • coopr.pysp/stable/2.3/scripts/ph_test_client

    r2361 r2434  
    44import Pyro.naming
    55
    6 from coopr.opt.base import SolverFactory
    7 from coopr.opt import *
     6import pickle
     7
     8import pyutilib.misc
     9import pyutilib.component.core
     10from coopr.opt.results import SolverResults
    811
    912import sys
     
    3841
    3942print "CALLING SOLVE METHOD"
    40 result = obj.solve(scenario_name)
     43encoded_result = obj.solve(scenario_name)
     44result = pickle.loads(encoded_result)
     45print "SUCCESSFULLY OBTAINED RESULT!"
    4146
    42 print "RESULT=",result
     47print "RESULTS:"
     48result.write(num=1)
    4349
    4450print "DONE"
Note: See TracChangeset for help on using the changeset viewer.