LINUX.ORG.RU

Есть ли смысл использовать для численных расчетов python?

 , ,


6

6

Есть ли смысл использовать для численных расчетов python (методы конечных элементов, математические расчеты, много циклов, большие данные)?

Или лучше использовать c++? Насколько медленнее код получается?

Плюсы питона:

  • простой
  • можно быстро создать алгоритм или переделать
  • много библиотек для численных расчетов, построения графиков
  • низкий порог вхождения, другим будет проще его изучить
  • можно потом куски на си переписать

Минусы питона:

  • медленнее плюсов
  • после c++ трудно переключится, кое-что по-другому (структуры, switch)
  • я его гораздо хуже знаю

Дал прогу на c++ одному, от так и не смог его осилить :(

Поделитесь историей успеха.

★★★★★

Ответ на: комментарий от dikiy

Код слинкуется с тем, с чем скажешь (или сказал тот, кто собирал код для тебя). numpy из репозитория ТСа - с blas'ом, предоставленным ему дистрибутивом (что оказалось дико тормознутым решением). numpy из репозитория pypi - с openblas'ом, который pip притащит вместе с numpy. numpy на Ломоносове - с atlas'ом, который предоставил Intel. numpy собранный лично тобой - с тем, с чем ему укажешь.

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

Не, я понимаю, ты бы матлаб расхваливал — там реально круто в последнее время стало: и тебе многопоточность, и CUDA!

Но говнопхытон...

Не понимаю. Пхытон — такой шлак, что с матлабом вообще не сравнить! Да даже у октавы он отсасывает по полной! Так на кой хрен пхытон юзать?

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

Что С, что пистон пользуются же одной и той же линамической либой.

Ну а я о чем.

Явно о чём-то глубоко своём, личном, потаённом.

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

numpy слинкован с blas, который идет в репах. При линковке программы на С, она линкуется с той же самой blas.

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

ТС вопрошал про питон, а не про октаву.

В октаве нет такой шняги, как matplotlib. Ну и вообще в программировании питон посильнее.

PS
пользователь octave

dikiy ★★☆☆☆
()

методы конечных элементов, математические расчеты, много циклов, большие данные

А можешь описать подробнее задачи?

Почему не рассматриваете Clojure или CL, не хватает каких-то конкретных библиотек? Скорость разработки / прототипирования в них выше, чем в Python - при наличии необходимых предметоспецифичных либ. А под капотом очень шустрые глобальнонадежные компиляторы - по завершению этапа прототипирования, ничего не надо будет переписывать на другие языки с левой семантикой, да и на этапе исследования алгоритмов скорость для вычислительных задач никогда лишней не бывает. Есть либы для численных методов и разного матана, machine learning, NLP, statistics, ... Clojure хорошо параллелится на кластеры метапарадигм^W. Много учебников для новичков, неудобные, но быстрые java-библиотеки подключаются и работают как нативные.

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

Мне понравилось. После C и подоных ему (Pascal) пистон предстает как калькулятор с человеческим лицом.

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

А можешь описать подробнее задачи?

Расчеты конструкций на прочность, МДТТ.

Почему не рассматриваете Clojure или CL

Потому что лисп я плохо знаю и он не сложнее питона. Задача у меня плохо параллелиться, поэтому этот вопрос не так важен.

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

Хочу оценить производительность, не подскажешь как лучше подсчитать тройной интеграл в питоне?
Не пойму как лучше это реализовать.

Интеграл вычисляется через квадратуру Гаусса (тройная сумма):

\sum_i\sum_j\sum_k {a_i a_j a_k transpose(B).D.B dV}

На плюсах схематично такой код:

for (int l=0; l<el; ++l) {
  dK = 0;
  for (int i=0; i<iInt; ++i)
    for (int j=0; j<iInt; ++j)
      for (int k=0; k<iInt; ++k) {
        matmul(D,B,Bt); 
        matmul(B,Bt,dK); // dK = B^T.D.B
        dK *= a_i*a_j*a_k*dV
        sum(dK,K); //K += dK
      }
  }

P.S. dK, K, B, D - матрицы

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

Расчеты конструкций на прочность, МДТТ.
Задача у меня плохо параллелиться

А что в МДТТ такого, что сложно параллелится? Интегралы (большинство), матрицы и прочая линейная алгебра - параллелятся на ура. А без параллелизации - какие там вообще могут быть вычисления на одном ядре :D

Потому что лисп я плохо знаю и он не сложнее питона

Он и сложнее, и мощнее питона. Но как я понял, вы выбираете ровно из двух языков. Всегда удивляло это в людях, которые профессионально занимаются каким-нибудь «матаном», «наукой» итд - вам ни инфраструктура там особо не нужна, ни фреймворки с формочками, ни ынтерпрайз и поддержка кода командами индусов из 200 человек, ни рынок заменяемых разрабочиков по доскам объявлений. Нужны - профессиональные либы со всякой банальщиной, а если что-то за пределами банальщины, ведь все равно придется все самим писать. И они выбирают, как ни жабу, так плюсы. Для исследовательских задач, ага.

Ох уже эти ученые... :)

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

