DX担当者になったら考えたいこと

「DX(デジタル・トランスフォーメーション)」というキーワードを目にすることが大変多くなった。

この5年間の間のだけでもIndustry 4.0, Society 5.0, ディープラーニング, AIといった様々なキーワードが入れ替わり立ち替わり似たようなメッセージを発し続けている。 それらの描く美しい未来像に触れ技術の要素の習得に力を入れつつも、 思い通りにならない会社組織を恨めしく思ったり、社会の重鈍さに腹を立てたりしてきたエンジニアは少なくないと思う。 あるいは一歩進んで、そういった現状を言語化・図示した資料に出会い、社内でのICT介護のやりとりを思い浮かべながら、 「よくぞ言ってくれた、これを檄文にみんなを決起させよう」と気持ちを高ぶらせているところなのかもしれない。

私自身そんなエンジニアの一人なのかもしれない。

だが、これは「無敵の『DX』でなんとかしてくださいよォーー」の前兆ではないのか?

「『DX』は手段です」と、したり顔で公演する識者の講演に首がもげそうなくらい頷くモブへと一直線に舗装された道の入り口に立っているのではないか?

そんな未来に到達する前に少しだけ悪あがきをしてみたい。

DX担当者に任命されたら…

ある時、あなたは上司に呼び出されたとしよう。

なんと「社長から『DX』を推進するように指示があった。『DX』は何やらパソコンを活用する話で若い人の柔軟な考え方が大事だそうだから君に任せたい。」という。

なんとも危ない匂いのする話だ。どこかでダメと書いてあった始まりそのものではないか。

社長の頭の中

さて、よくわかりもしないDXをまるっと若手のあなたに押し付けた社長は全くのアホで要高度ICT介護の上司なのだろうか。

だが、特に一部上場企業などのサラリーマン社長であれば、実績を積み、類まれな何らかの才能を発揮してその立場についていると考えるのが妥当だろうし、その立場から様々な形で一流の指南役も社内外についていると見るのが良いと思う。

私自身もしがないサラリーマンで、社長の経験もなくなったこともない立場の人間の経験は想像でしか語れないが、幸い、"CEO Job Description"とググるとCEO(≒社長)の役割についていくつもの記事があり、 ドラッカーなど著名な経営研究家が色々な名言を残してくれている。

また、CEOの直面する「課題」などもアンケートが取られていたりする。

そういった情報をまとめると

  • (サラリーマン)社長の最も重要なミッションは「会社を永続させること」
であり、そのために会社に関するメッセージを顧客・社員・投資家・政府などあらゆる方面とやりとりし、 それらをまとめ上げながら実行に移していく(指示を出す)といった役割を担っているようだ。

会社が永続できなくなる落とし穴は山ほどある。

もっとも当たり前のものは既存事業のシェアを落とすなどしてコストが売上より大きくなり赤字転落後そのまま運転資金が枯渇すること。「永続」すれば良いと言ってゾンビ飛行ギリギリとなってしまえば、成長が期待できない会社として投資家に愛想をつかされ、これもまた資金調達が苦しくなる。

だが、現在においては、人種や性別などの多様性対応を誤り、顧客離れ、製品の生産やサービスの提供が止まってしまったり、温暖化などの環境対策を狙いとした政策への対応に失敗し厳しい罰則を受けたり、 品質検査の不正行為で賠償に発展したり、昨今コロナ禍で顕在化したように市場が蒸発したりしてしまうこともある。

全部優秀なスタッフに任せることができれば良いが、意外と社長によるきめ細かなマニュアル運転が必要になってしまっている企業は多いと思う。

そんな中、投資家、つきあいのある政治家(官庁/政治団体)、社外取締役から「御社のDXはどうなっていますか?」とくると、何かしらの回答準備が必要になってくる。

社長は思う。

「Industry 4.0, Society 5.0, AI, 自分がポジションに付く前に似たような話はたくさんあったじゃないか。なんで今また夏休みの宿題を最後の1日でこなさないとならないような状況に…。とはいえ、なんとか形にしていかないと。」

デジタルネイティブな会社であれば、社長になったあとでデジタル技術に対する意識が急激に高まるという間抜けな話にはならないだろう。だが、生産や営業といった役割分化した組織の中で上がってくると、いくらゼネラリスト養成学校のJTCといえども壁の向こうを本当の本気で見るのは社長になったときくらいだ。

DX担当者として…

過小評価かもしれないがこうやって社長の頭の中を仮定すると、DX担当者として最初に着手することは「会社の永続」と「DX」を結びつけることだ。

そして、その方向性としては大きく分けると

  1. 既存事業のExploit(売上増加[アップセル]/コスト削減)
  2. 新規事業のExplore
の2つになる。

幸い、この方向でまとめている記事を見つけたので紹介する。

この記事と重複する部分もでてくるが、次回は、ここからさらに具体的な打ち手の選択肢へと掘り下げたい。

メモ

vueのプロジェクト設定メモ

メモの要略

vue-cliで作ったSingle Page Webアプリのプロジェクト設定でつまづいたところのメモ

  1. 環境変数(.envファイル)をvueのプロジェクトで読み込んで使う
  2. Vueファイルに対してESlintとPrettierを設定し、vscodeで使う

前提環境

  • vue
    • vue/cli: 4.5.9
    • vue: 2.6.12
  • eslint: 7.16.0
  • prettier: 2.2.1

1. 環境変数の読み込みをできるようにする

公式ドキュメント(Modes and Environment Variables | Vue CLI)によると、vue-cliでは.envファイルで定義した環境変数をプロジェクトやpublic/index.html内の埋め込み変数として使えるとのこと。 しかし、何が悪いのかはっきりしなかったが、デフォルトでは.envファイルを読み込まなかったので注意として記しておく。

解決策

プロジェクトにdotenv-webpackを追加すること

yarn add -D dotenv-webpack

.envファイルの作り方での注意

公式ドキュメント通り

  • .envファイル内で定義する変数名はPREFIXとしてVUE_APP_を与え、VUE_APP_SOMETHINGという名前にすること
  • .envファイルはpackage.jsonと同じ、プロジェクトのルートに置くこと
  • .envファイルを更新したら開発サーバ(yarn serve)を再起動すること

環境変数の使い方

.envファイルでVUE_APP_SOMETHINGという変数を定義したとして

  • .js/.vueではprocess.env.VUE_APP_SOMETHINGで参照する
  • public/index.htmlでは<%= VUE_APP_SOMETHING %>で参照する

参考

vue.js - Vue-cli 3 Environment Variables all undefined - Stack Overflow
同じような症状にあった方がいる様子。Dylan Nguyenさんがdotenv-webpackのインストールが必要だった旨に言及。

2. Vueファイルに対してESlintとPrettierを設定

2020年11月頃にPrettierのLintとの併用における公式設定が変ったらしいとのことで、Qiitaなどに転がっていた2020年6月頃の情報が古くなっていたことから試行錯誤。

これを読んでいるあなたも一旦公式ドキュメントを確認したほうが良いかも。情報古くなるの早すぎ…。

必要なパッケージ

プロジェクトに追加するものとして

yarn add -D eslint prettier eslint-plugin-vue eslint-config-prettier eslint-plugin-prettier
  • eslint: 7.16.0
  • prettier: 2.2.1
  • eslint-plugin-vue: 7.3.0
  • eslint-config-prettier: 7.1.0
  • eslint-plugin-prettier: 3.3.0

VSCodeに追加するものとして

設定ファイル

prettierrcは作らずにeslintrcに含めてしまいました

vscodeのsetting.json

veturを入れている多いのでvalidationを無効にする設定とeslintの有効対象設定. 心を込めた手作業にするためautosaveなどはonにしてない。

{
  "vetur.validation.template": false,
  "eslint.workingDirectories": ["./"],
  "eslint.validate": ["vue", "javascript", "html"]
}

.eslintrc.js

module.exports = {
  root: true,
  env: {
    node: true,
    browser: true,
    es6: true,
  },
  plugins: ['prettier'],
  extends: [
    'plugin:vue/recommended',
    'eslint:recommended',
    'prettier',
    'prettier/vue',
  ],
  rules: {
    'no-console': process.env.NODE_ENV === 'production' ? 'error' : 'off',
    'no-debugger': process.env.NODE_ENV === 'production' ? 'error' : 'off',
    'prettier/prettier': [
      'error',
      {
        endOfLine: 'auto',
        singleQuote: true,
        semi: false,
      },
    ],
  },
  parserOptions: { parser: 'babel-eslint' },
}

参考

GMail APIとPandasでメールボックスの大掃除

Googleアカウントの保存ポリシー変更の件、みなさんいかがされていますか?

Google ストレージの仕組みに関する今後の変更 - Gmail ヘルプ

Google Oneも契約しているので容量には余裕があるのですが、GMailには二度どころか一度も読むことのないようなメールが山と溜まっています。年末にもなりましたし掃除でもするかとポチポチ消しはじめ、 楽○さんやY○hooさんが熱心にお勧めしてくださっている人生を豊かにするサービスや商品に関するメールを心苦しながら読まぬうちに削除し、 Promotionカテゴリ下に28,000件あった未読メールのうち、14,000件まで減りましたが、 はて、それでも14,000件ある…残りがどこから来たのか気になったので分析することにしました。

GMail削除UI改善のお知らせ

これまでGMailの検索結果の「全部選択」はリスト表示されているページのメールだけがまとめて選択できたのはご存知でしょうか?本日付で試してみたところ「全部選択」後に選択範囲を「検索結果全体に拡張」できるようになっていたので、まとめて削除が大変捗るようになっています。

Work Package

  • GMail APIでメールリストを取得
  • Pandasで前処理
  • グラフ化

環境

  • Python 3.7
  • Windows 10
  • Scoop
    • wget

GMail APIでメールリストを取得

準備

Python Quickstart  |  Gmail API  |  Google Developersに従って何も迷うことはありません。 APIを有効化し、Desktopアプリ用のcredential.jsonを取得します。これはPythonコードと同じディレクトリに置きましょう。

