AnacondaのNumPyとPyPIのNumPyの速度を比較する

python

Anaconda Pythonで提供されるNumPyはIntelのMKLを利用しているため高速だという話を聞いたことがあります。実際どの程度違いがあるのか試してみました。

環境構築

実験は、自作PCに入れたUbuntu 16.04で行いました。環境構築にはAnaconda Pythonが提供する仮想環境を用います。

Anaconda版のNumPy

以下のコマンドで環境を作成しました。Cythonは必要ないかもしれませんが、chainerをインストールするときにmkl-random 1.0.1 requires cython, which is not installed. mkl-fft 1.0.2 requires cython, which is not installed.という警告が表示されたので、一応入れています。

$ conda create -n py37_conda_numpy python=3.7 numpy
$ source activate py37_conda_numpy
$ conda install cython
$ pip install chainer

インストールされたパッケージは以下のとおりです。

$ pip freeze
certifi==2018.4.16
chainer==4.3.0
Cython==0.28.3
filelock==3.0.4
mkl-fft==1.0.2
mkl-random==1.0.1
numpy==1.14.5
protobuf==3.6.0
six==1.11.0

PyPI版のNumPy

以下のコマンドで環境を作成しました。

$ conda create -n py37_pip_numpy python=3.7       
$ source activate py37_pip_numpy
$ pip install numpy
$ pip install chainer

インストールされたパッケージは以下のとおりです。

$ pip freeze
certifi==2018.4.16
chainer==4.3.0
filelock==3.0.4
numpy==1.14.5
protobuf==3.6.0
six==1.11.0

実験1. 行列積

1024x1024の行列二つの行列積を100回計算する以下のコードを使って、処理時間を比較しました。

import timeit
import numpy as np

def main():
    HEIGHT = 1024
    WIDTH = 1024
    REPEAT = 100

    np.random.seed(0)
    result = np.zeros((HEIGHT, WIDTH), dtype=np.float32)
    for _ in range(REPEAT):
        a = np.random.rand(HEIGHT, WIDTH)
        b = np.random.rand(HEIGHT, WIDTH)
        result += a @ b
    print(result.sum())

s = timeit.default_timer()
main()
print("elapsed time: {:.2f} sec.".format(timeit.default_timer() - s))

以下に結果を示します。

  • Anaconda版NumPy: 3.35 秒
  • PyPI版NumPy: 8.93 秒

Anaconda版のNumpyがPyPI版Pythonの38%程度の時間で処理を終えており、大差がつきました。プログラム実行中のCPUの負荷をhtopコマンドで眺めると、PyPI版のNumPyは8コアを全消費するのに比べてAnaconda版のNumPyは4コアしか消費しておらず、余裕を感じました。

実験2. ニューラルネット

手書きの数字を０から９に分類する分類器の学習chainer/train_mnist.py at master · chainer/chainer · GitHubを題材に、処理時間を比較しました。この分類器にはCNNは使われておらず、全結合層のみです。以下は実験の詳細です。

  • ソースコードには実験1と同様の時間計測コードを追加した。
  • 学習はCPUのみを用いる。GPUは用いない。
  • 最初の実行時は実験用データセットをダウンロードする処理が走るので、それは処理時間に含めない。
  • 引数に--epoch 5を指定し、5エポックで処理を打ち切る。

以下に結果を示します。

  • Anaconda版NumPy: 54.51 秒
  • PyPI版NumPy: 100.49 秒

Anaconda版のNumpyがPyPI版Pythonの54%程度の時間で処理を終えました。実験1ほどではないですが、やはりAnaconda版のNumpyは高速です。

まとめ

NumPyの速度がボトルネックになっているのであれば、Anaconda版のNumPyの導入は一考の価値があります。

AGC026 B rng_10s

programming contest

おとといのAGCに参加して1問しか解けず。Bに二時間以上考えて解ききれないというのはさすがに問題だと思い、解説を読み込んで自分なりに咀嚼しました。本番で書いていたひどすぎるコードが消えてすっきりしました。しかし、本当に思考に穴がないかはまだ自信がありません。

こういう数学的な問題、いちど袋小路に入ってしまうと思考のループから抜け出せず時間を浪費してしまうことが多く、どうにもいけません。

