numpy.amax() найдет максимальное значение в массиве, а numpy.amin() делает то же самое для значения min. Если я хочу найти как max, так и min, мне нужно вызвать обе функции, которые требуют прохождения над (очень большим) массивом дважды, что кажется медленным.
Есть ли функция в API-интерфейсе numpy, которая находит как max, так и min только с одним проходом через данные?
Есть ли функция в API-интерфейсе numpy, которая находит как max, так и min только с одним проходом через данные?
Нет. На момент написания этой статьи нет такой функции. (И да, если бы была такая функция, ее производительность была бы значительно лучше, чем вызов numpy.amin() и numpy.amax() последовательно на большом массиве.)
Я не думаю, что прохождение массива дважды является проблемой. Рассмотрим следующий псевдокод:
minval = array[0]
maxval = array[0]
for i in array:
if i < minval:
minval = i
if i > maxval:
maxval = i
Пока здесь всего 1 цикл, есть еще две проверки. (Вместо того, чтобы иметь 2 петли с 1 проверкой каждого). На самом деле единственное, что вы сохраняете, это накладные расходы на 1 цикл. Если массивы действительно велики, как вы говорите, эти накладные расходы малы по сравнению с фактической рабочей нагрузкой на петле. (Обратите внимание, что все это реализовано на C, поэтому петли более или менее свободны в любом случае).
РЕДАКТИРОВАТЬ Извините 4 из вас, кто поддержал меня и верю в меня. Вы определенно можете оптимизировать это.
Вот код fortran, который может быть скомпилирован в модуль python через f2py (возможно, гуру Cython может прийти и сравнить это с оптимизированной версией C...):
subroutine minmax1(a,n,amin,amax)
implicit none
!f2py intent(hidden) :: n
!f2py intent(out) :: amin,amax
!f2py intent(in) :: a
integer n
real a(n),amin,amax
integer i
amin = a(1)
amax = a(1)
do i=2, n
if(a(i) > amax)then
amax = a(i)
elseif(a(i) < amin) then
amin = a(i)
endif
enddo
end subroutine minmax1
subroutine minmax2(a,n,amin,amax)
implicit none
!f2py intent(hidden) :: n
!f2py intent(out) :: amin,amax
!f2py intent(in) :: a
integer n
real a(n),amin,amax
amin = minval(a)
amax = maxval(a)
end subroutine minmax2
Скомпилируйте его с помощью:
f2py -m untitled -c fortran_code.f90
И теперь мы находимся в месте, где мы можем его протестировать:
import timeit
size = 100000
repeat = 10000
print timeit.timeit(
'np.min(a); np.max(a)',
setup='import numpy as np; a = np.arange(%d, dtype=np.float32)' % size,
number=repeat), " # numpy min/max"
print timeit.timeit(
'untitled.minmax1(a)',
setup='import numpy as np; import untitled; a = np.arange(%d, dtype=np.float32)' % size,
number=repeat), '# minmax1'
print timeit.timeit(
'untitled.minmax2(a)',
setup='import numpy as np; import untitled; a = np.arange(%d, dtype=np.float32)' % size,
number=repeat), '# minmax2'
Результаты для меня немного ошеломляющие:
8.61869883537 # numpy min/max
1.60417699814 # minmax1
2.30169081688 # minmax2
Я должен сказать, я не совсем понимаю это. Сравнение только np.min против minmax1 и minmax2 по-прежнему проигрывает, поэтому это не просто проблема с памятью...
примечания. Увеличение размера в 10**a и уменьшение повторения с коэффициентом 10**a (постоянный размер задачи) меняют производительность, но не кажущимся последовательным способом, который показывает, что есть некоторые взаимодействие между производительностью памяти и служебными функциями вызова функции в python. Даже сравнение простой реализации min в fortran превосходит numpy примерно в 2...
Существует функция поиска (max-min), называемая numpy.ptp, если это полезно для вас:
>>> import numpy
>>> x = numpy.array([1,2,3,4,5,6])
>>> x.ptp()
5
но я не думаю, что есть способ найти как min, так и max с одним обходом.
Вы можете использовать Numba, динамический компилятор Python с поддержкой NumPy, использующий LLVM. Полученная реализация довольно проста и понятна:
import numpy
import numba
@numba.jit
def minmax(x):
maximum = x[0]
minimum = x[0]
for i in x[1:]:
if i > maximum:
maximum = i
elif i < minimum:
minimum = i
return (minimum, maximum)
numpy.random.seed(1)
x = numpy.random.rand(1000000)
print(minmax(x) == (x.min(), x.max()))
Это также должно быть быстрее, чем реализация Numpy min() & max(). И все это без необходимости писать одну строчку кода на C/Fortran.
Проведите собственные тесты производительности, так как они всегда зависят от вашей архитектуры, ваших данных, версий вашего пакета...
Это старый поток, но в любом случае, если кто-нибудь еще раз посмотрит на это...
При одновременном поиске мин и макс можно уменьшить количество сравнений. Если это поплавки, которые вы сравниваете (это, я думаю, это так), это может сэкономить вам некоторое время, хотя и не сложность вычислений.
Вместо (код Python):
_max = ar[0]
_min= ar[0]
for ii in xrange(len(ar)):
if _max > ar[ii]: _max = ar[ii]
if _min < ar[ii]: _min = ar[ii]
вы можете сначала сравнить два смежных значения в массиве, а затем сравнить меньшее значение с текущим минимумом, а большее - против текущего максимума:
## for an even-sized array
_max = ar[0]
_min = ar[0]
for ii in xrange(0, len(ar), 2)): ## iterate over every other value in the array
f1 = ar[ii]
f2 = ar[ii+1]
if (f1 < f2):
if f1 < _min: _min = f1
if f2 > _max: _max = f2
else:
if f2 < _min: _min = f2
if f1 > _max: _max = f1
Код здесь написан на Python, явно для скорости, которую вы использовали бы C или Fortran или Cython, но таким образом вы выполняете 3 сравнения на итерацию с len (ar)/2 итерациями, давая 3/2 * len (ar ) сравнения. В противоположность этому, делая сравнение "очевидным образом", вы делаете два сравнения на итерацию, что приводит к сопоставлениям 2 * len (ar). Экономит 25% от времени сравнения.
Возможно, кто-то однажды найдет это полезным.
В общем случае вы можете уменьшить количество сравнений для алгоритма minmax, обрабатывая одновременно два элемента и сравнивая меньший с временным минимумом, а другой - с временным максимумом. В среднем нужно всего 3/4 сравнений, чем наивный подход.
Это может быть реализовано в c или fortran (или любом другом низкоуровневом языке) и должно быть почти непревзойденным с точки зрения производительности. Я использую numba, чтобы проиллюстрировать принцип и получить очень быструю, независимую от dtype реализацию:
import numba as nb
import numpy as np
@nb.njit
def minmax(array):
# Ravel the array and return early if it empty
array = array.ravel()
length = array.size
if not length:
return
# We want to process two elements at once so we need
# an even sized array, but we preprocess the first and
# start with the second element, so we want it "odd"
odd = length % 2
if not odd:
length -= 1
# Initialize min and max with the first item
minimum = maximum = array[0]
i = 1
while i < length:
# Get the next two items and swap them if necessary
x = array[i]
y = array[i+1]
if x > y:
x, y = y, x
# Compare the min with the smaller one and the max
# with the bigger one
minimum = min(x, minimum)
maximum = max(y, maximum)
i += 2
# If we had an even sized array we need to compare the
# one remaining item too.
if not odd:
x = array[length]
minimum = min(x, minimum)
maximum = max(x, maximum)
return minimum, maximum
Это определенно быстрее, чем наивный подход, который Peque представил:
arr = np.random.random(3000000)
assert minmax(arr) == minmax_peque(arr) # warmup and making sure they are identical
%timeit minmax(arr) # 100 loops, best of 3: 2.1 ms per loop
%timeit minmax_peque(arr) # 100 loops, best of 3: 2.75 ms per loop
Как и ожидалось, новая реализация minmax займет примерно 3/4 времени наивности (2.1 / 2.75 = 0.7636363636363637)
Никто не упомянул numpy.percentile, поэтому я подумал, что буду. Если вы попросите процентили [0, 100], вы получите массив из двух элементов: минимальный (0-й процентиль) и максимальный (100-й процентиль).
Однако он не удовлетворяет цели OP: он не быстрее, чем min и max отдельно. Это, вероятно, из-за некоторого механизма, который позволил бы использовать неэкстремальные процентили (более сложная проблема, которая должна занять больше времени).
In [1]: import numpy
In [2]: a = numpy.random.normal(0, 1, 1000000)
In [3]: %%timeit
...: lo, hi = numpy.amin(a), numpy.amax(a)
...:
100 loops, best of 3: 4.08 ms per loop
In [4]: %%timeit
...: lo, hi = numpy.percentile(a, [0, 100])
...:
100 loops, best of 3: 17.2 ms per loop
In [5]: numpy.__version__
Out[5]: '1.14.4'
В будущей версии Numpy можно указать специальный случай, чтобы пропустить расчет нормального процентиля, если запрашивается только [0, 100]. Не добавляя ничего к интерфейсу, есть способ запросить у Numpy минимальное и максимальное значения за один вызов (вопреки тому, что было сказано в принятом ответе), но стандартная реализация библиотеки не использует этот случай, чтобы сделать его стоящим ,
На первый взгляд numpy.histogram, кажется, делает свое дело:
count, (amin, amax) = numpy.histogram(a, bins=1)
... но если вы посмотрите на источник для этой функции, он просто вызывает a.min() и a.max() независимо и поэтому не может избежать проблем с производительностью, рассмотренных в этом вопросе. :-(
Точно так же scipy.ndimage.measurements.extrema выглядит как возможность, но он также просто вызывает a.min() и a.max() независимо друг от друга.
В любом случае, это стоило мне усилий, поэтому я предложу самое трудное и наименее элегантное решение для тех, кто может быть заинтересован. Мое решение состоит в том, чтобы реализовать многопоточный алгоритм min-max за один проход в C++ и использовать его для создания модуля расширения Python. Эти усилия требуют некоторых дополнительных затрат для изучения того, как использовать API-интерфейсы Python и NumPy C/C++, и здесь я покажу код и дам несколько небольших пояснений и ссылок для тех, кто хочет идти по этому пути.
Здесь нет ничего слишком интересного. Массив разбит на куски размером length / workers. Мин./Макс. Рассчитываются для каждого фрагмента в future, который затем сканируется для поиска глобального мин./Макс.
// mt_np.cc
//
// multi-threaded min/max algorithm
#include <algorithm>
#include <future>
#include <vector>
namespace mt_np {
/*
* Get {min,max} in interval [begin,end)
*/
template <typename T> std::pair<T, T> min_max(T *begin, T *end) {
T min{*begin};
T max{*begin};
while (++begin < end) {
if (*begin < min) {
min = *begin;
continue;
} else if (*begin > max) {
max = *begin;
}
}
return {min, max};
}
/*
* get {min,max} in interval [begin,end) using #workers for concurrency
*/
template <typename T>
std::pair<T, T> min_max_mt(T *begin, T *end, int workers) {
const long int chunk_size = std::max((end - begin) / workers, 1l);
std::vector<std::future<std::pair<T, T>>> min_maxes;
// fire up the workers
while (begin < end) {
T *next = std::min(end, begin + chunk_size);
min_maxes.push_back(std::async(min_max<T>, begin, next));
begin = next;
}
// retrieve the results
auto min_max_it = min_maxes.begin();
auto v{min_max_it->get()};
T min{v.first};
T max{v.second};
while (++min_max_it != min_maxes.end()) {
v = min_max_it->get();
min = std::min(min, v.first);
max = std::max(max, v.second);
}
return {min, max};
}
}; // namespace mt_np
Вот где все начинает становиться ужасным... Одним из способов использования кода C++ в Python является реализация модуля расширения. Этот модуль может быть собран и установлен с использованием стандартного модуля distutils.core. Полное описание того, что это влечет за собой, описано в документации по Python: https://docs.python.org/3/extending/extending.html. ПРИМЕЧАНИЕ. Существуют и другие способы получения аналогичных результатов, например, https://docs.python.org/3/extending/index.html#extending-index:
Это руководство охватывает только основные инструменты для создания расширений, предоставляемых как часть этой версии CPython. Сторонние инструменты, такие как Cython, cffi, SWIG и Numba, предлагают как более простые, так и более сложные подходы к созданию расширений C и C++ для Python.
По сути, этот маршрут скорее академический, чем практический. С учетом сказанного, что я сделал дальше, придерживаясь довольно близко к уроку, создаю файл модуля. Это по сути шаблон для distutils, чтобы знать, что делать с вашим кодом, и создавать из него модуль Python. Прежде чем делать что-либо из этого, возможно, стоит создать виртуальную среду Python, чтобы не загрязнять системные пакеты (см. https://docs.python.org/3/library/venv.html#module-venv).
Вот файл модуля:
// mt_np_forpy.cc
//
// C++ module implementation for multi-threaded min/max for np
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <python3.6/numpy/arrayobject.h>
#include "mt_np.h"
#include <cstdint>
#include <iostream>
using namespace std;
/*
* check:
* shape
* stride
* data_type
* byteorder
* alignment
*/
static bool check_array(PyArrayObject *arr) {
if (PyArray_NDIM(arr) != 1) {
PyErr_SetString(PyExc_RuntimeError, "Wrong shape, require (1,n)");
return false;
}
if (PyArray_STRIDES(arr)[0] != 8) {
PyErr_SetString(PyExc_RuntimeError, "Expected stride of 8");
return false;
}
PyArray_Descr *descr = PyArray_DESCR(arr);
if (descr->type != NPY_LONGLTR && descr->type != NPY_DOUBLELTR) {
PyErr_SetString(PyExc_RuntimeError, "Wrong type, require l or d");
return false;
}
if (descr->byteorder != '=') {
PyErr_SetString(PyExc_RuntimeError, "Expected native byteorder");
return false;
}
if (descr->alignment != 8) {
cerr << "alignment: " << descr->alignment << endl;
PyErr_SetString(PyExc_RuntimeError, "Require proper alignement");
return false;
}
return true;
}
template <typename T>
static PyObject *mt_np_minmax_dispatch(PyArrayObject *arr) {
npy_intp size = PyArray_SHAPE(arr)[0];
T *begin = (T *)PyArray_DATA(arr);
auto minmax =
mt_np::min_max_mt(begin, begin + size, thread::hardware_concurrency());
return Py_BuildValue("(L,L)", minmax.first, minmax.second);
}
static PyObject *mt_np_minmax(PyObject *self, PyObject *args) {
PyArrayObject *arr;
if (!PyArg_ParseTuple(args, "O", &arr))
return NULL;
if (!check_array(arr))
return NULL;
switch (PyArray_DESCR(arr)->type) {
case NPY_LONGLTR: {
return mt_np_minmax_dispatch<int64_t>(arr);
} break;
case NPY_DOUBLELTR: {
return mt_np_minmax_dispatch<double>(arr);
} break;
default: {
PyErr_SetString(PyExc_RuntimeError, "Unknown error");
return NULL;
}
}
}
static PyObject *get_concurrency(PyObject *self, PyObject *args) {
return Py_BuildValue("I", thread::hardware_concurrency());
}
static PyMethodDef mt_np_Methods[] = {
{"mt_np_minmax", mt_np_minmax, METH_VARARGS, "multi-threaded np min/max"},
{"get_concurrency", get_concurrency, METH_VARARGS,
"retrieve thread::hardware_concurrency()"},
{NULL, NULL, 0, NULL} /* sentinel */
};
static struct PyModuleDef mt_np_module = {PyModuleDef_HEAD_INIT, "mt_np", NULL,
-1, mt_np_Methods};
PyMODINIT_FUNC PyInit_mt_np() { return PyModule_Create(&mt_np_module); }
В этом файле широко используется как Python, так и API NumPy, для получения дополнительной информации обратитесь к: https://docs.python.org/3/c-api/arg.html#c.PyArg_ParseTuple, а для NumPy: https://docs.scipy.org/doc/numpy/reference/c-api.array.html.
Следующее, что нужно сделать, это использовать distutils для установки модуля. Для этого требуется установочный файл:
# setup.py
from distutils.core import setup,Extension
module = Extension('mt_np', sources = ['mt_np_module.cc'])
setup (name = 'mt_np',
version = '1.0',
description = 'multi-threaded min/max for np arrays',
ext_modules = [module])
Чтобы наконец установить модуль, выполните python3 setup.py install из своей виртуальной среды.
Наконец, мы можем проверить, действительно ли реализация C++ превосходит наивное использование NumPy. Для этого вот простой тестовый скрипт:
# timing.py
# compare numpy min/max vs multi-threaded min/max
import numpy as np
import mt_np
import timeit
def normal_min_max(X):
return (np.min(X),np.max(X))
print(mt_np.get_concurrency())
for ssize in np.logspace(3,8,6):
size = int(ssize)
print('********************')
print('sample size:', size)
print('********************')
samples = np.random.normal(0,50,(2,size))
for sample in samples:
print('np:', timeit.timeit('normal_min_max(sample)',
globals=globals(),number=10))
print('mt:', timeit.timeit('mt_np.mt_np_minmax(sample)',
globals=globals(),number=10))
Вот результаты, которые я получил от всего этого:
8
********************
sample size: 1000
********************
np: 0.00012079699808964506
mt: 0.002468645994667895
np: 0.00011947099847020581
mt: 0.0020772050047526136
********************
sample size: 10000
********************
np: 0.00024697799381101504
mt: 0.002037393998762127
np: 0.0002713389985729009
mt: 0.0020942929986631498
********************
sample size: 100000
********************
np: 0.0007130410012905486
mt: 0.0019842900001094677
np: 0.0007540129954577424
mt: 0.0029724110063398257
********************
sample size: 1000000
********************
np: 0.0094779249993735
mt: 0.007134920000680722
np: 0.009129883001151029
mt: 0.012836456997320056
********************
sample size: 10000000
********************
np: 0.09471094200125663
mt: 0.0453535050037317
np: 0.09436299200024223
mt: 0.04188535599678289
********************
sample size: 100000000
********************
np: 0.9537652180006262
mt: 0.3957935369980987
np: 0.9624398809974082
mt: 0.4019058070043684
Они гораздо менее обнадеживающие, чем результаты, указанные ранее в теме, где указано ускорение примерно в 3,5 раза, и не было многопоточности. Достигнутые мною результаты довольно разумны, я ожидал бы, что накладные расходы на многопоточность будут доминировать во времени, пока массивы не станут очень большими, и в этот момент увеличение производительности начнет приближаться к std::thread::hardware_concurrency x увеличению.
Кажется, что для некоторого кода NumPy есть место для оптимизации под конкретные приложения, особенно в том, что касается многопоточности. Мне не ясно, стоит ли это усилий или нет, но это, безусловно, похоже на хорошее упражнение (или что-то в этом роде). Я думаю, что, возможно, изучение некоторых из этих "сторонних инструментов", таких как Cython, может лучше использовать время, но кто знает.