NumPy: berfungsi untuk maks () dan min () secara bersamaan

109

numpy.amax () akan menemukan nilai maksimal dalam sebuah array, dan numpy.amin () melakukan hal yang sama untuk nilai min. Jika saya ingin menemukan max dan min, saya harus memanggil kedua fungsi tersebut, yang mengharuskan melewati dua kali array (sangat besar), yang tampaknya lambat.

Apakah ada fungsi dalam numpy API yang menemukan max dan min hanya dengan satu kali pass melalui data?

Stuart Berg
sumber
1
Seberapa besar sangat besar? Jika saya punya waktu, saya akan menjalankan beberapa tes yang membandingkan implementasi fortran dengan amaxdanamin
mgilson
1
Saya akui "sangat besar" itu subjektif. Dalam kasus saya, saya berbicara tentang array yang berukuran beberapa GB.
Stuart Berg
itu cukup besar. Saya telah membuat kode contoh untuk menghitungnya di fortran (bahkan jika Anda tidak tahu fortran, seharusnya cukup mudah untuk memahami kodenya). Benar-benar membuat perbedaan menjalankannya dari fortran vs. berlari melalui numpy. (Agaknya, Anda harus bisa mendapatkan kinerja yang sama dari C ...) Saya tidak yakin - Saya kira kita memerlukan dev yang numpy untuk mengomentari mengapa fungsi saya bekerja jauh lebih baik daripada fungsi mereka ...
mgilson
Tentu saja, ini bukanlah ide baru. Misalnya, pustaka boost minmax (C ++) menyediakan implementasi algoritme yang saya cari.
Stuart Berg
3
Bukan jawaban atas pertanyaan yang diajukan, tetapi mungkin menarik bagi orang-orang di utas ini. Menanyakan NumPy tentang menambahkan minmaxke pustaka yang dipermasalahkan ( github.com/numpy/numpy/issues/9836 ).
jakirkham

Jawaban:

49

Apakah ada fungsi dalam numpy API yang menemukan max dan min hanya dengan satu kali pass melalui data?

Tidak. Pada saat tulisan ini dibuat, belum ada fungsi seperti itu. (Dan ya, jika ada yang fungsi seperti, kinerjanya akan secara signifikan lebih baik daripada menelepon numpy.amin()dan numpy.amax()berturut-turut pada array besar.)

Stuart Berg
sumber
31

Saya tidak berpikir bahwa melewati array dua kali adalah masalah. Pertimbangkan pseudo-code berikut:

minval = array[0]
maxval = array[0]
for i in array:
    if i < minval:
       minval = i
    if i > maxval:
       maxval = i

Meskipun hanya ada 1 loop di sini, masih ada 2 pemeriksaan. (Alih-alih memiliki 2 loop dengan masing-masing 1 centang). Sungguh satu-satunya hal yang Anda simpan adalah overhead 1 loop. Jika array benar-benar besar seperti yang Anda katakan, overhead itu kecil dibandingkan dengan beban kerja loop yang sebenarnya. (Perhatikan bahwa ini semua diimplementasikan di C, jadi loop lebih atau kurang gratis).


EDIT Maaf untuk 4 dari Anda yang memberikan suara positif dan percaya pada saya. Anda pasti bisa mengoptimalkan ini.

Berikut beberapa kode fortran yang dapat dikompilasi menjadi modul python melalui f2py(mungkin seorang Cythonguru dapat datang dan membandingkannya dengan versi C yang dioptimalkan ...):

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

Kompilasi melalui:

f2py -m untitled -c fortran_code.f90

Dan sekarang kami berada di tempat di mana kami dapat mengujinya:

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'

Hasilnya agak mengejutkan bagi saya:

8.61869883537 # numpy min/max
1.60417699814 # minmax1
2.30169081688 # minmax2

Saya harus mengatakan, saya tidak sepenuhnya memahaminya. Membandingkan hanya np.minversus minmax1dan minmax2masih merupakan pertarungan yang kalah, jadi ini bukan hanya masalah memori ...

catatan - Meningkatkan ukuran dengan faktor 10**adan mengurangi pengulangan dengan faktor 10**a(menjaga ukuran masalah konstan) memang mengubah kinerja, tetapi tidak dengan cara yang tampaknya konsisten yang menunjukkan bahwa ada beberapa interaksi antara kinerja memori dan overhead panggilan fungsi di python. Bahkan membandingkan minimplementasi sederhana di fortran beats numpy's dengan faktor sekitar 2 ...

