Merge branch 'use_old_rd'

src/pymor/algorithms/basisextension.py

@@ -51,7 +51,7 @@ def trivial_basis_extension(basis, U, U_ind=None, copy_basis=True, copy_U=True):
     new_basis = basis.copy() if copy_basis else basis
     new_basis.append(U, o_ind=U_ind, remove_from_other=(not copy_U))
 
-    return new_basis
+    return new_basis, {'hierarchic': True}
 
 
 def numpy_trivial_basis_extension(basis, U):
@@ -78,7 +78,7 @@ def numpy_trivial_basis_extension(basis, U):
     '''
     assert isinstance(U, NumpyVectorArray)
     if basis is None:
-        return U
+        return U, {'hierarchic': True}
     assert isinstance(basis, NumpyVectorArray)
     basis = basis.data
     U = U.data
@@ -90,7 +90,7 @@ def numpy_trivial_basis_extension(basis, U):
     new_basis[:-1, :] = basis
     new_basis[-1, :] = U
 
-    return NumpyVectorArray(new_basis)
+    return NumpyVectorArray(new_basis), {'hierarchic': True}
 
 
 def gram_schmidt_basis_extension(basis, U, U_ind=None, product=None, copy_basis=True, copy_U=True):
@@ -135,7 +135,7 @@ def gram_schmidt_basis_extension(basis, U, U_ind=None, product=None, copy_basis=
     if len(new_basis) <= basis_length:
         raise ExtensionError
 
-    return new_basis
+    return new_basis, {'hierarchic': True}
 
 
 def pod_basis_extension(basis, U, count=1, copy_basis=True, product=None):
@@ -167,7 +167,7 @@ def pod_basis_extension(basis, U, count=1, copy_basis=True, product=None):
         due to rounding errors ...
     '''
     if basis is None:
-        return pod(U, modes=count, product=product)
+        return pod(U, modes=count, product=product), {'hierarchic': True}
 
     basis_length = len(basis)
 
@@ -183,4 +183,4 @@ def pod_basis_extension(basis, U, count=1, copy_basis=True, product=None):
     if len(new_basis) <= basis_length:
         raise ExtensionError
 
-    return new_basis
+    return new_basis, {'hierarchic': True}

src/pymor/algorithms/greedy.py

