pydub.split_on_silence を使ってm4a,mp3,wavの無音部分を削除する

音声ファイルから無音部分をトリミングして保存するプログラムを書きました。

実装方法

pydubというライブラリを使います。 silencesplit_on_silenceという( str.split() のような)無音部分を削除する関数が用意されています。

ソースコード

recording-file.m4aというファイルから無音部分を削除してsilence_removed_-20.mp3などの名前で保存するプログラムです。

from pydub import AudioSegment
from pydub.silence import split_on_silence

SOURCE_FILE =  'recording-file.m4a'

sound = AudioSegment.from_file(SOURCE_FILE)

org_ms = len(sound)
print('original: {:.2f} [min]'.format(org_ms/60/1000))

for silence_thresh in range(-20, -60, -5):
    chunks = split_on_silence(sound, min_silence_len=100, silence_thresh=silence_thresh, keep_silence=100)
    if not chunks:
        continue
    cutted_sound = sum(chunks)
    cutted_ms = len(cutted_sound)
    print('silence_thresh = {}: {:.2f} [min]'.format(silence_thresh, cutted_ms/60/1000))
    cutted_sound.export('silence_removed_{}.mp3'.format(silence_thresh), format='mp3')

もともとがGoogle Pixelのレコードで録音したデータで試してみました。 無音のしきい値によって残る部分の時間は下記のとおりでした。耳で聞いた所しきい値は-35くらいがちょうど良いと感じました。

original: 35.37 [min]
silence_thresh = -20: 2.10 [min]
silence_thresh = -25: 5.58 [min]
silence_thresh = -30: 8.95 [min]
silence_thresh = -35: 11.40 [min]
silence_thresh = -40: 13.16 [min]
silence_thresh = -45: 15.91 [min]
silence_thresh = -50: 20.34 [min]
silence_thresh = -55: 25.17 [min]

real    5m56.119s
user    5m44.925s
sys     0m13.276s

平方根の値を計算する方法(二分法・ニュートン法・開平法)について

 \sqrt{5} = 2.2360679 \dots であることが知られていて「富士山麓オウム鳴く」のごろ合わせも有名です。実際に 2.2360679^{2} = 499999965341041 \lt 5 \lt 500000010062400 = 2.2360680^{2} であることから小数第6,7位が80ではなく79が正しいことが確認できます。

ところで平方根の値を計算する効率の良い方法を知っていますか?

この記事では平方根の値を求める方法として、計算方法を3つ紹介し、PythonとRustでの実装を示します。

  1. 二分法
  2. ニュートン法
  3. 開平法

二分法

求めたい平方根の値は  0 \le \sqrt{x} \le max(1.0, x) の範囲にあるため、二分法で範囲を絞っていくことで近似値を求めることができます。

1回のループごとに解の存在範囲が 1/2 倍になります。浮動小数点数型を使う場合は型の精度に気を付けましょう。

def sqrt_by_bisection(k, epslilon=0.001):
    """k-epsilon <= x*x < k であるxを返す。"""
    upper = 1 if k < 1 else k
    bottom = 0
    while upper - bottom >= epslilon:
        m = bottom + (upper - bottom) / 2
        if m**2 > k:
            upper = m-epslilon
        else:
            bottom = m
    return bottom
// k-epsilon <= x*x < k であるxを返す。
fn sqrt_by_bisection(k: f64, epsilon: f64) -> f64 {
    debug_assert!(k >= 0.0);
    let mut upper = if k < 1.0 { 1.0 } else { k };
    let mut bottom = 0.0;
    while upper - bottom >= epsilon {
        let m = bottom + (upper - bottom) / 2.0;
        if m * m > k {
            upper = m - epsilon;
        } else {
            bottom = m;
        }
    }
    bottom
}

ニュートン法

 \sqrt{x} の値を求めるという問題は、 f(x) = x^{2} - k として f(x) = 0 となる  x を求める問題と言い換えることができ、この問題はニュートン法で解の近似値を求めることができます。

