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

Это ты про «array notation» для сокращения кол-ва for в ряде случаев?

Технически фортран и есть современный. Все мировые ускорители на нём обсчитываются.

Спасибо за пример, стало понятно как писать на питоне.

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

Матрица B зависит от координат и каждый раз вычисляется

А что скажете на счёт Julia? Взлетит или нет?

на заданный тобой список вопросов:

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

с учётом указанным тобой достоинств и дефицитов отмеченных тобой лишь двух вариантов ответы в обратном порядке таковы:

• Есть!
• утвердительно
• смысл не только слово но и денотат.
• использовать само по себе не грех.
• целенаправленность деятельности есть признак осознанности деяния.
• а что есть не численное?
• всё есть расчёт числа
• для изготовления рецепта да. для эксплуатации уже окаменевшего способа погружаешь на плюсы либо на чистый си с возможным пипхолом на асм горячих частей для окначтельного выжимания последних процентов

По задаче: под интегралом матрица

размерность какая? По идее можно meshgrid уже по произведению матриц прогонять.

посмотри еще в сторону convolve. Это может еще ускорит

Неверно. Есть такой немаленький проект — QChem, написанный на плюсах. Ещё есть, mpqc, libint (на C) и, наверное, другие.

Чего на плюсах нет, это стандартной/стандартизованной библиотеки работы с матрицами/многомерными массивами/тензорами. Да, некоторые проекты используют, например, armadillo или Eigen, некоторые пишут свое (libtensor), но в целом именно это делает C++ менее удобным, чем Fortran 95/2003. Ну и ещё излишняя гибкость, не сдерживающая потуги отдельных товарищей ударится в дебри множественного наследования шаблонов и прочих контейнеров.

Если бы я сам знал

ну че, как успехи?

Норм, потихоньку пилю свою прогу. Сразу же выкинул свой костыль для работы с матрицами, в питоне все есть изкаробки))

Мозг ломается от питона, тут все по-другому.

Не совсем понятно как модуль из родительской папки импортировать в подпапку.

Скоро будут первые результаты, напишу про сравнение по скорости с сиплюсом

да, отпишись. У меня есть пару идей. Как обычно, одна — быстрая и память жрет, зато параллелится (с meshgrid). Но можно еще попробовать и без meshgrid на все твои матрицы делать (я сам еще не особо силен в этом деле, так что будет интересно)

Пока буду работать с обычными циклами (чтоб не парится), потом перепишу с помощью einsum итп.

Одна проблема осталась, вот структура проекта:

proj/
  ..b.py
    /fol1
    .. a.py

Могу загрузить из b.py модуль a.py. А наоборот не получается.

Для этого создал в /fol1 файл __init__.py: __all_ = ['a']

Что сделать для подключения снаружи?

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

Человеку нравится лисп и вангую, что он рад бы его использовать на работе, но обстоятельства не позволяют.

В общем, к совету прислушиваться стоит только если цель «изучить что-то новое/прикольное» важнее, чем «решить задачу побыстрее».

Я запускал сегодня через BLAS(от netlib.org) используя их бинарную сборку для Win32 - получилось медленнее чем реализация через циклы в лоб с оптимизациями при сборке (с предварительным транспонированием). Но, возможно, потому, что второй массив надо было всё же указать транспонированным. Да, ответ он мне выдаёт в виде транпонированной матрицы, видимо, сказывается, что библиотека фортрановская.

А вот с использованием бинарной сборки от OpenBLAS и флагом сборки -pthread выдаёт результат быстрее чем код с циклами (и транспонированиме) с использованием openmp (-fopenmp). Всё эьто для mingw64-4.9.3. Емнип, получилось в 2-2.5 раза быстрее кода с openmp.

P.S.

реализация BLAS от netlib.org захотела mingw64 именно с «exception handling method sjlj».

модуль из родительской папки

Должен быть добавлен глобальный путь sys.path до родительской папки. Читай https://leemendelowitz.github.io/blog/how-does-python-find-packages.html например. И не слушай эдди, по мне питон как плюсы без скобочек.

Спасибо. Это эдди из-под анонимуса пишет?

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

Не подскажу, у меня классы относительно небольшие. А IDE не пользуешься?

Ты можешь наследовать класс:

class foo():
    def aaa():
        dosomething

class bar(foo):
    def bbb():
        doanother

>>> nn = bar()
>>> nn.aaa()
>>> nn.bbb()