mgilson.dll
sumber
21
Keuntungan dari single pass adalah efisiensi memori. Terutama jika array Anda cukup besar untuk ditukar, ini bisa sangat besar.
Dougal
4
Itu tidak sepenuhnya benar, ini hampir setengahnya lebih cepat, karena dengan array semacam ini, kecepatan memori biasanya menjadi faktor pembatas, jadi bisa setengahnya lebih cepat ...
seberg
3
Anda tidak selalu membutuhkan dua pemeriksaan. Jika i < minvalbenar, maka i > maxvalselalu salah, jadi Anda hanya perlu melakukan 1,5 pemeriksaan rata-rata per iterasi ketika detik ifdiganti dengan elif.
Fred Foo
2
Catatan kecil: Saya ragu Cython adalah cara untuk mendapatkan modul C yang paling optimal yang dapat dipanggil Python. Tujuan Cython adalah menjadi semacam Python beranotasi tipe, yang kemudian diterjemahkan mesin ke C, sedangkan f2pyhanya membungkus Fortran dengan kode tangan sehingga dapat dipanggil oleh Python. Tes yang "lebih adil" mungkin adalah C coding tangan dan kemudian menggunakan f2py(!) Untuk membungkusnya untuk Python. Jika Anda mengizinkan C ++, maka Shed Skin mungkin merupakan tempat yang tepat untuk menyeimbangkan kemudahan pengkodean dengan kinerja.
John Y
4
pada numpy 1,8 menit dan maks di-vectorisasi pada platform amd64, pada core2duo numpy saya bekerja sebaik kode fortran ini. Tetapi satu pass akan menguntungkan jika array melebihi ukuran cache cpu yang lebih besar.
jtaylor
23

Ada fungsi untuk mencari (max-min) yang disebut numpy.ptp jika itu berguna untuk Anda:

>>> import numpy
>>> x = numpy.array([1,2,3,4,5,6])
>>> x.ptp()
5

tapi saya rasa tidak ada cara untuk menemukan min dan max dengan satu traversal.

EDIT: ptp hanya memanggil min dan max di bawah tenda

jterrace
sumber
2
Ini menjengkelkan karena kiranya cara ptp diimplementasikan harus melacak max dan min!
Andy Hayden
1
Atau mungkin hanya menelepon max dan min, tidak yakin
jterrace
3
@hayden ternyata ptp hanya menelepon max dan min
jterrace
1
Itu adalah kode array bertopeng; kode ndarray utama ada di C. Tapi ternyata kode C juga melakukan iterasi pada array dua kali: github.com/numpy/numpy/blob/… .
Ken Arnold
20

Anda dapat menggunakan Numba , yang merupakan kompiler Python dinamis yang sadar NumPy menggunakan LLVM. Implementasi yang dihasilkan cukup sederhana dan jelas:

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()))

Ini juga harus lebih cepat daripada min() & max()implementasi Numpy . Dan semuanya tanpa harus menulis satu baris kode C / Fortran.

Lakukan tes kinerja Anda sendiri, karena itu selalu bergantung pada arsitektur Anda, data Anda, versi paket Anda ...

