LINUX.ORG.RU

Численные расчеты с python

 , , ,


1

4

Сегодня будет продолжение темы: https://www.linux.org.ru/forum/development/13234056

С того времени переписал базовую часть программы расчета на python.
Что-то пока сделано чтобы только работало.
Часть реализовано один в один как на плюсах.

Делюсь некоторыми результатами, надеюсь вам будет интересно.
sudo cast: dikiy, grem, unanimous, Shadow, AlexVR.

Результаты представил в виде графика

https://ibb.co/e7N8La

Расчеты сделаны для разного количества конечных элементов.
Есть два варианта на питоне: 1 - совсем сырой, 2-ой после беглого прохода профайлером.

Выводы: плюсы все-таки быстрее, но на питоне много чего можно оптимизировать, код пока сырой.
Есть еще возможность воспользоваться компиляцией и переписыванием медленной части на си.

★★★★★

Что-то пока сделано чтобы только работало.
Часть реализовано один в один как на плюсах.

Вам, вроде бы, советовали numpy. Вы его используете? Если просто переписать C++ на Python - то он естественно будет медленнее. Поэтому так не делают.

RazrFalcon ★★★★★
()
Ответ на: комментарий от sad_but_true1

Конечно numpy (в том числе lapack) много где используется. Еще не все с помощью него сделал, но основная часть на нем.

Zodd ★★★★★
() автор топика

Вот для примера не оптимизированный код, еще написанный несколько лет назад.
Очень много раз вызывается и занимает существенное время в расчете.