Так в чем плюсы лиспа так я и не понял? Его мало кто знает.

Вот у питона куча библиотек, есть и cuda и даже распараллелить немного, также можно не только посчитать что-то, но и вывести через matplotlib результаты в виде графиков и сразу же вставить их в latex.

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

Скорей с жопомордой. Потому как это лицом называть никак язык не поворачивается! Пхытон — это же полное днище!

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

может тебе просто воспользоваться готовыми функциями из scipy для подсчета интегралов?

https://docs.scipy.org/doc/scipy-0.18.1/reference/tutorial/integrate.html

И я не до конца понимаю какой ты интеграл считаешь. Судя по приведенной тобой строчке там какая-то трансформация базиса. Напиши может просто формулу?

dikiy ★★☆☆☆
()
Последнее исправление: dikiy (всего исправлений: 1)
Ответ на: комментарий от Zodd

Да, ты прав.

В общем я бы на твоем месте попробовал порыть по ортогональным полиномам в многомерном пространстве - это раз. А по-твоему коду (насколько я вижу)

матрицу B^t D B можно вынести за скобки а остальное numpy.einsum

А если матрицы таки зависимы (ты просто не описал этого), то все равно einsum.

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

в общем твоя задача сделать так, чтобы функиця f, получая на вход ветора, возвращала тензор. Остальное сделает за тебя einsum (уже писал парой страниц выше).

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

Так в чем плюсы лиспа так я и не понял? Его мало кто знает.

Так а вы говорите, что и python не особо знаете, команда, как я подозреваю не знает его вообще. Какая тогда разница, что учить вместо С++?

Вот у питона куча библиотек

Да, потенциально либ там вагон. Но вопрос в другом - нужны ли вам все те либы, часто ли вы используете библиотеки «как есть», и как много приходится дописывать функционал существующих библиотек, изучая непосредственно их код.

В Java тоже либ достаточно много для расчетов чего угодно + скорость. И из clojure все это можно юзать нативно (!) - используешь java-методы и классы как родной код clojure-функций.

Вот у питона куча библиотек, есть и cuda и даже распараллелить немного

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

и cuda

Хз, я сам каким-то супер тяжелыми вычислениями на GPU не занимался. Но у тебя там была фраза про «мои задачи плохо параллеляться» (хотя я в это не верю), какое тогда нахрен вычисление на GPU? Если у тебя там из ресурсов или одна видео-карточка или одно хехе :) ядро проца.

Но чем не подойдет, допустим, OpenCL?

http://clojurecl.uncomplicate.org/articles/guides.html
http://neanderthal.uncomplicate.org/articles/tutorial_opencl.html

Gihub у этих проектов вполне себе живой и развивающийся. Судя по описанию - можно пользоваться обычной документацией для OpenCL.

Если нужна именно CUDA, наверняка ее можно использовать из Java. Но на этот счет утверждать не буду, хотя скорее всего так оно и есть.

также можно не только посчитать что-то, но и вывести через matplotlib

Матплотлиб это ж, при всем уважении, кусок говна. Куча кода кишками наружу, того, что надо нет, достаточно примитивная программка для рисования графиков в стиле «нарисовал и забыл», расширить эту хрень практические нереально без подсаживания на водку и средней тяжести опиаты. API убойное, куча кода депрекейтед.

Для clojure есть incanter - https://github.com/incanter/incanter#docs на сколько он жив неизвестно, последние активные коммиты 3 года назад, но существующий функционал был весьма нормальным, и главное - что там человеческий код, который можно поддерживать и расширять. Сделана поверх Parellel Colt https://github.com/rwl/ParallelColt и JFreeChart http://www.jfree.org/jfreechart/samples.html. clojure-api всегда можно расширить собственными (небольшими) усилиями. Еще есть такая вещь как Quil: http://quil.info/examples. Наверняка, если погуглить - найдется что-то еще. Даже если ничего не найдется, ты можешь просто из clojure рисовать на JFreeChart, уверен он ничем не хуже matplot-a.

результаты в виде графиков и сразу же вставить их в latex