Peque
sumber
2
> Ini juga harus lebih cepat dari implementasi min () & max () Numpy. Saya rasa ini tidak benar. numpy bukan python asli - itu C. `` x = numpy.random.rand (10000000) t = time () untuk i dalam range (1000): minmax (x) print ('numba', time () - t) t = time () untuk i dalam range (1000): x.min () x.max () print ('numpy', time () - t) `` Hasil di: ('numba', 10.299750089645386 ) ('numpy', 9.898081064224243)
Authman Apatira
1
@AuthmanApatira: Ya, tolok ukur selalu seperti itu, itulah mengapa saya mengatakan " harus " (lebih cepat) dan " lakukan uji kinerja Anda sendiri, karena selalu bergantung pada arsitektur Anda, data Anda ... ". Dalam kasus saya, saya mencoba dengan 3 komputer dan mendapatkan hasil yang sama (Numba lebih cepat dari Numpy), tetapi hasil di komputer Anda mungkin berbeda ... Apakah Anda mencoba menjalankan numbafungsi sekali sebelum benchmark untuk memastikan itu dikompilasi JIT ?. Juga, jika Anda menggunakan ipython, untuk kesederhanaan, saya akan menyarankan Anda menggunakan %timeit whatever_code()untuk mengukur eksekusi waktu.
Peque
3
@AuthmanApatira: Bagaimanapun yang saya coba tunjukkan dengan jawaban ini adalah bahwa terkadang kode Python (dalam hal ini JIT-dikompilasi dengan Numba) bisa secepat perpustakaan terkompilasi C tercepat (setidaknya kita berbicara tentang urutan yang sama besarnya), yang sangat mengesankan mengingat kami tidak menulis apa pun selain kode Python murni, bukankah Anda setuju? ^^
Peque
Saya setuju =) Juga, terima kasih atas tip di komentar sebelumnya tentang jupyter dan kompilasi fungsi sekali di luar kode waktu.
Authman Apatira
1
Hanya berlari melintasi ini, bukan itu penting dalam kasus praktis, tetapi elifmemungkinkan minimum Anda menjadi lebih besar dari maks. Misalnya, dengan larik dengan panjang 1, nilai maks adalah berapa pun nilainya, sedangkan min adalah + tak terhingga. Bukan masalah besar untuk satu kali saja, tapi bukan kode yang baik untuk dimasukkan jauh ke dalam perut monster produksi.
Mike Williamson
12

Secara umum, Anda dapat mengurangi jumlah perbandingan untuk algoritme minmax dengan memproses dua elemen sekaligus dan hanya membandingkan yang lebih kecil ke minimum sementara dan yang lebih besar dengan maksimum sementara. Rata-rata seseorang hanya membutuhkan 3/4 dari perbandingan daripada pendekatan yang naif.

Ini dapat diimplementasikan dalam c atau fortran (atau bahasa tingkat rendah lainnya) dan hampir tidak terkalahkan dalam hal kinerja. saya menggunakan untuk menggambarkan prinsip dan mendapatkan implementasi yang sangat cepat, tipe-independen:

import numba as nb
import numpy as np

@nb.njit
def minmax(array):
    # Ravel the array and return early if it's 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

Ini jelas lebih cepat daripada pendekatan naif yang disajikan 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

Seperti yang diharapkan, implementasi minmax baru hanya membutuhkan sekitar 3/4 dari waktu implementasi naif ( 2.1 / 2.75 = 0.7636363636363637)

MSeifert
sumber
1
Di komputer saya, solusi Anda tidak lebih cepat dari solusi Peque. Numba 0.33.
John Zwinck
@johnzwinck apakah Anda menjalankan benchmark dalam jawaban saya yang berbeda? Jika demikian, bisakah Anda membagikannya? Tapi itu mungkin: saya melihat beberapa regresi di versi yang lebih baru juga.
MSeifert
Saya menjalankan patokan Anda. Pengaturan waktu solusi Anda dan @Peque hampir sama (~ 2,8 md).
John Zwinck
@JohnZwinck Aneh, saya baru saja mengujinya lagi dan di komputer saya sudah pasti lebih cepat. Mungkin ada hubungannya dengan numba dan LLVM yang bergantung pada perangkat kerasnya.
MSeifert
Saya mencoba di komputer lain sekarang (workstation besar) dan mendapatkan 2,4 ms untuk komputer Anda vs 2,6 untuk Peque. Jadi, kemenangan kecil.
John Zwinck
11

Hanya untuk mendapatkan beberapa ide tentang angka yang diharapkan, dengan pendekatan berikut:

import numpy as np


def extrema_np(arr):
    return np.max(arr), np.min(arr)
import numba as nb


@nb.jit(nopython=True)
def extrema_loop_nb(arr):
    n = arr.size
    max_val = min_val = arr[0]
    for i in range(1, n):
        item = arr[i]
        if item > max_val:
            max_val = item
        elif item < min_val:
            min_val = item
    return max_val, min_val
import numba as nb


@nb.jit(nopython=True)
def extrema_while_nb(arr):
    n = arr.size
    odd = n % 2
    if not odd:
        n -= 1
    max_val = min_val = arr[0]
    i = 1
    while i < n:
        x = arr[i]
        y = arr[i + 1]
        if x > y:
            x, y = y, x
        min_val = min(x, min_val)
        max_val = max(y, max_val)
        i += 2
    if not odd:
        x = arr[n]
        min_val = min(x, min_val)
        max_val = max(x, max_val)
    return max_val, min_val