pip install --upgrade google-api-python-client google-auth-httplib2 google-auth-oauthlib
wget https://raw.githubusercontent.com/googleworkspace/python-samples/master/gmail/quickstart/quickstart.py
python quickstart.py

スクリプト実行の際、ブラウザが開いて認証をするとtoken.pickleが保存され、以降のアクセスではこれを使うことができるようになります。

メールリスト取得

コードの全体はこちら⇒Fetch GMail List

メールのリストを取得してmessage.jsonというファイルに保存します。

メールIDの取得

messages().listのAPIでは、GMail内部で割り振られたメールIDを取得できます。

広告に分類されたメールを引っこ抜きます。リファレンスではlistでは好きな件数だけ情報を引き抜けるような事になっているようですが、実態としては1回のリクエストで500件までのようです。

try:
    maxResults = 500
    request = service.users().messages().list(userId='me', labelIds=['CATEGORY_PROMOTIONS'], maxResults=maxResults)
    while request is not None:
        results = request.execute()
        consumed_token += 5
        messages = results.get('messages', [])
        db += messages
        request = service.users().messages().list_next(request, results)
finally:
    df = pd.DataFrame(db)
    df['from'] = pd.Series(dtype='string')
    df['subject'] = pd.Series(dtype='string')
    df['timestamp'] = pd.Series(dtype='string')

    df.to_json('messages.json', orient='records')

メールヘッダー情報の取得: 'From', 'Subject', 'Date'

APIを叩ける回数を調べると1秒間に50件のメール情報まで引き抜けるようでしたので、ちょっとSleepさせてます。(Usage limits  |  Gmail API  |  Google Developers)

try: 
    for i in tqdm.tqdm(df.index):
        if not pd.isnull(df['from'][i]):
            # Skip fetched data
            continue
        result = service.users().messages().get(userId='me', id=df['id'][i], format='metadata', 
                                                metadataHeaders=['From','Subject','Date']).execute()
        consumed_token += 5
        body = result.get('payload',[])

        if not body:
            print("Token Limit per minutes may be reached")
            break

        for el in body['headers']:
            if el['name'] == 'From':
                df.loc[i, 'from'] = el['value']
            elif el['name'] == 'Subject':
                df.loc[i,'subject'] = el['value']
            elif el['name'] == 'Date':
                df.loc[i, 'timestamp'] = pd.Timestamp(el['value'])
        time.sleep(1/40)
finally:
    df.to_json('messages.json', orient='records')
    print("Consumed Token {}".format(consumed_token))

Pandasで前処理

さて、ここから取得した情報をjupyter notebookで内容分析します。

import numpy as np
import seaborn as sns
import matplotlib.pylab as plt
from matplotlib.ticker import PercentFormatter
import pandas as pd

import re

df = pd.read_json('messages.json')

前処理

まず、From形式情報からemailアドレスを取り出します。

def domain(row):
    m = re.search('\<(.*)\>', row)
    if m is None:
        return row
    return m.group(1)
df['email'] = df['from'].apply(domain)
# HEAT MAPを作るときに使う変数
df['Year'] = df['timestamp'].apply(lambda x: x.strftime('%Y'))
df['Month'] = df['timestamp'].apply(lambda x: x.strftime('%m'))

メール送信者の集計

aggregations = {
  'email': ['count']
}
grouped = df.groupby(['email']).agg(aggregations)
grouped.columns = ["_".join(x) for x in grouped.columns.ravel()]
grouped = grouped.sort_values('email_count', ascending=False).reset_index()

グラフ化

メール送信者のトップ20とメールボックスに占める割合を算出します

grouped["cumpercentage"] = grouped["email_count"].cumsum()/grouped["email_count"].sum()*100

fig, ax = plt.subplots()
ax.bar(grouped.email[0:20], grouped["email_count"][0:20], color="C0")
ax2 = ax.twinx()
ax2.plot(grouped.email[0:20], grouped["cumpercentage"][0:20], color="C1", marker="D", ms=7)
ax2.yaxis.set_major_formatter(PercentFormatter())

ax.tick_params(axis="y", colors="C0")
ax2.tick_params(axis="y", colors="C1")
for tick in ax.get_xticklabels():
    tick.set_rotation(90)
plt.show()

ついでにどのくらい熱心にメールを下さっているのかチェック

p = pd.pivot_table(df[df['email']==grouped['email'][0]], values='id', columns=['Year'] , index=['Month'], aggfunc='count')
sns.heatmap(p, cmap=sns.light_palette("seagreen", as_cmap=True),  linewidths=.5)

4年以上メールくれてる

まとめ

  • トップ20のメール送信者を削除すれば50%くらいが削除できそう
  • 読んでない広告メールは購読解除しよう

遠隔なjupyter notebookに接続する話

ポートフォワーディングのメモ書き

localhost:8888localhost:8000192.168.1.100sshN f -L 8000:localhost:8888 user@192.168.1.100sshN f -R 8000:localhost:8888 user@192.168.1.3 192.168.1.3

ローカルの端末からリモートのjupyterを呼び寄せる場合

ssh –N –f -L 8000:localhost:8888 user@192.168.1.100

リモートの端末からローカルにjupyterを送りつける場合

ssh –N –f -R 8000:localhost:8888 user@192.168.1.2

Xをフォワードする場合

LinuxサーバにWindowsからPower ShellでOpenSSHをクライアントにしてアクセスする

$env:DISPLAY= 'localhost:0.0'
ssh -Y user@192.168.1.2

別途、VcXsrvを忘れないように。WSL上にXサーバをインストールしてGUIを実現する

ナイキストの安定判別法と小ゲイン定理

H∞制御の中心的役割を担う小ゲイン定理で求める制御安定(十分)条件はナイキストの安定判別法などから求まる必要十分な安定条件よりも保守的な条件として知られている。

一巡伝達関数
\[ P(s) = \frac{K_p}{\tau s+1}e^{-Ls} \]
で与えられる閉ループ系が安定であるための$K_p$の条件を求めよ

ナイキスト線図を用いて安定条件を求めます。PI制御チューニング : Ziegler Nicholsの限界感度法 - むだ時間ありの1次遅れ系への適用から、

\[ K_p \in \left[0, \sqrt{(\tau \omega_0)^2 + 1}\right) \qquad s.t. \arctan(\tau \omega_0) = \pi - L \omega_0 \]

小ゲイン定理では

\[ \Vert P(j\omega) \Vert_{\infty} < 1 \]
となる$ K_p $の範囲が求める範囲。線形時不変システムでは
\[ \Vert P(j\omega) \Vert_{\infty} = \sup_\omega \vert P(j\omega) \vert < 1 \]
です。これから
\[ \Vert P(j\omega) \Vert_{\infty} = \sup_\omega \left| \frac{K_p}{\sqrt{(\tau \omega)^2 + 1}} \right| = K_p \]
なので
\[ K_p \in [0, 1) \]
ここで、$\omega, \tau > 0$なので明らかに$[0,1)\subset{\left[0, \sqrt{(\tau \omega_0)^2 + 1}\right)}$であって、 保守的な条件になっていることがわかります。

ベンチャー投資のディールフロー

リモートワーク主体・給与体系が超透明で有名なGitlabですが、スタートアップ買収のプロセスまで透明になっています。

このディールフローは、

投資実行にいたらないまでも、VCが保有するディールフロー・投資候補の企業群リストそのものに価値があったりする。
deal flow, ディールフロー | VCレスキュー(仮)
等々言われていたりして、一子相伝だったり、門外不出だったりします。 実際のところ、ガイドラインはあるもののカチッと字にして標準化しろ, PDCAしろと言われると辛いものがあるせいかもしれません。

投資の対象

スタートアップへの投資の目的は通常Financial Investment(財務リターンの追求)とStrategic Investment(戦略リターンの追求)の2つに分類されております。

事業会社(メーカー等を中心とした投資専業じゃない会社)がスタートアップへ投資する場合は、CVC(Corporate Venture Capital)を組織し、そこから投資することもままあり、 特に日系のCVCは「事業提携を狙いとしたStrategic Investment」を目的として唄い、 財務部門や開発部門から人をローテーションで集めるなどして対応にあたることも少なくないように見えます。

対して、欧米系のCVCはスタートアップ投資への打席に立ち続けることを重視し、「Financial Investment」に比重を置くことが多いように見えます。 CVCのメンバーはVCなど金融系(の投資系)のキャリアが背景にあるメンバーを雇って、本体と独立で動くパターンが多いと思います。

もちろん、Financial Investmentに比重をおいているファンドであっても、金が儲かるならどこでも投資するというわけではなく、 ファンド自身の差別化をはかるためにも、事業部門との何らかの関係が見込める投資テーマを持っています。

例えば、シリコンバレーのベンチャーキャピタルの、5つの秘密と成功の秘訣 | 500 Startups Japanでは、ファンドに対し「1.何かの象徴として認識されるようになれ=ブランドになれ」ということと、 「2. 投資Thesis(ポリシー・理論)を持て=投資基準を明確にせよ」ということを挙げています。

Gitlabでは?

さて、この点、Gitlabは投資(買収先)スタートアップを

  • 社員数10名以下
  • Gitlabの長期開発戦略(3年)に合致すること

としており、

  • シードステージ対象
  • プラダクトの開発加速を狙いとしたAcqui-hiring(Strategic Investment)

を活動のガイドラインとして位置付けています。早い話が、スタートアップ買収を採用活動の一環としてやっているということです。

実施部門はCorporate developmentですね。

ディールフロー

通常スタートアップ投資は、たくさんの会社情報を集めて相対評価しながらふるいをかけながら進めます。東証と違って非公開株式を扱うので相場観の形成は属人的です。普段からたくさん見ておくことが大事なのでしょう。 さまざまな事情で「絶対ここに投資するんだい」と決めてから、一生懸命理屈をつけていく場合もあるようですが…。

ディールフローはどこも似ていてざっくり

  1. ソーシング/テーマ決め
  2. 候補リストの作成/声かけ(提案)
  3. デューデリジェンス; 会社の身体検査
  4. 投資実行

という流れです。こちらも様々な事情によりソーシングの時点でかなり的が絞られているパターンもあるようです。

