source: pyomo/trunk/pyomo/pysp/ph.py @ 9519

Last change on this file since 9519 was 9519, checked in by wehart, 4 years ago

Adding copyright headers for files that are missing them.

Updating some of the copyright assertions

File size: 191.1 KB
Line 
1#  _________________________________________________________________________
2#
3#  Pyomo: Python Optimization Modeling Objects
4#  Copyright (c) 2014 Sandia Corporation.
5#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
6#  the U.S. Government retains certain rights in this software.
7#  This software is distributed under the BSD License.
8#  _________________________________________________________________________
9
10_OLD_OUTPUT = True
11
12import copy
13import gc
14import logging
15import pickle
16import sys
17import traceback
18import types
19import time
20import types
21import inspect
22from operator import itemgetter
23import uuid
24
25from math import fabs, log, exp, sqrt
26from os import path
27
28from six import iterkeys, itervalues, iteritems, advance_iterator
29
30from pyomo.opt import SolverResults, SolverStatus, UndefinedData, ProblemFormat, undefined
31from pyomo.opt.base import SolverFactory
32from pyomo.opt.parallel import SolverManagerFactory
33from pyomo.core import *
34from pyomo.core.base import BasicSymbolMap, CounterLabeler
35from pyomo.pysp.phextension import IPHExtension
36from pyomo.pysp.ef import create_ef_instance
37from pyomo.pysp.generators import scenario_tree_node_variables_generator, \
38                                    scenario_tree_node_variables_generator_noinstances
39from pyomo.pysp.phsolverserverutils import *
40from pyomo.pysp.phsolverserverutils import TransmitType
41from pyomo.pysp.convergence import *
42from pyomo.pysp.phutils import *
43from pyomo.pysp.phobjective import *
44from pyomo.pysp.scenariotree import *
45from pyomo.pysp.dualphmodel import DualPHModel
46import pyomo.solvers.plugins.smanager.phpyro
47 
48from pyomo.util.plugin import ExtensionPoint
49
50import pyutilib.common
51
52try:
53    from guppy import hpy
54    guppy_available = True
55except ImportError:
56    guppy_available = False
57
58logger = logging.getLogger('pyomo.pysp')
59
60# PH iteratively solves scenario sub-problems, so we don't want to
61# waste a ton of time preprocessing unless some specific aspects of
62# the scenario instances change.  for example, a variable was fixed,
63# the objective was modified, or constraints were added. and if
64# instances do change, we only want to do the minimal amount of work
65# to get the instance back to a consistent "preprocessed" state.  the
66# following attributes are introduced to help perform the minimal
67# amount of work, and should be augmented in the future if we can
68# somehow do less.  these attributes are initially cleared, and are
69# re-set - following preprocessing, if necessary - at the top of the
70# PH iteration loop. this gives a chance for plugins and linearization
71# to get a chance at modification, and to set the appropriate
72# attributes so that the instances can be appropriately preprocessed
73# before solves for the next iteration commence. we assume (by
74# prefixing the attribute name with "instance") that modifications of
75# the indicated type have been uniformly applied to all instances.
76class ProblemStates(object):
77
78    #def __getstate__(self):
79    #    return dict(self.__slots__)
80
81    def __init__(self, instances):
82
83        # ph objects added to each model
84        self.has_ph_objective_weight_terms = dict.fromkeys(instances, False)
85        self.has_ph_objective_proximal_terms = dict.fromkeys(instances, False)
86        self.ph_objective_proximal_expressions = dict.fromkeys(instances, None)
87        self.ph_objective_weight_expressions = dict.fromkeys(instances, None)
88        self.ph_constraints = dict((inst_name,[]) for inst_name in instances)
89        self.ph_variables = dict((inst_name,[]) for inst_name in instances)
90
91        # TODO: Reconcile this new method with the persistent solver plugin
92        """
93        # keeps track of instances with recently fixed or freed variables
94        self.fixed_variables = dict.fromkeys(instances, False)
95        self.freed_variables = dict.fromkeys(instances, False)
96        """
97        # maps between instance name and a list of (variable-name, index) pairs
98        self.fixed_variables = dict((inst_name,[]) for inst_name in instances)
99        self.freed_variables = dict((inst_name,[]) for inst_name in instances)
100
101        # just coefficients modified
102        self.objective_updated = dict.fromkeys(instances, False)
103        self.ph_constraints_updated = dict.fromkeys(instances, False)
104        self.user_constraints_updated = dict.fromkeys(instances, False)
105
106    def clear_update_flags(self,name=None):
107        if name is not None:
108            self.objective_updated[name] = False
109            self.ph_constraints_updated[name] = False
110            self.user_constraints_updated[name] = False
111        else:
112            for key in iterkeys(self.objective_updated):
113                self.objective_updated[key] = False
114            for key in iterkeys(self.ph_constraints_updated):
115                self.ph_constraints_updated[key] = False
116            for key in iterkeys(self.user_constraints_updated):
117                self.user_constraints_updated[key] = False
118
119    # TODO
120    """
121    def has_fixed_variables(self,name=None):
122        if name is None:
123            for val in itervalues(self.fixed_variables):
124                if val:
125                    return True
126            return False
127        else:
128            return self.fixed_variables[name]
129
130    def has_freed_variables(self,name=None):
131        if name is None:
132            for val in itervalues(self.freed_variables):
133                if val:
134                    return True
135            return False
136        else:
137            return self.freed_variables[name]
138    """
139    def has_fixed_variables(self,name=None):
140        if name is None:
141            for val in itervalues(self.fixed_variables):
142                if len(val) > 0:
143                    return True
144            return False
145        else:
146            return len(self.fixed_variables[name]) > 0
147
148    def has_freed_variables(self,name=None):
149        if name is None:
150            for val in itervalues(self.freed_variables):
151                if len(val) > 0:
152                    return True
153            return False
154        else:
155            return len(self.freed_variables[name]) > 0
156
157    def has_ph_constraints(self,name=None):
158        if name is None:
159            for val in itervalues(self.ph_constraints):
160                if len(val) > 0:
161                    return True
162            return False
163        else:
164            return len(self.ph_constraints[name]) > 0
165
166    def has_ph_variables(self,name=None):
167        if name is None:
168            for val in itervalues(self.ph_variables):
169                if len(val) > 0:
170                    return True
171            return False
172        else:
173            return len(self.ph_variables[name]) > 0
174
175    # TODO
176    """
177    def clear_fixed_variables(self, name=None):
178        if name is None:
179            for key in self.fixed_variables:
180                self.fixed_variables[key] = False
181        else:
182            if name in self.fixed_variables:
183                self.fixed_variables[name] = False
184            else:
185                raise KeyError("KeyError: %s" % name)
186
187    def clear_freed_variables(self, name=None):
188        if name is None:
189            for key in self.freed_variables:
190                self.freed_variables[key] = False
191        else:
192            if name in self.freed_variables:
193                self.freed_variables[name] = False
194            else:
195                raise KeyError("KeyError: %s" % name)
196    """
197
198    def clear_fixed_variables(self, name=None):
199        if name is None:
200            for key in self.fixed_variables:
201                self.fixed_variables[key] = []
202        else:
203            if name in self.fixed_variables:
204                self.fixed_variables[name] = []
205            else:
206                raise KeyError("KeyError: %s" % name)
207
208    def clear_freed_variables(self, name=None):
209        if name is None:
210            for key in self.freed_variables:
211                self.freed_variables[key] = []
212        else:
213            if name in self.freed_variables:
214                self.freed_variables[name] = []
215            else:
216                raise KeyError("KeyError: %s" % name)
217
218    def clear_ph_variables(self, name=None):
219        if name is None:
220            for key in self.ph_variables:
221                self.ph_variables[key] = []
222        else:
223            if name in self.ph_variables:
224                self.ph_variables[name] = []
225            else:
226                raise KeyError("KeyError: %s" % name)
227
228    def clear_ph_constraints(self, name=None):
229        if name is None:
230            for key in self.ph_constraints:
231                self.ph_constraints[key] = []
232        else:
233            if name in self.ph_constraints:
234                self.ph_constraints[name] = []
235            else:
236                raise KeyError("KeyError: %s" % name)
237
238
239def assign_aggregate_data(ph,
240                          scenario_tree,
241                          scenario_tree_object,
242                          aggregate_data):
243    ph._aggregate_user_data = aggregate_data
244
245class _PHBase(object):
246
247    def __init__(self):
248
249        # PH solver information / objects.
250        self._solver = None
251        self._solver_type = "cplex"
252        self._solver_io = None
253        self._comparison_tolerance_for_fixed_vars = 1e-5
254
255        self._problem_states = None
256
257        # a flag indicating whether we preprocess constraints in our
258        # scenario instances when variables are fixed/freed, or
259        # whether we simply write the bounds while presenting the
260        # instances to solvers.
261        self._write_fixed_variables = True
262        self._fix_queue = {}
263
264        # For the users to modify as they please in the aggregate
265        # callback as long as the data placed on it can be serialized
266        # by Pyro
267        self._aggregate_user_data = {}
268
269        # maps scenario name to the corresponding model instance
270        self._instances = {}
271
272        # for various reasons (mainly hacks at this point), it's good
273        # to know whether we're minimizing or maximizing.
274        self._objective_sense = None
275        self._objective_sense_option = None
276
277        # maps scenario name (or bundle name, in the case of bundling)
278        # to the last gap reported by the solver when solving the
279        # associated instance. if there is no entry, then there has
280        # been no solve.
281        # NOTE: This dictionary could expand significantly, as we
282        #       identify additional solve-related information
283        #       associated with an instance.
284        self._gaps = {}
285
286        # maps scenario name (or bundle name, in the case of bundling)
287        # to the last solve time reported for the corresponding
288        # sub-problem.
289        # presently user time, due to deficiency in solver plugins. ultimately
290        # want wall clock time for PH reporting purposes.
291        self._solve_times = {}
292
293        # defines the stochastic program structure and links in that
294        # structure with the scenario instances (e.g., _VarData
295        # objects).
296        self._scenario_tree = None
297
298        # there are situations in which it is valuable to snapshot /
299        # store the solutions associated with the scenario
300        # instances. for example, when one wants to use a warm-start
301        # from a particular iteration solve, following a modification
302        # and re-solve of the problem instances in a user-defined
303        # callback. the following nested dictionary is intended to
304        # serve that purpose. The nesting is dependent on whether
305        # bundling and or phpyro is in use
306        self._cached_solutions = {}
307        self._cached_scenariotree_solutions = {}
308
309        # results objects from the most recent round of solves These
310        # may hold more information than just variable values and so
311        # can be useful to hold on to until the next round of solves
312        # (keys are bundle name or scenario name)
313        self._solver_results = {}
314
315        # PH reporting parameters
316
317        # do I flood the screen with status output?
318        self._verbose = False
319        self._output_times = False
320
321        # PH configuration parameters
322
323        # a default, global value for rho. 0 indicates unassigned.
324        self._rho = 0.0
325
326        # do I retain quadratic objective terms associated with binary
327        # variables? in general, there is no good reason to not
328        # linearize, but just in case, we introduced the option.
329        self._retain_quadratic_binary_terms = False
330
331        # do I linearize the quadratic penalty term for continuous
332        # variables via a piecewise linear approximation? the default
333        # should always be 0 (off), as the user should be aware when
334        # they are forcing an approximation.
335        self._linearize_nonbinary_penalty_terms = 0
336
337        # the breakpoint distribution strategy employed when
338        # linearizing. 0 implies uniform distribution between the
339        # variable lower and upper bounds.
340        self._breakpoint_strategy = 0
341
342        # PH default tolerances - for use in fixing and testing
343        # equality across scenarios, and other stuff.
344        self._integer_tolerance = 0.00001
345
346        # when bundling, we cache the extensive form binding instances
347        # to save re-generation costs.
348        # maps bundle name in a scenario tree to the binding instance
349        self._bundle_binding_instance_map = {}
350        # maps bundle name in a scenario tree to a name->instance map
351        # of the scenario instances in the bundle
352        self._bundle_scenario_instance_map = {}
353
354        # a simple boolean flag indicating whether or not this ph
355        # instance has received an initialization method and has
356        # successfully processed it.
357        self._initialized = False
358
359    def initialize(self, *args, **kwds):
360        raise NotImplementedError("_PHBase::initialize() is an abstract method")
361
362    #
363    # Creates a deterministic symbol map for variables on an
364    # instance. This allows convenient transmission of information to
365    # and from PHSolverServers and makes it easy to save solutions
366    # using a pickleable dictionary of symbols -> values
367    #
368    def _create_instance_symbol_maps(self, ctypes):
369
370        for instance in itervalues(self._instances):
371
372            create_block_symbol_maps(instance, ctypes)
373
374    def _setup_scenario_instances(self):
375
376        self._problem_states = \
377            ProblemStates([scen._name for scen in \
378                           self._scenario_tree._scenarios])
379
380        for scenario in self._scenario_tree._scenarios:
381
382            scenario_instance = scenario._instance
383
384            assert scenario_instance.name == scenario._name
385
386            if scenario_instance is None:
387                raise RuntimeError("ScenarioTree has not been linked "
388                                   "with Pyomo model instances")
389
390            self._problem_states.objective_updated[scenario._name] = True
391            self._problem_states.user_constraints_updated[scenario._name] = True
392
393            # IMPT: disable canonical representation construction
394            #       for ASL solvers.  this is a hack, in that we
395            #       need to address encodings and the like at a
396            #       more general level.
397            if self._solver.problem_format() == ProblemFormat.nl:
398                scenario_instance.skip_canonical_repn = True
399                # We will take care of these manually within
400                # _preprocess_scenario_instance This will also
401                # prevent regenerating the ampl_repn when forming
402                # the bundle_ef's
403                scenario_instance.gen_obj_ampl_repn = False
404                scenario_instance.gen_con_ampl_repn = False
405
406            self._instances[scenario._name] = scenario_instance
407
408    def solve(self, *args, **kwds):
409        raise NotImplementedError("_PHBase::solve() is an abstract method")
410
411    # restores the variable values for all of the scenario instances
412    # that I maintain.  restoration proceeds from the
413    # self._cached_solutions map. if this is not populated (via
414    # setting cache_results=True when calling solve_subproblems), then
415    # an exception will be thrown.
416
417    def restoreCachedSolutions(self, cache_id, release_cache):
418
419        cache = self._cached_scenariotree_solutions.get(cache_id,None)
420        if cache is None:
421            raise RuntimeError("PH scenario tree solution cache "
422                               "with id %s does not exist"
423                               % (cache_id))
424        if release_cache and (cache is not None):
425            del self._cached_scenariotree_solutions[cache_id]
426
427        for scenario in self._scenario_tree._scenarios:
428
429            scenario.update_current_solution(cache[scenario._name])
430
431        if (not len(self._bundle_binding_instance_map)) and \
432           (not len(self._instances)):
433            return
434
435        cache = self._cached_solutions.get(cache_id,None)
436        if cache is None:
437                raise RuntimeError("PH scenario tree solution cache "
438                                   "with id %s does not exist"
439                                   % (cache_id))
440
441        if release_cache and (cache is not None):
442            del self._cached_solutions[cache_id]
443
444        if self._scenario_tree.contains_bundles():
445
446            for bundle_name, bundle_ef_instance in iteritems(self._bundle_binding_instance_map):
447
448                results, fixed_results = cache[bundle_name]
449
450                for scenario_name, scenario_fixed_results in iteritems(fixed_results):
451                    scenario_instance = self._instances[scenario_name]
452                    bySymbol = scenario_instance._PHInstanceSymbolMaps[Var].bySymbol
453                    for instance_id, varvalue, stale_flag in scenario_fixed_results:
454                        vardata = bySymbol[instance_id]
455                        vardata.fix(varvalue)
456                        vardata.stale = stale_flag
457
458                bundle_ef_instance.load(
459                    results,
460                    allow_consistent_values_for_fixed_vars=self._write_fixed_variables,
461                    comparison_tolerance_for_fixed_vars=self._comparison_tolerance_for_fixed_vars)
462
463        else:
464            for scenario_name, scenario_instance in iteritems(self._instances):
465
466                results, fixed_results = cache[scenario_name]
467
468                bySymbol = scenario_instance._PHInstanceSymbolMaps[Var].bySymbol
469                for instance_id, varvalue, stale_flag in fixed_results:
470                    vardata = bySymbol[instance_id]
471                    vardata.fix(varvalue)
472                    vardata.stale = stale_flag
473
474                scenario_instance.load(
475                    results,
476                    allow_consistent_values_for_fixed_vars=self._write_fixed_variables,
477                    comparison_tolerance_for_fixed_vars=self._comparison_tolerance_for_fixed_vars)
478
479    def cacheSolutions(self, cache_id):
480
481        for scenario in self._scenario_tree._scenarios:
482            self._cached_scenariotree_solutions.\
483                setdefault(cache_id,{})[scenario._name] = \
484                    scenario.package_current_solution()
485
486        if self._scenario_tree.contains_bundles():
487
488            for bundle_name, scenario_map in iteritems(self._bundle_scenario_instance_map):
489
490                fixed_results = {}
491                for scenario_name, scenario_instance in iteritems(scenario_map):
492
493                    fixed_results[scenario_name] = \
494                        tuple((instance_id, vardata.value, vardata.stale) \
495                              for instance_id, vardata in \
496                              iteritems(scenario_instance.\
497                                        _PHInstanceSymbolMaps[Var].\
498                                        bySymbol) \
499                              if vardata.fixed)
500
501                self._cached_solutions.\
502                    setdefault(cache_id,{})[bundle_name] = \
503                        (self._solver_results[bundle_name],
504                         fixed_results)
505
506        else:
507
508            for scenario_name, scenario_instance in iteritems(self._instances):
509
510                fixed_results = \
511                    tuple((instance_id, vardata.value, vardata.stale) \
512                          for instance_id, vardata in \
513                          iteritems(scenario_instance.\
514                                    _PHInstanceSymbolMaps[Var].bySymbol) \
515                          if vardata.fixed)
516
517                self._cached_solutions.\
518                    setdefault(cache_id,{})[scenario_name] = \
519                        (self._solver_results[scenario_name],
520                         fixed_results)
521
522    #
523    # when bundling, form the extensive form binding instances given
524    # the current scenario tree specification.  unless bundles are
525    # dynamic, only needs to be invoked once, before PH iteration
526    # 0. otherwise, needs to be invoked each time the bundle structure
527    # is redefined.
528    #
529    # the resulting binding instances are stored in:
530    # self._bundle_extensive_form_map.  the scenario instances
531    # associated with a bundle are stored in:
532    # self._bundle_scenario_instance_map.
533    #
534
535    def _form_bundle_binding_instances(self, preprocess_objectives=True):
536
537        start_time = time.time()
538        if self._verbose:
539            print("Forming binding instances for all scenario bundles")
540
541        self._bundle_binding_instance_map.clear()
542        self._bundle_scenario_instance_map.clear()
543
544        if self._scenario_tree.contains_bundles() is False:
545            raise RuntimeError("Failed to create binding instances for scenario "
546                               "bundles - no scenario bundles are defined!")
547
548        for scenario_bundle in self._scenario_tree._scenario_bundles:
549
550            if self._verbose:
551                print("Creating binding instance for scenario bundle=%s"
552                      % (scenario_bundle._name))
553
554            self._bundle_scenario_instance_map[scenario_bundle._name] = {}
555            for scenario_name in scenario_bundle._scenario_names:
556                self._bundle_scenario_instance_map[scenario_bundle._name]\
557                    [scenario_name] = self._instances[scenario_name]
558
559            # IMPORTANT: The bundle variable IDs must be idential to
560            #            those in the parent scenario tree - this is
561            #            critical for storing results, which occurs at
562            #            the full-scale scenario tree.
563
564            # WARNING: THIS IS A PURE HACK - WE REALLY NEED TO CALL
565            #          THIS WHEN WE CONSTRUCT THE BUNDLE SCENARIO
566            #          TREE.  AS IT STANDS, THIS MUST BE DONE BEFORE
567            #          CREATING THE EF INSTANCE.
568
569            scenario_bundle._scenario_tree.linkInInstances(
570                self._instances,
571                create_variable_ids=False,
572                master_scenario_tree=self._scenario_tree,
573                initialize_solution_data=False)
574
575            bundle_ef_instance = create_ef_instance(
576                scenario_bundle._scenario_tree,
577                ef_instance_name = scenario_bundle._name,
578                verbose_output = self._verbose)
579
580            self._bundle_binding_instance_map[scenario_bundle._name] = \
581                bundle_ef_instance
582
583            # Adding the ph objective terms to the bundle
584            bundle_ef_objective_data = \
585                find_active_objective(bundle_ef_instance, safety_checks=True)
586
587            # augment the EF objective with the PH penalty terms for
588            # each composite scenario.
589            for scenario_name in scenario_bundle._scenario_names:
590                proximal_expression_component = \
591                    self._problem_states.\
592                    ph_objective_proximal_expressions[scenario_name][0]
593                weight_expression_component = \
594                    self._problem_states.\
595                    ph_objective_weight_expressions[scenario_name][0]
596                scenario = self._scenario_tree._scenario_map[scenario_name]
597                bundle_ef_objective_data.expr += \
598                    (scenario._probability / scenario_bundle._probability) * \
599                    proximal_expression_component
600                bundle_ef_objective_data.expr += \
601                    (scenario._probability / scenario_bundle._probability) * \
602                    weight_expression_component
603
604            if self._solver.problem_format == ProblemFormat.nl:
605                    ampl_preprocess_block_objectives(bundle_ef_instance)
606                    ampl_preprocess_block_constraints(bundle_ef_instance)
607            else:
608                var_id_map = {}
609                canonical_preprocess_block_objectives(bundle_ef_instance,
610                                                      var_id_map)
611                canonical_preprocess_block_constraints(bundle_ef_instance,
612                                                       var_id_map)
613
614        end_time = time.time()
615
616        if self._output_times:
617            print("Scenario bundle construction time=%.2f seconds"
618                  % (end_time - start_time))
619
620    def _destory_bundle_binding_instances(self):
621
622        for scenario in self._scenario_tree._scenarios:
623
624            if scenario._instance.parent_block() is not None:
625
626                scenario._instance.parent_block().del_component(scenario._instance)
627
628
629        self._bundle_binding_instance_map.clear()
630        self._bundle_scenario_instance_map.clear()
631
632    def add_ph_objective_proximal_terms(self):
633
634        start_time = time.time()
635
636        for instance_name, instance in iteritems(self._instances):
637
638            if not self._problem_states.\
639                  has_ph_objective_proximal_terms[instance_name]:
640                expression_component, proximal_expression = \
641                    add_ph_objective_proximal_terms(instance_name,
642                                                    instance,
643                                                    self._scenario_tree,
644                                                    self._linearize_nonbinary_penalty_terms,
645                                                    self._retain_quadratic_binary_terms)
646
647                self._problem_states.\
648                    ph_objective_proximal_expressions[instance_name] = \
649                        (expression_component, proximal_expression)
650                self._problem_states.\
651                    has_ph_objective_proximal_terms[instance_name] = True
652                # Flag the preprocessor
653                self._problem_states.objective_updated[instance_name] = True
654
655        end_time = time.time()
656
657        if self._output_times:
658            print("Add PH objective proximal terms time=%.2f seconds"
659                  % (end_time - start_time))
660
661    def activate_ph_objective_proximal_terms(self):
662
663        for instance_name, instance in iteritems(self._instances):
664
665            if not self._problem_states.\
666                  has_ph_objective_proximal_terms[instance_name]:
667                expression_component, expression = \
668                    self._problem_states.\
669                    ph_objective_proximal_expressions[instance_name]
670                expression_component.value = expression
671                self._problem_states.\
672                    has_ph_objective_proximal_terms[instance_name] = True
673                # Flag the preprocessor
674                self._problem_states.objective_updated[instance_name] = True
675
676    def deactivate_ph_objective_proximal_terms(self):
677
678        for instance_name, instance in iteritems(self._instances):
679
680            if self._problem_states.\
681                  has_ph_objective_proximal_terms[instance_name]:
682                self._problem_states.\
683                    ph_objective_proximal_expressions[instance_name][0].value = 0.0
684                self._problem_states.\
685                    has_ph_objective_proximal_terms[instance_name] = False
686                # Flag the preprocessor
687                self._problem_states.objective_updated[instance_name] = True
688
689    def add_ph_objective_weight_terms(self):
690
691        start_time = time.time()
692
693        for instance_name, instance in iteritems(self._instances):
694
695            if not self._problem_states.\
696                  has_ph_objective_weight_terms[instance_name]:
697                expression_component, expression = \
698                    add_ph_objective_weight_terms(instance_name,
699                                                  instance,
700                                                  self._scenario_tree)
701
702                self._problem_states.\
703                    ph_objective_weight_expressions[instance_name] = \
704                        (expression_component, expression)
705                self._problem_states.\
706                    has_ph_objective_weight_terms[instance_name] = True
707                # Flag the preprocessor
708                self._problem_states.objective_updated[instance_name] = True
709
710        end_time = time.time()
711
712        if self._output_times:
713            print("Add PH objective weight terms time=%.2f seconds"
714                  % (end_time - start_time))
715
716    def activate_ph_objective_weight_terms(self):
717
718        for instance_name, instance in iteritems(self._instances):
719
720            if not self._problem_states.\
721               has_ph_objective_weight_terms[instance_name]:
722                expression_component, expression = \
723                    self._problem_states.\
724                    ph_objective_weight_expressions[instance_name]
725                expression_component.value = expression
726                self._problem_states.\
727                    has_ph_objective_weight_terms[instance_name] = True
728                # Flag the preprocessor
729                self._problem_states.objective_updated[instance_name] = True
730
731    def deactivate_ph_objective_weight_terms(self):
732
733        for instance_name, instance in iteritems(self._instances):
734
735            if self._problem_states.\
736               has_ph_objective_weight_terms[instance_name]:
737                self._problem_states.\
738                    ph_objective_weight_expressions[instance_name][0].value = 0.0
739                self._problem_states.\
740                    has_ph_objective_weight_terms[instance_name] = False
741                # Flag the preprocessor
742                self._problem_states.objective_updated[instance_name] = True
743
744    def _push_w_to_instances(self):
745
746        for scenario in self._scenario_tree._scenarios:
747            scenario.push_w_to_instance()
748            # The objectives are always updated when the weight params are updated
749            # and weight terms exist
750            if self._problem_states.has_ph_objective_weight_terms[scenario._name]:
751                # Flag the preprocessor
752                self._problem_states.objective_updated[scenario._name] = True
753
754    def _push_rho_to_instances(self):
755
756        for scenario in self._scenario_tree._scenarios:
757            scenario.push_rho_to_instance()
758            # The objectives are always updated when the rho params are updated
759            # and the proximal terms exist
760            if self._problem_states.has_ph_objective_proximal_terms[scenario._name]:
761                # Flag the preprocessor
762                self._problem_states.objective_updated[scenario._name] = True
763
764    def _push_xbar_to_instances(self):
765
766        for stage in self._scenario_tree._stages[:-1]:
767            for tree_node in stage._tree_nodes:
768                tree_node.push_xbar_to_instances()
769        for scenario in self._scenario_tree._scenarios:
770            # The objectives are always updated when the xbar params are updated
771            # and proximal terms exist
772            if self._problem_states.has_ph_objective_proximal_terms[scenario._name]:
773                # Flag the preprocessor
774                self._problem_states.objective_updated[scenario._name] = True
775
776    def _push_fix_queue_to_instances(self):
777
778        for tree_node in self._scenario_tree._tree_nodes:
779
780            if len(tree_node._fix_queue):
781
782                """
783                some_fixed = tree_node.has_fixed_in_queue()
784                some_freed = tree_node.has_freed_in_queue()
785                # flag the preprocessor
786                if some_fixed or some_freed:
787                    for scenario in tree_node._scenarios:
788                        scenario_name = scenario._name
789                        self._problem_states.\
790                            fixed_variables[scenario_name] |= some_fixed
791                        self._problem_states.\
792                            freed_variables[scenario_name] |= some_freed
793                """
794                for scenario in tree_node._scenarios:
795                    scenario_name = scenario._name
796                    for variable_id, (fixed_status, new_value) in \
797                          iteritems(tree_node._fix_queue):
798                        variable_name, index = tree_node._variable_ids[variable_id]
799                        if fixed_status == tree_node.VARIABLE_FREED:
800                            self._problem_states.\
801                                freed_variables[scenario_name].\
802                                append((variable_name, index))
803                        elif fixed_status == tree_node.VARIABLE_FIXED:
804                            self._problem_states.\
805                                fixed_variables[scenario_name].\
806                                append((variable_name, index))
807                tree_node.push_fix_queue_to_instances()
808
809    def _push_all_node_fixed_to_instances(self):
810
811        for tree_node in self._scenario_tree._tree_nodes:
812
813            tree_node.push_all_fixed_to_instances()
814
815            # flag the preprocessor
816            for scenario in tree_node._scenarios:
817
818                for variable_id in tree_node._fixed:
819
820                    self._problem_states.\
821                        fixed_variables[scenario._name].\
822                        append(tree_node._variable_ids[variable_id])
823
824    #
825    # when linearizing the PH objective, PHQUADPENALTY* variables are
826    # introduced. however, the inclusion / presence of these variables
827    # in warm-start files leads to infeasible MIP starts. thus, we
828    # want to flag their value as None in all scenario instances prior
829    # to performing scenario sub-problem solves.
830    #
831
832    def _reset_instance_linearization_variables(self):
833
834       for scenario_name, scenario_instance in iteritems(self._instances):
835           if self._problem_states.has_ph_variables(scenario_name):
836               reset_linearization_variables(scenario_instance)
837
838    def form_ph_linearized_objective_constraints(self):
839
840        start_time = time.time()
841
842        for instance_name, instance in iteritems(self._instances):
843
844            if self._problem_states.has_ph_objective_proximal_terms[instance_name]:
845                new_attrs = form_linearized_objective_constraints(
846                    instance_name,
847                    instance,
848                    self._scenario_tree,
849                    self._linearize_nonbinary_penalty_terms,
850                    self._breakpoint_strategy,
851                    self._integer_tolerance)
852
853                self._problem_states.ph_constraints[instance_name].extend(new_attrs)
854                self._problem_states.ph_constraints_updated[instance_name] = True
855
856        end_time = time.time()
857
858        if self._output_times:
859            print("PH linearized objective constraint formation "
860                  "time=%.2f seconds" % (end_time - start_time))
861
862    #
863    # a utility to perform preprocessing on all scenario instances, on
864    # an as-needed basis.  queries the instance modification indicator
865    # attributes on the ProgressiveHedging (self) object. intended to
866    # be invoked before each iteration of PH, just before scenario
867    # solves.
868    #
869
870    def _preprocess_scenario_instances(self, ignore_bundles=False):
871
872        start_time = time.time()
873
874        if (not self._scenario_tree.contains_bundles()) or ignore_bundles:
875
876            for scenario_name, scenario_instance in iteritems(self._instances):
877
878                preprocess_scenario_instance(
879                    scenario_instance,
880                    self._problem_states.fixed_variables[scenario_name],
881                    self._problem_states.freed_variables[scenario_name],
882                    self._problem_states.user_constraints_updated[scenario_name],
883                    self._problem_states.ph_constraints_updated[scenario_name],
884                    self._problem_states.ph_constraints[scenario_name],
885                    self._problem_states.objective_updated[scenario_name],
886                    not self._write_fixed_variables,
887                    self._solver)
888
889                # We've preprocessed the instance, reset the relevant flags
890                self._problem_states.clear_update_flags(scenario_name)
891                self._problem_states.clear_fixed_variables(scenario_name)
892                self._problem_states.clear_freed_variables(scenario_name)
893
894        else:
895
896            for scenario_bundle_name, bundle_ef_instance in iteritems(
897                    self._bundle_binding_instance_map):
898
899                # Until proven otherwise
900                preprocess_bundle_objective = False
901                preprocess_bundle_constraints = False
902
903                for scenario_name in \
904                      self._bundle_scenario_instance_map[scenario_bundle_name]:
905
906                    scenario_instance = self._instances[scenario_name]
907                    fixed_vars = self._problem_states.fixed_variables[scenario_name]
908                    freed_vars = self._problem_states.freed_variables[scenario_name]
909                    objective_updated = \
910                        self._problem_states.objective_updated[scenario_name]
911
912                    if objective_updated:
913                        preprocess_bundle_objective = True
914                    # TODO
915                    """
916                    if (fixed_vars or freed_vars) and \
917                       (not self._write_fixed_variables):
918                    """
919                    if (len(fixed_vars) > 0 or len(freed_vars) > 0) and \
920                       (not self._write_fixed_variables):
921                        preprocess_bundle_objective = True
922                        preprocess_bundle_constraints = True
923
924                    preprocess_scenario_instance(
925                        scenario_instance,
926                        fixed_vars,
927                        freed_vars,
928                        self._problem_states.\
929                            user_constraints_updated[scenario_name],
930                        self._problem_states.ph_constraints_updated[scenario_name],
931                        self._problem_states.ph_constraints[scenario_name],
932                        objective_updated,
933                        not self._write_fixed_variables,
934                        self._solver)
935
936                    # We've preprocessed the instance, reset the relevant flags
937                    self._problem_states.clear_update_flags(scenario_name)
938                    self._problem_states.clear_fixed_variables(scenario_name)
939                    self._problem_states.clear_freed_variables(scenario_name)
940
941                if self._solver.problem_format == ProblemFormat.nl:
942                    if preprocess_bundle_objective:
943                        ampl_preprocess_block_objectives(bundle_ef_instance)
944                    if preprocess_bundle_constraints:
945                        ampl_preprocess_block_constraints(bundle_ef_instance)
946                else:
947                    var_id_map = {}
948                    if preprocess_bundle_objective:
949                        canonical_preprocess_block_objectives(bundle_ef_instance,
950                                                              var_id_map)
951                    if preprocess_bundle_constraints:
952                        canonical_preprocess_block_constraints(bundle_ef_instance,
953                                                               var_id_map)
954        end_time = time.time()
955
956        if self._output_times:
957            print("Scenario instance preprocessing time=%.2f seconds"
958                  % (end_time - start_time))
959
960    #
961    # create PH weight and xbar vectors, on a per-scenario basis, for
962    # each variable that is not in the final stage, i.e., for all
963    # variables that are being blended by PH. the parameters are
964    # created in the space of each scenario instance, so that they can
965    # be directly and automatically incorporated into the
966    # (appropriately modified) objective function.
967    #
968
969    def _create_scenario_ph_parameters(self):
970
971        create_nodal_ph_parameters(self._scenario_tree)
972
973        for instance_name, instance in iteritems(self._instances):
974            new_penalty_variable_names = create_ph_parameters(
975                instance,
976                self._scenario_tree,
977                self._rho,
978                self._linearize_nonbinary_penalty_terms)
979
980            if new_penalty_variable_names != []:
981                self._problem_states.ph_variables[instance_name].\
982                    extend(new_penalty_variable_names)
983
984    #
985    # a pair of utilities intended for folks who are brave enough to
986    # script rho setting in a python file.
987    #
988
989    def _rho_check(self, tree_node, variable_id):
990
991        if not variable_id in tree_node._standard_variable_ids:
992            # Generate a helpful error message
993            if variable_id in tree_node._derived_variable_ids:
994                variable_name, index = tree_node._variable_ids[variable_id]
995                raise ValueError("Cannot access rho for variable '%s' with scenario "
996                                 "tree id '%s' on tree node '%s'. The variable is "
997                                 "derived and therefore exclued from nonanticipativity "
998                                 "conditions." % (variable_name+indexToString(index),
999                                                  variable_id,
1000                                                  tree_node._name))
1001            # search the other tree nodes
1002            for other_tree_node in self._scenario_tree._tree_nodes:
1003                if variable_id in other_tree_node._variable_ids:
1004                    variable_name, index = other_tree_node._variable_ids[variable_id]
1005                    raise ValueError("Cannot access rho for variable '%s' with scenario "
1006                                     "tree id '%s' on tree node '%s' because the variable "
1007                                     "is tracked by a different tree node (%s)."
1008                                     % (variable_name+indexToString(index),
1009                                        variable_id,
1010                                        tree_node._name,
1011                                        other_tree_node._name))
1012            raise ValueError("Invalid scenario tree id '%s' for accessing rho. "
1013                             "No tree nodes were found with an associated "
1014                             "instance variable having that id." % (variable_id))
1015
1016
1017    # NOTE: rho_expression can be Pyomo expression, or a constant
1018    #       float/int. either way, the underlying value will be
1019    #       extracted via a value() call...
1020    def setRhoAllScenarios(self, tree_node, variable_id, rho_expression):
1021
1022        self._rho_check(tree_node, variable_id)
1023
1024        new_rho_value = value(rho_expression)
1025
1026        for scenario in tree_node._scenarios:
1027
1028            scenario._rho[tree_node._name][variable_id] = new_rho_value
1029
1030    def setRhoOneScenario(self, tree_node, scenario, variable_id, rho_expression):
1031
1032        self._rho_check(tree_node, variable_id)
1033
1034        scenario._rho[tree_node._name][variable_id] = value(rho_expression)
1035
1036    def getRhoOneScenario(self, tree_node, scenario, variable_id):
1037
1038        self._rho_check(tree_node, variable_id)
1039
1040        return scenario._rho[tree_node._name][variable_id]
1041
1042    #
1043    # a utility intended for folks who are brave enough to script
1044    # variable bounds setting in a python file.
1045    #
1046    def setVariableBoundsAllScenarios(self,
1047                                      tree_node,
1048                                      variable_id,
1049                                      lower_bound,
1050                                      upper_bound):
1051
1052        for scenario in tree_node._scenarios:
1053            vardata = scenario._instance.\
1054                      _ScenarioTreeSymbolMap.getObject(variable_id)
1055            vardata.setlb(lower_bound)
1056            vardata.setub(upper_bound)
1057
1058    def setVariableBoundsOneScenario(self,
1059                                     tree_node,
1060                                     scenario,
1061                                     variable_id,
1062                                     lower_bound,
1063                                     upper_bound):
1064
1065        vardata = scenario._instance._ScenarioTreeSymbolMap.getObject(variable_id)
1066        vardata.setlb(lower_bound)
1067        vardata.setub(upper_bound)
1068
1069    #
1070    # a utility intended for folks who are brave enough to script
1071    # variable bounds setting in a python file.  same functionality as
1072    # above, but applied to all indicies of the variable, in all
1073    # scenarios.
1074    #
1075    """
1076    def setVariableBoundsAllIndicesAllScenarios(self, variable_name, lower_bound, upper_bound):
1077
1078        if isinstance(lower_bound, float) is False:
1079            raise ValueError("Lower bound supplied to PH method setVariableBoundsAllIndiciesAllScenarios for variable="+variable_name+" must be a constant; value supplied="+str(lower_bound))
1080
1081        if isinstance(upper_bound, float) is False:
1082            raise ValueError("Upper bound supplied to PH method setVariableBoundsAllIndicesAllScenarios for variable="+variable_name+" must be a constant; value supplied="+str(upper_bound))
1083
1084        for instance_name, instance in iteritems(self._instances):
1085
1086            variable = instance.find_component(variable_name)
1087            for index in variable:
1088                variable[index].setlb(lower_bound)
1089                variable[index].setub(upper_bound)
1090    """
1091
1092class ProgressiveHedging(_PHBase):
1093
1094    def set_dual_mode(self):
1095
1096        self._dual_mode = True
1097
1098    def primal_mode(self):
1099
1100        self._set_dual_mode = False
1101
1102    def save_solution(self, label="ph"):
1103
1104        if self._solution_plugins is not None:
1105
1106            for plugin in self._solution_plugins:
1107
1108                plugin.write(self._scenario_tree, label)
1109
1110    def release_components(self):
1111
1112        if not self._initialized:
1113
1114            return
1115
1116        if isinstance(self._solver_manager,
1117                      pyomo.solvers.plugins.smanager.\
1118                      phpyro.SolverManager_PHPyro):
1119
1120            release_phsolverservers(self)
1121
1122        if self._solver is not None:
1123            self._solver.deactivate()
1124            self._solver = None
1125
1126        # cleanup the scenario instances for post-processing -
1127        # ideally, we want to leave them in their original state,
1128        # minus all the PH-specific stuff. we don't do all cleanup
1129        # (leaving things like rhos, etc), but we do clean up
1130        # constraints, as that really hoses up the ef writer.
1131        self._cleanup_scenario_instances()
1132        self._clear_bundle_instances()
1133
1134        self._initialized = False
1135
1136    def activate_ph_objective_proximal_terms(self):
1137
1138        start_time = time.time()
1139
1140        _PHBase.activate_ph_objective_proximal_terms(self)
1141
1142        if isinstance(self._solver_manager,
1143                      pyomo.solvers.plugins.smanager.\
1144                      phpyro.SolverManager_PHPyro):
1145
1146                activate_ph_objective_proximal_terms(self)
1147
1148        end_time = time.time()
1149
1150        if self._output_times:
1151            print("Activate PH objective proximal terms time=%.2f seconds"
1152                  % (end_time - start_time))
1153
1154    def deactivate_ph_objective_proximal_terms(self):
1155
1156        start_time = time.time()
1157
1158        _PHBase.deactivate_ph_objective_proximal_terms(self)
1159
1160        if isinstance(self._solver_manager,
1161                      pyomo.solvers.plugins.smanager.\
1162                      phpyro.SolverManager_PHPyro):
1163
1164            deactivate_ph_objective_proximal_terms(self)
1165
1166        end_time = time.time()
1167
1168        if self._output_times:
1169            print("Deactivate PH objective proximal terms time=%.2f seconds"
1170                  % (end_time - start_time))
1171
1172    def activate_ph_objective_weight_terms(self):
1173
1174        start_time = time.time()
1175
1176        _PHBase.activate_ph_objective_weight_terms(self)
1177
1178        if isinstance(self._solver_manager,
1179                      pyomo.solvers.plugins.smanager.\
1180                      phpyro.SolverManager_PHPyro):
1181
1182            activate_ph_objective_weight_terms(self)
1183
1184        end_time = time.time()
1185
1186        if self._output_times:
1187            print("Activate PH objective weight terms time=%.2f seconds"
1188                  % (end_time - start_time))
1189
1190    def deactivate_ph_objective_weight_terms(self):
1191
1192        start_time = time.time()
1193
1194        _PHBase.deactivate_ph_objective_weight_terms(self)
1195
1196        if isinstance(self._solver_manager,
1197                      pyomo.solvers.plugins.smanager.\
1198                      phpyro.SolverManager_PHPyro):
1199
1200            deactivate_ph_objective_weight_terms(self)
1201
1202        end_time = time.time()
1203
1204        if self._output_times:
1205            print("Deactivate PH objective weight terms time=%.2f "
1206                  "seconds" % (end_time - start_time))
1207
1208    #
1209    # checkpoint the current PH state via pickle'ing. the input iteration count
1210    # simply serves as a tag to create the output file name. everything with the
1211    # exception of the _ph_plugins, _solver_manager, and _solver attributes are
1212    # pickled. currently, plugins fail in the pickle process, which is fine as
1213    # JPW doesn't think you want to pickle plugins (particularly the solver and
1214    # solver manager) anyway. For example, you might want to change those later,
1215    # after restoration - and the PH state is independent of how scenario
1216    # sub-problems are solved.
1217    #
1218
1219    def checkpoint(self, iteration_count):
1220
1221        checkpoint_filename = "checkpoint."+str(iteration_count)
1222
1223        tmp_ph_plugins = self._ph_plugins
1224        tmp_solver_manager = self._solver_manager
1225        tmp_solver = self._solver
1226
1227        self._ph_plugins = None
1228        self._solver_manager = None
1229        self._solver = None
1230
1231        checkpoint_file = open(checkpoint_filename, "w")
1232        pickle.dump(self,checkpoint_file)
1233        checkpoint_file.close()
1234
1235        self._ph_plugins = tmp_ph_plugins
1236        self._solver_manager = tmp_solver_manager
1237        self._solver = tmp_solver
1238
1239        print("Checkpoint written to file="+checkpoint_filename)
1240
1241    def _report_bundle_objectives(self):
1242
1243        assert self._scenario_tree.contains_bundles()
1244
1245        max_name_len = max(len(str(_scenario_bundle._name)) \
1246                           for _scenario_bundle in \
1247                           self._scenario_tree._scenario_bundles)
1248        max_name_len = max((len("Scenario Bundle"), max_name_len))
1249        line = (("  %-"+str(max_name_len)+"s    ") % "Scenario Bundle")
1250        line += ("%-20s %-20s %-20s"
1251                 % ("PH Objective",
1252                    "Cost Objective",
1253                    "PH Objective Gap"))
1254        if self._output_times:
1255            line += (" %-10s" % ("Solve Time"))
1256        print(line)
1257        for scenario_bundle in self._scenario_tree._scenario_bundles:
1258
1259            bundle_gap = self._gaps[scenario_bundle._name]
1260            bundle_objective_value = 0.0
1261            bundle_cost_value = 0.0
1262            for scenario in scenario_bundle._scenario_tree._scenarios:
1263                # The objective must be taken from the scenario
1264                # objects on PH full scenario tree
1265                scenario_objective = \
1266                    self._scenario_tree.get_scenario(scenario._name)._objective
1267                scenario_cost = \
1268                    self._scenario_tree.get_scenario(scenario._name)._cost
1269                # And we need to make sure to use the
1270                # probabilities assigned to scenarios in the
1271                # compressed bundle scenario tree
1272                bundle_objective_value += scenario_objective * \
1273                                          scenario._probability
1274                bundle_cost_value += scenario_cost * \
1275                                     scenario._probability
1276
1277            line = ("  %-"+str(max_name_len)+"s    ")
1278            line += ("%-20.4f %-20.4f")
1279            if (not isinstance(bundle_gap, UndefinedData)) and \
1280               (bundle_gap is not None):
1281                line += (" %-20.4f")
1282            else:
1283                bundle_gap = "None Reported"
1284                line += (" %-20s")
1285            line %= (scenario_bundle._name,
1286                     bundle_objective_value,
1287                     bundle_cost_value,
1288                     bundle_gap)
1289            if self._output_times:
1290                solve_time = self._solve_times.get(scenario_bundle._name)
1291                if (not isinstance(solve_time, UndefinedData)) and \
1292                   (solve_time is not None):
1293                    line += (" %-10.2f"
1294                             % (solve_time))
1295                else:
1296                    line += (" %-10s" % "None Reported")
1297            print(line)
1298        print("")
1299
1300    def _report_scenario_objectives(self):
1301
1302        max_name_len = max(len(str(_scenario._name)) \
1303                           for _scenario in self._scenario_tree._scenarios)
1304        max_name_len = max((len("Scenario"), max_name_len))
1305        line = (("  %-"+str(max_name_len)+"s    ") % "Scenario")
1306        line += ("%-20s %-20s %-20s"
1307                 % ("PH Objective",
1308                    "Cost Objective",
1309                    "PH Objective Gap"))
1310        if self._output_times:
1311            line += (" %-10s" % ("Solve Time"))
1312        print(line)
1313        for scenario in self._scenario_tree._scenarios:
1314            objective_value = scenario._objective
1315            scenario_cost = scenario._cost
1316            gap = self._gaps.get(scenario._name)
1317            line = ("  %-"+str(max_name_len)+"s    ")
1318            line += ("%-20.4f %-20.4f")
1319            if (not isinstance(gap, UndefinedData)) and (gap is not None):
1320                line += (" %-20.4f")
1321            else:
1322                gap = "None Reported"
1323                line += (" %-20s")
1324            line %= (scenario._name,
1325                     objective_value,
1326                     scenario_cost,
1327                     gap)
1328            if self._output_times:
1329                solve_time = self._solve_times.get(scenario._name)
1330                if (not isinstance(solve_time, UndefinedData)) and \
1331                   (solve_time is not None):
1332                    line += (" %-10.2f"
1333                             % (solve_time))
1334                else:
1335                    line += (" %-10s" % "None Reported")
1336            print(line)
1337        print("")
1338
1339    def _push_w_to_instances(self):
1340
1341        if isinstance(self._solver_manager,
1342                      pyomo.solvers.plugins.smanager.\
1343                      phpyro.SolverManager_PHPyro):
1344
1345            transmit_weights(self)
1346
1347        else:
1348
1349            _PHBase._push_w_to_instances(self)
1350
1351    def _push_rho_to_instances(self):
1352
1353        if isinstance(self._solver_manager,
1354                      pyomo.solvers.plugins.smanager.\
1355                      phpyro.SolverManager_PHPyro):
1356
1357            transmit_rhos(self)
1358
1359        else:
1360
1361            _PHBase._push_rho_to_instances(self)
1362
1363    def _push_xbar_to_instances(self):
1364
1365        if isinstance(self._solver_manager,
1366                      pyomo.solvers.plugins.smanager.\
1367                      phpyro.SolverManager_PHPyro):
1368
1369            transmit_xbars(self)
1370
1371        else:
1372
1373            _PHBase._push_xbar_to_instances(self)
1374
1375    def _push_fix_queue_to_instances(self):
1376
1377        if isinstance(self._solver_manager,
1378                      pyomo.solvers.plugins.smanager.\
1379                      phpyro.SolverManager_PHPyro):
1380
1381            transmit_fixed_variables(self)
1382            for tree_node in self._scenario_tree._tree_nodes:
1383                # Note: If the scenario tree doesn't not have
1384                #       instances linked in this method will simply
1385                #       empty the queue into the tree node fixed list
1386                #       without trying to fix instance variables
1387                tree_node.push_fix_queue_to_instances()
1388        else:
1389            _PHBase._push_fix_queue_to_instances(self)
1390
1391    #
1392    # restores the current solutions of all scenario instances that I maintain.
1393    # Additionally, if running with PHPyro, asks solver servers to do the same.
1394    #
1395
1396    def restoreCachedSolutions(self, cache_id, release_cache=False):
1397
1398        if isinstance(self._solver_manager,
1399                      pyomo.solvers.plugins.smanager.\
1400                      phpyro.SolverManager_PHPyro):
1401
1402            restore_cached_scenario_solutions(self, cache_id, release_cache)
1403
1404        _PHBase.restoreCachedSolutions(self, cache_id, release_cache)
1405
1406    def cacheSolutions(self, cache_id=None):
1407
1408        if cache_id is None:
1409            cache_id = str(uuid.uuid4())
1410            while cache_id in self._cached_scenariotree_solutions:
1411                cache_id = str(uuid.uuid4())
1412
1413        if isinstance(self._solver_manager,
1414                      pyomo.solvers.plugins.smanager.\
1415                      phpyro.SolverManager_PHPyro):
1416
1417            cache_scenario_solutions(self, cache_id)
1418
1419        _PHBase.cacheSolutions(self, cache_id)
1420
1421        return cache_id
1422
1423    #
1424    # a simple utility to count the number of continuous and discrete
1425    # variables in a set of instances.  this is based on the number of
1426    # active, non-stale variables in the instances. returns a pair -
1427    # num-discrete, num-continuous.  IMPT: This should obviously only
1428    # be called *after* the scenario instances have been solved -
1429    # otherwise, everything is stale - and you'll get back (0,0).
1430    #
1431    # This method is assumed to be called ONCE immediately following
1432    # the iteration 0 solves. Otherwise, the fixing/freeing of
1433    # variables may interact with the stale flag in different and
1434    # unexpected ways depending on what option we choose for treating
1435    # fixed variable preprocessing.
1436    #
1437
1438    def compute_blended_variable_counts(self):
1439
1440        assert self._called_compute_blended_variable_counts is False
1441        self._called_compute_blended_variable_counts = True
1442
1443        num_continuous_vars = 0
1444        num_discrete_vars = 0
1445
1446        for stage, tree_node, variable_id, var_values, is_fixed, is_stale \
1447              in scenario_tree_node_variables_generator_noinstances(
1448                  self._scenario_tree,
1449                  includeDerivedVariables=False,
1450                  includeLastStage=False):
1451
1452            if not is_stale:
1453                if tree_node.is_variable_discrete(variable_id):
1454                    num_discrete_vars = num_discrete_vars + 1
1455                else:
1456                    num_continuous_vars = num_continuous_vars + 1
1457
1458        return (num_discrete_vars, num_continuous_vars)
1459
1460    #
1461    # ditto above, but count the number of fixed discrete and
1462    # continuous variables.
1463    #
1464
1465    def compute_fixed_variable_counts(self):
1466
1467        num_fixed_continuous_vars = 0
1468        num_fixed_discrete_vars = 0
1469
1470        for stage, tree_node, variable_id, var_values, is_fixed, is_stale \
1471              in scenario_tree_node_variables_generator_noinstances(
1472                  self._scenario_tree,
1473                  includeDerivedVariables=False,
1474                  includeLastStage=False):
1475
1476            if is_fixed:
1477                if tree_node.is_variable_discrete(variable_id):
1478                    num_fixed_discrete_vars = num_fixed_discrete_vars + 1
1479                else:
1480                    num_fixed_continuous_vars = num_fixed_continuous_vars + 1
1481
1482        return (num_fixed_discrete_vars, num_fixed_continuous_vars)
1483
1484    #
1485    # when the quadratic penalty terms are approximated via piecewise
1486    # linear segments, we end up (necessarily) "littering" the
1487    # scenario instances with extra constraints.  these need to and
1488    # should be cleaned up after PH, for purposes of post-PH
1489    # manipulation, e.g., writing the extensive form. equally
1490    # importantly, we need to re-establish the original instance
1491    # objectives.
1492    #
1493
1494    def _cleanup_scenario_instances(self):
1495
1496        for instance_name, instance in iteritems(self._instances):
1497
1498            # Eliminate references to ph constraints
1499            for constraint_name in self._problem_states.ph_constraints[instance_name]:
1500                instance.del_component(constraint_name)
1501            self._problem_states.clear_ph_constraints(instance_name)
1502
1503            # Eliminate references to ph variables
1504            for variable_name in self._problem_states.ph_variables[instance_name]:
1505                instance.del_component(variable_name)
1506            self._problem_states.clear_ph_variables(instance_name)
1507
1508            if hasattr(instance, "skip_canonical_repn"):
1509                del instance.skip_canonical_repn
1510            if hasattr(instance, "gen_obj_ampl_repn"):
1511                del instance.gen_obj_ampl_repn
1512            if hasattr(instance, "gen_con_ampl_repn"):
1513                del instance.gen_con_ampl_repn
1514
1515            if hasattr(instance, "_PHInstanceSymbolMaps"):
1516                del instance._PHInstanceSymbolMaps
1517
1518    #
1519    # a simple utility to extract the first-stage cost statistics, e.g., min, average, and max.
1520    #
1521
1522    def _extract_first_stage_cost_statistics(self):
1523
1524        maximum_value = 0.0
1525        minimum_value = 0.0
1526        sum_values = 0.0
1527        num_values = 0
1528        first_time = True
1529
1530        root_node = self._scenario_tree.findRootNode()
1531        first_stage_name = root_node._stage._name
1532        for scenario in root_node._scenarios:
1533            this_value = scenario._stage_costs[first_stage_name]
1534            # None means not reported by the solver.
1535            if this_value is not None:
1536                num_values += 1
1537                sum_values += this_value
1538                if first_time:
1539                    first_time = False
1540                    maximum_value = this_value
1541                    minimum_value = this_value
1542                else:
1543                    if this_value > maximum_value:
1544                        maximum_value = this_value
1545                    if this_value < minimum_value:
1546                        minimum_value = this_value
1547
1548        if num_values > 0:
1549            sum_values = sum_values / num_values
1550
1551        return minimum_value, sum_values, maximum_value
1552
1553    def __init__(self, options):
1554
1555        _PHBase.__init__(self)
1556
1557        self._options = options
1558
1559        self._solver_manager = None
1560
1561        self._phpyro_worker_jobs_map = {}
1562        self._phpyro_job_worker_map = {}
1563
1564        # (ph iteration, expected cost)
1565        self._cost_history = {}
1566        # (ph iteration with non-anticipative solution, expected cost)
1567        self._incumbent_cost_history = {}
1568        # key in the above dictionary of best solution
1569        self._best_incumbent_key = None
1570        # location in cache of best incumbent solution
1571        self._incumbent_cache_id = 'incumbent'
1572
1573
1574        # Make sure we don't call a method more than once
1575        self._called_compute_blended_variable_counts = False
1576
1577        # Augment the code where necessary to run the dual ph algorithm
1578        self._dual_mode = False
1579
1580        # Define the default configuration for what variable
1581        # values to include inside interation k phsolverserver
1582        # results. See phsolverserverutils for more information.
1583        # The default case sends the minimal set of information
1584        # needed to update ph objective parameters, which would
1585        # be nonleaf node ph variables that are not stale.
1586        # Plugins like phhistoryextension modify this flag to
1587        # include more information (e.g., leaf stage values)
1588        # ** NOTE **: If we do not preprocess fixed variables (default behavior),
1589        #             stale=True is equivalent to extraneous variables declared
1590        #             on the model and scenario tree that are never used (fixed or not).
1591        #             When we do preprocess fixed variables, the stale=True applies
1592        #             to the above case as well as fixed variables (which become stale).
1593        #             ...just something to keep in mind
1594        self._phpyro_variable_transmission_flags = \
1595            TransmitType.nonleaf_stages | \
1596            TransmitType.derived | \
1597            TransmitType.blended
1598
1599        self._overrelax = False
1600        # a default, global value for nu. 0 indicates unassigned.
1601        self._nu = 0.0
1602        # filename for the modeler to set rho on a per-variable or
1603        # per-scenario basis.
1604        self._rho_setter = None
1605        # filename for the modeler to collect aggregate scenario data
1606        self._aggregate_getter = None
1607        # filename for the modeler to set rho on a per-variable basis,
1608        # after all scenarios are available.
1609        self._bound_setter = None
1610        self._max_iterations = 0
1611        self._async = False
1612        self._async_buffer_len = 1
1613
1614        # it may be the case that some plugins think they can do a
1615        # better job of weight updates than PH - and it might even be
1616        # true! if so, set this flag to False and this class will not
1617        # invoke the update_weights() method.
1618        self._ph_weight_updates_enabled = True
1619
1620        # PH reporting parameters
1621        # do I report solutions after each PH iteration?
1622        self._report_solutions = False
1623        # do I report PH weights prior to each PH iteration?
1624        self._report_weights = False
1625        # do I report PH rhos prior to each PH iteration?
1626        self._report_rhos_each_iteration = False
1627        # do I report PH rhos prior to PH iteration 1?
1628        self._report_rhos_first_iteration = False
1629        # do I report only variable statistics when outputting
1630        # solutions and weights?
1631        self._report_only_statistics = False
1632        # do I report statistics (via pprint()) for all variables,
1633        # including those whose values equal 0?
1634        self._report_for_zero_variable_values = False
1635        # do I report statistics (via pprint()) for only non-converged
1636        # variables?
1637        self._report_only_nonconverged_variables = False
1638        # when in verbose mode, do I output weights/averages for
1639        # continuous variables?
1640        self._output_continuous_variable_stats = True
1641        self._output_solver_results = False
1642        self._output_scenario_tree_solution = False
1643
1644        #
1645        # PH performance diagnostic parameters and related timing
1646        # parameters.
1647        #
1648
1649        # indicates disabled.
1650        self._profile_memory = 0
1651        self._time_since_last_garbage_collect = time.time()
1652        # units are seconds
1653        self._minimum_garbage_collection_interval = 5
1654
1655        # PH run-time variables
1656        self._current_iteration = 0 # the 'k'
1657
1658        # options for writing solver files / logging / etc.
1659        self._keep_solver_files = False
1660        self._symbolic_solver_labels = False
1661        self._output_solver_logs = False
1662
1663        # string to support suffix specification by callbacks
1664        self._extensions_suffix_list = None
1665
1666        # PH convergence computer/updater.
1667        self._converger = None
1668
1669        # the checkpoint interval - expensive operation, but worth it
1670        # for big models. 0 indicates don't checkpoint.
1671        self._checkpoint_interval = 0
1672
1673
1674        self._ph_plugins = []
1675        self._solution_plugins = []
1676
1677        # PH timing statistics - relative to last invocation.
1678        self._init_start_time = None # for initialization() method
1679        self._init_end_time = None
1680        self._solve_start_time = None # for solve() method
1681        self._solve_end_time = None
1682        # seconds, over course of solve()
1683        self._cumulative_solve_time = 0.0
1684        # seconds, over course of update_xbars()
1685        self._cumulative_xbar_time = 0.0
1686        # seconds, over course of update_weights()
1687        self._cumulative_weight_time = 0.0
1688
1689        # do I disable warm-start for scenario sub-problem solves
1690        # during PH iterations >= 1?
1691        self._disable_warmstarts = False
1692
1693        # PH maintains a mipgap that is applied to each scenario solve
1694        # that is performed.  this attribute can be changed by PH
1695        # extensions, and the change will be applied on all subsequent
1696        # solves - until it is modified again. the default is None,
1697        # indicating unassigned.
1698        self._mipgap = None
1699
1700        # we only store these temporarily...
1701        scenario_solver_options = None
1702
1703        # process the keyword options
1704        self._flatten_expressions                 = options.flatten_expressions
1705        self._max_iterations                      = options.max_iterations
1706        self._overrelax                           = options.overrelax
1707        self._nu                                  = options.nu
1708        self._async                               = options.async
1709        self._async_buffer_len                    = options.async_buffer_len
1710        self._rho                                 = options.default_rho
1711        self._rho_setter                          = options.rho_cfgfile
1712        self._aggregate_getter                    = options.aggregate_cfgfile
1713        self._bound_setter                        = options.bounds_cfgfile
1714        self._solver_type                         = options.solver_type
1715        self._solver_io                           = options.solver_io
1716        scenario_solver_options                   = options.scenario_solver_options
1717        self._handshake_with_phpyro               = options.handshake_with_phpyro
1718        self._mipgap                              = options.scenario_mipgap
1719        self._write_fixed_variables               = options.write_fixed_variables
1720        self._keep_solver_files                   = options.keep_solver_files
1721        self._symbolic_solver_labels              = options.symbolic_solver_labels
1722        self._output_solver_results               = options.output_solver_results
1723        self._output_solver_logs                  = options.output_solver_logs
1724        self._verbose                             = options.verbose
1725        self._report_solutions                    = options.report_solutions
1726        self._report_weights                      = options.report_weights
1727        self._report_rhos_each_iteration          = options.report_rhos_each_iteration
1728        self._report_rhos_first_iteration         = options.report_rhos_first_iteration
1729        self._report_only_statistics              = options.report_only_statistics
1730        self._report_for_zero_variable_values     = options.report_for_zero_variable_values
1731        self._report_only_nonconverged_variables  = options.report_only_nonconverged_variables
1732        self._output_times                        = options.output_times
1733        self._output_instance_construction_times  = options.output_instance_construction_times
1734        self._disable_warmstarts                  = options.disable_warmstarts
1735        self._retain_quadratic_binary_terms       = options.retain_quadratic_binary_terms
1736        self._linearize_nonbinary_penalty_terms   = options.linearize_nonbinary_penalty_terms
1737        self._breakpoint_strategy                 = options.breakpoint_strategy
1738        self._checkpoint_interval                 = options.checkpoint_interval
1739        self._output_scenario_tree_solution       = options.output_scenario_tree_solution
1740        self._phpyro_transmit_leaf_stage_solution = options.phpyro_transmit_leaf_stage_solution
1741        # clutters up the screen, when we really only care about the
1742        # binaries.
1743        self._output_continuous_variable_stats = not options.suppress_continuous_variable_output
1744
1745        self._objective_sense = options.objective_sense
1746        self._objective_sense_option = options.objective_sense
1747        if hasattr(options, "profile_memory"):
1748            self._profile_memory = options.profile_memory
1749        else:
1750            self._profile_memory = False
1751
1752        if self._phpyro_transmit_leaf_stage_solution:
1753            self._phpyro_variable_transmission_flags |= TransmitType.all_stages
1754
1755        # Note: Default rho has become a required ph input. At this
1756        #       point it seems more natural to make the "-r" or
1757        #       "--default-rho" command-line option required (as
1758        #       contradictory as that sounds) rather than convert it
1759        #       into a positional argument. Unfortunately, optparse
1760        #       doesn't handle "required options", so the most natural
1761        #       location for checking default rho is here.
1762        if (self._rho == ""):
1763            raise ValueError("PH detected an invalid value for default rho: %s. "
1764                             "Use --default-rho=X to specify a positive number X for default rho. "
1765                             "A value of 1.0 is no longer assumed."
1766                             % (self._rho))
1767        if (self._rho == "None"):
1768            self._rho = None
1769            print("***WARNING***: PH is using a default rho value of "
1770                  "None for all blended scenario tree variables. This "
1771                  "will result in error during weight updates following "
1772                  "PH iteration 0 solves unless rho is changed. This "
1773                  "option indicates that a user intends to set rho for "
1774                  "all blended scenario tree variables using a PH extension.")
1775        else:
1776            self._rho = float(self._rho)
1777            if self._rho < 0:
1778                raise ValueError("PH detected an invalid value for default rho: %s. "
1779                                 "Use --default-rho=X to specify a positive number X for default rho. "
1780                                 "A value of 1.0 is no longer assumed."
1781                                 % (self._rho))
1782            elif self._rho == 0:
1783                print("***WARNING***: PH is using a default rho value of "
1784                      "0 for all blended scenario tree variables. This "
1785                      "will effectively disable non-anticipativity "
1786                      "for all variables unless rho is change using a "
1787                      "PH extension")
1788
1789        # cache stuff relating to scenario tree manipulation - the ph
1790        # solver servers may need it.
1791        self._scenario_bundle_specification = options.scenario_bundle_specification
1792        self._create_random_bundles = options.create_random_bundles
1793        self._scenario_tree_random_seed = options.scenario_tree_random_seed
1794
1795        # validate all "atomic" options (those that can be validated independently)
1796        if self._max_iterations < 0:
1797            raise ValueError("Maximum number of PH iterations must be non-negative; value specified=" + str(self._max_iterations))
1798        if self._nu <= 0.0 or self._nu >= 2:
1799            raise ValueError("Value of the nu parameter in PH must be on the interval (0, 2); value specified=" + str(self._nu))
1800        if (self._mipgap is not None) and ((self._mipgap < 0.0) or (self._mipgap > 1.0)):
1801            raise ValueError("Value of the mipgap parameter in PH must be on the unit interval; value specified=" + str(self._mipgap))
1802
1803        #
1804        # validate the linearization (number of pieces) and breakpoint
1805        # distribution parameters.
1806        #
1807        # if a breakpoint strategy is specified without linearization
1808        # enabled, halt and warn the user.
1809        if (self._breakpoint_strategy > 0) and \
1810           (self._linearize_nonbinary_penalty_terms == 0):
1811            raise ValueError("A breakpoint distribution strategy was "
1812                             "specified, but linearization is not enabled!")
1813        if self._linearize_nonbinary_penalty_terms < 0:
1814            raise ValueError("Value of linearization parameter for nonbinary penalty terms must be non-negative; value specified=" + str(self._linearize_nonbinary_penalty_terms))
1815        if self._breakpoint_strategy < 0:
1816            raise ValueError("Value of the breakpoint distribution strategy parameter must be non-negative; value specified=" + str(self._breakpoint_strategy))
1817        if self._breakpoint_strategy > 3:
1818            raise ValueError("Unknown breakpoint distribution strategy specified - valid values are between 0 and 2, inclusive; value specified=" + str(self._breakpoint_strategy))
1819
1820        # validate that callback functions exist in specified modules
1821        self._callback_function = {}
1822        self._mapped_module_name = {}
1823        for ph_attr, callback_name in (("_aggregate_getter","ph_aggregategetter_callback"),
1824                                       ("_rho_setter","ph_rhosetter_callback"),
1825                                       ("_bound_setter","ph_boundsetter_callback")):
1826            module_name = getattr(self,ph_attr)
1827            if module_name is not None:
1828                sys_modules_key, module = load_external_module(module_name)
1829                callback = None
1830                for oname, obj in inspect.getmembers(module):
1831                    if oname == callback_name:
1832                        callback = obj
1833                        break
1834
1835                if callback is None:
1836                    raise ImportError("PH callback with name '%s' could not be found in module file: "
1837                                      "%s" % (callback_name, module_name))
1838                self._callback_function[sys_modules_key] = callback
1839                setattr(self,ph_attr,sys_modules_key)
1840                self._mapped_module_name[sys_modules_key] = module_name
1841
1842        # validate the checkpoint interval.
1843        if self._checkpoint_interval < 0:
1844            raise ValueError("A negative checkpoint interval with value="
1845                             +str(self._checkpoint_interval)+" was "
1846                             "specified in call to PH constructor")
1847
1848        # construct the sub-problem solver.
1849        if self._verbose:
1850            print("Constructing solver type="+self._solver_type)
1851        self._solver = SolverFactory(self._solver_type, solver_io=self._solver_io)
1852        if self._solver == None:
1853            raise ValueError("Unknown solver type=" + self._solver_type + " specified in call to PH constructor")
1854        if self._keep_solver_files:
1855            self._solver.keepfiles = True
1856        self._solver.symbolic_solver_labels = self._symbolic_solver_labels
1857        if len(scenario_solver_options) > 0:
1858            if self._verbose:
1859                print("Initializing scenario sub-problem solver with options="+str(scenario_solver_options))
1860            self._solver.set_options("".join(scenario_solver_options))
1861        if self._output_times:
1862            self._solver._report_timing = True
1863
1864        # a set of all valid PH iteration indicies is generally useful for plug-ins, so create it here.
1865        self._iteration_index_set = Set(name="PHIterations")
1866        for i in range(0,self._max_iterations + 1):
1867            self._iteration_index_set.add(i)
1868
1869        #
1870        # construct the convergence "computer" class.
1871        #
1872        converger = None
1873        # go with the non-defaults first, and then with the default
1874        # (normalized term-diff).
1875        if options.enable_free_discrete_count_convergence:
1876            if self._verbose:
1877                print("Enabling convergence based on a fixed number of discrete variables")
1878            converger = NumFixedDiscreteVarConvergence(convergence_threshold=options.free_discrete_count_threshold)
1879        elif options.enable_termdiff_convergence:
1880            if self._verbose:
1881                print("Enabling convergence based on non-normalized term diff criterion")
1882            converger = TermDiffConvergence(convergence_threshold=options.termdiff_threshold)
1883        else:
1884            converger = NormalizedTermDiffConvergence(convergence_threshold=options.termdiff_threshold)
1885        self._converger = converger
1886
1887        # spit out parameterization if verbosity is enabled
1888        if self._verbose:
1889            print("PH solver configuration: ")
1890            print("   Max iterations="+str(self._max_iterations))
1891            print("   Async mode=" + str(self._async))
1892            print("   Async buffer len=" + str(self._async_buffer_len))
1893            print("   Default global rho=" + str(self._rho))
1894            print("   Over-relaxation enabled="+str(self._overrelax))
1895            if self._overrelax:
1896                print("   Nu=" + self._nu)
1897            if self._aggregate_getter is not None:
1898                print("   Aggregate getter callback file=" + self._aggregate_getter)
1899            if self._rho_setter is not None:
1900                print("   Rho setter callback file=" + self._rho_setter)
1901            if self._bound_setter is not None:
1902                print("   Bound setter callback file=" + self._bound_setter)
1903            print("   Sub-problem solver type='%s'" % str(self._solver_type))
1904            print("   Keep solver files? " + str(self._keep_solver_files))
1905            print("   Output solver results? " + str(self._output_solver_results))
1906            print("   Output solver log? " + str(self._output_solver_logs))
1907            print("   Output times? " + str(self._output_times))
1908            print("   Checkpoint interval="+str(self._checkpoint_interval))
1909
1910    """ Initialize PH with model and scenario data, in preparation for solve().
1911        Constructs and reads instances.
1912    """
1913    def initialize(self,
1914                   scenario_tree=None,
1915                   solver_manager=None,
1916                   ph_plugins=None,
1917                   solution_plugins=None):
1918
1919        self._init_start_time = time.time()
1920
1921        print("Initializing PH")
1922        print("")
1923
1924        if ph_plugins is not None:
1925            self._ph_plugins = ph_plugins
1926
1927        if solution_plugins is not None:
1928            self._solution_plugins = solution_plugins
1929
1930        # The first step in PH initialization is to impose an order on
1931        # the user-defined plugins. Invoking wwextensions and
1932        # phhistoryextensions last ensures that any any changes the
1933        # former makes in order to affect algorithm convergence will
1934        # not be seen by any other plugins until the next iteration
1935        # and so that the latter can properly snapshot any changes
1936        # made by other plugins. All remaining extensions are invoked
1937        # prior to these in random order. The reason is that we're
1938        # being lazy - ideally, this would be the user defined order
1939        # on the command-line.
1940        phboundextensions = \
1941            [plugin for plugin in self._ph_plugins \
1942             if isinstance(plugin,
1943                           pyomo.pysp.plugins.phboundextension.\
1944                           phboundextension)]
1945
1946        convexhullboundextensions = \
1947            [plugin for plugin in self._ph_plugins \
1948             if isinstance(plugin,
1949                           pyomo.pysp.plugins.convexhullboundextension.\
1950                           convexhullboundextension)]
1951
1952        wwextensions = \
1953            [plugin for plugin in self._ph_plugins \
1954             if isinstance(plugin,
1955                           pyomo.pysp.plugins.wwphextension.wwphextension)]
1956
1957        phhistoryextensions = \
1958            [plugin for plugin in self._ph_plugins \
1959             if isinstance(plugin,
1960                           pyomo.pysp.plugins.phhistoryextension.\
1961                           phhistoryextension)]
1962
1963        userdefinedextensions = []
1964        for plugin in self._ph_plugins:
1965            if not (isinstance(plugin,
1966                               pyomo.pysp.plugins.wwphextension.wwphextension) or \
1967                    isinstance(plugin,
1968                               pyomo.pysp.plugins.phhistoryextension.\
1969                               phhistoryextension) or \
1970                    isinstance(plugin,
1971                               pyomo.pysp.plugins.phboundextension.\
1972                               phboundextension) or \
1973                    isinstance(plugin,
1974                               pyomo.pysp.plugins.convexhullboundextension.\
1975                               convexhullboundextension)):
1976                userdefinedextensions.append(plugin)
1977
1978        # note that the order of plugin invocation is important. the history
1979        # extension goes last, to capture all modifications made by all plugins.
1980        # user-defined extensions should otherwise go after the built-in plugins.
1981        ph_plugins = []
1982        ph_plugins.extend(convexhullboundextensions)
1983        ph_plugins.extend(phboundextensions)
1984        ph_plugins.extend(wwextensions)
1985        ph_plugins.extend(userdefinedextensions)
1986        ph_plugins.extend(phhistoryextensions)
1987        self._ph_plugins = ph_plugins
1988
1989        # let plugins know if they care.
1990        if self._verbose:
1991            print("Invoking pre-initialization PH plugins")
1992        for plugin in self._ph_plugins:
1993            plugin.pre_ph_initialization(self)
1994
1995        if scenario_tree is None:
1996            raise ValueError("A scenario tree must be supplied to the "
1997                             "PH initialize() method")
1998
1999        if solver_manager is None:
2000            raise ValueError("A solver manager must be supplied to "
2001                             "the PH initialize() method")
2002
2003        # Eventually some of these might really become optional
2004        self._scenario_tree = scenario_tree
2005        self._solver_manager = solver_manager
2006
2007        self._converger.reset()
2008
2009        isPHPyro =  isinstance(self._solver_manager,
2010                               pyomo.solvers.plugins.\
2011                               smanager.phpyro.SolverManager_PHPyro)
2012
2013        initialization_action_handles = []
2014        if isPHPyro:
2015
2016            if self._verbose:
2017                print("Broadcasting requests to initialize PH solver servers")
2018
2019            initialization_action_handles.extend(initialize_ph_solver_servers(self))
2020
2021            if self._verbose:
2022                print("PH solver server initialization requests successfully transmitted")
2023
2024        else:
2025            # gather the scenario tree instances into
2026            # the self._instances dictionary and
2027            # tag appropriate preprocessing flags
2028            self._setup_scenario_instances()
2029
2030        # let plugins know if they care - this callback point allows
2031        # users to create/modify the original scenario instances
2032        # and/or the scenario tree prior to creating PH-related
2033        # parameters, variables, and the like.
2034        post_instance_plugin_callback_start_time = time.time()
2035        for plugin in self._ph_plugins:
2036            plugin.post_instance_creation(self)
2037        post_instance_plugin_callback_end_time = time.time()
2038        if self._output_times:
2039            print("PH post-instance plugin callback time=%.2f seconds"
2040                  % (post_instance_plugin_callback_end_time - \
2041                     post_instance_plugin_callback_start_time))
2042
2043        if not isPHPyro:
2044
2045            # create ph-specific parameters (weights, xbar, etc.) for
2046            # each instance.
2047            if self._verbose:
2048                print("Creating weight, average, and rho parameter "
2049                      "vectors for scenario instances")
2050            scenario_ph_parameters_start_time = time.time()
2051            self._create_scenario_ph_parameters()
2052            scenario_ph_parameters_end_time = time.time()
2053            if self._output_times:
2054                print("PH parameter vector construction time=%.2f seconds"
2055                      % (scenario_ph_parameters_end_time - \
2056                         scenario_ph_parameters_start_time))
2057
2058            # create symbol maps for easy storage/transmission of
2059            # variable values
2060            # NOTE: Not sure of the timing overhead that comes with
2061            #       this, but it's likely we can make this optional
2062            #       when we are not running parallel PH.
2063            if self._verbose:
2064                print("Creating deterministic SymbolMaps for scenario instances")
2065            scenario_ph_symbol_maps_start_time = time.time()
2066            # Define for what components we generate symbols
2067            symbol_ctypes = (Var, Suffix)
2068            self._create_instance_symbol_maps(symbol_ctypes)
2069            scenario_ph_symbol_maps_end_time = time.time()
2070            if self._output_times:
2071                print("PH SymbolMap creation time=%.2f seconds"
2072                      % (scenario_ph_symbol_maps_end_time - \
2073                         scenario_ph_symbol_maps_start_time))
2074
2075            # form the ph objective weight and proximal expressions
2076            # Note: The Expression objects for the weight and proximal
2077            #       terms will be added to the instances objectives
2078            #       but will be assigned values of 0.0, so that the
2079            #       original objective function form is maintained.
2080            #       The purpose is so that we can use this shared
2081            #       Expression object in the bundle binding instance
2082            #       objectives as well when we call
2083            #       _form_bundle_binding_instances a few lines down
2084            #       (so regeneration of bundle objective expressions
2085            #       is no longer required before each iteration k
2086            #       solve).
2087            self.add_ph_objective_weight_terms()
2088            self.deactivate_ph_objective_weight_terms()
2089            self.add_ph_objective_proximal_terms()
2090            self.deactivate_ph_objective_proximal_terms()
2091
2092            # if we have bundles and are not running with PH Pyro, we
2093            # need to create the binding instances - because we are
2094            # responsible for farming the instances out for solution.
2095            if self._scenario_tree.contains_bundles():
2096                self._form_bundle_binding_instances(preprocess_objectives=False)
2097
2098        # If specified, run the user script to collect aggregate
2099        # scenario data. This can slow down PH initialization as
2100        # syncronization across all phsolverservers is required
2101        if self._aggregate_getter is not None:
2102
2103            if isPHPyro:
2104
2105                # Transmit invocation to phsolverservers
2106                print("Transmitting user aggregate callback invocations "
2107                      "to phsolverservers")
2108                if self._scenario_tree.contains_bundles():
2109                    for scenario_bundle in self._scenario_tree._scenario_bundles:
2110                        ah = transmit_external_function_invocation_to_worker(
2111                            self,
2112                            scenario_bundle._name,
2113                            self._mapped_module_name[self._aggregate_getter],
2114                            "ph_aggregategetter_callback",
2115                            invocation_type=InvocationType.\
2116                                            PerScenarioChainedInvocation,
2117                            return_action_handle=True,
2118                            function_args=(self._aggregate_user_data,))
2119                        result = self._solver_manager.wait_for(ah)
2120                        assert len(result) == 1
2121                        self._aggregate_user_data = result[0]
2122
2123                else:
2124                    for scenario in self._scenario_tree._scenarios:
2125                        ah = transmit_external_function_invocation_to_worker(
2126                            self,
2127                            scenario._name,
2128                            self._mapped_module_name[self._aggregate_getter],
2129                            "ph_aggregategetter_callback",
2130                            invocation_type=InvocationType.SingleInvocation,
2131                            return_action_handle=True,
2132                            function_args=(self._aggregate_user_data,))
2133                        result = self._solver_manager.wait_for(ah)
2134                        assert len(result) == 1
2135                        self._aggregate_user_data = result[0]
2136
2137                # Transmit final aggregate state to phsolverservers
2138                print("Broadcasting final aggregate data to phsolverservers")
2139                transmit_external_function_invocation(
2140                    self,
2141                    "pyomo.pysp.ph",
2142                    "assign_aggregate_data",
2143                    invocation_type=InvocationType.SingleInvocation,
2144                    return_action_handles=False,
2145                    function_args=(self._aggregate_user_data,))
2146
2147            else:
2148
2149                print("Executing user aggregate getter callback function")
2150                for scenario in self._scenario_tree._scenarios:
2151                    result = self._callback_function[self._aggregate_getter](
2152                        self,
2153                        self._scenario_tree,
2154                        scenario,
2155                        self._aggregate_user_data)
2156                    assert len(result) == 1
2157                    self._aggregate_user_data = result[0]
2158
2159        # if specified, run the user script to initialize variable
2160        # rhos at their whim.
2161        if self._rho_setter is not None:
2162
2163            if isPHPyro:
2164
2165                # Transmit invocation to phsolverservers
2166                print("Transmitting user rho callback invocations "
2167                      "to phsolverservers")
2168                if self._scenario_tree.contains_bundles():
2169                    for scenario_bundle in self._scenario_tree._scenario_bundles:
2170                        transmit_external_function_invocation_to_worker(
2171                            self,
2172                            scenario_bundle._name,
2173                            self._mapped_module_name[self._rho_setter],
2174                            "ph_rhosetter_callback",
2175                            invocation_type=InvocationType.PerScenarioInvocation,
2176                            return_action_handle=False)
2177                else:
2178                    for scenario in self._scenario_tree._scenarios:
2179                        transmit_external_function_invocation_to_worker(
2180                            self,
2181                            scenario._name,
2182                            self._mapped_module_name[self._rho_setter],
2183                            "ph_rhosetter_callback",
2184                            invocation_type=InvocationType.SingleInvocation,
2185                            return_action_handle=False)
2186
2187                # NOTE: For the time being we rely on the
2188                #       gather_scenario_tree_data call at the end this
2189                #       initialize method in order to collect the
2190                #       finalized rho values
2191
2192            else:
2193
2194                print("Executing user rho setter callback function")
2195                for scenario in self._scenario_tree._scenarios:
2196                    self._callback_function[self._rho_setter](
2197                        self,
2198                        self._scenario_tree,
2199                        scenario)
2200
2201        # if specified, run the user script to initialize variable
2202        # bounds at their whim.
2203        if self._bound_setter is not None:
2204
2205            if isPHPyro:
2206
2207                # Transmit invocation to phsolverservers
2208                print("Transmitting user bound callback invocations to "
2209                      "phsolverservers")
2210                if self._scenario_tree.contains_bundles():
2211                    for scenario_bundle in self._scenario_tree._scenario_bundles:
2212                        transmit_external_function_invocation_to_worker(
2213                            self,
2214                            scenario_bundle._name,
2215                            self._mapped_module_name[self._bound_setter],
2216                            "ph_boundsetter_callback",
2217                            invocation_type=InvocationType.PerScenarioInvocation,
2218                            return_action_handle=False)
2219                else:
2220                    for scenario in self._scenario_tree._scenarios:
2221                        transmit_external_function_invocation_to_worker(
2222                            self,
2223                            scenario._name,
2224                            self._mapped_module_name[self._bound_setter],
2225                            "ph_boundsetter_callback",
2226                            invocation_type=InvocationType.SingleInvocation,
2227                            return_action_handle=False)
2228
2229            else:
2230
2231                print("Executing user bound setter callback function")
2232                for scenario in self._scenario_tree._scenarios:
2233                    self._callback_function[self._bound_setter](
2234                        self,
2235                        self._scenario_tree,
2236                        scenario)
2237
2238        # at this point, the instances are complete - preprocess them!
2239        # BUT: only if the phpyro solver manager isn't in use.
2240        if isPHPyro:
2241
2242            if self._verbose:
2243                print("Broadcasting requests to collect scenario tree "
2244                      "instance data from PH solver servers")
2245
2246            gather_scenario_tree_data(self, initialization_action_handles)
2247
2248            if self._verbose:
2249                print("Scenario tree instance data successfully collected")
2250
2251            if self._verbose:
2252                print("Broadcasting scenario tree id mapping"
2253                      "to PH solver servers")
2254
2255            transmit_scenario_tree_ids(self)
2256
2257            if self._verbose:
2258                print("Scenario tree ids successfully sent")
2259
2260        self._objective_sense = \
2261            self._scenario_tree._scenarios[0]._objective_sense
2262
2263        # indicate that we're ready to run.
2264        self._initialized = True
2265
2266        if self._verbose:
2267            print("PH successfully created model instances for all scenarios")
2268
2269        if self._verbose:
2270            print("PH is successfully initialized")
2271
2272        if self._output_times:
2273            print("Cumulative initialization time=%.2f seconds"
2274                  % (time.time() - self._init_start_time))
2275
2276        # let plugins know if they care.
2277        if self._verbose:
2278            print("Invoking post-initialization PH plugins")
2279        post_ph_initialization_plugin_callback_start_time = time.time()
2280        for plugin in self._ph_plugins:
2281            plugin.post_ph_initialization(self)
2282        post_ph_initialization_plugin_callback_end_time = time.time()
2283        if self._output_times:
2284            print("PH post-initialization plugin callback time=%.2f seconds"
2285                  % (post_ph_initialization_plugin_callback_end_time - \
2286                     post_ph_initialization_plugin_callback_start_time))
2287
2288        if self._output_times:
2289            print("Cumulative PH initialization time=%.2f seconds"
2290                  % (time.time() - self._init_start_time))
2291
2292        self._init_end_time = time.time()
2293        if self._output_times:
2294            print("Overall initialization time=%.2f seconds"
2295                  % (self._init_end_time - self._init_start_time))
2296            print("")
2297
2298    #
2299    # Transmits Solver Options, Queues Solves, and Collects/Loads
2300    # Results... nothing more. All subproblems are expected to be
2301    # fully preprocessed.
2302    #
2303    def solve_subproblems(self, warmstart=False):
2304
2305        iteration_start_time = time.time()
2306
2307        # Preprocess the scenario instances before solving we're
2308        # not using phpyro
2309        if not isinstance(self._solver_manager,
2310                          pyomo.solvers.plugins.smanager.\
2311                          phpyro.SolverManager_PHPyro):
2312            self._preprocess_scenario_instances()
2313
2314        # STEP -1: clear the gap and solve_time dictionaries - we
2315        #          don't have any results yet.
2316        if self._scenario_tree.contains_bundles():
2317            for scenario_bundle in self._scenario_tree._scenario_bundles:
2318                self._gaps[scenario_bundle._name] = undefined
2319                self._solve_times[scenario_bundle._name] = undefined
2320        else:
2321            for scenario in self._scenario_tree._scenarios:
2322                self._gaps[scenario._name] = undefined
2323                self._solve_times[scenario._name] = undefined
2324
2325        # STEP 0: set up all global solver options.
2326        if self._mipgap is not None:
2327            self._solver.options.mipgap = float(self._mipgap)
2328
2329        # if running the phpyro solver server, we need to ship the
2330        # solver options across the pipe.
2331        if isinstance(self._solver_manager,
2332                      pyomo.solvers.plugins.smanager.\
2333                      phpyro.SolverManager_PHPyro):
2334            solver_options = {}
2335            for key in self._solver.options:
2336                solver_options[key]=self._solver.options[key]
2337
2338        # STEP 1: queue up the solves for all scenario sub-problems
2339        # we could use the same names for scenarios and bundles, but
2340        # we are easily confused.
2341        scenario_action_handle_map = {} # maps scenario names to action handles
2342        action_handle_scenario_map = {} # maps action handles to scenario names
2343
2344        bundle_action_handle_map = {} # maps bundle names to action handles
2345        action_handle_bundle_map = {} # maps action handles to bundle names
2346
2347        common_kwds = {
2348            'tee':self._output_solver_logs,
2349            'keepfiles':self._keep_solver_files,
2350            'symbolic_solver_labels':self._symbolic_solver_labels,
2351            'output_fixed_variable_bounds':self._write_fixed_variables}
2352
2353        # TODO: suffixes are not handled equally for
2354        # scenario/bundles/serial/phpyro
2355        if isinstance(self._solver_manager,
2356                      pyomo.solvers.plugins.smanager.phpyro.SolverManager_PHPyro):
2357            common_kwds['solver_options'] = solver_options
2358            common_kwds['solver_suffixes'] = []
2359            common_kwds['warmstart'] = warmstart
2360            common_kwds['variable_transmission'] = \
2361                self._phpyro_variable_transmission_flags
2362        else:
2363            common_kwds['verbose'] = self._verbose
2364
2365        if self._scenario_tree.contains_bundles():
2366
2367            for scenario_bundle in self._scenario_tree._scenario_bundles:
2368
2369                if self._verbose:
2370                    print("Queuing solve for scenario bundle=%s"
2371                          % (scenario_bundle._name))
2372
2373                # and queue it up for solution - have to worry about
2374                # warm-starting here.
2375                new_action_handle = None
2376                if isinstance(self._solver_manager,
2377                              pyomo.solvers.plugins.smanager.\
2378                              phpyro.SolverManager_PHPyro):
2379                    new_action_handle = \
2380                        self._solver_manager.queue(action="solve",
2381                                                   name=scenario_bundle._name,
2382                                                   **common_kwds)
2383                else:
2384
2385                    if (self._output_times is True) and (self._verbose is False):
2386                        print("Solver manager queuing instance=%s"
2387                              % (scenario_bundle._name))
2388
2389                    if warmstart and self._solver.warm_start_capable():
2390                        new_action_handle = \
2391                            self._solver_manager.queue(
2392                                self._bundle_binding_instance_map[\
2393                                    scenario_bundle._name],
2394                                opt=self._solver,
2395                                warmstart=True,
2396                                **common_kwds)
2397                    else:
2398                        new_action_handle = \
2399                            self._solver_manager.queue(
2400                                self._bundle_binding_instance_map[\
2401                                    scenario_bundle._name],
2402                                opt=self._solver,
2403                                **common_kwds)
2404
2405                bundle_action_handle_map[scenario_bundle._name] = new_action_handle
2406                action_handle_bundle_map[new_action_handle] = scenario_bundle._name
2407
2408        else:
2409
2410            for scenario in self._scenario_tree._scenarios:
2411
2412                if self._verbose:
2413                    print("Queuing solve for scenario=%s" % (scenario._name))
2414
2415                # once past iteration 0, there is always a feasible
2416                # solution from which to warm-start.  however, you
2417                # might want to disable warm-start when the solver is
2418                # behaving badly (which does happen).
2419                new_action_handle = None
2420                if isinstance(self._solver_manager,
2421                              pyomo.solvers.plugins.smanager.\
2422                              phpyro.SolverManager_PHPyro):
2423
2424                    new_action_handle = \
2425                        self._solver_manager.queue(action="solve",
2426                                                   name=scenario._name,
2427                                                   **common_kwds)
2428                else:
2429
2430                    instance = scenario._instance
2431
2432                    if (self._output_times is True) and (self._verbose is False):
2433                        print("Solver manager queuing instance=%s"
2434                              % (scenario._name))
2435
2436                    if warmstart and self._solver.warm_start_capable():
2437
2438                        if self._extensions_suffix_list is not None:
2439                            new_action_handle = \
2440                                self._solver_manager.queue(
2441                                    instance,
2442                                    opt=self._solver,
2443                                    warmstart=True,
2444                                    suffixes=self._extensions_suffix_list,
2445                                    **common_kwds)
2446                        else:
2447                            new_action_handle = \
2448                                self._solver_manager.queue(instance,
2449                                                           opt=self._solver,
2450                                                           warmstart=True,
2451                                                           **common_kwds)
2452                    else:
2453
2454                        if self._extensions_suffix_list is not None:
2455                            new_action_handle = \
2456                                self._solver_manager.queue(
2457                                    instance,
2458                                    opt=self._solver,
2459                                    suffixes=self._extensions_suffix_list,
2460                                    **common_kwds)
2461                        else:
2462                            new_action_handle = \
2463                                self._solver_manager.queue(instance,
2464                                                           opt=self._solver,
2465                                                           **common_kwds)
2466
2467                scenario_action_handle_map[scenario._name] = new_action_handle
2468                action_handle_scenario_map[new_action_handle] = scenario._name
2469
2470        # STEP 3: loop for the solver results, reading them and
2471        #         loading them into instances as they are available.
2472        if self._scenario_tree.contains_bundles():
2473
2474            if self._verbose:
2475                print("Waiting for bundle sub-problem solves")
2476
2477            num_results_so_far = 0
2478
2479            while (num_results_so_far < \
2480                   len(self._scenario_tree._scenario_bundles)):
2481
2482                bundle_action_handle = self._solver_manager.wait_any()
2483                bundle_results = \
2484                    self._solver_manager.get_results(bundle_action_handle)
2485                bundle_name = action_handle_bundle_map[bundle_action_handle]
2486
2487                if isinstance(self._solver_manager,
2488                              pyomo.solvers.plugins.smanager.phpyro.\
2489                              SolverManager_PHPyro):
2490
2491                    if self._output_solver_results:
2492                        print("Results for scenario bundle=%s:"
2493                              % (bundle_name))
2494                        print(bundle_results)
2495
2496                    if len(bundle_results) == 0:
2497                        if self._verbose:
2498                            print("Solve failed for scenario bundle=%s; "
2499                                  "no solutions generated" % (bundle_name))
2500                        return bundle_name
2501
2502                    for scenario_name, scenario_solution in \
2503                                      iteritems(bundle_results[0]):
2504                        scenario = self._scenario_tree._scenario_map[scenario_name]
2505                        scenario.update_current_solution(scenario_solution)
2506
2507                    auxilliary_values = bundle_results[2]
2508                    if "gap" in auxilliary_values:
2509                        self._gaps[bundle_name] = auxilliary_values["gap"]
2510
2511                    if "user_time" in auxilliary_values:
2512                        self._solve_times[bundle_name] = \
2513                            auxilliary_values["user_time"]
2514
2515                else:
2516
2517                    # a temporary hack - if results come back from
2518                    # pyro, they won't have a symbol map attached. so
2519                    # create one.
2520                    if bundle_results._symbol_map is None:
2521                        bundle_results._symbol_map = \
2522                            symbol_map_from_instance(bundle_instance)
2523
2524                    bundle_instance = self._bundle_binding_instance_map[bundle_name]
2525
2526                    if self._verbose:
2527                        print("Results obtained for bundle=%s" % (bundle_name))
2528
2529                    if len(bundle_results.solution) == 0:
2530                        raise RuntimeError("Solve failed for scenario bundle=%s; no solutions "
2531                                           "generated\n%s" % (bundle_name, bundle_results))
2532
2533                    if self._output_solver_results:
2534                        print("Results for bundle=%s" % (bundle_name))
2535                        bundle_results.write(num=1)
2536
2537                    start_time = time.time()
2538                    bundle_instance.load(
2539                        bundle_results,
2540                        allow_consistent_values_for_fixed_vars=\
2541                            self._write_fixed_variables,
2542                        comparison_tolerance_for_fixed_vars=\
2543                            self._comparison_tolerance_for_fixed_vars)
2544
2545                    self._solver_results[bundle_name] = bundle_results
2546
2547                    solution0 = bundle_results.solution(0)
2548                    if hasattr(solution0, "gap") and \
2549                       (solution0.gap is not None):
2550                        self._gaps[bundle_name] = solution0.gap
2551
2552                    # if the solver plugin doesn't populate the
2553                    # user_time field, it is by default of type
2554                    # UndefinedData - defined in pyomo.opt.results
2555                    if hasattr(bundle_results.solver,"user_time") and \
2556                       (not isinstance(bundle_results.solver.user_time,
2557                                       UndefinedData)) and \
2558                       (bundle_results.solver.user_time is not None):
2559                        # the solve time might be a string, or might
2560                        # not be - we eventually would like more
2561                        # consistency on this front from the solver
2562                        # plugins.
2563                        self._solve_times[bundle_name] = \
2564                            float(bundle_results.solver.user_time)
2565
2566                    scenario_bundle = \
2567                        self._scenario_tree._scenario_bundle_map[bundle_name]
2568                    for scenario_name in scenario_bundle._scenario_names:
2569                        scenario = self._scenario_tree._scenario_map[scenario_name]
2570                        scenario.update_solution_from_instance()
2571
2572                    end_time = time.time()
2573                    if self._output_times:
2574                        print("Time loading results for bundle %s=%0.2f seconds"
2575                              % (bundle_name, end_time-start_time))
2576
2577                if self._verbose:
2578                    print("Successfully loaded solution for bundle=%s"
2579                          % (bundle_name))
2580
2581                num_results_so_far = num_results_so_far + 1
2582
2583        else:
2584
2585            if self._verbose:
2586                print("Waiting for scenario sub-problem solves")
2587
2588            num_results_so_far = 0
2589
2590            while (num_results_so_far < len(self._scenario_tree._scenarios)):
2591
2592                action_handle = self._solver_manager.wait_any()
2593                results = self._solver_manager.get_results(action_handle)
2594                scenario_name = action_handle_scenario_map[action_handle]
2595                scenario = self._scenario_tree._scenario_map[scenario_name]
2596
2597                if isinstance(self._solver_manager,
2598                              pyomo.solvers.plugins.smanager.\
2599                              phpyro.SolverManager_PHPyro):
2600
2601                    if self._output_solver_results:
2602                        print("Results for scenario= "+scenario_name)
2603
2604                    # TODO: Use these keywords to perform some
2605                    #       validation of fixed variable values in the
2606                    #       results returned
2607                    #allow_consistent_values_for_fixed_vars =\
2608                    #    self._write_fixed_variables,
2609                    #comparison_tolerance_for_fixed_vars =\
2610                    #    self._comparison_tolerance_for_fixed_vars
2611                    # results[0] are variable values
2612                    # results[1] are suffix values
2613                    # results[2] are auxilliary values
2614                    scenario.update_current_solution(results[0])
2615
2616                    auxilliary_values = results[2]
2617                    if "gap" in auxilliary_values:
2618                        self._gaps[scenario_name] = auxilliary_values["gap"]
2619
2620                    if "user_time" in auxilliary_values:
2621                        self._solve_times[scenario_name] = \
2622                            auxilliary_values["user_time"]
2623
2624                else:
2625
2626                    instance = scenario._instance
2627
2628                    # a temporary hack - if results come back from
2629                    # pyro, they won't have a symbol map attached. so
2630                    # create one.
2631                    if results._symbol_map is None:
2632                        results._symbol_map = symbol_map_from_instance(instance)
2633
2634                    if self._verbose:
2635                        print("Results obtained for scenario=%s" % (scenario_name))
2636
2637                    if len(results.solution) == 0:
2638                        raise RuntimeError("Solve failed for scenario=%s; no solutions "
2639                                           "generated\n%s" % (scenario_name, results))
2640
2641                    if self._output_solver_results:
2642                        print("Results for scenario="+scenario_name)
2643                        results.write(num=1)
2644
2645                    start_time = time.time()
2646
2647                    # TBD: Technically, we should validate that there
2648                    #      is only a single solution. Or at least warn
2649                    #      if there are multiple.
2650
2651                    instance.load(
2652                        results,
2653                        allow_consistent_values_for_fixed_vars=\
2654                            self._write_fixed_variables,
2655                        comparison_tolerance_for_fixed_vars=\
2656                            self._comparison_tolerance_for_fixed_vars)
2657
2658                    self._solver_results[scenario._name] = results
2659
2660                    scenario.update_solution_from_instance()
2661
2662                    solution0 = results.solution(0)
2663                    if hasattr(solution0, "gap") and \
2664                       (solution0.gap is not None):
2665                        self._gaps[scenario_name] = solution0.gap
2666
2667                    # if the solver plugin doesn't populate the
2668                    # user_time field, it is by default of type
2669                    # UndefinedData - defined in pyomo.opt.results
2670                    if hasattr(results.solver,"user_time") and \
2671                       (not isinstance(results.solver.user_time,
2672                                       UndefinedData)) and \
2673                       (results.solver.user_time is not None):
2674                        # the solve time might be a string, or might
2675                        # not be - we eventually would like more
2676                        # consistency on this front from the solver
2677                        # plugins.
2678                        self._solve_times[scenario_name] = \
2679                            float(results.solver.user_time)
2680
2681                    end_time = time.time()
2682
2683                    if self._output_times:
2684                        print("Time loading results into instance %s=%0.2f seconds"
2685                              % (scenario_name, end_time-start_time))
2686
2687                if self._verbose:
2688                    print("Successfully loaded solution for scenario=%s"
2689                          % (scenario_name))
2690
2691                num_results_so_far = num_results_so_far + 1
2692
2693        if len(self._solve_times) > 0:
2694            # if any of the solve times are of type
2695            # pyomo.opt.results.container.UndefinedData, then don't
2696            # output timing statistics.
2697            undefined_detected = False
2698            for this_time in itervalues(self._solve_times):
2699                if isinstance(this_time, UndefinedData):
2700                    undefined_detected=True
2701            if undefined_detected:
2702                print("At least one sub-problem solve time was "
2703                      "undefined - skipping timing statistics")
2704            else:
2705                mean = sum(self._solve_times.values()) / \
2706                        float(len(self._solve_times.values()))
2707                std_dev = sqrt(
2708                    sum(pow(x-mean,2.0) for x in self._solve_times.values()) /
2709                    float(len(self._solve_times.values())))
2710                if self._output_times:
2711                    print("Sub-problem solve time statistics - Min: "
2712                          "%0.2f Avg: %0.2f Max: %0.2f StdDev: %0.2f (seconds)"
2713                          % (min(self._solve_times.values()),
2714                             mean,
2715                             max(self._solve_times.values()),
2716                             std_dev))
2717
2718#                print "**** SOLVE TIMES:",self._solve_times.values()
2719#                print "*** GAPS:",sorted(self._gaps.values())
2720
2721        iteration_end_time = time.time()
2722        self._cumulative_solve_time += (iteration_end_time - iteration_start_time)
2723
2724        if self._output_times:
2725            print("Sub-problem solve time=%.2f seconds"
2726                  % (iteration_end_time - iteration_start_time))
2727
2728
2729    """ Perform the non-weighted scenario solves and form the initial w and xbars.
2730    """
2731    def iteration_0_solves(self):
2732
2733        # return None unless a sub-problem failure is detected, then
2734        # return its name immediately
2735
2736        if self._verbose:
2737            print("------------------------------------------------")
2738            print("Starting PH iteration 0 solves")
2739
2740        # scan any variables fixed prior to iteration 0, set up the
2741        # appropriate flags for pre-processing, and - if appropriate -
2742        # transmit the information to the PH solver servers.
2743        self._push_fix_queue_to_instances()
2744
2745        self.solve_subproblems(warmstart=False)
2746
2747        if self._verbose:
2748            print("Successfully completed PH iteration 0 solves\n"
2749                  "- solution statistics:\n")
2750            if self._scenario_tree.contains_bundles():
2751                self._report_bundle_objectives()
2752            self._report_scenario_objectives()
2753
2754    #
2755    # recompute the averages, minimum, and maximum statistics for all
2756    # variables to be blended by PH, i.e., not appearing in the final
2757    # stage. technically speaking, the min/max aren't required by PH,
2758    # but they are used often enough to warrant their computation and
2759    # it's basically free if you're computing the average.
2760    #
2761    # **When compute_xbars is False, the xbar is not assigned the
2762    #   calculated node averages. The instance parameters are still
2763    #   updated to the current value of xbar as usual. The dual ph
2764    #   algorithm uses both versions of this method
2765    #
2766
2767    def update_variable_statistics(self, compute_xbars=True):
2768
2769        start_time = time.time()
2770        # cache the lookups - don't want to do them deep in the index loop.
2771        overrelax = self._overrelax
2772        current_iteration = self._current_iteration
2773
2774        # skip the last stage, as there is only a single scenario there - no
2775        # meaningful statistics required!
2776        for stage in self._scenario_tree._stages[:-1]:
2777
2778            for tree_node in stage._tree_nodes:
2779
2780                xbars = tree_node._xbars
2781
2782                scenario_solutions = \
2783                    [(scenario._probability, scenario._x[tree_node._name]) \
2784                     for scenario in tree_node._scenarios]
2785
2786                for variable_id in tree_node._standard_variable_ids:
2787
2788                    stale = False
2789                    values = []
2790                    avg_value = 0.0
2791                    for probability, var_values in scenario_solutions:
2792                        val = var_values[variable_id]
2793                        if val is not None:
2794                            avg_value += probability * val
2795                            values.append(val)
2796                        else:
2797                            stale = True
2798
2799                    if not stale:
2800
2801                        avg_value /= tree_node._probability
2802                        tree_node._minimums[variable_id] = min(values)
2803                        tree_node._maximums[variable_id] = max(values)
2804                        tree_node._averages[variable_id] = avg_value
2805
2806                        if compute_xbars:
2807                            if (overrelax) and (current_iteration >= 1):
2808                                xbars[variable_id] = self._nu*avg_value + (1-self._nu)*var_xbar
2809                            else:
2810                                xbars[variable_id] = avg_value
2811
2812        end_time = time.time()
2813        self._cumulative_xbar_time += (end_time - start_time)
2814
2815        if self._output_times:
2816            print("Variable statistics compute time=%.2f seconds" % (end_time - start_time))
2817
2818    def update_weights(self):
2819
2820        start_time = time.time()
2821
2822        # because the weight updates rely on the xbars, and the xbars
2823        # are node-based, I'm looping over the tree nodes and pushing
2824        # weights into the corresponding scenarios.
2825        start_time = time.time()
2826
2827        # NOTE: the following code has some optimizations that are not
2828        #       normally recommended, in particular the direct access
2829        #       and manipulation of parameters via the .value
2830        #       attribute instead of the user-level-preferred value()
2831        #       method. this is justifiable in this particular
2832        #       instance because we are creating the PH parameters
2833        #       (and therefore can manipulate them safely), and this
2834        #       routine takes a non-trivial amount of the overall
2835        #       run-time.
2836
2837        # cache the lookups - don't want to do them deep in the index
2838        # loop.
2839        over_relaxing = self._overrelax
2840        objective_sense = self._objective_sense
2841
2842        # no blending over the final stage, so no weights to worry
2843        # about.
2844        for stage in self._scenario_tree._stages[:-1]:
2845
2846            for tree_node in stage._tree_nodes:
2847
2848                tree_node_xbars = None
2849                if self._dual_mode is True:
2850                    tree_node_xbars = tree_node._xbars
2851                else:
2852                    tree_node_xbars = tree_node._averages
2853                blend_values = tree_node._blend
2854
2855                # These will be updated inside this loop
2856                tree_node_wbars = tree_node._wbars = \
2857                    dict((var_id,0) for var_id in tree_node._variable_ids)
2858
2859                for scenario in tree_node._scenarios:
2860
2861                    instance = scenario._instance
2862
2863                    weight_values = scenario._w[tree_node._name]
2864                    rho_values = scenario._rho[tree_node._name]
2865                    var_values = scenario._x[tree_node._name]
2866
2867                    for variable_id in tree_node._standard_variable_ids:
2868
2869                        varval = var_values[variable_id]
2870
2871                        if varval is not None:
2872
2873                            # we are currently not updating weights if
2874                            # blending is disabled for a variable.
2875                            # this is done on the premise that unless
2876                            # you are actively trying to move the
2877                            # variable toward the mean, the weights
2878                            # will blow up and be huge by the time
2879                            # that blending is activated.
2880
2881                            nu_value = 1.0
2882                            if over_relaxing:
2883                                nu_value = self._nu
2884
2885                            if not self._dual_mode:
2886
2887                                if objective_sense == minimize:
2888                                    weight_values[variable_id] += \
2889                                        blend_values[variable_id] * \
2890                                        rho_values[variable_id] * \
2891                                        nu_value * \
2892                                        (varval - \
2893                                         tree_node_xbars[variable_id])
2894                                else:
2895                                    weight_values[variable_id] -= \
2896                                        blend_values[variable_id] * \
2897                                        rho_values[variable_id] * \
2898                                        nu_value * \
2899                                        (varval - \
2900                                         tree_node_xbars[variable_id])
2901                            else:
2902                                # **Adding these asserts simply
2903                                # **because we haven't thought about
2904                                # **what this means for other steps in
2905                                # **the code
2906                                assert blend_values[variable_id] == 1.0
2907                                assert nu_value == 1.0
2908                                assert objective_sense == minimize
2909                                weight_values[variable_id] = \
2910                                    blend_values[variable_id] * \
2911                                    (rho_values[variable_id]) * \
2912                                    nu_value * \
2913                                    (varval - \
2914                                     tree_node_xbars[variable_id])
2915
2916                            tree_node_wbars[variable_id] += \
2917                                scenario._probability * \
2918                                weight_values[variable_id] / tree_node._probability
2919
2920        end_time = time.time()
2921        self._cumulative_weight_time += (end_time - start_time)
2922
2923        if self._output_times:
2924            print("Weight update time=%.2f seconds" % (end_time - start_time))
2925
2926    def update_weights_for_scenario(self, scenario):
2927
2928        start_time = time.time()
2929
2930        # cache the lookups - don't want to do them deep in the index
2931        # loop.
2932        over_relaxing = self._overrelax
2933        objective_sense = self._objective_sense
2934
2935        nodeid_to_var_map = scenario._instance._ScenarioTreeSymbolMap.bySymbol
2936
2937        for tree_node in scenario._node_list[:-1]:
2938
2939            weight_values = scenario._w[tree_node._name]
2940            rho_values = scenario._rho[tree_node._name]
2941            var_values = scenario._x[tree_node._name]
2942
2943            tree_node_xbars = None
2944            if self._dual_mode is True:
2945                tree_node_xbars = tree_node._xbars
2946            else:
2947                tree_node_xbars = tree_node._averages
2948            blend_values = tree_node._blend
2949
2950            # Note: This does not update wbar
2951            for variable_id in tree_node._standard_variable_ids:
2952
2953                varval = var_values[variable_id]
2954
2955                if varval is not None:
2956
2957                    # we are currently not updating weights if
2958                    # blending is disabled for a variable.  this is
2959                    # done on the premise that unless you are actively
2960                    # trying to move the variable toward the mean, the
2961                    # weights will blow up and be huge by the time
2962                    # that blending is activated.
2963
2964                    nu_value = 1.0
2965                    if over_relaxing:
2966                        nu_value = self._nu
2967
2968                    if self._dual_mode is False:
2969                        if objective_sense == minimize:
2970                            weight_values[variable_id] += \
2971                                blend_values[variable_id] * \
2972                                rho_values[variable_id] * \
2973                                nu_value * \
2974                                (varval - \
2975                                 tree_node_xbars[variable_id])
2976                        else:
2977                            weight_values[variable_id] -= \
2978                                blend_values[variable_id] * \
2979                                rho_values[variable_id] * \
2980                                nu_value * \
2981                                (varval - \
2982                                 tree_node_xbars[variable_id])
2983                    else:
2984                        # **Adding these asserts simply because we
2985                        # **haven't thought about what this means for
2986                        # **other steps in the code
2987                        assert blend_values[variable_id] == 1.0
2988                        assert nu_value == 1.0
2989                        assert objective_sense == minimize
2990                        weight_values[variable_id] = \
2991                            blend_values[variable_id] * \
2992                            rho_values[variable_id] * \
2993                            nu_value * \
2994                            (varval - \
2995                             tree_node_xbars[variable_id])
2996
2997        end_time = time.time()
2998        self._cumulative_weight_time += (end_time - start_time)
2999
3000    def iteration_k_solves(self):
3001
3002        if self._verbose:
3003            print("------------------------------------------------")
3004            print("Starting PH iteration %s solves"
3005                  % (self._current_iteration))
3006
3007        # scan any variables fixed/freed, set up the appropriate flags
3008        # for pre-processing, and - if appropriate - transmit the
3009        # information to the PH solver servers.
3010        self._push_fix_queue_to_instances()
3011
3012        # update parameters on instances (transmitting to ph solver
3013        # servers when appropriate)
3014        self._push_xbar_to_instances()
3015        self._push_w_to_instances()
3016        # NOTE: We now transmit rho at every iteration as
3017        #       it can be adjusted for better convergence
3018        self._push_rho_to_instances()
3019
3020        # STEP -1: if using a PH solver manager, propagate current
3021        #          weights/averages to the appropriate solver servers.
3022        #          ditto the tree node statistics, which are necessary
3023        #          if linearizing (so an optimization could be
3024        #          performed here).
3025        if isinstance(self._solver_manager,
3026                      pyomo.solvers.plugins.smanager.phpyro.SolverManager_PHPyro):
3027
3028            # we only transmit tree node statistics if we are
3029            # linearizing the PH objectives.  otherwise, the
3030            # statistics are never referenced by the PH solver
3031            # servers, so don't want the time.
3032            if (self._linearize_nonbinary_penalty_terms > 0) or \
3033               (self._scenario_tree.contains_bundles()):
3034                transmit_tree_node_statistics(self)
3035
3036        else:
3037
3038            # GAH: We may need to redefine our concept of
3039            #      warmstart. These values could be helpful in the
3040            #      nonlinear case (or could be better than 0.0, the
3041            #      likely default used by the solver when these
3042            #      initializations are not placed in the NL
3043            #      file. **Note: Initializations go into the NL file
3044            #      independent of the "warmstart" keyword
3045
3046            # STEP 0.85: if linearizing the PH objective, clear the
3047            #            values for any PHQUADPENALTY* variables -
3048            #            otherwise, the MIP starts are likely to be
3049            #            infeasible.
3050            if self._linearize_nonbinary_penalty_terms > 0:
3051                self._reset_instance_linearization_variables()
3052
3053            if self._scenario_tree.contains_bundles():
3054                # clear non-converged variables and stage cost
3055                # variables, to ensure feasible warm starts.
3056                reset_nonconverged_variables(self._scenario_tree, self._instances)
3057                reset_stage_cost_variables(self._scenario_tree, self._instances)
3058            else:
3059                # clear stage cost variables, to ensure feasible warm starts.
3060                reset_stage_cost_variables(self._scenario_tree, self._instances)
3061
3062        self.solve_subproblems(warmstart=not self._disable_warmstarts)
3063
3064        if self._verbose:
3065            print("Successfully completed PH iteration %s solves\n"
3066                  "- solution statistics:\n" % (self._current_iteration))
3067            if self._scenario_tree.contains_bundles():
3068                self._report_bundle_objectives()
3069            self._report_scenario_objectives()
3070
3071    def async_iteration_k_plus_solves(self):
3072        # note: this routine retains control until a termination
3073        # criterion is met modified nov 2011 by dlw to do async with a
3074        # window-like paramater
3075
3076        if self._scenario_tree.contains_bundles():
3077            raise RuntimeError("Async PH does not currently support bundling")
3078
3079        if (self._async_buffer_len <= 0) or \
3080           (self._async_buffer_len > len(self._scenario_tree._scenarios)):
3081            raise RuntimeError("Async buffer length parameter is bad: %s"
3082                               % (self._async_buffer_len))
3083        if self._verbose:
3084            print("Starting PH iteration k+ solves - running async "
3085                  "with buffer length=%s" % (self._async_buffer_len))
3086
3087        # we are going to buffer the scenario names
3088        ScenarioBuffer = []
3089
3090        # things progress at different rates - keep track of what's
3091        # going on.
3092        total_scenario_solves = 0
3093        # a map of scenario name to the number of sub-problems solved
3094        # thus far.
3095        scenario_ks = {}
3096        for scenario in self._scenario_tree._scenarios:
3097            scenario_ks[scenario._name] = 0
3098
3099        # keep track of action handles mapping to scenarios.
3100        action_handle_instance_map = {}
3101
3102        # only form the PH objective components once, before we start
3103        # in on the asychronous sub-problem solves
3104        self.activate_ph_objective_weight_terms()
3105        self.activate_ph_objective_proximal_terms()
3106
3107        # if linearizing, form the necessary terms to compute the cost
3108        # variables.
3109        if self._linearize_nonbinary_penalty_terms > 0:
3110            self.form_ph_linearized_objective_constraints()
3111
3112        # scan any variables fixed/freed, set up the appropriate flags
3113        # for pre-processing, and - if appropriate - transmit the
3114        # information to the PH solver servers.
3115        self._push_fix_queue_to_instances()
3116
3117        # update parameters on instances (transmitting to ph solver
3118        # servers when appropriate)
3119        self._push_xbar_to_instances()
3120        self._push_w_to_instances()
3121        # NOTE: We aren't currently propagating rhos, as they
3122        #       generally don't change - we need to have a flag,
3123        #       though, indicating whether the rhos have changed, so
3124        #       they can be transmitted if needed.
3125        self._push_rho_to_instances()
3126
3127        self._preprocess_scenario_instances()
3128
3129        if self._verbose or self._report_rhos_first_iteration or self._report_rhos_each_iteration:
3130            print("Async starting rhos:")
3131            self.pprint(False, False, False, False, True,
3132                        output_only_statistics=\
3133                        self._report_only_statistic,
3134                        output_only_nonconverged=\
3135                        self._report_only_nonconverged_variables,
3136                        report_stage_costs=False)
3137
3138        # STEP 1: queue up the solves for all scenario sub-problems.
3139
3140        for scenario in self._scenario_tree._scenarios:
3141
3142            instance = self._instances[scenario._name]
3143
3144            if self._verbose == True:
3145                print("Queuing solve for scenario=%s" % (scenario._name))
3146
3147            # once past iteration 0, there is always a feasible
3148            # solution from which to warm-start.
3149            if (not self._disable_warmstarts) and \
3150               (self._solver.warm_start_capable()):
3151                new_action_handle = \
3152                    self._solver_manager.queue(instance,
3153                                               opt=self._solver,
3154                                               warmstart=True,
3155                                               tee=self._output_solver_logs,
3156                                               verbose=self._verbose)
3157            else:
3158                new_action_handle = \
3159                    self._solver_manager.queue(instance,
3160                                               opt=self._solver,
3161                                               tee=self._output_solver_logs,
3162                                               verbose=self._verbose)
3163
3164            action_handle_instance_map[new_action_handle] = scenario._name
3165
3166        # STEP 2: wait for the first action handle to return, process
3167        #         it, and keep chugging.
3168
3169        while(True):
3170
3171            this_action_handle = self._solver_manager.wait_any()
3172            solved_scenario_name = action_handle_instance_map[this_action_handle]
3173            solved_scenario = self._scenario_tree.get_scenario(solved_scenario_name)
3174
3175            scenario_ks[solved_scenario_name] += 1
3176            total_scenario_solves += 1
3177
3178            if int(total_scenario_solves / len(scenario_ks)) > \
3179               self._current_iteration:
3180                new_reportable_iteration = True
3181                self._current_iteration += 1
3182            else:
3183                new_reportable_iteration = False
3184
3185            if self._verbose:
3186                print("Solve for scenario=%s completed - new solve count for "
3187                      "this scenario=%s"
3188                      % (solved_scenario_name,
3189                         scenario_ks[solved_scenario_name]))
3190
3191            instance = self._instances[solved_scenario_name]
3192            results = self._solver_manager.get_results(this_action_handle)
3193
3194            if len(results.solution) == 0:
3195                raise RuntimeError("Solve failed for scenario=%s; no solutions "
3196                                   "generated" % (solved_scenario_name))
3197
3198            if self._verbose:
3199                print("Solve completed successfully")
3200
3201            if self._output_solver_results == True:
3202                print("Results:")
3203                results.write(num=1)
3204
3205            # in async mode, it is possible that we will receive
3206            # values for variables that have been fixed due to
3207            # apparent convergence - but the outstanding scenario
3208            # solves will obviously not know this. if the values are
3209            # inconsistent, we have bigger problems - an exception
3210            # will be thrown, and we currently lack a recourse
3211            # mechanism in such a case.
3212            instance.load(results,
3213                          allow_consistent_values_for_fixed_vars=\
3214                              self._write_fixed_variables,
3215                          comparison_tolerance_for_fixed_vars=\
3216                              self._comparison_tolerance_for_fixed_vars)
3217
3218            self._solver_results[solved_scenario._name] = results
3219
3220            solved_scenario.update_solution_from_instance()
3221
3222            if self._verbose:
3223                print("Successfully loaded solution")
3224
3225            if self._verbose:
3226                print("%20s       %18.4f     %14.4f"
3227                      % (solved_scenario_name,
3228                         solved_scenario._objective,
3229                         0.0))
3230
3231            # changed 19 Nov 2011 to support scenario buffers for async
3232            ScenarioBuffer.append(solved_scenario_name)
3233            if len(ScenarioBuffer) == self._async_buffer_len:
3234                if self._verbose:
3235                    print("Processing async buffer.")
3236
3237                # update variable statistics prior to any output.
3238                self.update_variable_statistics()
3239                for scenario_name in ScenarioBuffer:
3240                    scenario = self._scenario_tree.get_scenario(scenario_name)
3241                    self.update_weights_for_scenario(scenario)
3242                    scenario.push_w_to_instance()
3243
3244                # we don't want to report stuff and invoke callbacks
3245                # after each scenario solve - wait for when each
3246                # scenario (on average) has reported back a solution.
3247                if new_reportable_iteration:
3248
3249                    # let plugins know if they care.
3250                    for plugin in self._ph_plugins:
3251                        plugin.post_iteration_k_solves(self)
3252
3253                    # update the fixed variable statistics.
3254                    self._total_fixed_discrete_vars,\
3255                        self._total_fixed_continuous_vars = \
3256                            self.compute_fixed_variable_counts()
3257
3258                    if self._report_rhos_each_iteration:
3259                        print("Async Reportable Iteration Current rhos:")
3260                        self.pprint(False, False, False, False, True,
3261                                    output_only_statistics=\
3262                                    self._report_only_statistic,
3263                                    output_only_nonconverged=\
3264                                    self._report_only_nonconverged_variables,
3265                                    report_stage_costs=False)
3266
3267                    if self._verbose or self._report_weights:
3268                        print("Async Reportable Iteration Current variable "
3269                              "averages and weights:")
3270                        self.pprint(True, True, False, False, False,
3271                                    output_only_statistics=\
3272                                    self._report_only_statistics,
3273                                    output_only_nonconverged=\
3274                                    self._report_only_nonconverged_variables)
3275
3276                    # check for early termination.
3277                    self._converger.update(self._current_iteration,
3278                                           self,
3279                                           self._scenario_tree,
3280                                           self._instances)
3281                    first_stage_min, first_stage_avg, first_stage_max = \
3282                        self._extract_first_stage_cost_statistics()
3283                    print("Convergence metric=%12.4f  First stage cost avg=%12.4f "
3284                          "Max-Min=%8.2f" % (self._converger.lastMetric(),
3285                                             first_stage_avg,
3286                                             first_stage_max-first_stage_min))
3287
3288                    expected_cost = self._scenario_tree.findRootNode().computeExpectedNodeCost()
3289                    if not _OLD_OUTPUT: print("Expected Cost=%14.4f" % (expected_cost))
3290                    self._cost_history[self._current_iteration] = expected_cost
3291                    if self._converger.isConverged(self):
3292
3293                        if self._total_discrete_vars == 0:
3294                            print("PH converged - convergence metric is below "
3295                                  "threshold="
3296                                  +str(self._converger._convergence_threshold))
3297                        else:
3298                            print("PH converged - convergence metric is below "
3299                                  "threshold="
3300                                  +str(self._converger._convergence_threshold)+
3301                                  " or all discrete variables are fixed")
3302
3303                        if (len(self._incumbent_cost_history) == 0) or \
3304                           ((self._objective_sense == minimize) and \
3305                            (expected_cost < min(self._incumbent_cost_history))) or \
3306                           ((self._objective_sense == maximize) and \
3307                            (expected_cost > max(self._incumbent_cost_history))):
3308                            if not _OLD_OUTPUT: print("Caching results for new incumbent solution")
3309                            self.cacheSolutions(self._incumbent_cache_id)
3310                            self._best_incumbent_key = self._current_iteration
3311                        self._incumbent_cost_history[self._current_iteration] = expected_cost
3312
3313                        plugin_convergence = True
3314                        for plugin in self._ph_plugins:
3315                            if hasattr(plugin,"ph_convergence_check"):
3316                                if not plugin.ph_convergence_check(self):
3317                                    plugin_convergence = False
3318
3319                        if plugin_convergence:
3320                            break
3321
3322                # see if we've exceeded our patience with the
3323                # iteration limit.  changed to be based on the average
3324                # on July 10, 2011 by dlw (really, it should be some
3325                # combination of the min and average over the
3326                # scenarios)
3327                if total_scenario_solves / len(self._scenario_tree._scenarios) >= \
3328                   self._max_iterations:
3329                    return
3330
3331                # update parameters on instances
3332                self._push_xbar_to_instances()
3333                self._push_w_to_instances()
3334
3335                # we're still good to run - re-queue the instance,
3336                # following any necessary linearization
3337                for  scenario_name in ScenarioBuffer:
3338                    instance = self._instances[scenario_name]
3339
3340                    # if linearizing, form the necessary terms to
3341                    # compute the cost variables.
3342                    if self._linearize_nonbinary_penalty_terms > 0:
3343                        new_attrs = \
3344                            form_linearized_objective_constraints(
3345                                scenario_name,
3346                                instance,
3347                                self._scenario_tree,
3348                                self._linearize_nonbinary_penalty_terms,
3349                                self._breakpoint_strategy,
3350                                self._integer_tolerance)
3351                        self._problem_states.ph_constraints[scenario_name].\
3352                            extend(new_attrs)
3353                        # Flag the preprocessor
3354                        self._problem_states.\
3355                            ph_constraints_updated[scenario_name] = True
3356
3357                    # preprocess prior to queuing the solve.
3358                    preprocess_scenario_instance(
3359                        instance,
3360                        self._problem_states.fixed_variables[scenario_name],
3361                        self._problem_states.freed_variables[scenario_name],
3362                        self._problem_states.\
3363                            user_constraints_updated[scenario_name],
3364                        self._problem_states.ph_constraints_updated[scenario_name],
3365                        self._problem_states.ph_constraints[scenario_name],
3366                        self._problem_states.objective_updated[scenario_name],
3367                        not self._write_fixed_variables,
3368                        self._solver)
3369
3370                    # We've preprocessed the instance, reset the
3371                    # relevant flags
3372                    self._problem_states.clear_update_flags(scenario_name)
3373                    self._problem_states.clear_fixed_variables(scenario_name)
3374                    self._problem_states.clear_freed_variables(scenario_name)
3375
3376                    # once past the initial sub-problem solves, there
3377                    # is always a feasible solution from which to
3378                    # warm-start.
3379                    if (not self._disable_warmstarts) and \
3380                       self._solver.warm_start_capable():
3381                        new_action_handle = \
3382                            self._solver_manager.queue(
3383                                instance,
3384                                opt=self._solver,
3385                                warmstart=True,
3386                                tee=self._output_solver_logs,
3387                                verbose=self._verbose)
3388                    else:
3389                        new_action_handle = \
3390                            self._solver_manager.queue(
3391                                instance,
3392                                opt=self._solver,
3393                                tee=self._output_solver_logs,
3394                                verbose=self._verbose)
3395
3396                    action_handle_instance_map[new_action_handle] = scenario_name
3397
3398                    if self._verbose:
3399                        print("Queued solve k=%s for scenario=%s"
3400                              % (scenario_ks[scenario_name]+1,
3401                                 solved_scenario_name))
3402
3403                    if self._verbose:
3404                        for sname, scenario_count in iteritems(scenario_ks):
3405                            print("Scenario=%s was solved %s times"
3406                                  % (sname, scenario_count))
3407                        print("Cumulative number of scenario solves=%s"
3408                              % (total_scenario_solves))
3409                        print("PH Iteration Count (computed)=%s"
3410                              % (self._current_iteration))
3411
3412                    if self._verbose:
3413                        print("Variable values following scenario solves:")
3414                        self.pprint(False, False, True, False, False,
3415                                    output_only_statistics=\
3416                                    self._report_only_statistics,
3417                                    output_only_nonconverged=\
3418                                        self._report_only_nonconverged_variables)
3419
3420                if self._verbose is True:
3421                    print("Emptying the asynch scenario buffer.")
3422                # this is not a speed issue, is there a memory issue?
3423                ScenarioBuffer = []
3424
3425    def solve(self):
3426        # return None unless a solve failure was detected in iter0,
3427        # then immediately return the iter0 solve return value (which
3428        # should be the name of the scenario detected)
3429
3430        self._solve_start_time = time.time()
3431        self._cumulative_solve_time = 0.0
3432        self._cumulative_xbar_time = 0.0
3433        self._cumulative_weight_time = 0.0
3434        self._current_iteration = 0;
3435
3436        print("Starting PH")
3437
3438        if self._initialized == False:
3439            raise RuntimeError("PH is not initialized - cannot invoke solve() method")
3440
3441        # garbage collection noticeably slows down PH when dealing with
3442        # large numbers of scenarios. fortunately, there are well-defined
3443        # points at which garbage collection makes sense (and there isn't a
3444        # lot of collection to do). namely, after each PH iteration.
3445        re_enable_gc = gc.isenabled()
3446        gc.disable()
3447
3448        print("")
3449        print("Initiating PH iteration=" + str(self._current_iteration))
3450
3451        iter0retval = self.iteration_0_solves()
3452
3453        if iter0retval is not None:
3454            if self._verbose:
3455                print("Iteration zero reports trouble with scenario: "
3456                      +str(iter0retval))
3457            return iter0retval
3458
3459        # now that we have scenario solutions, compute and cache the
3460        # number of discrete and continuous variables.  the values are
3461        # of general use, e.g., in the converger classes and in
3462        # plugins. this is only invoked once, after the iteration 0
3463        # solves.
3464        (self._total_discrete_vars,self._total_continuous_vars) = \
3465            self.compute_blended_variable_counts()
3466
3467        if self._verbose:
3468            print("Total number of non-stale discrete instance variables="
3469                  +str(self._total_discrete_vars))
3470            print("Total number of non-stale continuous instance variables="
3471                  +str(self._total_continuous_vars))
3472
3473        # very rare, but the following condition can actually happen...
3474        if (self._total_discrete_vars + self._total_continuous_vars) == 0:
3475            raise RuntimeError("***ERROR: The total number of non-anticipative "
3476                               "discrete and continuous variables equals 0! "
3477                               "Did you set the StageVariables set(s) in "
3478                               "ScenarioStructure.dat")
3479
3480        # update variable statistics prior to any output, and most
3481        # importantly, prior to any variable fixing by PH extensions.
3482        self.update_variable_statistics()
3483
3484        if (self._verbose) or (self._report_solutions):
3485            print("Variable values following scenario solves:")
3486            self.pprint(False, False, True, False, False,
3487                        output_only_statistics=\
3488                            self._report_only_statistics,
3489                        output_only_nonconverged=\
3490                            self._report_only_nonconverged_variables)
3491
3492        # let plugins know if they care.
3493        for plugin in self._ph_plugins:
3494            plugin.post_iteration_0_solves(self)
3495
3496        # update the fixed variable statistics.
3497        self._total_fixed_discrete_vars, \
3498            self._total_fixed_continuous_vars = \
3499                self.compute_fixed_variable_counts()
3500
3501        print("Number of discrete variables fixed="
3502              +str(self._total_fixed_discrete_vars)+
3503              " (total="+str(self._total_discrete_vars)+")")
3504        print("Number of continuous variables fixed="
3505              +str(self._total_fixed_continuous_vars)+
3506              " (total="+str(self._total_continuous_vars)+")")
3507
3508        # always output the convergence metric and first-stage cost
3509        # statistics, to give a sense of progress.
3510        self._converger.update(self._current_iteration,
3511                               self,
3512                               self._scenario_tree,
3513                               self._instances)
3514        first_stage_min, first_stage_avg, first_stage_max = \
3515            self._extract_first_stage_cost_statistics()
3516        print("Convergence metric=%12.4f  "
3517              "First stage cost avg=%12.4f  "
3518              "Max-Min=%8.2f"
3519              % (self._converger.lastMetric(),
3520                 first_stage_avg,
3521                 first_stage_max-first_stage_min))
3522
3523        expected_cost = self._scenario_tree.findRootNode().computeExpectedNodeCost()
3524        if not _OLD_OUTPUT: print("Expected Cost=%14.4f" % (expected_cost))
3525        self._cost_history[self._current_iteration] = expected_cost
3526        if self._converger.isConverged(self):
3527            if not _OLD_OUTPUT: print("Caching results for new incumbent solution")
3528            self.cacheSolutions(self._incumbent_cache_id)
3529            self._best_incumbent_key = self._current_iteration
3530            self._incumbent_cost_history[self._current_iteration] = expected_cost
3531
3532        # let plugins know if they care.
3533        for plugin in self._ph_plugins:
3534            plugin.post_iteration_0(self)
3535
3536        # IMPT: update the weights after the PH iteration 0 callbacks;
3537        #       they might compute rhos based on iteration 0
3538        #       solutions.
3539        if self._ph_weight_updates_enabled:
3540            self.update_weights()
3541
3542        # checkpoint if it's time - which it always is after iteration 0,
3543        # if the interval is >= 1!
3544        if (self._checkpoint_interval > 0):
3545            self.checkpoint(0)
3546
3547        # garbage-collect if it wasn't disabled entirely.
3548        if re_enable_gc:
3549            if (time.time() - self._time_since_last_garbage_collect) >= self._minimum_garbage_collection_interval:
3550               gc.collect()
3551               self._time_last_garbage_collect = time.time()
3552
3553        # everybody wants to know how long they've been waiting...
3554        print("Cumulative run-time=%.2f seconds" % (time.time() - self._solve_start_time))
3555
3556        # gather memory statistics (for leak detection purposes) if specified.
3557        # XXX begin debugging - commented
3558        #if (pympler_available) and (self._profile_memory >= 1):
3559        #    objects_last_iteration = muppy.get_objects()
3560        #    summary_last_iteration = summary.summarize(objects_last_iteration)
3561        # XXX end debugging - commented
3562
3563        ####################################################################################################
3564        # major logic branch - if we are not running async, do the usual PH - otherwise, invoke the async. #
3565        ####################################################################################################
3566        if self._async is False:
3567
3568            # Note: As a first pass at our implementation, the solve method
3569            #       on the DualPHModel actually updates the xbar dictionary
3570            #       on the ph scenario tree.
3571            if (self._dual_mode is True):
3572                WARM_START = True
3573                dual_model = DualPHModel(self)
3574                print("Warm-starting dual-ph weights")
3575                if WARM_START is True:
3576                    self.activate_ph_objective_proximal_terms()
3577
3578                    self.iteration_k_solves()
3579                    dual_model.add_cut(first=True)
3580                    dual_model.solve()
3581                    self.update_variable_statistics(compute_xbars=False)
3582                    self.update_weights()
3583
3584                    self.deactivate_ph_objective_weight_terms()
3585                    self.activate_ph_objective_proximal_terms()
3586                else:
3587                    dual_model.add_cut()
3588                    dual_model.solve()
3589                    self.update_variable_statistics(compute_xbars=False)
3590                    self.update_weights()
3591
3592                    self.activate_ph_objective_proximal_terms()
3593            else:
3594                self.activate_ph_objective_weight_terms()
3595                self.activate_ph_objective_proximal_terms()
3596
3597
3598            ####################################################################################################
3599
3600            # there is an upper bound on the number of iterations to execute -
3601            # the actual bound depends on the converger supplied by the user.
3602            for i in xrange(1, self._max_iterations+1):
3603
3604                # XXX begin debugging
3605                #def muppetize(self):
3606                #    if (pympler_available) and (self._profile_memory >= 1):
3607                #        objects_this_iteration = muppy.get_objects()
3608                #        summary_this_iteration = summary.summarize(objects_this_iteration)
3609                #        summary.print_(summary_this_iteration)
3610                #        del summary_this_iteration
3611                # XXX end debugging
3612
3613                self._current_iteration = self._current_iteration + 1
3614
3615                print("")
3616                print("Initiating PH iteration=" + str(self._current_iteration))
3617
3618                # let plugins know if they care.
3619                for plugin in self._ph_plugins:
3620                    plugin.pre_iteration_k_solves(self)
3621
3622                if not _OLD_OUTPUT:
3623                    if self._report_rhos_each_iteration or \
3624                       ((self._verbose or self._report_rhos_first_iteration) and \
3625                        (self._current_iteration == 1)):
3626                        print("Rhos prior to scenario solves:")
3627                        self.pprint(False, False, False, False, True,
3628                                    output_only_statistics=self._report_only_statistics,
3629                                    output_only_nonconverged=self._report_only_nonconverged_variables,
3630                                    report_stage_costs=False)
3631
3632                if (self._verbose) or (self._report_weights):
3633                    print("Variable averages and weights prior to scenario solves:")
3634                    self.pprint(True, True, False, False, False,
3635                                output_only_statistics=self._report_only_statistics,
3636                                output_only_nonconverged=self._report_only_nonconverged_variables)
3637
3638                # with the introduction of piecewise linearization, the form of the
3639                # penalty-weighted objective is no longer fixed. thus, when linearizing,
3640                # we need to construct (or at least modify) the constraints used to
3641                # compute the linearized cost terms.
3642                if self._linearize_nonbinary_penalty_terms > 0:
3643                    self.form_ph_linearized_objective_constraints()
3644
3645                try:
3646                    self.iteration_k_solves()
3647                except SystemExit:
3648                    print("")
3649                    print(" ** Caught SystemExit exception. "
3650                          "Attempting to gracefully exit PH")
3651                    print(" Signal: "+str(sys.exc_info()[1]))
3652                    print("")
3653                    self._current_iteration -= 1
3654                    break
3655                except pyutilib.common._exceptions.ApplicationError:
3656                    print("")
3657                    print(" ** Caught ApplicationError exception. "
3658                          "Attempting to gracefully exit PH")
3659                    print(" Signal: "+str(sys.exc_info()[1]))
3660                    print("")
3661                    self._current_iteration -= 1
3662                    break
3663
3664                # update variable statistics prior to any output.
3665                if self._dual_mode is False:
3666                    self.update_variable_statistics()
3667                else:
3668                    dual_rc = dual_model.add_cut()
3669                    dual_model.solve()
3670                    self.update_variable_statistics(compute_xbars=False)
3671
3672                # we don't technically have to do this at the last
3673                # iteration, but with checkpointing and re-starts,
3674                # you're never sure when you're executing the last
3675                # iteration.
3676                if self._ph_weight_updates_enabled:
3677                    self.update_weights()
3678
3679                # let plugins know if they care.
3680                for plugin in self._ph_plugins:
3681                    plugin.post_iteration_k_solves(self)
3682
3683                if (self._verbose) or (self._report_solutions):
3684                    print("Variable values following scenario solves:")
3685                    self.pprint(False, False, True, False, False,
3686                                output_only_statistics=self._report_only_statistics,
3687                                output_only_nonconverged=self._report_only_nonconverged_variables)
3688
3689                # update the fixed variable statistics.
3690                self._total_fixed_discrete_vars,self._total_fixed_continuous_vars = \
3691                    self.compute_fixed_variable_counts()
3692
3693                print("Number of discrete variables fixed="
3694                      +str(self._total_fixed_discrete_vars)+" "
3695                      "(total="+str(self._total_discrete_vars)+")")
3696                print("Number of continuous variables fixed="
3697                      +str(self._total_fixed_continuous_vars)+" "
3698                      "(total="+str(self._total_continuous_vars)+")")
3699
3700                # update the convergence statistic - prior to the
3701                # plugins callbacks; technically, computing the
3702                # convergence metric is part of the iteration k work
3703                # load.
3704
3705                self._converger.update(self._current_iteration,
3706                                       self,
3707                                       self._scenario_tree,
3708                                       self._instances)
3709                first_stage_min, first_stage_avg, first_stage_max = \
3710                    self._extract_first_stage_cost_statistics()
3711                print("Convergence metric=%12.4f  First stage cost avg=%12.4f  "
3712                      "Max-Min=%8.2f" % (self._converger.lastMetric(),
3713                                         first_stage_avg,
3714                                         first_stage_max-first_stage_min))
3715
3716                expected_cost = self._scenario_tree.findRootNode().computeExpectedNodeCost()
3717                if not _OLD_OUTPUT: print("Expected Cost=%14.4f" % (expected_cost))
3718                self._cost_history[self._current_iteration] = expected_cost
3719                if self._converger.isConverged(self):
3720                    if (len(self._incumbent_cost_history) == 0) or \
3721                       ((self._objective_sense == minimize) and \
3722                        (expected_cost < min(self._incumbent_cost_history.values()))) or \
3723                       ((self._objective_sense == maximize) and \
3724                        (expected_cost > max(self._incumbent_cost_history.values()))):
3725                        if not _OLD_OUTPUT: print("Caching results for new incumbent solution")
3726                        self.cacheSolutions(self._incumbent_cache_id)
3727                        self._best_incumbent_key = self._current_iteration
3728                    self._incumbent_cost_history[self._current_iteration] = expected_cost
3729
3730                # let plugins know if they care.
3731                for plugin in self._ph_plugins:
3732                    plugin.post_iteration_k(self)
3733
3734                # at this point, all the real work of an iteration is
3735                # complete.
3736
3737                # checkpoint if it's time.
3738                if (self._checkpoint_interval > 0) and \
3739                   (i % self._checkpoint_interval is 0):
3740                    self.checkpoint(i)
3741
3742                # everybody wants to know how long they've been waiting...
3743                print("Cumulative run-time=%.2f seconds"
3744                      % (time.time() - self._solve_start_time))
3745
3746                # check for early termination.
3747                if self._dual_mode is False:
3748                    if self._converger.isConverged(self):
3749                        if self._total_discrete_vars == 0:
3750                            print("PH converged - convergence metric is below "
3751                                  "threshold="
3752                                  +str(self._converger._convergence_threshold))
3753                        else:
3754                            print("PH converged - convergence metric is below "
3755                                  "threshold="
3756                                  +str(self._converger._convergence_threshold)+" "
3757                                  "or all discrete variables are fixed")
3758
3759                        plugin_convergence = True
3760                        for plugin in self._ph_plugins:
3761                            if hasattr(plugin,"ph_convergence_check"):
3762                                if not plugin.ph_convergence_check(self):
3763                                    plugin_convergence = False
3764
3765                        if plugin_convergence:
3766                            break
3767                else:
3768                    # This is ugly. We can fix later when we figure
3769                    # out convergence criteria for the dual ph
3770                    # algorithm
3771                    if dual_rc is True:
3772                        break
3773
3774                # if we're terminating due to exceeding the maximum
3775                # iteration count, print a message indicating so -
3776                # otherwise, you get a quiet, information-free output
3777                # trace.
3778                if i == self._max_iterations:
3779                    print("Halting PH - reached maximal iteration count="
3780                          +str(self._max_iterations))
3781
3782                # garbage-collect if it wasn't disabled entirely.
3783                if re_enable_gc:
3784                    if (time.time() - self._time_since_last_garbage_collect) >= \
3785                         self._minimum_garbage_collection_interval:
3786                       gc.collect()
3787                       self._time_since_last_garbage_collect = time.time()
3788
3789                # gather and report memory statistics (for leak
3790                # detection purposes) if specified.
3791                if (guppy_available) and (self._profile_memory >= 1):
3792                    print(hpy().heap())
3793
3794                    #print "New (persistent) objects constructed during PH iteration "+str(self._current_iteration)+":"
3795                    #memory_tracker.print_diff(summary1=summary_last_iteration,
3796                    #                          summary2=summary_this_iteration)
3797
3798                    ## get ready for the next iteration.
3799                    #objects_last_iteration = objects_this_iteration
3800                    #summary_last_iteration = summary_this_iteration
3801
3802                    # XXX begin debugging
3803                    #print "Current type: {0} ({1})".format(type(self), type(self).__name__)
3804                    #print "Recognized objects in muppy:", len(muppy.get_objects())
3805                    #print "Uncollectable objects (cycles?):", gc.garbage
3806
3807                    ##from pympler.muppy import refbrowser
3808                    ##refbrowser.InteractiveBrowser(self).main()
3809
3810                    #print "Referents from PH solver:", gc.get_referents(self)
3811                    #print "Interesting referent keys:", [k for k in gc.get_referents(self)[0].keys() if type(gc.get_referents(self)[0][k]).__name__ not in ['list', 'int', 'str', 'float', 'dict', 'bool']]
3812                    #print "_ph_plugins:", gc.get_referents(self)[0]['_ph_plugins']
3813                    #print "_converger:", gc.get_referents(self)[0]['_converger']
3814                    # XXX end debugging
3815
3816            ####################################################################################################
3817
3818        else:
3819            # TODO: Eliminate this option from the rundph command
3820            if self._dual_mode is True:
3821                raise NotImplementedError("The 'async' option has not been implemented for dual ph.")
3822            ####################################################################################################
3823            self.async_iteration_k_plus_solves()
3824            ####################################################################################################
3825
3826        # re-enable the normal garbage collection mode.
3827        if re_enable_gc:
3828            gc.enable()
3829
3830        print("Number of discrete variables fixed "
3831              "before final plugin calls="
3832              +str(self._total_fixed_discrete_vars)+" "
3833              "(total="+str(self._total_discrete_vars)+")")
3834        print("Number of continuous variables fixed "
3835              "before final plugin calls="
3836              +str(self._total_fixed_continuous_vars)+" "
3837              "(total="+str(self._total_continuous_vars)+")")
3838
3839        if (self._best_incumbent_key is not None) and \
3840           (self._best_incumbent_key != self._current_iteration):
3841            if not _OLD_OUTPUT:
3842                print("")
3843                print("Restoring scenario tree solution "
3844                      "to best incumbent solution with "
3845                      "expected cost=%14.4f"
3846                      % self._incumbent_cost_history[self._best_incumbent_key])
3847            self.restoreCachedSolutions(self._incumbent_cache_id)
3848
3849        if isinstance(self._solver_manager,
3850                      pyomo.solvers.plugins.smanager.phpyro.SolverManager_PHPyro):
3851            collect_full_results(self,
3852                                 TransmitType.all_stages | \
3853                                 TransmitType.blended | \
3854                                 TransmitType.derived | \
3855                                 TransmitType.fixed)
3856
3857        # let plugins know if they care. do this before
3858        # the final solution / statistics output, as the plugins
3859        # might do some final tweaking.
3860        for plugin in self._ph_plugins:
3861            plugin.post_ph_execution(self)
3862
3863        # update the fixed variable statistics - the plugins might have done something.
3864        (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
3865
3866        self._solve_end_time = time.time()
3867
3868        print("PH complete")
3869
3870        if _OLD_OUTPUT:
3871            print("Convergence history:")
3872            self._converger.pprint()
3873        else:
3874            print("")
3875            print("Algorithm History: ")
3876            label_string = "  "
3877            label_string += ("%10s" % "Iteration")
3878            label_string += ("%14s" % "Metric Value")
3879            label_string += ("%17s" % "Expected Cost")
3880            label_string += ("%17s" % "Best Converged")
3881            print(label_string)
3882            best_incumbent_cost = None
3883            for i in xrange(self._current_iteration+1):
3884                row_string = "{0} "
3885                row_string += ("%10d" % i)
3886                metric = self._converger._metric_history[i]
3887                if self._converger._convergence_threshold >= 0.0001:
3888                    row_string += ("%14.4f" % metric)
3889                else:
3890                    row_string += ("%14.3e" % metric)
3891                expected_cost = self._cost_history[i]
3892                row_string += ("%17.4f" % (expected_cost))
3893                if i in self._incumbent_cost_history:
3894                    row_string = row_string.format('*')
3895                    expected_cost = self._incumbent_cost_history[i]
3896                    updated_best = False
3897                    if best_incumbent_cost is None:
3898                        best_incumbent_cost = expected_cost
3899                        updated_best = True
3900                    else:
3901                        if self._objective_sense == minimize:
3902                            if expected_cost < best_incumbent_cost:
3903                                best_incumbent_cost = expected_cost
3904                                updated_best = True
3905                        else:
3906                            if expected_cost > best_incumbent_cost:
3907                                best_incumbent_cost = expected_cost
3908                                updated_best = False
3909                    row_string += ('%17.4f' % (best_incumbent_cost))
3910                    if updated_best:
3911                        row_string += "  (new incumbent)"
3912                else:
3913                    row_string = row_string.format(' ')
3914                    if best_incumbent_cost is not None:
3915                        row_string += ('%17.4f' % (best_incumbent_cost))
3916                    else:
3917                        row_string += ('%17s' % ('-'))
3918                print(row_string)
3919
3920        # print *the* metric of interest.
3921        print("")
3922        print("***********************************************************************************************")
3923        root_node = self._scenario_tree._stages[0]._tree_nodes[0]
3924        print(">>>THE EXPECTED SUM OF THE STAGE COST VARIABLES="+str(root_node.computeExpectedNodeCost())+"<<<")
3925        print(">>>***CAUTION***: Assumes full (or nearly so) convergence of scenario solutions at each node in the scenario tree - computed costs are invalid otherwise<<<")
3926        print("***********************************************************************************************")
3927
3928        print("Final number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")")
3929        print("Final number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")")
3930
3931        # populate the scenario tree solution from the instances - to ensure consistent state
3932        # across the scenario tree instance and the scenario instances.
3933        self._scenario_tree.snapshotSolutionFromScenarios()
3934
3935        print("Final variable values:")
3936        self.pprint(False, False, True, True, False,
3937                    output_only_statistics=self._report_only_statistics,
3938                    output_only_nonconverged=self._report_only_nonconverged_variables)
3939
3940        print("Final costs:")
3941        self._scenario_tree.pprintCosts()
3942
3943        if self._output_scenario_tree_solution:
3944            print("Final solution (scenario tree format):")
3945            self._scenario_tree.pprintSolution()
3946
3947        if (self._verbose) and (self._output_times):
3948            print("Overall run-time=%.2f seconds" % (self._solve_end_time - self._solve_start_time))
3949
3950        for tree_node in self._scenario_tree._tree_nodes:
3951            if tree_node.has_fixed_in_queue() or \
3952               tree_node.has_freed_in_queue():
3953                print("***WARNING***: PH exiting with fixed or freed variables "
3954                      "in the scenario tree queue whose statuses have not been "
3955                      "pushed to the scenario instances. This is because these "
3956                      "requests were placed in the queue after the most recent "
3957                      "solve. If these statuses are pushed to the scenario "
3958                      "instances, a new solution should be obtained in order to "
3959                      "reflect accurate solution data within the scenario tree. "
3960                      "This warning can be safely ignored in most cases.")
3961                break
3962
3963        # Activate the original objective form Don't bother
3964        # transmitting these deactivation signals to the ph solver
3965        # servers as this function is being called at the end of ph
3966        # (for now)
3967        if not isinstance(self._solver_manager,
3968                          pyomo.solvers.plugins.smanager.\
3969                          phpyro.SolverManager_PHPyro):
3970            self.deactivate_ph_objective_weight_terms()
3971            self.deactivate_ph_objective_proximal_terms()
3972
3973    def _clear_bundle_instances(self):
3974
3975        for bundle_name, bundle_instance in iteritems(self._bundle_binding_instance_map):
3976            for scenario_name in self._bundle_scenario_instance_map[bundle_name]:
3977                bundle_instance.del_component(scenario_name)
3978
3979        self._bundle_binding_instance_map = {}
3980        self._bundle_scenario_instance_map = {}
3981    #
3982    # prints a summary of all collected time statistics
3983    #
3984
3985    def print_time_stats(self):
3986
3987        print("PH run-time statistics:")
3988
3989        print("Initialization time=  %.2f seconds" % (self._init_end_time - self._init_start_time))
3990        print("Overall solve time=   %.2f seconds" % (self._solve_end_time - self._solve_start_time))
3991        print("Scenario solve time=  %.2f seconds" % self._cumulative_solve_time)
3992        print("Average update time=  %.2f seconds" % self._cumulative_xbar_time)
3993        print("Weight update time=   %.2f seconds" % self._cumulative_weight_time)
3994
3995    #
3996    # a utility to determine whether to output weight / average / etc. information for
3997    # a variable/node combination. when the printing is moved into a callback/plugin,
3998    # this routine will go there. for now, we don't dive down into the node resolution -
3999    # just the variable/stage.
4000    #
4001
4002    def should_print(self, stage, variable):
4003
4004        if self._output_continuous_variable_stats is False:
4005
4006            variable_type = variable.domain
4007
4008            if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):
4009
4010                return False
4011
4012        return True
4013
4014    #
4015    # pretty-prints the state of the current variable averages, weights, and values.
4016    # inputs are booleans indicating which components should be output.
4017    #
4018
4019    def pprint(self,
4020               output_averages,
4021               output_weights,
4022               output_values,
4023               output_fixed,
4024               output_rhos,
4025               output_only_statistics=False,
4026               output_only_nonconverged=False,
4027               report_stage_costs=True):
4028
4029        if self._initialized is False:
4030            raise RuntimeError("PH is not initialized - cannot invoke "
4031                               "pprint() method")
4032
4033        # print tree nodes and associated variable/xbar/ph information
4034        # in stage-order we don't blend in the last stage, so we don't
4035        # current care about printing the associated information.
4036        for stage in self._scenario_tree._stages[:-1]:
4037
4038            print("   Stage: %s" % (stage._name))
4039
4040            # tracks the number of outputs on a per-index basis.
4041            num_outputs_this_stage = 0
4042
4043            for variable_name, index_template in sorted(iteritems(stage._variables), key=itemgetter(0)):
4044
4045                # track, so we don't output the variable names unless
4046                # there is an entry to report.
4047                num_outputs_this_variable = 0
4048
4049                for tree_node in stage._tree_nodes:
4050
4051                    if not output_only_statistics:
4052                        sys.stdout.write("          (Scenarios: ")
4053                        for scenario in tree_node._scenarios:
4054                            sys.stdout.write(str(scenario._name)+"  ")
4055                            if scenario == tree_node._scenarios[-1]:
4056                                sys.stdout.write(")\n")
4057
4058                    variable_indices = tree_node._variable_indices[variable_name]
4059
4060                    # this is moderately redundant, but shouldn't show
4061                    # up in profiles - printing takes more time than
4062                    # computation. determine the maximimal index
4063                    # string length, so we can output readable column
4064                    # formats.
4065                    max_index_string_length = 0
4066                    for index in variable_indices:
4067                        if index != None:
4068                            this_index_length = len(indexToString(index))
4069                            if this_index_length > max_index_string_length:
4070                                max_index_string_length = this_index_length
4071
4072                    for index in sorted(variable_indices):
4073
4074                        # track, so we don't output the variable index
4075                        # more than once.
4076                        num_outputs_this_index = 0
4077
4078                        # determine if the variable/index pair is used
4079                        # across the set of scenarios (technically, it
4080                        # should be good enough to check one
4081                        # scenario). ditto for "fixed" status. fixed
4082                        # does imply unused (see note below), but we
4083                        # care about the fixed status when outputting
4084                        # final solutions.
4085
4086                        # should be consistent across scenarios, so
4087                        # one "unused" flags as invalid.
4088                        variable_id = \
4089                            tree_node._name_index_to_id[variable_name,index]
4090                        is_fixed = tree_node.is_variable_fixed(variable_id)
4091
4092                        is_not_stale = \
4093                            all((not scenario.is_variable_stale(tree_node,
4094                                                                variable_id)) \
4095                                for scenario in tree_node._scenarios)
4096
4097                        # IMPT: this is far from obvious, but
4098                        #       variables that are fixed will -
4099                        #       because presolve will identify them as
4100                        #       constants and eliminate them from all
4101                        #       expressions - be flagged as "unused"
4102                        #       and therefore not output.
4103
4104                        if ((output_fixed) and (is_fixed)) or \
4105                           ((is_not_stale) and (not is_fixed)):
4106
4107                            minimum_value = tree_node._minimums[variable_id]
4108                            average_value = tree_node._averages[variable_id]
4109                            maximum_value = tree_node._maximums[variable_id]
4110
4111                            # there really isn't a default need to
4112                            # output variables whose values are equal
4113                            # to 0 across-the-board. and there is good
4114                            # reason not to, i.e., the volume of
4115                            # output.
4116                            if ((fabs(minimum_value) > self._integer_tolerance) or \
4117                                (fabs(maximum_value) > self._integer_tolerance)) or\
4118                               (self._report_for_zero_variable_values is True):
4119
4120                                if (fabs(maximum_value - minimum_value) <= \
4121                                    self._integer_tolerance) and \
4122                                   (output_only_nonconverged == True):
4123                                    pass
4124                                else:
4125                                    num_outputs_this_stage = \
4126                                        num_outputs_this_stage + 1
4127                                    num_outputs_this_variable = \
4128                                        num_outputs_this_variable + 1
4129                                    num_outputs_this_index = \
4130                                        num_outputs_this_index + 1
4131
4132                                    if num_outputs_this_variable == 1:
4133                                        sys.stdout.write("      Variable: "
4134                                                         + variable_name+'\n')
4135
4136                                    if num_outputs_this_index == 1:
4137                                        if index is not None:
4138                                            format_string = \
4139                                                ("         Index: %"
4140                                                 +str(max_index_string_length)+"s")
4141                                            sys.stdout.write(format_string
4142                                                             % indexToString(index))
4143
4144                                    if len(stage._tree_nodes) > 1:
4145                                        sys.stdout.write("\n")
4146                                        sys.stdout.write("         Tree Node: %s"
4147                                                         % (tree_node._name))
4148
4149                                    if output_values:
4150                                        if output_only_statistics is False:
4151                                            sys.stdout.write("\tValues:  ")
4152                                        last_scenario = tree_node._scenarios[-1]
4153                                        for scenario in tree_node._scenarios:
4154                                            scenario_probability = \
4155                                                scenario._probability
4156                                            this_value = \
4157                                                scenario._x[tree_node._name]\
4158                                                           [variable_id]
4159                                            if not output_only_statistics:
4160                                                sys.stdout.write("%12.4f"
4161                                                                 % this_value)
4162                                            if scenario is last_scenario:
4163                                                if output_only_statistics:
4164                                                    # there
4165                                                    # technically
4166                                                    # isn't any good
4167                                                    # reason not to
4168                                                    # always report
4169                                                    # the min and max;
4170                                                    # the only reason
4171                                                    # we're not doing
4172                                                    # this currently
4173                                                    # is to avoid
4174                                                    # updating our
4175                                                    # regression test
4176                                                    # baseline output.
4177                                                    sys.stdout.write(
4178                                                        "    Min:  %12.4f"
4179                                                        % (minimum_value))
4180                                                    sys.stdout.write(
4181                                                        "    Avg:  %12.4f"
4182                                                        % (average_value))
4183                                                    sys.stdout.write(
4184                                                        "    Max:  %12.4f"
4185                                                        % (maximum_value))
4186                                                else:
4187                                                    sys.stdout.write(
4188                                                        "    Max-Min:  %12.4f"
4189                                                        % (maximum_value - \
4190                                                           minimum_value))
4191                                                    sys.stdout.write(
4192                                                        "    Avg:  %12.4f"
4193                                                        % (average_value))
4194                                                sys.stdout.write("\n")
4195                                    if output_weights:
4196                                        sys.stdout.write("         Weights:  ")
4197                                        for scenario in tree_node._scenarios:
4198                                            sys.stdout.write(
4199                                                "%12.4f"
4200                                                % scenario._w[tree_node._name]\
4201                                                             [variable_id])
4202                                    if output_rhos:
4203                                        sys.stdout.write("         Rhos:  ")
4204                                        for scenario in tree_node._scenarios:
4205                                            sys.stdout.write(
4206                                                "%12.4f"
4207                                                % scenario._rho[tree_node._name]\
4208                                                               [variable_id])
4209
4210                                    if output_averages:
4211                                        sys.stdout.write("   Average:  %12.4f"
4212                                                         % (average_value))
4213                                    if output_weights or output_rhos or output_values:
4214                                        sys.stdout.write("\n")
4215            if num_outputs_this_stage == 0:
4216                print("\t\tNo non-converged variables in stage")
4217
4218            if not report_stage_costs:
4219                continue
4220
4221            # cost variables aren't blended, so go through the gory
4222            # computation of min/max/avg.  we currently always print
4223            # these.
4224            cost_variable_name = stage._cost_variable[0]
4225            cost_variable_index = stage._cost_variable[1]
4226            print("      Cost Variable: "
4227                  +cost_variable_name+indexToString(cost_variable_index))
4228            for tree_node in stage._tree_nodes:
4229                sys.stdout.write("         Tree Node: %s"
4230                                 % (tree_node._name))
4231                if output_only_statistics is False:
4232                    sys.stdout.write("      (Scenarios:  ")
4233                    for scenario in tree_node._scenarios:
4234                        sys.stdout.write(str(scenario._name)+" ")
4235                        if scenario == tree_node._scenarios[-1]:
4236                            sys.stdout.write(")\n")
4237                maximum_value = 0.0
4238                minimum_value = 0.0
4239                prob_sum_values = 0.0
4240                sum_prob = 0.0
4241                num_values = 0
4242                first_time = True
4243                if output_only_statistics is False:
4244                    sys.stdout.write("         Values:  ")
4245                else:
4246                    sys.stdout.write("         ")
4247                for scenario in tree_node._scenarios:
4248                    this_value = scenario._stage_costs[stage._name]
4249                    if output_only_statistics is False:
4250                        if this_value is not None:
4251                            sys.stdout.write("%12.4f" % this_value)
4252                        else:
4253                            # this is a hack, in case the stage cost
4254                            # variables are not returned. ipopt does
4255                            # this occasionally, for example, if stage
4256                            # cost variables are constrained to a
4257                            # constant value (and consequently
4258                            # preprocessed out).
4259                            sys.stdout.write("%12s" % "Not Reported")
4260                    if this_value is not None:
4261                        num_values += 1
4262                        prob_sum_values += scenario._probability * this_value
4263                        sum_prob += scenario._probability
4264                        if first_time:
4265                            first_time = False
4266                            maximum_value = this_value
4267                            minimum_value = this_value
4268                        else:
4269                            if this_value > maximum_value:
4270                                maximum_value = this_value
4271                            if this_value < minimum_value:
4272                                minimum_value = this_value
4273                    if scenario == tree_node._scenarios[-1]:
4274                        if num_values > 0:
4275                            if output_only_statistics:
4276                                sys.stdout.write("    Min:  %12.4f"
4277                                                 % (minimum_value))
4278                                if (sum_prob > 0.0):
4279                                    sys.stdout.write("    Avg:  %12.4f"
4280                                                     % (prob_sum_values/sum_prob))
4281                                sys.stdout.write("    Max:  %12.4f"
4282                                                 % (maximum_value))
4283                            else:
4284                                sys.stdout.write("    Max-Min:  %12.4f"
4285                                                 % (maximum_value-minimum_value))
4286                                if (sum_prob > 0.0):
4287                                    sys.stdout.write("   Avg:  %12.4f"
4288                                                     % (prob_sum_values/sum_prob))
4289                        sys.stdout.write("\n")
Note: See TracBrowser for help on using the repository browser.