@@ -77,10 +77,12 @@ def greedy(discretization, reductor, samples, initial_data=None, use_estimator=T
     extensions = 0
     max_errs = []
     max_err_mus = []
+    hierarchic = False
 
     while True:
         logger.info('Reducing ...')
-        rd, rc = reductor(discretization, data)
+        rd, rc, reduction_data = reductor(discretization, data) if not hierarchic \
+            else reductor(discretization, data, extends=(rd, rc, reduction_data))
 
         logger.info('Estimating errors ...')
         if use_estimator:
@@ -102,22 +104,28 @@ def greedy(discretization, reductor, samples, initial_data=None, use_estimator=T
         logger.info('Extending with snapshot for mu = {}'.format(max_err_mu))
         U = discretization.solve(max_err_mu)
         try:
-            data = extension_algorithm(data, U)
+            data, extension_data = extension_algorithm(data, U)
         except ExtensionError:
             logger.info('Extension failed. Stopping now.')
             break
         extensions += 1
+        if not 'hierarchic' in extension_data:
+            logger.warn('Extension algorithm does not report if extension was hierarchic. Assuming it was\'nt ..')
+            hierarchic = False
+        else:
+            hierarchic = extension_data['hierarchic']
 
         logger.info('')
 
         if max_extensions is not None and extensions >= max_extensions:
             logger.info('Maximal number of {} extensions reached.'.format(max_extensions))
             logger.info('Reducing once more ...')
-            rd, rc = reductor(discretization, data)
+            rd, rc, reduction_data = reductor(discretization, data) if not hierarchic \
+                else reductor(discretization, data, extends=(rd, rc, reduction_data))
             break
 
     tictoc = time.time() - tic
     logger.info('Greedy search took {} seconds'.format(tictoc))
     return {'data': data, 'reduced_discretization': rd, 'reconstructor': rc, 'max_err': max_err,
             'max_err_mu': max_err_mu, 'max_errs': max_errs, 'max_err_mus': max_err_mus, 'extensions': extensions,
-            'time': tictoc}
+            'time': tictoc, 'reduction_data': reduction_data}

src/pymor/demos/burgers_ei.py

@@ -167,8 +167,8 @@ def burgers_demo(args):
 
     print('RB generation ...')
 
-    def reductor(discretization, rb):
-        return reduce_generic_rb(ei_discretization, rb)
+    def reductor(discretization, rb, extends=None):
+        return reduce_generic_rb(ei_discretization, rb, extends=extends)
 
     extension_algorithm = partial(pod_basis_extension)
 
@@ -187,7 +187,7 @@ def reductor(discretization, rb):
 
     def error_analysis(N, M):
         print('N = {}, M = {}: '.format(N, M), end='')
-        rd, rc = reduce_to_subbasis(rb_discretization, N, reconstructor)
+        rd, rc, _ = reduce_to_subbasis(rb_discretization, N, reconstructor)
         rd = rd.with_(operator=rd.operator.projected_to_subbasis(dim_collateral=M))
         l2_err_max = -1
         mumax = None

src/pymor/demos/thermalblock.py

@@ -153,7 +153,7 @@ def error_analysis(d, rd, rc, mus):
         Ns = np.linspace(1, real_rb_size, N_count).astype(np.int)
     else:
         Ns = np.array([real_rb_size])
-    rd_rcs = [reduce_to_subbasis(rb_discretization, N, reconstructor) for N in Ns]
+    rd_rcs = [reduce_to_subbasis(rb_discretization, N, reconstructor)[:2] for N in Ns]
     mus = list(discretization.parameter_space.sample_randomly(args['--test']))
 
     errs, err_mus, ests, est_mus, conds, cond_mus = zip(*(error_analysis(discretization, rd, rc, mus) for rd, rc in rd_rcs))

src/pymor/demos/thermalblock_pod.py

@@ -101,7 +101,7 @@ def thermalblock_demo(args):
 
     print('Reducing ...')
     reductor = reduce_generic_rb
-    rb_discretization, reconstructor = reductor(discretization, rb)
+    rb_discretization, reconstructor, _ = reductor(discretization, rb)
 
     toc = time.time()
     t_offline = toc - tic

src/pymor/reductors/basic.py

@@ -25,7 +25,8 @@ def restricted_to_subbasis(self, dim):
         return GenericRBReconstructor(self.RB.copy(ind=range(dim)))
 
 
-def reduce_generic_rb(discretization, RB, product=None, disable_caching=True):
+def reduce_generic_rb(discretization, RB, product=None, disable_caching=True,
+                      extends=None):
     '''Generic reduced basis reductor.
 
     Reduces a discretization by applying `operators.rb_project_operator` to
@@ -52,6 +53,7 @@ def reduce_generic_rb(discretization, RB, product=None, disable_caching=True):
         The reconstructor providing a `reconstruct(U)` method which reconstructs
         high-dimensional solutions from solutions U of the reduced discretization.
     '''
+    assert extends is None or len(extends) == 3
 
     if RB is None:
         RB = discretization.type_solution.empty(discretization.dim_solution)
@@ -72,7 +74,7 @@ def reduce_generic_rb(discretization, RB, product=None, disable_caching=True):
     rd.disable_logging()
     rc = GenericRBReconstructor(RB)
 
-    return rd, rc
+    return rd, rc, {}
 
 
 class SubbasisReconstructor(core.BasicInterface):
@@ -130,4 +132,4 @@ def estimate(self, U, mu=None, discretization=None):
         rc = SubbasisReconstructor(next(discretization.operators.itervalues()).dim_source, dim,
                                    old_recontructor=reconstructor)
 
-    return rd, rc
+    return rd, rc, {}

src/pymor/reductors/linear.py

@@ -5,6 +5,7 @@
 from __future__ import absolute_import, division, print_function
 
 import types
+from itertools import izip
 
 import numpy as np
 
@@ -15,7 +16,8 @@
 from pymor.reductors.basic import reduce_generic_rb
 
 
-def reduce_stationary_affine_linear(discretization, RB, error_product=None, disable_caching=True):
+def reduce_stationary_affine_linear(discretization, RB, error_product=None, disable_caching=True,
+                                    extends=None):
     '''Reductor for stationary linear problems whose `operator` and `rhs` are affinely decomposed.
 
     We simply use reduce_generic_rb for the actual RB-projection. The only addition
@@ -54,9 +56,16 @@ def reduce_stationary_affine_linear(discretization, RB, error_product=None, disa
     if discretization.rhs.parametric:
         assert isinstance(discretization.rhs, LincombOperatorInterface)
         assert all(not op.parametric for op in discretization.rhs.operators)
+    assert extends is None or len(extends) == 3
 
     d = discretization
-    rd, rc = reduce_generic_rb(d, RB, product=None, disable_caching=disable_caching)
+    rd, rc, data = reduce_generic_rb(d, RB, product=None, disable_caching=disable_caching,
+                                     extends=extends)
+    if extends:
+        old_data = extends[2]
+        old_RB_size = len(extends[1].RB)
+    else:
+        old_RB_size = 0
 
     # compute data for estimator
     space_dim = d.operator.dim_source
@@ -81,7 +90,9 @@ def append_vector(U, R, RR):
     if RB is None:
         RB = discretization.type_solution.empty(discretization.dim_solution)
 
-    if not d.rhs.parametric:
+    if extends:
+        R_R, RR_R = old_data['R_R'], old_data['RR_R']
+    elif not d.rhs.parametric:
         R_R = space_type.empty(space_dim, reserve=1)
         RR_R = space_type.empty(space_dim, reserve=1)
         append_vector(d.rhs.as_vector(), R_R, RR_R)
@@ -92,24 +103,29 @@ def append_vector(U, R, RR):
             append_vector(op.as_vector(), R_R, RR_R)
 
     if len(RB) == 0:
-        R_O = space_type.empty(space_dim)
-        RR_O = space_type.empty(space_dim)
+        R_Os = [space_type.empty(space_dim)]
+        RR_Os = [space_type.empty(space_dim)]
     elif not d.operator.parametric:
-        R_O = space_type.empty(space_dim, reserve=len(RB))
-        RR_O = space_type.empty(space_dim, reserve=len(RB))
+        R_Os = [space_type.empty(space_dim, reserve=len(RB))]
+        RR_Os = [space_type.empty(space_dim, reserve=len(RB))]
         for i in xrange(len(RB)):
-            append_vector(-d.operator.apply(RB, ind=i), R_O, RR_O)
+            append_vector(-d.operator.apply(RB, ind=i), R_Os[0], RR_Os[0])
     else:
-        R_O = space_type.empty(space_dim, reserve=len(d.operator.operators) * len(RB))
-        RR_O = space_type.empty(space_dim, reserve=len(d.operator.operators) * len(RB))
-        for op in d.operator.operators:
-            for i in xrange(len(RB)):
+        R_Os = [space_type.empty(space_dim, reserve=len(RB)) for _ in xrange(len(d.operator.operators))]
+        RR_Os = [space_type.empty(space_dim, reserve=len(RB)) for _ in xrange(len(d.operator.operators))]
+        if old_RB_size > 0:
+            for op, R_O, RR_O, old_R_O, old_RR_O in izip(d.operator.operators, R_Os, RR_Os,
+                                                         old_data['R_Os'], old_data['RR_Os']):
+                R_O.append(old_R_O)
+                RR_O.append(old_RR_O)
+        for op, R_O, RR_O in izip(d.operator.operators, R_Os, RR_Os):
+            for i in xrange(old_RB_size, len(RB)):
                 append_vector(-op.apply(RB, [i]), R_O, RR_O)
 
     # compute Gram matrix of the residuals
     R_RR = RR_R.dot(R_R, pairwise=False)
-    R_RO = RR_R.dot(R_O, pairwise=False)
-    R_OO = RR_O.dot(R_O, pairwise=False)
+    R_RO = np.hstack([RR_R.dot(R_O, pairwise=False) for R_O in R_Os])
+    R_OO = np.vstack([np.hstack([RR_O.dot(R_O, pairwise=False) for R_O in R_Os]) for RR_O in RR_Os])
 
     estimator_matrix = np.empty((len(R_RR) + len(R_OO),) * 2)
     estimator_matrix[:len(R_RR), :len(R_RR)] = R_RR
@@ -121,8 +137,9 @@ def append_vector(U, R, RR):
 
     estimator = StationaryAffineLinearReducedEstimator(estimator_matrix)
     rd = rd.with_(estimator=estimator)
+    data.update(R_R=R_R, RR_R=RR_R, R_Os=R_Os, RR_Os=RR_Os)
 
-    return rd, rc
+    return rd, rc, data
 
 
 class StationaryAffineLinearReducedEstimator(ImmutableInterface):