先日読んだ「経営パワーの危機」では、

リスク投資ばかりを行う米国のベンチャーキャピタルはディールフローと呼んで広い情報網から案件を集め、一〇〇件に一件くらいの割で厳選する。
三枝 匡. 経営パワーの危機 日本経済新聞社. Kindle 版.

とあって、100件はどこの数字かなと思いながら読んでいました。ちなみに、本書では

面白い投資先があまり見つからなくて数少ない投資案件から選ばざるを得ませんでした
三枝 匡. 経営パワーの危機 日本経済新聞社. Kindle 版.

と、投資失敗の風景が描かれています。1の時点で頑なになりすぎるのも良くないようですね。

Gitlabでは?

量が多いので超訳とともに、ざっくり抜粋します。

ディールフロー

これに対してGitlabのディールフローは以下のようです。

  1. Pipeline Building
    • (CrunchbaseのDBなど)マスターリストに掲載されている会社数: 15,000+社
    • プロジェクトで投資検討対象にする会社数: 1,000社
      • 社外データ: 800社
      • 社内情報: 100社
      • Cooperate Development部で作成する候補: 100社
    • 投資検討優先会社: 400社
  2. Exploratory
    • 投資提案会社数: 400社
    • 意思確認・会社の状況確認: 300社
    • 絞り込み(優先づけ): 100社
    • 事業部門同席との絞り込み・投資先への技術的質疑応答: 50社
  3. Early Diligence
    • 投資検討: 30社
  4. Confirmatory Due Diligence
    • 投資条件交渉: 15社
  5. Integration
    • 投資実行: 10社

1. Pipeline Buildingで何をするか

商品企画チームとM&Aを通した強化が必要そうな分野を特定する。活動を通して買収検討リストを作成する。

  1. Crunchbaseなどスタートアップ情報が載っているデータベースからどういった分野がホットか眺めて絞る
  2. Gitlab社社員や知り合いなどにコンタクトがあったスタートアップ情報を集める
  3. Gitlab社の長期ビジョンや戦略に応じて調査をかける

2. Exploratory

お互いの方向性や買収そのものがマッチするかどうかの確認。統合のイメージ感も検討。

3. Early Diligence

買収検討着手。検討期間中はコードネームをつけるらしい。

  1. NDA締結
  2. チームと責任者(リード)の設定
  3. Preliminary diligence
    各種情報の整理
    1. Financials
      • 財務3表
      • 税務状況
    2. Employees
      • 社員名簿: 名前・役職・人気・在籍年数・勤務地・給与・LinkedInプロファイル・プログラミング能力など
      • 社員の履歴書
      • 雇用契約書・知財権譲渡契約
    3. 取引先情報
      • 顧客情報: 名前・月間売上・契約日
      • サプライヤ情報: 月間支出
    4. 資産
      • 買収で取り込む有形・無形資産
      • 買収では除外される有形・無形資産
    5. Tecnical BOM / 技術資産
      ソフトウェア、データ、SNSアカウント等々
  4. Early technical diligence
    コードレビューやOSS利用状況、開発規約等のレビュー
  5. 従業員のレビュー(キーパーソンの特定等)
  6. 給与など待遇のレビュー
  7. 開発状況に関するインタビュー
    対象: キーパーソン
  8. 製品統合のマイルストーン検討
  9. ROIの算定
  10. 企画の作成・社内承認
  11. 買収条件の調整(Term Sheet)

4. Confirmatory Due Diligence

本格的なデューディリジェンス. 以下略.

参考文献

  1. ベンチャーキャピタルの仕組みとディールソーシングの重要性:とむさとうのコラム集
  2. CVC : ベンチャーキャピタリストへの道
  3. Acquisition Process | GitLab

カルマンフィルタの誤差共分散の意味

先日、逐次ベイズ推定とカルマンフィルタにて、ベイズ推定とカルマンフィルタの関係を考察したが、「誤差の共分散」が何を表すのか深く考えないでいた。68–95–99.7則から「共分散すなわち $ \sigma $ の多変量版」だから、68%確率円(気象庁|台風情報の種類と表現方法) くらいに思っていたが、甘かったのでした。

おかげで、カルマンフィルタにおける誤差楕円の計算方法 - MyEnigmaを拝見した時に、意気がったツイートをしてしまったのですが、私のほうが違っていたのでお詫び申し上げます。すみませんでした!!

正規分布と確率

1次元

確率変数Xがスカラーで、1次元の場合を考えます。

\[ X \sim \mathcal{N}(\mu, \sigma^2) \qquad \sigma \in \mathbb{R}^+ \]

このとき、$ X $ が閉区間 $ [\mu -a \sigma, \mu + a \sigma] $ (ただし、$ a \in \mathbb{R}^+ $) に収まる確率は、標準正規分布の累積密度関数 $ \Phi(x) $ を用いて

\[ \text{Pr}(|X - \mu |\leq a \sigma) = \Phi(a) - \Phi(-a) = erf(\frac{a}{\sqrt{2}}) \]
(厳密には半開区間だが、下限の等号の有無は気にしない方が多いようなので…)

あとは標準正規分布表やらerfやらを用いると「$ 3 \sigma $ の範囲に99.7%に収まる」などが求まります。「何が範囲に収まるんですか?」というところについては、 カルマンフィルタの予測ステップなりフィルタステップなりで得られた事後分布 $ \mathcal{N}(\hat{\mu}, \sigma) $ に対して、確信区間という形で「真値が収まる確率」という解釈で良いでしょう。

erfの逆関数を用いる形で任意の確率について、aを求めることができます。

多次元

さて、k次元の場合に移ります。まず、おさらいで、多変量正規分布の確率密度分布は、