Ну здесь я вообще никак проблем не вижу. Можно даже latex-описания к графикам херачить: https://data-sorcery.org/category/latex/.

Так в чем плюсы лиспа так я и не понял?

* интерактивная разработка в REPL (самый большой плюс)
* макросы (!)
* скорость, быстрая JVM под капотом со всеми плюшками - для вычислительных задач плюс несомненный, сравни циферки
* быстрые java-либы для вычислений нашару без необходимости ковыряться в java, C, C++ и прочем xxx
* халявный параллелизм и многопоточность, встроенные в язык
* функциональный язык - не поощряет неявных мутабельных состояний -> простота отладки, идиоматичный код
* удовольствие от процесса разработки

Для исследовательских задач я бы тоже выбирал из языков: python, clojure, CL, racket. Но на python остановился бы _только_ если бы в его инфраструктуре были супер критичные задачеспецифичные либы, которых больше нигде нет. При всех прочих равных, естественно бы, не выбирал.

Как минимум, проведите хотя бы небольшой рисерч по либам, которые вам реально нужны. А вдруг в clojure + java все это есть?

Clojure бы не взял в одном случае - если у вас уже известны алгоритмы и большая часть их императивщина по техническим причинам (серьезные требования по скорости на одном вычислительном узле или специфика самих алгоритмов) с постоянныеми инплейс-заменами и прочими мутабельными состояниями (речь не о циклах и не о матрицах - «немутабельные» матрицы прекрасно изменяются без копирования, а для циклов есть какие-то либы в стиле CL-макроса iterate, который раскрывается в обычный кложуровский loop-recur) и их сложно / неестественно выражать в функциональном стиле (рекурсия, неизменяемые данные), но в таком случае python вам тоже не подойдет - только java, C, C++ (и то не факт), а может и исключительно фортран использовать надо. И дальше руками все параллелить. При таком раскладе можно посмотреть еще на Common Lisp. Но батареек там будет поменьше, чем для clojure. Но опять же, смотря, что вам надо: http://www.cliki.net/Mathematics, http://cliki.net/plotting - здесь, возможно, уголок некрофилии, смотреть по обновлениям в dvcs-репах.

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

Если тебе самому нужно реализовывать такие циклы, то в самом питоне они медленные очень.

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

Сделал тест на с++ и на питоне. Питон проигрывает пока в 4 раза.

Ниже привел код (убрал все лишнее, оставил основную часть).

#!/usr/bin/env python
from __future__ import print_function

import numpy as np
import random
import time

el = 10000
node = 20
iInt = 3
a = [-1., 0., 1.]

B = [[random.uniform(0, 1.1) for x in range(node)] for y in range(node)]
D = [[random.uniform(0, 1.1) for x in range(node)] for y in range(node)]

dK = np.array([node, node])
K = np.array([node, node])
B = np.array(B)
Bt = np.array([node, node])
D = np.array(D)

K = 0

start_time = time.time()

for i in xrange(el):
    for m in xrange(iInt):
        for n in xrange(iInt):
            for l in xrange(iInt):
                koef = a[m]*a[n]*a[l]
                Bt = np.matmul(D, B)
                dK = np.matmul(koef*B, Bt)
                K += dK

print(" Time %s seconds:" % (time.time() - start_time))

$ python ktest.py
Time 3.4896569252 seconds:
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <ctime>

using namespace std;

const int el=10000;
const int node=20;
const int iInt = 3;
double dK[node*node], K[node*node], B[node*node], Bt[node*node], D[node*node], a[iInt];

void matmul(int n, double *a, double *b, double *c){
  double cc;
  for (int i=0; i<n; i++)
    for (int j=0; j<n; j++){
      cc = 0;
      for (int k=0; k<n; k++)
        cc += a[n*i+k]*b[n*j+k];
      c[n*i+j] = cc;
    }
  }

void sum(double koef, int n, double *a, double *b) {
  for (int i=0; i<n; i++)
    for (int j=0; j<n; j++)
      b[n*i+j]=b[n*i+j]+koef*a[n*i+j];
  }

int main() {
  double koef;
  a[0] = -1.;
  a[1] = 0.;
  a[2] = 1.;

  for (int i=0; i<node; i++)
    for (int j=0; j<node; j++)
      D[node*i+j] = (rand() % 100)/100.0+1.0;

    for (int i=0; i<node; i++)
      for (int j=0; j<node; j++)
        B[node*i+j] = (rand() % 100)/100.0+1.0;

  for (int i=0; i<node*node; ++i)
    K[i] = 0.;

  clock_t start;
  start = clock();

  for (int i=0; i<el; ++i)
    for (int m=0; m<iInt; ++m)
      for (int n=0; n<iInt; ++n)
        for (int l=0; l<iInt; ++l) {
          matmul(node,D,B,Bt);
          matmul(node,B,Bt,dK);
          koef = a[m]*a[n]*a[l];
          sum(koef,node,dK,K);
        }

  cout << (clock()-start)/(double) CLOCKS_PER_SEC << endl;
  return 0;
}

