Introduce VectorArray views

docs/source/technical_overview.rst

@@ -16,30 +16,30 @@ operating on objects of the following types:
     |copied| to a new array, |appended| to an existing array or |removed| from the
     array. Basic linear algebra operations can be performed on the vectors of the
     array: vectors can be |scaled| in-place, the BLAS |axpy| operation is
-    supported and |scalar products| between vectors can be formed. Linear
+    supported and |inner products| between vectors can be formed. Linear
     combinations of vectors can be formed using the |lincomb| method. Moreover,
     various norms can be computed and selected |components| of the vectors can
     be extracted for :mod:`empirical interpolation <pymor.algorithms.ei>`.
-
-    Each of these methods takes optional `ind` parameters to specify the subset
-    of vectors on which to operate. If the parameter is not specified, the whole
-    array is selected for the operation. 
+    To act on subsets of vectors of an array, arrays can be |indexed|, returning
+    a new |VectorArray| which acts as a view onto the respective vectors in the
+    original array.
 
     New vector arrays can be created using the |empty| and |zeros| method. As a
     convenience, many of Python's math special methods are implemented in terms
     of the interface methods.
 
     Note that there is not the notion of a single vector in pyMOR. The main
     reason for this design choice is to take advantage of vectorized
-    implementations like |NumpyVectorArray| which internally store the
-    vectors as two-dimensional |NumPy| arrays. As an example, the application of
-    a linear matrix based operator to an array via the |apply| method boils down
-    to a call to |NumPy|'s optimized :meth:`~numpy.ndarray.dot` method. If there
-    were only lists of vectors in pyMOR, the above matrix-matrix multiplication
-    would have to be expressed by a loop of matrix-vector multiplications. However,
-    when working with external solvers, vector arrays will often be just lists
-    of vectors. For this use-case we provide |ListVectorArray|, a vector array
-    based on a Python list of vectors.
+    implementations like |NumpyVectorArray| which internally store the vectors
+    as two-dimensional |NumPy| arrays. As an example, the application of a
+    linear matrix based operator to an array via the |apply| method boils down
+    to a call to |NumPy|'s optimized :meth:`~numpy.ndarray.dot` method. If
+    there were only lists of vectors in pyMOR, the above matrix-matrix
+    multiplication would have to be expressed by a loop of matrix-vector
+    multiplications.  However, when working with external solvers, vector
+    arrays will often be given as lists of indiviual vector objects. For this
+    use-case we provide |ListVectorArray|, a |VectorArray| based on a Python
+    list of vectors.
 
     Associated to each vector array is a |VectorSpace|. A Vector space in pyMOR
     is simply the combination of a |VectorArray| subclass and an appropriate
@@ -62,10 +62,11 @@ operating on objects of the following types:
     .. |copied|           replace:: :meth:`copied <pymor.vectorarrays.interfaces.VectorArrayInterface.copy>`
     .. |dimension|        replace:: :attr:`dimension <pymor.vectorarrays.interfaces.VectorArrayInterface.dim>`
     .. |empty|            replace:: :meth:`~pymor.vectorarrays.interfaces.VectorArrayInterface.empty`
+    .. |indexed|          replace:: :meth:`indexed <pymor.vectorarrays.interfaces.VectorArrayInterface.__getitem__>`
+    .. |inner products|   replace:: :meth:`inner products <pymor.vectorarrays.interfaces.VectorArrayInterface.dot>`
     .. |lincomb|          replace:: :meth:`~pymor.vectorarrays.interfaces.VectorArrayInterface.lincomb`
     .. |make_array|       replace:: :meth:`~pymor.vectorarrays.interfaces.VectorArrayInterface.make_array`
-    .. |removed|          replace:: :meth:`deleted <pymor.vectorarrays.interfaces.VectorArrayInterface.remove>`
-    .. |scalar products|  replace:: :meth:`scalar products <pymor.vectorarrays.interfaces.VectorArrayInterface.dot>`
+    .. |removed|          replace:: :meth:`deleted <pymor.vectorarrays.interfaces.VectorArrayInterface.__delitem__>`
     .. |scaled|           replace:: :meth:`scaled <pymor.vectorarrays.interfaces.VectorArrayInterface.scal>`
     .. |subtype|          replace:: :attr:`~pymor.vectorarrays.interfaces.VectorSpace.subtype`
     .. |zeros|            replace:: :meth:`~pymor.vectorarrays.interfaces.VectorArrayInterface.zeros`

