Python Numpyの代替としてのD std.ndslice

はじめに： 私は Pythonで6年以上執筆しており、自分をこの言語の専門家と呼ぶことができます。 最近、私は彼についての本さえ書いた。 ただし、過去8か月間Dに切り替え、4か月間、標準のPhobosライブラリの拡張に関するこの言語の開発に積極的に関与しました。 また、議論されるstd.ndsliceモジュールのコードレビューにも参加しました。

Numpyのようなstd.ndsliceは、多次元配列で動作するように設計されています。 ただし、Numpyとは異なり、ndsliceは通常のライブラリのあらゆる場所で使用される範囲（範囲）に基づいているため、オーバーヘッドが非常に低くなっています。 範囲を使用すると、不必要なコピー手順を回避でき、遅延計算を美しく整理できます。

この記事では、stump.ndsliceがNumpyと比較してどのような利点があるかについて説明します。

なぜPythonプログラマーはDに目を向けるのですか？

Dを使用すると、Pythonとほぼ同じシンプルで明快なコードを記述できますが、このコードはPythonコードよりもはるかに高速に動作します。

次の例では、 iota関数（Pythonのアナログはxrangeと呼ばれxrange ）を使用して0から999までの数値の範囲を作成し、5x5x40の次元の配列を返します。

 import std.range : iota; import std.experimental.ndslice; void main() { auto slice = sliced(iota(1000), 5, 5, 40); } 

Dは静的に型付けされた言語であり、型を明示的に指定するとコードの可読性が向上するため、Pythonプログラマーはautoキーワードを使用して自動型付けを使用する方が簡単です。

sliced関数は、多次元スライスを返します。 入力では、単純な配列とranges両方を受け入れることができranges 。 その結果、0〜999の数字で構成される5x5x40のキューブが得られました。

Rangesとは何かという言葉。 それらを範囲としてロシア語に翻訳する方がより正確です。 範囲を使用すると、クラスまたは構造にかかわらず、あらゆるタイプのデータを列挙するためのルールを記述することができます。 関数を実装するだけで十分ですpopFrontは、範囲の最初の要素popFrontを返しますpopFrontは、次の要素に移動し、空です。検索対象のシーケンスが空であることを示すブール値empty返します。 Ranges使用すると、必要に応じてデータにアクセスして、遅延反復を実行できます。 Rangesの概念については、 ここで詳しく説明します 。

空のメモリ割り当てがないことに注意してください！ これは、 iotaでも遅延範囲を生成できるため、遅延モードでslicedすると、 iotaからデータを取得して、到着時に処理するためです。

ご覧のとおり、std.ndsliceの動作はNumpyとは少し異なります。 Numpyはデータ用に独自のタイプを作成し、std.ndsliceは既存のデータを操作する方法を提供します。 これにより、無駄なメモリ割り当てにリソースを浪費することなく、プログラム全体で同じデータを使用できます！ これがソリューションのパフォーマンスに非常に深刻な影響を与えることは容易に推測できます。

