'DiscreteField' object has no attribute 'reshape'

[gatling-nrl]

This leads to confusing behavior in forms.

def some_functional():
    def _form(w):
        # arrays are shaped to group by facet, but we don't care for that here so reshape(2,-1)
        ii = w['ind_p0'].reshape(2,-1).T
        gg = grad(w['ind_p1']).reshape(2,-1).T
        pp = w.x.reshape(2,-1).T  
        nn = w.n.reshape(2,-1).T
        ...
        return w.x[0] # must return *something* of the right shape
    return skfem.Functional(_form)

skfem.asm(
    some_functional(),
    skfem.InteriorFacetBasis(mesh, skfem.ElementTriP1(), side=0),
    ind_p1=ind_p1,
    ind_p0=basis_p1.project(basis_p0.interpolate(ind_p0)),
)

gg, pp, nn all work, but ii returns 'DiscreteField' object has no attribute 'reshape'.

The "work around" or maybe the intended usage appears to be

ii = w['ind_p0'].value.reshape(2,-1).T

and this might be okay, except that

gg = grad(w['ind_p1']).reshape(2,-1).T

strongly implies (imo) that w['ind_p1'] is the "value". Moreover w[ind_p0]**2 works, and returns (w[ind_p0].value)**2 making me think maybe w[ind_p0] supports any ufunc.

I think reshape should work.

[kinnala]

You have it correct, the idea has been to have DiscreteField behave like DiscreteField.value, mainly for convenience. Would it be reasonable to simply add DiscreteField.reshape which returns reshape applied to the value?

[gatling-nrl]

That works, but there might be other ndarray methods that come up. Is duck typing DiscreteField so it acts like an ndarray the best implementation? Or should it also subclass ndarray. What is the difference between the return value from solve and a DiscreteField?

What's the difference between skfem.helpers.grad(w['thing']) and w['thing'].grad? Is w['thing'].grad a lazy evaluation (i.e. not created until it is used somewhere?

I think there is a common abstraction that could tie interpolate, project, solve, and Functional all together in a really intuitive way. Conceptually, aren't all these things discrete fields? If DiscreteField subclassed ndarray, could it keep its existing functionality, and add its basis as a property? Does that make it easier to fix project/interpolate?

[kinnala]

DiscreteField is a tuple of ndarrays so I'm not sure what you mean by subclassing ndarray. Everything is precomputed and DiscreteField.grad == grad(DiscreteField). DiscreteField represents a function and its derivatives evaluated at quadrature points and it does not necessarily have any connection to a Basis.

[kinnala]

We used to define forms like this:

def laplace(u, ux, uy, uz, v, vx, vy, vz):
     return ux * vx + uy * vy + uz * vz

So I invented DiscreteField to group all those starting with u into a single object.

[kinnala]

Then after making U = (u, ux, uy, uz) we had the form definition

def laplace(U, V):
    ...

but then I had to figure out how to access u so that it doesn't look like U[0]. So I made U a special tuple with U.toarray().

[gatling-nrl]

This might be some gap in my skfem / fem in general understanding then...

DiscreteField represents a function

That idea, of a function somehow represented in FEM space, I thought was very closely coupled to the idea of a Basis. Are you saying that the return value of solve very much has a connection to a Basis, but that DiscreteField is literally just collection of evaluations at fixed (xyz) tuples (the quadrature points)?

Or maybe another way of saying it is DiscreteField samples a function at discrete points and stores those samples in an ndarray whereas solve, basis.zeros, project all give you back a more abstract array containing (the coefficients of?) a function represented by the given Basis?