この方法は二分法よりも収束が速いことが知られています。

2次収束と呼ばれ,実際に計算してみるとわかるように,非常に速い. 大ざっぱにいって「誤差が毎回2乗で減る」あるいは「正しい桁数が反復1回ごとに倍に増える」といえる.

『数値計算の常識』P60より

def sqrt_by_newton(k, epsilon=0.001):
    def f(x): return x - (x * x - k) / (2.0 * x)
    if k == 0.0:
        return 0.0
    x = k
    next_x = f(x)
    while abs(x - next_x) >= epsilon:
        x = next_x
        next_x = f(x)
    return next_x
fn sqrt_by_newton(k: f64, epsilon: f64) -> f64 {
    let f = |x| x - (x * x - k) / (2.0 * x);
    debug_assert!(k >= 0.0);
    if k == 0.0 {
        return 0.0;
    }

    let mut x = k;
    let mut next_x = f(x);
    while (x - next_x).abs() >= epsilon {
        x = next_x;
        next_x = f(x);
    }

    next_x
}

実測:ニュートン法は二分法よりも収束が速いのか

ニュートン法のほうが速いし実装も簡単です。

上記のPython実装を使って精度を変更しながらループ回数との関係を調べました。浮動小数点数型の誤差の問題でi=14以降計算することが出来ませんでした。Decimalモジュールを使うとより多い桁数を正しく求めることができるかも知れません。

f:id:matsu7874:20210109035543p:plain

開平法

手計算(筆算)で行いやすい方法として開平法という方法があります。

 10a + b \le \sqrt{x} であるような  b (10a + b)^{2} = 100a^{2} + (2 \times 10a + b)bという形に変形して求めていく方法です。詳しい方法は参考のリンクを参照してください。

実装上、整数で演算するので多倍長整数が使える言語(Pythonなど)ではかんたんに多くの桁数を正確に計算することが出来ます。

def sqrt_by_extraction(n, prec=16):
    chunks = []
    k = n
    while k > 0:
        chunks.append(k % 100)
        k //= 100
    sqrt = []
    a = 0
    b = 0
    for _ in range(prec):
        if chunks:
            c = chunks.pop()
        else:
            c = 0
        a *= 10
        b = b * 100 + c
        for i in range(10)[::-1]:
            t = (a + i) * i
            if t <= b:
                b -= t
                a += 2 * i
                sqrt.append(i)
                break
    return sqrt
fn sqrt_by_extraction(n: i64, prec: usize) -> Vec<i64> {
    debug_assert!(n >= 0);
    debug_assert!(prec <= 16);
    let mut chunks = vec![];
    let mut k = n;
    while k > 0 {
        chunks.push(k % 100);
        k /= 100;
    }
    eprintln!("{:?}", chunks);
    let mut sqrt = vec![];
    let mut a = 0;
    let mut b = 0;
    for _ in 0..prec {
        let c = chunks.pop().unwrap_or(0);
        a *= 10;
        b = b * 100 + c;
        for i in (0..10).rev() {
            let t = (a + i) * i;
            if t <= b {
                b -= t;
                a += 2 * i;
                sqrt.push(i as i64);
                break;
            }
        }
    }

    sqrt
}

まとめ

平方根の値を計算する3種類の方法(二分法・ニュートン法・開平法)について紹介し、それぞれPythonとRustでの実装例を示しました。 大きな桁数が必要であれば開平法を、そうでなければ高速なニュートン法を使うと良さそうです。

SRM788 Div2

2020/07/24にTopcoderで行われたSRM788の参加記です。

公式解説はTopcoder Single Round Match 788 Editorials | Topcoderです。

Easy NextOlympics

"YYYY.MM.DD" 形式で今日の日付が与えられるので東京オリンピックの開催予定日である"July 23rd, 2021"までの日数を出力する問題です。

実装するだけですが、datetime ライブラリを使うより独自実装したほうが早いと思ったため、datetime ライブラリを使わずに実装しました。 ある月のX日からその先月のX日までの日数は先月が何月かに依存して決まるのですね。