問題：B - rng_10s 回答：ソースコード 

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# 本番中には解けなかった。解説pdfにある3つの場合分けはできた。
# B個の購入を可能なかぎりしたあとの個数が、「Cより大きくBより小さい」値に
# なることがあるとNoになる、ということには気づいたのだが、
# どう解いてよいか混乱して手がつけられなかった
#
#
# 解説pdf https://img.atcoder.jp/agc026/editorial.pdf を見た
#
# B個の購入を行った直後の個数の遷移を、mod Bの世界で考える。
#   最初は A mod B個。
#   B個の購入を可能な限り行った後も A mod B個で変わらず。
#   D個入荷すると    (A + D) mod B個。
# つまり、入荷するたびに、
# A からスタートして (A + D) mod B, (A + 2D) mod B, ... と遷移していく。
# 
# 遷移を無限に繰り返すと、mod Bの世界でのとりうる値は、g = gcd(B, D)として
#   (A % g), g + (A % g), 2g + (A % g), ..., kg + (A % g) (kは、kg + (A % g)がB未満となる最大の整数) のどれかとなる。
# ここでkを直接求めるのは難しい。
# かわりに、mを整数として、mg + (A % g)が最初にB以上となるのを求めるのは簡単で、それは B + (A % g)。
# これよりひとつ分小さいのは、B + (A % g) - gで、これが kg + (A % g)に等しい。
# 
# (A % g) + B - g が、C を超える場合、次に入荷が行われない。次もBを引くと個数が負になってしまうので No。
# (A % g) + B - g が、C 以下の場合、次に入荷が行われる。Yes。
#
#
# 解説PDFの、「また，(B/g − 1) × inv(D/g, B/g) 回 D を足したときにこの上界が達成できます」
# については理解できていない。
#



import array
from bisect import *
from collections import *
import fractions
import heapq 
from itertools import *
import math
import random
import re
import string
import sys


def gcd(a, b):
    while b:
        a, b = b, a % b
    return a


def solve(A, B, C, D):
    if B > A:
        return "No"
    if B > D:
        return "No"
    if B <= C:
        return "Yes"

    g = gcd(B, D)
    if B - g + (A % g) > C:
        return "No"
    else:
        return "Yes"


T = int(input())
for t in range(T):
    A, B, C, D = map(int, input().split())
    print(solve(A, B, C, D))

自分のブログに付与されたはてなブックマークの総数を公式APIで調べる

hatena

表題のAPIがはてなブックマーク件数取得API - Hatena Developer Centerに掲載されていました。2018年06月05日に追加されたばかりのhttp://api.b.st-hatena.com/entry.countというAPIを使えばよいそうです。