%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True


import numpy as np


cdef void _extrema_loop_cy(
        long[:] arr,
        size_t n,
        long[:] result):
    cdef size_t i
    cdef long item, max_val, min_val
    max_val = arr[0]
    min_val = arr[0]
    for i in range(1, n):
        item = arr[i]
        if item > max_val:
            max_val = item
        elif item < min_val:
            min_val = item
    result[0] = max_val
    result[1] = min_val


def extrema_loop_cy(arr):
    result = np.zeros(2, dtype=arr.dtype)
    _extrema_loop_cy(arr, arr.size, result)
    return result[0], result[1]
%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True


import numpy as np


cdef void _extrema_while_cy(
        long[:] arr,
        size_t n,
        long[:] result):
    cdef size_t i, odd
    cdef long x, y, max_val, min_val
    max_val = arr[0]
    min_val = arr[0]
    odd = n % 2
    if not odd:
        n -= 1
    max_val = min_val = arr[0]
    i = 1
    while i < n:
        x = arr[i]
        y = arr[i + 1]
        if x > y:
            x, y = y, x
        min_val = min(x, min_val)
        max_val = max(y, max_val)
        i += 2
    if not odd:
        x = arr[n]
        min_val = min(x, min_val)
        max_val = max(x, max_val)
    result[0] = max_val
    result[1] = min_val


def extrema_while_cy(arr):
    result = np.zeros(2, dtype=arr.dtype)
    _extrema_while_cy(arr, arr.size, result)
    return result[0], result[1]

( extrema_loop_*()pendekatannya mirip dengan yang diusulkan di sini , sedangkan extrema_while_*()pendekatan didasarkan pada kode dari sini )

Pengaturan waktu berikut:

bm

menunjukkan bahwa extrema_while_*()yang tercepat, dengan extrema_while_nb()yang tercepat. Bagaimanapun, solusi extrema_loop_nb()dan extrema_loop_cy()juga mengungguli pendekatan NumPy saja (menggunakan np.max()dan np.min()secara terpisah).

Terakhir, perhatikan bahwa tidak ada yang sefleksibel np.min()/ np.max()(dalam hal dukungan n-dim, axisparameter, dll.).

(kode lengkap tersedia di sini )

norok2
sumber
2
Sepertinya Anda bisa mendapatkan kecepatan ekstra 10% jika menggunakan @njit (fastmath = True)extrema_while_nb
argenisleon
10

Tidak ada yang menyebutkan numpy.percentile , jadi saya pikir saya akan melakukannya. Jika Anda meminta [0, 100]persentil, itu akan memberi Anda larik dua elemen, min (persentil ke-0) dan maks (persentil ke-100).

Namun, itu tidak memenuhi tujuan OP: itu tidak lebih cepat dari min dan max secara terpisah. Itu mungkin karena beberapa mesin yang memungkinkan persentil non-ekstrim (masalah yang lebih sulit, yang seharusnya membutuhkan waktu lebih lama).

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'

Versi Numpy yang akan datang dapat dimasukkan ke dalam kasus khusus untuk melewati penghitungan persentil normal jika hanya [0, 100]diminta. Tanpa menambahkan apa pun ke antarmuka, ada cara untuk meminta Numpy min dan max dalam satu panggilan (bertentangan dengan apa yang dikatakan dalam jawaban yang diterima), tetapi implementasi standar pustaka tidak memanfaatkan kasus ini untuk membuatnya bermanfaat.

Jim Pivarski
sumber
9

Ini adalah utas lama, tapi bagaimanapun, jika ada yang melihat ini lagi ...

Saat mencari min dan max secara bersamaan, adalah mungkin untuk mengurangi jumlah perbandingan. Jika float yang Anda bandingkan (yang menurut saya memang demikian), ini mungkin menghemat waktu Anda, meskipun bukan kompleksitas komputasi.

Alih-alih (kode 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]

Anda dapat membandingkan dua nilai yang berdekatan dalam larik terlebih dahulu, lalu hanya membandingkan nilai yang lebih kecil dengan nilai minimum saat ini, dan nilai yang lebih besar dengan nilai maksimum saat ini:

