Pythonのコード改善のためのツール5つを試してみた

python

Pythonのコードを改善するためのツールについて一通り試してみました。各ツールのインストール方法や使い方については Pythonのスタイルガイドとそれを守るための各種Lint・解析ツール5種まとめ！ - Sider Blog に詳細にまとまっているのでおすすめです。

サンプルコード

以下のサンプルコードを対象に、各ツールの出力を確かめてみます。

import time
import sys
import fractions

def func1(varA,varB):
    '''return sum of a and b'''
    varC = 42
    return (varA + varB)


print(func1(fractions.Fraction(1, 2), fractions.Fraction(1, 3)))
3 + 5
sys.exit(0)

このスクリプトをsample.pyという名前で保存したのち、以下のスクリプトで実験をしました。各ツールはすべてデフォルトです。

for name in pycodestyle \
            pydocstyle \
            pyflakes \
            flake8 \
            pylint
do
    echo "================="
    echo ${name}
    echo "================="
    echo "version:" $(${name} --version)
    ${name} sample.py
    echo ""
done

各ツールの紹介

pycodestyle

  • 用途：スタイルチェック
    • PEP 8で規定されているPython公式のスタイルガイドに違反する箇所を見つけてくれる。
  • 詳細

出力は以下のとおりです。空白行やスペースの不足が指摘されています。

version: 2.0.0
sample.py:5:1: E302 expected 2 blank lines, found 1
sample.py:5:15: E231 missing whitespace after ','

pydocstyle

出力は以下のとおりです。docstringに関するエラーのみが指摘されています。

version: 2.1.1
sample.py:1 at module level:
        D100: Missing docstring in public module
