Jak skutecznie łączyć projektowanie klas i matematykę macierzową?

Question

Od pewnego czasu jestem psychicznie nękany przez zderzenie dwóch filozofii projektowania modeli fizycznych systemów i zastanawiam się, jakie rozwiązania wymyśliła społeczność.

Dla złożonych (er) symulacji, uwielbiam abstrakcję tworzenia klas dla obiektów i jak instancje obiektów klas można identyfikować z rzeczywistymi obiektami, które chcę studiować i jak pewne atrybuty obiektu reprezentują fizyczne cechy przedmioty prawdziwego życia. Weźmy balistyczne systemy cząsteczkowe jako prosty przykład:

class Particle(object): 
    def __init__(self, x=0, y=0, z=0): 
     self.x = x 
     self.y = y 
     self.z = z 
    def __repr__(self): 
     return "x={}\ny={}\nz={}".format(self.x, self.y, self.z) 
    def apply_lateral_wind(self, dx, dy): 
     self.x += dx 
     self.y += dy 

Gdybym zainicjować to z milion wartości mogę to zrobić:

start_values = np.random.random((int(1e6),3)) 

particles = [Particle(*i) for i in start_values] 

Teraz załóżmy, że muszę zrobić konkretnej rzeczy do wszystkich moich cząsteczek, takich jak dodanie bocznego wektora wiatru, tylko powoduje siekierę, przesunięcie y dla tej konkretnej operacji, ponieważ mam tylko garść (listę) wszystkich moich cząstek, musiałem pętnąć wszystkie moje cząstki, aby to zrobić i zajmuje to tyle czasu:

%timeit _ = [p.apply_lateral_wind(0.5, 1.2) for p in particles] 
1 loop, best of 3: 551 ms per loop 

Teraz przeciwstawnych oczywiste paradygmat dla tego, co jest oczywiście bardziej wydajny, ma pozostać na numpy poziomie i po prostu zrobić operację matematyczną bezpośrednio na tablicy, która jest ponad 10-krotnie szybszy:

%timeit start_values[...,:2] += np.array([0.5,1.2]) 
10 loops, best of 3: 20.3 ms per loop 

Moje pytanie brzmi: czy istnieją jakieś wzorce projektowe, które skutecznie łączą te dwa podejścia, aby nie stracić tak dużej wydajności? Uważam, że osobiście naprawdę łatwiej jest myśleć w kategoriach metod i atrybutów obiektów, jest to o wiele jaśniejsze w mojej głowie, a także dla mnie również podstawową przyczyną sukcesu programowania obiektowego (lub jego użycia w (fizycznym) modelowaniu) . Ale wady są oczywiste. Czy chciałbyś, jeśli jest jakaś elegancka wymiana zdań pomiędzy tymi podejściami?

Answer 1

Można zdefiniować klasę, która dba o wielu cząstek:

class Particles(object): 
    def __init__(self, coords): 
     self.coords = coords 

    def __repr__(self): 
     return "Particles(coords={})".format(self.coords) 

    def apply_lateral_wind(self, dx, dy): 
     self.coords[:, 0] += dx 
     self.coords[:, 1] += dy 

start_values = np.random.random((int(1e6), 3))   
particles = Particles(start_values) 

Te czasy na moich systemów pokazuje, że jest to rzeczywiście szybciej niż wersja numpy, przypuszczalnie dlatego, że nie budować dodatkową tablicę i robi nie trzeba się martwić o radiofonii:

%timeit particles.apply_lateral_wind(0.5, 1.2) 
100 loops, best of 3: 3.17 ms per loop 

natomiast przy użyciu numpy przykład daje

%timeit start_values[..., :2] += np.array([0.5, 1.2]) 
10 loops, best of 3: 21.1 ms per loop 

Answer 2

Cython extension types może być interesującą opcją, jeśli naprawdę chcesz użyć obiektu zamiast czystej tablicy NumPy. (. Mimo, że również pochodzić z pewnymi ograniczeniami, które są opisane w tej połączonej stronie dokumentacji)

edytowany kod przykładowy nieco, aby go Cythonize i podszedł z tym:

cdef class Particle(object): 
    cdef double x, y, z 

    def __init__(self, double x=0, double y=0, double z=0): 
     self.x = x 
     self.y = y 
     self.z = z 
    def __repr__(self): 
     return "x={}\ny={}\nz={}".format(self.x, self.y, self.z) 
    cpdef apply_lateral_wind(self, double dx, double dy): 
     self.x += dx 
     self.y += dy 

Z tego podstawowego Cython wersja, mam następujący harmonogram gdy zapętlenie ponad 1 milion cząsteczek:

>>> %timeit _ = [p.apply_lateral_wind(0.5, 1.2) for p in particles] 
10 loops, best of 3: 102 ms per loop 

Chodzi o czynnik 5 szybciej niż zwykły wersji Pythona, który odbył 532 ms na tej samej maszynie.Mimo to, jest wyraźnie wolniejszy o rząd wielkości niż przy użyciu czystej tablicy NumPy, ale myślę, że to cena czytelności.