\[ f(\mathbf{x}) = \mathcal{N}_k(\textbf{x}|\mathbf{\mu}, \Sigma) =\left(\frac{1}{2\pi}\right)^{k/2}|\Sigma|^{-1/2}\exp\{-\frac{1}{2}(\textbf{x}-\mathbf{\mu})'\Sigma^{-1}(\textbf{x}-\mathbf{\mu})\} \]
のように定義されます。1次元のときと同様に閉区間にあたるものを定義したいですね。 expの中の項に注目すると $ (\textbf{x}-\mathbf{\mu})'\Sigma^{-1}(\textbf{x}-\mathbf{\mu}) $ が目に付きます。 この値が等しいときは確率密度 $ f(\cdot) $ も等しくなるわけですね。 $ \sqrt{(\textbf{x}-\mathbf{\mu})'\Sigma^{-1}(\textbf{x}-\mathbf{\mu})} $ はマハラノビス距離と呼ばれていて、 マハラノビス距離が($ \mathbf{\mu}$から)等距離のベクトルはk次元の楕円の表面上にあります。 ここから、
\[ g(d) = \text{Pr}\{(\textbf{x}-\mathbf{\mu})'\Sigma^{-1}(\textbf{x}-\mathbf{\mu}) \le d\} \]
と、dと確率を結びつける関数g(1次元の例ではerf)を知りたくなります。 さて、これについては証明は参考文献[1]に讓りますが、自由度kのカイ二乗分布を用いて下記の関係があります。
\[ \text{Pr}\{(\textbf{x}-\mathbf{\mu})'\Sigma^{-1}(\textbf{x}-\mathbf{\mu}) \le \chi^2_{k,\alpha}\}=1-\alpha \]
したがって、カイ二乗分布の累積密度関数の逆関数を使えば、k次元の正規分布について確率 $ \alpha $とマハラノビス距離を結び付けられそうですね。 よって、カルマンフィルタで求めた正規分布$ p(x_n \vert Y_n) = \mathcal{N}(x_n \vert \hat{x}_{n|n}, P_{n|n}) , x_n \in \mathbb{R}^k $ などに対しては、 自由度kのカイ二乗分布を使って誤差楕円を求められそうです。 具体的な方法は、こちらの神サイトを御覧ください。カルマンフィルタにおける誤差楕円の計算方法 - MyEnigma

カルマンフィルタ誤差の共分散の意味

というわけで、68–95–99.7則の感覚で意味づけすることを考えると

\[ \text{Pr}\{(x-\hat{x})'P^{-1}(x-\hat{x}) \le n^2 \} \quad n = 1\ or\ 3 \]
などと一般化した上で、 $ P_{t|t} $ (あるいは $ P_{t|t-1} $) の次元に応じて「真値がXX%の確率で収まる確信区間」という意味でした。 XX%の下りは参考に以下を御覧ください(次元の呪いっぽいものも見えます。)

次元   1σ(?) 3σ(?)
1 68.3% 99.7%
2 39.3% 98.9%
3 19.9% 97.1%
4 9.0% 93.9%

というか、ひょっとしてモデルの次元が高まるほど真値が楕円の中心に来る確率は限りなく低くなってくるのか…?

参考

  1. N-DIMENSIONAL CUMULATIVE FUNCTION, AND OTHER USEFUL FACTS ABOUT GAUSSIANS AND NORMAL DENSITIES
  2. 4.3 - Exponent of Multivariate Normal Distribution | STAT 505
  3. カイ2乗分布 - 高精度計算サイト
  4. 高次元空間中の正規分布は超球面状に分布する - Qiita

HyperX Alloy FPS Pro

スイスでは3月中旬にCOVID19の感染爆発に伴うロックダウンが実施されました。私もWFHに入ったのですが、ノートパソコン・スタンドを使った環境となったので 外付けキーボードが必要となり、KingstonのHyperX FPS Proを購入しました。 使いはじめて4ヶ月ほど経ったのでレビューします。

購入の動機

  • スイスで手に入るキーボードのうち、QWERTZなCHキーボードではなく、
    QWERTYなキーボードがほしかった
  • メカニカルキーボードに興味あった
  • キーボードは卓上を占拠しないコンパクトなものが好き
  • 予算100CHF(約1万円以内) / 納期1週間以内

使ってみて

  • 購入の動機は達成。テンキーなしにも慣れたので、さらに小さいキーボードに移行できるかも。
  • 一番スタンダードな赤軸(Cherry)ということですが、やはりペコペコキーボードより音が大きい。慣れてきた結果、底打ち(※)はしなくなりましたが、 静かなオフィスで使うのは同僚に遠慮して躊躇するかも。
    ※ どうやら「キーを必要以上に押し込んでしまい、大きな接触音がすること」を言うらしい
  • 購入した製品に不具合があるのか、たまにキーが押下固定されてしまい、"dddddddd...."などと虚無を量産することがある。redditでも何件か同じような不具合のコンプレインを目にした。
    • 別のキーを押すと開放されるが、その後、そのキーはキーボードにを再接続するまで使えない
    • 再現条件不明。同じキーを連打している時に起きやすい。2週間に1回程度の頻度。
    • 最新ファームウェア(2.1.1.2) はこちらから. アップデートのご利益なし。
      ->HyperX Alloy FPS Pro Firmware Update / Alloy FPS Pro | HyperX
    • しばらくは我慢してましたが、キーの誤動作で我慢ならん自体が発生したので、サポートに連絡したら代替品を出荷してくれた。到着待ちです。
      誤動作の原因が電気信号っぽい気がする。代替品作戦、地球環境に悪いし付録ケーブル含めてノイズとか確認したほうが良いと思うんですがどうですか>kingstonさん
  • 付属のUSBケーブルが硬くて長め。自分でUSB-mini B用意したほうが良いかも。

逐次ベイズ推定とカルマンフィルタ

カルマンフィルタは逐次ベイズフィルタの特殊系として言われますが、式を眺めていても結びつきがさっぱりイメージできなかったので数式いじりをしました。 もっとしっかりした入門は他でたくさん見かけますので、 参考文献として記事のリンクしておきます。

数式操作に関しては須山さんの本を参考にしました。

集合・位相・圏 数学の言葉への最短コース
ベイズ推論による機械学習入門
  • 作者: 須山敦志
  • 出版社/メーカー: 講談社
  • 発売日: 2017/10/21

逐次的ベイズ推定から見たカルマン・フィルタ

セオリーと違ってカルマンフィルタで得たい出力値 状態量 $ \hat{x}_{k|k} $ の位置付けから入ります。理解を簡単にするため、システムとしてD次元の酔歩 $ x_k \in \mathbb{R}^D $と 酔歩する粒子の観測誤差を含む観測値 $ y_{k} \in \mathbb{R}^D $ で考察します。

\[ \begin{array}{rll} x_{n+1} &= x_{n} + w_n & \qquad w_n \sim \mathcal{N}(0, Q) \quad Q \in \mathbb{R}^{D \times D} \\ y_{n} &= x_{n} + v_n & \qquad v_n \sim \mathcal{N}(0, R) \quad R \in \mathbb{R}^{D \times D} \\ \end{array} \]
観測値 $ y_n $ は真の平均値のわからない正規分布 $ \mathcal{N}(x_n, R) $ から取り出されたもので、 このn=1のデータからは、せいぜい $ \hat{x}_n = y_n $ とするのが関の山です。 しかし、うまく $ x_{n+1} - x_{n} \sim \mathcal{N}(0, Q) $ と、 観測したデータ $ Y_n = \{Y_{n-1}, y_{n}\} $ を使って、より高い精度で $ x_n $ を推定することを考えます。

カルマンフィルタの式

システムをざっくり簡略化したので、カルマンフィルタの式もザクッと簡単になります。この式は逐次ベイズ推定から求めた結果とあとから比較するために掲載しました。

予測のステップ
\[ \begin{array}{rl} \hat{x}_{k|k-1} &= \hat{x}_{k-1|k-1} \\ P_{k|k-1} &= P_{k-1|k-1} + Q \\ \end{array} \]
フィルタリングのステップ
\[ \begin{array}{rl} e_k &= y_k - \hat{x}_{k|k-1} \\ S_k &= R + P_{k|k-1} \\ K_k &= P_k S_k^{-1} = P_k (R + P_{k|k-1} )^{-1}\\ \hat{x}_{k|k} &= \hat{x}_{k|k-1} + K_k e_k\\ P_{k|k} &= (I - K_k) P_{k|k-1} \\ \end{array} \]

フィルタステップ

観測データ $ Y_n $から状態量 $ x_n $ の確率密度分布 $ p(x_n|Y_n) $ を求めます。確率密度分布にはベイズの定理から下記の性質があります。赤字では $ y_n $ と $ Y_{n-1} $ の独立性(マルコフ性)を用いました。

\[ \begin{array}{rll} p(x_n \vert Y_n) & = p(x_n | y_n, Y_{n-1}) \\ & \propto \color{red}{p(y_n | x_n, Y_{n-1})} p(x_n | Y_{n-1})\\ & = \color{red}{p(y_n | x_n)} p(x_n | Y_{n-1})\\ \end{array} \]

ここで、事前分布として $ p(x_n | Y_{n-1}) = \mathcal{N}(x_n | \hat{x}_{n \vert n-1}, P_{n \vert n-1}) $ をカルマンフィルタでは仮定する。このとき、$ p(y_n \vert x_n) = \mathcal{N}(y_n | x_n, R) $ であるので、

\[ \begin{array}{rll} \ln \left(p(y_n | x_n) p(x_n | Y_{n-1}) \right) &= \ln \mathcal{N}(y_n | x_n, k) + \ln \mathcal{N}(x_n | \hat{x}_{n|n-1}, P_{n|n-1}) \\ &= - \frac{1}{2} \left\{ x_n^T (\color{green}{R^{-1}+P_{n|n-1}^{-1}}) x_n - 2 x_n^T \color{blue}{(R^{-1} y_n + P_{n|n-1}^{-1} \hat{x}_{n|n-1}}) \right\} + \text{const.} \end{array} \]
この数式の形状により $ p(x_n \vert Y_n) $ がガウス分布であることがわかるので (須山本の手法です)、
\[ p(x_n \vert Y_n) = \mathcal{N}(x_n \vert \hat{x}_{n|n}, P_{n|n}) \\ \]
とおくと
\[ \begin{array}{rl} \ln p(x_n \vert Y_n) &= \left( -\frac{1}{2}(x_n - \hat{x}_{n|n})^T P_{n|n}^{-1} (x_n - \hat{x}_{n|n}) - \ln | P_{n|n}| \right) + \text{const}\\ & = -\frac{1}{2} \left( x_n^T \color{green}{P_{n|n}^{-1}} x_n - 2 x_n^T \color{blue}{P_{n|n}^{-1} \hat{x}_{n|n}} \right) + \text{const.} \end{array} \]
ここから各パラメータの比較をします。すると
\[ \begin{array}{crl} &P^{-1}_{n|n} &= \color{green}{R^{-1} + P_{n|n-1}^{-1}}\\ \Leftrightarrow & P_{n|n} &= \left( 1- P_{n|n-1}(P_{n|n-1}+R)^{-1} \right) P_{n|n-1} \\ && = \left( 1- K_n \right) P_{n|n-1} \\ \end{array} \]
と誤差共分散の更新式が求まります。さらに目当ての推定値が
\[ \begin{array}{crl} & P_{n|n}^{-1} \hat{x}_{n|n} &= \color{blue}{ R^{-1} y_n + P_{n|n-1}^{-1} \hat{x}_{n|n-1} }\\ \Leftrightarrow & \hat{x}_{n|n} &= P_{n|n} \left( R^{-1} y_n + P_{n|n-1}^{-1} \hat{x}_{n|n-1} \right) \\ && = \hat{x}_{n|n-1} + P_{n|n-1}(P_{n|n-1}+R)^{-1} (y_n - \hat{x}_{n|n-1}) \\ && = \hat{x}_{n|n-1} + K_n e_n \\ \text{where} & K_n &= P_{n|n-1}(P_{n|n-1}+R)^{-1} \\ & e_n &= y_n - \hat{x}_{n|n-1} \end{array} \]

と求まり、$ K $が、このシステムにおける最適カルマンゲインになっているなど、カルマンフィルタの更新ステップの式になっています。 $ \hat{x}_{n|n} $は、条件付き確率密度分布 $ p(x_n \vert Y_n) $ に従う、$ x_n $ の条件付き期待値 $ \mathbb{E}(x_n|Y_n) $ ですので、 これ自体が最小分散推定になっています(証明は省略. 参考文献[4]参照)。

予測のステップ

先ほどのフィルタリングのステップでは、事前分布として$ p(x_n | Y_{n-1}) = \mathcal{N}(x_n | \hat{x}_{n \vert n-1}, P_{n \vert n-1}) $を用いていました。この事前分布のパラメータ、$ \hat{x}_{n|n-1} $と$ P_{n|n-1} $を求めます。 $ x_{n+1} - x_{n} \sim \mathcal{N}(0, Q) $ であることから、 その1ステップ前のフィルタリングステップで求めた粒子位置に関する確率密度分布 $ p(x_{n-1}|Y_{n-1}) = \mathcal{N}(x_{n-1}|\hat{x}_{n-1|n-1}, P_{n-1|n-1}) $ に対して、 正規分布の加法性から

\[ \begin{array}{rl} \hat{x}_{n|n-1} &= \hat{x}_{n-1|n-1}\\ P_{n|n-1} &= P_{n-1|n-1} + Q\\ \end{array} \]
となります。そのままですね。

初期値

さて、同じようにしてフィルタリングステップ→予測ステップ→フィルタリングと、ステップをどんどんさかのぼっていくと、最初の位置 $ \hat{x}_{0|0} $, そしてその位置に関する共分散 $ P_{0|0} $ に行き着きます。 これはハイパーパラメータというか、えいやと決めるものになります。

システムがもう少し複雑な場合

前節では酔歩を仮定しましたが、一般の線形システムとして、

\[ \begin{array}{rl} x_{n+1} &= A x_n + B u_n + w_n \\ y_n &= C x_n + v_n \end{array} \]
について考えます。このとき、
\[ \begin{array}{rl} p(x_{n+1}|x_{n}) &= \mathcal{N}(A x_n + B u_n, Q) \\ p(y_n|x_n) &= \mathcal{N}(C x_n, R) \end{array} \]
であって、他の確率分布についても正規分布の再生性をふまえて
\[ \begin{array}{rl} x &\sim \mathcal{N}(\mu, \Sigma) \\ \Rightarrow G x &\sim \mathcal{N}(G \mu, G^T \Sigma G) \end{array} \]
に注意すれば、もとのカルマンフィルタの式になるのは自明ですね。

参考文献

  1. シンプルなモデルとイラストでカルマンフィルタを直観的に理解してみる - Qiita
  2. ベイズの定理からカルマンフィルタを導出する - ssk tech blog
  3. データ同化|カルマンフィルタと尤度 - ari23の研究ノート
  4. 知能制御システム学 画像追跡 (3) ― ベイズ推定とパーティクルフィルタ ― 東北大学 大学院情報科学研究科 鏡 慎吾 准教授
  5. カルマンフィルタの誤差共分散の意味 - 夜間飛行

Numpy使いのためのEigenチートシート

Eigen

v3.4からスライスが便利になるそうだが、最新Stable版の利用を前提にする。

時期が来たらこっそりエントリー書き換えるかも。

git clone https://gitlab.com/libeigen/eigen.git
git checkout tags/3.3.7

導入

#include <Eigen/Core>
using namespace Eigen;

基本

ベクトル - 1D Array - Vector

  • numpyのデフォルトはfloat64(=double)らしいので、それに倣う。
  Numpy Eigen
配列確保 a = np.empty(size) VectorXd a(size);
ゼロベクトル a = np.zeros(size) VectorXd a = VectorXd::Zero(size);
1ベクトル a = np.ones(size) VectorXd a = VectorXd::Ones(size);
fill a.fill(1) a.fill(1);
ベクトル長 a.shape[0] a.size();
要素アクセス a[0]
a[-1]
a(0)
a(a.size()-1)
スライス・アクセス a[a:a+n]
a[:n]
a[-n:]
a.segment(a,n)
a.head(n)
a.tail(n)
要素べき乗 a**n
a**2
a.array().pow(n)
a.array().square()
要素ルート a**0.5 a.array().sqrt()
要素総和 np.sum(a) a.sum()

行列 - 2D Array - Matrix

  Numpy Eigen
配列確保 a = np.empty((r,c)) MatrixXd a(r,c);
ゼロ行列 a = np.zeros((r,c)) MatrixXd a = MatrixXd::Zero(r,c);
1s 行列 a = np.ones((r,c)) MatrixXd a = MatrixXd::Ones(r,c);
サイズ a.shape a.rows(), a.cols()
要素アクセス a[r,c] a(r,c)
スライス・アクセス a[r:r+n,c:c+m] a.block<n,m>(r,c)
a.block(r,c,n,m)
列アクセス a[r] a.row(r)
行アクセス a[:,c] a.col(c)
flatten#1 a.flatten('F') Map<const VectorXd>(a.data(), a.size())
flatten#2 a.flatten() MatrixXd b = a.transpose();
Map<const VectorXd>(b.transpose().data(), b.size());
reshape
- MatrixXdのみ
a.reshape((r,c)) a.resize(r,c)
// 非破壊でreshapeするなら a.reshaped(r,c)
reshape a.reshape((r,c)) MatrixXd x = Map<Matrix<double, r, c>>(a.data());
transepose a.T a.transpose()
逆行列 np.linalg.inv(a) a.inverse()

配列構築・スニペット系

Eigenにはないnumpy/scipyの便利関数たち

  Numpy Eigen
vstack np.vstack([a, b])
MatrixXd c(a.rows()+b.rows(), a.cols());
c << a, b;
hstack np.hstack([a, b])
MatrixXd c(a.rows(), a.cols()+b.cols());
c << a, b;
対角行列 np.diag([0, 1, 2, 3])
MatrixXd v(4,1);
v << 0.,1.,2.,3;
MatrixXd m = v.asDiagonal();
block Diag # Scipy # numpycpp参照

参考

SympyでC++/Eigenな運動方程式の生成

ちょっと複雑な運動方程式をCやPythonコード(numba)に落とすとき、Sympyが便利です。Sympyのコード生成にちょっと癖があったのでメモ代わりに残します。 例として倒立振子をモデルにした。

モデル

欲しい式

非線形な運動方程式に対して、

\[ \dot{q} = f(q, u) \]

テーラー展開で一次線形近似した式の赤いところと青いところを計算してくれる式が欲しい。

\[ \dot{q} = \color{red}{f(q_0, u_0)} + \color{blue}{\frac{\partial f}{\partial q} \Big\vert_{q_0 u_0}} q +\color{blue}{ \frac{\partial f}{\partial u}\Big\vert_{q_0 u_0}} u \]

以降簡単のため、

\[ \begin{array}{rl} A_c(q_0, u_0) &= \frac{\partial f}{\partial q} \Big\vert_{q_0 u_0} \\ B_c(q_0, u_0) &= \frac{\partial f}{\partial u}\Big\vert_{q_0 u_0} \\ g(q_0, u_0) &= f(q_0, u_0) \end{array} \]

運動方程式

一般化座標 $ q $を使って、カートの駆動力 $ u $, カート質量 $ M $, 振り子長 $ l $, 振り子先端質量 $ m $ について運動方程式を立てます。

\[ q = [x, \theta, \dot{x}, \dot{\theta}] \]

あとの導出方法等は他所様の記事を引用して終わります

(7節まで、8節はスキップ. あと回転座標のとり方が若干違う)

Sympyモデルの準備

さて、Qiitaの記事ではラグランジアンを使っていましたが、別ルートで導出したので、こちらになります。

import sympy as sy
class cart_pole():
    def gen_rhe_sympy(self):
        g = sy.symbols('g')
        l = sy.symbols('l')
        M = sy.symbols('M')
        m = sy.symbols('m')

        q  = sy.symbols('q:{0}'.format(4))
        qd = q[2:4]
        u  = sy.symbols('u')
        
        I = sy.Matrix([[1, 0, 0, 0], 
                      [0, 1, 0, 0], 
                      [0, 0, M + m, l*m*sy.cos(q[1])], 
                      [0, 0, l*m*sy.cos(q[1]), l**2*m]])
        f = sy.Matrix([
                      qd[0], 
                      qd[1],
                      l*m*sy.sin(q[1])*qd[1]**2 + u,
                      -g*l*m*sy.sin(q[1])])
        return sy.simplify(I.inv()*f)
    
    def gen_lmodel(self):
        mat = self.gen_rhe_sympy()
        q = sy.symbols('q:{0}'.format(4))
        u = sy.symbols('u')
        
        A = mat.jacobian(q)
        #B = mat.jacobian(u)
        B = mat.diff(u)
        
        return A,B

Python

Pythonの場合はlamdifyを使って

c = CartPole()
q = sy.symbols('q:{0}'.format(4))
u = sy.symbols('u')
# calc_rhe(q, u) -> np.array
calc_rhe = sy.lambdify([q,u], c.gen_rhe_sympy())

などとしたほうが自然ですが、numbaに放り込みたいなどの事情があることもあるでしょう。そこでリファレンス通りにやってみます

from sympy.printing.pycode import pycode
c = CartPole()
pycode(c.gen_rhe_sympy())

すると、あんまりありがたくない形の出力が出ます。

'ImmutableDenseMatrix([[q2], [q3], [((1/2)*g*m*math.sin(2*q1) + l*m*q3**2*math.sin(q1) + u)/(M + m*math.sin(q1)**2)], [-(g*(M + m)*math.sin(q1) + (l*m*q3**2*math.sin(q1) + u)*math.cos(q1))/(l*(M + m*math.sin(q1)**2))]])'
  • ImmutableDenseMatrixというSympyの型そのままの出力
  • sin/cosがmath? numpyがいいな!
  • q2, q3じゃなくてq[2], q[3]と出して欲しい


さらに、リファレンスには使い方が詳しくない裏メニューのNumPyPrinterを使うとこうなります

from sympy.printing.pycode import NumPyPrinter
c = cart_pole()
NumPyPrinter().doprint(c.gen_rhe_sympy())

出力

'numpy.array([[q2], [q3], [((1/2)*g*m*numpy.sin(2*q1) + l*m*q3**2*numpy.sin(q1) + u)/(M + m*numpy.sin(q1)**2)], [-(g*(M + m)*numpy.sin(q1) + (l*m*q3**2*numpy.sin(q1) + u)*numpy.cos(q1))/(l*(M + m*numpy.sin(q1)**2))]])'

かなり良くなりましたが、

  • q2, q3じゃなくてq[2], q[3]と出して欲しい

が残ります。掘ってみるとSympyの立式で

q = sy.symbols('q:{0}'.format(4))

と定義すると、q[0] が q0 というsympyシンボルとして定義されているところから来ている様子です。

Matrixやindexingを駆使するとできるようなことも書いてある気がしないでもないですが、面倒なのでこうしました。

from sympy.printing.pycode import NumPyPrinter
c = cart_pole()
class NumPyPrinterR(NumPyPrinter):
  def _print_Symbol(self, expr):
      name = super(NumPyPrinter, self)._print_Symbol(expr)

      # 変数名が他とかぶらないのでマッチングせずに頭文字だけで判断
      if name[0] == 'q':
          name = name[0] + '[' + name[1] + ']'
      elif name[0] == 'u':
          name = 'u[0]'
      return name
NumPyPrinterR().doprint(np.squeeze(c.gen_rhe_sympy()))

これでめでたし。

C++11

Pythonで検証したあとは同じようにするだけです。C++ではEigenを使いたいのとImmutableDenseMatrixの処理が定義されていないのでもう一捻りします。

別途

Eigen::Matrix<double, 4, 1> g;
Eigen::Matrix<double, 4, 4> A;
Eigen::Matrix<double, 4, 1> B;

と定義されていることとしましょう

c = cart_pole()
class CXX11CodePrinterR(CXX11CodePrinter):
    def _print_Symbol(self, expr):
        name = super(CXX11CodePrinter, self)._print_Symbol(expr)
        
        if name[0] == 'q':
            name = name[0] + '(' + name[1] + ', 0)'
        elif name[0] == 'u':
            name = 'u(0, 0)'
        return name

print('g << ')
for expr in np.squeeze(c.gen_rhe_sympy()).tolist():
    print(CXX11CodePrinterR().doprint(expr), end=',\n')

A, B = np.squeeze(c.gen_lmodel())

print('\nA << ', end='')
for expr in A:
    print(CXX11CodePrinterR().doprint(expr), end=',\n')

print('\nB << ', end='')
for expr in B:
    print(CXX11CodePrinterR().doprint(expr), end=',\n')

あとはセミコロンとか処理してお終い。

その他

Sympyには他にも

  • julia
  • Rust

などのプリンタも用意されているそうです!

SOCPのソルバーでQPを解く話

PythonRoboticsの一例でQPのソルバーとしてECOSが使われていました。ECOSは二次錐計画問題 (SOCP)のソルバーです。 QP ⊂ QCQP(Quadratically constrained quadratic program)⊂ SOCP であって、
制約なしのQP形式の問題

\[ \begin{array}{ll} \min_{x} & \frac{1}{2} x^T Q x + q^T x + r \\ \end{array} \]

は、変数tを追加して

\[ \begin{array}{ll} \min_{x,t} & t \\ {\rm s.t.} & \frac{1}{2} x^T Q x + q^T x + r \leq t \\ \end{array} \]

と変換すればSOCPとして解ける(制約ありも同様)とのことでしたので、ちょっと実装してみます。

まず改めて、$ x \in \mathbb{R}^n $ のとき、

\[ \begin{aligned} y &\equiv \left[ x, t\right] \\ c &\equiv [0_n,1] \\ Q &= U^T U \\ Q_s &\equiv \begin{bmatrix} U & 0 \end{bmatrix} \\ q_s &\equiv [q,0] \\ \end{aligned} \]

とすると表現がこう変わります

\[ \begin{array}{ll} \min_{y} & c^T \cdot y \\ {\rm s.t.} & \frac{1}{2} y^T Q_s^T Q_s y + q_s^T y + r \leq c^T \cdot y \\ \end{array} \]

制約式をノルムの形式に変換すると

\[ \left\| \begin{matrix} \left( 1 + (q_s+c)^{T}y+r \right) / \sqrt{2} \\ Q_s y \end{matrix} \right\|_{2} \leq \frac{1-(q_s+c)^{T}y-r}{\sqrt{2}} \]

最後にCVXOPTなどでも使われるスラック変数を用いた表現に変形します。

\[ \begin{bmatrix} - (q_s + c)/\sqrt{2} \\ (q_s + c)/\sqrt{2} \\ Q_s \end{bmatrix} \cdot y + \begin{bmatrix} (1-r)/\sqrt{2} \\ (1+r)/\sqrt{2} \\ 0 \end{bmatrix} = \left[ \begin{array}{c} s_0 \\ s_1 \\s_2 \end{array} \right], s_0 \geq \left\Vert \begin{array}{c} s_1 \\ s_2 \end{array} \right\Vert_2, s_0, s_1 \in \mathbb{R}, s_2 \in \mathbb{R}^{n} \]

続けて $ G \cdot x \preceq_K h $ という形式にあわせて、

\[ \begin{bmatrix} (q_s + c)/\sqrt{2} \\ - (q_s + c)/\sqrt{2} \\ - Q_s \end{bmatrix} \cdot y + \begin{bmatrix} s_0 \\ s_1 \\ s_2 \end{bmatrix} = \begin{bmatrix} (1-r)/\sqrt{2} \\ (1+r)/\sqrt{2} \\ 0 \end{bmatrix} \]

QPをSOCPとして解く

くどいですが、以上をまとめると

\[ \begin{array}{ll} \min_{x} & \frac{1}{2} x^T Q x + q^T x + r \\ where & x \in \mathbb{R}^n \end{array} \]
というQPは、
\[ \begin{array}{ll} \min_{y} & c^T \cdot y \\ \text{s.t.} & G \cdot y \preceq_{\mathcal{K}} h \quad \left( \Leftrightarrow G \cdot y + s = h, s \in \mathcal{K} \right) \\ \text{where} & y \equiv [x\ t], t \in \mathbb{R} \\ & c = [0_{1 \times n}\ 1] \\ & Q = U^T U \\ & \mathcal{K} \equiv \{ (u_0, u_1) \in \mathbb{R} \times \mathbb{R}^{n+1} \ |\ u_0 \geq \Vert u_1 \Vert_2 \} \\ & G = \left[ \begin{array}{cc} q/\sqrt{2} & 1/\sqrt{2} \\ -q/\sqrt{2} & -1/\sqrt{2} \\ U & 0 \end{array} \right], h = \left[ \begin{array}{cc} (1-r)/\sqrt{2} \\ (1+r)/\sqrt{2} \\ 0_{n \times 1} \end{array} \right] \\ \end{array} \]
というSOCPとして解くことができます。

余談

ECOSに投げるときは、Usage from C · embotech/ecos Wikiを参照しつつ

// Nは元となるQPのxの次元
idxint l = 0; // もし制約付QPを解く場合は l>0, 配列Gに対してはconeの定義にあたる係数よりも前方に配置する
idxint ncones = 1;
idxint q[] = {N+2};

とします。

補足

$ Q = U^T U $ は、コレスキー分解です。

コレスキー分解でベンチマークを取った先日のポストも広告:Pythonと関数単位での高速化(cython, numba, swig) - 夜間飛行

参考

脱PDFからのシングルHTMLのすすめ

レポートやメモの作成に定番のMicrosoft Wordや同時編集が魅力的なgoogle documentsを使ってきたが、gitで管理できるフォーマットを試しながら、gulp+pugで作るシングルページHTMLに行き着いたのでメモする。

これまで試したもの

やりたかったこと

脱ペーパー
これまでPDF化するなど印刷前提のフォーマットを貫いてきたけどペーパーレスの時代に「紙ページに収める」という縛られ方をする必要はないのではないだろうか。 挿絵の位置で改ページがあっちいったり、こっち行ったりをどうこうしたり、 横に平べったく長い表の貼り付けたり、 章タイトルの直後に改ページが来てpandocフィルタ追加したり… そういう調節とお別れしたい。
メール添付やクラウドストレージを通した配布
クラウドストレージのプレビュー機能の充実とブラウザ間の互換性向上によって、シングルHTMLファイルがPDFの代替になりそうな気がした。
インタラクティブなコンテンツ
いったん脱PDFし、HTMLをフォーマットとして選ぶとJSとかSVGとかの採用が見えてくる。

環境の準備

Webpackも考えたけどentry pointがjsファイルになる点が今ひとつコンセプトに合わなくてgulpにした。

npm install --global gulp-cli
npm init -y
npm i -D gulp pug gulp-pug
npm i -D gulp-inline-source gulp-inline-images

構成

root
├── _template
├── dest <-- ここに出来上がったレポートが出力される
├── gulpfile.js
├── package-lock.json
├── package.json
├── img <-- レポートに使う画像を放り込んでおく
└── src <-- ここにレポートの素を置く

gulpfile

const gulp = require('gulp')
const pug = require('gulp-pug')
const inlineSource = require('gulp-inline-source')
const inlineImages = require('gulp-inline-images')

gulp.task('pug', () => {
  return gulp.src(['./src/*.pug', '!./_**/*', '!./node_modules/**/*'])
    .pipe(pug({ pretty: true }))
    .pipe(inlineSource())
    .pipe(inlineImages({ basedir: './img/' }))
    .pipe(gulp.dest('./dest/'))
})
gulp.task('watch', () => {
  return gulp.watch(['./src/*.pug', './_template/*.pug'], gulp.series(['pug']))
})
gulp.task('default', () =>  {
  return gulp.parallel(['pug'])
})

あとはgulpでもgulp watchでも…. 画像, CSSそれにjavascriptはinline化してhtmlに埋め込んでしまいます。

レポートのテンプレート

コンセプトにも関わらず、紙に印刷したいとなったらpaper.css、プレゼンしたいならreveal.jsと思いつつ、Bootstrapをcssのみで使う選択肢で作ってみた。

ライブラリを追加する場合inlineのアノテーションを忘れないように!!

doctype html
html(lang="ja")
  block config
  head
    meta(charset="UTF-8")
    meta(name="viewport", content="width=device-width, initial-scale=1")
    link(rel="stylesheet", href="https://stackpath.bootstrapcdn.com/bootstrap/4.5.0/css/bootstrap.min.css", inline)
    style.
      body {
        font-family: "Helvetica Neue", 
              "Hiragino Kaku Gothic ProN", Meiryo, sans-serif;
      }
    block css
    block scripts

    title #{title} 
body
  .container
    block content

本文の書き方

カラム数などによっては40emで折り返すようにしておくと読みやすいかも。

extends ../_template/template.pug

block config
  - var title='私のレポート'

block css
  style.
    .alert-trim {
      display: inline-block;
    }
    .col p {
      width: 40em;
    }

block content
  .row
    .col
      h1 報告
      h2 はじめに
      p 今月は…

Boilerplate

ボイラープレート用意した

Chachay/pug_boilerplate: My Boilerplate of Single HTML reports

Pythonと関数単位での高速化(cython, numba, swig)

Python、計算が遅いです。javascriptよりも遅いです。

しかし、Python as グルー言語として「データの読み込み」や「計算処理」、「可視化」を橋渡しする機能に注目すると本当に便利なプログラミング言語だと思います。 橋渡し先の「計算処理」の定番としてはnumpyがあったり、scipyがあったりしますが、 ここをオリジナルの計算処理に置き換える手段として、 pybind11, swig, numbaやcythonがあると思っています。 使い勝手や速度の面で気になっていたので比較しました。

環境

  • Windows 10
  • Visual Studio 2019
  • miniconda - python 3.7
    • swig
    • numba
    • cython

やったこと

LU分解がベンチマークでは良く使われるそうですが、実装を簡単にするため、コレスキー分解を選びました。 ベンチマークではランダムに生成した実対称正方行列(の正定行列)を使っています。

scipyの例をあげると

# 128 x 128の行列を生成(numpy)
A128 =  gen_sym_matrix(128)

from scipy.linalg import cholesky
# コレスキー分解の時間を計測
%timeit l = cholesky(A128, lower=True)

結果

リポジトリ

python_bench/cholesky at master · Chachay/python_bench

実行時間

行列のサイズで2水準計測。

python_bench/main.ipynb at master · Chachay/python_bench

  方法 128x128
[msec]
256x256
[msec]
コメント
1 scipy 0.09 0.37 内部でLAPACKを呼び出し. 参考値.
2 Swig - CPP - Eigen 0.17 0.91 ホットスポットをCPPで書きたい場合、バランスが良い
インターフェース・ファイル(.i)さえ用意すれば、cppっぽく書ける。 コンパイラの差し替え(Clangとか)で更に性能が伸ばせるかも。
3 Numba - numpy 0.35 1.79 ほとんど、そのままpythonなのに魔術的な力を感じる。Pythonでプロトタイプした直後に処理速度が気になったときの第一候補.
4 Swig - CPP 0.25 2.27 スタンダードC++のみ。もっと速くできるかも。
5 Cython - Numpy 11.6 45.8 お前にはいつもがっかりしている. 俺のことが嫌いなのか。
6 Python - Numpy 21.8 87 いわゆる標準的なアプローチ
7 javascript - 295 飛び入り参加。Pythonとは関係ない。
8 Pure Python 44.4 327 List in Listで行列を表現し、組み込み関数のみで作成。
List in List形式はnumbaやcythonで上手に取り扱えないので、高速化はなし。

実行メモ

Numba

Pythonのnumpy版関数にjitデコレータをつけただけ。numpyの定義でdtypeなどで型を指定してやる必要がある。引数やローカル変数は型推論でうまく処理したのかな。

from numba import jit
@jit(nopython=True)
def cholesky_numpy_numba(A):
    M = A.shape[0]
    L = np.zeros_like(A, dtype=np.float64)

    for i in range(M):
        for j in range(0, i+1):
            s = np.dot(L[i,:j],L[j,:j])
            L[i, j] = np.sqrt(A[i,i]-s) if(i==j) else ((A[i,j]-s)/L[j,j])
    return L

Swig - Numpy

公式資料: numpy.i: a SWIG Interface File for NumPy — NumPy v1.18 Manualにいろいろ詳しい。闇についても若干言及がある。

ファイル構成

swig
├── numpy.i
├── setup.py
├── swig_mod.cpp
├── swig_mod.h
└── swig_mod.i

setup.py

ビルドの条件を書きます。モジュール名の頭にアンダースコア(_)をつけないとリンカが悲鳴をあげます。これは超重要です。

コンパイラの最適化オプションとかコンパイラの差し替えとかもここで定義できた気がします。思い出したら追記します。

SwigMod = Extension("_swig_eigen_mod", # ここのモジュール名にアンダースコア
                  ["swig_mod.i","swig_mod.cpp"],
                  swig_opts=['-c++', '-py3'],
                  include_dirs = [numpy_include]
                  )

swig_mod.cpp & swig_mod.h

本体。Swigのマクロで2次元配列のnumpyの配列を(行サイズ、列サイズ、ポインタ)に置き換えるので、

import numpy as np
from SwigMod import cholesky_swig
A = np.zeros((3,3))
L = cholesky_swig(A)

のような使い方を予定しているモジュールを作るときには、引数をint mm1, int mn1, double* mat1のように定義しておき、 戻り値も引数としてdouble** outmat, int* mm, int* mnのように定義しておく。

// 関数の定義としては戻り値はvoid
void cholesky_swig(
          // 引数Aの成れの果て
          int mm1, int mn1, double* mat1,
          // 返り値Lの素
          double** outmat, int* mm, int* mn){
  // メモリの確保はC++で行い、そのまま管理をpythonに引き継ぐ
  double* arr = (double*)calloc(mn1 * mm1, sizeof(double));

  // ---- 中略 ------

  // 返り値はポインタで渡す。returnしない。
  *mm = mm1;
  *mn = mn1;
  *outmat = arr;
}

戻り値についてはSwigがラッパーを作るときに、引数を戻り値に置き換えるらしい。黒魔術なので深入りしてない。公式資料に詳しい。numpyの1次元配列、多次元配列についても同様。

後述のSwig - Eigenのようにnumpyをtypedefして、きちんとインターフェースを書くとこんな気持ち悪いことにはならないのかもしれない。

callocで確保したヒープが最後どのようになるのか杳として知れず。

swig_mod.i

ヘッダーファイルとpythonをつなぐラッパーを生成するためのスクリプト。numpyを開くところについては下記でマクロ定義。変数名はヘッダーで使ったものと一致するように。

%apply (double** ARGOUTVIEWM_ARRAY2, int* DIM1, int* DIM2){(double** outmat, int* mm, int* mn)}
%apply (int DIM1, int DIM2, double* IN_ARRAY2){(int mm1, int mn1, double* mat1)}
// 複数定義する場合
// %apply (int DIM1, int DIM2, double* IN_ARRAY2){(int mm1, int mn1, double* mat1), (int mm2, int mn2, double* mat2)}

モジュールのビルド

python setup.py build_ext --inplace

するとsetup.pyと同じ階層にモジュールが作られるので、jupyterから呼び出しやすいところへコピーして使う。

Swig - Eigen

ありがたいベストプラクティスに出会えたので追加。

eigen.iが秘伝のタレです。numpyとeigenを橋渡ししてくれます。46行目の#include "Eigen/Core"だけ環境によっては書き換え必要。

インターフェースファイルの書き方も詳しいので熟読するべき

swig-eigen-numpy/inverter_wrapper.i

eigen.iがしっかりしているので、hppもこれがpythonモジュールの定義ファイルなのか疑うレベルできれいなC++になります。

#include <exception>
#include <vector>
#include "eigen/Eigen/Core"

Eigen::MatrixXd cholesky_swig_eigen(const Eigen::MatrixXd &M);

Cython

お前にはがっかりしている。どうしたら良いか私も知りたい。

%%cython
import numpy as np
cimport numpy as np
cimport cython

ctypedef double DTYPE_t

@cython.boundscheck(False)
@cython.wraparound(False)
def cholesky_numpy_cython(np.ndarray[DTYPE_t, ndim=2]  A):
    cdef int i
    cdef int j
    cdef double s
    cdef int M = A.shape[0]
    cdef np.ndarray[DTYPE_t, ndim=2] L = np.zeros_like(A, dtype=np.float64)

    for i in range(M):
        for j in range(0, i+1):
            s = np.dot(L[i,:j],L[j,:j])
            L[i, j] = (A[i,i]-s)**0.5 if(i==j) else ((A[i,j]-s)/L[j,j])
    return L

参考文献

WindowsにSWIG導入 - Volo di notte
いまやconda install swigで楽にインストールできるのですが、anacondaを使いたくない場合などに参照ください
Write C++ extensions for Python - Visual Studio | Microsoft Docs
今後よみたい
SWIG and Python
ちょっと古い話もある

d3.js v5でコンター図を描く

d3.js v5でコンター図を描く

以下のコードでobjcet_functionでコンター図を書きたい対象の関数を定義しています。

x,yはlinspace関数を別途定義してすべての定義域をベクトル化した関数(vect_obj_function)で計算しています。

const margin = {
  top: 20,
  right: 30,
  bottom: 30,
  left: 40
}
const width = 700 - margin.left - margin.right
const height = 400 - margin.top - margin.bottom

const linspace = (start, stop, step = 1) => Array(Math.ceil((stop - start) / step)).fill(start).map((v, i) => v + i * step)
const meshgrid = (x, y) => y.map((y_) => x.map((x_) => [x_, y_])).flat()

// 描画したい関数
const object_function = (x, y) => (x ** 2 + y - 11) ** 2 + (x + y ** 2 - 7) ** 2
const vect_obj_function = x => x.map((v) => object_function(v[0], v[1]))

// 描画したい定義域 linspace(min, max, 刻み値)で変数を作ります
// d3.rangeと機能がかぶってしまっていますが独自定義しています
const x = linspace(-6, 6, 0.1)
const y = linspace(-6, 6, 0.1)

const grid = meshgrid(x, y)

const z = vect_obj_function(grid)

// '等高線'の定義. 等間隔じゃなくて2乗則で決めています
const thresholds = [0.001, ...d3.range(0, 10).map(i => Math.pow(2, i))]

const xScale = d3.scaleLinear(d3.extent(x), [0, width]).nice()
const yScale = d3.scaleLinear(d3.extent(y), [height, 0]).nice()

const transform = ({ type, value, coordinates }) => {
  return {
    type, value,
    coordinates: coordinates.map(rings => {
      return rings.map(points => {
        return points.map(([x, y]) => ([ xScale(x), yScale(y) ]))
      })
    })
  }
}

const contours = d3.contours()
  .size([x.length, y.length])
  .thresholds(thresholds)(z)
  .map(transform)

const color = d3.scaleSequentialLog([d3.max(thresholds), d3.min(thresholds)], d3.interpolateGreys)
const div = d3.select("body")
  .append("div")
  .attr("class", "contour_tooltip")
  .style("font-size", "12px")
  .style("position", "absolute")
  .style("text-align", "center")
  .style("width", "128px")
  .style("height", "34px")
  .style("background", "#333")
  .style("color", "#ddd")
  .style("padding", "0px")
  .style("border", "0px")
  .style("border-radius", "8px")
  .style("opacity", "0");

const svg = d3.select("#fig01")
  .attr("width", width + margin.left + margin.right)
  .attr("height", height + margin.top + margin.bottom)
  .append("g")
  .attr("transform",
    `translate(${margin.left},${margin.top})`)


svg.append("g")
  .attr("fill", "none")
  .attr("stroke", "#fff")
  .selectAll("path")
  .data(contours)
  .join("path")
  .attr("fill", d => color(d.value))
  .attr("stroke", 'white')
  .attr("d", d3.geoPath().projection(d3.geoIdentity()
    .fitSize([width, height], contours[0])))
  .style("stroke-width", 2)

d3.select("#fig01").on("mouseover", () => {
    div.transition().duration(400).style("opacity", 0.9);
    div.style("z-index", "");
  })
  .on("mousemove", () => {
    const point = d3.mouse(d3.event.target)
    const x = xScale.invert(point[0])
    const y = yScale.invert(point[1])
    const z = object_function(x, y)

    div.html(`x = ${x.toFixed(2)} y = ${y.toFixed(2)}<br>f(x,y) = ${z.toFixed(2)}`)
      .style("left", (d3.event.pageX + 20) + "px")
      .style("top", (d3.event.pageY - 35) + "px");
  })
  .on("mouseout", () => {
    div.transition()
      .duration(500)
      .style("opacity", 0);
  })

svg.append("g")
  .attr("transform", `translate(0,${height})`)
  .call(d3.axisBottom(xScale))
svg.append("g")
  .call(d3.axisLeft(yScale))

参考文献

Vue CLIでjavascript/cssをひとつのhtmlに埋め込む

クラウドストレージなどでHTMLプレビュー機能が強化され、Vueで作ったSingle Page Applicationもホストしたりできそうという期待から、 JavascriptやCSSの全てをひとつのHTMLファイルに埋め込む方法を試しました。

環境

VueのVersion等

いわゆるデフォルトのプロジェクトを使います。

html-webbpack-inline-source-pluginの導入

javascriptとcssの展開・埋め込みにhtml-webpack-inline-source-plugin - npmを使います。

yarn add --dev html-webpack-inline-source-plugin

Buildの設定(vue.config.js)

プロジェクトのルートにvue.config.jsを作成します。

// vue.config.js
module.exports = {
  css: {
    extract: false
  },
  chainWebpack: config => {
    config.plugin('inline-source')
      .use(require('html-webpack-inline-source-plugin'))

    config.module
      .rule('images')
      .use('url-loader')
      .loader('url-loader')
      .tap(options => Object.assign(options, { limit: 10240 }))

    config
      .plugin('html')
      .tap(args => {
        args[0].inlineSource = '.(js|css)$'
        return args
      })
  }
}

あとはいつもどおり

yarn build

で/dist以下にindex.htmlが作成されます。

参考

[javascript] Array of ObjectsとObject of Arraysの相互変換

JavascriptでArray of ObjectsとObject of Arraysの変換をたまに使うのでメモ。

はじめに

Array of Objects

const array_objects = [
  { age: 40, first_name: 'Dickerson', last_name: 'Macdonald' },
  { age: 21, first_name: 'Larsen', last_name: 'Shaw' },
  { age: 89, first_name: 'Geneva', last_name: 'Wilson' },
  { age: 38, first_name: 'Jami', last_name: 'Carney' }
]

Object of Arrays

const object_arrays = {
  age: Array [40, 21, 89, 38], 
  first_name: Array ["Dickerson", "Larsen", "Geneva", "Jami"],
  last_name: Array ["Macdonald", "Shaw", "Wilson", "Carney"] 
}

Array of Objects → Object of Array

使いみち

  • グラフ描画系ライブラリの利用時

コード

const new_object_arrays = Object.assign(
  ...Object
    .keys(array_objects[0])
    .map( key => ({ [key]: array_objects.map( o => o[key] ) }))
  );

Object of Array → Array of Objects

使いみち

コード

const new_array_object = Object.keys(object_arrays).reduce(
  (r, key) =>{
    object_arrays[key].forEach(
      (value, i) => {
        r[i] = r[i] || {}
        r[i][key] = value
      })
      return r
  }, [])
または
const new_array_objects = object_arrays[Object.keys(object_arrays)[0]].reduce(
  (r, _, i) =>{
    r.push(
       Object.assign(...Object
      .keys(object_arrays)
      .map( key => ({ [key]: object_arrays[key][i]})))
    )
    return r
  },[])

Chrome拡張の'Create Link'設定メモ

Chromeでさっとリンクを作りたいときに使っているExtention Create Link の設定メモ

ソースコード@githubはku/CreateLink: Make Link alternative to chrome

Formats

いつも使いの基本のフォーマット。markdownはデフォルトで定義されているけど、特にエスケープすると逆に困る場面に遭遇しているので置き換えている。

markdown
[%title%](%url%)
Restructured Text(reST)
`%title% <%url%>`_
pug (jade)
a(href="%url%") %title%

Filter

Amazonのリンクを簡略化したいときなんかにFilter Columnに追加する内容。

正規表現に自信がなければ、RegExr: Learn, Build, & Test RegExの利用がおすすめ。

Filterの文法が下記で定義されているのでスラッシュ(Forward Slash)はHEXでエスケープ(\x2F)する必要がある。

var m = def.filter.match(/^s\/(.+?)\/(.*?)\/(\w*)$/);
Amazonリンク クリーナー(For HTML, Pug, Markdown)
s/(https.*amazon\.co\.jp\x2F).*(?:d|dp|ASIN|product)\x2F(\w{10})(?:\x2F|\?).*?(?="|>|\))/$1dp/$2//g
Amazonリンク クリーナー(For Direct Text)
s/(https.*amazon\.co\.jp\x2F).*(?:d|dp|ASIN|product)\x2F(\w{10}).*/$1dp/$2//g

Power BIでkintoneのレコードを取り出す(REST API)

最近Kintoneを使い始めデータが溜まってきたので分析をしたくなりました。Power BI Desktopの標準機能「Power Query M言語」を使って、Kintoneからデータを取り出します。

用意するもの

  • URI: https://(サブドメイン名).cybozu.com/k/v1/records.json
  • アプリID
  • APIトークン
  • PowerBIに読み込みたい変数名

クエリ編集画面を開く

PowerBI Desktopを起動後、「データの取得」 から「空のクエリ」を開きます。

「クエリエディタ」が開いたら「詳細エディタ」へ。以降が詳細エディタに貼り付ける内容です。

スクリプト

丸括弧で書かれているところは書き直してください。COLUMNSは好きなだけ並べられます。取り出したい変数名のリストを作るのが大変なときは、Excel 2016 や Power BI から kintone のレコードを REST API で取り出すで紹介されている方法でColumn1を展開して、ソースコードで1つ目の'Table.ExpandRecordColumn'をコピペすると楽だと思います。

let 
  // Author: Chachay https://blog.chachay.org/2020/04/kintone-power-bi.html
  BaseUrl         = "https://(サブドメイン名).cybozu.com/k/v1/records.json?",
  APPID           = "app=(アプリID)",
  Token           = "(APIトークン)",
  COLUMNS         = {"(取り出したい変数1)", "(取り出したい変数2)", "(取り出したい変数3)"},
  EntitiesPerPage = 500,

  GetJson = (Url) =>
      let Options = [Headers=[ #"X-Cybozu-API-Token" =  Token ]],
          RawData = Web.Contents(Url, Options),
          Source  = Json.Document(RawData),
          records = Source[records]
      in  records,

  Table.GenerateByPage = (getNextPage as function) as table =>
      let        
          listOfPages = List.Generate(
              () => getNextPage(null),
              (lastPage) => lastPage <> null,
              (lastPage) => getNextPage(lastPage)
          ),
          tableOfPages = Table.FromList(listOfPages, Splitter.SplitByNothing(), {"Column1"})
      in
          Table.ExpandTableColumn(tableOfPages, "Column1", COLUMNS),

  CheckNextPage = (response, Index) as number => 
      if (response=EntitiesPerPage) then Index+1 else -1,

  GetPage = (Index) =>
      let Skip  = "offset " & Text.From(Index * EntitiesPerPage),
          Limit = "limit " & Text.From(EntitiesPerPage),
          Url   = BaseUrl & APPID & "&query=" & Limit & " " & Skip,
          Json  = GetJson(Url),
          data  = Table.FromRecords(Json),
          next  = CheckNextPage(Table.RowCount(data), Index)
      in  
          data meta [Next = next],
  
  GetAllPages = () as table =>
      Table.GenerateByPage((previous) => 
          let
              next_index = if (previous = null) then 0 else Value.Metadata(previous)[Next]?,
              page = if (next_index > -1) then GetPage(next_index) else null
          in
              page
      ),

  ExpandRecordColumns = (Table, columns, fieldName) =>
      if List.Count(columns)=0 then Table
      else 
          let tmp    = Table.ExpandRecordColumn(Table, List.First(columns), {fieldName}, {List.First(columns)}),
              Result = @ExpandRecordColumns(tmp, List.Skip(columns,1), fieldName)
          in Result,

  Table       = GetAllPages(),
  Data        = ExpandRecordColumns(Table, COLUMNS, "value")
in
  Data  

類似手法

Power BI でのデータ取得からビジュアライズまで - Qiita
KintoneからCSVエクスポートを使った方法とKintone ODBC Driverを使った方法の紹介。後述のCDATAさんのプラグイン利用‥?
Excel 2016 や Power BI から kintone のレコードを REST API で取り出す
Power Queryを使ってKintoneのREST APIを使って取り出す方法の紹介。取り出し件数は500件まで。
Power BI から R スクリプトを介して kintone データにアクセスする - R3 Cloud Journey
Visual Studioを導入した上で、Rスクリプトを利用する方法の紹介。500件以上のデータも取り出せるが、ソフトウェアのインストール規模が大きくなるのが難点。
kintone Power BI Connector - Power BI | Connectors | Direct Query
Kintoneのデータを吸い出してPower BIにわたすプラグイン CDATAさん。ローカルからAPIをたたいてPower BI Desktopに読み出すソフトと クラウド上にデータマートを構築してKintoneからデータを複製するサービスの二段構え。

参考文献

レコードの取得(GET) – cybozu developer network
Kintone REST APIのリファレンス。取り出す内容に一工夫したいときはGetPage関数の中でqueryを工夫しましょう。
Power Query M 言語仕様 - PowerQuery M | Microsoft Docs
Power Query M言語の仕様
TripPin 5 - Paging | Microsoft Docs
複数ページにまたがるJSONの取得例. Pagenation × PowerQuery.
How To Do Pagination In Power Query - Mark Tiedemann - Medium
クエリ設計の参考にした