Я в sublime сижу пока. Pycharm слишком уж тормозной.

В общем, к совету прислушиваться стоит если цель «изучить что-то новое/другое» более интересна, чем «побыстрее завалить задачу кривыми костылями».

fixed

fixed

Осталось доказать, что в лиспе костыли менее кривые, чем в питоне.

Ну и ты правда думаешь, что автор на совершенно незнакомом языке родит решение более «прямое», чем на худо-бедно знакомом питоне?

Наркоманский синтаксис лиспа может даже посоперничать с пхытоновским!

Хотя, подозреваю, что лисп таки юзабельней, чем пхытоноговнище.

Тогда попробуй последний Eclipse с Pydev. Только отключи проверку орфографии, поставь тему darkest и color themes.

А можешь привести наивный код на плюсах?

...........

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

З.Ы.: в общем задело это сообщение, как результат:

#include <chrono>
#include <iostream>
#include <Eigen/Dense>

using std::chrono::duration;
using std::chrono::steady_clock;
using std::cout;
using std::endl;
using Eigen::MatrixXd;

int main()
{
  const size_t N = 4096;
  cout << "N: " << N << endl;

  MatrixXd a = 0.5*MatrixXd::Random(N,N)+MatrixXd::Ones(N,N);
  MatrixXd b = 0.5*MatrixXd::Random(N,N)+MatrixXd::Ones(N,N);

  auto t_start = steady_clock::now();
  MatrixXd c = a*b;
  auto t_end = steady_clock::now();

  cout << "Time: " << duration<double>(t_end - t_start).count()
       << " seconds" << endl;

  return 0;
}

$ export OMP_NUM_THREADS=4
$ g++ -O2 -DEIGEN_USE_BLAS -fopenmp -std=c++11 -I ~/opt/Eugen-3.3/include/eigen3 -lopenblaso test.cpp -o test
$ for i in {1..5} ; do time ./test ; done
N: 4096
Time: 1.65822 seconds

real    0m2.022s
user    0m6.753s
sys     0m0.158s
N: 4096
Time: 1.7515 seconds

real    0m2.067s
user    0m6.942s
sys     0m0.295s
N: 4096
Time: 1.53577 seconds

real    0m1.852s
user    0m6.322s
sys     0m0.068s
N: 4096
Time: 1.692 seconds

real    0m2.009s
user    0m6.802s
sys     0m0.206s
N: 4096
Time: 1.69986 seconds

real    0m2.016s
user    0m6.866s
sys     0m0.174s

import time
import numpy as np

N = 4096
print("N: %d" % N)

A = np.random.uniform(1.0, 2.0, (N,N))
B = np.random.uniform(1.0, 2.0, (N,N))

t_start = time.time()
C = np.matmul(A, B)
t_end = time.time()

print(" Time %s seconds :" % (t_end - t_start))

$ for i in {1..5} ; do time python test.py ; done
N: 4096
 Time 1.80444788933 seconds :

real    0m2.356s
user    0m7.511s
sys     0m0.468s
N: 4096
 Time 1.87723398209 seconds :

real    0m2.385s
user    0m7.672s
sys     0m0.553s
N: 4096
 Time 1.63046002388 seconds :

real    0m2.139s
user    0m6.906s
sys     0m0.337s
N: 4096
 Time 2.06297397614 seconds :

real    0m2.572s
user    0m7.214s
sys     0m0.861s
N: 4096
 Time 1.66444993019 seconds :

real    0m2.172s
user    0m6.972s
sys     0m0.405s

Я уже и не помню сколько ядер нагружала функция из numpy и не уверен (нужно проверить), что больше одного. Реализации BLAS сильно разные по скорости бывают. Референсная от netlib.org не самая быстрая и не умеет параллельность, емнип. А вот уже собранная для windows от OpenBLAS, похоже, всегда запускается в параллельном режиме, хотя в их виках в примере писали, что для параллельности при линковке нужен флаг "-lpthread", но его отсутствие никак не влияет на результат - всё так же грузит все ядра, включая виртуальные.

Реализации BLAS сильно разные по скорости бывают.

Так и есть.

Вот так не в какие ворота не лезет:

$ g++ -O2 -DEIGEN_USE_BLAS -fopenmp -std=c++11 -I ~/opt/Eugen-3.3/include/eigen3 -lblas test.cpp -o test-libblas
$ time ./test-libblas 
N: 4096
Time: 63.8347 seconds

