論文

【タイトル】 Enhancing the Spatial Resolution of Stereo Images using a Parallax Prior (CVPR2018)
【概要】

• Stereo Imageを用いた超解像手法の提案
• 白黒画像の生成→カラー画像の生成のtwo stage学習
• 視差を表現したシフト画像の利用

開発環境

• Windows10
  
• Pytorch

データセット

Flickr1024: A Dataset for Stereo Image Super-resolution (https://yingqianwang.github.io/Flickr1024/)

ソースコード

https://github.com/kyoukuntaro/StereoSR_Jeon_CVPR2018

論文

【タイトル】 Enhancing the Spatial Resolution of Stereo Images using a Parallax Prior (CVPR2018)
【概要】

• Stereo Imageを用いた超解像手法の提案
• 白黒画像の生成→カラー画像の生成のtwo stage学習
• 視差を表現したシフト画像の利用

開発環境

• Windows10
  
• Pytorch

データセット

Flickr1024: A Dataset for Stereo Image Super-resolution (https://yingqianwang.github.io/Flickr1024/)

ソースコード

https://github.com/kyoukuntaro/StereoSR_Jeon_CVPR2018

はじめに

突然ですが、なんの脈絡もなく、マンデルブロ集合を描きたくなりました。
マンデルブロ集合 (Wikipedia)

↓こんなやつ↓

ただ描くだけだと面白くないので、Google Mapみたいにスクロールしたり拡大したり…というのをやってみよう…と思いちょっと手を動かしてみました。

先に完成品…

先に完成品をリンクしておきます。目標としては、このリンク先のようなことをやりたい、ということになります。

技術選定

Google Mapみたいな地図のアプリを作る、となると、いくつか選択肢はありますが、今回はタイル画像を用意してLeafletを使ってブラウザ上に描画する、という方法をとってみることにしました。

タイル画像

地図アプリで一般的なのが、縮尺、緯度、経度ごとに地図の画像を用意して、それをつなぎ合わせて表示する、という方法。国土地理院のサイトでわかりやすく説明されていますので詳細はそちらを参照ください。
地理院地図 | 地理院タイル仕様

上記リンク先で図示されているような形でズームレベルごとのタイル分割してみると、こんな感じになります。

一般的にタイル1枚の画像のサイズは256x256ピクセルにします。(Retina用だと512x512ピクセルとかを使うこともあるようです)

この画像を{z}/{x}/{y}.pngというパスに配置します。(zはズームレベル、x,yはそのズームレベルでの画像の位置)

Leaflet

詳細はLeafletのサイトを参照のこと。

日本語での説明は、Wikipediaによると

Leaflet は広く使われているWeb地図のためのJavaScriptライブラリである。 2011年に最初にリリースされた。 モバイルとデスクトップのプラットフォームのほとんどに対応し、HTML5とCSS3に対応している。 OpenLayersやGoogle Maps API（英語版）とともに最も人気のあるJavaScript地図ライブラリの一つであり、FourSquare、Pinterest、Flickrなどの有名なサイトで使われている。

ということになります。
上記のようなタイル画像を使ってWeb上で地図表示をするのによく使われています。(地理院タイルを使って地図表示とかも簡単)

地図表示用に使われるLeafletを今回のように地図でない画像表示に使う方法に関しては、この記事の続きで記載したいと思います。

マンデルブロ集合のタイル画像の作成

で、マンデルブロ集合のタイル画像の作成に関してです。

-2-2i〜2+2iの範囲をズームレベル0で描画するとすると、-2〜2の4の範囲を256分割で計算するわけなので、1ピクセルあたり4/256=1/64単位で計算することになります。
ズームレベル1だとタイル1枚あたり、2の範囲を256分割になりますので、1ピクセルあたり2/256=1/128、ズームレベル2だと1/256…という形で計算単位を細かくしていきます。

実際にPythonで計算してPNGに保存するソースはこちら。githubにもおいてあります。

mandelbrot.py

"""マンデルブロ集合を描画する."""
from argparse import ArgumentParser
from itertools import product
from os import makedirs, path

import numpy
from PIL import Image
from tqdm import tqdm


def main():
    """メイン関数."""
    parser = ArgumentParser()
    parser.add_argument('outdir')
    parser.add_argument('level', type=int)
    args = parser.parse_args()
    outpath = path.join(args.outdir, str(args.level))
    tile_num = 2 ** args.level
    tile_wid = 4.0 / tile_num

    lefttop = -2.0+2.0j
    delta = tile_wid / 256

    xylist = list(product(range(tile_num), range(tile_num)))
    for x, y in tqdm(xylist):
        img = get_image(lefttop + x * tile_wid - (y * tile_wid) * 1j, delta)
        outpath = path.join(args.outdir, str(args.level), str(x))
        makedirs(outpath, exist_ok=True)
        imgname = path.join(outpath, f'{y}.png')
        img.save(imgname)


def get_image(lefttop, delta):
    """マンデルブロ集合のpillowイメージを取得する

    Arguments:
        lefttop {complex} -- 左上
        delta {float} -- 1pxの差

    Returns:
        Image -- 作成したイメージ
    """
    imgarray = numpy.zeros((256, 256), dtype=numpy.uint8)

    for rstep, istep in product(range(256), range(256)):
        c = lefttop + rstep * delta - (istep * delta) * 1j
        z = 0
        for i in range(256):
            z = z ** 2 + c
            if abs(z) > 2:
                imgarray[istep][rstep] = 255 - i
                break

    img = Image.fromarray(imgarray)
    return img


if __name__ == "__main__":
    main()

1ピクセルごと愚直に計算する処理にしてます。

python mandelbrot.py docs/tiles/grayscale 2

とすることで、ズームレベル2のタイルをdocs/tiles/grayscaleディレクトリの下に出力する…という動作になります。

最後に

まずはタイル画像の生成までをこの記事に書きました。
そのタイル画像をLeafletを使って実際に表示するところに関しては続きの記事に記載したいと思います。
↓
続きを書きました。

概要

• 様々なプログラミング言語でクラス (あるいはクラスもどき) を書く

Java

MyCounter.java

public class MyCounter {

  public static MyCounter globalCounter = null;

  private String name = "No Name";
  private int count;

  public MyCounter(int initValue) {
    count = initValue;
  }

  public String getName() {
    return name;
  }

  public void setName(String name) {
    this.name = name;
  }

  public int addCount(int number) {
    count += number;
    return count;
  }

  public void print() {
    System.out.println(format());
  }

  private String format() {
    return name + ": " + count;
  }

  public static void createGlobalCounter(String name, int initValue) {
    globalCounter = new MyCounter(initValue);
    globalCounter.setName(name);
  }
}

MyMain.java

public class MyMain {

  public static void main(String[] args) {

    MyCounter counter = new MyCounter(0);
    counter.print();

    counter.setName("John Doe");
    String currentName = counter.getName();
    System.out.println("Current Name=" + currentName);

    int currentValue = counter.addCount(10);
    System.out.println("Current Value=" + currentValue);

    counter.print();

    MyCounter.createGlobalCounter("CafeBabe! Write once, run anywhere", 256);
    MyCounter.globalCounter.print();
  }
}

実行結果

No Name: 0
Current Name=John Doe
Current Value=10
John Doe: 10
CafeBabe! Write once, run anywhere: 256

JavaScript

my_counter.js

class MyCounter {

  constructor(initValue) {
    this._name = 'No Name';
    this.count = initValue;
  }

  get name() {
    return this._name;
  }

  set name(name) {
    this._name = name;
  }

  addCount(number) {
    this.count += number;
    return this.count;
  }

  print() {
    console.log(this._format());
  }

  _format() {
    return this._name + ': ' + this.count;
  }

  static createGlobalCounter(name, init_value) {
    MyCounter.global_counter = new MyCounter(init_value);
    MyCounter.global_counter.name = name;
  }
}

module.exports = MyCounter;

my_main.js

const MyCounter = require('./my_counter');

counter = new MyCounter(0);
counter.print();

counter.name = 'John Doe';
const currentName = counter.name;
console.log('Current Name=' + currentName);

const currentValue = counter.addCount(10);
console.log('Current Value=' + currentValue);

counter.print();

MyCounter.createGlobalCounter('JavaScript was born from Netscape Navigator 2.0', 256);
MyCounter.global_counter.print();

実行結果

$ node my_main.js
No Name: 0
Current Name=John Doe
Current Value=10
John Doe: 10
JavaScript was born from Netscape Navigator 2.0: 256

Python

my_counter.py

class MyCounter:

  global_counter = None

  def __init__(self, init_value):
    self.__name = "No Name"
    self.__count = init_value

  @property
  def name(self):
    return self.__name

  @name.setter
  def name(self, name):
    self.__name = name

  def add_count(self, number):
    self.__count += number
    return self.__count

  def print(self):
    print(self.__format())

  def __format(self):
    return "{0}: {1}".format(self.__name, self.__count)

  @classmethod
  def create_global_counter(cls, name, init_value):
    cls.global_counter = MyCounter(init_value)
    cls.global_counter.name = name

my_main.py

from my_counter import MyCounter

counter = MyCounter(0)
counter.print()

counter.name = "John Doe"
current_name = counter.name
print("Current Name={}".format(current_name))

current_value = counter.add_count(10)
print("Current Value={}".format(current_value))

counter.print()

MyCounter.create_global_counter("Batteries included. Beautiful is better than ugly.", 256)
MyCounter.global_counter.print()

実行結果

$ python my_main.py
No Name: 0
Current Name=John Doe
Current Value=10
John Doe: 10
Batteries included. Beautiful is better than ugly.: 256

Ruby

my_counter.rb

class MyCounter
  attr_accessor :name

  @@global_counter = nil

  def initialize(init_value)
    @name = 'No Name'
    @count = init_value
  end

  def add_count(number)
    @count += number
  end

  def print
    puts format
  end

  private

  def format
    "#{@name}: #{@count}"
  end

  def self.create_global_counter(name, init_value)
    @@global_counter = self.new(init_value)
    @@global_counter.name = name
  end

  def self.global_counter
    @@global_counter
  end
end

my_main.rb

require './my_counter'

counter = MyCounter.new(0)
counter.print

counter.name = 'John Doe'
current_name = counter.name
puts "Current Name=#{current_name}"

current_value = counter.add_count(10)
puts "Current Value=#{current_value}"

counter.print

MyCounter.create_global_counter('Everything is an object in Ruby', 256)
MyCounter.global_counter.print

$ ruby my_main.rb
No Name: 0
Current Name=John Doe
Current Value=10
John Doe: 10
Everything is an object in Ruby: 256

はじめに

以前に混乱した[i for i in l]の書き方の名前が「リスト内包表記」であることが判明した。
GPSデータの処理でもリストを使うので、内包表記について慣れていく。

実装

l = [i for i in [1, 2, 3, 4, 5, 6, 7, 8, 9] if i % 2 == 0]
print( l )

結果

[2, 4, 6, 8]

びっくりするほど簡単に書けた。
Cからすると恐ろしいほどだ。
とはいえ、少し長くはなってきた。
これを見ても「簡単」と思えるほどにはPythonに親しんできたのだと思いたい。

クローラーとは

まずクローリングとは、スクレイピングとセットで扱われ、自動的にインターネットを巡回し、
様々なWebサイトからコンテンツを収集・保存していく処理
それを行うソフトウェアをクローラーと呼ぶ
スクレイピング
webページから取得したコンテンツから必要な情報を抜き出したり、整形したりすることを指す
クローリング
ソフトウェアが自動的にWebコンテンツを収集・保存していくことを指す
コンテンツ内のハイパーリンクを次々辿っていったりもできる

Scrapy

役割的にはRequests、BeautifulSoupと同じ
コンテンツのダウンロード、解析、保存をより簡単にできる、かつ多機能なのがScrapy（先に言ってよ...w）
scrapyのプロジェクトを作る必要があるみたい

$ pip install scrapy

いつも通りインストール　なんかちょっと長かった

• Create new Scrapy Project

$ scrapy startproject sample

なるほど〜

Spider

• spiderというクラスを定義してscrapyコマンドを実行するだけでクローリングしてくれる（ずいぶん簡単に言いますな）

sample/spiders/sample_spider.py

import scrapy

class SampleSpider(scrapy.Spider):
    name = 'sample' # このクローラーの名前
    allowed_domains = ['google.com'] # クローリングを許可するドメイン
    start_urls = ['https://www.google.com/?hl=ja'] # 起点となるURL（もちろん許可されたドメインでないとNG）

    def parse(self, response):
        print(response.text)

上記のようにいくつか設定値があるみたい
starts_urls・・・ 順不同で実行される
parseメソッド ・・・ コンテンツ取得後のコールバックとして実行される、つまりHTML取ってきたあとはここが動く

HTMLを取得してみる

あとは実行するだけでHTMLを取得して表示してくれるとのこと、、、ほんまかいな、、前より簡単すぎぃ
（どうでもいいけど allow を ｱﾛｳ って発音されると気になって一瞬思考停止する）

$ scrapy crawl sample
...
2019-05-29 14:40:16 [scrapy.core.engine] INFO: Spider opened
2019-05-29 14:40:16 [scrapy.extensions.logstats] INFO: Crawled 0 pages (at 0 pages/min), scraped 0 items (at 0 items/min)
2019-05-29 14:40:16 [scrapy.extensions.telnet] INFO: Telnet console listening on 127.0.0.1:6023
2019-05-29 14:40:16 [scrapy.core.engine] DEBUG: Crawled (200) <GET https://www.google.com/robots.txt> (referer: None)
2019-05-29 14:40:16 [scrapy.downloadermiddlewares.robotstxt] DEBUG: Forbidden by robots.txt: <GET https://www.google.com/?hl=ja>
2019-05-29 14:40:16 [scrapy.core.engine] INFO: Closing spider (finished)
...

おー？なんかSpiderがopenしてcloseしてるし、crawledでURLが表示されてる！！
でも Forbidden が出てる...　うまくはいってないっぽい
robots.txtに設定が足りてないってこと？robots.txtないけど。。。httpsも許可する設定とかがあるのかしら。
今は置いておいて

sample/spiders/sample_spider.py

# （略）

allowed_domains = ['quotes.toscrape.com']
start_urls = ['http://quotes.toscrape.com/page/1/'] # Scrapy公式チュートリアルで使われてたサンプルURL

# （略）

こう変えたらゴチャゴチャ返ってくる中にHTMLが！！

...
2019-05-29 20:26:02 [scrapy.core.engine] INFO: Spider opened
2019-05-29 20:26:02 [scrapy.extensions.logstats] INFO: Crawled 0 pages (at 0 pages/min), scraped 0 items (at 0 items/min)
2019-05-29 20:26:02 [scrapy.extensions.telnet] INFO: Telnet console listening on 127.0.0.1:6023
2019-05-29 20:26:03 [scrapy.core.engine] DEBUG: Crawled (404) <GET http://quotes.toscrape.com/robots.txt> (referer: None)
2019-05-29 20:26:03 [scrapy.core.engine] DEBUG: Crawled (200) <GET http://quotes.toscrape.com/page/1/> (referer: None)
<!DOCTYPE html>
<html lang="en">
<head>
        <meta charset="UTF-8">
        <title>Quotes to Scrape</title>

...

</footer>
</body>
</html>
2019-05-29 20:26:03 [scrapy.core.engine] INFO: Closing spider (finished)
2019-05-29 20:26:03 [scrapy.statscollectors] INFO: Dumping Scrapy stats:
{'downloader/request_bytes': 453,
 'downloader/request_count': 2,
 'downloader/request_method_count/GET': 2,
 'downloader/response_bytes': 2701,
 'downloader/response_count': 2,
 'downloader/response_status_count/200': 1,
 'downloader/response_status_count/404': 1,
 'finish_reason': 'finished',
 'finish_time': datetime.datetime(2019, 5, 29, 11, 26, 3, 683008),
 'log_count/DEBUG': 2,
 'log_count/INFO': 9,
 'memusage/max': 49848320,
 'memusage/startup': 49848320,
 'response_received_count': 2,
 'robotstxt/request_count': 1,
 'robotstxt/response_count': 1,
 'robotstxt/response_status_count/404': 1,
 'scheduler/dequeued': 1,
 'scheduler/dequeued/memory': 1,
 'scheduler/enqueued': 1,
 'scheduler/enqueued/memory': 1,
 'start_time': datetime.datetime(2019, 5, 29, 11, 26, 2, 333924)}
2019-05-29 20:26:03 [scrapy.core.engine] INFO: Spider closed (finished)

いい感じ！！　成功したから長めに結果を載せちゃう

取得したHTMLを解析、保存してみる

より実戦で使えそうな感じにするため、上記の方法で取得したHTMLから必要な情報だけ抜き出して、CSVに出してみる
さっきのチュートリアルのサイトは、名言？とその主がダーっと載っているので、名言＆主CSVを出す
さっきprintしてたコールバックの部分でやりくりすれば良さそうよね

sample/spiders/sample_spider.py

import scrapy

class SampleSpider(scrapy.Spider):
    name = 'sample'
    allowed_domains = ['quotes.toscrape.com']
    start_urls = ['http://quotes.toscrape.com/page/1/']

    def parse(self, response):
        result_list = []
        div_quotes = response.css('div.quote')
        for div_quote in div_quotes:
            span_text = div_quote.css('span.text::text').extract_first() # css同様class指定だと複数取れるため1件目
            quote = span_text.strip()
            author = div_quote.css('small.author::text').extract_first().strip()
            # 属性で抜きたい場合は ::attr(href) とかできるので覚えといたほうがよさそう

            result_list.append({
                'quote': quote,
                'author': author
            })

        print(result_list)

• ポイント１：必要なデータのlistを作ることを目指す　各要素は辞書にしておく
• ポイント２：response.css('~~~')のようにするとCSSセレクターを利用して欲しい情報まで辿れる
  

printの結果（１行なので見やすく整形しました）.txt

[
  {
    'quote': '“The world as we have created it is a process of our thinking. It cannot be changed without changing our thinking.”', 
    'author': 'Albert Einstein'
  }, 
  {
    'quote': '“It is our choices, Harry, that show what we truly are, far more than our abilities.”', 
    'author': 'J.K. Rowling'
  }, 
  ...
]

１発でうまくいってちょっと感動！
さてここまでくればもう簡単、printではなくcsvに出すぞ
printの箇所をreturnに変える

sample/spiders/sample_spider.py

        # （略）
        return result_list

そして実行するコマンドを以下のようにする

$ scrapy crawl sample -o quates_list.csv

おおおおおおででで出てきた！！

完璧すぎる
そう、
辞書のlistにしておくだけで、ファイル種類（今回はcsv）に合わせて良しなに出してくれる
という涙級の優しさ

なのでコマンドを変えるだけでjsonにもなる（泣）

もっと深掘り：リンクを辿ってクローリング

さっきのサイトのAboutリンクの先のお誕生日まで抜いてきちゃおう
はじめにポイント
- 別のページを取得するにはscrapy.Request()を使う
- その別のページ取得後に実行したいcallbackを引数で渡す

sample/spiders/sample_spider.py

    def parse(self, response):
        result_list = []
        div_quotes = response.css('div.quote')
        for div_quote in div_quotes:
            span_text = div_quote.css('span.text::text').extract_first()
            quote = span_text.strip()
            author = div_quote.css('small.author::text').extract_first().strip()
            about_link = div_quote.css('a::attr(href)').extract_first().strip()

            # 初めはこう書いてたけど・・・
            # yield scrapy.Request(f'http://quotes.toscrape.com{about_link}', callback=self.parse_about_page)
            # 公式に倣えばもっとスマートに書けた urljoin() でURL調整
            # yield scrapy.Request(response.urljoin(about_link), callback=self.parse_about_page)
            # もっともっとスマートに　follow() を使えばurljoinすら要らない
            yield response.follow(about_link, callback=self.parse_about_page)

            result_list.append({
                'quote': quote,
                'author': author,
            })

        return result_list

    def parse_about_page(self, response):
        birthday = response.css('span.author-born-date::text').extract_first().strip()

        return { 'birthday':birthday }

一応詳細ページのクローリングには成功したけど、一覧と詳細を組み合わせることはできなかった・・・
（quote, author, birthday のCSVが出したかった）
yield、returnを書いた時点でそこで結果が出力されてしまう
イメージ的にはこういうことがしたかったのにな〜　やりたいことは伝わるはずw

image.py

            # ここでクローリングして解析してきた結果を変数に保持して
            birthday = yield response.follow(about_link, callback=self.parse_about_page)

            # 一緒に出す
            result_list.append({
                'quote': quote,
                'author': author,
                'birthday': birthday
            })

        return result_list

    def parse_about_page(self, response):
        birthday = response.css('span.author-born-date::text').extract_first().strip()

        return { 'birthday':birthday }

一旦別々に出力しておいて、あとでうまいことすれば結果的には同じことできなくはないけど、
もっとスマートにできそう。。。一旦断念

感想

色々やり方、機能があるみたいだけど、基本的なことするだけなら超お手軽
一覧で取得したデータと詳細のデータを組み合わせてってのができそうでできなかったけど、これらを活かして自分用のクローラーとか作ってみたいな〜
うまくいかなかったのはもう少し調べてみる
わかったら更新します

クローラーとは

まずクローリングとは、スクレイピングとセットで扱われ、自動的にインターネットを巡回し、
様々なWebサイトからコンテンツを収集・保存していく処理
それを行うソフトウェアをクローラーと呼ぶ
スクレイピング
webページから取得したコンテンツから必要な情報を抜き出したり、整形したりすることを指す
クローリング
ソフトウェアが自動的にWebコンテンツを収集・保存していくことを指す
コンテンツ内のハイパーリンクを次々辿っていったりもできる

Scrapy

役割的にはRequests、BeautifulSoupと同じ
コンテンツのダウンロード、解析、保存をより簡単にできる、かつ多機能なのがScrapy（先に言ってよ...w）
scrapyのプロジェクトを作る必要があるみたい

$ pip install scrapy

いつも通りインストール　なんかちょっと長かった

• Create new Scrapy Project

$ scrapy startproject sample

なるほど〜

Spider

• spiderというクラスを定義してscrapyコマンドを実行するだけでクローリングしてくれる（ずいぶん簡単に言いますな）

sample/spiders/sample_spider.py

import scrapy

class SampleSpider(scrapy.Spider):
    name = 'sample' # このクローラーの名前
    allowed_domains = ['google.com'] # クローリングを許可するドメイン
    start_urls = ['https://www.google.com/?hl=ja'] # 起点となるURL（もちろん許可されたドメインでないとNG）

    def parse(self, response):
        print(response.text)

上記のようにいくつか設定値があるみたい
starts_urls・・・ 順不同で実行される
parseメソッド ・・・ コンテンツ取得後のコールバックとして実行される、つまりHTML取ってきたあとはここが動く

HTMLを取得してみる

あとは実行するだけでHTMLを取得して表示してくれるとのこと、、、ほんまかいな、、前より簡単すぎぃ
（どうでもいいけど allow を ｱﾛｳ って発音されると気になって一瞬思考停止する）

$ scrapy crawl sample
...
2019-05-29 14:40:16 [scrapy.core.engine] INFO: Spider opened
2019-05-29 14:40:16 [scrapy.extensions.logstats] INFO: Crawled 0 pages (at 0 pages/min), scraped 0 items (at 0 items/min)
2019-05-29 14:40:16 [scrapy.extensions.telnet] INFO: Telnet console listening on 127.0.0.1:6023
2019-05-29 14:40:16 [scrapy.core.engine] DEBUG: Crawled (200) <GET https://www.google.com/robots.txt> (referer: None)
2019-05-29 14:40:16 [scrapy.downloadermiddlewares.robotstxt] DEBUG: Forbidden by robots.txt: <GET https://www.google.com/?hl=ja>
2019-05-29 14:40:16 [scrapy.core.engine] INFO: Closing spider (finished)
...

おー？なんかSpiderがopenしてcloseしてるし、crawledでURLが表示されてる！！
でも Forbidden が出てる...　うまくはいってないっぽい
robots.txtに設定が足りてないってこと？robots.txtないけど。。。httpsも許可する設定とかがあるのかしら。
今は置いておいて

sample/spiders/sample_spider.py

# （略）

allowed_domains = ['quotes.toscrape.com']
start_urls = ['http://quotes.toscrape.com/page/1/'] # Scrapy公式チュートリアルで使われてたサンプルURL

# （略）

こう変えたらゴチャゴチャ返ってくる中にHTMLが！！

...
2019-05-29 20:26:02 [scrapy.core.engine] INFO: Spider opened
2019-05-29 20:26:02 [scrapy.extensions.logstats] INFO: Crawled 0 pages (at 0 pages/min), scraped 0 items (at 0 items/min)
2019-05-29 20:26:02 [scrapy.extensions.telnet] INFO: Telnet console listening on 127.0.0.1:6023
2019-05-29 20:26:03 [scrapy.core.engine] DEBUG: Crawled (404) <GET http://quotes.toscrape.com/robots.txt> (referer: None)
2019-05-29 20:26:03 [scrapy.core.engine] DEBUG: Crawled (200) <GET http://quotes.toscrape.com/page/1/> (referer: None)
<!DOCTYPE html>
<html lang="en">
<head>
        <meta charset="UTF-8">
        <title>Quotes to Scrape</title>

...

</footer>
</body>
</html>
2019-05-29 20:26:03 [scrapy.core.engine] INFO: Closing spider (finished)
2019-05-29 20:26:03 [scrapy.statscollectors] INFO: Dumping Scrapy stats:
{'downloader/request_bytes': 453,
 'downloader/request_count': 2,
 'downloader/request_method_count/GET': 2,
 'downloader/response_bytes': 2701,
 'downloader/response_count': 2,
 'downloader/response_status_count/200': 1,
 'downloader/response_status_count/404': 1,
 'finish_reason': 'finished',
 'finish_time': datetime.datetime(2019, 5, 29, 11, 26, 3, 683008),
 'log_count/DEBUG': 2,
 'log_count/INFO': 9,
 'memusage/max': 49848320,
 'memusage/startup': 49848320,
 'response_received_count': 2,
 'robotstxt/request_count': 1,
 'robotstxt/response_count': 1,
 'robotstxt/response_status_count/404': 1,
 'scheduler/dequeued': 1,
 'scheduler/dequeued/memory': 1,
 'scheduler/enqueued': 1,
 'scheduler/enqueued/memory': 1,
 'start_time': datetime.datetime(2019, 5, 29, 11, 26, 2, 333924)}
2019-05-29 20:26:03 [scrapy.core.engine] INFO: Spider closed (finished)

いい感じ！！　成功したから長めに結果を載せちゃう

取得したHTMLを解析、保存してみる

より実戦で使えそうな感じにするため、上記の方法で取得したHTMLから必要な情報だけ抜き出して、CSVに出してみる
さっきのチュートリアルのサイトは、名言？とその主がダーっと載っているので、名言＆主CSVを出す
さっきprintしてたコールバックの部分でやりくりすれば良さそうよね

sample/spiders/sample_spider.py

import scrapy

class SampleSpider(scrapy.Spider):
    name = 'sample'
    allowed_domains = ['quotes.toscrape.com']
    start_urls = ['http://quotes.toscrape.com/page/1/']

    def parse(self, response):
        result_list = []
        div_quotes = response.css('div.quote')
        for div_quote in div_quotes:
            span_text = div_quote.css('span.text::text').extract_first() # css同様class指定だと複数取れるため1件目
            quote = span_text.strip()
            author = div_quote.css('small.author::text').extract_first().strip()
            # 属性で抜きたい場合は ::attr(href) とかできるので覚えといたほうがよさそう

            result_list.append({
                'quote': quote,
                'author': author
            })

        print(result_list)

• ポイント１：必要なデータのlistを作ることを目指す　各要素は辞書にしておく
• ポイント２：response.css('~~~')のようにするとCSSセレクターを利用して欲しい情報まで辿れる
  

printの結果（１行なので見やすく整形しました）.txt

[
  {
    'quote': '“The world as we have created it is a process of our thinking. It cannot be changed without changing our thinking.”', 
    'author': 'Albert Einstein'
  }, 
  {
    'quote': '“It is our choices, Harry, that show what we truly are, far more than our abilities.”', 
    'author': 'J.K. Rowling'
  }, 
  ...
]

１発でうまくいってちょっと感動！
さてここまでくればもう簡単、printではなくcsvに出すぞ
printの箇所をreturnに変える

sample/spiders/sample_spider.py

        # （略）
        return result_list

そして実行するコマンドを以下のようにする

$ scrapy crawl sample -o quates_list.csv

おおおおおおででで出てきた！！

完璧すぎる
そう、
辞書のlistにしておくだけで、ファイル種類（今回はcsv）に合わせて良しなに出してくれる
という涙級の優しさ

なのでコマンドを変えるだけでjsonにもなる（泣）

もっと深掘り：リンクを辿ってクローリング

さっきのサイトのAboutリンクの先のお誕生日まで抜いてきちゃおう
はじめにポイント
- 別のページを取得するにはscrapy.Request()を使う
- その別のページ取得後に実行したいcallbackを引数で渡す

sample/spiders/sample_spider.py

    def parse(self, response):
        result_list = []
        div_quotes = response.css('div.quote')
        for div_quote in div_quotes:
            span_text = div_quote.css('span.text::text').extract_first()
            quote = span_text.strip()
            author = div_quote.css('small.author::text').extract_first().strip()
            about_link = div_quote.css('a::attr(href)').extract_first().strip()

            # 初めはこう書いてたけど・・・
            # yield scrapy.Request(f'http://quotes.toscrape.com{about_link}', callback=self.parse_about_page)
            # 公式に倣えばもっとスマートに書けた urljoin() でURL調整
            # yield scrapy.Request(response.urljoin(about_link), callback=self.parse_about_page)
            # もっともっとスマートに　follow() を使えばurljoinすら要らない
            yield response.follow(about_link, callback=self.parse_about_page)

            result_list.append({
                'quote': quote,
                'author': author,
            })

        return result_list

    def parse_about_page(self, response):
        birthday = response.css('span.author-born-date::text').extract_first().strip()

        return { 'birthday':birthday }

一応詳細ページのクローリングには成功したけど、一覧と詳細を組み合わせることはできなかった・・・
（quote, author, birthday のCSVが出したかった）
yield、returnを書いた時点でそこで結果が出力されてしまう
イメージ的にはこういうことがしたかったのにな〜　やりたいことは伝わるはずw

image.py

            # ここでクローリングして解析してきた結果を変数に保持して
            birthday = yield response.follow(about_link, callback=self.parse_about_page)

            # 一緒に出す
            result_list.append({
                'quote': quote,
                'author': author,
                'birthday': birthday
            })

        return result_list

    def parse_about_page(self, response):
        birthday = response.css('span.author-born-date::text').extract_first().strip()

        return { 'birthday':birthday }

一旦別々に出力しておいて、あとでうまいことすれば結果的には同じことできなくはないけど、
もっとスマートにできそう。。。一旦断念

感想

色々やり方、機能があるみたいだけど、基本的なことするだけなら超お手軽
一覧で取得したデータと詳細のデータを組み合わせてってのができそうでできなかったけど、これらを活かして自分用のクローラーとか作ってみたいな〜
うまくいかなかったのはもう少し調べてみる
わかったら更新します

TL;DR

$ . your-venv/bin/activate
$ pip install -U boto
$ sudo cat > /etc/boto.cfg <<'^Z'
[s3]
use-sigv4 = True
host = s3.ap-northeast-1.amazonaws.com
^Z

はじめに

間もなく 2019/6/24 に AWS S3 の署名バージョン２が廃止になり、現在署名バージョン２を使って S3 アクセスを行っているシステムに影響が出ます。

私が前任者から受け継がされたシステムも古い Boto2 ライブラリ (Python 2.7) を利用し、アップデートもせずに放置されているので、当然署名バージョン２を利用しています。

このシステムは再デプロイのやり方も不明（というかデプロイされた後にサーバー上でいじられた形跡がある）なため、なんとか動作したままライブで署名バージョン４に対応したいところです。

なぜか Boto2 を署名バージョン４に対応させる情報がググってもあまり出てこず、ややはまったのでここに書き記しておきます。

Boto2 と Boto3

この記事で説明する方法は Boto2 が対象です。Boto3 には当てはまらないので注意してください。

署名バージョン２を使っていることを確認する

実際にこのシステムが署名バージョン２を使っていることを確かめるために「【Amazon S3を利用しているすべての人が安心して2019/06/24を迎えるために】CloudTrailとAthenaを用いたS3の署名バージョンの確認方法 ｜ DevelopersIO」を参考に署名バージョン２を利用した S3 アクセスを検出しました。

上の記事の最後の Athena に対するクエリを実行すると CSV がダウンロードできるようになります。この CSV（巨大です）を data.csv として、以下のようにすると簡単にインスタンスが特定できます。

$ tail -n +2 data.csv | sort -t, -k1r,1 | sort -sut, -k6,6 -k5,5 > data.uniq.csv
$ tail -n +2 data.csv | sort -t, -k1r,1 | sort -st, -k6,6 -k5,5 > data.sorted.csv

data.uniq.csv は (S3 バケット名, インスタンスの IP アドレス）が同一のもののうち、最も最近にあったアクセスだけが抽出されます。data.sorted.csv は同一のものを時間順にソートしただけのものです。data.uniq.csv は小さいので見やすいです。詳細（バケット名に加えてキーを見たい場合など）を見たい場合は data.sorted.csv を参照します。

この方法で、修正対象のシステムのインスタンスの IP から確かに署名バージョン２でアクセスがされており、そのユーザーエージェントが [Boto/2.25.0 Python/2.7.4 Linux/3.4.73-64.112.amzn1.x86_64] であることが確認できました。Boto2 の 2.25.0 が使われているのが観察できますね。

修正方法を調べる

AWS SDK の種類ごとの修正方法が「署名バージョン 2 から署名バージョン 4 への移行」に掲載されています。これによると：

SDK	Upgrade	コード変更が必要
Boto2	Boto2 v2.49.0 にアップグレード	はい

う～む「コード変更が必要」は困る...

さらに「リクエスト認証での署名バージョンの指定」を見ると：

boto デフォルト設定ファイルに以下を指定します。

[s3] use-sigv4 = True

とあります。コード変更といっても設定ファイルの変更だけで済みそうですが、この boto デフォルト設定ファイルってどこにあるの？しかも何この変なフォーマットは

[s3]
use-sigv4 = True

じゃないの？

Boto2 の設定ファイル

Boto2 の設定ファイルは「Boto Config — boto v2.49.0」に書いてあります。 /etc/boto.cfg とか ~/.boto に書けばいいみたい。

あと、やはりフォーマットは ini フォーマットで、[s3] の後に改行が必要みたいです。

修正作業

材料が揃ったので修正しましょう。今回修正対象のシステムは cron で動作しているので、まずそれを止めます。virtualenv を利用して動いている（よかった...）ので virtualenv 内の boto をアップグレードします。

$ cp -a your-venv your-venv-backup  # ← 非常時ロールバック用
$ . your-venv/bin/activate
$ pip install -U boto
$ pip show boto
---
Name: boto
Version: 2.49.0
Location: /home/user/your-venv/lib/python2.7/site-packages
Requires: 

Boto2 は依存先が無く他のパッケージを巻き添えにせずにアップグレードできてレガシーシステムにやさしいです。

次に Boto2 設定ファイルを変更します。今回は既存の設定ファイルは無さそうだったので新規で作成します。

$ sudo cat > /etc/boto.cfg <<'^Z'
[s3]
use-sigv4 = True
^Z

さてテスト実行してみると... S3 に全くアクセスできず！！！... 困った。

さらなる修正

原因をググってもなかなか探し当てられませんでしたが、それらしきものを発見: python - Using boto for AWS S3 Buckets for Signature V4 - Stack Overflow

この記事には API エンドポイントに s3.amazonaws.com ではなく s3.eu-central-1.amazonaws.com というリージョン毎のエンドポイントを指定しなくてはならない、と書いてあります。う～ん、これ boto.cfg に書けるのか？

まず東京リージョンの S3 API エンドポイントを調べます: AWS のリージョンとエンドポイント - アマゾン ウェブ サービス によると s3.ap-northeast-1.amazonaws.com でした。ということで...

$ sudo cat > /etc/boto.cfg <<'^Z'
[s3]
use-sigv4 = True
host = s3.ap-northeast-1.amazonaws.com
^Z

テスト実行は... うまくいきました。cron を元通りにして終了。

この [s3] セクションの設定の公式ドキュメントは一体どこにあるのだろうか？見つからない。

OneMax問題

OneMax問題とは, [0, 1, 1, 1, 0, 1, 0, 1, 1, 0]のように0/1からなる数列の和を最大化する問題です. 和が最大となるには, すべての要素が1となることを目指します.
今回は遺伝的アルゴリズムを用いて, OneMax問題を解いてみます.

遺伝的アルゴリズム

遺伝的アルゴリズムは, 生物の進化の過程を真似て作られたアルゴリズムで確率的探索手法の1つです. 交叉や突然変異という操作を繰り返しながら, 個体を進化させていきます.
基本的な流れは以下のようになります.

選択

次世代個体群を生成するために, 親個体を選択します.
親個体の選択方法には以下のようなものがあります.
-ルーレット方式
-ランキング方式
-エリート方式

交叉

選択した親個体を交叉させて, 子個体を生成します.
交叉方法には以下のようなものがあります.
-一様交叉
-一点交叉
-多点交叉

突然変異

生成した個体に突然変異を加えます.

実践

それでは実際にGAでOneMax問題を解いてみます.
実験の設定は
遺伝子長:100
世代数:100
選択:トーナメント選択(トーナメントサイズ:2)
交叉:一様交叉
突然変異確率:1/遺伝子長
親個体数:100
子個体数:50
個体更新:親個体の優良個体50+子個体

# -*- coding: utf-8 -*-
import random
import copy

# 各種パラメータ
length = 100  # 遺伝子長
population = 100  # 親個体数
offspring_n = 50    # 子個体数
generation = 100     # 世代数
mutation_rate = 1.0/100.0    # 突然変異確率(1/遺伝子長)

# 初期個体群生成
def initialize():
    gene = [[random.randint(0,1) for i in range(length)] for j in range(population)]
    return gene

# 個体評価
def evaluate(gene):
    fitness = []
    for i in range(population):
        fitness.append(sum(gene[i])/length)
    return fitness

# minを見つける
def find_min(fitness):
    min = 10
    for i in range(population):
        if fitness[i] < min:
            min = fitness[i]
    return min

# maxを見つける
def find_max(fitness):
    max = 0
    for i in range(population):
        if fitness[i] > max:
            max = fitness[i]
    return max

# 親選択
def choice_parents(gene, fitness):
    father_index = random.randint(0,population-1)
    mother_index = random.randint(0,population-1)
    if fitness[father_index] > fitness[mother_index]:
        parents = gene[father_index]
    else:
        parents = gene[mother_index]
    return parents

# 交叉
def crossover(father, mother):
    offspring = []
    for i in range(length):
        p = random.random()
        if p < 0.5:
            offspring.append(father[i])
        else:
            offspring.append(mother[i])
    return offspring

# 突然変異
def mutation(offspring):
    for i in range(length):
        p = random.random()
        if p < mutation_rate:
            if offspring[i] == 0:
                offspring[i] = 1
            else:
                offspring[i] = 0
    return offspring

# 個体更新(親個体エリート50 + 子個体50 = 計100)
def elite(gene, fitness, next_gene):
    sort_fitness = sorted(fitness, reverse=True)
    gen_tmp = []
    for i in range(int(population/2)):
        index = fitness.index(sort_fitness[i])
        gen_tmp.append(gene[index])
    gen_tmp.extend(next_gene)
    return gen_tmp

def main():
    generation_count = 1
    next_gene = []

    # 初期個体群生成
    gene = initialize()

    while generation_count <= generation:
        next_gene.clear()
        #print(gene)
        # 個体評価
        fitness = evaluate(gene)
        min = find_min(fitness)
        max = find_max(fitness)
        print("generation:", generation_count)
        print("min", min)
        print("max", max)
        if min == 100:
            break
        # 次世代個体群生成
        for i in range(offspring_n):
            # 個体選択
            father = choice_parents(gene, fitness)
            mother = choice_parents(gene, fitness)
            # 交叉
            offspring = crossover(father, mother)
            # 突然変異
            offspring = mutation(offspring)
            next_gene.append(offspring)
        # 個体更新
        gene = elite(gene, fitness, next_gene)
        generation_count += 1
        print("--------------------------------")

if __name__ == "__main__":
    main()

結果は以下のようになりました.
最初の世代ではmin, maxともに低い値となっています.
それが最終世代ではmaxは最大値に達し, minも大きな値を獲得できています.
遺伝的アルゴリズムもより, 個体が進化していることを確認できました.

generation: 1
min 0.37
max 0.62
--------------------------------
generation: 2
min 0.39
max 0.63
--------------------------------

.
.
.

generation: 99
min 0.96
max 1.0
--------------------------------
generation: 100
min 0.95
max 1.0
--------------------------------

はじめに

こんにちは。
今回は「追いかけっこゲーム」を作っていきたいと思います。
完成するとこんな感じに動きます！

赤い円が自分(player)で青い円が敵(enemy)です。
マウスカーソルを操作して、赤い円を動かすことができます。青い円に触れたら止まります。つまりゲームオーバーといった感じです！

環境

• Windows 10 home
• Python 3.7.1

numpyを使うので以下のコマンドでインストールしておきましょう。

pip install numpy

「追いかけっこゲーム」の制作

インポート

使うライブラリは以下の2つです。

import tkinter as tk
import numpy as np

ウィンドウの作成

import tkinter as tk
import numpy as np

class Application(tk.Frame):
    def __init__(self,master):
        super().__init__(master)
        self.pack()

        self.width=self.height=500
        master.geometry(str(self.width)+"x"+str(self.height))
        master.title("追いかけっこゲーム")


def main():
    win = tk.Tk()
    app = Application(master = win)
    app.mainloop()


if __name__ == "__main__":

    main()

このプログラムでウィンドウを作っていきます。
サイズは500×500です。タイトル名は追いかけっこゲームです。
詳しくはPythonによるTkinterを使った雛形（クラス化手法）の記事をご参照ください。

敵の動きのベクトル取得に関して

敵の動きのベクトルは以下のenemyVec関数で取得しています。少しややこしい部分なので、解説を入れようと思います。

def enemyVec(self,player,enemy,speed):#敵の動き(x,y)のベクトルを返す
    rad=np.arctan((player.y-enemy.y)/(player.x-enemy.x))#向き（角度の計算）
    if player.x-enemy.x >= 0:
        vx=np.cos(rad)*speed
        vy=np.sin(rad)*speed
    else:
        vx=np.cos(rad)*(-1*speed)
        vy=np.sin(rad)*(-1*speed)
    return [vx,vy]

引数として、playerオブジェクトとenemyオブジェクトを渡しています。それぞれの中心座標を計算で使うためです。ここで、この2つのオブジェクトの向きの計算が少しややこしくなると思います。つまり、2点間の角度の計算です。
2点の座標が分かっているので、ベクトルのx成分とy成分を求めることができます。

(例) A(x1,y1),B(x2,y2)の場合

\vec{AB} =
\begin{matrix}
x2 - x1 \\
y2 - y1 
\end{matrix}

そして、x軸とベクトルがなす角度θを以下のtanの逆関数を使って表すことができます。

θ = \tan^-1(\frac{y}{x})

tanθの範囲が-π/2～π/2なので、playerとenemyのx座標の差の正負で場合分けをし、speedに-1をかけて、方向を変えています。

完成したプログラム

こちらが完成したプログラムになります。

追いかけっこゲーム.py

import tkinter as tk
import numpy as np

class Circle():
    def __init__(self,canvas,x,y,r,color,tag=None):
        self.canvas = canvas
        self.x = x #円の中心のx座標
        self.y = y #円の中心のy座標
        self.r = r #円の半径
        self.color = color
        self.tag = tag

    def createCircle(self):
        self.canvas.create_oval(self.x-self.r,self.y-self.r,self.x+self.r,self.y+self.r,fill=self.color,tag=self.tag)

class Application(tk.Frame):
    def __init__(self,master):
        super().__init__(master)
        self.pack()

        self.width=self.height=500
        master.geometry(str(self.width)+"x"+str(self.height))
        master.title("追いかけっこゲーム")

        self.canvas = tk.Canvas(master,width=self.width,height=self.height,bg="black")
        self.canvas.pack()

        self.player = Circle(self.canvas,250,250,30,"red","player") #インスタンスplayerの生成
        self.enemy = Circle(self.canvas,0,0,30,"blue","enemy") #インスタンスenemyの生成

        self.canvas.bind("<Motion>",self.mouseEvent)
        self.master.after(50,self.update)

    def update(self):
        if self.judgeflag(self.player,self.enemy):
            eV = self.enemyVec(self.player,self.enemy,10)
            self.enemy.x += eV[0]
            self.enemy.y += eV[1]
            self.canvas.delete("player")
            self.canvas.delete("enemy")
            self.player.createCircle()
            self.enemy.createCircle()

        self.master.after(50,self.update)

    def enemyVec(self,player,enemy,speed): #敵の動き(x,y)のベクトルを返す
        rad=np.arctan((player.y-enemy.y)/(player.x-enemy.x)) #向き（角度の計算）
        if player.x-enemy.x >= 0:
            vx=np.cos(rad)*speed
            vy=np.sin(rad)*speed
        else:
            vx=np.cos(rad)*(-1*speed)
            vy=np.sin(rad)*(-1*speed)
        return [vx,vy]

    def judgeflag(self,player,enemy): #当たり判定
        if np.sqrt((player.x-enemy.x)**2+(player.y-enemy.y)**2) > player.r+enemy.r: return True
        else: return False

    def mouseEvent(self,event):
        self.player.x = event.x
        self.player.y = event.y

def main():
    win = tk.Tk()
    app = Application(master=win)
    app.mainloop()

if __name__ == "__main__":
    main()

以上で追いかけっこゲームは完成になります。
ここまで読んでいただき、ありがとうございました。

はじめに

今ではお馴染みとなったGradient Boostingですが、きちんと説明を読んだことがなかったので調べてみました。

Boostingというと色々ありますが、今回まとめたのはFriedmanのGradient Boostingです。

参考

• https://en.wikipedia.org/wiki/Gradient_boosting
  
• Gradientboostingmachines,atutorial(PDF)

手順

定義

説明変数を$\mathbf{x}_i(i=1,...,N)$とし、目的変数を$t_i (i=1,...,N)$とします。説明変数の次元は$d$とします。

予測器を$f$と書き、予測値を$y$で表します。我々が求めたいものは、損失関数$\Psi$を最小にする$\hat f$です。

\begin{align}
y_i &= \hat f(\mathbf{x_i}) \\
\hat f &= \arg \min_{f \in \mathcal{F}} \sum_{i=1}^N \Psi(t_i, f(\mathbf{x_i}))
\end{align}

ここで$\mathcal{F}$は関数$f$の空間です。

関数推定の問題を扱いやすくするために、関数$f$を関数空間がパラメータ$\theta$で張られる族$f(\mathbf{x}, \theta)$に制限します。すると上の式は次のように書き直せます。

\begin{align}
y_i &= \hat f(\mathbf{x_i}) \\
\hat f(\mathbf{x}) &= f(\mathbf{x}, \hat\theta) \\
\hat \theta &= \arg \min_{\theta \in D} \sum_{i=1}^N \Psi(t_i, f(\mathbf{x_i}, \theta))
\end{align}

勾配降下法

$\theta$を導入したことで全データの損失関数を次のように表せます。

J(\theta) = \sum_{i=1}^N \Psi(t_i, f(\mathbf{x_i}, \theta))

$J$が最小になる$\theta$が$\hat \theta$です。解析的には$J$のgradを計算し、大域的に極小となる点を求めればよいです。しかし、損失関数の形によっては解析的に解を求めるのが難しいこともあります。そこで、$J$の微分を計算し、$J$が小さくなる方向に$\theta$を変化させる、ということを逐次的に行います。

\theta_{\mathrm{new}} =\theta - \eta \nabla_\theta J(\theta)

上記の式では更新ごとにすべての訓練用データが含まれていますが、データ一つごとに$\theta$を更新すればオンライン学習に、少量ごとに更新すればミニバッチ法になります。ちなみに、訓練データをすべて使い切るのに必要な回数をイテレーション数と呼びます。

\mathrm{イテレーション数} = \mathrm{データ数} / \mathrm{バッチサイズ}

Gradient Boosting

前章での$\theta$の更新ですが、勾配法との比較のために次のように書き直してみましょう。

イテレーションを$M$回行うとし、各イテレーションを$m=1,...,M$でラベルします。更新の式を次のように書き換えます。

\begin{align}
\theta^{(m)} &= \theta^{(m-1)} + \delta\theta^{(m)}\\
\delta \theta^{(m)} &= - \eta \nabla_\theta J(\theta^{(m-1)})
\end{align}

ここで$\theta^{(0)}$は初期値で、適当に与えます。すると$M$回のイテレーションで$\theta$が次のように更新されることが分かります。

\theta^{(M)} = \theta^{(0)} + \sum_{m=1}^{M}\delta\theta^{(m)}

式からわかるように、初期パラメータ$\theta^{(0)}$に修正$\delta\theta^{(m)}$を加えていくことで最適な$\hat\theta$を求める仕組みになっています。

そこで、この発想を関数空間に適用することを考えます。すなわちパラメータを修正するのではなく、初期予測器に対して修正予測器を加えていく仕組みを考えます。

f^{(M)} = f^{(0)} + \sum_{m=1}^{M}\delta f^{(m)}

勾配法とのアナロジーで更新は次のようにすれば良いとわかります。

\begin{align}
f^{(m)}(\mathbf{x}) &= f^{(m-1)}(\mathbf{x}) + \delta f^{(m)}(\mathbf{x})\\
\delta f^{(m)} &= -\rho \nabla_f J(f^{(m-1)})
\end{align}

上の式は素朴なアナロジーでしたが、関数空間の場合はこの更新をブラッシュアップできます。

どういうことかというと、$\theta$の空間の場合、どれくらい$\theta$をシフトさせれば良いかわかりませんでしたので、学習率$\eta$を手動でチューニングする必要がありました。ところが、関数空間の場合は手元に予測器があるので、予測がより良くなるように、すなわち損失関数が小さくなるように関数のシフトを決めることができます。

\begin{align}
f^{(m)}(\mathbf{x}) &= f^{(m-1)}(\mathbf{x}) - \rho^{(m)}\nabla_f J\left(f^{(m-1)}(\mathbf{x})\right) \\
\rho^{(m)} &=  \arg \min_{\rho} \sum_{i=1}^N \Psi\left(t_i, f^{(m-1)}(\mathbf{x_i}) - \rho \nabla_f J\left(f^{(m-1)}(\mathbf{x_i})\right)\right)
\end{align}

この発想で考えれば、もはやシフトに使用する関数も自由に選べます。したがってより一般には次のように書けます。

\begin{align}
f^{(m)}(\mathbf{x}) &= f^{(m-1)}(\mathbf{x}) +g^{(m)}(\mathbf{x}) \\
g^{(m)} &=  \arg \min_{g\in \mathcal{G}} \sum_{i=1}^N \Psi\left(t_i, f^{(m-1)}(\mathbf{x_i}) + g(\mathbf{x_i})\right)
\end{align}

このようにしてGradient Boostingでは、ひとつ前のイテレーションの予測器に新たな予測器を加えていくことで、経験損失を小さくしていきます。

Friedmanの論文の説明について補足

Friedmanの論文は実装のことも考えて書かれており、説明が少し異なるので、これについて補足を残します。

前章の最後で、関数のシフトに任意の$g$を使用しました。Friedmanの論文では、この$g$（と$f^{(0)}$）にパラメータ$\theta$を持つ弱学習器$h(\mathbf{x}, \theta)$を考えます。すなわち、前章の最後の式は次のようになります。

\begin{align}
f^{(m)}(\mathbf{x}) &= f^{(m-1)}(\mathbf{x}) + \beta^{(m)} h(\mathbf{x}, \theta^{(m)}) \\
\left(\beta^{(m)}, \theta^{(m)}\right) &=  \arg \min_{\beta, \theta} \sum_{i=1}^N \Psi\left(t_i, f^{(m-1)}(\mathbf{x_i}) + \beta h(\mathbf{x}_i, \theta)\right)
\end{align}

$g$をこのように置くのは、任意の損失関数に対応するためです。損失関数が複雑な場合、$\beta^{(m)}, \theta^{(m)}$を求めるのは大変です。そこで、gradient boosting近似を行います。

まず、損失関数の勾配に予測器をフィットさせます。

\begin{align}
\theta^{(m)} &= \arg \min_{\rho, \theta}\left[ \sum_{i=1}^N\left(\tilde y_{i}^{(m)} - \rho h(\mathbf{x}_i, \theta) \right)^2\right] \\
\tilde y_{i}^{(m)} &= -\left[\frac{\partial\Psi(t_i, f)}{\partial f}\right]_{f=f^{(m-1)}, \mathbf{x}=\mathbf{x}_i}
\end{align}

次に、この弱学習機の係数を最適化します。

\beta^{(m)} =  \arg \min_{\beta} \sum_{i=1}^N \Psi\left(t_i, f^{(m-1)}(\mathbf{x_i}) + \beta h(\mathbf{x}_i, \theta^{(m)})\right)

この手順によって、最適化問題を最二乗問題と損失関数$\Psi$の1変数最適化問題に分解しました。

予測器の具体例として、弱学習器である決定木を考えます。この決定木は葉を$L$個持つとします。また、この決定木は回帰木です。イテレーション$m$の決定木により、説明変数$\mathbf{x}$の空間が$L$個に分割されます。この分割された領域を$R_m^l (m=1,...,M;l=1,...,L)$でラベルします。gradient boosting近似によって、決定木は次の予測値を出力するように学習します。

\begin{align}
h(\mathbf{x}, \{R_m^l\}_{l=1}^L) &= \sum_{l=1}^L \bar y_l^{(m)} 1(\mathbf{x} \in R_m^l)\\
\bar y_l^{(m)} &= \frac{1}{\#(\mathbf{x}_i \in R_m^l)}\sum_{\mathbf{x}_i \in R_m^l}\tilde y_{i}^{(m)}
\end{align}

あとは$\beta$を決めるのですが、領域ごとに$\beta$を最適化できます。予測器が出力する値は$R_m^l$の各領域内で定数なので、$\beta$との積にしてしまって、これを$\gamma$とすると次のように書けます。

\gamma_{l}^{(m)} = \arg \min_{\gamma} \sum_{\mathbf{x}_i\in R_m^l} \Psi\left(t_i, f^{(m-1)}(\mathbf{x_i}) + \gamma\right)

これによって学習器が更新されます。

f^{(m)}(\mathbf{x}) = f^{(m-1)}(\mathbf{x}) +  \gamma_{l}^{(m)}1(\mathbf{x} \in R_m^l)

オーバーフィッティングを避けるために、学習率$\nu(0<\nu<1)$を導入して次のようにする場合もあります。

f^{(m)}(\mathbf{x}) = f^{(m-1)}(\mathbf{x}) +  \nu \gamma_{l}^{(m)}1(\mathbf{x} \in R_m^l)

これをShrinkageと呼びます。Friedmanによると$\nu<0.1$とすると汎化性能が上がるそうです。

おわりに

Boostingアルゴリズムは残差を予測する予測器を加えていくというアルゴリズムですが、他にもいくつか手法があり、それぞれアルゴリズムが結構違うので調べてみると面白いです。

• XGBoost -> XGBoostの概要
• AdaBoost -> AdaBoostについて調べたのでまとめる

買ったきっかけ

Django REST Framework（DRF）を使った開発を行うことになり、
ネット等で情報収集やサンプルプログラムの作成をしていたのですが、
なかなか欲しい情報がなく、まとまった知識を得る為に購入しました。
（ちなみにDjangoはいくつか本があるのですが、DRFはこの本以外にたぶんない）
現場で使える Django REST Framework の薄い本【紙の本】（技術書典6バージョン） - - BOOTH
で買いました。

内容と感想

第1章： Django REST Framework（DRF）の概要

REST APIの基本から、DRFの全体構成などの説明など
図があってわかりやすい。

第2章： モデル・シリアライザ・ビュー・URLconf の基本的な使い方

使い方や処理の流れが記載されてます。
他サイトのチュートリアルは、ModelViewSetを使ったものばかりで、
少し複雑なAPIを作ろうとする時にどうすればいいかわからなかったのですが、
ここでは、APIViewにも言及されていて、こちらの方が汎用的であることがわかりました。
開発ではAPIViewを使うことが多くなりそうです。

第3章： DRF とセキュリティ

基本的なセキュリティーに関する知識で、DRF特有の情報が少ないです。

第4章： DRF と認証

Basic認証、Cokkie認証、トークン認証、リモートユーザー認証、JWT認証の各方式に関しての紹介とどれを使うべきが記載されています。
今回使う予定だった認証方式についても詳しく紹介されていて参考になりました。

第5章： チュートリアル その1、第6章チュートリアル その2

よくある本では、前半はチュートリアルで、後半に詳細な機能やリファレンスが記載されていることが多いのですが、この本は逆でした。
書籍にもソースコードが記載されていますが、全部ではないので、サンプルコードをダウンロードした上で読み進めるとわかりやすいです。
チュートリアル その1は比較的簡単にできますが、チュートリアル その2はVueの知識が少し必要です。
（dockerを使ったローカルの動作環境が特殊でその2はうまく動かなかった・・・。）
今回、APIのみ実装予定だったので、その1を作り、JWT認証に変えたものを試しました。
JWT認証は第6章よりも、第4章の方が参考になります。

第7章： 現場で使える Tips 集

今後、実際に開発する中で適宜読む予定。

全体を通して

DRFは、Djangoの知識が前提になるので、Djangoをうろ覚えだったらまずDjangoを覚えた方がいいです。
ネットのチュートリアルの方が環境の作り方からわかりやすく載ってたりするので、この書籍が難しければそちらを試してから読む方がいいです。
全体通してとてもわかりやすく詳しく書かれていて、最初からこれを読んだらよかったと思いました。

今後は

書籍にも書かれていましたが、
[Django REST Framework] Serializer の 使い方 をまとめてみた
[Django REST Framework] View の使い方をまとめてみた
を見ていく予定です。

概要

前回の記事で取り上げた、scikit-learnのCountVectorizerのanalyzerキーワード引数について、変えると具体的にどう結果が変化するんだろう？
と思ったので一式コードを書いて確認。

環境

• Windows10 Pro
• Python 3.7.3
• scikit-learn 0.21.2

コード

一通りのanalyzerキーワード引数を使ったインスタンスを作成する

from sklearn.feature_extraction.text import CountVectorizer

# 単語でN-gram
cv_word = CountVectorizer(analyzer = 'word')
# 文字でN-gram
cv_char = CountVectorizer(analyzer = 'char')
# 'char_wb'でN-gram
cv_wb = CountVectorizer(analyzer = 'char_wb')
# 呼び出し可能オブジェクト
cv_callable = CountVectorizer(analyzer = lambda x: [x])

対象にする単語と文章の配列を作成する

即興で書いたので、英文の不味さには目をつぶってください……

words = ['Tokyo', 'Osaka', 'Sendai', 'Fukuoka', 'Sapporo', 'Kouchi']
texts = ['Tokyo is the biggest city in japan.', 
    'Osaka is the secondary biggest city in japan.', 
    'Fukuoka is most nearly city from Koria, China in japan.',
    'Sendai called we ''Mori no Miyako''.',
    'Sapporo is most nearly city from Russia in japan.',
    'Kouchi is a japanese city, in japanese region in ''Shikoku''.']

単語の配列で、fit_transformでの特徴抽出と特徴を出力

print("CountVectorizer: target is words")

X = cv_word.fit_transform(words)
print("CountVectorizer: analyzer is word")
print("results array:")
print(str(X.toarray()))
print("feature_names:")
print(str(cv_word.get_feature_names()))

X = cv_char.fit_transform(words)
print("CountVectorizer: analyzer is char")
print("results array:")
print(str(X.toarray()))
print("feature_names:")
print(str(cv_char.get_feature_names()))

X = cv_wb.fit_transform(words)
print("CountVectorizer: analyzer is char_wb")
print("results array:")
print(str(X.toarray()))
print("feature_names:")
print(str(cv_wb.get_feature_names()))

X = cv_callable.fit_transform(words)
print("CountVectorizer: analyzer is callable object")
print("results array:")
print(str(X.toarray()))
print("feature_names:")
print(str(cv_callable.get_feature_names()))

単語で試した出力結果

CountVectorizer: target is words
CountVectorizer: analyzer is word
results array:
[[0 0 0 0 0 1]
 [0 0 1 0 0 0]
 [0 0 0 0 1 0]
 [1 0 0 0 0 0]
 [0 0 0 1 0 0]
 [0 1 0 0 0 0]]
feature_names:
['fukuoka', 'kouchi', 'osaka', 'sapporo', 'sendai', 'tokyo']
CountVectorizer: analyzer is char
results array:
[[0 0 0 0 0 0 0 1 0 2 0 0 0 1 0 1]
 [2 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0]
 [1 0 1 1 0 0 1 0 1 0 0 0 1 0 0 0]
 [1 0 0 0 1 0 0 2 0 1 0 0 0 0 2 0]
 [1 0 0 0 0 0 0 0 0 2 2 1 1 0 0 0]
 [0 1 0 0 0 1 1 1 0 1 0 0 0 0 1 0]]
feature_names:
['a', 'c', 'd', 'e', 'f', 'h', 'i', 'k', 'n', 'o', 'p', 'r', 's', 't', 'u', 'y']
CountVectorizer: analyzer is char_wb
results array:
[[2 0 0 0 0 0 0 0 1 0 2 0 0 0 1 0 1]
 [2 2 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0]
 [2 1 0 1 1 0 0 1 0 1 0 0 0 1 0 0 0]
 [2 1 0 0 0 1 0 0 2 0 1 0 0 0 0 2 0]
 [2 1 0 0 0 0 0 0 0 0 2 2 1 1 0 0 0]
 [2 0 1 0 0 0 1 1 1 0 1 0 0 0 0 1 0]]
feature_names:
[' ', 'a', 'c', 'd', 'e', 'f', 'h', 'i', 'k', 'n', 'o', 'p', 'r', 's', 't', 'u', 'y']
CountVectorizer: analyzer is callable object
results array:
[[0 0 0 0 0 1]
 [0 0 1 0 0 0]
 [0 0 0 0 1 0]
 [1 0 0 0 0 0]
 [0 0 0 1 0 0]
 [0 1 0 0 0 0]]
feature_names:
['Fukuoka', 'Kouchi', 'Osaka', 'Sapporo', 'Sendai', 'Tokyo']

結果から見えるもの

• fit_transformした結果の配列をDataFrameの値部分と考えると、Indexを元の配列、Columnを特徴（feature_name）とすることが出来る。
  例として、一番最初の結果を表とすると以下のようになる。

Index	fukuoka	kouchi	osaka	sapporo	sendai	tokyo
Tokyo	0	0	0	0	0	1
Osaka	0	0	1	0	0	0
Sendai	0	0	0	0	1	0
Fukuoka	1	0	0	0	0	0
Sapporo	0	0	0	1	0	0
Kouchi	0	1	0	0	0	0

• 文字列で指定した場合は特徴（feature_name）は基本的にlower case
  • 公式を見ると、キーワード引数としてlowercaseが存在し、デフォルト値がTrueとなっているので、ここの設定を変えることで変わるかも
• 特徴（feature_name）は文字の昇順で格納される
• analyzerにchar_wbを指定した場合は、特徴（feature_name）に' 'が含まれる
  • おそらくは単語境界を' 'としている、2回カウントされるのは、前後を1回ずつ単語境界とカウントしているから？
• カテゴリ変数のエンコーディングに使う場合、lower caseな単語のみが格納されていて、単純にOne-hotエンコードの結果を取得するだけだけなら、analyzer = 'word'でいけるのでは？

文章の配列で、fit_transformでの特徴抽出と特徴を出力

print("CountVectorizer: target is texts")

X = cv_word.fit_transform(texts)
print("CountVectorizer: analyzer is word")
print("results array:")
print(str(X.toarray()))
print("feature_names:")
print(str(cv_word.get_feature_names()))

X = cv_char.fit_transform(texts)
print("CountVectorizer: analyzer is char")
print("results array:")
print(str(X.toarray()))
print("feature_names:")
print(str(cv_char.get_feature_names()))

X = cv_wb.fit_transform(texts)
print("CountVectorizer: analyzer is char_wb")
print("results array:")
print(str(X.toarray()))
print("feature_names:")
print(str(cv_wb.get_feature_names()))

X = cv_callable.fit_transform(texts)
print("CountVectorizer: analyzer is callable object")
print("results array:")
print(str(X.toarray()))
print("feature_names:")
print(str(cv_callable.get_feature_names()))

文章で試した出力結果

CountVectorizer: target is texts
CountVectorizer: analyzer is word
results array:
[[1 0 0 1 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0]
 [1 0 0 1 0 0 1 1 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0]
 [0 0 1 1 1 1 1 1 1 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0]
 [0 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 0 1 0 0 0 1]
 [0 0 0 1 1 0 1 1 1 0 0 0 0 0 1 1 0 0 0 1 1 0 0 0 0 0 0]
 [0 0 0 1 0 0 2 1 0 2 0 1 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0]]
feature_names:
['biggest', 'called', 'china', 'city', 'from', 'fukuoka', 'in', 'is', 'japan', 'japanese', 'koria', 'kouchi', 'miyako', 'mori', 'most', 'nearly', 'no', 'osaka', 'region', 'russia', 'sapporo', 'secondary', 'sendai', 'shikoku', 'the', 'tokyo', 'we']
CountVectorizer: analyzer is char
results array:
[[6 0 1 2 1 1 0 2 0 2 1 4 1 1 0 0 2 2 1 0 2 4 0 0 2]
 [7 0 1 5 1 2 1 3 0 2 1 4 1 1 0 0 3 2 1 1 4 3 0 0 2]
 [9 1 1 6 0 2 0 1 2 0 1 5 1 3 1 2 4 4 1 3 2 2 2 0 2]
 [5 0 1 3 0 1 2 3 0 0 0 3 0 1 2 2 2 3 0 1 1 0 0 1 1]
 [8 0 1 5 0 1 0 1 1 0 0 4 1 0 1 2 3 4 3 4 5 2 1 0 2]
 [9 1 1 5 0 2 0 5 0 1 2 7 2 3 0 0 5 3 2 1 4 1 2 0 1]]
feature_names:
[' ', ',', '.', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'r', 's', 't', 'u', 'w', 'y']
CountVectorizer: analyzer is char_wb
results array:
[[14  0  1  2  1  1  0  2  0  2  1  4  1  1  0  0  2  2  1  0  2  4  0  0
   2]
 [16  0  1  5  1  2  1  3  0  2  1  4  1  1  0  0  3  2  1  1  4  3  0  0
   2]
 [20  1  1  6  0  2  0  1  2  0  1  5  1  3  1  2  4  4  1  3  2  2  2  0
   2]
 [12  0  1  3  0  1  2  3  0  0  0  3  0  1  2  2  2  3  0  1  1  0  0  1
   1]
 [18  0  1  5  0  1  0  1  1  0  0  4  1  0  1  2  3  4  3  4  5  2  1  0
   2]
 [20  1  1  5  0  2  0  5  0  1  2  7  2  3  0  0  5  3  2  1  4  1  2  0
   1]]
feature_names:
[' ', ',', '.', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'r', 's', 't', 'u', 'w', 'y']
CountVectorizer: analyzer is callable object
results array:
[[0 0 0 0 0 1]
 [0 0 1 0 0 0]
 [1 0 0 0 0 0]
 [0 0 0 0 1 0]
 [0 0 0 1 0 0]
 [0 1 0 0 0 0]]
feature_names:
['Fukuoka is most nearly city from Koria, China in japan.', 'Kouchi is a japanese city, in japanese region in Shikoku.', 'Osaka is the secondary biggest city in japan.', 'Sapporo is most nearly city from Russia in japan.', 'Sendai called we Mori no Miyako.', 'Tokyo is the biggest city in japan.']

結果から見えるもの

• analyzerにcharとchar_wbを指定した場合の結果を比較すると、analyzer = 'char'の場合の' 'の抽出回数をxとしたとき、analyzer = 'char_wb'の場合の' 'の抽出回数は2x+2回になる
  • analyzer = 'char'の場合は単純に文章中に含まれる' 'の数、analyzer = 'char_wb'の場合は加えて単語の前後が' 'ではない場合も境界として' 'に計上している
• analyzer = lambda x: [x]としたインスタンスで実行した場合、入力した文章が特徴（feature_name）としてそのまま戻ってくる。

予測された延滞率が実際どれだけ当たっているか。予測と結果の合致性の評価

はじめに

• 今回は、説明可能性についての3回目、モデルを評価するための指標の1つであるブライアスコアについてまとめています。
• MLで作成されたモデルをどう評価するか、どう説明するかについてさまざまな指標・方法を用いて行うことを目的としており、主にPoCなどで利用していくことを想定しています。
  
• 説明可能性についてのまとめはこちら。
  
  • POC作成のために、機械学習したモデルをどう評価し説明するかのまとめ。

Brier Score（ブライアスコア）

-　Brier Scoreとは

• 日本語ではWIKIもないブライアスコアですが、天気予報のWIKIの中にこんな記述があります。

  予報確率Fi、実際値Ai（発生した場合1、発生しなかった場合0などとする）としたときの確率誤差の総和を予報回数Nで割ったもの。値が0に近いほど精度が高い。

  $$\frac{1}{N} \sum_{i=1}^N(F_i - A_i)^2$$

• 英語版のWikiにある例を引用すると

  • 延滞予測が100%（$F_i$=1）で結果も延滞（$A_i$=1）の場合は、Brier scoreは0となり、最高値となります。
    
  • 延滞予測が100%（$F_i$=1）で結果は延滞しない（$A_i$=0）の場合は、Brier scoreは1となり、最低値となります。
    
  • 延滞予測が70%（$F_i$=0.7）で結果も延滞（$A_i$=1）の場合は、Brier scoreは$(0.70-1)^2$=0.09となります。
    
  • 延滞予測が30%（$F_i$=0.3）で結果も延滞（$A_i$=1）の場合は、Brier scoreは$(0.30-1)^2$=0.49となります。

• 予測値と結果に乖離があるほどスコアが下がる（Bier scoreの値が高くなる）ことになります。　

• 上記のように、確率は0から1の値となり、結果も0と1の2項分類となるときには機能しますが、それ以外の場合には別な方法を使用することとなります。別な方法についてはWikiのOriginal definition by Brierを参照してください。

• 検索して出てくる文献のほとんどは天気予報関連でした。たぶんそちらの業界ではかなり一般的な指標と思われますが、あまり金融系では使われていないようです。

• こちらでも書いた通りに、Forecastに拘りたい私としては、予測した確率の精度は重要な指標の1つとなります。

• クレジットスコアで置き換えて考えると、機械学習で得られた延滞確率が実際にはどれだけの確率で発生しているのか、その確率の誤差を評価できます。

- 利用方法

• 上記のように計算も割と単純ですが、Sklernにも準備されているので容易に使うことができます。
  

• サンプルは下記のようになります。

  import numpy as np
  from sklearn.metrics import brier_score_loss
  
  # Probabilityの値
  y_prob = np.array([0.1, 0.9, 0.8, 0.3])
  
  # 2項分類の場合  
  # 結果の値。
  y_true = np.array([0, 1, 1, 0])
  brier_score_loss(y_true, y_prob)
  0.037...
  

• ちなみに、多項分類なども準備されています。

  # 多項分類の場合  
  # 結果の値。
  y_true_categorical = np.array(["spam", "ham", "ham", "spam"])
  brier_score_loss(y_true_categorical, y_prob, pos_label="ham")
  0.037...
  

まとめ

• 0と1の予測結果の精度を評価するのがAUCやジニ係数などになります。今回のブライアスコアは何%の確率と予測して、その結果がどうだったのかを評価します。
• クレジットスコアでは、延滞する・しないの予測のほかに、何%で延滞するのかを適切に予測する必要があります。その確率が正しければ、延滞確率に見合う融資金額や利率、全体のポートフォリオでの許容アセットなどを適切にすることで、しっかりとリスクをコントロールしていくことができます。
• 学習する分類器によっては、予測値の値をそのまま使用すると正確ではない数値を利用することになります。それを防ぐため適切な予測値に変換するprobability calibrationを行ってから利用しましょう。たとえばランダムフォレストの場合は予測値が >0.1 <0.7の値になっていますので、他の分類器と比較する際は注意が必要です。
  

関連・参考記事

• 金融リスクと統計的学習
  
• sklearn.metrics.brier_score_loss

MFCCとは

MFCCは聴覚フィルタに基づく音響分析手法で、主に音声認識の分野で使われることが多いです。

最近だとニューラルネットワークに学習させる音声特徴量としてよく使われていますね。
2019.5.29訂正
Deep Learning for Audio Signal ProcessingによるとDeep Learningにおいては必要な情報が失われるためMFCCは使わずに、最後の計算ステップである離散コサイン変換を省いたメルスペクトラム(log-mel spectrum)が使われるそうです。MFCCは従来手法である隠れマルコフモデル、混合ガウスモデル、サポートベクターマシンで使われることが多いです。

今回はMFCC「メル周波数」や「ケプストラム」についても説明し、具体的なMFCCの実装方法も見ていきたいと思います。

メル尺度

心理学者のStanley Smith Stevensらによって提案された、人間の音高知覚が考慮された尺度です。

1000Hzの純音の高さの感覚を1000メルと決めた上で、1000メルの半分の高さに感じた音を500メル、1000メルの2倍の高さに感じた音を2000メルという容量で定めたものです。

実験により、可聴域の下限に近い音は高め、上限に近い音は低めに聞こえることがわかっています。

メル尺度は基底膜上の座標とほとんど一致します。1000メルは基底膜の1cmに相当し、臨海帯域幅は約137メルに相当します。

周波数$f$をメル尺度$m$に変換する式は

m = m_0\log(\frac{f}{f_0} + 1)

逆変換は

f = f_0(\exp\frac{m}{m_0} - 1)

で与えられます。

$f_0$は自由パラメータの周波数パラメータであり、$m_0$は「1000Hzは1000メル」という成約から導かれる

m_0 = \frac{1000}{\log(\frac{1000\mathrm{Hz}}{f_0} + 1)}

であらわされる縦続パラメータです。

よく用いられるパラメータとして

f_0 = 700, m_0 = 2595

があります。

周波数[Hz]とメル周波数[mel]の関係をグラフに表すとこのようになります。

ケプストラム

音声信号のフーリエ変換の絶対値に対数をかけ、さらに逆フーリエ変換（フーリエ変換とする場合もある）したもの。

信号$x(n)$のフーリエ変換を$X(e^{j\omega}$とすると、ケプストラム$c(m)$は

c(m) = \frac{1}{2\pi}\int_{\pi}^{\pi}\log|X(e^{j\omega})|e^{j\omega n}d\omega

高ケフレンシ領域は微細振動（音声でいうと声帯の振動）、低ケフレンシ領域はスペクトル包絡（音声でいうと声道フィルタの特性）を表します。

リフタリングという処理を施すことによってこれらを分離することができます。

変換手順

この記事を参考にしながら実際に変換してみます。

コードはほぼそのままですが、Python3対応させたり、変数名やパラメータをちょこちょこ変えたりしてます。

今回は「あ」の音声波形の定常部分のMFCCを求めてみましょう。（音声ファイルは各自用意してください。）

まずは必要なモジュールのインポートから

import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf

1. 波形を適当な長さに分割し、窓関数をかけてFFTを行う

# 音声の読み込み
master, fs = sf.read('wavfiles/a_1.wav')

t = np.arange(0, len(master) / fs, 1/fs)

# 音声波形の中心部分（定常部）を切り出す
center = len(master)//2 #　中心のサンプル番号
cuttime = 0.04 # 秒
x = master[int(center - cuttime / 2 * fs):int(center + cuttime / 2 * fs)]
time = t[int(center - cuttime/2*fs): int(center + cuttime/2*fs)]

plt.plot(time * 1000, x)
plt.xlabel("time [ms]")
plt.ylabel("amplitude")
plt.show()

上のコードを実行するとこのようなグラフが描けます。

次に窓処理を行って振幅スペクトルを求めます。

# ハミング窓をかける
hamming = np.hamming(len(x))
x = x * hamming

# 振幅スペクトルを求める
N = 2048 # FFTのサンプル数
spec = np.abs(np.fft.fft(x, N))[:N//2]
fscale = np.fft.fftfreq(N, d = 1.0 / fs)[:N//2]

plt.plot(fscale, spec)
plt.xlabel("frequency [Hz]")
plt.ylabel("amplitude spectrum")
plt.show()

2.メルフィルタバンクをかける

フィルタバンクは複数のバンドパスフィルタを並べたものです。

MFCCでは三角形のバンドパスフィルタをオーバラップさせながら並べ、フィルタバンクとしたメルフィルタバンクを用います。

メルフィルタバンクは周波数領域での三角窓がメル尺度上で等間隔になるように並べたものであり、低周波ほど幅が狭く、高周波ほど幅が広くなります。

並べられたバンドパスフィルタの数をチャンネル数といいます。

まずは周波数とメル尺度を変換する関数を実装します。

def hz2mel(f):
    """Hzをmelに変換"""
    return 2595 * np.log(f / 700.0 + 1.0)

def mel2hz(m):
    """melをhzに変換"""
    return 700 * (np.exp(m / 2595) - 1.0)

次にメルフィルタバンクを作る関数を実装します。

def melFilterBank(fs, N, numChannels):
    """メルフィルタバンクを作成"""
    # ナイキスト周波数（Hz）
    fmax = fs / 2
    # ナイキスト周波数（mel）
    melmax = hz2mel(fmax)
    # 周波数インデックスの最大数
    nmax = N // 2
    # 周波数解像度（周波数インデックス1あたりのHz幅）
    df = fs / N
    # メル尺度における各フィルタの中心周波数を求める
    dmel = melmax / (numChannels + 1)
    melcenters = np.arange(1, numChannels + 1) * dmel
    # 各フィルタの中心周波数をHzに変換
    fcenters = mel2hz(melcenters)
    # 各フィルタの中心周波数を周波数インデックスに変換
    indexcenter = np.round(fcenters / df)
    # 各フィルタの開始位置のインデックス
    indexstart = np.hstack(([0], indexcenter[0:numChannels - 1]))
    # 各フィルタの終了位置のインデックス
    indexstop = np.hstack((indexcenter[1:numChannels], [nmax]))
    filterbank = np.zeros((numChannels, nmax))
    print(indexstop)
    for c in range(0, numChannels):
        # 三角フィルタの左の直線の傾きから点を求める
        increment= 1.0 / (indexcenter[c] - indexstart[c])
        for i in range(int(indexstart[c]), int(indexcenter[c])):
            filterbank[c, i] = (i - indexstart[c]) * increment
        # 三角フィルタの右の直線の傾きから点を求める
        decrement = 1.0 / (indexstop[c] - indexcenter[c])
        for i in range(int(indexcenter[c]), int(indexstop[c])):
            filterbank[c, i] = 1.0 - ((i - indexcenter[c]) * decrement)

    return filterbank, fcenters

# メルフィルタバンクを作成
numChannels = 20  # メルフィルタバンクのチャネル数
df = fs / N   # 周波数解像度（周波数インデックス1あたりのHz幅）
filterbank, fcenters = melFilterBank(fs, N, numChannels)

# メルフィルタバンクのプロット
for c in np.arange(0, numChannels):
    plt.plot(np.arange(0, N / 2) * df, filterbank[c])

plt.title('Mel filter bank')
plt.xlabel('Frequency[Hz]')
plt.show()

次に振幅スペクトルにメルフィルタバンクをかけます。

振幅スペクトルに対して作成した各フィルタバンクをかけ、フィルタリングされた振幅を足し合わせて対数を取ります。

forループを使って書いてもいいですが、内積の形でも書け、こっちの方がエレガントです。

# 振幅スペクトルにメルフィルタバンクを適用
mspec = np.dot(spec, filterbank.T)

元の振幅スペクトルとフィルタバンクをかけた後のスペクトルを比べてみましょう。

# 元の振幅スペクトルとフィルタバンクをかけて圧縮したスペクトルを表示
plt.figure(figsize=(13, 5))

plt.plot(fscale, 10* np.log10(spec), label='Original Spectrum')
plt.plot(fcenters, 10 * np.log10(mspec), "o-", label='Mel Spectrum')
plt.xlabel("frequency[Hz]")
plt.ylabel('Amplitude[dB]')
plt.legend()
plt.show()

3.離散コサイン変換を行う

メルフィルタバンクをかけた後のスペクトルをケフレンシ空間に戻すために離散コサイン変換を行います。

離散コサイン変換を行う理由として、係数間の相関が減り、特徴量としての性能が向上するなどがあります。

離散コサイン変換の結果から低次の係数を取り出すと、MFCCが得られます。大体12次まで取ることが多いです。

from scipy.fftpack.realtransforms import dct

ceps = dct(mspec, type=2, norm="ortho", axis=-1)
nceps = 12
mfcc = ceps[:nceps]
print(mfcc)

[185.60782114 200.56668729  92.15345128   1.16538975 -47.4746488
 -60.60365351 -43.46459693   3.15481979  50.37136343  64.6038092
  48.42430931  27.41620554]

LibrosaでMFCCを求める

上では音声データ全体の中の１フレームのみを用いてMFCCを求めましたが、Librosaを使うと簡単に各フレームごとのMFCCを求めることができます。

import librosa

x, fs = sf.read('wavfiles/a_1.wav')
mfccs = librosa.feature.mfcc(x, sr=fs)
print(mfccs.shape)
print(mfccs[0])

(20, 26)
[-299.48698558 -306.19008252 -281.77019349 -205.61561205 -165.07338147
 -182.40750356 -172.10454081 -175.49027956 -191.51493931 -210.77326972
 -217.29547751 -223.91082068 -232.23909119 -247.65589292 -256.91669445
 -263.66001621 -275.46132992 -299.72892089 -288.90436044 -262.39782466
 -252.39945491 -256.61596037 -260.47453287 -274.17839766 -299.54783037
 -342.87809203]

カラーマップで可視化することも出来ます。

import librosa.display

librosa.display.specshow(mfccs, sr=fs, x_axis='time')
plt.colorbar()

おわりに

今回はMFCCに入門してみました。
GitHubにJupyter Notebookを置いていますので、よかったらこちらもどうぞ。

https://github.com/tmtakashi/signal_processing/blob/master/mfcc.ipynb

参考

リアルタイム音声認識（安藤彰夫　著）

音楽と機械学習 前処理編 MFCC ~ メル周波数ケプストラム係数

メル周波数ケプストラム係数（MFCC） - 人工知能に関する断創録

メル尺度 (Wikipedia)

http://abcpedia.acoustics.jp/acoustic_feature_2.pdf

ATTENTION

• チャンネルIDの取得方法は適当にググる (チャンネル名ではない)
• Slack App を登録して Tokenを得て、適切な Permission を設定しておく
• メンバー数が多い場合には対応していない、ページングに対応していない
• 重複招待でエラーが返ってくる問題にも対応していない
• API制限に引っかかったら適当に sleep 入れるなりして対処してください
• 例では無意味に jq コマンドかましてます

all-member-list.py

#!/usr/bin/env python3

# https://api.slack.com/methods/users.list

import requests, os, json

url = "https://slack.com/api/users.list"

params = {
  "token" : os.environ.get('SLACK_TOKEN'),
}

headers = {
  'Content-type': 'application/json'
}

res = requests.get(url, headers=headers, params=params)

print(json.dumps(res.json()))

invite-to-channel.py

#!/usr/bin/env python3

# https://api.slack.com/methods/channels.invite
# https://slack.com/api/channels.invite

import requests, os, json, sys

url = "https://slack.com/api/channels.invite"

members = json.loads(sys.stdin.read())

for member in members:
  params = {
    "token" : os.environ.get('SLACK_TOKEN'),
    "channel" : os.environ.get('CHANNEL'),
    "user" : member.get('id'),
  }

  headers = {
    'Content-type': 'application/json'
  }

  res = requests.post(url, headers=headers, params=params)

  print(json.dumps(res.json()))

exe

export SLACK_TOKEN=xxxxxxxxxxxxxxxx && ./all-member-list.py | jq -r '.members' | CHANNEL=xxxxxxx ./invite-to-channel.py

ワークスペースのデフォルトチャンネルを設定

一度全員招待しても、WorkSpaceのメンバーが増えたらやり直しなので、admin設定でデフォルトチャンネルも設定しておくと良いかも

https://example.slack.com/admin/settings

有料プランなら

• 全メンバーをSlackグループに入れておく
• グループに対して invite する

っていう手順で、API叩かなくてもできるかもしれません。

Original by Github issue

https://github.com/YumaInaura/YumaInaura/issues/2067

初カキコ…ども…

表題の通りですが、簡単に使い方を説明していこうと思います

simplekmlとは

簡単に言えばKMLをシンプルに作成できるpythonのパッケージですね。
あんまりsimplekmlに関する記事がないようでしたので、自分のコードとかを備忘録的に記事にしておこうと思いました。

記事作成時点でのバージョンは1.3.1になります。
https://simplekml.readthedocs.io/en/latest/index.html

めちゃ簡単なsample

simplekml.py

kml = simplekml.Kml()
kml.newpoint(name="Kirstenbosch", coords=[(18.432314,-33.988862)])  # lon, lat, optional height
kml.save("botanicalgarden.kml")

ざっくりこれだけでkmlが作れます。
sampleでは点だけを作ってますが、これ以外にもpolygonやlineが描画できます。
またフォルダを作成して、その中にpointやlineを作ることももちろん可能です。

ポリゴンのスタイル

polygonやline,pointにはそれぞれスタイルが設定できます。
設定できるものはオブジェクトによって違います。詳細は公式Documentを参照してください

polygonStyle.py

pol = kml.newpolygon(name = testName)
pol.outerboundaryis = [a, b, c, d, a]
pol.style.linestyle.color = "bf00ffff"
pol.style.polystyle.color = "3f00ff00"
pol.style.linestyle.width = 3

上のサンプルではまずポリゴンの外形境界線(outerBoundary)を設定し
境界線の色、ポリゴンの色、境界線の太さを設定しています。
また、このほかに内側の境界線(innerBoundary)を設定することもできます。

ポイントのスタイル

pointStyle.py

point = kml.newpoint(name = testName)
point.coords = [(a,b)]
point.style.iconstyle.icon.href = None
point.style.labelstyle.scale = 0.8

pointのスタイルになります。
割と使いそうなのはiconstyleだと思います。
こう記述するとiconなしで表示することができます。

注意点

最後に注意したい点なんですが、kmlの座標の表記は度分秒ではなく十進数表記になります。
若干めんどうなんですけど、そこだけ注意してください。

おわり

いかがでしたか！？
初めての記事なのでちょっと締めの言葉がわかりません。
転職やめ太郎師父の記事がめちゃおもしろいので皆さんも読んでください
https://qiita.com/Yametaro

以上です。

極めてニッチで一体誰の参考になるのかという小ネタだけど、自分用のメモとして。

前提

np で NumPy を使えるものとする。

import numpy as np

やりたいこと

以下のような要素のリスト (例: 1から6までの整数) があったとする。

arr = np.arange(1, 7)  # array([1, 2, 3, 4, 5, 6])

このリストから2つの要素を選ぶ全ての組み合わせについて、何らかの計算 (例: 2要素の積) をした行列を求めたとする。

X, Y = np.meshgrid(arr, arr)
matrix = X * Y
# array([[ 1,  2,  3,  4,  5,  6],
#        [ 2,  4,  6,  8, 10, 12],
#        [ 3,  6,  9, 12, 15, 18],
#        [ 4,  8, 12, 16, 20, 24],
#        [ 5, 10, 15, 20, 25, 30],
#        [ 6, 12, 18, 24, 30, 36]])

(補足) 上記のコードは、下記のコードでも同じ結果が得られる。

prod = lambda X, Y: X * Y
matrix = prod(*np.meshgrid(arr, arr))
# 上記と同じ行列が得られる

ここで、リストに対して要素が追加されることになった、とする。

arr_new = np.arange(7, 10)  # array([7, 8, 9])

もともとの要素 arr と追加要素 arr_new を結合した新しいリストについて、先ほど計算したような行列を再び求めたい。

愚直に書くとこうなる。

arr_all = np.concatenate([arr, arr_new])
matrix_new = prod(*np.meshgrid(arr_all, arr_all))
# array([[ 1,  2,  3,  4,  5,  6,  7,  8,  9],
#        [ 2,  4,  6,  8, 10, 12, 14, 16, 18],
#        [ 3,  6,  9, 12, 15, 18, 21, 24, 27],
#        [ 4,  8, 12, 16, 20, 24, 28, 32, 36],
#        [ 5, 10, 15, 20, 25, 30, 35, 40, 45],
#        [ 6, 12, 18, 24, 30, 36, 42, 48, 54],
#        [ 7, 14, 21, 28, 35, 42, 49, 56, 63],
#        [ 8, 16, 24, 32, 40, 48, 56, 64, 72],
#        [ 9, 18, 27, 36, 45, 54, 63, 72, 81]])

しかし、もともとの要素 arr に関する計算済みの行列 matrix があるので、この行列を利用して新しい要素に関する拡張部分だけを計算することにして、計算量を抑えたい。

この差分計算を NumPy で書きたい。

解法

このような関数を作ると実現できる。

def extend_matrix(matrix_base, arr_base, arr_new, func):
    part0 = func(*np.meshgrid(arr, arr_new))
    part1 = func(*np.meshgrid(arr_new, arr_new))
    return np.concatenate([
        np.concatenate([matrix_base, part0], axis=0),
        np.concatenate([part0.T, part1], axis=0),
    ], axis=1)

行列の新しく増える部分だけを計算して、既存の行列と結合している。

使うときは、こう。

matrix_new = extend_matrix(matrix, arr, arr_new, prod)
# array([[ 1,  2,  3,  4,  5,  6,  7,  8,  9],
#        [ 2,  4,  6,  8, 10, 12, 14, 16, 18],
#        [ 3,  6,  9, 12, 15, 18, 21, 24, 27],
#        [ 4,  8, 12, 16, 20, 24, 28, 32, 36],
#        [ 5, 10, 15, 20, 25, 30, 35, 40, 45],
#        [ 6, 12, 18, 24, 30, 36, 42, 48, 54],
#        [ 7, 14, 21, 28, 35, 42, 49, 56, 63],
#        [ 8, 16, 24, 32, 40, 48, 56, 64, 72],
#        [ 9, 18, 27, 36, 45, 54, 63, 72, 81]])

追加要素についての差分計算だけで全体の行列を求めることができた。

概要

Pythonでmultiprocessingを使って並列処理を実行するときに、それぞれのプロセス毎に進捗を表示したい時があるかと思います。
そこで、tqdm を用いてプロセス毎にバーを表示してみました。

参考 : tqdm

目標

プロセス毎に独立してプログレスバーを表示！

Code

from time import sleep
import random
from tqdm import tqdm
from multiprocessing import Pool, freeze_support, RLock

L = 8


def long_time_process(p):
    info = f'#{p:>2} '  # 進捗バーの左側に表示される文字列
    for _ in tqdm(range(20), desc=info, position=p+1):
        sleep(random.random())
    return p * 2


if __name__ == '__main__':
    freeze_support()  # Windows のみ必要
    with Pool(L,
             # Windows のみ必要
             initializer=tqdm.set_lock, initargs=(RLock(),)) as p:
             result = p.map(long_time_process, range(L))
    print("\n" * L)  # tqdm終了後のカーソル位置を最下部に持ってくる
    print(result)

解説

long_time_process: 並列実行したい関数

freeze_support, initializer, initargs: Windowsでのみ必要になるそうです。Macではこれを省いても特に問題ありませんでした。

position=p+1: 本来はpositionは0から始まるため、画像を見るとわかるように、プログラムを実行したときに上に1行分のスペースができています。position=0から表示させると、一番上のプログレスバーが終了したタイミングで全体の表示がずれてしまうため、ここでは1から始めています。

結果

レイアウトが崩れることなく、プログレスバーの表示、及びその後のprint出力もできています。

TurnoverをQuantX上で算出

• positive returns：収益が正である
• Low exposure to sector risk and style risk：Quantopian独自の Risk model の条件を満たしている
• Beta-to-spy：
• Leverage：レバレッジ0.8-1.1倍（QuantXには関係ない？）想定資本に対して８０％〜１１０％の範囲で投資
• turnover : 中程度のturnoverを維持。
• アルゴリズムに用いる株はある程度の流動性を持っていないといけない
• position concentration ：アルゴリズムは１つの銘柄に対して全資産の５％以上投資してはいけない
• Long and short positions ：買いポジションと売りポジションの差額がトータル資金に対して１０％以上になってはいけない。

以下補足

• Risk model は大きく分けて２つ。セクター要因とスタイル要因
• セクター要因とは株式会社を１１のセクターに業界で分別して（ex,技術、金融、不動産...）、各セクターに２０％以上の資金を投資していはいけない
• スタイル要因とは？？？momentum size volatility とかで計算していく
• turnover ＝（１日の投資資金の絶対値）/（全ポートフォリオの価値）
• コンテストに使用できる株が流動性のある株に限定するために、Quantopianは独自で流動性があるとみなした株式銘柄群（株の集合）を作成し、コンテストに使用する株はその銘柄群の中から選ぶ必要がある。

この中で投資比率制約、turnover制約、を実際に実装してみてこれらを満たしているかを確認する。

turnover制約の実装

turnoverの具体的な計算方法としては（１日の投資資金の絶対値）/（全ポートフォリオの価値）を毎日計算して、それの６３日営業日移動平均をとる。この値が５〜６５％の間に入っていることが条件。

X日目のturnover = \sum_{銘柄} \frac{X日目の、ある銘柄の投資量}{全資産}
= \sum_{銘柄} {|(X日目の、ある銘柄の全資産に対する保有割合)-(X-1日目の、ある銘柄の全資産に対する保有割合)|}

で算出できる。

以下でQuantX上でturnoverを表示する

X日目の、ある銘柄の全資産に対する保有割合 current_portfolio_ratio は、symで銘柄を指定して

current_portfolio_ratio = ctx.portfolio.positions[sym]["portfolio_ratio"]

で取得できる。
さらにあらかじめctx.localStorageを利用して保存しておいた、X-1日目のある銘柄の全資産に対する保有割合 previous_portfolio_ratio も用いて

X日目のturnoverは

X日目のturnover = \sum_{銘柄} |current_portfolio_ratio - previous_portfolio_ratio|

で求まる。ここでQuantopian上での制約はturnoverの63日移動平均を計算してそれが5~63%に入ることが条件である。
ゆえに、上で算出したX日目のturnoverを毎日localStorageに保管して、実際にはlocalStorageに保管されている過去63日分のturnoverの平均を数値評価するのである。

以上を実際にQuantX上で実装してみる。
とりあえず、取り扱う通過などを設定する

import pandas as pd
import talib as ta
import numpy as np

def initialize(ctx):
    ctx.flag_profit = False #利益確定売りを用いるかTrueなら用いるFalseなら用いない
    ctx.flag_loss = False #損切りを用いるかTrueなら用いるFalseなら用いない
    ctx.loss_cut = -0.03 #損切りのボーダーマイナス％
    ctx.profit_taking = 0.05 #利益確定売りのボーダープラス％
    ctx.codes = [
      4502,
      4506,
      ]
    ctx.symbol_list  = ["jp.stock.{}".format(code) for code in ctx.codes]
    ctx.configure(
      channels={          # 利用チャンネル
        "jp.stock": {
          "symbols":ctx.symbol_list,
          "columns": [
            "close_price",        # 終値
            "close_price_adj",    # 終値(株式分割調整後) 
                     ]}})

次にmy_signal内でbuy_sigとsell_sigを設定する
しかし今回はわかりやすいようにするためにbuysigの値を全てTrueに、sell_sigの値を全てFalseに設定する。

     def _my_signal(data):
      # この部分に作成するアルゴの指標を書き込んで下さい。
      #シグナル定義
      #各銘柄の終値(株式分割調整後)を取得、欠損データの補完 
      cp = data["close_price_adj"].fillna(method="ffill")

      #単純移動平均線(SMA)の設定
      #データの入れ物を用意
      s_sma = pd.DataFrame(data=0,columns=[], index=cp.index)
      m_sma = pd.DataFrame(data=0,columns=[], index=cp.index)
      l_sma = pd.DataFrame(data=0,columns=[], index=cp.index)
      #TA-Libによる計算
      for (sym,val) in cp.items():
        s_sma[sym] = ta.SMA(cp[sym].values.astype(np.double), timeperiod=5)
        m_sma[sym] = ta.SMA(cp[sym].values.astype(np.double), timeperiod=25)
        l_sma[sym] = ta.SMA(cp[sym].values.astype(np.double), timeperiod=75)
      #SMAの売買シグナルの定義
      s_sma_ratio = s_sma/m_sma
      l_sma_ratio = m_sma/l_sma

      buy_sig = (s_sma_ratio > 1.00) & (s_sma_ratio < 1.10) & (l_sma_ratio > 1.00) & (l_sma_ratio < 1.10)
      sell_sig = (s_sma_ratio < 0.99) & (l_sma_ratio < 0.99) 

       #実験的にbuy_sigの値を全てtrueに、sell_sigの値を全てFalseに設定。
      buy_sig.iloc[:,:] = True
      sell_sig.iloc[:,:] = False


      return {
        "sma5":s_sma,
        "sma25":m_sma,
        "sma75":l_sma,
        "buy:sig": buy_sig,
        "sell:sig": sell_sig,
      }

    # シグナル登録
    ctx.regist_signal("my_signal", _my_signal)

さらにhandle_signal内で売買を行い、その売買結果を利用してturnoverの値を算出する。
なお設定としてbuy_sigにヒットすると、一回の取引で全資産の5%分の株を買うように設定している。
実際には全資産の５%ちょうど売買することは不可能で、5%未満でキリのいいunit数売買することになる。
また、quantopianではturnoverの63日の移動平均をとって数値評価しているが、今回の試験的なQuantX上でのturnoverの実装なので、簡単にturnoverの5日移動平均をとって数値評価する。

def handle_signals(ctx, date, current):
    '''
    current: pd.DataFrame
    '''
    df = current.copy()
    done_syms = set([])
#X日目の各銘柄の全資産に対する保有割合の変化を足し合わせるためにsum_rate_changeを作成
    sum_rate_change = 0

    for (sym,val) in ctx.portfolio.positions.items():
        current_portfolio_rate = ctx.portfolio.positions[sym]["portfolio_ratio"]

        '''
        localStorageにsymbolのkeyがあれば先に進み、なければsymbolをkey、
        valueは空の辞書とする、辞書を作る
        '''

        if sym in ctx.localStorage:
          pass
        else:
          ctx.localStorage[sym] = {}

        '''
        localStorage[sym]に"portfolio_ratio"のkeyがあれば、
        前日と当日のportfolio_ratioの差の絶対値を計算しsum_rate_changeにたす
        なければ、新しくlocalStorage[sym]にからの辞書を作る
        '''

        if "portfolio_ratio" in ctx.localStorage[sym]:

          previous_portfolio_rate = ctx.localStorage[sym]["portfolio_ratio"]
          abs_rate_change =  abs(current_portfolio_rate - previous_portfolio_rate)
          sum_rate_change += abs_rate_change
          ctx.localStorage[sym]["portfolio_ratio"] = current_portfolio_rate
        else:
          ctx.localStorage[sym]["portfolio_ratio"] = current_portfolio_rate

        returns = val["returns"]
        if (ctx.flag_loss) & (returns < ctx.loss_cut):
          sec = ctx.getSecurity(sym)
          sec.order(-val["amount"], comment="損切り(%f)" % returns)
          done_syms.add(sym)
        elif (ctx.flag_profit) & (returns > ctx.profit_taking):
          sec = ctx.getSecurity(sym)
          sec.order(-val["amount"], comment="利益確定売(%f)" % returns)
          done_syms.add(sym)

    if "turnover" in ctx.localStorage.keys():
      ctx.localStorage["turnover"] = np.append(ctx.localStorage["turnover"], float(sum_rate_change))
      sma_turnover = ta.SMA(ctx.localStorage["turnover"], timeperiod=5)
      ctx.logger.debug(sma_turnover[-1])
    elif "turnover" not in ctx.localStorage.keys():
      ctx.localStorage["turnover"] = np.zeros(1)
    ctx.logger.debug(ctx.localStorage['turnover'])



    # 買いシグナル
    df_long = df[df["buy:sig"]]
    # ctx.logger.debug(df_long)
    if not df_long.empty:
      for (sym, val) in df_long.iterrows(): 
        if sym in done_syms:
          continue

        sec = ctx.getSecurity(sym)
        sec.order_percent(0.05, comment="SIGNAL BUY")

    # 売りシグナル
    df_sell = df[df["sell:sig"]]
    if not df_sell.empty:
      for (sym, val) in df_sell.iterrows():
        if sym in done_syms:
          continue

        sec = ctx.getSecurity(sym)
        sec.order_percent(0, comment= "SIGNAL SELL")

これを実行するとdebug画面に
array型のturnoverと、そのarrayの５日移動平均をとった実際に数値評価するsma_turnoverが得られる。

ところで今回、二銘柄に対してsignalが発火するごとに5％づつ取引するように設定しているので、１日の取引量の全資産に対する割合はおおよそ８〜１０％であると予想される。

実際の結果は以下である。（以下はn日目のturnoverのarrayと）

# turnoverのarray
[0.         0.         0.         0.08652656 0.08774072 0.08964686
 0.08793647 0.08958893 0.08674776 0.08879918 0.09168621 0.08819629]

turnoverのarrayの最新５日分の平均をとった値。
0.0888118524395898

となり予想して結果と一致している。

以上よりturnoverのQuantX上での計算が完了

はじめに

Openpyxlを使って、エクセルを読み込み、マージされたセルをどのように扱うのかについてメモを残しておきます。

このようなデータを使用。

Excelデータを読み込む

import openpyxl

# エクセエルファイルを指定し、読み込みたいシート名を指定
excel_file = 'sample.xlsx'
wb = openpyxl.load_workbook(filename = excel_file)
sheet = wb['Sheet4']

セル範囲を指定してセルを取得する

# 先ほど読み込んだシートの中で読み込みたい範囲を指定する
for 行変数 in シート変数[開始セル:終了セル]:
 for セル変数 in 行変数:
  print(セル変数.value ) 

例えば、A1~D7に書かれたデータを読み込みたい場合は以下のように範囲を指定する

# １回目のfor文では、行単位で、2回目のfor文で行のなかにあるcell単位で出力することができる
for row in sheet['A1:D7']:
 for cell in row:
  print(cell, cell.value) 

出力結果

<Cell 'Sheet4'.A1> None
<Cell 'Sheet4'.B1> サンプル
<MergedCell 'Sheet4'.C1> None
<MergedCell 'Sheet4'.D1> None
<MergedCell 'Sheet4'.A2> None
<Cell 'Sheet4'.B2> レコードヘッダー1
<Cell 'Sheet4'.C2> レコードヘッダー2
<Cell 'Sheet4'.D2> レコードヘッダー3
<MergedCell 'Sheet4'.A3> None
<MergedCell 'Sheet4'.B3> None
<MergedCell 'Sheet4'.C3> None
<MergedCell 'Sheet4'.D3> None
<MergedCell 'Sheet4'.A4> None
<MergedCell 'Sheet4'.B4> None
<MergedCell 'Sheet4'.C4> None
<MergedCell 'Sheet4'.D4> None
<Cell 'Sheet4'.A5> カラムヘッダー1
<Cell 'Sheet4'.B5> サンプルセル1
<Cell 'Sheet4'.C5> サンプルセル2
<Cell 'Sheet4'.D5> サンプルセル3
<Cell 'Sheet4'.A6> カラムヘッダー2
<Cell 'Sheet4'.B6> サンプルセル4
<Cell 'Sheet4'.C6> サンプルセル5
<Cell 'Sheet4'.D6> サンプルセル6
<Cell 'Sheet4'.A7> カラムヘッダー3
<Cell 'Sheet4'.B7> サンプルセル7
<Cell 'Sheet4'.C7> サンプルセル8
<Cell 'Sheet4'.D7> サンプルセル9

cellオブジェクトから得られる情報は
マージされたセルか(Cell or MergedCell)
セルの属するシート名 (シート名)
セルの座標 (アルファベット+数字)

cell.valueから得られる情報は
セルの中に書かれたテキスト
(* mergeされたcellの場合はNoneになる)

mergeされたセルの値

セルが縦方向または横方向にマージされている場合、そのセルの一番左上が代表のセルとなり、
(縦方向でのマージでは一番上にあるセル、横方向のマージでは一番左にあるセル)
それ以外のセルはMergedCellとなる。

具体的には、
今回のデータでは、サンプルと書かれたセルが、B1~D1にかけてマージされている。

この時、先ほどの実行結果を見てみると、

<Cell 'Sheet4'.B1> サンプル
<MergedCell 'Sheet4'.C1> None
<MergedCell 'Sheet4'.D1> None
<MergedCell 'Sheet4'.A2> None
<Cell 'Sheet4'.B2> レコードヘッダー1
<Cell 'Sheet4'.C2> レコードヘッダー2
<Cell 'Sheet4'.D2> レコードヘッダー3

左上のB1のセルは'Cell'として、"サンプル"というテキストをもち、
それ以外のC1, D1は'MergedCell'となり、valueがNoneとなる。

同様に、レコードヘッダー1, レコードヘッダー2, レコードヘッダー3と書かれた
縦にマージされたセルについて、一番上のB2, C2, D2のみテキストが入り、
それ以外のセルは'MergedCell'とされる。

<Cell 'Sheet4'.B2> レコードヘッダー1
<Cell 'Sheet4'.C2> レコードヘッダー2
<Cell 'Sheet4'.D2> レコードヘッダー3
<MergedCell 'Sheet4'.A3> None
<MergedCell 'Sheet4'.B3> None
<MergedCell 'Sheet4'.C3> None
<MergedCell 'Sheet4'.D3> None
<MergedCell 'Sheet4'.A4> None
<MergedCell 'Sheet4'.B4> None
<MergedCell 'Sheet4'.C4> None
<MergedCell 'Sheet4'.D4> None

mergeされたセルの範囲を取得

sheetの中のマージされたセルの範囲を取得する関数が用意されていたのでメモ。

sheet.merged_cells.ranges

• merged_cell_rangesが廃止されて、バージョンいくつからmerged_cells.rangesを使用することになったそう。 https://openpyxl.readthedocs.io/en/stable/api/openpyxl.worksheet.worksheet.html#openpyxl.worksheet.worksheet.Worksheet.merged_cell_ranges

出力結果

[<CellRange B1:D1>,
<CellRange B2:B4>,
<CellRange C2:C4>,
<CellRange A1:A4>,
<CellRange D2:D4>]


これを使えば、どのセルがどこからどこまでマージされているのかリストで取得できるため、
先ほどの実行結果と組み合わせれば、色々使えそう！