class NextOlympics:
    def countDays(self, today):
        y, m, d = today.split('.')
        y = int(y)
        m = int(m)
        d = int(d)

        ty = 2021
        tm = 7
        td = 23
        total = 0
        if y < 2021:
            total += 365 * (ty-y)
        while m < tm:
            total += [0, 31, 28, 31, 30, 31, 30][m]
            m += 1
        while tm < m:
            total -= [0, 0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31][m]
            m -= 1
        if d < td:
            total += (td-d)
        else:
            total -= (d-td)

        return total

Medium StarsInTheSky

2次元平面上の点集合が与えられるので、軸に平行な矩形で、異なる点集合を含むものを数え上げる問題です。

点の数 N <= 15 なので、総当りを行いました。矩形の判別は左下、右上の頂点対で行いました。

import itertools


class StarsInTheSky:
    def countPictures(self, N, X, Y):
        if N <= 2:
            return [0, 1, 3][N]
        count = set()
        for i in range(2, N+1):
            for p in itertools.combinations(range(N), r=i):
                xs = [X[c] for c in p]
                ys = [Y[c] for c in p]
                count.add(self.getPic(xs, ys))
        return len(count) + N

    def getPic(self, xs, ys):
        max_x, max_y = -float('inf'), -float('inf')
        min_x, min_y = float('inf'), float('inf')
        for x in xs:
            max_x = max(max_x, x)
            min_x = min(min_x, x)
        for y in ys:
            max_y = max(max_y, y)
            min_y = min(min_y, y)
        return (max_x, max_y, min_x, min_y)

Hard largestGCD

2次元グリッド上で (0,0) から (N-1,N-1) までの最短経路のうち、経路上の各頂点の最大公約数(GCD: Greatest Common Divisor)の最大値を求める問題です。

左上から右下に DPするだけ に見えます。

class SquareCityWalking:
    def largestGCD(self, N, A):
        best = [1]*(N*N)
        best[0] = A[0]
        for i in range(1, N):
            best[i] = self.gcd(best[i-1], A[i])
        for i in range(1, N):
            best[N*i] = self.gcd(best[N*i-N], A[N*i])
        for i in range(1, N):
            for j in range(1, N):
                a = A[N*i+j]
                best[N*i+j] = max(self.gcd(best[N*i+j-N], a),
                                  self.gcd(best[N*i+j-1], a))
        return best[N*N-1]


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

です。解説にある通り、途中で小さい方の値を保持する必要があるケースでWAします。

import unittest
class Test(unittest.TestCase):
    def test(self):
        scw = SquareCityWalking()
        N = 3
        A = [6, 2, 1, 3, 6, 2, 1, 2, 6]
        expected = 2
        self.assertEqual(scw.largestGCD(N, A), expected)

解になる可能性のある値を大きい順にループしながら、 その値を約数に持つマスだけを通って(0,0) から (N-1,N-1) まで最短経路で到達できるか検証するのが想定解でした。 下記がPythonによる実装です。

class SquareCityWalking:
    def largestGCD(self, N, A):
        divisors = []
        for i in range(1, int(A[0]**0.5)+1):
            if A[0] % i == 0:
                divisors.append(i)
                divisors.append(int(A[0]/i))
        divisors.sort(reverse=True)

        for d in divisors:
            visitable = [False]*(N*N)
            visitable[0] = True
            for i in range(1, N):
                visitable[i] = visitable[i-1] and A[i] % d == 0
            for i in range(1, N):
                visitable[N*i] = visitable[N*i-N] and A[N*i] % d == 0
            for i in range(1, N):
                for j in range(1, N):
                    a = A[N*i+j]
                    visitable[N*i+j] |= visitable[N*i + j - N] and a % d == 0
                    visitable[N*i+j] |= visitable[N*i + j - 1] and a % d == 0
            if visitable[N*N-1]:
                return d
        assert False, "unreachable"

確認として求めた値で経路が復元できることを確かめたのですが、それ以上の値で経路が存在しないことを確認する必要がありました。