sample.py:5 in public function `func1`:
        D300: Use """triple double quotes""" (found '''-quotes)
sample.py:5 in public function `func1`:
        D400: First line should end with a period (not 'b')
sample.py:5 in public function `func1`:
        D403: First word of the first line should be properly capitalized ('Return', not 'return')

pyflakes

  • 用途：エラー解析
  • 詳細
    • ソースファイルを解析して、importしたけど使っていないモジュールや、未使用の変数などを報告。
    • pycodestyleとは違い、コードのスタイルについては分析しない。
    • 高速さに重点を置く。Pylintよりも高速。
    • 正しいコードを過ってエラーと報告しないように細心の注意が払われている。
    • 解析はファイル単位で行われるので、importで読み込むファイルにまたいだ解析はできない。
    • もっと詳しい分析が必要な場合は、Flake8を使う。

出力は以下のとおりです。未使用のimport文、未使用の変数が報告されています。

version: 1.2.3
sample.py:1: 'time' imported but unused
sample.py:7: local variable 'varC' is assigned to but never used

Flake8

  • 用途：スタイルチェック + エラー解析 + 複雑度チェック
  • 詳細

出力は以下のとおりです。ちょうど出力がpycodestyleとpyflakesの和になっていることがわかります。

version: 2.6.2 (pycodestyle: 2.0.0, mccabe: 0.5.3, pyflakes: 1.2.3) CPython 3.6.4 on Windows
sample.py:1:1: F401 'time' imported but unused
sample.py:5:1: E302 expected 2 blank lines, found 1
sample.py:5:15: E231 missing whitespace after ','
sample.py:7:5: F841 local variable 'varC' is assigned to but never used

Pylint

  • 用途：スタイルチェック + エラー解析
  • 詳細
    • 老舗のツール。effective pythonで推されていた。
    • Plugin機能がある。

出力は以下のとおりです。flake8では報告されなかった、変数名の規則違反(varAなど)や、意味のない文('3 + 5')が報告されています。 最終行にこのコードの得点が表示されるのが面白いです。この場合10点満点中0点という意味です。さらには負の点数になることもあります。

No config file found, using default configuration
 Python 3.6.4 |Anaconda custom (64-bit)| (default, Jan 16 2018, 10:22:32) [MSC v.1900 64 bit (AMD64)]
No config file found, using default configuration
************* Module sample
C:  5, 0: Exactly one space required after comma
def func1(varA,varB):
              ^ (bad-whitespace)
C:  8, 0: Unnecessary parens after 'return' keyword (superfluous-parens)
C:  1, 0: Missing module docstring (missing-docstring)
C:  5, 0: Argument name "varA" doesn't conform to snake_case naming style (invalid-name)
C:  5, 0: Argument name "varB" doesn't conform to snake_case naming style (invalid-name)
C:  7, 4: Variable name "varC" doesn't conform to snake_case naming style (invalid-name)
W:  7, 4: Unused variable 'varC' (unused-variable)
W: 12, 0: Statement seems to have no effect (pointless-statement)
W:  1, 0: Unused import time (unused-import)

------------------------------------------------------------------
Your code has been rated at 0.00/10 (previous run: 0.00/10, +0.00)

flake8のプラグイン

flake8のプラグインを別途インストールすることで、flake8にルールを追加することが可能です。pip search flake8-するとわかるのですが、大量のプラグインが存在しています。ここでは2つだけ紹介します。

hacking

  • 用途：Flake8にOpenstack社のルールを追加するプラグイン
  • 詳細
    • pip install hackingしたあとにflake8を呼び出すと、flake8に、Openstack社が独自に追加したルールが追加される

出力は以下のとおりです。Hで始まる行が、flake8に追加された検出された項目になります。この例だとimport文の順番がアルファベット順になっていなかったことで怒られています。

version: 2.6.2 (pycodestyle: 2.0.0, mccabe: 0.5.3, ProxyChecker: 0.0.1, hacking.core: 0.0.1, pyflakes: 1.2.3) CPython 3.6.4 on Windows
sample.py:1:1: F401 'time' imported but unused
sample.py:2:1: H306  imports not in alphabetical order (time, sys)
sample.py:3:1: H306  imports not in alphabetical order (sys, fractions)
sample.py:5:1: E302 expected 2 blank lines, found 1
sample.py:5:15: E231 missing whitespace after ','
sample.py:7:5: F841 local variable 'varC' is assigned to but never used

flake8-docstrings

  • 用途：Flake8にpydocstyleのルールを追加するプラグイン
  • 詳細
    • pip install pydocstyleしたあとにflake8を呼び出すと、flake8に、pydocstyleのルールが追加される

出力は以下のとおりです。Dで始まる行がdocstringに関する報告です。

version: 2.6.2 (pycodestyle: 2.0.0, mccabe: 0.5.3, pyflakes: 1.2.3, flake8-docstrings: 1.3.0, pydocstyle: 2.1.1) CPython 3.6.4 on Windows
sample.py:1:1: D100 Missing docstring in public module
sample.py:1:1: F401 'time' imported but unused
sample.py:5:1: D300 Use """triple double quotes"""
sample.py:5:1: D400 First line should end with a period
sample.py:5:1: D403 First word of the first line should be properly capitalized
sample.py:5:1: E302 expected 2 blank lines, found 1
sample.py:5:15: E231 missing whitespace after ','
sample.py:7:5: F841 local variable 'varC' is assigned to but never used

サンプルコードの修正例

すべてのエラーや警告が消えるまで修正したのが以下のコードです。良いコードになったでしょうか？

"""Sample for hatena diary."""

import fractions
import sys


def func1(var_a, var_b):
    """Return sum of a and b."""
    return var_a + var_b


print(func1(fractions.Fraction(1, 2), fractions.Fraction(1, 3)))
sys.exit(0)

その他のツール

mypy

型ヒントのためのツールです。自分はまだ理解できていません。

  • 用途： Optional Static Typing for Python
    • Python 3.6に導入された型ヒントを使っている人向けのツール。型の誤りを静的に解析する。

PyChecker

  • 用途：スタイルチェック + エラー解析
  • 詳細
    • Python 2.x系が対象。開発が止まっているようなので今から使う理由はないでしょう。

個人的な結論

  • まずはflake8をデフォルトで使う
  • flake8で足りないところがあれば、Pylintを併用するか、flake8の適当なプラグインを見つけて導入する

flake8についてはこちらの記事により詳しい使い方をまとめました。

参考

KUPC 2018 C - 七目

programming contest

本番中にパニクって解けなかったのを反省して解き直しました。DFSくらいさっと書けないとだめですね…。

  • 問題

    • C - 七目
    • 9x9のマスが白で埋められている。白マスが縦・横・斜めのいずれにも白マスが7個以上連続しないように、11個のマスを黒く塗る塗り方を求める

  • 解法

    • 素直にDFS。黒く塗るしかないマスは黒く塗り、そうでない箇所では黒く塗る場合・塗らない場合の両方を試す。
    • 各行に少なくとも一つは黒マスが必要なことを利用して枝刈りすると、自分の古いPCで1分20秒ほどですべての解が出ました。

  • コード

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


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

const int N = 9;
const int W = N;
const int H = N;

void debug_print(vector<vector<char>>& cells, int h_now, int w_now, bool print_cursor=true) {
    cout << "=================" << endl;
    REP(h, H) {
        REP(w, W) {
            if (h == h_now && w == w_now && print_cursor) {
                cout << "[";
            }
            else {
                cout << " ";
            }
            cout << (cells[h][w] == 1 ? '#' : '.');
            if (h == h_now && w == w_now && print_cursor) {
                cout << "]";
            }
            else {
                cout << " ";
            }
        }
        cout << endl;
    }
    cout << "=================" << endl;
}


pair<int,int> idx_to_h_w(int idx) {
    return {idx / N, idx % N};
}

// 今いるマスの左・上・右上・左上を見て、このマスを黒く塗る必要があるかを判定
bool should_put_stone(int idx, vector<vector<char>>& cells) {
    auto hw = idx_to_h_w(idx);
    int h, w;
    tie(h, w) = hw;

    int dx[] = {-1, 0, -1, 1};
    int dy[] = {0, -1, -1, -1};
    REP(d, 4) {
        int seq = 0;
        FOR(i, 1, 7) {
            int h2 = h + dy[d] * i;
            int w2 = w + dx[d] * i;
            if (w2 < 0 || w2 >= W || h2 < 0 || h2 >= H) {
                break;
            }
            if (cells[h2][w2] == 0) ++seq;
        }
        if (seq == 6) return true;
    }
    return false;
}

bool dfs(int idx, vector<vector<char>>& cells, int remain) {
    auto hw = idx_to_h_w(idx);
    int h, w;
    tie(h, w) = hw;

    // 今いるマスをとりあえず白くする
    // これを忘れると正しくdfsできない！
    cells[h][w] = 0;

//    cout << "remain = << " << remain << endl;
//    debug_print(cells, h, w);

    // 各行に少なくとも1つ黒マスが必要なことから
    // もう条件を満たせないならfalseを返す
    int minimum_need = H - h - 1;
    if (remain < minimum_need) return false;

    // 最後の判定
    if (idx == H * W - 1) {
        if (should_put_stone(idx, cells)) {
            if (remain > 0) {
                cells[h][w] = 1;
                remain--;
            }
            else {
                return false;
            }
        }
        // 条件を満たすのでprint
        debug_print(cells, h, w, false);
        return true;
    }

    // 上・左・右上・左上いずれかに6連続で白マスなら現マスは黒く塗るしかない
    if (should_put_stone(idx, cells)) {
        // もう残数がなければおわり
        if (remain == 0) {
            return false;
        }
        else {
            cells[h][w] = 1;
            return dfs(idx + 1, cells, remain - 1);
        }
    }
    // 石を置いても置かなくてもよい
    else {
        // 石を置かない場合
        auto ret = dfs(idx + 1, cells, remain);
        // 石を置く場合
        if (remain > 0) {
            cells[h][w] = 1;
            ret |= dfs(idx + 1, cells, remain - 1);
        }
        return ret;
    }
}

int main(void)
{
    vector<vector<char>> cells(N, vector<char>(N));
    dfs(0, cells, 11);
    return 0;
}

出力されたのは以下の16盤面です。

=================
......#..
....#....
..#......
#......#.
.....#...
...#.....
.#......#
......#..
....#....
=================
=================
......#..
...#.....
#......#.
....#....
.#......#
.....#...
..#......
......#..
...#.....
=================
=================
......#..
..#......
.....#...
.#......#
....#....
#......#.
...#.....
......#..
..#......
=================
=================
......#..
..#......
#......#.
.....#...
....#....
...#.....
.#......#
......#..
..#......
=================
=================
......#..
..#......
#......#.
...#.....
....#....
.....#...
.#......#
......#..
..#......
=================
=================
.....#...
...#.....
.#......#
......#..
....#....
..#......
#......#.
.....#...
...#.....
=================
=================
.....#...
..#......
......#..
...#.....
#......#.
....#....
.#......#
.....#...
..#......
=================
=================
....#....
......#..
.#......#
...#.....
.....#...
#......#.
..#......
....#....
......#..
=================
=================
....#....
..#......
#......#.
.....#...
...#.....
.#......#
......#..
....#....
..#......
=================
=================
...#.....
......#..
..#......
.....#...
.#......#
....#....
#......#.
...#.....
......#..
=================
=================
...#.....
.....#...
#......#.
..#......
....#....
......#..
.#......#
...#.....
.....#...
=================
=================
..#......
......#..
...#.....
#......#.
....#....
.#......#
.....#...
..#......
......#..
=================
=================
..#......
......#..
.#......#
.....#...
....#....
...#.....
#......#.
..#......
......#..
=================
=================
..#......
......#..
.#......#
...#.....
....#....
.....#...
#......#.
..#......
......#..
=================
=================
..#......
.....#...
.#......#
....#....
#......#.
...#.....
......#..
..#......
.....#...
=================
=================
..#......
....#....
......#..
.#......#
...#.....
.....#...
#......#.
..#......
....#....
=================

TCO19 Single Round Match 737 - Div1 Easy (AliceAndBobEasy)

programming contest

2018-09-19のTopcoderに３ヶ月ぶりで出場。実験を執念で繰り返してかろうじてEasyが解けました。あわててNimゲームの勝利条件を蟻本で確認しましたが、もはやこれが前提知識となっているとは。

以下、Easyの解法とPythonコードです。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# [問題]
#   Nが与えられる。
#   以下の条件を満たす、N個の山を構成せよ。
#     Nimゲームで先攻が勝つような先攻のmove (winning move) 数を最大化する
#     各山の石の個数は1 ~ 737373737個
#   Nは1~37。
#
# [解法]
#   まずNimの勝利条件を蟻本で確認。
#   これによると、N個の山の石の数でxorをとった結果が0なら負け、非0なら勝ちの状態である。
#
#   まずNが奇数の場合を考える。
#   例えば山が5個あるとする。
#   山を順にa, b, c, d, e,
#   各山の石の数をn(a), n(b), ...とする。
#
#   １つ目の山から任意の個数の石を取り除いたあと、
#   勝ちの状態を作れるか、言い換えると、残りの山のXORをちょうど0にできるか、を考えてみる。
#
#   石の個数をa, b, c, d, eだとすると、
#   山aから石を取り除いた後の値が、n(b) ^ n(c) ^ n(d) ^ n(e)に等しくなることが、
#   この山から石を取る場合の唯一の勝利条件である。
#
#   必要な条件は、n(a) > n(b) ^ n(c) ^ n(d) ^ n(e)。
#
#   ここで、n(a) = 0b100001
#          n(b) = 0b100010
#          n(c) = 0b100011
#          n(d) = 0b100100
#          n(e) = 0b100101
#   としてみると、n(b) ^ n(c) ^ n(d) ^ n(e)の結果、最上位ビットは0になる。
#   よって n(a) > (n(b) ^ n(c) ^ n(d) ^ n(e)) とできるので、
#   n(a)から石をとることで相手を負けにできる。
#   Nが奇数のときはこの構成法でN種類のwinning moveができる。
#
#   Nが偶数のときは、実はN-1種類のwinning moveが常に最大値となる。
#   N = 4の場合を考えると、
#          n(a) = 0b100001
#          n(b) = 0b100010
#          n(c) = 0b100011
#          n(d) = 0b000001
#   と、最後の山を1にすることで、N-1種類のwinning moveができる。
#
#   では、N種類のwinning moveは可能なのか。
#   本番中では、実験的に不可能でありそうなことを確かめたが、
#   証明はできなかった。https://www.topcoder.com/blog/single-round-match-737-editorials/ 
#   を読んで、証明を理解できた。
#   
#   証明:
#   今、Nが偶数で、winning positionにいるとする。(つまり先手が石を取る直前である）
#   すべての石のxorは非ゼロなので、
#   すべての石のxorをとった値のうち、もっとも左側の1が立つビットの位置をxとする。
#   山の数は偶数だが、x番目のビットが立っている山の数は奇数である。
#   ということは、少なくとも1つの山(Pとおく)は、そのx番目のビットが立っていない。
#   山Pについては、どう石をとっても、先手は勝てない。
#   よってNが偶数の場合は、最大でもN-1種類のwinning moveしかできない。
#
# [本番]
#   ブサイクだが上記解法とほぼ同じ実装ができた。
#   生成された配列が条件を満たすかチェックする関数を作っていたので2個撃墜できた


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

class AliceAndBobEasy:
    def make(self, N):
        if N == 1:
            return (737,)
        elif N == 2:
            return (737, 373)
        elif N % 2 == 1:
            ret = []
            ret.append(1 << 28)
            ret.append((1 << 28) + (1 << 27))
            ret.append((1 << 28) + (1 << 26))
            for i in range(N - 3):
                ret.append((1 << 28) + i + 1)
            return tuple(ret)
        else:
            ret = []
            ret.append(1)
            ret2 = list(self.make(N - 1))
            ret += ret2
            return tuple(ret)

AGC027 B - Garbage Collector

programming contest

2018-09-15のAGCに出場。Aと、Bの部分点が解けて325位。Bは考え方は正しかったのですが、計算途中のオーバーフローに泣きました。上位陣でもBはスルーしている人が多かった中で、DPをやりたくなる衝動を抑えつつ部分点が取れたので、少し自信が付きました。もっとも証明できてない箇所があるのですが。

以下、Bの解法とコードです。g++の拡張(__int128)を使っているため、Visual Studioではビルドが通りません。

// 本番時、解説PDFと同様に考えられたが、オーバーフローでlargeがWA。
// 証明できていない箇所がある。
//
// [解法]
// 以下の場合を例にとって考える
//    |------*-----*-----*------*-----*------*-----*
//    0      a     b     c      d     e      f     g
// 原点から距離iの場所にあるゴミを、ゴミiと呼ぶことにする。
//
// ロボットがゴミ箱にゴミを捨てる回数をK回と固定して考える。例えばK=3とする。
// 
// このとき、
// 1回目のゴミ捨ては、「ゴミgまで行ってゴミgを拾い、ゴミa ~ dのどれかを拾いながら原点に戻る」
// 2回目のゴミ捨ては、「ゴミfまで行ってゴミfを拾い、ゴミa ~ dのどれかを拾いながら原点に戻る」
// 3回目のゴミ捨ては、「ゴミeまで行ってゴミeを拾い、ゴミa ~ dのどれかを拾いながら原点に戻る」
// のが最適となる。本番中も直感的にそう思ったが、証明できていない。解説pdfにも「明らか」としか書かれていない。
//
// 上記を信じると、ゴミを拾いに行く行きのコストは常に定数となる。
// したがって残りの問題は、「3回のゴミ捨ての帰りのうち、ゴミa, b, c, dをどのタイミングで拾うか」
// を、帰りのコストが最小になるように決める問題となる。
//
// まずゴミa~cはないものとして、ゴミdについてついて考える。するとゴミdは3回のうち
// どのタイミングで拾っても、かかるコストは同じになることがわかる。よって1回目で拾うことにする。
//
// 次、ゴミcについて考える。選択肢として、1回目で拾うか、2回目で拾うかの2通りがある。
// (3回目に拾うのは、2回目に拾うのとコストが同じであるので考える必要がない)
// これは2回目に拾うのがよい。
// なぜなら、もし1回目で拾うと、ゴミbからゴミcまでの区間のコストが9になるのに対し、
// 2回目で拾うと、同区間のコストが4になるから。
// （ここも証明が怪しい）
//
// 同様に、遠くにあるゴミから順に、1回目のゴミ捨て、2回目のゴミ捨て、3回目のゴミ捨て、
// また最初に戻って1回目のゴミ捨て、...と拾うのが最適。
//
// この場合だと、
// 1回目のゴミ捨ては、「ゴミgまで行ってゴミgを拾い、ゴミd、ゴミaを拾いながら原点に戻る」
// 2回目のゴミ捨ては、「ゴミfまで行ってゴミgを拾い、ゴミcを拾いながら原点に戻る」
// 3回目のゴミ捨ては、「ゴミeまで行ってゴミgを拾い、ゴミbを拾いながら原点に戻る」
// が最適となる。
//
// 累積和を使うとこれを効率的に計算でき、largeが解ける。

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

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

int main(void)
{
    cin.sync_with_stdio(false);
    ll N, X;
    cin >> N >> X;
    vector<ll> Xs(N);
    REP(n, N) cin >> Xs[n];

    // 原点から遠い順に距離をソート
    auto Ys = Xs;
    reverse(ALL(Ys));

    // 累積和を計算
    vector<ll> Y_sum(N);
    Y_sum[0] = Ys[0];
    FOR(i, 1, N) Y_sum[i] = Y_sum[i - 1] + Ys[i];

    __int128 ans = 1e18;  // ★ 最終的な解は1e18以下に収まるのだが、他の解候補がlong long intの範囲を超えることが
                          // あるため、__int128を使う必要がある！！！
    FOR(k, 1, N + 1) { // ゴミ箱に行く回数をkとする

        // ★この変数がlong long intを超える!!
        // ゴミを捨てる回数に応じたコストをまず計算
        __int128 local_ans = X * k;

        // 距離dにあるゴミを取りに行くコストはd, 取ったあと帰るコストは2^2 * d = 4d。
        // よって足して5d。
        // k回のゴミ捨てで、原点から遠いゴミk個をそれぞれ取って帰るコストは、
        // 累積和を使って以下のように計算できる。
        local_ans += 5 * Y_sum[k-1];

        ll i = 2 * k - 1;
        ll n = 1;

        // まだ取っていないゴミのうち、遠くにあるゴミから順番にk個のゴミを、
        // k回のゴミ捨てでの帰路にそれぞれ1個ずつ取っていくところをシミュレート。
        // もともとn個のゴミを持って帰るところをn+1個のゴミを持って帰ることにすることで
        // 増えるコストは、((n+2)^2 - (n+1)^2 = 2n + 3)。
        while (i <= N - 1) {
            local_ans += (2LL * n + 3) * (Y_sum[i] - Y_sum[i - k]);
            n += 1;
            i += k;
        }
        // 端数処理
        local_ans += (2LL * n + 3) * (Y_sum[N - 1] - Y_sum[i - k]);

        ans = min(ans, local_ans);
    }

    // __int128は出力に失敗するのでllにキャストする
    ll ans2 = (ll)ans;
    // ゴミを拾うコストは常に定数なので最後に足す
    cout << (ans2 + N * X) << endl;

    return 0;
}

ARC102 D - All Your Paths are Different Lengths

programming contest

全然分かりませんでしたが、解説を見て理解できました。思いつけば簡単な問題なので悔しいです。

問題：D: All Your Paths are Different Lengths - AtCoder Regular Contest 102 | AtCoder

問題概要：パスの長さが0..L-1のちょうどL通りになるグラフを生成

解法： 例としてL=37の場合を考えてみます。

まず、長さが0から31までの32通りのパスからなるグラフは、以下のように構築できることに気付く必要があります。 f:id:minus9d:20180902003543p:plain

L=37ですから、上記のグラフを改造して、長さが32, 33, 34, 35, 36のパスをさらに追加する必要があります。

ここで、頂点1から頂点3までは、ちょうど4種類の長さ0, 1, 2, 3が取れることに着目します。頂点3から頂点6までコスト32の辺をショートカットとして加えると、 ちょうど長さ32, 33, 34, 35のパスが追加できることがわかります。 f:id:minus9d:20180902004155p:plain

残るは長さ36のパスです。今度は頂点1から頂点6までコスト36の辺を加えればよいです。よって以下が最終形です。 f:id:minus9d:20180902004345p:plain

終了後ACしたコードを以下に記します。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# 終了後に解説を見て解けた
# パス数が２のべき乗であれば簡単にグラフを構築できることに気付きたかった。
#
# 本番中は、10進法で考えて 頂点1と頂点2の間を辺{0, 100, 200, ..., 900}で繋ぎ、頂点2と頂点3の間を辺{0, 10, ..., 90}で繋ぎ、
# などという事を考えていたが、それ以上進展させることが難しかった。


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

L = int(input())

v = 1
edges = []  # (from, to, cost)

# L = 37の場合を例にとりコメントを付与
# まず長さ0 ~ 31までのパスを作成する。
#   頂点1 -> 頂点2を0の辺と1の辺で結ぶ
#   頂点2 -> 頂点3を0の辺と2の辺で結ぶ
#   頂点3 -> 頂点4を0の辺と4の辺で結ぶ
#   頂点4 -> 頂点5を0の辺と8の辺で結ぶ
#   頂点5 -> 頂点6を0の辺と16の辺で結ぶ
n = 2
cost = 1
while n <= L:
    edges.append((v, v+1, 0))
    edges.append((v, v+1, cost))
    cost *= 2
    v += 1
    n *= 2

# 頂点6をv_sink(最終到達点)とする
v_sink = v

# ここから先、長さ32, 33, 34, 35, 36のパスを生成することになる。
# 頂点3と頂点6を、32の辺で直でつなぐと、長さ32, 33, 34, 35のパスが作れる。
# 頂点1と頂点6を、36の辺で直でつなぐと、長さ36のパスが作れる。
n //= 2
N = n
L -= N

while L > 0:
    if L >= n:
        edges.append((v, v_sink, N))
        L -= n
        N += n
    n //= 2
    v -= 1

print(v_sink, len(edges))
for e in edges:
    print(*e)

SoundHound Programming Contest 2018 Masters Tournament 本戦 B - Neutralize

programming contest

問題：B - Neutralize

問題概要：N個の薬品が並んでいる。各薬品の効用は-109 ~ 109。連続したK個の薬品の効用を0にする処理を任意回おこなったあとの薬品の効用の和の最大値を求める。

400点問題なのに2時間考えてまったくわかりませんでした。実行時間制限がなかったとしても解けてなかったと思います。sugim48さんの回答 ほぼそのまま参考にして本番後ACしたコードを最後に示します。また具体例をもとにアルゴリズムの動きを書き下した説明を以下に示します。

例としてK = 3で薬品が3 4 -8 2 -7 -5と並んだ場合を考えます。この解法では左側の薬品から順に部分配列を伸ばしていき、最後に増やした薬品を右端としたK個の薬品の効用を0にすべきかどうかを決めます。以下でその手順を見ていきましょう。

まず1個目を取り出して3。これではK個の薬品を取り出しようがないのでそのまま。この部分配列に対する答えは3。これをdp[1] = 3として記録します。

次に2個目を取り出して3 4。同様の理由でそのまま。この部分配列に対する答は7。これをdp[2] = 7として記録します。

次に3個目を取り出して3 4 -8。この部分配列に対する答は、「3 4 -8の3つの効用をゼロにするか、しないかのうち、良い方」です。今回はゼロにしたほうが得です。よってこの部分配列に対する答は0。これをdp[3] = 0として記録します。

次に4個目を取り出して3 4 -8 2。この部分配列に対する答は、「4 -8 2の3つの効用をゼロにするか、しないかのうち、良い方」です。4 -8 2の効用をゼロにしない場合は簡単で、dp[2]に2を足した値なので、2になります。4 -8 2の効用をゼロにする場合はちょっとややこしくて、4 -8 2だけでなくて4 -8 2のすぐ左側にある任意の数の薬品をもゼロにしてよい、と考えなければいけません。それを求めるにはこれまでのdpの値の最大値を記録する配列maを用意しておきます（詳細はコードを参照）。この場合は4 -8 2を0とするのがベストなので、3になります。2と3のうち良い方である3が答なので、dp[4] = 3として記録します。

次に5個目を取り出して3 4 -8 2 -7。この部分配列に対する答は、「-8 2 -7の3つの効用をゼロにするか、しないかのうち、良い方」です。上と同様の手順により、-8 2 -7をゼロにするのがベストなので、dp[5] = 7として記録します。

最後に6個目を取り出して3 4 -8 2 -7 -5。この部分配列に対する答は、「2 -7 -5の3つの効用をゼロにするか、しないかのうち、良い方」です。今回の場合は2 -7 -5の3つの効用をゼロにし、さらにその左側にある-8も巻き込んでゼロにするのがベストです。よってdp[6] = 7として記録します。

以上のアルゴリズムによるdp, maの値の変化を以下に示します。Bの値とdp, maの値との比較をしやすくするため、あえてdp, maの値は本来の位置よりひとつ左にずらしています。

idx 0 1 2 3 4 5
B 3 4 -8 2 -7 -5
dp 3 7 0 3 7 7
ma 3 7 7 7 7 7

しつこく考えてなんとか理解できたつもりですが、本当に正しいか、また人が読んでわかる説明になっているかは自信ありません。

続きを読む

CUDAを使って2枚の画像の平均をとるプログラム

cuda C++

CUDAの勉強がてら書いてみました。2枚の同サイズの画像を読み込んで、平均画像を作成して保存するだけのプログラムです。環境はUbuntu 16.04です。

// 2枚の画像の平均を取るプログラム
// 画像は以下から6000x4000のものを取得
//     https://www.pexels.com/photo/black-plants-photo-943903/
//     https://www.pexels.com/photo/person-riding-on-zip-line-1090551/

#include <chrono>
#include <cstdio>
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>

const int THREADS_PER_BLOCK = 16;

// https://developer.nvidia.com/cuda-example のサンプルから引用
static void HandleError( cudaError_t err, const char *file, int line ) {
    if (err != cudaSuccess) {
        printf( "%s in %s at line %d\n", cudaGetErrorString( err ), file, line );
        exit( EXIT_FAILURE );
    }
}
#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))

__global__ void use_gpu_sub(unsigned char *a, unsigned char *b, unsigned char *c, int N) {
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    if (tid < N) {
        c[tid] = (a[tid] + b[tid]) / 2;
    }
}

void use_gpu(unsigned char* a, unsigned char* b, unsigned char* c, int height, int width, int channel) {
    unsigned char *dev_a, *dev_b, *dev_c;
    const int N = height * width * 3;

    // GPUのためのメモリ確保
    HANDLE_ERROR( cudaMalloc( (void**)&dev_a, N * sizeof(unsigned char) ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_b, N * sizeof(unsigned char) ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_c, N * sizeof(unsigned char) ) );

    // a, bがさす画像データをGPUのメモリにコピー
    HANDLE_ERROR( cudaMemcpy( dev_a, a, N * sizeof(unsigned char),
                            cudaMemcpyHostToDevice ) );
    HANDLE_ERROR( cudaMemcpy( dev_b, b, N * sizeof(unsigned char),
                            cudaMemcpyHostToDevice ) );

    // (N + THREADS_PER_BLOCK - 1)/THREADS_PER_BLOCK個のブロックをGPU上で実行。
    // 各ブロックでは THREADS_PER_BLOCK 個のスレッドが立つ
    use_gpu_sub<<<(N + THREADS_PER_BLOCK - 1)/THREADS_PER_BLOCK, THREADS_PER_BLOCK>>>( dev_a, dev_b, dev_c, N );

    // GPのメモリにある計算結果をCPU側にコピー
    HANDLE_ERROR( cudaMemcpy( c, dev_c, N * sizeof(unsigned char),
                            cudaMemcpyDeviceToHost ) );

    // free the memory allocated on the GPU
    HANDLE_ERROR( cudaFree( dev_a ) );
    HANDLE_ERROR( cudaFree( dev_b ) );
    HANDLE_ERROR( cudaFree( dev_c ) );
}

void use_cpu(unsigned char* a, unsigned char* b, unsigned char* c, int height, int width, int channel) {
    int N = height * width * channel;
    for (int i = 0; i < N; ++i) {
        c[i] = (a[i] + b[i]) / 2;
    }
}

void average_two_images() {
    cv::Mat img1 = cv::imread("adventure-asia-bali-1090551.jpg");
    cv::Mat img2 = cv::imread("beautiful-close-up-color-943903.jpg");

    assert(img1.size() == img2.size());
    assert(img1.type() == img2.type());
    std::cout << "image size: " << img1.size() << std::endl;

    cv::Mat img3_gpu = cv::Mat(img1.size(), img1.type());
    cv::Mat img3_cpu = cv::Mat(img1.size(), img1.type());

    // gpuを使って
    auto s1 = std::chrono::system_clock::now();
    use_gpu(img1.data, img2.data, img3_gpu.data, img1.rows, img1.cols, img1.channels());
    auto e1 = std::chrono::system_clock::now();
    std::cout << "gpu time: " << std::chrono::duration_cast<std::chrono::milliseconds>(e1 - s1).count() << " ms" << std::endl;

    // cpuを使って
    auto s2 = std::chrono::system_clock::now();
    use_cpu(img1.data, img2.data, img3_cpu.data, img1.rows, img1.cols, img1.channels());    
    auto e2 = std::chrono::system_clock::now();
    std::cout << "cpu time: " << std::chrono::duration_cast<std::chrono::milliseconds>(e2 - s2).count() << " ms" << std::endl;

    imwrite("averaged.jpg", img3_gpu);
    
    // 結果が同じことを確認
    assert(cv::countNonZero(img3_gpu != img3_cpu) == 0);
}

int main()
{
    average_two_images();
    return 0;
}

これを以下のコマンドでビルドします。

nvcc -O2 -std=c++11 average_two_images.cu $(pkg-config --cflags --libs opencv)

以下は実行結果です。GPUは転送速度がネックになっていてCPUより遅いです。

image size: [6000 x 4000]
gpu time: 180 ms
cpu time: 51 ms

コード中で定数にしている THREADS_PER_BLOCKの値を変えると処理時間が変わりました。どうやって最適な値を決めればよいのでしょう？

    176ms: 512スレッド
    172ms: 256スレッド
    158ms: 128スレッド
    158ms: 64スレッド
    180ms: 32スレッド
    180ms: 16スレッド