30

I'm trying to implement a simulation for a lattice model (lattice boltzmann) in Python. Each site of the lattice has a number of properties, and interact with neighboring sites according to certain rules. I figured that it might be clever to make a class with all the properties and make a grid of instances of that class. (As I'm inexperienced with Python, this might not be a good idea at all, so feel free to comment on my approach.)

Here is a toy example of what I'm doing

class site:
    def __init__(self,a,...):
        self.a = a
        .... other properties ...
    def set_a(self, new_a):
        self.a = new_a

Now I want to deal with a 2D/3D lattice (grid) of such sites so I tried to do the following (here is a 2D 3x3 grid as an example, but in simulation I would need the order of >1000x1000X1000)

lattice = np.empty( (3,3), dtype=object)
lattice[:,:] = site(3)

Now, the problem is that each lattice point refer to the same instance, for example

lattice[0,0].set_a(5)

will also set the value of lattice[0,2].a to 5. This behavior is unwanted. To avoid the problem i can loop over each grid point and assign the objects element by element, like

for i in range(3):
    for j in range(3):
        lattice[i,j] = site(a)

But is there a better way (not involving the loops) to assign objects to a multidimensional array?

Thanks

jonalm
  • 935
  • 2
  • 11
  • 21
  • 8
    If you're dealing with a >1000x1000X1000 array, _don't_ use an object array!! It will use egregious amounts of memory compared to using a "normal" numpy array. Object arrays aren't what you want here... – Joe Kington Feb 02 '11 at 18:40
  • 5
    by simulation I guess you mean fluid simulation? If so, then I'll recommend you to rethink your approach. Perhaps elements of your arrays should be numerical ones, so you could harness all the power of linear algebra ;-). Particle propagation and collision processes must be done globally! No local object lattice is able to handle that in any reasonable computation time. Just toughts, don't know really what you are aiming to ;-). Thanks – eat Feb 02 '11 at 21:57
  • @eat: I am doing fluid simulation. I wanted to code a generic grid of sites, where all the local properties were collected in a class (collision is local in my book, not propagation tho), but I guess you're right after all. At least bpowah taught me how to vectorize the __init__ function :) – jonalm Feb 02 '11 at 23:17
  • 1
    Incidentally, have you seen sailfish? http://sailfish.us.edu.pl/index.html It's a GPU accelerated fluid simulation package using a Lattice-Boltzman method implemented in numpy and pyopencl/pycuda. It's pretty slick from what I've seen (which are just the demo videos...). At any rate, thought you might find it relevant. – Joe Kington Feb 03 '11 at 01:06

3 Answers3

29

You can vectorize the class's __init__ function:

import numpy as np

class Site:
    def __init__(self, a):
        self.a = a
    def set_a(self, new_a):
        self.a = new_a

vSite = np.vectorize(Site)

init_arry = np.arange(9).reshape((3,3))

lattice = np.empty((3,3), dtype=object)
lattice[:,:] = vSite(init_arry)

This may look cleaner but has no performance advantage over your looping solution. The list comprehension answers create an intermediate python list which would cause a performance hit.

Paul
  • 42,322
  • 15
  • 106
  • 123
  • 1
    What if `Site` takes more than one parameter? – Neil G Apr 12 '13 at 00:47
  • 1
    @NeilG You'd have to encapsulate all the initialization information in `a` as a higher dimension array or recarray. The `vectorize`'d function only assumes that the first parameter is a vector (unless of course you used the `excluded` parameter to define a different one). If you insist on having `Site` multi-parameter, you could write a wrapper for `Site` that expands the recarray custom dtype to Site's `(*args, **vargs)` and vectorize the wrapper. Keep in mind that this is all just for fun. There is no performance advantage here over the OP's looping solution. – Paul Apr 12 '13 at 02:18
  • 2
    @Paul: Yes, I understand your point. I've never seen `vectorize` applied to a constructor before. It's a neat solution for this problem, but it is not a general solution since it makes adding parameters to `Site`'s constructor cumbersome. – Neil G Apr 12 '13 at 06:47
9

The missing piece for you is that Python treats everything as a reference. (There are some "immutable" objects, strings and numbers and tuples, that are treated more like values.) When you do

lattice[:,:] = site(3)

you are saying "Python: make a new object site, and tell every element of lattice to point to that object." To see that this is really the case, print the array to see that the memory addresses of the objects are all the same:

array([[<__main__.Site object at 0x1029d5610>,
        <__main__.Site object at 0x1029d5610>,
        <__main__.Site object at 0x1029d5610>],
       [<__main__.Site object at 0x1029d5610>,
        <__main__.Site object at 0x1029d5610>,
        <__main__.Site object at 0x1029d5610>],
       [<__main__.Site object at 0x1029d5610>,
        <__main__.Site object at 0x1029d5610>,
        <__main__.Site object at 0x1029d5610>]], dtype=object)

The loop way is one correct way to do it. With numpy arrays, that may be your best option; with Python lists, you could also use a list comprehension:

lattice = [ [Site(i + j) for i in range(3)] for j in range(3) ]

You can use a list comprehension with the numpy.array construction:

lattice = np.array( [ [Site(i + j) for i in range(3)] for j in range(3) ],
                    dtype=object)

Now when you print lattice, it's

array([[<__main__.Site object at 0x1029d53d0>,
        <__main__.Site object at 0x1029d50d0>,
        <__main__.Site object at 0x1029d5390>],
       [<__main__.Site object at 0x1029d5750>,
        <__main__.Site object at 0x1029d57d0>,
        <__main__.Site object at 0x1029d5990>],
       [<__main__.Site object at 0x1029d59d0>,
        <__main__.Site object at 0x1029d5a10>,
        <__main__.Site object at 0x1029d5a50>]], dtype=object)

so you can see that every object in there is unique.

You should also note that "setter" and "getter" methods (e.g., set_a) are un-Pythonic. It's better to set and get attributes directly, and then use the @property decorator if you REALLY need to prevent write access to an attribute.

Also note that it's standard for Python classes to be written using CamelCase, not lowercase.

Seth Johnson
  • 14,762
  • 6
  • 59
  • 85
  • Thank you so much for the input. But I don't know if I understand you correctly; For a class instance A = site(4,5), Is it better (more Pythonic) to update a variable with 'A.a = 6' than 'A.set_a(6)'? Any references on good guides to pythonic coding? – jonalm Feb 02 '11 at 21:19
  • Be careful with this is `Site` implements `__len__`. – Neil G Apr 12 '13 at 00:46
8

I don't know about better, but as an alternative to an explicit set of loops, you could write

lattice = np.empty( (3,3), dtype=object)
lattice.flat = [site(3) for _ in lattice.flat]

which should work whatever the shape of the lattice.

DSM
  • 342,061
  • 65
  • 592
  • 494