感想

SRMには久しぶりに参加しましたが、とても楽しいコンテストでした。どの問題もサンプルが強く、これ通ったらchallengeやること無いのでは?と思っているところに上記の落とし穴がありました。

今回のコンテストでレーティングは1197->1229と上昇しました。

ちなみにTopcoderには下記リンクのARENAから参加することが出来ます。興味がある方はぜひ参加してみて下さい。 www.topcoder.com

Pythonでシフトを自動で組む

2015年にシフトを自動で組むプログラムをPythonで書いた - matsu7874のブログを書いているのですが、これはランダムで100案作って制約を満たすか確認する雑な実装なので、もう少し勉強して書き直してみるシリーズの第一回目です。

問題設定

  • N日分のシフト枠をM人のワーカーで埋める

入力

N,Mと各ワーカーのNG日程のリスト

例:5日3人の場合、下記の入力でworker0はday2がNG、worker3はday0,1,2がNGということを表します。

[
  [2,],
  [3,],
  [4,],
  [0,1,2,],
  [2,3,4,],
]

満たすべき制約

  • 各ワーカーのNGを優先する。
  • 毎日capacity人以下の割当を行う。

なるべく実現したいこと

  • 毎日の割当はなるべくcapacity人にする。
  • ワーカーの割当回数はなるべく均等にする。

方針

まずはシンプルな貪欲法で実装したいと思います。

  • リストの先頭から働ける人に働いてもらう。
  • 割当回数を均等化するために、割当日数の少ない人から割り当てていく。

最悪時は毎日全員が働けるかをチェックするため、NGが重ならない限りは割当が可能です。

実装

import sys
import random
import heapq


class Worker:
    def __init__(self, name, ng):
        self.name = name
        self.ng = ng

    def __str__(self):
        return str((self.name, str(self.ng)))


class Shift:
    def __init__(self, days, capacity):
        self.days = days
        self.capacity = capacity
        self.shift = [['UNASSIGNED']*capacity for i in range(days)]

    def assign(self, workers):
        hq = [(0, i) for i in range(len(workers))]
        used_hq = []
        for d in range(self.days):
            day_assigned_cnt = 0
            while hq and day_assigned_cnt < self.capacity:
                worker_assigned_cnt, worker_id = heapq.heappop(hq)
                if d in workers[worker_id].ng:
                    heapq.heappush(used_hq, (worker_assigned_cnt, worker_id))
                else:
                    self.shift[d][day_assigned_cnt] = workers[worker_id].name
                    day_assigned_cnt += 1
                    heapq.heappush(used_hq, (worker_assigned_cnt+1, worker_id))
            while used_hq:
                heapq.heappush(hq, heapq.heappop(used_hq))
            self.shift[d].sort()
        worker_assigned_cnt = [0 for i in range(len(workers))]
        for wa_cnt, worker_id in hq:
            worker_assigned_cnt[worker_id] = wa_cnt
        print(worker_assigned_cnt, file=sys.stderr)

    def __str__(self):
        return ',\n'.join([str((i, d)) for i, d in enumerate(self.shift)])