What is the intended way to move back and forth between DiscreteField and basis.zeros? Does basis.zeros have a name? Is it a Projection? (I know that class doesn't actually exist in skfem, just conceptually.)

[gatling-nrl]

DiscreteField is a tuple of ndarrays so I'm not sure what you mean by subclassing ndarray.

I meant that when you do

x = DiscreteField()

what you get back is a child of ndarray. The .value property goes away... x literally is the .value. The other properties, e.g. .grad remain basically as is. Most of the methods on the class would go away because they are already implemented in numpy using this approach.

Essentially DiscreteField is a ndarray with some extra data attached. I think this was essentially what you built using the NamedTuple and I think what I am suggesting is transparent to user code, so long as we support backward compatibility with a @property to return self when value is accessed.

This issue is very similar to #846 and perhaps the two should be merged.

[kinnala]

Ok, now I see what you mean. I'll experiment with your idea of subclassing ndarray to see if it works. It might even simplify the implementation of DiscreteField.

[gatling-nrl]

I'm doing that now, I'll have a PR for you in a few.

[gatling-nrl]

scikit-fem/skfem/element/discrete_field.py

Lines 81 to 89 in 84ef5ad

	def zeros_like(self) -> 'DiscreteField':
	"""Return zero :class:`~skfem.element.DiscreteField` with same size."""
	def zero_or_none(x):
	if x is None:
	return None
	return np.zeros_like(x)
	return DiscreteField(*[zero_or_none(field) for field in self])

Is not exactly how zeros_like works. I think this might be a break in backward compatibility, but... I think you shouldn't override zeros_like and should accept the standard numpy behavior. Thoughts?

[gatling-nrl]

scikit-fem/skfem/element/discrete_field.py

Lines 71 to 75 in 84ef5ad

	def _split(self):
	"""Split all components based on their first dimension."""
	return [DiscreteField(*[f[i] for f in self if f is not None])
	for i in range(self.value.shape[0])]

What's the goal of this? Is the idea when you slice a DiscreteField you want to slice all of the derived arrays as well? Basically that

x = DiscreteField(...)
x[10:20,:,:-5].grad

should work as expected?

[kinnala]

No no, there is something known as ElementComposite which will complicate matters. I didn't remember split is done that way. I think Basis.interpolate also assumes DiscreteField being iterable.

[gatling-nrl]

So

for a in DiscreteField():
    ...

should iterate over the properties in the order they're listed in the NamedTuple? Which would obviously be painfully confusing :(

Because anyone would be expecting that to iterate over the elements of .value

[kinnala]

Unless you want to redo Basis.interpolate and all those split thingies.

[gatling-nrl]

It certainly seems like the way to go... but I think DiscreteField(np.ndarray) is going to be a breaking change for this reason. Because its unexpected to the extreme for an ndarray to iterate over something other than its values. But that is the current implementation.

I guess its somewhat semantics, but you were trying to make w['thing'] act like w['thing'].value (and I like that abstraction a lot) but in the current implementation .reshape didn't work and if I had tried for a in w['thing'] I would have been even more confused. So I think the current implementation is a little too inconsistent and something should be done?

[kinnala]

I can try it out. This is very core part of the lib so it requires some careful consideration to not break too many things.

[gdmcbain]

Possible alternative to deriving DiscreteField from ndarray: define its __getattr__ method so that when someone tries to call a method that isn't thete, e.g. reshape, it returns the result of calling that merhod on the .value attribute. This is a popular pattern in OOP: delegation.

[gatling-nrl]

@gdmcbain delegation isn't going to fix what I think is the more fundamental issue. If f is a DiscreteField then

f**2 # works, f seems like an ndarray
f + f # works, f seems like an ndarray 
[a for a in f] # !! f doesn't act like ndarray at all !! 

Whether inherit or delegate, I think f should adopt a single abstraction, and my vote is that f always act like the ndarray f.value. I say this not arbitrarily, but because that's how w.x and w.n work, so it seems like a consistent pattern if w.my_df behaves the same way. And certainly

for df, x in zip(w.my_df, w.x):

shouldn't trick me the way it does now.

Re: delegate vs inherit, I have a personal preference for inherit, but I also think this is a semantic/syntactic debate. One of the trickier parts of implementation I think is supporting something like f[1:2].grad or perhaps f.reshape().grad. And this will eventually come up because someone will pass a reshaped or sliced DiscreteField in to some function that accesses .grad and there's no real reason that can't work. Maybe choose delegate vs inherit on whichever makes that easier to implement.

[gatling-nrl]

Also, after delegation, .zeros_like and .ones_like do two different things under certain circumstances. I'm not sure what other sneaky bugs are lurking here.

[gdmcbain]

Yes I see. Anyway i've published a draft branch in #847 to experiment with in case. It's only two lines and probably backwards compatible.

[gatling-nrl]

Egads, w.my_df and w['my_df'] aren't the same thing!!

[gatling-nrl]

After #847 len(w['my_df']) and w['my_df'].shape are confusing since len() returns length of the tuple and .shape gets delegated. I think ideally len() would act on .value and return whatever.

[kinnala]

Then again many examples given in this thread are really confusing to me. Why would you want to reshape a DiscreteField? Why would you want to know it's length?

It's really about defining forms conveniently and if you don't need those things inside a form definition then we can forget about it. No need to make such drastic changes for no benefit.

[kinnala]

Re: delegate vs inherit, I have a personal preference for inherit, but I also think this is a semantic/syntactic debate. One of the trickier parts of implementation I think is supporting something like f[1:2].grad or perhaps f.reshape().grad. And this will eventually come up because someone will pass a reshaped or sliced DiscreteField in to some function that accesses .grad and there's no real reason that can't work. Maybe choose delegate vs inherit on whichever makes that easier to implement.

Slicing and then taking grad cannot work. In general, value and grad have different shapes. So after slicing grad would need to be dropped.

[kinnala]

My try on this: #849

[kinnala]

I went as far as I could in #849. Waiting for you to try out the branch and give feedback. The remaining irritation regarding slicing is explained in the PR.

[gatling-nrl]

Then again many examples given in this thread are really confusing to me. Why would you want to know it's length?

I admit I got a little distracted by the computer sciency aspects yesterday. I don't have good use cases for the examples I gave. I was more pointing out how the DiscreteField was not behaving like a consistent abstraction.

Why would you want to reshape a DiscreteField?

Well, that's how I got here. I was still abusing the Functional to make interesting plots while I worked on making a "current source" boundary condition in skfem. My main argument of these examples is that if you're going to tamper with DiscreteField, might as well make it as well behaved as possible, which I think means making it act as much like an ordinary ndarray as possible.

[kinnala]

My main argument of these examples is that if you're going to tamper with DiscreteField, might as well make it as well behaved as possible, which I think means making it act as much like an ordinary ndarray as possible.

Fair enough. I think #849 somehow makes sense and seems to solve some minor pains we had with the existing implementation such as

• inability to do w.fun.grad
• sometimes in complicated forms I still had to go for u.value for some reasons I don't fully understand

If you agree it is a reasonable approach I think I will continue by

• making a nice __repr__ for DiscreteField
• fixing the constructor as pointed out in your review comment
• fixing .value to return ndarray and not DiscreteField
• removing uses of .value under skfem
• adding a separate test module for DiscreteField
• thinking about the slicing; can we somehow consistently slice over all fields at once? it would somehow make sense, e.g., to slice over -2 dimension to restrict to a subset of elements.

[gatling-nrl]

I like your plan. Here are my thoughts.

fixing .value to return ndarray and not DiscreteField

I looked around at what other packages do. Both pandas and xarray use .to_numpy() to get from their respective data types to a bare array. I think that's what skfem should do too, and .value should just be a deprecated way to call .to_numpy(). I think pandas in particular is popular enough that it can form a sort of standard calling convention.

removing uses of .value under skfem
adding a separate test module for DiscreteField
I have the tendency of leaving around deprecated stuff so that I remember the state of affairs and support those as long as it's reasonable.

Test case for .value! :)