## 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

Kode di sini ditulis dengan Python, jelas untuk kecepatan Anda akan menggunakan C atau Fortran atau Cython, tetapi dengan cara ini Anda melakukan 3 perbandingan per iterasi, dengan iterasi len (ar) / 2, memberikan perbandingan 3/2 * len (ar). Berbeda dengan itu, melakukan perbandingan "dengan cara yang jelas" Anda melakukan dua perbandingan per iterasi, yang mengarah ke perbandingan 2 * len (ar). Menghemat 25% waktu perbandingan.

Mungkin seseorang suatu hari akan menganggap ini berguna.

Bennet
sumber
6
sudahkah Anda membandingkan ini? pada perangkat keras x86 modern Anda memiliki instruksi mesin untuk min dan max seperti yang digunakan dalam varian pertama, ini menghindari kebutuhan akan cabang sementara kode Anda menempatkan ketergantungan kontrol yang mungkin tidak dipetakan dengan baik ke perangkat keras.
jtaylor
Sebenarnya aku belum. Akan dilakukan jika saya mendapat kesempatan. Saya pikir cukup jelas bahwa kode python murni akan kehilangan implementasi terkompilasi yang masuk akal, tetapi saya bertanya-tanya apakah percepatan dapat dilihat di Cython ...
Bennet
13
Ada implementasi minmax di numpy, under the hood, digunakan oleh np.bincount, lihat di sini . Itu tidak menggunakan trik yang Anda tunjukkan, karena ternyata hingga 2x lebih lambat dari pendekatan naif. Ada tautan dari PR ke beberapa tolok ukur komprehensif dari kedua metode tersebut.
Jaime
5

Pada pandangan pertama, tampaknya untuk melakukan trik:numpy.histogram

count, (amin, amax) = numpy.histogram(a, bins=1)