src/pymor/algorithms/basic.py

@@ -12,7 +12,7 @@
 
 
 @defaults('rtol', 'atol')
-def almost_equal(U, V, U_ind=None, V_ind=None, product=None, norm=None, rtol=1e-14, atol=1e-14):
+def almost_equal(U, V, product=None, norm=None, rtol=1e-14, atol=1e-14):
     """Compare U and V for almost equality.
 
     The vectors of `U` and `V` are compared in pairs for almost equality.
@@ -23,16 +23,14 @@ def almost_equal(U, V, U_ind=None, V_ind=None, product=None, norm=None, rtol=1e-
     The norm to be used can be specified via the `norm` or `product`
     parameter.
 
-    If the length of `U` (`U_ind`) resp. `V` (`V_ind`) is 1, the single
-    specified vector is compared to all vectors of the other array.
+    If the length of `U`  resp. `V`  is 1, the single specified
+    vector is compared to all vectors of the other array.
     Otherwise, the lengths of both indexed arrays have to agree.
 
     Parameters
     ----------
     U, V
         |VectorArrays| to be compared.
-    U_ind, V_ind
-        Indices of the vectors that are to be compared (see |VectorArray|).
     product
         If specified, use this inner product |Operator| to compute the norm.
         `product` and `norm` are mutually exclusive.
@@ -57,16 +55,15 @@ def almost_equal(U, V, U_ind=None, V_ind=None, product=None, norm=None, rtol=1e-
         norm_str = norm
         norm = lambda U: getattr(U, norm_str + '_norm')()
 
-    X = V.copy(V_ind)
+    X = V.copy()
     V_norm = norm(X)
 
     # broadcast if necessary
     if len(X) == 1:
-        len_U = U.len_ind(U_ind)
-        if len_U > 1:
-            X.append(X, o_ind=np.zeros(len_U - 1, dtype=np.int))
+        if len(U) > 1:
+            X.append(X[np.zeros(len(U) - 1, dtype=np.int)])
 
-    X.axpy(-1, U, x_ind=U_ind)
+    X -= U
     ERR_norm = norm(X)
 
     return ERR_norm <= atol + V_norm * rtol

src/pymor/algorithms/basisextension.py

@@ -65,11 +65,11 @@ def trivial_basis_extension(basis, U, copy_basis=True, copy_U=True):
     old_basis_length = len(basis)
     remove = set()
     for i in range(len(U)):
-        if np.any(almost_equal(U, basis, U_ind=i)):
+        if np.any(almost_equal(U[i], basis)):
             remove.add(i)
 
     new_basis = basis.copy() if copy_basis else basis
-    new_basis.append(U, o_ind=[i for i in range(len(U)) if i not in remove],
+    new_basis.append(U[[i for i in range(len(U)) if i not in remove]],
                      remove_from_other=(not copy_U))
 
     if len(new_basis) == old_basis_length:

src/pymor/algorithms/ei.py

@@ -139,7 +139,7 @@ def ei_greedy(U, error_norm=None, atol=None, rtol=None, max_interpolation_dofs=N
             break
 
         # compute new interpolation dof and collateral basis vector
-        new_vec = U.copy(ind=max_err_ind)
+        new_vec = U[max_err_ind].copy()
         new_dof = new_vec.amax()[0][0]
         if new_dof in interpolation_dofs:
             logger.info('DOF {} selected twice for interplation! Stopping extension loop.'.format(new_dof))
@@ -160,8 +160,8 @@ def ei_greedy(U, error_norm=None, atol=None, rtol=None, max_interpolation_dofs=N
         if projection == 'orthogonal':
             onb_collateral_basis.append(new_vec)
             gram_schmidt(onb_collateral_basis, offset=len(onb_collateral_basis) - 1, copy=False)
-            coeffs = ERR.dot(onb_collateral_basis, o_ind=len(onb_collateral_basis) - 1)
-            ERR.axpy(-coeffs[:, 0], onb_collateral_basis, x_ind=len(onb_collateral_basis) - 1)
+            coeffs = ERR.dot(onb_collateral_basis[len(onb_collateral_basis) - 1])
+            ERR.axpy(-coeffs[:, 0], onb_collateral_basis[len(onb_collateral_basis) - 1])
         errs = ERR.l2_norm() if error_norm is None else error_norm(ERR)
         max_err_ind = np.argmax(errs)
         max_err = errs[max_err_ind]
@@ -229,12 +229,12 @@ def deim(U, modes=None, error_norm=None, product=None):
 
         if len(interpolation_dofs) > 0:
             coefficients = np.linalg.solve(interpolation_matrix,
-                                           collateral_basis.components(interpolation_dofs, ind=i).T).T
-            U_interpolated = collateral_basis.lincomb(coefficients, ind=list(range(len(interpolation_dofs))))
-            ERR = collateral_basis.copy(ind=i)
+                                           collateral_basis[i].components(interpolation_dofs).T).T
+            U_interpolated = collateral_basis[:len(interpolation_dofs)].lincomb(coefficients)
+            ERR = collateral_basis[i].copy()
             ERR -= U_interpolated
         else:
-            ERR = collateral_basis.copy(ind=i)
+            ERR = collateral_basis[i].copy()
 
         err = np.max(ERR.l2_norm() if error_norm is None else error_norm(ERR))
 
@@ -248,13 +248,13 @@ def deim(U, modes=None, error_norm=None, product=None):
             break
 
         interpolation_dofs = np.hstack((interpolation_dofs, new_dof))
-        interpolation_matrix = collateral_basis.components(interpolation_dofs, ind=list(range(len(interpolation_dofs)))).T
+        interpolation_matrix = collateral_basis[:len(interpolation_dofs)].components(interpolation_dofs).T
         errs.append(err)
 
         logger.info('')
 
     if len(interpolation_dofs) < len(collateral_basis):
-        collateral_basis.remove(ind=list(range(len(interpolation_dofs), len(collateral_basis))))
+        del collateral_basis[len(interpolation_dofs):len(collateral_basis)]
 
     logger.info('Finished.')
 
@@ -444,7 +444,7 @@ def _parallel_ei_greedy_initialize(U=None, error_norm=None, copy=None, data=None
 
 
 def _parallel_ei_greedy_get_vector(data=None):
-    return data['U'].copy(ind=data['max_err_ind'])
+    return data['U'][data['max_err_ind']].copy()
 
 
 def _parallel_ei_greedy_update(new_vec=None, new_dof=None, data=None):

src/pymor/algorithms/genericsolvers.py

@@ -170,7 +170,7 @@ def apply_inverse(op, rhs, options=None):
 
     if options['type'] == 'generic_lgmres':
         for i in range(len(rhs)):
-            r, info = lgmres(op, rhs.copy(i),
+            r, info = lgmres(op, rhs[i],
                              tol=options['tol'],
                              maxiter=options['maxiter'],
                              inner_m=options['inner_m'],
@@ -181,7 +181,7 @@ def apply_inverse(op, rhs, options=None):
             R.append(r)
     elif options['type'] == 'least_squares_generic_lsmr':
         for i in range(len(rhs)):
-            r, info, itn, _, _, _, _, _ = lsmr(op, rhs.copy(i),
+            r, info, itn, _, _, _, _, _ = lsmr(op, rhs[i],
                                                damp=options['damp'],
                                                atol=options['atol'],
                                                btol=options['btol'],
@@ -195,7 +195,7 @@ def apply_inverse(op, rhs, options=None):
             R.append(r)
     elif options['type'] == 'least_squares_generic_lsqr':
         for i in range(len(rhs)):
-            r, info, itn, _, _, _, _, _, _ = lsqr(op, rhs.copy(i),
+            r, info, itn, _, _, _, _, _, _ = lsqr(op, rhs[i],
                                                   damp=options['damp'],
                                                   atol=options['atol'],
                                                   btol=options['btol'],

src/pymor/algorithms/gram_schmidt.py

@@ -59,17 +59,17 @@ def gram_schmidt(A, product=None, atol=1e-13, rtol=1e-13, offset=0, find_duplica
     for i in range(offset, len(A)):
         # first calculate norm
         if product is None:
-            initial_norm = A.l2_norm(ind=i)[0]
+            initial_norm = A[i].l2_norm()[0]
         else:
-            initial_norm = np.sqrt(product.pairwise_apply2(A, A, V_ind=i, U_ind=i))[0]
+            initial_norm = np.sqrt(product.pairwise_apply2(A[i], A[i]))[0]
 
         if initial_norm < atol:
             logger.info("Removing vector {} of norm {}".format(i, initial_norm))
             remove.append(i)
             continue
 
         if i == 0:
-            A.scal(1/initial_norm, ind=0)
+            A[0].scal(1/initial_norm)
 
         else:
             first_iteration = True
@@ -88,16 +88,16 @@ def gram_schmidt(A, product=None, atol=1e-13, rtol=1e-13, offset=0, find_duplica
                     if j in remove:
                         continue
                     if product is None:
-                        p = A.pairwise_dot(A, ind=i, o_ind=j)[0]
+                        p = A[i].pairwise_dot(A[j])[0]
                     else:
-                        p = product.pairwise_apply2(A, A, V_ind=i, U_ind=j)[0]
-                    A.axpy(-p, A, ind=i, x_ind=j)
+                        p = product.pairwise_apply2(A[i], A[j])[0]
+                    A[i].axpy(-p, A[j])
 
                 # calculate new norm
                 if product is None:
-                    old_norm, norm = norm, A.l2_norm(ind=i)[0]
+                    old_norm, norm = norm, A[i].l2_norm()[0]
                 else:
-                    old_norm, norm = norm, np.sqrt(product.pairwise_apply2(A, A, V_ind=i, U_ind=i))[0]
+                    old_norm, norm = norm, np.sqrt(product.pairwise_apply2(A[i], A[i])[0])
 
                 # remove vector if it got too small:
                 if norm / initial_norm < rtol:
@@ -106,16 +106,16 @@ def gram_schmidt(A, product=None, atol=1e-13, rtol=1e-13, offset=0, find_duplica
                     break
 
             if norm > 0:
-                A.scal(1 / norm, ind=i)
+                A[i].scal(1 / norm)
 
     if remove:
-        A.remove(remove)
+        del A[remove]
 
     if check:
         if product:
-            error_matrix = product.apply2(A, A, V_ind=list(range(offset, len(A))))
+            error_matrix = product.apply2(A[offset:len(A)], A)
         else:
-            error_matrix = A.dot(A, ind=list(range(offset, len(A))))
+            error_matrix = A[offset:len(A)].dot(A)
         error_matrix[:len(A) - offset, offset:len(A)] -= np.eye(len(A) - offset)
         if error_matrix.size > 0:
             err = np.max(np.abs(error_matrix))

src/pymor/algorithms/image.py

@@ -208,11 +208,11 @@ def estimate_image_hierarchical(operators=(), vectors=(), domain=None, extends=N
     if extends:
         image = extends[0]
         image_dims = extends[1]
-        ind_range = list(range(len(image_dims) - 1, len(domain))) if operators else list(range(len(image_dims) - 1, 0))
+        ind_range = range(len(image_dims) - 1, len(domain)) if operators else range(len(image_dims) - 1, 0)
     else:
         image = image_space.empty()
         image_dims = []
-        ind_range = list(range(-1, len(domain))) if operators else [-1]
+        ind_range = range(-1, len(domain)) if operators else [-1]
 
     for i in ind_range:
         logger.info('Estimating image for basis vector {} ...'.format(i))
@@ -221,7 +221,7 @@ def estimate_image_hierarchical(operators=(), vectors=(), domain=None, extends=N
                                        orthonormalize=False, product=product,
                                        riesz_representatives=riesz_representatives)
         else:
-            new_image = estimate_image(operators, [], domain.copy(i), extends=True,
+            new_image = estimate_image(operators, [], domain[i], extends=True,
                                        orthonormalize=False, product=product,
                                        riesz_representatives=riesz_representatives)