Vector and tensor fields

[kohr-h]

Suggestion:
Since vector fields, or more general, tensor fields, are always power spaces of an underlying scalar function space, we can safely assume that it is homogeneous in data type. So we could use Numpy's advanced data type functionality to represent elements of such a space:

>>> dt = np.dtype((float, (100, 100)))
>>> dt
dtype(('<f8', (100, 100)))
>>> a = np.ones((2, 3), dtype=dt)
>>> a.shape
(2, 3, 100, 100)

The only problem I see is when we get non-contiguous input data, i.e. one component at a time from some location we don't have control of. The simple way out here would be to just make it contiguous if necessary - which probably means always if the data is not one big array.

Out of curiosity: Is this feature supported in gpuarray?

[kohr-h]

Okay, now that we have a plan for this functionality (see #731) how are we going to go about it? Somehow both this issue and #731 have to be solved simultaneously since they both implement multi-indexed spaces, this one on the CPU, the other one on the GPU.

That change is going to be massive, a bloody hell to review as a whole. Is there any way to make this transition more smooth and digestible? To me, there seem to be at least two options:

1. Change everything at once. That would mean a single point of transition from purely linear data containers to multi-dimensional ones. Somehow the ideal scenario, but as I wrote, it's going to be quite an amount of code at once.
2. Introduce tensors as add-on and make the full transition with the GPU arrays. The disadvantage of this is that both variants (tensors and ntuples) co-exist for a while, which can be confusing. But since it's not the most public part of the library, that shouldn't be that bad.
3. Replace only the CPU ntuples with tensors and keep the downstream classes (like DiscreteLp) "in the dark" about the change. They still see storage as linear. The convenience functions stay (which they will anyway), and the tensor classes print themselves exactly the same way as the NtuplesBase classes do now. The additional capabilities are "unlocked" when the GPU stuff goes in.

In any case, I will put up "WIP" PRs for both branches, and I think the best is to have some kind of continuous development, no big changes in one chunk. Preliminary design reviews would make much sense.

[adler-j]

If at all possible I'd go for option 1 since any of the other ones would perhaps imply we need to make design choices that will be sub-optimal down the line.

I'll have a look at the WIP stuff and review.

[adler-j]

Out of curiosity: Is this feature supported in gpuarray?

Absolutely no idea, but I think the idea is that it should. Anyway dtypes like this are only "hacks" in a sense, and don't "really" exist, see for example

>>> arr = np.zeros(5, dtype=(float, (3, 3)))
>>> arr2 = np.zeros([5, 3, 3], dtype=float)
>>> arr.dtype == arr2.dtype
True
>>> arr.shape == arr2.shape
True

So an equivalent way to approach it is to simply give a shape.

[adler-j]

The only problem I see is when we get non-contiguous input data, i.e. one component at a time from some location we don't have control of.

This should i.m.o. (if the user wants to explicitly have non-contiguous data) be handled with a productspace.

[kohr-h]

Absolutely no idea, but I think the idea is that it should. Anyway dtypes like this are only "hacks" in a sense, and don't "really" exist, see for example

arr = np.zeros(5, dtype=(float, (3, 3)))
arr2 = np.zeros([5, 3, 3], dtype=float)
arr.dtype == arr2.dtype
True
arr.shape == arr2.shape
True

So an equivalent way to approach it is to simply give a shape.

Totally agree. I now have a much better picture of what you can do in libgpuarray. Using dtypes with shapes is only a hack and not necessary. The space should have a shape.

The only problem I see is when we get non-contiguous input data, i.e. one component at a time from some location we don't have control of.

This should i.m.o. (if the user wants to explicitly have non-contiguous data) be handled with a productspace.

That may not even be necessary, but I'll investigate. A situation where this could occur would be slicing with [::2] or similar. The question is whether you can still use those arrays with kernels in an efficient way or not.

[kohr-h]

Closed by #1088