Add capability for naming fields of q and aux

[ketch]

Example usage:

state.q_names = ['density','momentum','energy']

state.q.density[i,j] # gives state.q[0,i,j] state.q.momentum[i,j] # gives state.q[1,i,j]

[mandli]

ndarrays allow field names to be built into them and referenced although I have never used that particular functionality. I might be a nice way to implement this rather than overriding the setter and getters.

[ketch]

Thanks for the suggestion, Kyle. It seems easy to do:

http://docs.scipy.org/doc/numpy/user/basics.rec.html

http://www.scipy.org/Cookbook/Recarray

It seems there are two approaches. They mention that one of them may incur overhead, but say nothing about the performance of the other. Although my guess is that it won't be significant as long as it doesn't cause arrays to be copied when passing to Fortran (Amal has just eliminated that issue in the current code). We would need to try it and check timings.

[mandli]

So I have tried implementing this using the first methodology above but have run into some issues that may not allow this to really happen (at least via this approach). So the following works:

>>> q = numpy.empty( (10,10), dtype={'names':['density', 'momentum', 'energy'], 'formats':[float64, float64, float64]})

but when asked its shape says:

>>> q.shape
(10, 10)
>>> q['density'].shape
(10, 10)

and unfortunately

>>> q[0,:,:]
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-56-de2ba28b1d31> in <module>()
----> 1 q[0,:,:]

Subclassing ndarray is strongly discouraged and I do not see an easy way to leverage a __getattr__ approach to implement this. In light of this, unless there's an alternative way to build the dtype, I don't think we should do this as it does not allow for the traditional indexing.

[ketch]

I recall finding the same problem when I tried it.

I agree that we shouldn't implement anything that is backwards-incompatible -- the current way of indexing must continue to work.

What about creating a dictionary of views into the q array? I tried it out:

http://nbviewer.ipython.org/urls/dl.dropboxusercontent.com/u/656693/Named%20fields%20as%20a%20dict%20of%20references.ipynb

Seems to work just fine. I think it would be nice to incorporate this.

[ketch]

Sorry, here's a better link: http://nbviewer.ipython.org/urls/dl.dropboxusercontent.com/u/656693/Named%20fields%20as%20a%20dict%20of%20references.ipynb

[mandli]

I thought about this approach but I think it might be difficult to ensure that copies of q do not get made and ensuring that the fields are truly pointing at the same thing. The obvious alternative to this is clearly to subclass the ndarray object which we should not (cannot) do. Otherwise I am not sure how to do this.

[ketch]

You're probably right. But just to humor me, could you give an example of something that one might normally do that would create a copy?

[mandli]

The code

import numpy

q = numpy.empty((3,4,4))

fields = {'density':q[0,:,:], 'momentum':q[1,:,:], 'energy':q[2,:,:]}

for n in xrange(3):
    q[n,:,:] = float(n)

fields['energy'] = numpy.exp(fields['energy'])

print q[2,:,:]
print fields['energy']

gives the output

[[ 2.  2.  2.  2.]
 [ 2.  2.  2.  2.]
 [ 2.  2.  2.  2.]
 [ 2.  2.  2.  2.]]
[[ 7.3890561  7.3890561  7.3890561  7.3890561]
 [ 7.3890561  7.3890561  7.3890561  7.3890561]
 [ 7.3890561  7.3890561  7.3890561  7.3890561]
 [ 7.3890561  7.3890561  7.3890561  7.3890561]]

Clearly this is somewhat contrived but it may lead people to unexpected things.

[ketch]

Alright, I'm convinced. But let me propose something different.

We could implement solution.fields as a read-only property. This would guarantee that it's always an up-to-date view of the appropriate part of q. We could:

• Automatically pick up the field names based on the choice of Riemann solver (for now, this can be done in a hacky way like how we get num_eqn; later it could be f90 module data)
• Write the field names as metadata to the output files (we should start doing this with HDF5 files anyway)
• Read the field names back in when loading data and (for instance) automatically set visclaw plot titles with them

I would like to implement this in a way that doesn't require field names to be given for a new Riemann solver, if possible. I.e., the code shouldn't break if there are no field names (also for backward compatibility).

[mandli]

I agree on all points. I had a simple list implementation that I was thinking might at the very least be written out and read in from file formats that supported such things (similar to units in the dimensions).

Making this even fancier (i.e. a read-only object that would behave properly when called) will take a bit of experimentation but I think this will be a lot easier than what we were trying to do before.