... tetapi jika Anda melihat sumber untuk fungsi itu, itu hanya memanggil a.min()dan a.max()secara independen, dan karena itu gagal untuk menghindari masalah kinerja yang dibahas dalam pertanyaan ini. :-(

Demikian pula, scipy.ndimage.measurements.extrematampak seperti kemungkinan, tetapi itu juga hanya panggilan a.min()dan a.max()mandiri.

Stuart Berg
sumber
3
np.histogramtidak selalu berfungsi untuk ini karena nilai yang dikembalikan (amin, amax)adalah untuk nilai minimum dan maksimum nampan. Jika saya memiliki, misalnya a = np.zeros(10), np.histogram(a, bins=1)pengembalian (array([10]), array([-0.5, 0.5])). Pengguna mencari (amin, amax)= (0, 0) dalam kasus itu.
eclark
3

Itu sepadan dengan usaha saya, jadi saya akan mengusulkan solusi paling sulit dan paling tidak elegan di sini untuk siapa pun yang mungkin tertarik. Solusi saya adalah mengimplementasikan multi-threaded min-max dalam algoritma one pass di C ++, dan menggunakannya untuk membuat modul ekstensi Python. Upaya ini membutuhkan sedikit overhead untuk mempelajari cara menggunakan Python dan NumPy C / C ++ API, dan di sini saya akan menunjukkan kode dan memberikan beberapa penjelasan dan referensi kecil untuk siapa pun yang ingin menempuh jalur ini.

Min / Maks multi-utas

Tidak ada yang terlalu menarik di sini. Array dipecah menjadi potongan-potongan ukuran length / workers. Min / max dihitung untuk setiap potongan di a future, yang kemudian dipindai untuk min / max global.

    // 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

Modul Ekstensi Python

Di sinilah segalanya mulai menjadi jelek ... Salah satu cara untuk menggunakan kode C ++ dengan Python adalah dengan menerapkan modul ekstensi. Modul ini dapat dibangun dan dipasang menggunakan distutils.coremodul standar. Penjelasan lengkap tentang apa yang diperlukan tercakup dalam dokumentasi Python: https://docs.python.org/3/extending/extending.html . CATATAN: tentunya ada cara lain untuk mendapatkan hasil yang serupa, dengan mengutip https://docs.python.org/3/extending/index.html#extending-index :

Panduan ini hanya mencakup alat dasar untuk membuat ekstensi yang disediakan sebagai bagian dari versi CPython. Alat pihak ketiga seperti Cython, cffi, SWIG dan Numba menawarkan pendekatan yang lebih sederhana dan lebih canggih untuk membuat ekstensi C dan C ++ untuk Python.

Pada dasarnya, rute ini mungkin lebih bersifat akademis daripada praktis. Dengan itu, apa yang saya lakukan selanjutnya adalah, menempel cukup dekat dengan tutorial, membuat file modul. Ini pada dasarnya adalah boilerplate untuk distutils untuk mengetahui apa yang harus dilakukan dengan kode Anda dan membuat modul Python darinya. Sebelum melakukan semua ini, mungkin bijaksana untuk membuat lingkungan virtual Python sehingga Anda tidak mencemari paket sistem Anda (lihat https://docs.python.org/3/library/venv.html#module-venv ).

Ini file modulnya:

// 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); }

Dalam file ini terdapat penggunaan signifikan dari Python serta NumPy API, untuk informasi lebih lanjut lihat: https://docs.python.org/3/c-api/arg.html#c.PyArg_ParseTuple , dan untuk NumPy : https://docs.scipy.org/doc/numpy/reference/c-api.array.html .

Memasang Modul

Hal berikutnya yang harus dilakukan adalah memanfaatkan distutils untuk menginstal modul. Ini membutuhkan file setup:

# 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])

Untuk akhirnya menginstal modul, jalankan python3 setup.py installdari lingkungan virtual Anda.

Menguji Modul

Terakhir, kita dapat menguji untuk melihat apakah implementasi C ++ benar-benar mengungguli penggunaan NumPy yang naif. Untuk melakukannya, berikut ini skrip pengujian sederhana:

# 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))

Inilah hasil yang saya dapat dari melakukan semua ini:

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  

Ini jauh kurang menggembirakan daripada hasil yang ditunjukkan sebelumnya di utas, yang menunjukkan sekitar 3,5x percepatan, dan tidak menyertakan multi-threading. Hasil yang saya capai agak masuk akal, saya berharap bahwa overhead threading dan akan mendominasi waktu sampai array menjadi sangat besar, di mana peningkatan kinerja akan mulai mendekati std::thread::hardware_concurrencypeningkatan x.

Kesimpulan

Jelas ada ruang untuk pengoptimalan khusus aplikasi untuk beberapa kode NumPy, tampaknya, khususnya yang berkaitan dengan multi-threading. Apakah itu sepadan atau tidak, tidak jelas bagi saya, tetapi itu jelas terlihat seperti latihan yang baik (atau sesuatu). Saya pikir mungkin mempelajari beberapa "alat pihak ketiga" seperti Cython mungkin penggunaan waktu yang lebih baik, tapi siapa tahu.

Nathan Chappell
sumber
1
Saya mulai mempelajari kode Anda, mengetahui beberapa C ++ tetapi masih belum menggunakan std :: future dan std :: async. Pada fungsi template 'min_max_mt' Anda, bagaimana cara mengetahui bahwa setiap pekerja telah selesai antara mengaktifkan dan mengambil hasil? (Meminta hanya untuk memahami, tidak mengatakan bahwa ada yang salah dengan itu)
ChrCury78
Garis v = min_max_it->get();. The getMetode blok sampai hasilnya siap dan kembali itu. Karena perulangan melewati setiap masa depan, itu tidak akan selesai sampai semuanya selesai. future.get ()
Nathan Chappell
0

Cara terpendek yang saya temukan adalah ini:

mn, mx = np.sort(ar)[[0, -1]]

Tapi karena itu mengurutkan array, itu bukan yang paling efisien.

Cara singkat lainnya adalah:

mn, mx = np.percentile(ar, [0, 100])

Ini seharusnya lebih efisien, tetapi hasilnya dihitung, dan float dikembalikan.

Israel Unterman
sumber
Sayangnya, keduanya adalah solusi paling lambat dibandingkan dengan yang lain di halaman ini: m = np.min (a); M = np.max (a) -> 0,54002 ||i> m, M = f90_minmax1 (a) -> 0,72134 ||i> m, M = numba_minmax (a) -> 0,77323 ||i> m, M = np.sort (a) [[0, -1]] -> 12.01456 ||i> m, M = np. persentil (a, [0, 100]) -> 11.09418 ||i> dalam detik untuk 10.000 pengulangan untuk larik 100k elemen
Isaías