real    1m4.207s
user    1m4.121s
sys     0m0.068s

А поменялась только слинкованная библиотека.

Даже с CentOS идёт несколько реализаций:

$ yum search blas
...
blas-devel.x86_64 : BLAS development libraries
blas-static.x86_64 : BLAS static libraries
blas64-devel.x86_64 : BLAS development libraries
blas64-static.x86_64 : BLAS static libraries (64bit INTEGER)
openblas.x86_64 : An optimized BLAS library based on GotoBLAS2
openblas-devel.x86_64 : Development headers and libraries for openblas-openmp.x86_64 : An optimized BLAS library based on GotoBLAS2, OpenMP version
openblas-openmp64.x86_64 : An optimized BLAS library based on GotoBLAS2, OpenMP version
openblas-openmp64_.x86_64 : An optimized BLAS library based on GotoBLAS2, OpenMP version
openblas-serial64.x86_64 : An optimized BLAS library based on GotoBLAS2, serial version
openblas-serial64_.x86_64 : An optimized BLAS library based on GotoBLAS2, serial version
openblas-static.x86_64 : Static version of OpenBLAS
openblas-threads.x86_64 : An optimized BLAS library based on GotoBLAS2, pthreads version
openblas-threads64.x86_64 : An optimized BLAS library based on GotoBLAS2, pthreads version
openblas-threads64_.x86_64 : An optimized BLAS library based on GotoBLAS2, pthreads version
...

А ещё MKL, сборка numpy с которым под Win-у легко находится.

Я уже и не помню сколько ядер нагружала функция из numpy и не уверен (нужно проверить), что больше одного.

Кроме меня тут есть и у других ситуации, когда time выдал user>real, так что это нормально, что numpy использует openmp по умолчанию.

Системная реализация openblas с потоками показала более не стабильные результаты, плюс общее время работы программы увеличилось

$ g++ -O2 -DEIGEN_USE_BLAS -pthread -std=c++11 -I ~/opt/Eugen-3.3/include/eigen3 -lopenblasp test.cpp -o test-openblas-pthreads
$ for i in {1..5}; do time ./test-openblas-pthreads ; done
N: 4096
Time: 1.78025 seconds

real    0m2.566s
user    0m7.713s
sys     0m0.383s
N: 4096
Time: 1.49981 seconds

real    0m2.437s
user    0m6.863s
sys     0m0.276s
N: 4096
Time: 1.48917 seconds

real    0m2.238s
user    0m6.626s
sys     0m0.276s
N: 4096
Time: 1.85577 seconds

real    0m2.603s
user    0m7.943s
sys     0m0.421s
N: 4096
Time: 1.77573 seconds

real    0m2.522s
user    0m7.591s
sys     0m0.456s

З.Ы.: надо проверять не на рабочей станции, а то результаты слишком различаются от запуска к запуску.

Хотя, подозреваю, что лисп таки юзабельней, чем пхытоноговнище.

А чем тебе именно питон так не угодил? В смысле, я-то думал, что ты «упоротыйёртый сишник», которому всё остальное не по душе.

exec(open('..\\b.py').read())

?

Оставьте напрастные попытки, Эдуарда вразумить.

Эдуард - ответ на все Вопросы (без лишних !?Ъ)

А что скажете на счёт Julia? Взлетит или нет?

Терпит крушение над водами Атлантического океана в течение 325 серий.

Грубо Денег дали им НА РАЗВИТИЕ не мало наблюдаем

Честно, не понимаю, почему julia не сольётся с Python на пользу человетчества

Конечная цель у них ОДНА?

Оставьте напрастные попытки, Эдуарда вразумить.

Вразумлять не пытаюсь, интересно логику понять.

Тем, что пхытон — абсолютно неюзабельное говнище! А на лиспе хоть какие-то программки вполне полезные бывают...

Осталось доказать, что в лиспе костыли менее кривые, чем в питоне.

Ну вот как там дела у питона с нативными тредами?

Ну и ты правда думаешь, что автор на совершенно незнакомом языке родит решение более «прямое», чем на худо-бедно знакомом питоне?

Никто автору не предлагает писать на языке, не изучив его хотя бы на уровне быстрого использования Hyperspec'а.

