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

4.2

не? А как тогда процесс происходит?

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

Ну а я о чем. Если брать из репов, то numpy берет тот же blas, что и программа на С.

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

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

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

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

Ну а я о чем.

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

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

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

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

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

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

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

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

Питон и программирование — штуки совершенно несовместимые!

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

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

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

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

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

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

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

\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 - матрицы

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

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

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

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

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

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

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

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

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

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

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

Мне нужен fixed_quad, только трехмерный)

не, ну блин, неужели не видно там сразу ссылки на tplquad?

Там другое, это не квадратура гаусса

Да, ты прав.

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

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

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

у тебя изначально область интегрирования какая?

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

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

Так а вы говорите, что и 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-репах.

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

Сделал тест на с++ и на питоне. Питон проигрывает пока в 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

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

ты рукожоп, с таким кодом тебе никакой 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

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

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

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)

я не понимаю чего ты хочешь. f возвращает у тебя матрицу, которую тупо за скобки можно вынести

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

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

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

\[K\] = \int\limits_{-1}^1 \int\limits_{-1}^1 \int\limits_{-1}^1 { \[B\]^T \[D\] \[B\] } dV,

\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],

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

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

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

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

То есть в предложенных автором языках (C++ и Python) они есть?

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

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

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

Именно поэтому на нём написано тонны мат. кода? Всякие eigen, armadillo и тд.

На фортране написано ещё больше мат. кода, и он хоть как-то может в многомерные массивы.

А C++ не может? Я не спорю, что фортран под это больше заточен, но мы говорим про современные языки.

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

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

Как будто что-то плохое.

Таки плохое, причём не только из-за тотального уродства и убогости подобного кода.

Эстеты в треде. Всё, что не haskell - гавно. Я правильно понимаю?

У меня скромные требования. Даже FORTRAN может чуть лучше, не говоря уже об APL'е.
Языки полувековой давности.

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

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

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

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

В нём нет зависимых типов, но зато в нём есть некоторые средства работы с многомерными массивами из коробки.

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