making a nice repr for DiscreteField

pandas and xarray and even numpy itself point the way. Let the repr show, for each property that isn't None, a subset of the values with ... in between. Dont bother including anything about the attributes that are None.

thinking about the slicing;

One brute force approach I had is that then you slice, re-compute the derived attributes. This would probably work a lot better if things were lazy evaluated, but maybe it is still worth considering some version of this idea.

---

While we're talking about all this, what are your thoughts on adding class Projection(np.ndarray). It would be returned by Basis.zeros and Basis.project. (and others maybe?) It would be very analogous to what you're doing with DiscreteField. Its extra attribute would be its basis. This is very nearly the bare ndarray like you want, but facilitates all sorts of automated stuff (including helpful exceptions) and makes it clear in user code when a particular array is actually a projection.

[gatling-nrl]

I forgot to add, since the projection is currently an ndarray, making class Projection(np.ndarray) should be 100% backward compatible, fully compatible with scipy, and current skfem users could just ignore the extra properties on the object. At least, this is how I was picturing the idea I was trying to describe.

[kinnala]

pandas and xarray and even numpy itself point the way. Let the repr show, for each property that isn't None, a subset of the values with ... in between. Dont bother including anything about the attributes that are None.

But there is something important about DiscreteField that I want to message.

1. The last dimension of the numpy array is always over elementwise quadrature points
2. The second to last dimension of the numpy array is always over elements
3. The remaining dimensions, if they exist, are related to the spatial dimension of the problem