Почему чтобы изменить в питоне код метода надо перезапускать весь проект с интерпретатором? Почему такому якобы динамичному языку, позиционируемому как экспериментальная программная лаборатория, по динамизму далеко до всяких lua как раком до луны?

У лиспа тоже есть проблемы. Одна из этих проблем - после его изучение почти всё АЙТИ начинает казаться просто огромным стадом говноедов от которого потом хочется держаться подальше.

наличие большого количества циклов в реализаци мат. рассчетов обычно говорит о плохом уровне владения инструментом.

Вообще то это говорит о большом числе сумм/произведений в вычисляемом выражение, или большой размерности пространства;-) У меня напр 8й час крутится софтинка на ноуте - там пять вложенных циклов и никуда не денешься...

В итоге, если все выразить через матформулы и соответственно скормить в scipy, то все чики-пуки.

Мне бы очень хотелось увидеть решение банального скалярного волнового у-я в scipy для области 1024х1024х1024. Или системы у-й ландау-лифшица для атомистического моделирования магнетика, атомов так хотя бы на 10^4, 10^5 шагов. Сколько недель придется ждать окончания расчета?

Про многолучевую сейсмическую миграцию в асимптотическом приближении на scipy я уж промолчу...

Мне бы очень хотелось увидеть решение банального скалярного волнового у-я в scipy для области 1024х1024х1024.

Хватило бы памяти, а так — не вопрос.

Или системы у-й ландау-лифшица для атомистического моделирования магнетика, атомов так хотя бы на 10^4, 10^5 шагов. Сколько недель придется ждать окончания расчета?

Именно это мы как-то и делали в рамках курса вычислительной физики. Я лично делал задачу по распределению растительности на поверхности (это нелинейное уравнение второго порядка). По сути параболическое уравнение с нелинейным элементом.

А волновое уравнение — это вообще ерунда. Тупо convolve и шагаем по времени. Да и то, у меня есть подозрения, что по времени тоже можно без явных циклов сделать.

Мдя... Про волновое уравнение - пропускная способность памяти, memory bound/compute bound, roofline model - ни о чем не говорят слова?

Про магнетик - то же самое.

Как только возникает необходимость сделать что то стандартное с явной численной схемой для большого объема данных, пусть и помещающегося в память, но если задача мемори боунд (а она чем проще схема тем скорее будет мемори боунд - все эти numpy/scipy начинают радикально сливать обычному коду.

Как только возникает необходимость сделать что то не ложащееся по шерсти в рамках линейной алгебры - то же самое.

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

Мдя... Про волновое уравнение - пропускная способность памяти, memory bound/compute bound, roofline model - ни о чем не говорят слова?

ни о чем не говорят. Волновое уравнение — говоирт. bound — мало о чем говорит. Ты уравнение покажи лучше.

Про магнетик - то же самое.

про магнетик мой одногруппник делал.

Как только возникает необходимость сделать что то стандартное с явной численной схемой для большого объема данных, пусть и помещающегося в память, но если задача мемори боунд (а она чем проще схема тем скорее будет мемори боунд - все эти numpy/scipy начинают радикально сливать обычному коду.

да неужели?

Как только возникает необходимость сделать что то не ложащееся по шерсти в рамках линейной алгебры - то же самое.

если не ложится по шерсти линейной алгебры, значит надо ее туда положить.

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

Аааа, я понял. Ты из воинствующих «практиков»?

Почему я заинтересовался питоном - есть у нас в городе одна фирма, которая делает расчеты на элементах 1000х1000х1000 по времени (сколько там шагов не помню). Как я понял у них ядро на си, а поверх питон.

Это очень годная связка.

ни о чем не говорят. Волновое уравнение — говоирт. bound — мало о чем говорит. Ты уравнение покажи лучше.

Ну волновое уравнение я надеюсь Вы видели. Лучше Вы погуглите по ключевым словам которые я привел, м.б. поймете в чем именно Вы заблуждаетесь...

да неужели?

Гуглите, коллега, гуглите...

Аааа, я понял. Ты из воинствующих «практиков»?

и это говорит человек, выдавший

если не ложится по шерсти линейной алгебры, значит надо ее туда положить.

Скорее тут Вы - воинствующий сторонник линейной алгебры, ссылающийся на свои студенческие поделки и поделки своих однокурсников. Давайте я тогда тоже сошлюсь на задачи решенные моими студентами/аспирантами и на свои 20 лет опыта занятий числодробилками, если это у Вас уже катит за аргументы?;-)