EmacsでPythonを書く設定 2021

emacs python

はじめに

PythonまわりのEmacsの設定を整理したということ。Company-modeの設定は言語共通の部分が多く長くなるので省略した。Language serverの紹介がメイン。

※最新の設定は2022の記事に

tam5917.hatenablog.com

Language server

Emacsからlanguage serverを使うためにeglotを入れる。MELPAからインストール可能。さらにPython用のpython-language-serverをインストールしておく。ターミナルから以下を実行する。

pip3 install python-language-server

関連パッケージをまとめて入れる場合は以下がラク。

pip install 'python-language-server[all]'

Flymake

flycheckはOFFにして、flymakeを使うようにした。flymake-diagnostic-at-pointはなかなか便利なのでおすすめ。MELPAからインストール可能。

python-black

Emacsからblackを呼び出せるようにするため、python-blackを導入する。
https://github.com/wbolster/emacs-python-black

バッファ保存時に「常に」blackをかけたいときは、以下のように設定可能である。

(declare-function python-black-on-save-mode "python-black")
(add-hook 'python-mode-hook
          #'(lambda ()
              (python-black-on-save-mode)))

python-black-on-save-modeを持ち出さない場合は、以下でもOK。

(add-hook 'python-mode-hook
          #'(lambda ()
              (add-hook 'before-save-hook
                        'python-black-buffer nil t)))

前者のminor-modeを使うとlighterの文字列が増えてゴチャゴチャするので、個人的には後者が好み。
projectがblackによるフォーマットを必要とするかに応じて振る舞いを変えたい場合はpython-black-on-save-mode-enable-dwimを使う。このようにユーザに選択の余地が残されているのがよい。

設定

;; Python
(setq auto-mode-alist (cons '("\\.py$" . python-mode) auto-mode-alist))
(setq interpreter-mode-alist (cons '("python" . python-mode)
                                   interpreter-mode-alist))

;; pip3 install python-language-server
;; https://github.com/palantir/python-language-server
;;
;; 以下は追加で（お好みで）インストール
;; pip install 'python-language-server[all]' としてしまうのが早い
;; 
;; pip3 install rope ; completions and renaming
;; pip3 install pyflakes ; linter to detect various errors
;; pip3 install mccabe ; linter for complexity checking
;; pip3 install pycodestyle ; linter for style checking
;; pip3 install pydocstyle ; linter for docstring style checking
;; pip3 install pyls-mypy ; type checker (heavy)
;; pip3 install pyls-isort ; import sort code formatting
;; pip3 install pyls-black ; code formatting using Black

(require 'eglot)
(add-hook 'python-mode-hook 'eglot-ensure)

;; server として起動
;; (add-to-list 'eglot-server-programs
;;              `(python-mode . ("pyls" "-v" "--tcp" "--host"
;;                               "localhost" "--port" :autoport)))

;; 保存時に自動整形 (black)
;; python-blackをmelpaからインストールしておく
(add-hook 'python-mode-hook
          #'(lambda ()
              (add-hook 'before-save-hook
                        'python-black-buffer nil t)))

;; 以下の設定でもOKだが、minor-modeの分、lighterが増えてごちゃごちゃする
;; (declare-function python-black-on-save-mode "python-black")
;; (add-hook 'python-mode-hook
;;           #'(lambda ()
;;               (python-black-on-save-mode)))

;; 保存時にインポート順を自動修正
;; py-isortをmelpaからインストールしておく
(add-hook 'python-mode-hook #'(lambda ()
                                (add-hook 'before-save-hook
                                          'py-isort-buffer nil t)))

;; point位置にメッセージを表示する。eldocと干渉しないのが利点
(with-eval-after-load "flymake"
  (require 'flymake-diagnostic-at-point)
  (add-hook 'flymake-mode-hook #'flymake-diagnostic-at-point-mode))

その他

ほか、以下は好みで設定した。エコーエリアが複数行になるのを防ぐ効果がある。

(custom-set-variables
  '(eldoc-echo-area-use-multiline-p nil))

またsmart-jumpやivy-xrefも便利なので使っている。

(require 'smart-jump)
(smart-jump-setup-default-registers)  ;; pythonはデフォルトで有効になる
;; M-. (smart-jump-go; xref-find-definitions)
;; M-, (smart-jump-back; xref-pop-marker-stack)
;; M-? (smart-jump-references; xref-find-references)

(require 'ivy-xref)
(setq xref-show-xrefs-function #'ivy-xref-show-xrefs)

もし、blackユーザならばpycodestyleの設定をしておくと良い。

[pycodestyle]
ignore = E203, W503, W504
max-line-length = 88

おわりに

eglotとpython-language-serverの組み合わせがなかなか快適だったので紹介した。lsp-modeを設定するのも悪くはないし、elpyとflycheckの組み合わせで書くのも悪くはないので（特に軽快性が優れる）、このあたりは好みだろう。

参考

Griffin-Limアルゴリズムの実行時間をlibrosaとtorchaudioで比較してみた話

音声 python

はじめに

音声の振幅スペクトルから位相を推定し、元の音声を復元するためのGriffin-Limアルゴリズムが知られている。
Griffin-Limアルゴリズムはlibrosaパッケージとtorchaudioパッケージの両方に実装されている。

  • librosa

librosa.org

  • torchaudio

pytorch.org

今回はGriffin-Limアルゴリズムに関してCPUとGPUとで実行時間を計測し比較してみた。

環境

  • OS、ハードウェア、ドライバ
  • ソフトウェア
    • Python 3.6.9
    • librosa 0.8.0
    • torch 1.7.1
    • torchaudio 0.7.2
    • torchlibrosa 0.0.9

スクリプト

"""
PyTorch script for verification of execution time.
"""

# Commentary:
# torch 1.8.1
# torchaudio 0.8.1
# torchlibrosa 0.0.9

# GeForce RTX 2070
# CUDA 11.2
# CUDA Driver 460.32.03

import time

import librosa
import numpy as np
import torch
import torchaudio
import torchaudio.transforms as T

N_FFT = 1024
WIN_LENGTH = None
HOP_LENGTH = 512
N_ITER = 50

if __name__ == "__main__":

    waveform, sample_rate = librosa.load(librosa.ex("trumpet"))

    # librosa
    start = time.time()
    specgram = np.abs(
        librosa.stft(
            waveform, n_fft=N_FFT, hop_length=HOP_LENGTH, win_length=WIN_LENGTH
        )
    )
    y_inv = librosa.griffinlim(
        specgram, n_iter=N_ITER, hop_length=HOP_LENGTH, win_length=WIN_LENGTH
    )
    elapsed_time = time.time() - start
    print("elapsed_time (librosa): {0:.6f}".format(elapsed_time) + "[sec]")

    # load waveform as torch.Tensor
    waveform, sample_rate = torchaudio.load(librosa.ex("trumpet"))
    waveform_cuda = waveform.cuda()

    # torchaudio (CPU)
    specgram = T.Spectrogram(
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
    )
    griffin_lim = T.GriffinLim(
        n_fft=N_FFT,
        n_iter=N_ITER,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
    )
    start = time.time()
    reconstructed = griffin_lim(specgram(waveform))
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; CPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchaudio (GPU)
    specgram = T.Spectrogram(
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
    ).cuda()
    griffin_lim = T.GriffinLim(
        n_fft=N_FFT,
        n_iter=N_ITER,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
    ).cuda()
    torch.cuda.synchronize()
    start = time.time()
    reconstructed = griffin_lim(specgram(waveform_cuda))
    torch.cuda.synchronize()
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; GPU): {0:.6f}".format(elapsed_time) + "[sec]")

実行結果

GPUを使ったtorchaudioが一番早い結果となった。

elapsed_time (librosa): 0.264543[sec]
elapsed_time (torchaudio; CPU): 0.147203[sec]
elapsed_time (torchaudio; GPU): 0.022305[sec]

まとめ

位相復元に関しても、torchaudioのCPU版を使うだけで実行時間の短縮が望めそうである。GPU版も悪くない結果である。今回は内部的に実装されているアルゴリズムを考慮したフェアな比較というよりは、いちユーザが関数を外から触るときの体感を重視した。

音声認識結果を音声合成するPythonスクリプトをSpeechRecognitionとPyOpenJTalkで書いたみた話

音声 python

はじめに

かつて、音声認識と音声合成を組み合わせて遊んでみるという主旨の記事を書いたことがある。

tam5917.hatenablog.com

音声合成には、コマンドラインから音声合成できるOpenJTalkパッケージを用いたのだった。これをPythonから動かす場合には、専用の関数を定義して呼び出す必要があり、ひと手間かかる。

最近、山本りゅういち氏がPyOpenJTalkなるパッケージを整備し、単体でTTS（Text-to-Speech、テキストを音声に変換する）機能を実現した。

これを使って前回記事のPythonスクリプトを少しきれいにしたというのが本記事の主旨である。

インストール

音声認識にはSpeechRecognitionパッケージを使う。PyOpenJTalkが今回のキモなので忘れずにインストールする。最新版をインストールするという意味で--upgradeオプションもつけておく。

pip3 install SpeechRecognition
pip3 install --upgrade pyopenjtalk

スクリプト

以下のスクリプトを適当な名前で保存し、実行する。 「話しかけてみましょう！」と表示されたのち、適当な言葉をしゃべると、認識された言葉をオウム返しで音声合成して終了するというシンプルなものである。

import subprocess

import numpy as np
import pyopenjtalk
import speech_recognition as sr
from scipy.io import wavfile

OUT_WAV = "/tmp/tts.wav"

if __name__ == "__main__":

    # マイクからの音声入力
    r = sr.Recognizer()
    with sr.Microphone() as source:
        print("話しかけてみましょう！")
        audio = r.listen(source)

    try:
        # 日本語でGoogle音声認識
        text = r.recognize_google(audio, language="ja-JP")
    except sr.UnknownValueError:
        print("Google音声認識は音声を理解できませんでした。")
    except sr.RequestError as e:
        print("Google音声認識サービスからの結果を要求できませんでした;" " {0}".format(e))
    else:
        print(text)

        # 音声認識結果を音声合成する
        x, sr = pyopenjtalk.tts(text)

        # wavファイルに保存
        wavfile.write(OUT_WAV, sr, x.astype(np.int16))

        # 保存したwavファイルを再生
        subprocess.run("afplay " + OUT_WAV, shell=True)

おわりに

PyOpenJTalkのおかげでTTSが手軽に実現でき、楽しい。PyOpenJTalkには話速変換やピッチ変換のオプションもあるので、調べて使ってみるのが良いだろう。

おまけ

話速変換とピッチ変換のオプションを使ったnotebookはこちらから。 nbviewer.jupyter.org

注意点

実行時のエラー

File "pyopenjtalk/htsengine.pyx", line 1, in init pyopenjtalk.htsengine
ValueError: numpy.ndarray size changed, may indicate binary incompatibility. Expected 88 from C header, got 80 from PyObject

このようなエラーが出たら、numpyを再インストールしてみるのが吉。

Google音声認識API

ただし、本記事で利用しているGoogle音声認識APIはあくまでテスト/デモ用であり、いつGoogle側にAPI利用を打ち切られても文句は言えない。Google製の音声認識APIを使いたい場合はCould Speech-to-Textの認証情報を取得して使うのが本筋である。ひと月あたり60分の音声までなら利用無料であるが、それ以上を利用するにはお金がかかる。

cloud.google.com

torchaudioとtorchlibrosaの実行速度に違いはあるのか？

音声 python

はじめに

PyTorchには音声系データを処理するのに便利なtorchaudioというライブラリが存在する。
pytorch.org

一方、音声系データの処理に便利なlibrosaというパッケージが存在する。
librosa.org

さらにtorchlibrosaという、librosa内部の行列計算まわりをPyTorchで置き換えたパッケージが存在する。
github.com

ここで一つ疑問：
「で、結局どれ使えばいいの？（どれが早いの？）」

この疑問が気になったので、CPU実行とGPU実行における実行時間を比較検証するための簡易的なスクリプトを書いて調べてみたということである。

環境

  • OS、ハードウェア、ドライバ
  • ソフトウェア
    • Python 3.6.9
    • librosa 0.8.0
    • torch 1.7.1
    • torchaudio 0.7.2
    • torchlibrosa 0.0.9

スクリプト

メルスペクトログラムを計算するスクリプトである。

  • 3月24日追記：GPUの計測方法に誤りがあったため、コメントに従ってtorch.cuda.synchronize()を用いるように修正した。

スクリプトを表示する

"""
PyTorch script for verification of execution time.
"""

# Commentary:
# torch 1.7.1
# torchaudio 0.7.2
# torchlibrosa 0.0.9

# GeForce RTX 2070
# CUDA 11.2
# CUDA Driver 460.32.03

import time

import librosa
import librosa.feature
import torch
import torchaudio
import torchaudio.transforms as T
import torchlibrosa as tl

N_FFT = 1024
WIN_LENGTH = None
HOP_LENGTH = 512
N_MELS = 128

if __name__ == "__main__":

    waveform, sample_rate = librosa.load(librosa.ex("nutcracker"))
    # librosa
    start = time.time()
    mel_spectrogram = librosa.feature.melspectrogram(
        y=waveform,
        sr=sample_rate,
        n_fft=N_FFT,
        hop_length=HOP_LENGTH,
        n_mels=N_MELS,
        power=2.0,
    )
    elapsed_time = time.time() - start
    print("elapsed_time (librosa): {0:.6f}".format(elapsed_time) + "[sec]")

    waveform, sample_rate = torchaudio.load(librosa.ex("nutcracker"))
    waveform_cuda = waveform.cuda()

    # torchaudio (CPU)
    start = time.time()
    mel_spectrogram = T.MelSpectrogram(
        sample_rate=sample_rate,
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
        power=2.0,
        n_mels=N_MELS,
    )
    melspec = mel_spectrogram(waveform)
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; CPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchaudio (GPU)
    torch.cuda.synchronize()
    start = time.time()
    mel_spectrogram = T.MelSpectrogram(
        sample_rate=sample_rate,
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
        power=2.0,
        n_mels=N_MELS,
    )
    mel_spectrogram = mel_spectrogram.cuda()
    melspec = mel_spectrogram(waveform_cuda)
    torch.cuda.synchronize()
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; GPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchlibrosa (CPU)
    # TorchLibrosa feature extractor the same as librosa.feature.melspectrogram()
    start = time.time()
    feature_extractor = torch.nn.Sequential(
        tl.Spectrogram(
            n_fft=N_FFT,
            win_length=WIN_LENGTH,
            hop_length=HOP_LENGTH,
            power=2.0,
        ),
        tl.LogmelFilterBank(
            n_fft=N_FFT,
            sr=sample_rate,
            n_mels=N_MELS,
            is_log=False,  # Default is true
        ),
    )
    logmel = feature_extractor(waveform)
    elapsed_time = time.time() - start
    print("elapsed_time (torchlibrosa; CPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchlibrosa (GPU)
    # TorchLibrosa feature extractor the same as librosa.feature.melspectrogram()
    torch.cuda.synchronize()
    start = time.time()
    feature_extractor = torch.nn.Sequential(
        tl.Spectrogram(
            n_fft=N_FFT,
            win_length=WIN_LENGTH,
            hop_length=HOP_LENGTH,
            power=2.0,
        ),
        tl.LogmelFilterBank(
            n_fft=N_FFT,
            sr=sample_rate,
            n_mels=N_MELS,
            is_log=False,  # Default is true
        ),
    )
    feature_extractor = feature_extractor.cuda()
    logmel = feature_extractor(waveform_cuda)
    torch.cuda.synchronize()
    elapsed_time = time.time() - start
    print("elapsed_time (torchlibrosa; GPU): {0:.6f}".format(elapsed_time) + "[sec]")


  • 5月8日追記：koukyo1213様のコメントに従い、torch.nn.Sequentialやtorchaudio.transforms.MelSpectrogramのインスタンス化の部分を計測対象から外し、計算処理のみにフォーカスしたスクリプトは以下である。

スクリプトを表示する

"""
PyTorch script for verification of execution time.
"""

# Commentary:
# torch 1.8.1
# torchaudio 0.8.1
# torchlibrosa 0.0.9

# GeForce RTX 2070
# CUDA 11.2
# CUDA Driver 460.32.03

import time

import librosa
import librosa.feature
import torch
import torchaudio
import torchaudio.transforms as T
import torchlibrosa as tl

N_FFT = 1024
WIN_LENGTH = None
HOP_LENGTH = 512
N_MELS = 128

if __name__ == "__main__":

    torchaudio.set_audio_backend("sox_io")

    waveform, sample_rate = librosa.load(librosa.ex("nutcracker"))
    # librosa
    start = time.time()
    mel_spectrogram = librosa.feature.melspectrogram(
        y=waveform,
        sr=sample_rate,
        n_fft=N_FFT,
        hop_length=HOP_LENGTH,
        n_mels=N_MELS,
        power=2.0,
    )
    elapsed_time = time.time() - start
    print("elapsed_time (librosa): {0:.6f}".format(elapsed_time) + "[sec]")

    waveform, sample_rate = torchaudio.load(librosa.ex("nutcracker"))
    waveform_cuda = waveform.cuda()

    # torchaudio (cpu)
    mel_spectrogram = T.MelSpectrogram(
        sample_rate=sample_rate,
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
        power=2.0,
        n_mels=N_MELS,
    )
    start = time.time()
    melspec = mel_spectrogram(waveform)
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; CPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchaudio (GPU)
    mel_spectrogram = T.MelSpectrogram(
        sample_rate=sample_rate,
        n_fft=N_FFT,
        win_length=WIN_LENGTH,
        hop_length=HOP_LENGTH,
        power=2.0,
        n_mels=N_MELS,
    )
    torch.cuda.synchronize()
    start = time.time()
    mel_spectrogram = mel_spectrogram.cuda()
    melspec = mel_spectrogram(waveform_cuda)
    torch.cuda.synchronize()
    elapsed_time = time.time() - start
    print("elapsed_time (torchaudio; GPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchlibrosa (CPU)
    # TorchLibrosa feature extractor the same as librosa.feature.melspectrogram()
    feature_extractor = torch.nn.Sequential(
        tl.Spectrogram(
            n_fft=N_FFT,
            win_length=WIN_LENGTH,
            hop_length=HOP_LENGTH,
            power=2.0,
        ),
        tl.LogmelFilterBank(
            n_fft=N_FFT,
            sr=sample_rate,
            n_mels=N_MELS,
            is_log=False,  # Default is true
        ),
    )
    start = time.time()
    logmel = feature_extractor(waveform)
    elapsed_time = time.time() - start
    print("elapsed_time (torchlibrosa; CPU): {0:.6f}".format(elapsed_time) + "[sec]")

    # torchlibrosa (GPU)
    # TorchLibrosa feature extractor the same as librosa.feature.melspectrogram()
    feature_extractor = torch.nn.Sequential(
        tl.Spectrogram(
            n_fft=N_FFT,
            win_length=WIN_LENGTH,
            hop_length=HOP_LENGTH,
            power=2.0,
        ),
        tl.LogmelFilterBank(
            n_fft=N_FFT,
            sr=sample_rate,
            n_mels=N_MELS,
            is_log=False,  # Default is true
        ),
    )
    torch.cuda.synchronize()
    start = time.time()
    feature_extractor = feature_extractor.cuda()
    logmel = feature_extractor(waveform_cuda)
    torch.cuda.synchronize()
    elapsed_time = time.time() - start
    print("elapsed_time (torchlibrosa; GPU): {0:.6f}".format(elapsed_time) + "[sec]")


実行結果

その1（3月24日版）

GPUを使ったtorchaudioが一番早い結果となった。

elapsed_time (librosa): 0.042655[sec]
elapsed_time (torchaudio; CPU): 0.021293[sec]
elapsed_time (torchaudio; GPU): 0.002648[sec]
elapsed_time (torchlibrosa; CPU): 0.237388[sec]
elapsed_time (torchlibrosa; GPU): 0.209118[sec]

その2（5月8日版）

koukyo1213様のコメントに従って、torch.nn.Sequentialやtorchaudio.transforms.MelSpectrogramのインスタンス化部分を計測対象から外した場合の結果を以下に示す。

elapsed_time (librosa): 0.040965[sec]
elapsed_time (torchaudio; CPU): 0.014024[sec]
elapsed_time (torchaudio; GPU): 0.001927[sec]
elapsed_time (torchlibrosa; CPU): 0.027673[sec]
elapsed_time (torchlibrosa; GPU): 0.006538[sec]

今回もGPUを使ったtorchaudioが一番早い結果となった。

考察とまとめ

最初のスクリプトでtorchlibrosaの実行時間がより多くかかったのは、koukyo1213様ご指摘の通り、インスタンス化のオーバーヘッドがかなり大きいことに起因すると考えられる。その影響を除いた新たなスクリプトでは、GPU使用のtorchlibrosaは、同じくGPU使用のtorchaudioと比較して3倍ほど低速という結果が得られた。とはいえ実行時間の桁は揃い、どちらも十分実用に供すると考える。

black用py-isortの設定

emacs python

私はblackをpythonコードのformatterに用いているので、blackの仕様に合うようにisortを設定して動かしたいというわけである。isortをEmacsから触るためにはpy-isortを用いる。以下の記事は参考になる。

qiita.com

isortのマニュアルはこちらから： pycqa.github.io

まずisortをインストールする。

pip3 install isort

py-isortはEmacsのパッケージシステムから容易にインストール可能である。その後、Emacs設定ファイルに以下を追記する。こうするとバッファ保存時にisortが走ってくれるので便利である（＝つまりインポート順を自分で編集する必要がなくなる！）。

;; py-sort: isortを使ってライブラリのimport順をソート
(eval-when-compile (require 'py-isort))
;; for black
(setq py-isort-options '("--multi-line=3"
                         "--line-length=88"
                         "--trailing-comma"
                         "--use-parentheses"
                         "--ensure-newline-before-comments"
                         "--force-grid-wrap=0"))
;; バッファ保存時に自動ソート
(add-hook 'before-save-hook 'py-isort-before-save)

なおこれらのオプション値はblack公式のドキュメントを参考にした。 black.readthedocs.io

numpy.arange関数におけるstepに0.1を指定したときの振る舞いとその原因を追ってみた話

python

はじめに

先日、接点QB氏による以下のツイートを見かけた：

興味深いので少し調べてみたというわけである（少々長い記事）。

numpyのarange関数の振る舞い

numpyのarange関数のステップ数に0.1を指定したときの振る舞いについて、

np.arange(1.0, 1.5, 0.1)

とすると、0.1刻みで1.0, 1.1, 1.2, 1.3, 1.4までのarrayが得られた。array要素の最大値1.4は引数stopの1.5より小さいというのが、想定される挙動である。

一方、

np.arange(1.0, 1.6, 0.1)

とすると、0.1刻みで1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6のように、array要素の最大値が1.5を超えて1.6に達したので困惑したということである（「1.0, 1.1, 1.2, 1.3, 1.4, 1.5」かと思ったのに！）。

原因を探るため、numpyのマニュアルのarangeの項を見てみた。 numpy.org

すると、戻り値のarrayについて以下の記述が見つかった：

For floating point arguments, the length of the result is ceil((stop - start)/step).

つまり、浮動小数点数をarange関数の引数に与えた場合、array長がceil（＝小数点以下切り上げ）を伴う式で計算されるということだった。

手元で挙動を試してみると、以下のようになった。

つまり、(1.5 - 1.0) / 0.1の計算結果が5.0となる一方で、(1.6 - 1.0) / 0.1の計算結果が6.000000000000001となり6.0にならない。

したがってceil関数を通したとき、

ceil((1.5 - 1.0) / 0.1) は5になるが、ceil((1.5 - 1.0) / 0.1) は7になった

というわけである。

問題は「なぜこのような計算結果が得られたのか」ということである。原因を一言で言えば「浮動小数点数演算に由来する計算誤差」ということに尽きるが、本記事ではこれを詳細に追ってみた（浮動小数点数の復習問題！）。

浮動小数点数の表現形式と演算

さて、浮動小数点数のIEEE754倍精度形式は次式で表される。 $$(-1)^{sgn} \times 2^{(exp - 1023)} \times (1+frac×2^{-52})$$ $sgn$は符号部、$exp$は指数部、$frac$は仮数部である（ここでは10進数表現）。

この表現形式を元にして、計算を追っていこうという方針である。 ここでは以下の記事を参考にした。 qiita.com

浮動小数点数の倍精度表現

検証にあたり、接点QB氏の具体例に基づいて、対象とする浮動小数点数を

0.1, 1.0, 1.5, 1.6

の4つに限定する。

0.1の倍精度表現

0.1の倍精度表現は

符号部 = 0
指数部 = 01111111011
仮数部 = 1001100110011001100110011001100110011001100110011010

となり、これらを10進数に変換すると

符号部 = 0
指数部 = 1019
仮数部 = 2702159776422298

となる。0.1の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1019 - 1023} \times (1 + 2702159776422298 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \times 2^{-4} \times (2^{52} + 2702159776422298) \times 2^{-52} \\ &= (4503599627370496 + 2702159776422298) \times 2^{-56} \\ &= 7205759403792794 / 72057594037927936 \end{aligned} $$

1.0の倍精度表現

1.0の倍精度表現は

符号部 = 0
指数部 = 01111111111
仮数部 = 0000000000000000000000000000000000000000000000000000

であり、10進数に変換すると

符号部 = 0
指数部 = 1023
仮数部 = 0

となる。1.0の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1023 - 1023} \times (1 + 0 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \\ &= \frac{2^{52}}{2^{52}} \\ &= 4503599627370496 / 4503599627370496 \end{aligned} $$

1.5の倍精度表現

1.5の倍精度表現は

符号部 = 0
指数部 = 01111111111
仮数部 = 1001100110011001100110011001100110011001100110011010

となり、これらを10進数に変換すると

符号部 = 0
指数部 = 1023
仮数部 = 2251799813685248

となる。 1.5の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1023 - 1023} \times (1 + 2251799813685248 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \times 2^{0} \times (2^{52} + 2251799813685248) \times 2^{-52} \\ &= (4503599627370496 + 2251799813685248) \times 2^{-52} \\ &= 6755399441055744 / 4503599627370496 \end{aligned} $$

1.6の倍精度表現

1.6の倍精度表現は

符号部 = 0
指数部 = 01111111111
仮数部 = 1001100110011001100110011001100110011001100110011010

となり、これらを10進数に変換すると

符号部 = 0
指数部 = 1023
仮数部 = 2702159776422298

となる。 1.6の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1023 - 1023} \times (1 + 2702159776422298 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \times 2^{0} \times (2^{52} + 2702159776422298) \times 2^{-52} \\ &= (4503599627370496 + 2702159776422298) \times 2^{-52} \\ &= 7205759403792794 / 4503599627370496 \end{aligned} $$

倍精度表現に基づく浮動小数点数演算

以上で倍精度表現が揃ったので、検証に入ることができる（前置き長すぎワロタ）。 表に示しておく。

浮動小数点数 倍精度表現
0.1 7205759403792794 / 72057594037927936
1.0 1 (=4503599627370496 / 4503599627370496)
1.5 6755399441055744 / 4503599627370496
1.6 7205759403792794 / 4503599627370496

以上の表に基づいて、(1.5 - 1.0) / 0.1 を計算すると

((6755399441055744  - 4503599627370496) / 4503599627370496) / (7205759403792794 / 72057594037927936)
= 5

が得られた。ぴったし5である。

一方、(1.6 - 1.0) / 0.1 については

((7205759403792794 - 4503599627370496) / 4503599627370496) / (7205759403792794 / 72057594037927936)
= 6.000000000000001

と計算できた。

以上の計算結果に基づいてceilは小数点以下を切り上げるため、冒頭で示した挙動が説明できたわけである。

まとめ

本記事ではarange関数の引数のひとつであるstep数に浮動小数点数を適用した場合に、 戻り値のarrayの要素数がどのように変わりうるかを、具体的な数値例を元に調べてみた。

Numpyのマニュアルによれば、arangeによる戻り値のarrayの要素数はceil((stop - start)/step)で計算され、これらstart, stop, stepはarange関数の引数として与えられる。

接点QB氏のツイートをもとに、arangeの引数が具体的に

  • (1) start=1.0, stop=1.5, step=0.1
  • (2) start=1.0, stop=1.6, step=0.1

の場合に、浮動小数点数の倍精度表現に基づいて要素数を計算し、arange関数の振る舞いを検証した。 検証の結果、確かにPythonの計算結果は正しかったことが確かめられた。

Numpyマニュアルにも浮動小数点数が引数の場合はlinspace使ったほうがいいよ的な記述もあり、arangeを使うときには思い込みによるバグを埋め込まないように気をつけたい。

参考

vscodeのメモ

  • ターミナルの表示およびカーソル切り替えはCommand + Shift + @
  • ターミナルを閉じる場合はControl + Shift + @ （ターミナルにフォーカスされている場合）
  • サイドバーを閉じる場合はControl + Option + Space (Toggle Sidebar Visibility)