$ g++ -march=native -Ofast -o ktest ktest.cpp && ./ktest
0.72372
Zodd ★★★★★
() автор топика
Ответ на: комментарий от Zodd

Как можно код на питоне ускорить? Хотяб в 2 раза.

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

ты рукожоп, с таким кодом тебе никакой blas с numpy не поможет :D

этот код бессмысленный. Так как на выходе всегда ноль ибо сумма по a_m a_n a_l вычисляется (при твоих вводных) в ноль.

for i in xrange(el):
    for m in xrange(iInt):
        for n in xrange(iInt):
            for l in xrange(iInt):
                koef = a[m]*a[n]*a[l]
                Bt = np.matmul(D, B)
                dK = np.matmul(koef*B, Bt)
                K += dK


И еще раз повторяю: используй einsum.

dikiy ★★☆☆☆
()
Последнее исправление: dikiy (всего исправлений: 1)
Ответ на: комментарий от dikiy

Там другие коэффициенты, для простоты написал такие.

Подскажи, что у меня тут неправильно. Не до конца въехал в питон.

def f(B, D):
    return np.matmul(koef*B, np.matmul(D, B))

for i in xrange(el):
    K = np.einsum('i,j,k', a, a, a, f)
Zodd ★★★★★
() автор топика
Ответ на: комментарий от Zodd

f должен вернуть тензор размерности i_max на j_max на k_max

А у тебя матрица возвращается 20x20.

dikiy ★★☆☆☆
()
Последнее исправление: dikiy (всего исправлений: 1)
Ответ на: комментарий от dikiy

Мне нужно вычислить следующее:

\[K\] = \int\limits_{-1}^1 \int\limits_{-1}^1 \int\limits_{-1}^1 { \[B\]^T \[D\] \[B\] } dV,
где матрицы B, D, K - матрицы. Это приближенно можно вычислить так
\sum_i\sum_j\sum_k {a(i)*a(j)*a(k)* [B(x(i),x(j),x(k))]^T [D] [B(x(i),x(j),x(k))] * det[J],
a(i) - вектор весовых коэффициентов, x(i) - вектор координат, i,j,k=1..3

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

def f(x,y,z):
    return np.matmul(B(x,y,z),np.matmul(D,B(x,y,z))) # твоя функция. Насколько я понимаю, матрицы B функциональные


xx,yy,zz = np.meshgrid(x,y,z)  #здесь подразумевается, что x,y,z уже содержат нужные наборты опорных точек значения

f_meshed=f(xx,yy,zz)

np.einsum('i,j,k,ijk',a,a,a,f_meshed)
dikiy ★★☆☆☆
()
Последнее исправление: dikiy (всего исправлений: 1)
Ответ на: комментарий от Zodd

на выходе ты что хочешь получить? Матрицу или скаляр?

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

При всём моём уважении к расту, годным языком для численных расчётов он не является.
По крайней мере, пока в него не завезут зависимые типы.

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

В общем, покажи математически, чего тебе надо. Только точно. Начиная от исходной задачи.

Интеграл от матрицы? Имеет собсно смысл, только если у тебя матрица из функций. Иначе просто за интеграл можно вынести.

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

В питоне динамическая типизация, там это не критично. В С++ их можно сделать, хотя в случае с шаблонной магией «можно» обычно означает «лучше не надо». Но плюсы и не являются хорошим языком для численных расчётов.

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

Но плюсы и не являются хорошим языком для численных расчётов.

Как интересно. А что вообще является таким языком?

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

хотя в случае с шаблонной магией «можно» обычно означает «лучше не надо»


В С++ принято делать for for for, хотя и есть теоретическая возможность так не делать.

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

Даже FORTRAN может чуть лучш

Но ведь в нем нет зависимых типов. Чем же он лучше?

не говоря уже об APL'е.

А. Тогда вопросов нет.

tailgunner ★★★★★
()
Последнее исправление: tailgunner (всего исправлений: 1)
Ответ на: комментарий от RazrFalcon

О боже, фанаты раста уже срутся между собой. Как предсказуемо. Конечно же хачкелята и крестогеи там должны столкнуться и аннигилировать на радость всем.

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