次の例を見てみましょう。 その中で、 stdinからデータを受け取り、一意の行のみをフィルター処理し、並べ替えてからstdout戻ります。

 import std.stdio; import std.array; import std.algorithm; void main() { stdin // get stdin as a range .byLine(KeepTerminator.yes) .uniq // stdin is immutable, so we need a copy .map!(a => a.idup) .array .sort // stdout.lockingTextWriter() is an output range, meaning values can be // inserted into to it, which in this case will be sent to stdout .copy(stdout.lockingTextWriter()); } 

遅延範囲の生成をより詳細に処理したい場合は、 この記事を読むことをお勧めします。

sliceは3つの次元があるため、これは範囲の範囲を返す範囲です。 これは、次の例ではっきりとわかります。

 import std.range : iota; import std.stdio : writeln; import std.experimental.ndslice; void main() { auto slice = sliced(iota(1000), 5, 5, 40); foreach (item; slice) { writeln(item); } } 

彼の仕事の結果は次のようになります（省略記号）。

 [[0, 1, ... 38, 39], [40, 41, ... 78, 79], [80, 81, ... 118, 119], [120, 121, ... 158, 159], [160, 161, ... 198, 199]] ... [[800, 801, ... 838, 839], [840, 841, ... 878, 879], [880, 881, ... 918, 919], [920, 921, ... 958, 959], [960, 961, ... 998, 999]] 

foreachは、Pythonのfor foreachとほぼ同じように機能for foreach 。 ただし、Dでは、CスタイルとPythonのループの両方で使用できますが、 enumerateやxrange手間はxrange 。

UFCS（Uniform Function Call Syntax）を使用すると、次の方法でコードを書き換えることができます。

 import std.range : iota; import std.experimental.ndslice; void main() { auto slice = 1000.iota.sliced(5, 5, 40); } 

UFCSを使用すると、メソッド呼び出しをチェーンで記録し、次のように記述できます 。

 a.func(b) 

代わりに：

 func(a, b) 

ダブパッケージマネージャーを使用して空のプロジェクトを生成しましょう。 dub initコマンドおよび\source\app.d次のように記述します。

 import std.experimental.ndslice; void main() { } 

現在std.experimental.ndslice; std.experimentalセクションにあります。 これは、それが生であることを意味しません。 これは、彼が主張する必要があることを意味します。

次のコマンドでプロジェクトを組み立てます。

 dub 

D ndsliceモジュールはNumpyに非常に似ています：

 a = numpy.arange(1000).reshape((5, 5, 40)) b = a[2:-1, 1, 10:20] 

同等：

 auto a = 1000.iota.sliced(5, 5, 40); auto b = a[2 .. $, 1, 10 .. 20]; 

次に、2次元配列を作成して、各列の中央を取得します。

Python：

 import numpy data = numpy.arange(100000).reshape((100, 1000)) means = numpy.mean(data, axis=0) 

D：

 import std.range; import std.algorithm.iteration; import std.experimental.ndslice; import std.array : array; void main() { auto means = 100_000.iota .sliced(100, 1000) .transposed .map!(r => sum(r) / r.length) .array; } 

このコードが遅延モードで動作しないようにするには、 arrayメソッドを呼び出す必要がありました。 ただし、実際のアプリケーションでは、結果はプログラムの別の部分で使用されている間は計算されません。

現在、Phobosには組み込みの統計モジュールはありません。 したがって、この例では単純なラムダを使用して平均値を見つけます。 map!機能map! 最後に感嘆符が付いています。 これは、これがテンプレート関数であることを意味します。 本体で指定された式に基づいて、コンパイル段階でコードを生成できます。 Dでテンプレート自体がどのように機能するかについての良い記事があります。

DのコードはPythonよりも少し冗長であることがわかりましたが、 map!おかげですmap! コードは範囲であるすべての入力で機能します。 PythonコードはNumpyの特別な配列でのみ機能します。

ここで、このテストの後、Pythonが時々Dを失い、Hacker Newsで多くの議論を行った後、間違いを犯し、比較が完全に正しくないことに気付いたと言わなければなりません。 iotaは、 sliced関数sliced取得するデータを動的に作成します。 したがって、最後の再配置の瞬間までメモリに触れません。 Dは、データ型がlong配列も返しますが、Numpyはdoubleから返します。 その結果、コードを書き直して、配列の値を10000ではなく1000 000にしました。次に何が起こったかを示します。

 import std.range : iota; import std.array : array; import std.algorithm; import std.datetime; import std.conv : to; import std.stdio; import std.experimental.ndslice; enum test_count = 10_000; double[] means; int[] data; void f0() { means = data .sliced(100, 10000) .transposed .map!(r => sum(r, 0L) / cast(double) r.length) .array; } void main() { data = 1_000_000.iota.array; auto r = benchmark!(f0)(test_count); auto f0Result = to!Duration(r[0] / test_count); f0Result.writeln; } 

2.9 GHz Intel Core Broadwell i5を搭載した2015 MacBook Proで実施したテスト。 Pythonでは、速度を測定するためにD std.datetime.benchmark %timeit関数を使用しました。 LDC v0.17で次のキーを使用してすべてをコンパイルしました： ldmd2 -release -inline -boundscheck=off -O または、ダブを使用する場合、オプションdub --build=release-nobounds --compiler=ldmd2はこれらのキーに類似しています。

最初のテストの結果は次のとおりです。

 Python: 145 µs LDC: 5 µs D is 29x faster 

修正バージョンのテストは次のとおりです。

 Python: 1.43 msec LDC: 628 μs D is 2.27x faster 

NumpyはCで書かれており、Dでは誰もがガベージコレクターをscるので、それは悪い違いではないことに同意します。

DはどのようにNumpyの問題を回避しますか

はい、Numpyは高速ですが、単純なPython配列と比較した場合のみ高速です。 しかし、問題は、これらの単純な配列と部分的にしか互換性がないことです。

Numpyライブラリは、Pythonの残りの部分のどこかにあります。 彼女は自分の人生を生きています。 独自の機能を使用し、そのタイプで動作します。 たとえば、PythonでNumPyで作成された配列を使用する必要がある場合、新しい変数にコピーするnp.asarrayを使用する必要があります。 迅速なgithubの検索により、何千ものプロジェクトがこの松葉杖を使用していることが明らかになりました。 データは、これらの空のコピーがなくても、ある機能から別の機能に簡単に転送できます。

 import numpy as np a = [[0.2,0.5,0.3], [0.2,0.5,0.3]] p = np.asarray(a) y = np.asarray([0,1]) 

彼らは、Python標準ライブラリの一部を書き換えてNumpy型を使用することで、この問題を解決しようとしています。 ただし、これはまだ完全な解決策ではないため、作成時に素晴らしいジョークが発生します。

 sum(a) 

代わりに：

 a.sum() 

速度が10倍低下します。

Dは、単に設計上このような問題を抱えていません。 これは、コンパイル済みの静的に型指定された言語です。 コード生成中に、すべてのタイプの変数が認識されます。 std.ndslice自体では、たとえばstd.algorithmやstd.rangeなどのすばらしい機能など、Phobosライブラリのすべての機能に完全にアクセスできます。 そう、コンパイル段階でDテンプレートを使用してコードを生成できます。

以下に例を示します。

 import std.range : iota; import std.algorithm.iteration : sum; import std.experimental.ndslice; void main() { auto slice = 1000.iota.sliced(5, 5, 40); auto result = slice // sum expects an input range of numerical values, so to get one // we call std.experimental.ndslice.byElement to get the unwound // range .byElement .sum; } 

あなたは、 sum関数を取り、使用するだけで、基本ライブラリの他の関数と同様に、取り、機能します。

Python自体で、特定の値で初期化された定義済みの長さのリストを取得するには、次のように記述する必要があります。

 a = [0] * 1000 

Numpyは完全に異なります。

 a = numpy.zeros((1000)) 

また、Numpyを使用しない場合、メモリを消費する不要なコピー操作よりもコードが4倍遅くなります。 範囲はDの助けになります。これにより、空のコピー操作をせずに、同じ操作をすばやく実行できます。

 auto a = repeat(0, 1000).array; 

また、必要に応じて、すぐにndsliceを呼び出すことができます。

 auto a = repeat(0, 1000).array.sliced(5, 5, 40); 

現在のNumpyの主な利点は、その普及率です。 現在では、銀行システムから機械学習まで非常に広く使用されています。 多くの本、例、記事があります。 ただし、Dの数学的可能性は明らかにすぐに拡張されるでしょう。 そのため、ndsliceの作者は、現在phosのBLAS（Basic Linear Algebra Subprograms）に取り組んでおり、これもndsliceおよび標準ライブラリと完全に統合されると述べています。

強力な数学的サブシステムにより、ビッグデータを扱う必要がある多くの問題を非常に効率的に解決できます。 たとえば、コンピュータービジョンシステム。 これらのシステムの1つのプロトタイプはすでに開発中で、 DCVと呼ばれています。

Dの画像処理

次の例は、メディアンフィルターが画像のノイズを除去する方法を示しています。 movingWindowByChannel関数は、スライディングウィンドウが必要な他のフィルターでも使用できます。 movingWindowByChannel使用movingWindowByChannelと、スライディングウィンドウを使用して画像内を移動できmovingWindowByChannel 。 このような各ウィンドウは、選択したゾーンに基づいてピクセル値を計算するフィルターに渡されます。

この関数は、部分的に重複した領域を処理しません。 ただし、その助けを借りて、それらの値も計算できます。 これを行うには、元の画像の境界線を反映したエッジを持つ拡大画像を作成してから処理します。

 /** Params: filter = unary function. Dimension window 2D is the argument. image = image dimensions `(h, w, c)`, where  is the number of channels in the image nr = number of rows in the window n = number of columns in the window Returns: image dimensions `(h - nr + 1, w - nc + 1, c)`, where  is the number of channels in the image. Dense data layout is guaranteed. */ Slice!(3, C*) movingWindowByChannel(alias filter, C) (Slice!(3, C*) image, size_t nr, size_t nc) { // local imports in D work much like Python's local imports, // meaning if your code never runs this function, these will // never be imported because this function wasn't compiled import std.algorithm.iteration: map; import std.array: array; // 0. 3D // The last dimension represents the color channel. auto wnds = image // 1. 2D composed of 1D // Packs the last dimension. .pack!1 // 2. 2D composed of 2D composed of 1D // Splits image into overlapping windows. .windows(nr, nc) // 3. 5D // Unpacks the windows. .unpack // 4. 5D // Brings the color channel dimension to the third position. .transposed!(0, 1, 4) // 5. 3D Composed of 2D // Packs the last two dimensions. .pack!2; return wnds // 6. Range composed of 2D // Gathers all windows in the range. .byElement // 7. Range composed of pixels // 2D to pixel lazy conversion. .map!filter // 8. `C[]` // The only memory allocation in this function. .array // 9. 3D // Returns slice with corresponding shape. .sliced(wnds.shape); } 

次の関数は、オブジェクトの中央値を計算する方法の例です。 読みやすくするために、この機能は大幅に簡素化されています。

 /** Params: r = input range buf = buffer with length no less than the number of elements in `r` Returns: median value over the range `r` */ T median(Range, T)(Range r, T[] buf) { import std.algorithm.sorting: sort; size_t n; foreach (e; r) { buf[n++] = e; } buf[0 .. n].sort(); immutable m = n >> 1; return n & 1 ? buf[m] : cast(T)((buf[m - 1] + buf[m]) / 2); } 

さて、今メイン自体：

 void main(string[] args) { import std.conv: to; import std.getopt: getopt, defaultGetoptPrinter; import std.path: stripExtension; // In D, getopt is part of the standard library uint nr, nc, def = 3; auto helpInformation = args.getopt( "nr", "number of rows in window, default value is " ~ def.to!string, &nr, "nc", "number of columns in window, default value is equal to nr", &nc); if (helpInformation.helpWanted) { defaultGetoptPrinter( "Usage: median-filter [<options...>] [<file_names...>]\noptions:", helpInformation.options); return; } if (!nr) nr = def; if (!nc) nc = nr; auto buf = new ubyte[nr * nc]; foreach (name; args[1 .. $]) { import imageformats; // can be found at code.dlang.org IFImage image = read_image(name); auto ret = image.pixels .sliced(cast(size_t)image.h, cast(size_t)image.w, cast(size_t)image.c) .movingWindowByChannel !(window => median(window.byElement, buf)) (nr, nc); write_image( name.stripExtension ~ "_filtered.png", ret.length!1, ret.length!0, (&ret[0, 0, 0])[0 .. ret.elementsCount]); } } 

すべての例が明確に見えない場合は、すばらしい本「 Programming in D 」の無料版を読むことをお勧めします。

PSこの出版物を「翻訳」のステータスに翻訳する方法を知っているなら、プライベートで書いてください。