def main():
    n_days = 10
    n_workers = 5
    # 各workerはランダムで最大n_days/2日分までNGを出す。
    workers = [Worker('worker {:02}'.format(i), sorted(random.sample(
        range(n_days), random.randint(0, n_days//2)))) for i in range(n_workers)]
    print(*workers, sep='\n', file=sys.stderr)
    shift = Shift(n_days, 2)
    shift.assign(workers)
    for d, ws in enumerate(shift.shift):
        print(d, '\t'.join(ws), sep='\t')


if __name__ == "__main__":
    main()

実行結果

上記のプログラムを実行すると、乱数を使っているので同じ結果にはならないと思いますが、下記のような結果を得ることができます。

('worker 00', '[]')
('worker 01', '[5, 9]')
('worker 02', '[0, 9]')
('worker 03', '[1, 2, 4, 7, 9]')
('worker 04', '[1, 3, 7, 8, 9]')
[5, 4, 4, 4, 2]
0       worker 00       worker 01
1       worker 00       worker 02
2       worker 01       worker 04
3       worker 02       worker 03
4       worker 00       worker 04
5       worker 02       worker 03
6       worker 01       worker 03
7       worker 00       worker 01
8       worker 02       worker 03
9       UNASSIGNED      worker 00

4/5がNGを出しているday9以外は割当が行えていることがわかります。 一方で各ワーカーの割当回数には2から5のばらつきがあり、NGを5箇所出しているworker03,04間で割当回数が2,4と差がある点は課題です。 また、day0にworker00の代わりにworker04を割り当てることでより均等な割当が実現可能です。

最後に

ナース・スケジューリング:問題把握とモデリングを読んでいるのですが、いろいろなアルゴリズムが応用されていて非常に面白いです。 次回以降、内容を紹介しながらプログラムを改善していきたいと思います。

0-1ナップサック問題入門(if文とfor文と配列を前提とする)

この記事ではナップサック問題を解くプログラムをPythonで実装します。
再帰関数も動的計画法も使わず、if文とfor文と配列を使って解くことを目指します。

ナップサック問題とは

ナップサック問題とは次のような問題です。

いくつかのアイテムと1つのナップサックが与えられる。それぞれのアイテムには重さと価値があり、ナップサックには容量が設定されている。ナップサックの容量の範囲内でアイテム何種類か選び、選んだアイテムの価値の合計を最大化せよ。

各アイテムについてナップサックに入れられる数が0個か1個であることから0-1ナップサック問題と呼ばれることがあります。

具体例

具体例を示します。

ナップサックには 24kgまで入れられます。 アイテムの重さと価値は下記表のとおりです。

名前 重さ(kg) 価値(万円)
A 1 1
B 7 6
C 8 7
D 9 8
E 19 19

{B,C,D}の3つを選ぶと重さ24kg・価値21万円となり、これは価値の最大値になっています。

ナップサック問題の解き方

この問題は、動的計画法という高度なアルゴリズムを使うことで、商品が増えた場合も高速に解を求められることが知られていますが、今回はアイテムがたかだか20個としてしまってif文とfor文を使って素直に実装します。 ループは深いですが、やっている事自体は単純です。

capacity = 24  # ナップサックには 24kgまで入れられる
weights = [1, 7, 8, 9, 19]  # アイテムA-Eの重さ
values = [1, 6, 7, 8, 19]  # アイテムA-Eの価値
total_weight = 0
total_value = 0
max_value = 0

# それぞれ選ばない・選ぶに対応している
for use_A in (False, True):
    if use_A:  # アイテムAを選ぶなら、Aの重さと価値をナップサックに加える
        total_weight += weights[0]
        total_value += values[0]
    for use_B in (False, True):
        if use_B:  # アイテムBを選ぶなら、Bの重さと価値をナップサックに加える
            total_weight += weights[1]
            total_value += values[1]
        for use_C in (False, True):
            if use_C:  # アイテムCを選ぶなら、Cの重さと価値をナップサックに加える
                total_weight += weights[2]
                total_value += values[2]
            for use_D in (False, True):
                if use_D:
                    total_weight += weights[3]
                    total_value += values[3]
                for use_E in (False, True):
                    if use_E:
                        total_weight += weights[4]
                        total_value += values[4]

                    if total_weight <= capacity:
                        # 重さが容量以下かつ、これまでの価値の最大値を現在の価値が上回っている場合は価値の最大値を更新する
                        max_value = max(max_value, total_value)

                    if use_E:  # アイテムEを選んでいたら、Eの重さと価値をナップサックから除く
                        total_weight -= weights[4]
                        total_value -= values[4]
                if use_D:  # アイテムDを選んでいたら、Dの重さと価値をナップサックから除く
                    total_weight -= weights[3]
                    total_value -= values[3]
            if use_C:  # アイテムCを選んでいたら、Cの重さと価値をナップサックから除く
                total_weight -= weights[2]
                total_value -= values[2]
        if use_B:
            total_weight -= weights[1]
            total_value -= values[1]
    if use_A:
        total_weight -= weights[0]
        total_value -= values[0]

print(max_value)

それぞれのアイテムについて、選ばない・選ぶの各2通りのパターンをfor文をネストすることで網羅しています。2*2*2*2*2=32通りのtotal_valueを計算し、価値の最大値を更新しています。

プログラムを動的に生成する

※ここからちょっと難しくなります。

特定の問題(ナップサックの容量やアイテムの重さ・価値が固定)を解くのであれば上記で十分だろうと思いますが、上記のプログラムはアイテムの個数分のfor文のネストを深くする必要があり、アイテム数の変化に対応できません。 上記のアルゴリズムをプログラム自身に実装させてみることにします。つまりアイテムが与えられるとアイテムの個数に応じてループを増やし、全探索を行い、最大値を求める関数を作ります。 ここでのポイントは、プログラムとして成立している文字列を作成するロジックを作りexec関数を使って、文字列をプログラムとみて実行する関数を使うことです。

下記が「ナップサック問題を解くプログラムを作成して実行するプログラム」です。

import textwrap


def search_max_value(capacity: int, weights: list, values: list):
    item_count = len(weights)

    prog = [
        'total_weight = 0',
        'total_value = 0',
        'max_value = 0',
    ]

    item_loop_template = textwrap.dedent("""
        for use_{item_index} in (False, True):
            if use_{item_index}:
                total_weight += weights[{item_index}]
                total_value += values[{item_index}]
    """)
    for i in range(item_count):
        indent = ' '*4*i
        row = textwrap.indent(item_loop_template.format(item_index=i), indent)
        prog.append(row)

    indent = ' '*4*item_count
    prog.append(indent + 'if total_weight <= capacity:')
    prog.append(indent + '    max_value = max(max_value, total_value)')

    item_unloop_template = textwrap.dedent("""
        if use_{item_index}:
            total_weight -= weights[{item_index}]
            total_value -= values[{item_index}]
    """)
    for i in reversed(range(item_count)):
        indent = ' '*4*(i+1)
        row = textwrap.indent(item_unloop_template.format(item_index=i), indent)
        prog.append(row)

    local_vars = {
        'capacity': capacity,
        'weights': weights,
        'values': values,
    }
    exec('\n'.join(prog), None, local_vars)
    # print('\n'.join(prog))  # この行のコメントアウトを外すと実行されているプログラムが標準出力に出力されます
    return local_vars['max_value']

if __name__ == '__main__':
    # 具体的な問題を与える
    capacity = 24
    weights = [1, 7, 8, 9, 19]
    values = [1, 6, 7, 8, 19]
    max_value = search_max_value(capacity, weights, values)
    print(max_value)

textwrapはテキストの折り返しと詰め込みを行う標準ライブラリです。indent関数は3.3で導入されました。

※Pythonの言語仕様上の問題でアイテム数が20を超える問題は上記のプログラムで最大値を求めることができません。より大きなサイズの問題を解きたい場合は再帰関数や動的計画法を勉強しましょう。

ナップサック問題など動的計画法の勉強ように良い書籍があるので、ご興味のある方は見てみてください。

最強最速アルゴリズマー養成講座 プログラミングコンテストTopCoder攻略ガイド

最強最速アルゴリズマー養成講座 プログラミングコンテストTopCoder攻略ガイド

まとめ

この記事ではナップサック問題の典型的な導入である再帰関数・動的計画法という概念を使わず、if文とfor文を使った価値の合計の最大値を求めるプログラムを示しました。 またメタプログラミングを行うことでアイテム数が20以下のケースについて、ナップサック問題の解を求めまるプログラムを実装するプログラムを示しました。

f:id:matsu7874:20190217032530p:plain

2次元の凸包を求めるアルゴリズムと応用について

2次元の凸包(convex hull)を求めるアルゴリズムについてまとめました。また、凸包の応用先を列挙し、凸包を使って解ける競プロ問題を集めました。

この記事はデータ構造とアルゴリズム Advent Calendar 2018の17日目の記事です。 前日はゆらふな(@yurahuna)さんSentential Decision Diagramについてです。

凸包(convex hull)とは

凸包の定義

凸包(convex hull)とは,与えられた点をすべて包含する最小の凸多角形(凸多面体)のこと.

Microsoft PowerPoint - ppt-3.pptx

与えられた点集合が凸集合であるとは、その集合に属する点の任意の対を結ぶ線分がその集合に含まれることを言うのであった。与えられた集合 X に対して、その凸包は以下の同値な条件:

  1. X を含む(唯一の)最小の凸集合、
  2. X を含む凸集合全ての交わり、
  3. X に属する点から得られる凸結合全体の成す集合、
  4. X に属する点を頂点とする単体全ての合併

の何れか一つ(従って全て)を満たす集合として定義される。

凸包 - Wikipedia

凸包の特徴

  • 点集合Sに含まれる点のうち、ある軸で最小値・最大値を持つ点は凸包に含まれる
  • 点集合Sについて、凸包上の点を内部に含む3点が存在しない
  • 2次元の凸包は上部凸包(upper hull)と下部凸包(lower hull)に分けることができる

凸包を求めるアルゴリズム

ギフト包装法(1970)

ギフト包装法(Gift wrapping algorithm)やJarvisの行進法(Jarvis's march)。
半直線pqに対して、半直線qrのなす角が正かつ最小になるように全探索を行う。
n点からh点の凸包を構成する際の時間計算量はO(nh)になる。※hが小さいなら後述の分割統治より速い!

分割統治法で求める

凸包2つをマージする作業は上側接線と下側接線を求めて、内側を削除すれば良い。
n点からh点の凸包を構成する際の時間計算量はO(nlogn)になる。3次元への拡張が可能。

QuickHull(1977)

平均O(nlogn)、最悪O(n2)。4次元以上でも使えるのが強み。

グラハムスキャン(Graham's scan)(1972)

前処理としてソートを行うことで、ギフト包装法で次の点を決める際に貪欲に処理できるように工夫している。
n点からh点の凸包を構成する際の時間計算量はO(nlogn)になる。 2次元でしか使えない。

Monotone Chain(1979)

Andrew's algorithmとも。 グラハムスキャンは偏角でソートしていたが、Monotone Chainは座標でソートする。定数倍の高速化になっている。
蟻本のP233にはグラハムスキャンとして紹介されている方法はこれ。
グラハムスキャン同様に2次元でしか使えない。

def det(p, q):
    return p[0]*q[1] - p[1]*q[0]


def sub(p, q):
    return (p[0]-q[0], p[1]-q[1])


def get_convex_hull(points):
    # どの3点も直線上に並んでいないと仮定する。
    n = len(points)
    points.sort()
    size_convex_hull = 0
    ch = []
    for i in range(n):
        while size_convex_hull > 1:
            v_cur = sub(ch[-1], ch[-2])
            v_new = sub(points[i], ch[-2])
            if det(v_cur, v_new) > 0:
                break
            size_convex_hull -= 1
            ch.pop()
        ch.append(points[i])
        size_convex_hull += 1

    t = size_convex_hull
    for i in range(n-2, -1, -1):
        while size_convex_hull > t:
            v_cur = sub(ch[-1], ch[-2])
            v_new = sub(points[i], ch[-2])
            if det(v_cur, v_new) > 0:
                break
            size_convex_hull -= 1
            ch.pop()
        ch.append(points[i])
        size_convex_hull += 1

    return ch[:-1]

Chan's algorithm

O(n logh) - Chan's algorithm - Wikipedia - A Minimalist’s Implementation of the 3-d Divide-and-Conquer Convex Hull Algorithm

凸包を利用して解ける問題

凸包を利用して解ける問題の具体例

最遠点対(点集合の直径)を求める

最遠点対は凸包上の点である。 ※凸包の内部の点p,qが最遠点対だと仮定すると、半直線pqと交わる凸包上の辺abについてpa,pbの少なくともどちらかはpqより長くなることと矛盾する。

凸包上の点をいい感じに走査することで凸包の構築とは別にO(n)で最遠点対を見つけることができる。

他には

  • 当たり判定のコストを下げる

TODO

コンピュータ・ジオメトリ―計算幾何学:アルゴリズムと応用

コンピュータ・ジオメトリ―計算幾何学:アルゴリズムと応用

  • 作者: M.ドバーグ,M.ファンクリベルド,M.オーバマーズ,O.チョン,Mark De Berg,Mark Overmars,Mark van Kreveld,Otfried Cheong,浅野哲夫
  • 出版社/メーカー: 近代科学社
  • 発売日: 2010/02/01
  • メディア: 単行本
  • 購入: 4人 クリック: 34回
  • この商品を含むブログ (5件) を見る


明日はjkr_2255(すんぶ)さんが「AST」についてです。よろしくお願いします。

Pythonで文章から複合名詞や「〇〇の〇〇」といったフレーズを抽出する

複合名詞や「〇〇の〇〇」といったフレーズを、文章から取り出したい場面は多くあります。 この記事ではPythonを使ってそれらを抽出する方法を2つ紹介します。

  1. janomeのTokenFilterを使う
  2. negimaを使う

janomeのTokenFilterを使う

Janomeは@moco_betaさんが開発したPython製の形態素解析器です。 Janomeの後処理機能であるTokenFilterを使用する方法を説明します。

複合名詞用に関してはすぐに利用が可能です。 それ以外の品詞を指定してフレーズを探す場合には、TokenFilterクラスを継承して新しいFilterを自前で実装する必要があります。

利用が簡単な複合名詞を得る例を紹介します。CompoundNounFilterを使います。

from janome.tokenizer import Tokenizer
from janome.analyzer import Analyzer
from janome.tokenfilter import CompoundNounFilter


def main():
    tokenizer = Tokenizer()
    token_filter = [CompoundNounFilter()]
    analyzer = Analyzer(tokenizer=tokenizer, token_filters=token_filter)

    with open('sample.txt') as lines:
        for line in lines:
            for token in analyzer.analyze(line):
                print(token)


if __name__ == '__main__':
    main()

negimaを使う

negimaは@cocodripsさんさんが開発している品詞パターンマッチングライブラリです。

下記のようなマッチングのルールを記載したrule.csvを用意します。

id,min,max,pos0,pos1,pos2,pos3,pos4,pos5
0_xxない,1,5,名詞,,,,,
,1,5,形容詞,自立,,,,
5_xxがyy,1,5,名詞,,,,,
,1,5,助詞,格助詞,,,,
,1,5,名詞,,,,,
6_xxのyy,0,2,接頭詞,名詞接続,,,,
,1,5,名詞,,,,,
,1,5,助詞,連体化,,,,
,1,5,名詞,,,,,
999_複合名詞,2,6,名詞,,,,,

上記のルールファイルを保存したら、下記のスクリプトで該当するフレーズを抽出してみましょう。 ※まだnegimaをインストールしていない場合は、実行前にpip install -U negimaが必要です。

from negima import MorphemeMerger
mm = MorphemeMerger()
mm.set_rule_from_csv('rule.csv', sep=',')

with open('sample.txt') as lines:
    for line in lines:
        words, _ = mm.get_rule_pattern(line)
        print(words)

複合名詞や「〇〇の〇〇」といったフレーズが抽出できていることがわかります。

'暴君ディオニス', '娘さん',
'邪智暴虐の王', '村の牧人',  'メロスの裸体'

negimaは品詞を列記するだけで、文字列中から該当するフレーズを簡単に抽出することができます。

まとめ

  • 後処理を細かく実装したいならJanomeのTokenFilterが選択肢になる
  • 品詞のパターンマッチングに関してはnegimaが圧倒的に使いやすい