So somehow I want to transfer this information through __repr__.

thinking about the slicing;

One brute force approach I had is that then you slice, re-compute the derived attributes. This would probably work a lot better if things were lazy evaluated, but maybe it is still worth considering some version of this idea.

This is not exactly how you can think about it. Even though mathematically grad is derived from value, it is not the case here. Both are separately defined and simultaneously computed inside Element classes for the given basis functions. This is mainly due to performance.

Splitting evaluation of value and grad would be technically possible but such a large refactoring that I'm not up for it any time soon. :-)

[gatling-nrl]

x = skfem.DiscreteField(np.arange(10*11*12*3).reshape((10,11,12,3)))
print(x)

DiscreteField([[[[0, 1, 2, 3, ..., 17157, 17158, 17159]]]])

Quadrature points per element: 3
Elements: 12
Spatial Dimensions: (10, 11)
Available Attributes: grad, hess, grad3

print(x.shape)

(10, 11, 12, 3)

Ok so this is contrived and I don't have any code. But it's pretty :)

[kinnala]

I forgot to add, since the projection is currently an ndarray, making class Projection(np.ndarray) should be 100% backward compatible, fully compatible with scipy, and current skfem users could just ignore the extra properties on the object. At least, this is how I was picturing the idea I was trying to describe.

I still need to think about it but in other codes this seems to be typically called GridFunction:

• MFEM: https://mfem.github.io/doxygen/html/classmfem_1_1GridFunction.html
• NGSolve: https://ngsolve.org/forum/ngspy-forum/1456-gridfunction-to-numpy-array-to-new-gridfunction
• Bem++: http://bempp.com/handbook/api/grid_functions.html

[kinnala]

I forgot to add, since the projection is currently an ndarray, making class Projection(np.ndarray) should be 100% backward compatible, fully compatible with scipy, and current skfem users could just ignore the extra properties on the object. At least, this is how I was picturing the idea I was trying to describe.

How would this work in solve? Because those matrices and load vectors don't have any info about Basis. Should they also be subclasses?

[gatling-nrl]

Re: slicing... what if you didn't actually slice anything. Suppose you just recorded the start, stop, and stride.

Then when you do

a = DiscreteField(...)
s = slice(...)
b = a[s]

b is ... I'm hand waving a bit here, but b is a "view" into a that saved the slice info as a property. Then when you do

b.grad

you can access that slice info to properly slice grad, hess, or whatever, in a way that make sense in the context. This would require changing .grad and friends to be @property though so you can make functions to work with that slice property and the underlying ._grad ndarray.

[kinnala]

I forgot to add, since the projection is currently an ndarray, making class Projection(np.ndarray) should be 100% backward compatible, fully compatible with scipy, and current skfem users could just ignore the extra properties on the object. At least, this is how I was picturing the idea I was trying to describe.

How would this work in solve? Because those matrices and load vectors don't have any info about Basis. Should they also be subclasses?

Since this library is really about "finite element assembly" I sometimes even think if having solve is really too much. There are plenty of other libraries for that purpose. That's why it's inside utils.

[gatling-nrl]

In Bempp, data on a given grid is represented as a grid function object. A grid function consists of a set of basis function coefficients and a corresponding function space.

This does sound exactly like what we're talking about. However, they subclassed object

https://github.com/bempp/bempp-cl/blob/0eb59db86baecb16239337200a489333710b0bbe/bempp/api/assembly/grid_function.py#L107

You said you like the way skfem doesn't use such opaque objects, relying instead on ndarrays and scipy. I agree with you about that.

Since this library is really about "finite element assembly" I sometimes even think if having solve is really too much. There are plenty of other libraries for that purpose. That's why it's inside utils.

I agree with that too. I am happy... more happy, in fact, to import scipy and use that solver. Keeping skfem data objects compatible with scipy makes them easier to understand.

I don't think solve (the utils one or the scipy one) needs to know about the extra functionality of Projection (or GridFunction if you want). Since it subclassed ndarray, if it is coded properly it will act exactly like an ndarray. The gain comes in other places. Particularly, when you get an incompatible grid function, you can say A LOT more about it than IndexError or some generic numpy nonsense. I think it also allows syntax like

x_in_p1 = x_in_p0.project(basis_p1)
x_discrete_field = x_in_p0.interpolate(basis_p1.quadrature_points)

which feels intuitive and readable to me.