/
interfaces.py
634 lines (508 loc) · 20.3 KB
/
interfaces.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
# -*- coding: utf-8 -*-
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2017 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
from numbers import Number
import sys
from packaging.version import Version
import numpy as np
from pymor.core.config import config
from pymor.core.interfaces import BasicInterface, ImmutableInterface, abstractmethod
_INDEXTYPES = (Number,) if Version(np.__version__) >= Version('1.9') else (Number, np.intp)
class VectorArrayInterface(BasicInterface):
"""Interface for vector arrays.
A vector array should be thought of as a list of (possibly high-dimensional) vectors.
While the vectors themselves will be inaccessible in general (e.g. because they are
managed by an external PDE solver code), operations on the vectors like addition can
be performed via this interface.
It is assumed that the number of vectors is small enough such that scalar data
associated to each vector can be handled on the Python side. As such, methods like
:meth:`~VectorArrayInterface.l2_norm` or :meth:`~VectorArrayInterface.gramian` will
always return |NumPy arrays|.
An implementation of the `VectorArrayInterface` via |NumPy arrays| is given by
|NumpyVectorArray|. In general, it is the implementors decision how memory is
allocated internally (e.g. continuous block of memory vs. list of pointers to the
individual vectors.) Thus, no general assumptions can be made on the costs of operations
like appending to or removing vectors from the array. As a hint for 'continuous block
of memory' implementations, :meth:`~VectorSpaceInterface.zeros` provides a `reserve`
keyword argument which allows to specify to what size the array is assumed to grow.
As with |Numpy array|, |VectorArrays| can be indexed with numbers, slices and
lists or one-dimensional |NumPy arrays|. Indexing will always return a new
|VectorArray| which acts as a view into the original data. Thus, if the indexed
array is modified via :meth:`~VectorArrayInterface.scal` or :meth:`~VectorArrayInterface.axpy`,
the vectors in the original array will be changed. Indices may be negative, in
which case the vector is selected by counting from the end of the array. Moreover
indices can be repeated, in which case the corresponding vector is selected several
times. The resulting view will be immutable, however.
.. note::
It is disallowed to append vectors to a |VectorArray| view or to remove
vectors from it. Removing vectors from an array with existing views
will lead to undefined behavior of these views. As such, it is generally
advisable to make a :meth:`~VectorArrayInterface.copy` of a view for long
term storage. Since :meth:`~VectorArrayInterface.copy` has copy-on-write
semantics, this will usually cause little overhead.
Attributes
----------
data
Implementors can provide a `data` property which returns a |NumPy array| of
shape `(len(v), v.dim)` containing the data stored in the array. Access should
be assumed to be slow and is mainly intended for debugging / visualization
purposes or to once transfer data to pyMOR and further process it using NumPy.
In the case of |NumpyVectorArray|, an actual view of the internally used
|NumPy array| is returned, so changing it, will alter the |VectorArray|.
Thus, you cannot assume to own the data returned to you, in general.
dim
The dimension of the vectors in the array.
is_view
`True` if the array is a view obtained by indexing another array.
space
The |VectorSpace| the array belongs to.
"""
is_view = False
def zeros(self, count=1, reserve=0):
"""Create a |VectorArray| of null vectors of the same |VectorSpace|.
This is a shorthand for `self.space.zeros(count, reserve)`.
Parameters
----------
count
The number of vectors.
reserve
Hint for the backend to which length the array will grow.
Returns
-------
A |VectorArray| containing `count` vectors whith each component
zero.
"""
return self.space.zeros(count, reserve=reserve)
def empty(self, reserve=0):
"""Create an empty |VectorArray| of the same |VectorSpace|.
This is a shorthand for `self.space.zeros(0, reserve)`.
Parameters
----------
reserve
Hint for the backend to which length the array will grow.
Returns
-------
An empty |VectorArray|.
"""
return self.space.zeros(0, reserve=reserve)
@property
def dim(self):
return self.space.dim
@abstractmethod
def __len__(self):
"""The number of vectors in the array."""
pass
@abstractmethod
def __getitem__(self, ind):
"""Return a |VectorArray| view onto a subset of the vectors in the array."""
pass
@abstractmethod
def __delitem__(self, ind):
"""Remove vectors from the array."""
pass
@abstractmethod
def append(self, other, remove_from_other=False):
"""Append vectors to the array.
Parameters
----------
other
A |VectorArray| containing the vectors to be appended.
remove_from_other
If `True`, the appended vectors are removed from `other`.
For list-like implementations this can be used to prevent
unnecessary copies of the involved vectors.
"""
pass
@abstractmethod
def copy(self, deep=False):
"""Returns a copy of a subarray.
All |VectorArray| implementations in pyMOR have copy-on-write semantics:
if not specified otherwise by setting `deep` to `True`, the returned
copy will hold a handle to the same array data as the original array,
and a deep copy of the data will only be performed when one of the arrays
is modified.
Note that for |NumpyVectorArray|, a deep copy is always performed when only
some vectors in the array are copied.
Parameters
----------
deep
Ensure that an actual copy of the array data is made (see above).
Returns
-------
A copy of the |VectorArray|.
"""
pass
@abstractmethod
def scal(self, alpha):
"""BLAS SCAL operation (in-place scalar multiplication).
This method calculates ::
self = alpha*self
If `alpha` is a scalar, each vector is multiplied by this scalar. Otherwise, `alpha`
has to be a one-dimensional |NumPy array| of the same length as `self`
containing the factors for each vector.
Parameters
----------
alpha
The scalar coefficient or one-dimensional |NumPy array| of coefficients
with which the vectors in `self` are multiplied.
"""
pass
@abstractmethod
def axpy(self, alpha, x):
"""BLAS AXPY operation.
This method forms the sum ::
self = alpha*x + self
If the length of `x` is 1, the same `x` vector is used for all vectors
in `self`. Otherwise, the lengths of `self` and `x` have to agree.
If `alpha` is a scalar, each `x` vector is multiplied with the same factor `alpha`.
Otherwise, `alpha` has to be a one-dimensional |NumPy array| of the same length as
`self` containing the coefficients for each `x` vector.
Parameters
----------
alpha
The scalar coefficient or one-dimensional |NumPy array| of coefficients with which
the vectors in `x` are multiplied.
x
A |VectorArray| containing the x-summands.
"""
pass
@abstractmethod
def dot(self, other):
"""Returns the inner products between |VectorArray| elements.
Parameters
----------
other
A |VectorArray| containing the second factors.
Returns
-------
A |NumPy array| `result` such that:
result[i, j] = ( self[i], other[j] ).
"""
pass
def inner(self, other, product=None):
"""Inner products w.r.t. given product |Operator|.
Equivalent to `self.dot(other)` if `product` is None,
else equivalent to `product.apply2(self, other)`.
"""
if product is None:
return self.dot(other)
else:
return product.apply2(self, other)
@abstractmethod
def pairwise_dot(self, other):
"""Returns the pairwise inner products between |VectorArray| elements.
Parameters
----------
other
A |VectorArray| containing the second factors.
Returns
-------
A |NumPy array| `result` such that:
result[i] = ( self[i], other[i] ).
"""
pass
def pairwise_inner(self, other, product=None):
"""Pairwise inner products w.r.t. given product |Operator|.
Equivalent to `self.pairwise_dot(other)` if `product` is None,
else equivalent to `product.pairwise_apply2(self, other)`.
"""
if product is None:
return self.pairwise_dot(other)
else:
return product.pairwise_apply2(self, other)
@abstractmethod
def lincomb(self, coefficients):
"""Returns linear combinations of the vectors contained in the array.
Parameters
----------
coefficients
A |NumPy array| of dimension 1 or 2 containing the linear
coefficients. `coefficients.shape[-1]` has to agree with
`len(self)`.
Returns
-------
A |VectorArray| `result` such that:
result[i] = ∑ self[j] * coefficients[i,j]
in case `coefficients` is of dimension 2, otherwise
`len(result) == 1` and
result[0] = ∑ self[j] * coefficients[j].
"""
pass
def norm(self, product=None):
"""Norm w.r.t. given inner product |Operator|.
Equivalent to `self.l2_norm()` if `product` is None,
else equivalent to `np.sqrt(product.pairwise_apply2(self, self))`.
"""
if product is None:
return self.l2_norm()
else:
return np.sqrt(product.pairwise_apply2(self, self))
def norm2(self, product=None):
"""Squared norm w.r.t. given inner product |Operator|.
Equivalent to `self.l2_norm2()` if `product` is None,
else equivalent to `product.pairwise_apply2(self, self)`.
"""
if product is None:
return self.l2_norm2()
else:
return product.pairwise_apply2(self, self)
@abstractmethod
def l1_norm(self):
"""The l1-norms of the vectors contained in the array.
Returns
-------
A |NumPy array| `result` such that `result[i]` contains the norm
of `self[i]`.
"""
pass
@abstractmethod
def l2_norm(self):
"""The l2-norms of the vectors contained in the array.
Returns
-------
A |NumPy array| `result` such that `result[i]` contains the norm
of `self[i]`.
"""
pass
@abstractmethod
def l2_norm2(self):
"""The squared l2-norms of the vectors contained in the array.
Returns
-------
A |NumPy array| `result` such that `result[i]` contains the norm
of `self[i]`.
"""
pass
def sup_norm(self):
"""The l-infinity--norms of the vectors contained in the array.
Returns
-------
A |NumPy array| `result` such that `result[i]` contains the norm
of `self[i]`.
"""
if self.dim == 0:
return np.zeros(len(self))
else:
_, max_val = self.amax()
return max_val
@abstractmethod
def components(self, component_indices):
"""Extract components of the vectors contained in the array.
Parameters
----------
component_indices
List or 1D |NumPy array| of indices of the vector components that are to
be returned.
Returns
-------
A |NumPy array| `result` such that `result[i, j]` is the `component_indices[j]`-th
component of the `i`-th vector of the array.
"""
pass
@abstractmethod
def amax(self):
"""The maximum absolute value of the vectors contained in the array.
Returns
-------
max_ind
|NumPy array| containing for each vector an index at which the maximum is
attained.
max_val
|NumPy array| containing for each vector the maximum absolute value of its
components.
"""
pass
def gramian(self, product=None):
"""Shorthand for `self.inner(self, product)`."""
return self.inner(self, product)
def __add__(self, other):
"""The pairwise sum of two |VectorArrays|."""
if isinstance(other, Number):
assert other == 0
return self.copy()
result = self.copy()
result.axpy(1, other)
return result
def __iadd__(self, other):
"""In-place pairwise addition of |VectorArrays|."""
self.axpy(1, other)
return self
__radd__ = __add__
def __sub__(self, other):
"""The pairwise difference of two |VectorArrays|."""
result = self.copy()
result.axpy(-1, other)
return result
def __isub__(self, other):
"""In-place pairwise difference of |VectorArrays|."""
self.axpy(-1, other)
return self
def __mul__(self, other):
"""Product by a scalar."""
result = self.copy()
result.scal(other)
return result
__rmul__ = __mul__
def __imul__(self, other):
"""In-place product by a scalar."""
self.scal(other)
return self
def __neg__(self):
"""Product by -1."""
result = self.copy()
result.scal(-1)
return result
def check_ind(self, ind):
"""Check if `ind` is an admissable list of indices in the sense of the class documentation."""
l = len(self)
return (type(ind) is slice or
isinstance(ind, _INDEXTYPES) and -l <= ind < l or
isinstance(ind, (list, np.ndarray)) and all(-l <= i < l for i in ind))
def check_ind_unique(self, ind):
"""Check if `ind` is an admissable list of non-repeated indices in the sense of the class documentation."""
l = len(self)
return (type(ind) is slice or
isinstance(ind, _INDEXTYPES) and -l <= ind < l or
isinstance(ind, (list, np.ndarray)) and len(set(i if i >= 0 else l+i for i in ind if -l <= i < l)) == len(ind))
def len_ind(self, ind):
"""Return the number of given indices."""
l = len(self)
return (len(range(*ind.indices(l))) if type(ind) is slice else
1 if not hasattr(ind, '__len__') else
len(ind))
def len_ind_unique(self, ind):
"""Return the number of specified unique indices."""
l = len(self)
return (len(range(*ind.indices(l))) if type(ind) is slice else
1 if isinstance(ind, _INDEXTYPES) else
len(set(i if i >= 0 else l+i for i in ind)))
def normalize_ind(self, ind):
"""Normalize given indices such that they are independent of the array length."""
if type(ind) is slice:
return slice(*ind.indices(len(self)))
elif not hasattr(ind, '__len__'):
ind = ind if 0 <= ind else len(self)+ind
return slice(ind, ind+1)
else:
l = len(self)
return [i if 0 <= i else l+i for i in ind]
def sub_index(self, ind, ind_ind):
"""Return indices corresponding to the view `self[ind][ind_ind]`"""
if type(ind) is slice:
ind = range(*ind.indices(len(self)))
if type(ind_ind) is slice:
if config.PY2: # make ind a list since is xrange is not sliceable
return list(ind)[ind_ind]
else:
result = ind[ind_ind]
return slice(result.start, result.stop, result.step)
elif hasattr(ind_ind, '__len__'):
return [ind[i] for i in ind_ind]
else:
return [ind[ind_ind]]
else:
if not hasattr(ind, '__len__'):
ind = [ind]
if type(ind_ind) is slice:
return ind[ind_ind]
elif hasattr(ind_ind, '__len__'):
return [ind[i] for i in ind_ind]
else:
return [ind[ind_ind]]
class VectorSpaceInterface(ImmutableInterface):
"""Class describing a vector space.
Vector spaces act as factories for |VectorArrays| of vectors
contained in them. As such, they hold all data necessary to
create |VectorArrays| of a given type (e.g. the dimension of
the vectors, or a socket for communication with an external
PDE solver).
New |VectorArrays| of null vectors are created via
:meth:`~VectorSpaceInterface.zeros`. The
:meth:`~VectorSpaceInterface.make_array` method builds a new
|VectorArray| from given raw data of the underlying linear algebra
backend (e.g. a |Numpy array| in the case of |NumpyVectorSpace|).
Some vector spaces can create new |VectorArrays| from a given
|Numpy array| via the :meth:`~VectorSpaceInterface.from_data`
method.
Each vector space has a string :attr:`~VectorSpaceInterface.id`
to distinguish mathematically different spaces appearing
in the formulation of a given problem.
Vector spaces can be compared for equality via the `==` and `!=`
operators. To test if a given |VectorArray| is an element of
the space, the `in` operator can be used.
Attributes
----------
id
None, or a string describing the mathematical identity
of the vector space (for instance to distinguish different
components in an equation system).
dim
The dimension (number of degrees of freedom) of the
vectors contained in the space.
is_scalar
Equivalent to `isinstance(space, NumpyVectorSpace) and space.dim == 1`.
"""
id = None
dim = None
is_scalar = False
@abstractmethod
def make_array(*args, **kwargs):
"""Create a |VectorArray| from raw data.
This method is used in the implementation of |Operators|
and |Discretizations| to create new |VectorArrays| from
raw data of the underlying solver backends. The ownership
of the data is transferred to the newly created array.
The exact signature of this method depends on the wrapped
solver backend.
"""
pass
@abstractmethod
def zeros(self, count=1, reserve=0):
"""Create a |VectorArray| of null vectors
Parameters
----------
count
The number of vectors.
reserve
Hint for the backend to which length the array will grow.
Returns
-------
A |VectorArray| containing `count` vectors with each component zero.
"""
pass
def empty(self, reserve=0):
"""Create an empty |VectorArray|
This is a shorthand for `self.zeros(0, reserve)`.
Parameters
----------
reserve
Hint for the backend to which length the array will grow.
Returns
-------
An empty |VectorArray|.
"""
return self.zeros(0, reserve=reserve)
def from_data(self, data):
"""Create a |VectorArray| from a |NumPy array|
Note that this method will not be supported by all vector
space implementations.
Parameters
----------
data
|NumPy| array.
Returns
-------
A |VectorArray| with `data` as data.
"""
raise NotImplementedError
def __eq__(self, other):
return other is self
def __ne__(self, other):
return not (self == other)
def __contains__(self, other):
return self == getattr(other, 'space', None)
def __hash__(self):
return hash(self.id)
def __repr__(self):
return '{}({})'.format(self.__class__.__name__, self.id)