当ブログの場合、api.b.st-hatena.com/entry.total_count?url=http%3A%2F%2Fminus9d.hatenablog.com%2F にアクセスすると以下の結果が得られました (2018-07-04現在）。595件ものブックマークありがとうございます。

{"url":"http://minus9d.hatenablog.com/","total_bookmarks":595}

SRM735 Div. 2 1000 MajoritySubarray

programming contest

Div.2に落とされた直後のSRMで、自己最高の4位を取ることができました。本当は同時開催のTCOに出るつもりだったところ操作ミスでSRMにエントリーしてしまっていたのですが、結果オーライでした。記念のスクリーンショット。 f:id:minus9d:20180628224923p:plain

本番では解けなかった1000のMajoritySubarrayが勉強になったので解法をメモしておきます。

問題概要

  • TopCoder Statistics - Problem Statement
    • 長さN(最大100000)で、M(最大50)以下の整数からなる数列が与えられる。この数列のうち、"majority element”が存在するような部分数列の数を求める。
    • ここで"majority element"とは、数列の過半数を占める数のこと。例えば数列{2, 2, 3, 3, 3}では3がmajority element。

解法

数列の長さが100000なことから、すべての部分文字列でシミュレーションするのは計算量的に無理です。

公式ブログSingle Round Match 735 Editorials - Topcoder 2.0 にわかりやすい解法があります。まず、ある数字に着目して、その数字がmajority elementである部分数列を数えることだけを考えます。数列中に着目する数字が現れれば+1, 現れなければ-1して新たな数列を作成すると、この数列中の任意の2要素の大小関係を見るだけで、その部分数列に"majority element”が存在するかを判定できるようになります。これができれば、BITやセグメントツリーなどのデータ構造を使って、"majority element”が存在する部分数列の数を高速に数えられます。詳しくは以下に貼るソースコード中の説明を見てください。

この「ある数列の部分数列が条件を満たすかどうか、2要素の大小関係を見るだけで判定できるような数列を作れれば、BITで数え上げができる」という部分は他にも応用ができそうだと思いました。

// 公式解答を見た
//   https://www.topcoder.com/blog/single-round-match-735-editorials/

// 配列Aの構成要素はせいぜい50なので、各構成要素ごとに考えればよい。
// 次のような配列Aで、'o'が支配的になる部分配列の数を考える
// 
//   oxxooxo
//   0123456 (idx)
//
// 配列Bを以下のように構成する。
//   B[0]=7  (すべての要素が0以上になるよう下駄を履かせる)
//   a[i]が'o'であればb[i+1]=b[i]+1, 'o'以外であればb[i+1]=b[i]-1
// すると配列Bは以下のようになる
//  78767878
//  01234567 (idx)
// b[j] < b[i]となるj < iのペアを探す。このとき aの部分文字列[j, i)は'o'が支配的になっていることがわかる。
// 列挙すると以下のようになる。
// (j, i) = (0, 1), (0, 5), (0, 7),
//          (2, 5), (2, 7),
//          (3, 4), (3, 5), (3, 6), (3, 7),
//          (4, 5), (4, 7),
//          (6, 7)
//
// このようなペアの数は、BITを使うと効率的に求められる。
// たとえば B[5]に注目すると、これより左にある、B[5]より小さい値であるようなidxの数が、BITにより効率的に求められる

#include <iostream>
#include <sstream>
#include <string>
#include <cassert>
#include <cmath>
#include <climits>
#include <cstdio>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <deque>
#include <algorithm>
#include <functional>
#include <numeric>
#include <iomanip>
#include <cstring>
#include <fstream>

using namespace std;
typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;

#define REP(i,n) for(int i = 0; i < (int)(n); ++i)
#define FOR(i,a,b) for(int i = (a); i < (int)(b); ++i)
#define ALL(c) (c).begin(), (c).end()
#define SIZE(v) ((int)v.size())

#define pb push_back
#define mp make_pair
#define mt make_tuple

using namespace std;

class MajoritySubarray {
public:
    vector<ll> F;
    int fenwick_sum(int i) { // BITのこと。実装は蟻本と同じ
        int s = 0;
        while (i > 0) {
            s += F[i];
            i -= i & -i;  // i & -i: 1, 2, 4, 8, 16, ...のうち、iを割れる最大の数
                          //         言い換えると、iの一番右に立っているビットに等しい
        }
        return s;
    }
    void fenwick_add(int i, int val) {
        while(i < (1<<18)) {
            F[i] += val;
            i += i & -i;
        }
    }

    long long solve(ll N, ll seed, ll MOD) {
        vector<ll> As(N);
        REP(i, N) {
            As[i] = (seed / 65536) % MOD;
            seed = (seed * 1103515245 + 12345) % 2147483648;
        }

        F.resize(1<<18);
        ll ans = 0;
        REP(x, MOD) {
            fill(ALL(F), 0);
            int bal = N;
            REP(i, N) {
                fenwick_add(bal, 1);
                if (As[i] == x) {
                    ++bal;
                } else {
                    --bal;
                }
                ans += fenwick_sum(bal - 1);
            }
        }
        return ans;
    }

    long long getCount(int n_, int seed_, int m_) {
        ll n = n_;
        ll seed = seed_;
        ll MOD = m_;

        return solve(n, seed, MOD);
    }
};

TopCoderのプラグインGreedで自動生成したコードで「あいまいなシンボルです」エラーを回避

C++ programming contest

TopCoderのプラグインGreed（詳細はTopCoder の強欲プラグイン、Greed を使う！）で自動生成したコードをVisual Studioでビルドすると以下のようなエラーが出ます。

error C2872: 'data': あいまいなシンボルです。
note: 'std::ifstream data' である可能性があります。
note: または 'data'

これが指すコードは以下です。

using namespace std;

ifstream data("MajoritySubarray.sample");

string next_line() {
    string s;
    getline(data, s);
    return s;
}

エラーの原因は、<string>ヘッダで提供されるstd::dataという関数と名前がバッティングしていることでした。コード中のdataをdata2と変えるとビルドが通りました。

当たり前ですが、ちゃんとしたコードではそもそもusing namespace std;しないようにしましょう。

Linuxにおける2つのクリップボード - clipboard bufferとprimary selection

linux

Linuxには「2種類のクリップボードのようなもの」があるということを今更明確に理解しました。

一つ目の「クリップボードのようなもの」は、Windowsなどでおなじみの"clipboard buffer"です。たとえばChromeにて、文字列を選択してCtrl + Cでコピー、Ctrl + Vでペーストという例のやつです。端末(Terminal)の場合はCtrl + C, Ctrl + Vがほかの用途で予約されているので、かわりにCtrl + Shift + Cでコピー、Ctrl + Shift + Vでペーストになります。

二つ目の「クリップボードのようなもの」は、"primary selection"とよばれるものです。例えば、端末にてマウスで文字列を選択したあと、Chromeの検索窓にてホイールクリックをするだけで、選択した文字列がコピーされます。

"primary selection"の方を意識して使うようにするとわずかながら作業速度が上がるかもしれません。

参考

シェルスクリプトでの'--'の意味は、オプションの打ち止め

linux

シェルスクリプトでハイフンを2つ繋げた'--'の意味を最近はじめて知りました。これはオプションの打ち止めを意味します。

たとえばtouchコマンドで-hという空ファイルを作成しようとして

$ touch -h

とすると、-hがtouchコマンドのオプションとして解釈されてしまうため、下記のようにエラーが出てしまいます。

$ touch -h
touch: ファイルオペランドがありません
Try 'touch --help' for more information.

ここで--の出番です。--の後ろにある文字列はtouchのオプションではないことを明示できるので、以下のようにすることで-hという名前の空ファイルを生成できます。

$ touch -- -h