# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
from numbers import Number
import numpy as np
from scipy.sparse import issparse
from pymor.core.base import classinstancemethod
from pymor.tools.random import get_random_state
from pymor.vectorarrays.interface import VectorArray, VectorSpace, _create_random_values
[docs]class NumpyVectorArray(VectorArray):
"""|VectorArray| implementation via |NumPy arrays|.
This is the default |VectorArray| type used by all |Operators|
in pyMOR's discretization toolkit. Moreover, all reduced |Operators|
are based on |NumpyVectorArray|.
This class is just a thin wrapper around the underlying
|NumPy array|. Thus, while operations like
:meth:`~pymor.vectorarrays.interface.VectorArray.axpy` or
:meth:`~pymor.vectorarrays.interface.VectorArray.dot`
will be quite efficient, removing or appending vectors will
be costly.
.. warning::
This class is not intended to be instantiated directly. Use
the associated :class:`VectorSpace <NumpyVectorSpace>` instead.
"""
def __init__(self, array, space):
self._array = array
self.space = space
self._refcount = [1]
self._len = len(array)
[docs] def to_numpy(self, ensure_copy=False):
if ensure_copy:
return self._array[:self._len].copy()
if self._refcount[0] > 1:
# While changing the returned NumPy array may change the data of `self`,
# it should not change the data of other arrays.
self._deep_copy()
return self._array[:self._len]
@property
def _data(self):
"""Return NumPy Array view on data for hacking / interactive use."""
return self.to_numpy()
@property
def real(self):
return NumpyVectorArray(self.to_numpy().real.copy(), self.space)
@property
def imag(self):
return NumpyVectorArray(self.to_numpy().imag.copy(), self.space)
[docs] def conj(self):
if np.isrealobj(self.to_numpy()):
return self.copy()
return NumpyVectorArray(np.conj(self.to_numpy()), self.space)
[docs] def __len__(self):
return self._len
[docs] def __getitem__(self, ind):
return NumpyVectorArrayView(self, ind)
[docs] def __delitem__(self, ind):
assert self.check_ind(ind)
if self._refcount[0] > 1:
self._deep_copy()
if type(ind) is slice:
ind = set(range(*ind.indices(self._len)))
elif not hasattr(ind, '__len__'):
ind = {ind if 0 <= ind else self._len + ind}
else:
l = self._len
ind = {i if 0 <= i else l+i for i in ind}
remaining = sorted(set(range(len(self))) - ind)
self._array = self._array[remaining]
self._len = len(self._array)
if not self._array.flags['OWNDATA']:
self._array = self._array.copy()
[docs] def copy(self, deep=False, *, _ind=None):
if _ind is None and not deep:
C = NumpyVectorArray(self._array, self.space)
C._len = self._len
C._refcount = self._refcount
self._refcount[0] += 1
return C
else:
new_array = self._array[:self._len] if _ind is None else self._array[_ind]
if not new_array.flags['OWNDATA']:
new_array = new_array.copy()
return NumpyVectorArray(new_array, self.space)
[docs] def append(self, other, remove_from_other=False):
assert self.dim == other.dim
assert not remove_from_other or (other is not self and getattr(other, 'base', None) is not self)
if self._refcount[0] > 1:
self._deep_copy()
other_array = other.to_numpy()
len_other = len(other_array)
if len_other == 0:
return
if len_other <= self._array.shape[0] - self._len:
if self._array.dtype != other_array.dtype:
self._array = self._array.astype(np.promote_types(self._array.dtype, other_array.dtype))
self._array[self._len:self._len + len_other] = other_array
else:
self._array = np.append(self._array[:self._len], other_array, axis=0)
self._len += len_other
if remove_from_other:
if other.is_view:
del other.base[other.ind]
else:
del other[:]
[docs] def scal(self, alpha, *, _ind=None):
if _ind is None:
_ind = slice(0, self._len)
assert isinstance(alpha, Number) \
or isinstance(alpha, np.ndarray) and alpha.shape == (self.len_ind(_ind),)
if self._refcount[0] > 1:
self._deep_copy()
if type(alpha) is np.ndarray:
alpha = alpha[:, np.newaxis]
alpha_type = type(alpha)
alpha_dtype = alpha.dtype if alpha_type is np.ndarray else alpha_type
if self._array.dtype != alpha_dtype:
self._array = self._array.astype(np.promote_types(self._array.dtype, alpha_dtype))
self._array[_ind] *= alpha
[docs] def axpy(self, alpha, x, *, _ind=None):
if _ind is None:
_ind = slice(0, self._len)
assert self.dim == x.dim
assert isinstance(alpha, Number) \
or isinstance(alpha, np.ndarray) and alpha.shape == (self.len_ind(_ind),)
if self._refcount[0] > 1:
self._deep_copy()
B = x.base._array[x.ind] if x.is_view else x._array[:x._len]
assert self.len_ind(_ind) == len(B) or len(B) == 1
alpha_type = type(alpha)
alpha_dtype = alpha.dtype if alpha_type is np.ndarray else alpha_type
if self._array.dtype != alpha_dtype or self._array.dtype != B.dtype:
dtype = np.promote_types(self._array.dtype, alpha_dtype)
dtype = np.promote_types(dtype, B.dtype)
self._array = self._array.astype(dtype)
if type(alpha) is np.ndarray:
alpha = alpha[:, np.newaxis]
self._array[_ind] += B * alpha
[docs] def dot(self, other, *, _ind=None):
if _ind is None:
_ind = slice(0, self._len)
assert self.dim == other.dim
A = self._array[_ind]
B = other.base._array[other.ind] if other.is_view else other._array[:other._len]
# .conj() is a no-op on non-complex data types
return A.conj().dot(B.T)
[docs] def pairwise_dot(self, other, *, _ind=None):
if _ind is None:
_ind = slice(0, self._len)
assert self.dim == other.dim
A = self._array[_ind]
B = other.base._array[other.ind] if other.is_view else other._array[:other._len]
assert len(A) == len(B)
# .conj() is a no-op on non-complex data types
return np.sum(A.conj() * B, axis=1)
[docs] def lincomb(self, coefficients, *, _ind=None):
if _ind is None:
_ind = slice(0, self._len)
assert 1 <= coefficients.ndim <= 2
if coefficients.ndim == 1:
coefficients = coefficients[np.newaxis, ...]
return NumpyVectorArray(coefficients.dot(self._array[_ind]), self.space)
[docs] def l1_norm(self, *, _ind=None):
if _ind is None:
_ind = slice(0, self._len)
return np.linalg.norm(self._array[_ind], ord=1, axis=1)
[docs] def l2_norm(self, *, _ind=None):
if _ind is None:
_ind = slice(0, self._len)
return np.linalg.norm(self._array[_ind], axis=1)
[docs] def l2_norm2(self, *, _ind=None):
if _ind is None:
_ind = slice(0, self._len)
A = self._array[_ind]
return np.sum((A * A.conj()).real, axis=1)
[docs] def sup_norm(self, *, _ind=None):
if self.dim == 0:
if _ind is None:
_ind = slice(0, self._len)
return np.zeros(self.len_ind(_ind))
else:
_, max_val = self.amax(_ind=_ind)
return max_val
[docs] def dofs(self, dof_indices, *, _ind=None):
if _ind is None:
_ind = slice(0, self._len)
assert isinstance(dof_indices, list) and (len(dof_indices) == 0 or min(dof_indices) >= 0) \
or (isinstance(dof_indices, np.ndarray) and dof_indices.ndim == 1
and (len(dof_indices) == 0 or np.min(dof_indices) >= 0))
# NumPy 1.9 is quite permissive when indexing arrays of size 0, so we have to add the
# following check:
assert self._len > 0 \
or (isinstance(dof_indices, list)
and (len(dof_indices) == 0 or max(dof_indices) < self.dim)) \
or (isinstance(dof_indices, np.ndarray) and dof_indices.ndim == 1
and (len(dof_indices) == 0 or np.max(dof_indices) < self.dim))
return self._array[:, dof_indices][_ind, :]
[docs] def amax(self, *, _ind=None):
if _ind is None:
_ind = slice(0, self._len)
assert self.dim > 0
if self._array.shape[1] == 0:
l = self.len_ind(_ind)
return np.ones(l) * -1, np.zeros(l)
A = np.abs(self._array[_ind])
max_ind = np.argmax(A, axis=1)
max_val = A[np.arange(len(A)), max_ind]
return max_ind, max_val
[docs] def __str__(self):
return self._array[:self._len].__str__()
def _format_repr(self, max_width, verbosity):
return super()._format_repr(max_width, verbosity, override={'array': str(self._array[:self._len].__str__())})
def __del__(self):
self._refcount[0] -= 1
def _deep_copy(self):
self._array = self._array.copy() # copy the array data
self._refcount[0] -= 1 # decrease refcount for original array
self._refcount = [1] # create new reference counter
def __add__(self, other):
if isinstance(other, Number):
assert other == 0
return self.copy()
assert self.dim == other.dim
return NumpyVectorArray(self._array[:self._len]
+ (other.base._array[other.ind] if other.is_view else other._array[:other._len]),
self.space)
def __iadd__(self, other):
assert self.dim == other.dim
if self._refcount[0] > 1:
self._deep_copy()
other_dtype = other.base._array.dtype if other.is_view else other._array.dtype
common_dtype = np.promote_types(self._array.dtype, other_dtype)
if self._array.dtype != common_dtype:
self._array = self._array.astype(common_dtype)
self._array[:self._len] += other.base._array[other.ind] if other.is_view else other._array[:other._len]
return self
__radd__ = __add__
def __sub__(self, other):
assert self.dim == other.dim
return NumpyVectorArray(self._array[:self._len]
- (other.base._array[other.ind] if other.is_view else other._array[:other._len]),
self.space)
def __isub__(self, other):
assert self.dim == other.dim
if self._refcount[0] > 1:
self._deep_copy()
other_dtype = other.base._array.dtype if other.is_view else other._array.dtype
common_dtype = np.promote_types(self._array.dtype, other_dtype)
if self._array.dtype != common_dtype:
self._array = self._array.astype(common_dtype)
self._array[:self._len] -= other.base._array[other.ind] if other.is_view else other._array[:other._len]
return self
def __mul__(self, other):
assert isinstance(other, Number) \
or isinstance(other, np.ndarray) and other.shape == (len(self),)
if type(other) is np.ndarray:
other = other[:, np.newaxis]
return NumpyVectorArray(self._array[:self._len] * other, self.space)
def __imul__(self, other):
assert isinstance(other, Number) \
or isinstance(other, np.ndarray) and other.shape == (len(self),)
if type(other) is np.ndarray:
other = other[:, np.newaxis]
if self._refcount[0] > 1:
self._deep_copy()
other_dtype = other.dtype if isinstance(other, np.ndarray) else type(other)
common_dtype = np.promote_types(self._array.dtype, other_dtype)
if self._array.dtype != common_dtype:
self._array = self._array.astype(common_dtype)
self._array[:self._len] *= other
return self
def __neg__(self):
return NumpyVectorArray(-self._array[:self._len], self.space)
[docs]class NumpyVectorSpace(VectorSpace):
"""|VectorSpace| of |NumpyVectorArrays|.
Parameters
----------
dim
The dimension of the vectors contained in the space.
id
See :attr:`~pymor.vectorarrays.interface.VectorSpace.id`.
"""
def __init__(self, dim, id=None):
self.dim = dim
self.id = id
[docs] def __eq__(self, other):
return type(other) is type(self) and self.dim == other.dim and self.id == other.id
[docs] def __hash__(self):
return hash(self.dim) + hash(self.id)
[docs] def zeros(self, count=1, reserve=0):
assert count >= 0
assert reserve >= 0
va = NumpyVectorArray(np.empty((0, 0)), self)
va._array = np.zeros((max(count, reserve), self.dim))
va._len = count
return va
[docs] def full(self, value, count=1, reserve=0):
assert count >= 0
assert reserve >= 0
va = NumpyVectorArray(np.empty((0, 0)), self)
va._array = np.full((max(count, reserve), self.dim), value)
va._len = count
return va
[docs] def random(self, count=1, distribution='uniform', random_state=None, seed=None, reserve=0, **kwargs):
assert count >= 0
assert reserve >= 0
assert random_state is None or seed is None
random_state = get_random_state(random_state, seed)
va = self.zeros(count, reserve)
va._array[:count] = _create_random_values((count, self.dim), distribution, random_state, **kwargs)
return va
@classinstancemethod
def make_array(cls, obj, id=None):
return cls._array_factory(obj, id=id)
@make_array.instancemethod
def make_array(self, obj):
return self._array_factory(obj, space=self)
@classinstancemethod
def from_numpy(cls, data, id=None, ensure_copy=False):
return cls._array_factory(data.copy() if ensure_copy else data, id=id)
@from_numpy.instancemethod
def from_numpy(self, data, ensure_copy=False):
return self._array_factory(data.copy() if ensure_copy else data, space=self)
@classinstancemethod
def from_file(cls, path, key=None, single_vector=False, transpose=False, id=None):
assert not (single_vector and transpose)
from pymor.tools.io import load_matrix
array = load_matrix(path, key=key)
assert isinstance(array, np.ndarray)
assert array.ndim <= 2
if array.ndim == 1:
array = array.reshape((1, -1))
if single_vector:
assert array.shape[0] == 1 or array.shape[1] == 1
array = array.reshape((1, -1))
if transpose:
array = array.T
return cls.make_array(array, id=id)
@from_file.instancemethod
def from_file(self, path, key=None, single_vector=False, transpose=False):
return type(self).from_file(path, key=key, single_vector=single_vector, transpose=transpose, id=self.id)
@classmethod
def _array_factory(cls, array, space=None, id=None):
if type(array) is np.ndarray:
pass
elif issparse(array):
array = array.toarray()
else:
array = np.array(array, ndmin=2)
if array.ndim != 2:
assert array.ndim == 1
array = np.reshape(array, (1, -1))
if space is None:
return NumpyVectorArray(array, cls(array.shape[1], id))
else:
assert array.shape[1] == space.dim
return NumpyVectorArray(array, space)
@property
def is_scalar(self):
return self.dim == 1 and self.id is None
[docs]class NumpyVectorArrayView(NumpyVectorArray):
is_view = True
def __init__(self, array, ind):
assert array.check_ind(ind)
self.base = array
self.ind = array.normalize_ind(ind)
self.space = array.space
[docs] def to_numpy(self, ensure_copy=False):
result = self.base.to_numpy()[self.ind]
if ensure_copy and not result.flags['OWNDATA']:
result = result.copy()
return result
[docs] def __len__(self):
return self.base.len_ind(self.ind)
[docs] def __getitem__(self, ind):
return self.base[self.base.sub_index(self.ind, ind)]
[docs] def __delitem__(self):
raise ValueError('Cannot remove from NumpyVectorArrayView')
[docs] def append(self, other, remove_from_other=False):
raise ValueError('Cannot append to NumpyVectorArrayView')
[docs] def copy(self, deep=False):
return self.base.copy(_ind=self.ind, deep=deep)
[docs] def scal(self, alpha):
assert self.base.check_ind_unique(self.ind)
self.base.scal(alpha, _ind=self.ind)
[docs] def axpy(self, alpha, x):
assert self.base.check_ind_unique(self.ind)
self.base.axpy(alpha, x, _ind=self.ind)
[docs] def dot(self, other):
return self.base.dot(other, _ind=self.ind)
[docs] def pairwise_dot(self, other):
return self.base.pairwise_dot(other, _ind=self.ind)
[docs] def lincomb(self, coefficients):
return self.base.lincomb(coefficients, _ind=self.ind)
[docs] def l1_norm(self):
return self.base.l1_norm(_ind=self.ind)
[docs] def l2_norm(self):
return self.base.l2_norm(_ind=self.ind)
[docs] def l2_norm2(self):
return self.base.l2_norm2(_ind=self.ind)
[docs] def sup_norm(self):
return self.base.sup_norm(_ind=self.ind)
[docs] def dofs(self, dof_indices):
return self.base.dofs(dof_indices, _ind=self.ind)
[docs] def amax(self):
return self.base.amax(_ind=self.ind)
def __add__(self, other):
if isinstance(other, Number):
assert other == 0
return self.copy()
assert self.dim == other.dim
return NumpyVectorArray(self.base._array[self.ind]
+ (other.base._array[other.ind] if other.is_view else other._array[:other._len]),
self.space)
def __iadd__(self, other):
assert self.dim == other.dim
assert self.base.check_ind_unique(self.ind)
if self.base._refcount[0] > 1:
self._deep_copy()
other_dtype = other.base._array.dtype if other.is_view else other._array.dtype
common_dtype = np.promote_types(self.base._array.dtype, other_dtype)
if self.base._array.dtype != common_dtype:
self.base._array = self.base._array.astype(common_dtype)
self.base.array[self.ind] += other.base._array[other.ind] if other.is_view else other._array[:other._len]
return self
__radd__ = __add__
def __sub__(self, other):
assert self.dim == other.dim
return NumpyVectorArray(self.base._array[self.ind]
- (other.base._array[other.ind] if other.is_view else other._array[:other._len]),
self.space)
def __isub__(self, other):
assert self.dim == other.dim
assert self.base.check_ind_unique(self.ind)
if self.base._refcount[0] > 1:
self._deep_copy()
other_dtype = other.base._array.dtype if other.is_view else other._array.dtype
common_dtype = np.promote_types(self.base._array.dtype, other_dtype)
if self.base._array.dtype != common_dtype:
self.base._array = self.base._array.astype(common_dtype)
self.base._array[self.ind] -= other.base._array[other.ind] if other.is_view else other._array[:other._len]
return self
def __mul__(self, other):
assert isinstance(other, Number) \
or isinstance(other, np.ndarray) and other.shape == (len(self),)
if type(other) is np.ndarray:
other = other[:, np.newaxis]
return NumpyVectorArray(self.base._array[self.ind] * other, self.space)
def __imul__(self, other):
assert isinstance(other, Number) \
or isinstance(other, np.ndarray) and other.shape == (len(self),)
if type(other) is np.ndarray:
other = other[:, np.newaxis]
assert self.base.check_ind_unique(self.ind)
if self.base._refcount[0] > 1:
self._deep_copy()
other_dtype = other.dtype if isinstance(other, np.ndarray) else type(other)
common_dtype = np.promote_types(self.base._array.dtype, other_dtype)
if self.base._array.dtype != common_dtype:
self.base._array = self.base._array.astype(common_dtype)
self.base._array[self.ind] *= other
return self
def __neg__(self):
return NumpyVectorArray(-self.base._array[self.ind], self.space)
def __del__(self):
return
[docs] def __repr__(self):
return f'NumpyVectorArrayView({self.to_numpy()}, {self.space})'