def DNx(self, Xi, Eta, Zeta, i):
        '''DNx - производная от функция формы fNi по x в лок. коорд.'''
        return self.JcbInv[0, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[0, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[0, 2] * self.DNZeta(Xi, Eta, i)

    def DNy(self, Xi, Eta, Zeta, i):
        '''DNy - производная от функция формы fNi по y в лок. коорд.'''
        return self.JcbInv[1, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[1, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[1, 2] * self.DNZeta(Xi, Eta, i)

    def DNz(self, Xi, Eta, Zeta, i):
        '''DNz - производная от функция формы fNi по z в лок. коорд.'''
        return self.JcbInv[2, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[2, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[2, 2] * self.DNZeta(Xi, Eta, i)

    def BMat(self, Xi, Eta, Zeta):
        '''Подпрограмма для вычисления матрицы связи (B) между перемещениями и
        деформациями (eps=B*u).'''
        nFree = self.nFree
        self.B.fill(0)
        for i in xrange(self.eNode):
            dnx = self.DNx(Xi, Eta, Zeta, i)
            dny = self.DNy(Xi, Eta, Zeta, i)
            dnz = self.DNz(Xi, Eta, Zeta, i)
            self.B[0, nFree * i] = dnx
            self.B[1, nFree * i + 1] = dny
            self.B[2, nFree * i + 2] = dnz
            self.B[3, nFree * i] = dny
            self.B[3, nFree * i + 1] = dnx
            self.B[4, nFree * i + 1] = dnz
            self.B[4, nFree * i + 2] = dny
            self.B[5, nFree * i] = dnz
            self.B[5, nFree * i + 2] = dnx

Есть варианты как ускорить, или на си переносить?

Zodd ★★★★★
() автор топика
Последнее исправление: Zodd (всего исправлений: 2)

Выводы: плюсы все-таки быстрее,

На них быстрее писать или быстрее выполняется ? Если второе, то что, были в этом сомнения ?

Dred ★★★★★
()
Ответ на: комментарий от Dred

Ну думал будет хуже, но меня больше расстраивает вторая часть графика. При решение больших задач, в основном упор должен идти на решение системы алгебраических уравнений (и там и там lapack), т.о. графики для питона и c++ должны были сближаться.

Zodd ★★★★★
() автор топика
Ответ на: комментарий от Zodd

Скорее всего это немного можно решить, улучшив код на питоне.

peregrine ★★★★★
()

Надеюсь, ты в курсе про Cython, и не станешь бездумно переписывать кучу кода на чистый С и подключать с помощью, прости Господи, ctypes (если совсем не терпится, то хоть воспользуйся cffi). И, конечно, тебе верно сказали, что стандартные сишные подходы (типа итерирования по массиву) даже в связке с numpy будут адово тормозить.

lu4nik ★★★
()
Ответ на: комментарий от Zodd

можно. Скорей всего x,y,z можно вычислять по отдельности. Тоже относится и к функции BMat.

А вообще неплохо бы полный код фстудию.

dikiy ★★☆☆☆
()
Ответ на: комментарий от ZERG

В этом и суть, что никаких. По возможности итерирование по массиву должно производиться в numpy (фактически это делается в сишном коде библиотеки):

# Скользящее среднее с окном размером 2.

a = numpy.random.random(10000)
res = numpy.empty(len(a) - 1)

# Медленно:
for i in range(len(a) - 1):
    res[i] = 0.5*(a[i] + a[i+1])
# Быстрee:
res = 0.5*(a[1:] + a[:-1])

На моей машине первый вариант выполняется за 2.96 мс, второй — за 156 мкс.

Конечно, в numpy есть средства для чуть более быстрого итерирования, если очень надо (numpy.nditer), но ими не стоит злоупотреблять. Если быстрый for позарез нужен, то Cython ускорит его в 12-15 раз, если для i, a и res отметить статические типы. Предполагаю, что ТС скопировал подобные циклы как есть из С++ и получил вышеприведенные тормоза.

lu4nik ★★★
()
Ответ на: комментарий от ZERG

https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/calculat...

Если я правильно понимаю, то чисто питоновские циклы так тормозят из-за динамичности переменной-счетчика. Видимо, разнообразные валидации типа не проходят даром. В numpу, насколько знаю, циклы реализованы обычным для С образом, но с оптимизациями вроде избегания cache miss и кэширования результатов.

lu4nik ★★★
()
Ответ на: комментарий от Zodd

Откуда эта прорва self?

Сдаётся мне, ты мыслишь в терминах циклов, а не матрично-векторных операций.

Ты когда-нибудь в Matlab/Octave/R программировал?

unanimous ★★★★★
()
Ответ на: комментарий от lu4nik

Оптимизации питоновские или на С? А то я планирую к минорным оптимизациям перейти сам, вот и ищу, с чего начать.

ZERG ★★★★★
()

Отвечая скорее на пред. тему: для числодробилок мы используем связку С++ и питон - вычислительное ядро на плюсах (то что должно считаться быстро), интерфейс на питоне (верхний управляющий слой, чтение конфигов и т.д.). Для биндинга используем SWIG (см напр. http://a-iv.ru/pyart/cpp2py.pdf).

В итоге время на разработку сокращается в разы без потери производительности полученного кода. Единственная проблема - развертывание этого на кластерах/суперкомпьютерах, не всегда есть свежий SWIG, и не всегда такая связка взлетает из под MPI. Но как правило когда дошло до кластера, это значит что код настолько отлажен и стабилизирован, что можно питон и открутить малой кровью (иногда).

AIv ★★★★★
()
Ответ на: комментарий от ZERG

Как я понял, единственная питоновская оптимизация там — это то, что по возможности дольше управление не передается в питон.

Об оптимизациях: возможно имеет смысл начать с Cython — раскидать статические типы и воспользоваться арифметикой из libc. По крайней мере, это не особо трудозатратно. Еще numexpr может помочь в длинных выражениях с большими ndarray — он использует JIT и может распараллеливать операции (тоже весьма легко попробовать).

lu4nik ★★★
()
Ответ на: комментарий от unanimous

Матлабами не пользовался. Это все еще старый код, тупо переведенный на питон. Понятно, что надо менять парадигму и переходить на матрично-векторный подход.

Может подскажешь что почитать на эту тему, либо примеры кода?

Zodd ★★★★★
() автор топика
Ответ на: комментарий от dikiy

Прошу сильно не пинать, т.к. этот кусок кода был написан давно, а также питона пока еще плохо знаю.

class Solid8():
    def __init__(self):
        self.nDef = 6
        self.nFree = 3
        self.eNode = 8
        self.eInt = 3
        self.Xi0 = [-1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0]
        self.Eta0 = [-1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0]
        self.Zeta0 = [-1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0]
        self.Jcb = np.empty((self.nFree, self.nFree))
        self.JcbInv = np.empty((self.nFree, self.nFree))
        self.B = np.zeros((self.nDef, self.nFree * self.eNode))

    def fNi(self, Xi, Eta, Zeta, i):
        return (1.0 + Xi * self.Xi0[i]) * (1.0 + Eta * self.Eta0[i])\
            * (1.0 + Zeta * self.Zeta0[i]) / 8.0

    def DNXi(self, Eta, Zeta, i):
        return self.Xi0[i] * (1.0 + Eta * self.Eta0[i])\
            * (1.0 + Zeta * self.Zeta0[i]) / 8.0

    def DNEta(self, Xi, Zeta, i):
        return self.Eta0[i] * (1.0 + Xi * self.Xi0[i])\
            * (1.0 + Zeta * self.Zeta0[i]) / 8.0

    def DNZeta(self, Xi, Eta, i):
        return self.Zeta0[i] * (1.0 + Xi * self.Xi0[i])\
            * (1.0 + Eta * self.Eta0[i]) / 8.0

    def Jacobian(self, Xi, Eta, Zeta, x0, y0, z0):
        self.Jcb = np.zeros((self.nFree, self.nFree))

        # Вычисление Jcb
        for i in xrange(self.eNode):
            self.Jcb[0, 0] += x0[i] * self.DNXi(Eta, Zeta, i)
            self.Jcb[0, 1] += y0[i] * self.DNXi(Eta, Zeta, i)
            self.Jcb[0, 2] += z0[i] * self.DNXi(Eta, Zeta, i)
            self.Jcb[1, 0] += x0[i] * self.DNEta(Xi, Zeta, i)
            self.Jcb[1, 1] += y0[i] * self.DNEta(Xi, Zeta, i)
            self.Jcb[1, 2] += z0[i] * self.DNEta(Xi, Zeta, i)
            self.Jcb[2, 0] += x0[i] * self.DNZeta(Xi, Eta, i)
            self.Jcb[2, 1] += y0[i] * self.DNZeta(Xi, Eta, i)
            self.Jcb[2, 2] += z0[i] * self.DNZeta(Xi, Eta, i)

        # Вычисление определителя
        det = np.linalg.det(self.Jcb)

        # Вычисление обратной матрицы JcbInv
        self.JcbInv = np.linalg.inv(self.Jcb)
        return det

    def DNx(self, Xi, Eta, Zeta, i):
        return self.JcbInv[0, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[0, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[0, 2] * self.DNZeta(Xi, Eta, i)

    def DNy(self, Xi, Eta, Zeta, i):
        return self.JcbInv[1, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[1, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[1, 2] * self.DNZeta(Xi, Eta, i)

    def DNz(self, Xi, Eta, Zeta, i):
        return self.JcbInv[2, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[2, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[2, 2] * self.DNZeta(Xi, Eta, i)

    def BMat(self, Xi, Eta, Zeta):
        nFree = self.nFree
        self.B.fill(0)
        for i in xrange(self.eNode):
            dnx = self.DNx(Xi, Eta, Zeta, i)
            dny = self.DNy(Xi, Eta, Zeta, i)
            dnz = self.DNz(Xi, Eta, Zeta, i)
            self.B[0, nFree * i] = dnx
            self.B[1, nFree * i + 1] = dny
            self.B[2, nFree * i + 2] = dnz
            self.B[3, nFree * i] = dny
            self.B[3, nFree * i + 1] = dnx
            self.B[4, nFree * i + 1] = dnz
            self.B[4, nFree * i + 2] = dny
            self.B[5, nFree * i] = dnz
            self.B[5, nFree * i + 2] = dnx

elem = Solid8()
x0 = [0, 1, 1, 0, 0, 1, 1, 0]
y0 = [0, 0, 1, 1, 0, 0, 1, 1]
z0 = [0, 0, 0, 0, 1, 1, 1, 1]
det = elem.Jacobian(0, 0, 0, x0, y0, z0)
elem.BMat(1, 2, 3)
print(det)
print(elem.B)
Zodd ★★★★★
() автор топика
Последнее исправление: Zodd (всего исправлений: 1)
Ответ на: комментарий от Zodd

весь код пока не смотрел, но якобиан так считать в пистоне не очень хорошо. Используй convolve (я уже упоминал это несколько раз). Рассчет градиента — это по сути свертка функции с определенной «маской».

Или numdifftools заюзай. Еще глянь в сторону gradient

dikiy ★★☆☆☆
()
Ответ на: комментарий от lu4nik

Так мне не нужен питон :) Я ищу способы свой код на С отполировать без использования openmp/mpi на данный момент.

ZERG ★★★★★
()
Последнее исправление: ZERG (всего исправлений: 1)
Ответ на: комментарий от ei-grad

Тоже была мысль предложить theano для якобиана.
Ещё есть Numdifftools и SymPy

Bell
()
Ответ на: комментарий от ism

GPU поможет, если оверхед на инициализацию и общение с ним компенсируется объемом и распараллеливанием вычислений. Это не всегда так. Кстати, у theano который вспомнили выше, большой оверхед на прекомпиляцию, и по-моему она не зависит от бэкенда - CPU/GPU.

Bell
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.