皆様、在宅ワーク万歳な日々を過ごしでしょうか？
先週は久々に現場に駆り出されましたが、
電車内はコロナ禍なんて嘘のような様子でした

コロナとは…

セキュリティ需要の高まり

それはさておき、在宅の増えた昨今では端末のセキュリティについて
考えることも増えたのではないかと思います。

最近だとEMOTETなどのようなマルウェアなどのニュースを見聞きして
怖いなと警戒心が一時的に高まったりすることもあるかとおもいます。

でも日時が立つとそんな警戒心なども薄れて、変なサイトにアクセスしてしまったり、
行きつけのサイトに変なものを埋め込まれたりする可能性も
無きにしも非ずです。

そんなマルウェアが私たちの端末で悪さをするには
自動起動登録する必要があります。
なので、日々自動起動されるプログラムを確認することで、
うちに忍び込む不審な輩に気づくきっかけになったら…！

というわけで「autorunsc.exe」という便利なツールを使用して、
自動起動プログラムの際をログ収集するスクリプトを書いてみました。

import部分

import subprocess
import zipfile
import os
import sys
import urllib.request as req
import pandas as pd
from glob import glob
from plyer import notification
from alittleuseful import loglotate

# pip install pandas
# pip install plyer
# pip install git+https://github.com/ardnico/main

alittleusefulは私が個人的にgithubで公開してる、
ログを書き込むライブラリです
他にも変な機能を公開していますが、よければ使ってくれたら
嬉しいかなって。。。

命名定義

csv_file = f'{os.getcwd()}\\out.csv'
rcsv_file = f'{os.getcwd()}\\out_old.csv'
enc = "utf-16"
URL = "https://download.sysinternals.com/files/Autoruns.zip"
zip_file = "A.zip"
path='.'
logger = loglotate(
    logname = "StartUpSec",
    outputdir = [os.getcwd()],
    lsize = 100000,
    num = 20,
    timestanp = 1 # 1:on other:off
)

静的命名や関数の呼び出し部分です

関数：autorunsc.exeのダウンロード

def download_tool(tf:bool):
    # file download
    if tf == False:
        logger.write('[INFO]Because the tool has not existed, the one will download')
        req.urlretrieve(URL,zip_file)
        with zipfile.ZipFile(zip_file, 'r') as z_file:
            try:
                z_file.extractall(path=path)
                logger.write("[SUCCESS]Tool download succeeded")
            except Exception as e:
                logger.write('[ERROR]Failed to download or unzip autorunsc.exe')
                logger.write(f'[ERROR]{e}')
                sys.exit(0)

「autorunsc.exe」が存在しなかった場合、requestを使用して
ツールをダウンロードします。
ダウンロードがZIP解凍まで行います。

関数：「autorunsc.exe」の実行から自動起動プログラムの比較まで

def get_log():
    if os.path.exists(rcsv_file) == True:
        try:
            os.remove(rcsv_file)
        except Exception as e:
            logger.write('[ERROR]Failed to remove oldcsvfile')
            logger.write(f'[ERROR]{e}')
            sys.exit(0)
    if os.path.exists(csv_file) == True:
        try:
            os.rename(csv_file,rcsv_file)
            df_old = pd.read_csv(rcsv_file,encoding=enc)
        except Exception as e:
            logger.write('[ERROR]Failed to rename oldcsvfile')
            logger.write(f'[ERROR]{e}')
            sys.exit(0)
    else:
        df_old = ''
    with open(csv_file, mode='w', encoding=enc) as fp:
        cp = subprocess.run([f'{os.getcwd()}\\autorunsc.exe','-nobanner','-c','-a','*'], encoding=enc, stdout=fp)
    try:
        if df_old == '':
            flag = 0
        else:
            flag = 2
    except:
        if len(df_old.index) <= 0:
            flag = 0
        else:
            flag = 1
    if flag==0 or flag==2:
        logger.write("[INFO]StartUp Program's log has created")
    else:
        with open(csv_file,encoding=enc) as f:
            data = f.read().split('\n')
        with open(rcsv_file,encoding=enc) as f:
            data2 = f.read().split('\n')
        l_diff = list(set(data)^set(data2))
        if len(l_diff) > 0:
            logger.write("[DIFF INFO]The difference of the startup program has existed")
            for i in l_diff:
                logger.write(f"[DIFF]{i}")
                notification.notify(
                    title='The difference of startup program has existed',
                    message=i,
                    app_name='Diff notify'
                )
        else:
            logger.write("[INFO]The difference did not exsist")

CSVファイルの扱いで少々長くなってしまいましたが、
動きとしては以下の通りです。

1. CSVファイルをローテーション
2. 最新のCSVと一世代前のCSVファイルの内容を比較
3. 差異が存在する場合ログに書き込み＋ポップ通知

実行部分

if __name__ == "__main__":
    os.chdir(r"C:\python\notebooks\StartUpProgramSec")
    tooltf = os.path.exists(f"{os.getcwd()}\\autorunsc.exe")
    download_tool(tooltf)
    get_log()
    logger.write("[INFO]The process completed")

以上です。

ちょっとでも在宅ワークを安全なものにして、
普及されることを祈って…

AtCoder Beginner Contest 174 参戦記

9ペナの大爆死. 終了1分半前にギリギリE問題を通したおかげで、レーティングが小ダウンで済んだ.

ABC174A - Air Conditioner

1分で突破. 書くだけ.

X = int(input())

if X >= 30:
    print('Yes')
else:
    print('No')

ABC174B - Distance

2分半で突破. 書くだけ.

N, D = map(int, input().split())

result = 0
for _ in range(N):
    X, Y = map(int, input().split())
    if X * X + Y * Y <= D * D:
        result += 1
print(result)

ABC174D - Alter Altar

20分くらい?で突破. RE×1, WA×4. 一番左にある白い石を一番右にある赤い石とスワップするを繰り返して、白い石の右側に赤い石がない状態にすれば良い. 言うのは簡単なのに、実装でハマりました…….

N = int(input())
c = input()

d = list(c)

i = 0
j = N - 1
result = 0
while True:
    while i < N and d[i] != 'W':
        i += 1
    while j > 0 and d[j] != 'R':
        j -= 1
    if i == N or j == -1 or i >= j:
        break
    d[i] = 'R'
    d[j] = 'W'
    result += 1
print(result)

ABC174C - Repsept

28分くらい?で突破. C問題だからナイーブに書いても行けるかと思ったら、入力例3で TLE 食らって真っ青. 多分最大値は K - 1 だろうとヤマを張って出したら通った…….

K = int(input())

t = 7
for i in range(K):
    if t % K == 0:
        print(i + 1)
        break
    t = (t * 10 + 7) % K
else:
    print(-1)

追記: K で割った余りの0以外は K - 1 種類なので、最大値は K - 1 に決まっていた. 途中で同じ値がでてくるとループなので0になることはない.

ABC174E - Logs

47分で突破. TLE×1、RE×3. heapq で正答が出るコードを書いて提出してから K≤10⁹ に気づくチョンボ. 24分間 heapq のまま高速化の努力をした後に、これは駄目だとなった瞬間ににぶたんを閃いて、13分でコードを書いたらもう終了1分半前. 通った、助かった. しかし、 abs(ng - ok) > 1 でよく通ったな.

from math import ceil

N, K, *A = map(int, open(0).read().split())

def is_ok(n):
    t = 0
    for a in A:
        if a <= n:
            continue
        t += ceil(a / n) - 1
    return t <= K

ok = 1000000000
ng = 0.0000000001
while abs(ng - ok) > 1:
    m = (ok + ng) // 2
    if is_ok(m):
        ok = m
    else:
        ng = m

print(ceil(ok))

こちらの記事は、Rhea Moutafis 氏により2020年5月に公開された『 Bye-bye Python. Hello Julia! 』の和訳です。
本記事は原著者から許可を得た上で記事を公開しています。

バイバイ Python。 ハロー Julia！

Pythonの勢いに歯止めがかかると同時に新しい競争相手の登場だ

---

Juliaがまだあなたにとって未知であっても、心配しないでほしい。 Photo by Julia Caesar on Unsplash

誤解しないでほしい。 Pythonの人気は、コンピュータ科学者、データサイエンティスト、AIスペシャリストといった堅固なコミュニティによって支えられている。

しかし、これらの人々と一緒に夕食をともにしたことがあれば、彼らがPythonの弱点についてどれほどわめき散らしているのかも知っているだろう。 速度が遅いことに始まり過度のテストが必要になること、以前のテストにもかかわらずランタイムエラーが発生することまで、うんざりすることはたくさんある。

そのため、ますます多くのプログラマーが他の言語を使うようになってきている。トッププレイヤーはJulia、Go、Rustだ。 Juliaは数学的・技術的なタスクに最適である一方、Goはモジュラープログラムに最適であり、Rustはシステムプログラミングの最高の選択肢だ。

データサイエンティストとAIスペシャリストは多くの数学的問題を扱っているため、彼らにとってJuliaはその中の勝者だ。 そして、厳しく見てみても、JuliaにはPythonが勝つことのできない利点が存在する。

• Pythonが将来のプログラミング言語ではない理由

Pythonの禅 vs. Juliaの貪欲

人々が新しいプログラミング言語を作り出すのは古い言語の優れた機能を維持し、悪い部分を修正したいという理由があるときだ。

このような理由から、Guido van Rossumは1980年代後半にABCを改善するためにPythonを生み出した。 ABCはプログラミング言語としてはあまりに完璧すぎた。その柔軟性のない構造のため教えやすかったが、現実的な使用には適していなかったのだ。

対照的に、Pythonは非常に実用的だ。 作者たちの意図を反映した Zen of Python ¹の中に、それを見ることができる。

醜いより美しいほうがいい。

暗示するより明示するほうがいい。

複雑であるよりは平易であるほうがいい。

それでも、込み入っているよりは複雑であるほうがまし。

ネストは浅いほうがいい。

密集しているよりは隙間があるほうがいい。

読みやすいことは善である。

特殊であることはルールを破る理由にならない。

しかし、実用性を求めると純粋さが失われることがある。

...

それでも、PythonはABCの優れた点を引き継いでいる。例えば、読みやすさ、シンプルさ、初心者にとって優しい点などだ。 しかし、PythonはABCに比べてはるかに堅牢で、現実的な使用に適している。

ABCは、Pythonへの道を切り開き、PythonはJuliaへの道を切り開いた Photo by David Ballew on Unsplash

同様に、Juliaの作者たちは他の言語の良い部分を引き継いで、悪い部分を捨てたいと考えている。 しかし、Juliaの方がはるかに野心的だ。1つの言語を置き換えるのではなく、すべての言語を打ち負かしたいと考えているのだ。

Juliaの作者たちが言うにはこうだ：

つまり僕たちは欲張りなんだ。それじゃ物足りないんだ。

僕たちが欲しい言語はこんな感じだ。まず、ゆるいライセンスをもったオープン ソース言語。 そこに、Cの速さとRubyの動的さが欲しい。 Lispのような同図像性があって真のマクロを使えるけど、Matlabのように分かりやすくて自然な数学的記述もできる言語だ。 そのうえ、Python並みに汎用プログラミングに使えて、R並みに簡単に統計分析ができて、Perl並みに文字列処理を自然にできて、Matlab並みに線形代数計算に強くて、シェル並みにうまくプログラムをつなぎ合わせることができる。 さらに、すっごい簡単に覚えられるけど、凄腕ハッカーも満足させられるものだ。 インタラクティブに使えて、かつコンパイルもできる。

Juliaは、現在存在するすべての利点をブレンドしつつも、それを他の言語に存在する欠点と引き換えにしようとしない。 Juliaは若い言語であるものの、作者たちが設定した目標の多くはすでに達成している。

Julia開発者たちがJuliaを好んで使う理由

汎用性

Juliaは、単純な機械学習アプリケーションから巨大なスーパーコンピュータシミュレーションまで、あらゆる用途に使用可能だ。 ある程度はPythonでも実現可能だ。ところが、Pythonは何とかしてこれをお得意のものにした。

それとは対照に、Juliaはまさにこのために作られた。 １からだ。

速度

Juliaの作者たちはCと同じくらい高速な言語を作成したいと考えていたが、できたものはさらに高速だった 。 近年、Pythonの高速化は容易になってきたが、そのパフォーマンスはJuliaができることとはかけ離れている。

2017年にJuliaは Petaflop Club にまで加わった。それはピークパフォーマンスで毎秒1ペタフロップの速度を超えることができる言語の小さなクラブだ。 Juliaを除くと現在そのクラブにいるのはC、C++、Fortranだけだ。

• Pythonコードを高速化するための10のトリック

コミュニティ

Pythonには30年以上の歴史があり、巨大で協力的なコミュニティがある。 1回のGoogle検索では回答を得られないPython関連の質問はほとんどない。

対照的に、Juliaコミュニティは非常に小さい。 これは、答えを見つけるためにより調べるのに時間がかかることを意味する一方で、同じ人と何度も関わるかもしれないということだ。 そして、これはかけがいのないプログラマー関係になる可能性があるということだ。

コード変換

Juliaでコーディングするために、Juliaコマンドを１つも知る必要さえない。 Julia内でPythonおよびCコードを使うことができるのみならず、 Python内でJulia を使うことさえもできてしまう！

言うまでもなく、これにより、Pythonコードの弱点を極めて簡単に修正することができる。 言いかえれば、Juliaを学びながら生産性を維持することもできるということだ。

ライブラリは依然としてPythonの強みだ Photo by Susan Yin on Unsplash

ライブラリ

これは、Pythonの最も優れた点の1つだ。膨大な数のよくメンテナンスされたライブラリ群を抱えているということだ。 Juliaには多くのライブラリがなく、ユーザーはそれらが（まだ）驚くほどメンテナンスされていないと不満を述べている。

しかし、Juliaが限られたリソースを持つ非常に若い言語であることを考慮すると、既に存在するライブラリの数はかなり印象的だ。 Juliaのライブラリの量が増加しているという事実は別としても、Juliaはプロットを処理するためにCおよびFortranのライブラリとインターフェースを持つこともできる。

動的かつ静的な型

Pythonは100％の動的型付け言語だ。 これは、例をあげるとプログラムは実行時に変数が浮動小数点か整数かを決定しているということを意味する。

これは初心者に非常にやさしい一方で、多くのバグをもらしうる。 つまり、考えられるすべてのシナリオでPythonコードをテストする必要があるということだ。これは、多くの時間を要する非常にばかげたタスクだ。

Juliaの作者たちもJuliaを簡単に学習できるようにしたかったため、Juliaは動的型付けを完全にサポートしている。 しかし、Pythonとは対照的に、必要に応じて静的な型を導入することができる。たとえば、CやFortranでの静的な型と同じようにだ。

これにより、時間を大幅に節約できる。コードをテストしない言い訳を見つける代わりに、可能なところであればどこででも型を指定することができる。

• JuliaがPythonよりも優れている5つの点

データ：小さなときに投資する

StackOverflow でJulia（上）およびPython（下）のタグが付けられた質問の数。

これらすべてはかなり素晴らしいように聞こえるが、次のことを忘れないことが大切だ。JuliaはPythonに比べてまだ小さいということだ。

StackOverflowでの質問の数は、かなり良い指標の1つだ。現時点で、Pythonのタグ付けは、Juliaよりも約20倍多くなっている。

これは、Juliaが不人気であることを意味しているのではない。むしろ、多くのプログラマーに新しい言語が使われるようになるまでに時間がかかるのは当然のことだ。

次のことを考えてみてみよう。コード全体を別の言語で本当に書きたいだろうか？ いいや、そうではなくあなたは将来のプロジェクトで新しい言語を試してみたいと考えているのではないだろうか。 これがすべてのプログラミング言語でリリースから使われるようになるまでの間にタイムラグが生じる理由だ。

しかし、今すぐに使うことにした場合はどうだろうか。Juliaでは非常に多くの言語変換が可能なのでそれは難しくはない。そう、あなたは将来に投資しているのだ。 ますます多くの人々がJuliaを採用するときには、あなたはすでに彼らの質問に答えるのに十分な経験を積んでいることだろう。 また、より多くのPythonコードをJuliaに置き換えるほど、あなたのコードの耐久性は向上していくだろう。

Juliaに愛を示す時がきた Photo by Alexander Sinn on Unsplash

最後に伝えたいこと：Juliaを学び、あなたの強みにしよう

40年前、人工知能はニッチな現象にすぎなかった。 業界と投資家はそれを信じていなかった、そして多くの技術は洗練されたものではなく実用に耐えなかった。 しかし、当時それを学んだ者は今日の巨人となっている。需要が非常に高く、給与は NFLプレーヤーに相当する。

同様に、Juliaはまだ非常にニッチだ。 しかし、Juliaが大きくなったとき、勝者はJuliaを早期に採用した人だろう。

私はあなたが今Juliaを採用すれば、10年以内にうんとお金を稼ぐことが保証されると言っているわけではない。 しかし、チャンスの数は増えるだろう。

考えてみてほしい。たいていの巷のプログラマーは、職務経歴書にPythonを載せている。 そして今後数年で、さらに多くのPythonプログラマーが求人市場に出てくるだろう。 しかし、企業のPythonに対する需要が鈍化した場合、Pythonプログラマーの将来性は低下するだろう。 それは最初はゆっくりだが、必然的なものだ。

一方、Juliaを職務経歴書に載せることができれば、あなたは本当の優位性を持つことができる。 正直言って、他のPythonistaとの違いはどれほどあるだろうか。 大きな差ではないだろう。 しかし、3年後であっても、それほど多くのJuliaプログラマーは存在しないだろう。

Juliaスキルがあれば、仕事の要件を超えて興味があることを示すだけでなく、 学習意欲のある人間であり、プログラマーであることの意味をより広く理解していることも示していることになる。 言いかえれば、あなたはその仕事に適任だ。

あなたと他のジュリアプログラマーは未来のロックスターで、あなたはそれを知っている。 そういえば、 Juliaの作者たちは2012年にこのように言っている：

僕たちは言い訳できないほど欲張りだと認識しているけど、それでも全てが欲しいんだ。2年半ほど前に、この欲張りな言語を生み出すための仕事に取りかかった。 完璧ではないけど、バージョン1.0リリースまでもう少しだ。出来上がった言語は Juliaと名付けた。 すでに僕たちの無理な要求の90％を実現してくれてるけど、より完成度の高いものにするには他の人からの無理な要求も聞かないといけない。 だから、もし君が欲張りで理不尽で要求の厳しいプログラマーでもあるなら、Juliaを試してもらいたいんだ。

Pythonは今でも異常なほど人気がある。 しかし、あなたが今ジュリアを学ぶなら、それは将来もしかするとゴールデンチケットになるかもしれない。 そういう意味で、『バイバイ Python。 ハロー Julia!』

翻訳協力

Original Author: Rhea Moutafis
Original Article: Bye-bye Python. Hello Julia!
Thank you for letting us share your knowledge!

この記事は以下の方々のご協力により公開する事が出来ました。
改めて感謝致します。
選定担当: @_masa_u
翻訳担当: @_masa_u
監査担当: @r_pg10
公開担当: @_masa_u

ご意見・ご感想をお待ちしております

今回の記事は、いかがだったでしょうか？
・こうしたら良かった、もっとこうして欲しい、こうした方が良いのではないか
・こういったところが良かった
などなど、率直なご意見を募集しております。
いただいたお声は、今後の記事の質向上に役立たせていただきますので、お気軽にコメント欄にてご投稿ください。Twitterでもご意見を受け付けております。
みなさまのメッセージをお待ちしております。

---

1. 日本語訳はこちらから引用させていただきました： プログラマが持つべき心構え (The Zen of Python) ↩

Youtube

動画解説もしています。

問題

P-019: レシート明細データフレーム（df_receipt）に対し、1件あたりの売上金額（amount）が高い順にランクを付与し、先頭10件を抽出せよ。項目は顧客ID（customer_id）、売上金額（amount）、付与したランクを表示させること。なお、売上金額（amount）が等しい場合は同一順位を付与するものとする。

解答

コード

df_amount_rank = pd.concat([df_receipt[['customer_id', 'amount']] \
,df_receipt['amount'].rank(method='min', ascending=False)], axis=1)

df_amount_rank.columns = ['customer_id', 'amount', 'amount_ranking']

df_amount_rank.sort_values('amount_ranking', ascending=True).head(10)

出力

        customer_id     amount  amount_ranking
1202    CS011415000006  10925   1.0
62317   ZZ000000000000  6800    2.0
54095   CS028605000002  5780    3.0
4632    CS015515000034  5480    4.0
72747   ZZ000000000000  5480    4.0
10320   ZZ000000000000  5480    4.0
97294   CS021515000089  5440    7.0
28304   ZZ000000000000  5440    7.0
92246   CS009415000038  5280    9.0
68553   CS040415000200  5280    9.0 

解説

・PandasのDataFrame/Seriesにて、ランク列を新たに作成し、列を連結して、データを順位付けする方法です。
・数字情報をランキング形式で見たい時に使用します。
・'concat('<列名A>','<列名B>','<列名C>')'は、指定した列を連結する関数です。concatenate(=連結する)という動詞が元です。
・'rank(method(average/min/max/first),'ascending=True/False')'は、指定した列のランキングを表示する関数です。
・今回の場合、amount のランキングを rank で表示し、指定された列と連結するために concat を使用しています。
・'columns('<列名A>','<列名B>','<列名C>')'は、列名を再指定する関数です。
・P-017,018で使用した sort_values を用いて、ランキング順に並び直して表示しています。

※rankの参考記事はこちら

心がけたこと

• 型とか変数とか宣言とか、お作法の話は後回し
• とにかくプログラムを実行してもらってエラーに遭遇してもらう
• 一度にたくさん詰め込まない

前提

職場

• 製造業
• 色々あって僕も初心者も違う会社だけど同じチーム

ぼく（ざっくりいうと技術サイド）

• プログラミング経験は学生時代７年、業務で５年
• Python経験は学生時代２年、業務で１年
• 社内研修でプログラミング講師を２ヶ月経験済み

初心者（ざっくりいうと営業サイド）

• プログラミング経験なし（そもそも異業種）
• Pythonは自身では書かないが、業務で実行する機会はたまにある。
• 元会社からプログラミングを学ぶように言われているものの手も足も出ない。

はじめに

Fusion 360 に パスに沿った点 というコマンドがあります。これを Fusion 360 API で実現しようと紆余曲折した記録です

リファレンスマニュアルを調べてみる

Fusion 360 API で構築の点を作成する方法の指定は ConstructionPointInput Object を使います。リファレンスマニュアルのメソッドの中にある setBy で始まるものが UI での各コマンドに対応しているのでこの中から「パスに沿った点」を探します

setByCenter > 円/球/トーラスの中心点
setByEdgePlane > エッジおよび平面にある点
setByPoint > 頂点の点
setByThreePlanes > 3つの平面が通過する点
setByTwoEdges > 2つのエッジの通過点
あれ・・・？？？パスに沿った線に対応するメソッドがない？？？

「パスに沿った点」はないけど「パスに沿った平面」はあった！

パスに沿った点によく似たコマンドで パスに沿った平面 というのがあります。

なぜかこちらはちゃんと ConstructionPlaneInput.setByDistanceOnPath メソッドがあるようです。平面は有るのに点は無いのはなぜだー!といっても無いものは無いので開き直って代替案を考えることにしました。

パスに沿った平面をつくる > その平面とパスとの交点をつくる という手順を Fusion 360 API でやってみます

手順

スケッチ内の曲線状に点を作成する手順です

最初と最後のおまじないから途中まで

import adsk.core, adsk.fusion, traceback

def run(context):
    ui = None
    try:
        app = adsk.core.Application.get()
        ui  = app.userInterface
        design = app.activeProduct

        # Get the root component of the active design.
        rootComp = design.rootComponent

        # 
        # ここにコードを追加していく 
        # 

    except:
        if ui:
            ui.messageBox('Failed:\n{}'.format(traceback.format_exc()))

スケッチの中のスプラインを取得する

スケッチ1のなかのスプラインを取得して crvPath に代入しています

        # スケッチの中のスケッチ1を取得 > スケッチ1の中のスプラインカーブを取得
        sketch1 = rootComp.sketches.item(0)
        crvPath = sketch1.sketchCurves.sketchFittedSplines.item(0)

パスに沿った平面を作成する

setByDistanceOnPath を使ってパスに沿った面を作成します。始点を0、終点を1として0～1の値を指定するのですが、直接数値を入力するのではなく、ValueInput オブジェクトで指定する必要があります。

        # constructionPlanes オブジェクトを取得
        planes = rootComp.constructionPlanes

        # ConstructionPlaneInput オブジェクトを取得
        planeInput = planes.createInput()

        # パスに沿った平面を作成
        distance = adsk.core.ValueInput.createByReal(0.4)
        planeInput.setByDistanceOnPath(crvPath, distance)
        plane = planes.add(planeInput)

平面とパスの交点を作る

Fusion 360 には平面とパスとの交点を作るコマンドはありませんが、Sketch.intersectWithSketchPlane メソッド でスケッチ平面と曲線との交点を作ることができます。指定するパスは曲線一本だとしても配列にしないといけないようです。ここでは entities という配列に crvPath を追加して intersectWithSketchPlane メソッドの引数に入れています。

        # 新しいスケッチを追加
        sketch2 = rootComp.sketches.add(plane)

        # スケッチ平面とパスとの交点を作成
        entities = []
        entities.append(crvPath)
        skPoints = sketch2.intersectWithSketchPlane(entities)

このスクリプトを実行するとこうなります

まとめ

できたのはいいけど手数が多くて結構面倒でしたね。 pointInput.setByDistanceOnPath はリファレンスに載っていないだけで実は使えるのではと淡い期待を抱いてこんなコードを描いて実行してみた・・・けどエラーが出てしまってだめでした！

        # Get construction points
        points = rootComp.constructionPoints

        # Create construction point input
        pointInput = points.createInput()

        # Add construction point by distance on path
        distance = adsk.core.ValueInput.createByReal(0.4)
        pointInput.setByDistanceOnPath(crvPath, distance)
        point = points.add(pointInput)

早く API に実装されてほしいー！

importする

import scipy as sp
import sklearn.base
from sklearn.manifold import TSNE
import matplotlib.pyplot as plt
import pandas as pd

dataの読み込み

data = 'table_rdf.csv'
df = pd.read_csv(data,index_col=0).dropna(axis=1)

t-SNEの実行

X_reduced = TSNE(n_components=2, random_state=0).fit_transform(df)

色分けの設定
dataの読み込み

data_2 = 'y_rdf.csv'
df_2 = pd.read_csv(data_2,index_col=0).dropna(axis=1)

列方向（縦方向）にソート

df2_s = df_2.sort_index()
df2_s.head()

mergeする

df_merge = pd.merge(df, df_2, how='left',right_index=True,left_index=True)
df_merge.head()

プリントする

plt.scatter(X_reduced[:, 0], X_reduced[:, 1], c=df_merge['target'])
plt.colorbar()

PCAをする

from sklearn.decomposition import PCA
pca = PCA(n_components=2).fit(df)
values = pca.transform(df)
fig = plt.figure()
plt.scatter(values[:,0], values[:,1], c=df_merge['target'])
plt.xlabel("value1")
plt.ylabel("value2")
plt.colorbar()
plt.show()

はじめに

2021/2/19・20に実施される日本ディープラーニング協会(JDLA)Ｅ資格合格を目指して、ラビットチャレンジを受講した際の学習記録です。

ラビットチャレンジは「現場で潰しが効くディープラーニング講座」の通学講座録画ビデオを編集した教材を活用したコースです。
質問等のサポートはありませんが、E資格受験のための格安（2020年6月時点での最安値）の講座です。

詳細は以下のリンクからご確認ください。

科目一覧

応用数学
機械学習
深層学習（day1）
深層学習（day2）
深層学習（day3）
深層学習（day4）

Section1：実装演習

Tensorflow

線形回帰

• $d=3x+2$の関数に対して、ノイズの値を変更
  

• 関数$d$を$d=-6x+3$に変更して、ノイズの値も変更
  

→ ノイズの値が大きいほど、誤差が大きくなり、予測精度が下がることが確認された。

非線形回帰

• $d=-0.4x^3+1.6x^2-2.8x+1$の関数に対して、ノイズの値を変更
  

• 関数$d$を$d=2.1x^3-1.2x^2-4.8x+2$に変更して、ノイズの値も変更
  

→ ノイズの値が大きいほど、誤差が大きくなり、予測精度が下がることが確認された。

分類1層（mnist）

分類3層（mnist）

• 隠れ層のサイズを変更
  
  → 隠れ層のサイズが大きいほど、予測精度が上がることが確認された。

• optimizerを変更
  

分類CNN（mnist）

keras

線形回帰

単純パーセプトロン

Simple_perceptron

# モジュール読み込み
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Activation
from keras.optimizers import SGD

# 乱数を固定値で初期化
np.random.seed(0)

# シグモイドの単純パーセプトロン作成
model = Sequential()
model.add(Dense(input_dim=2, units=1))
model.add(Activation('sigmoid'))
model.summary()

model.compile(loss='binary_crossentropy', optimizer=SGD(lr=0.1))

# トレーニング用入力 X と正解データ T
X = np.array( [[0,0], [0,1], [1,0], [1,1]] )
T = np.array( [[0], [1], [1], [1]] )

# トレーニング
model.fit(X, T, epochs=30, batch_size=1)

# トレーニングの入力を流用して実際に分類
Y = model.predict_classes(X, batch_size=1)

print("TEST")
print(Y == T)

• np.random.seed(0)をnp.random.seed(1)に変更

発生する乱数が下記のような異なる値で固定されることになるので、学習結果も変わる。

np.random.seed()

import numpy as np

# 乱数を固定値で初期化
np.random.seed(0)
print(np.random.randn()) -> 1.764052345967664
# 乱数を固定値で初期化
np.random.seed(1)
print(np.random.randn()) -> 1.6243453636632417

• エポック数を100に変更

エポック数を増加させると、lossは減少することが確認された。

• AND回路, XOR回路に変更

OR回路、AND回路は線形分離可能だが、XORは線形分離不可能なので学習できないことが確認された。

• OR回路にしてバッチサイズを10に変更

バッチサイズを増加させると、lossは増加することが確認された。

• エポック数を300に変更

エポック数を増加させると、lossは減少することが確認された。

分類（iris）

classifiy_iris

import matplotlib.pyplot as plt
from sklearn import datasets
iris = datasets.load_iris()
x = iris.data
d = iris.target

from sklearn.model_selection import train_test_split
x_train, x_test, d_train, d_test = train_test_split(x, d, test_size=0.2)

from keras.models import Sequential
from keras.layers import Dense, Activation
# from keras.optimizers import SGD

#モデルの設定
model = Sequential()
model.add(Dense(12, input_dim=4))
model.add(Activation('relu'))
# model.add(Activation('sigmoid'))
model.add(Dense(3, input_dim=12))
model.add(Activation('softmax'))
model.summary()

model.compile(optimizer='sgd', loss='sparse_categorical_crossentropy', metrics=['accuracy'])

history = model.fit(x_train, d_train, batch_size=5, epochs=20, verbose=1, validation_data=(x_test, d_test))
loss = model.evaluate(x_test, d_test, verbose=0)

#Accuracy
plt.plot(history.history['accuracy'])
plt.plot(history.history['val_accuracy'])
plt.ylabel('accuracy', fontsize=14)
plt.xlabel('epoch', fontsize=14)
plt.legend(['train', 'test'], loc='lower right', fontsize=14)
plt.ylim(0, 1.0)
plt.show()

• 中間層の活性関数をsigmoidに変更

→　accuracyの減少が確認された。

• SGDをimportし、optimizerをSGD(lr=0.1)に変更

→　epochごとのaccuracyのばらつきが増加したように思われる。

分類（mnist）

classify_mnist

# 必要なライブラリのインポート
import sys, os
sys.path.append(os.pardir)  # 親ディレクトリのファイルをインポートするための設定
import keras
import matplotlib.pyplot as plt
from data.mnist import load_mnist

(x_train, d_train), (x_test, d_test) = load_mnist(normalize=True, one_hot_label=True)

# 必要なライブラリのインポート、最適化手法はAdamを使う
from keras.models import Sequential
from keras.layers import Dense, Dropout
from keras.optimizers import Adam

# モデル作成
model = Sequential()
model.add(Dense(512, activation='relu', input_shape=(784,)))
model.add(Dropout(0.2))
model.add(Dense(512, activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(10, activation='softmax'))
model.summary()

# バッチサイズ、エポック数
batch_size = 128
epochs = 20

model.compile(loss='categorical_crossentropy', 
              optimizer=Adam(lr=0.001, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.0, amsgrad=False), 
              metrics=['accuracy'])

history = model.fit(x_train, d_train, batch_size=batch_size, epochs=epochs, verbose=1, validation_data=(x_test, d_test))
loss = model.evaluate(x_test, d_test, verbose=0)
print('Test loss:', loss[0])
print('Test accuracy:', loss[1])

# Accuracy
plt.plot(history.history['accuracy'])
plt.plot(history.history['val_accuracy'])
plt.ylabel('accuracy')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='upper left')
# plt.ylim(0, 1.0)
plt.show()

• load_mnistのone_hot_labelをFalseに変更、または誤差関数をsparse_categorical_crossentropyに変更
  →　誤差関数として'categorical_crossentropy'を使う場合はデータの形式をone_hot_label形式に、'sparse_categorical_crossentropy'を使う場合はone_hot_labelではない形式にする必要がある。

• Adamの引数の値を変更
  
  →　学習率を0.0001から0.001に変更すると、予測精度はほぼ変わらないが、学習の早さは向上した。しかし、学習率を0.001からさらに0.01に変更すると、予測精度が低下した。

CNN分類（mnist）

classify_mnist_with_CNN

# 必要なライブラリのインポート
import sys, os
sys.path.append(os.pardir)  # 親ディレクトリのファイルをインポートするための設定
import keras
import matplotlib.pyplot as plt
from data.mnist import load_mnist

(x_train, d_train), (x_test, d_test) = load_mnist(normalize=True, one_hot_label=True)


# 行列として入力するための加工
batch_size = 128
num_classes = 10
epochs = 20

img_rows, img_cols = 28, 28

x_train = x_train.reshape(x_train.shape[0], img_rows, img_cols, 1)
x_test = x_test.reshape(x_test.shape[0], img_rows, img_cols, 1)
input_shape = (img_rows, img_cols, 1)


# 必要なライブラリのインポート、最適化手法はAdamを使う
from keras.models import Sequential
from keras.layers import Dense, Dropout, Flatten
from keras.layers import Conv2D, MaxPooling2D
from keras.optimizers import Adam

model = Sequential()
model.add(Conv2D(32, kernel_size=(3, 3),
                 activation='relu',
                 input_shape=input_shape))
model.add(Conv2D(64, (3, 3), activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))
model.add(Flatten())
model.add(Dense(128, activation='relu'))
model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='softmax'))
model.summary()

# バッチサイズ、エポック数
batch_size = 128
epochs = 20

model.compile(loss='categorical_crossentropy', optimizer=Adam(), metrics=['accuracy'])
history = model.fit(x_train, d_train, batch_size=batch_size, epochs=epochs, verbose=1, validation_data=(x_test, d_test))

# Accuracy
plt.plot(history.history['accuracy'])
plt.plot(history.history['val_accuracy'])
plt.title('model accuracy')
plt.ylabel('accuracy')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='upper left')
# plt.ylim(0, 1.0)
plt.show()

→　CNNにより、99%以上という非常に高い予測精度が得られることが確認された。

CIFAR-10

classify_cifar10

#CIFAR-10のデータセットのインポート
from keras.datasets import cifar10
(x_train, d_train), (x_test, d_test) = cifar10.load_data()

#CIFAR-10の正規化
from keras.utils import to_categorical

# 特徴量の正規化
x_train = x_train/255.
x_test = x_test/255.

# クラスラベルの1-hotベクトル化
d_train = to_categorical(d_train, 10)
d_test = to_categorical(d_test, 10)

# CNNの構築
import keras
from keras.models import Sequential
from keras.layers.convolutional import Conv2D, MaxPooling2D
from keras.layers.core import Dense, Dropout, Activation, Flatten
import numpy as np

model = Sequential()

model.add(Conv2D(32, (3, 3), padding='same',input_shape=x_train.shape[1:]))
model.add(Activation('relu'))
model.add(Conv2D(32, (3, 3)))
model.add(Activation('relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))

model.add(Conv2D(64, (3, 3), padding='same'))
model.add(Activation('relu'))
model.add(Conv2D(64, (3, 3)))
model.add(Activation('relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))

model.add(Flatten())
model.add(Dense(512))
model.add(Activation('relu'))
model.add(Dropout(0.5))
model.add(Dense(10))
model.add(Activation('softmax'))

# コンパイル
model.compile(loss='categorical_crossentropy',optimizer='adam',metrics=['accuracy'])

#訓練
history = model.fit(x_train, d_train, epochs=20)
tf.executing_eagerly()

#評価 & 評価結果出力
print(model.evaluate(x_test, d_test))

# Accuracy
plt.plot(history.history['accuracy'])
plt.title('model accuracy')
plt.ylabel('accuracy')
plt.xlabel('epoch')
plt.show()

# モデルの保存
# model.save('./CIFAR-10.h5')

RNN（2進数足し算の予測）

RNN（2進数足し算の予測）

# import tensorflow as tf
import tensorflow.compat.v1 as tf
tf.disable_v2_behavior()

# logging levelを変更
tf.logging.set_verbosity(tf.logging.ERROR)

import numpy as np
import matplotlib.pyplot as plt

import keras
from keras.models import Sequential
from keras.layers.core import Dense, Dropout,Activation
from keras.layers.wrappers import TimeDistributed
from keras.optimizers import SGD
from keras.layers.recurrent import SimpleRNN, LSTM, GRU


# データを用意
# 2進数の桁数
binary_dim = 8
# 最大値 + 1
largest_number = pow(2, binary_dim)

# largest_numberまで2進数を用意
binary = np.unpackbits(np.array([range(largest_number)], dtype=np.uint8).T,axis=1)[:, ::-1]


# A, B初期化 (a + b = d)
a_int = np.random.randint(largest_number/2, size=20000)
a_bin = binary[a_int] # binary encoding
b_int = np.random.randint(largest_number/2, size=20000)
b_bin = binary[b_int] # binary encoding

x_int = []
x_bin = []
for i in range(10000):
    x_int.append(np.array([a_int[i], b_int[i]]).T)
    x_bin.append(np.array([a_bin[i], b_bin[i]]).T)

x_int_test = []
x_bin_test = []
for i in range(10001, 20000):
    x_int_test.append(np.array([a_int[i], b_int[i]]).T)
    x_bin_test.append(np.array([a_bin[i], b_bin[i]]).T)

x_int = np.array(x_int)
x_bin = np.array(x_bin)
x_int_test = np.array(x_int_test)
x_bin_test = np.array(x_bin_test)


# 正解データ
d_int = a_int + b_int
d_bin = binary[d_int][0:10000]
d_bin_test = binary[d_int][10001:20000]

model = Sequential()

model.add(SimpleRNN(units=16,
               return_sequences=True,
               input_shape=[8, 2],
               go_backwards=False,
               activation='relu',
               # dropout=0.5,
               # recurrent_dropout=0.3,
               # unroll = True,
            ))
# 出力層
model.add(Dense(1, activation='sigmoid', input_shape=(-1,2)))
model.summary()
model.compile(loss='mean_squared_error', optimizer=SGD(lr=0.1), metrics=['accuracy'])
# model.compile(loss='mse', optimizer='adam', metrics=['accuracy'])

history = model.fit(x_bin, d_bin.reshape(-1, 8, 1), epochs=5, batch_size=2)

# テスト結果出力
score = model.evaluate(x_bin_test, d_bin_test.reshape(-1,8,1), verbose=0)
print('Test loss:', score[0])
print('Test accuracy:', score[1])

Epoch 1/5
10000/10000 [==============================] - 28s 3ms/step - loss: 0.0810 - accuracy: 0.9156
Epoch 2/5
10000/10000 [==============================] - 28s 3ms/step - loss: 0.0029 - accuracy: 1.0000
Epoch 3/5
10000/10000 [==============================] - 28s 3ms/step - loss: 9.7898e-04 - accuracy: 1.0000
Epoch 4/5
10000/10000 [==============================] - 28s 3ms/step - loss: 4.8765e-04 - accuracy: 1.0000
Epoch 5/5
10000/10000 [==============================] - 28s 3ms/step - loss: 3.3379e-04 - accuracy: 1.0000
Test loss: 0.0002865186012966387
Test accuracy: 1.0

• RNNの出力ノード数を128に変更
  Epoch 1/5
  10000/10000 [==============================] - 29s 3ms/step - loss: 0.0741 - accuracy: 0.9194
  Epoch 2/5
  10000/10000 [==============================] - 29s 3ms/step - loss: 0.0020 - accuracy: 1.0000
  Epoch 3/5
  10000/10000 [==============================] - 29s 3ms/step - loss: 6.9808e-04 - accuracy: 1.0000
  Epoch 4/5
  10000/10000 [==============================] - 28s 3ms/step - loss: 4.0773e-04 - accuracy: 1.0000
  Epoch 5/5
  10000/10000 [==============================] - 28s 3ms/step - loss: 2.8269e-04 - accuracy: 1.0000
  Test loss: 0.0002439875359702302
  Test accuracy: 1.0
  →　出力ノード数を16から128に変更した場合、予測精度が向上した。

• RNNの出力活性化関数をsigmoidに変更
  Epoch 1/5
  10000/10000 [==============================] - 28s 3ms/step - loss: 0.2498 - accuracy: 0.5131
  Epoch 2/5
  10000/10000 [==============================] - 28s 3ms/step - loss: 0.2487 - accuracy: 0.5302
  Epoch 3/5
  10000/10000 [==============================] - 28s 3ms/step - loss: 0.2469 - accuracy: 0.5481
  Epoch 4/5
  10000/10000 [==============================] - 28s 3ms/step - loss: 0.2416 - accuracy: 0.6096
  Epoch 5/5
  10000/10000 [==============================] - 27s 3ms/step - loss: 0.2166 - accuracy: 0.7125
  Test loss: 0.18766544096552618
  Test accuracy: 0.7449744939804077
  →　出力活性化関数をReLUからsigmoidに変更した場合、予測精度が低下した。

• RNNの出力活性化関数をtanhに変更
  Epoch 1/5
  10000/10000 [==============================] - 28s 3ms/step - loss: 0.1289 - accuracy: 0.8170
  Epoch 2/5
  10000/10000 [==============================] - 28s 3ms/step - loss: 0.0022 - accuracy: 1.0000
  Epoch 3/5
  10000/10000 [==============================] - 28s 3ms/step - loss: 7.1403e-04 - accuracy: 1.0000
  Epoch 4/5
  10000/10000 [==============================] - 27s 3ms/step - loss: 4.1603e-04 - accuracy: 1.0000
  Epoch 5/5
  10000/10000 [==============================] - 28s 3ms/step - loss: 2.8925e-04 - accuracy: 1.0000
  Test loss: 0.00024679134038564263
  Test accuracy: 1.0
  →　出力活性化関数をReLUからTanhに変更すると、予測精度が向上した。

• 最適化方法をadamに変更
  Epoch 1/5
  10000/10000 [==============================] - 31s 3ms/step - loss: 0.0694 - accuracy: 0.9385
  Epoch 2/5
  10000/10000 [==============================] - 32s 3ms/step - loss: 0.0012 - accuracy: 1.0000
  Epoch 3/5
  10000/10000 [==============================] - 31s 3ms/step - loss: 5.4037e-05 - accuracy: 1.0000
  Epoch 4/5
  10000/10000 [==============================] - 31s 3ms/step - loss: 3.3823e-06 - accuracy: 1.0000
  Epoch 5/5
  10000/10000 [==============================] - 31s 3ms/step - loss: 2.4213e-07 - accuracy: 1.0000
  Test loss: 5.572893907904208e-08
  Test accuracy: 1.0
  →　最適化方法をsgdからadamに変更すると、予測精度が向上した。

• RNNの入力Dropoutを0.5に設定
  Epoch 1/5
  10000/10000 [==============================] - 30s 3ms/step - loss: 0.2324 - accuracy: 0.5875
  Epoch 2/5
  10000/10000 [==============================] - 31s 3ms/step - loss: 0.2106 - accuracy: 0.6230
  Epoch 3/5
  10000/10000 [==============================] - 31s 3ms/step - loss: 0.2046 - accuracy: 0.6264
  Epoch 4/5
  10000/10000 [==============================] - 31s 3ms/step - loss: 0.2032 - accuracy: 0.6244
  Epoch 5/5
  10000/10000 [==============================] - 30s 3ms/step - loss: 0.1995 - accuracy: 0.6350
  Test loss: 0.15461890560702712
  Test accuracy: 0.8619986772537231
  →　入力Dropoutを0.5に設定すると、学習が進まなくなった。

• RNNの再帰Dropoutを0.3に設定
  Epoch 1/5
  10000/10000 [==============================] - 32s 3ms/step - loss: 0.1459 - accuracy: 0.8306
  Epoch 2/5
  10000/10000 [==============================] - 32s 3ms/step - loss: 0.0947 - accuracy: 0.9058
  Epoch 3/5
  10000/10000 [==============================] - 32s 3ms/step - loss: 0.0895 - accuracy: 0.9099
  Epoch 4/5
  10000/10000 [==============================] - 31s 3ms/step - loss: 0.0881 - accuracy: 0.9112
  Epoch 5/5
  10000/10000 [==============================] - 31s 3ms/step - loss: 0.0871 - accuracy: 0.9120
  Test loss: 0.09418297858566198
  Test accuracy: 0.9016276597976685
  →　再帰Dropoutを0.3に設定すると、予測精度が低下した。

• RNNのunrollをTrueに設定
  Epoch 1/5
  10000/10000 [==============================] - 20s 2ms/step - loss: 0.0938 - accuracy: 0.8871
  Epoch 2/5
  10000/10000 [==============================] - 20s 2ms/step - loss: 0.0032 - accuracy: 0.9999
  Epoch 3/5
  10000/10000 [==============================] - 20s 2ms/step - loss: 9.4733e-04 - accuracy: 1.0000
  Epoch 4/5
  10000/10000 [==============================] - 21s 2ms/step - loss: 5.3410e-04 - accuracy: 1.0000
  Epoch 5/5
  10000/10000 [==============================] - 20s 2ms/step - loss: 3.6424e-04 - accuracy: 1.0000
  Test loss: 0.0003125413816745379
  Test accuracy: 1.0
  →　unrollをTrueに設定（ループせずに計算を行う）した場合、予測精度がわずかに低下した。メモリが集約されたため、計算時間は早くなった。

Section2：強化学習

強化学習とは

長期的に報酬を最大化できるように環境の中で行動を選択できるエージェントを作ることを目標とする機械学習の一分野。
→　行動の結果として与えられる利益（報酬）を基に、行動を決定する原理を改善していく仕組み。

教師あり学習および教師なし学習では、データから予測すること、データに含まれるパターンを見つけ出すのがそれぞれの目標である。
一方、強化学習では、優れた方策を見つけることが目標。
環境について事前に完璧な知識があれば、最適な行動を予測し決定することは可能であるが、強化学習では不完全な知識を基に行動しながら、データを収集し最適な行動を見つけていく。

• Q学習
  行動する毎に行動価値関数を更新することにより学習を進める手法。

• 関数近似法
  価値関数や方策関数を関数近似する手法。

探索と利用のトレードオフ

過去のデータでベストとされる行動のみを常にとり続ければ、他のさらにベストな行動を見つけることはできない。（探索が足りない状態）
$\longleftrightarrow$　未知の行動のみを常にとり続ければ、過去の経験が活かせない、（利用が足りない状態）

強化学習のイメージ

行動価値関数

価値を表す関数として、状態価値関数（ある状態に注目する場合）と行動価値関数（状態と価値を組み合わせた価値に注目する場合）の2種類がある。

方策関数

方策関数とは、方策ベースの強化学習手法において、ある状態でどのような行動をとるのかの確率を与える関数のこと。

方策勾配法

$$ \theta^{(t+1)} = \theta^{(t)}+\epsilon \nabla J(\theta) $$

はじめに

2021/2/19・20に実施される日本ディープラーニング協会(JDLA)Ｅ資格合格を目指して、ラビットチャレンジを受講した際の学習記録です。

ラビットチャレンジは「現場で潰しが効くディープラーニング講座」の通学講座録画ビデオを編集した教材を活用したコースです。
質問等のサポートはありませんが、E資格受験のための格安（2020年6月時点での最安値）の講座です。

詳細は以下のリンクからご確認ください。

科目一覧

応用数学
機械学習
深層学習（day1）
深層学習（day2）
深層学習（day3）
深層学習（day4）

Section1：再帰型ニューラルネットワークの概念

Recurrent Neural Network（RNN）とは、時系列データに対応可能なニューラルネットワークのことである。
時系列データは時間的順序を追って一定間隔ごとに観察され、しかも相互に統計的依存関係が認められるようなデータの系列のことであり、音声データやテキストデータがある。

RNNは時系列データを扱うため、初期の状態と過去の時間 $t-1$ の状態を保持し、次の時間 $t$ の状態を求める再帰的構造が必要になる。

$$ u^t = W_{(in)}x^t+Wz^{(t-1)}+b $$
$$ z^t = f(W_{(in)}x^t+Wz^{(t-1)}+b) $$
$$ v^t = W_{(out)}z^t+c $$
$$ y^t = g(W_{(out)}z^t+c) $$

RNNにおいてのパラメータ調整方法として、誤差逆伝播の一種であるBPTTがある。

【パラメータの更新式】
$$ W_{(in)}^{t+1} = W_{(in)}^{t}-\epsilon\frac{\partial E}{\partial W_{(in)}} = W_{(in)}^{t}-\epsilon \sum_{z=0}^{T_t}\delta^{t-z}[x^{t-z}]^T $$
$$ W_{(out)}^{t+1} = W_{(out)}^{t}-\epsilon\frac{\partial E}{\partial W_{(out)}} = W_{(out)}^{t}-\epsilon \delta^{out,t}[z^{t}]^T $$
$$ W^{t+1} = W^{t}-\epsilon\frac{\partial E}{\partial W} = W_{(in)}^{t}-\epsilon \sum_{z=0}^{T_t}\delta^{t-z}[x^{t-z-1}]^T $$
$$ b^{t+1} = b^t-\epsilon\frac{\partial E}{\partial b} = b^t-\epsilon \sum_{z=0}^{T_t}\delta^{t-z} $$
$$ c^{t+1} = c^t-\epsilon\frac{\partial E}{\partial c} = c^t-\epsilon \delta^{out,t} $$

simple_RNN（バイナリ加算）

import sys, os
sys.path.append(os.pardir)  # 親ディレクトリのファイルをインポートするための設定
import numpy as np
from common import functions
import matplotlib.pyplot as plt


def d_tanh(x):
    return 1/(np.cosh(x) ** 2)

# データを用意
# 2進数の桁数
binary_dim = 8
# 最大値 + 1
largest_number = pow(2, binary_dim)
# largest_numberまで2進数を用意
binary = np.unpackbits(np.array([range(largest_number)],dtype=np.uint8).T,axis=1)

input_layer_size = 2
hidden_layer_size = 16
output_layer_size = 1

weight_init_std = 1
learning_rate = 0.1

iters_num = 10000
plot_interval = 100

# ウェイト初期化 (バイアスは簡単のため省略)
W_in = weight_init_std * np.random.randn(input_layer_size, hidden_layer_size)
W_out = weight_init_std * np.random.randn(hidden_layer_size, output_layer_size)
W = weight_init_std * np.random.randn(hidden_layer_size, hidden_layer_size)
# Xavier
# W_in = np.random.randn(input_layer_size, hidden_layer_size) / (np.sqrt(input_layer_size))
# W_out = np.random.randn(hidden_layer_size, output_layer_size) / (np.sqrt(hidden_layer_size))
# W = np.random.randn(hidden_layer_size, hidden_layer_size) / (np.sqrt(hidden_layer_size))
# He
# W_in = np.random.randn(input_layer_size, hidden_layer_size) / (np.sqrt(input_layer_size)) * np.sqrt(2)
# W_out = np.random.randn(hidden_layer_size, output_layer_size) / (np.sqrt(hidden_layer_size)) * np.sqrt(2)
# W = np.random.randn(hidden_layer_size, hidden_layer_size) / (np.sqrt(hidden_layer_size)) * np.sqrt(2)


# 勾配
W_in_grad = np.zeros_like(W_in)
W_out_grad = np.zeros_like(W_out)
W_grad = np.zeros_like(W)

u = np.zeros((hidden_layer_size, binary_dim + 1))
z = np.zeros((hidden_layer_size, binary_dim + 1))
y = np.zeros((output_layer_size, binary_dim))

delta_out = np.zeros((output_layer_size, binary_dim))
delta = np.zeros((hidden_layer_size, binary_dim + 1))

all_losses = []

for i in range(iters_num):

    # A, B初期化 (a + b = d)
    a_int = np.random.randint(largest_number/2)
    a_bin = binary[a_int] # binary encoding
    b_int = np.random.randint(largest_number/2)
    b_bin = binary[b_int] # binary encoding

    # 正解データ
    d_int = a_int + b_int
    d_bin = binary[d_int]

    # 出力バイナリ
    out_bin = np.zeros_like(d_bin)

    # 時系列全体の誤差
    all_loss = 0    

    # 時系列ループ
    for t in range(binary_dim):
        # 入力値
        X = np.array([a_bin[ - t - 1], b_bin[ - t - 1]]).reshape(1, -1)
        # 時刻tにおける正解データ
        dd = np.array([d_bin[binary_dim - t - 1]])

        u[:,t+1] = np.dot(X, W_in) + np.dot(z[:,t].reshape(1, -1), W)
        z[:,t+1] = functions.sigmoid(u[:,t+1])
#         z[:,t+1] = functions.relu(u[:,t+1])
#         z[:,t+1] = np.tanh(u[:,t+1])    
        y[:,t] = functions.sigmoid(np.dot(z[:,t+1].reshape(1, -1), W_out))


        #誤差
        loss = functions.mean_squared_error(dd, y[:,t])

        delta_out[:,t] = functions.d_mean_squared_error(dd, y[:,t]) * functions.d_sigmoid(y[:,t])        

        all_loss += loss

        out_bin[binary_dim - t - 1] = np.round(y[:,t])


    for t in range(binary_dim)[::-1]:
        X = np.array([a_bin[-t-1],b_bin[-t-1]]).reshape(1, -1)        

        delta[:,t] = (np.dot(delta[:,t+1].T, W.T) + np.dot(delta_out[:,t].T, W_out.T)) * functions.d_sigmoid(u[:,t+1])
#         delta[:,t] = (np.dot(delta[:,t+1].T, W.T) + np.dot(delta_out[:,t].T, W_out.T)) * functions.d_relu(u[:,t+1])
#         delta[:,t] = (np.dot(delta[:,t+1].T, W.T) + np.dot(delta_out[:,t].T, W_out.T)) * d_tanh(u[:,t+1])    

        # 勾配更新
        W_out_grad += np.dot(z[:,t+1].reshape(-1,1), delta_out[:,t].reshape(-1,1))
        W_grad += np.dot(z[:,t].reshape(-1,1), delta[:,t].reshape(1,-1))
        W_in_grad += np.dot(X.T, delta[:,t].reshape(1,-1))

    # 勾配適用
    W_in -= learning_rate * W_in_grad
    W_out -= learning_rate * W_out_grad
    W -= learning_rate * W_grad

    W_in_grad *= 0
    W_out_grad *= 0
    W_grad *= 0


    if(i % plot_interval == 0):
        all_losses.append(all_loss)        
        print("iters:" + str(i))
        print("Loss:" + str(all_loss))
        print("Pred:" + str(out_bin))
        print("True:" + str(d_bin))
        out_int = 0
        for index,x in enumerate(reversed(out_bin)):
            out_int += x * pow(2, index)
        print(str(a_int) + " + " + str(b_int) + " = " + str(out_int))
        print("------------")

lists = range(0, iters_num, plot_interval)
plt.plot(lists, all_losses, label="loss")
plt.xlabel("iters_num", fontsize=14)
plt.ylabel("loss", fontsize=14)

plt.show()

• weight_init_stdを変更
  
  →　weight_init_stdを1.0から0.5に減少または2.0に増加させると、どちらも収束が遅くなった。

• learning_rateを変更
  
  →　learning_rateを0.1から0.01に減少させると学習が遅くなったが、1.0に増加させると学習が早くなり、さらに3.0に増加させると学習が進まなかった。

• hidden_layer_sizeを変更
  
  →　hidden_layer_sizeを16から8に減少させると学習が遅くなったが、32に増加させると学習が早くなり、さらに128に増加させると学習が遅くなった。

• 重みの初期化方法を変更
  ・Xavierに変更
  
  ・Heに変更
  
  
  →　どちらの初期化方法も学習が遅くなった。

• 中間層の活性化関数を変更

・ReLUに変更

・tanhに変更

Section2：LSTM（Long Short-Term Memory）

単純なRNNでは、誤差逆伝播を行う際の勾配消失のため、長距離依存関係を上手く学習することができない。
LSTMはその問題点を解消する手法として用いられている。
- Input gate:$ i_t = \sigma(W^{(i)}x_t + U^{(i)}h_{t-1} + b^{(i)} ) $
- Forget gate:$ f_t = \sigma(W^{(f)}x_t + U^{(f)}h_{t-1} + b^{(f)} ) $
- output gate:$ o_t = \sigma(W^{(o)}x_t + U^{(o)}h_{t-1} + b^{(o)} ) $
- メモリーセル：$ \tilde c_t = tanh(W^{(\tilde c)}x_t+U^{(\tilde c)}h_{t-1}+b^{(\tilde c)}, \quad c_t = i_t \circ \tilde c_t+ f_t \circ c_{t-1} $
- 状態の更新：$ h_t = o_t \circ tanh(c_t) $

• CEC(Constant Error Carousel)
  勾配消失および勾配爆発の解決方法として、勾配が1であれば解決できる。
  課題：入力データについて、時間依存度に関係なく重みが一律である。
  →　ニューラルネットワークの学習特性がないということ。

• 入力ゲートと出力ゲート
  入力ゲートと出力ゲートを追加し、それぞれのゲートへの入力値の重みを重み行列W、Uで可変可能とすることにより、CECの課題を解決。

• 忘却ゲート
  CECは過去の情報がすべて保管されているが、過去の情報が不要になった場合、削除することはできず保管され続ける。
  そこで、過去の情報が不要になった場合、そのタイミングで情報を忘却する昨日として、忘却ゲートが誕生した。

• 覗き穴結合
  CECの保存されている過去の情報を任意のタイミングで他のノードに伝播させたり、忘却させたい。
  CEC自身の値は、ゲート制御に影響を与えていない。
  →　CEC自身の値に重み行列を介して伝播可能にした構造として覗き穴結合。

Section3：GRU（Gated Recurrent Unit）

従来のLSTMでは、パラメータが多数存在していたため、計算付加が大きかった。
そこで、GRUではそのパラメータを大幅に削減し、精度は同等またはそれ以上が望めるようになった構造をしている。
- Update gate:$ z_t = \sigma(W^{(z)}x_t + U^{(z)}h_{t-1}) $
- Reset gate:$ r_t = \sigma(W^{(r)}x_t + U^{(r)}h_{t-1}) $
- 状態の更新：$ \tilde h_t = tanh(Wx_t+r_t \circ Uh_{t-1}), \quad h_t = z_t \circ h_{t-1}) + (1-z_t) \circ \tilde h_t $
Reset gateの値が0の場合、前の状態を無視。
Update gateの値が1の場合、前の状態をそのままコピー。

Section4：双方向RNN（Bidirectional RNN）

過去の情報だけでなく、未来の情報を加味することで、精度を向上させるためのモデル。
文章の推敲や機械翻訳に使用される。

Section5：Seq2Seq

通常のRNNは入力と出力の長さや順序が同じでなければならない。
一方、Seq2SeqはEncoder-Decoderモデルの一種で、入力側と出力側で別々のRNNを使う。
機械対話や機械翻訳に使用される。

• Encoder RNN
  ユーザーがインプットしたテキストデータを単語等のトークンに区切って渡す構造。
  vec1をRNNに入力し、hidden stateを出力。
  このhidden stateと次の入力vec2をまたRNNに入力してきたhidden state
  最後のvecを入れたときのhidden stateをfinal stateをしてとっておく。
  このfinal stateがthought vectorと呼ばれ、入力した文の意味を表すベクトルとなる。

• Decoder RNN
  システムがアウトプットデータを単語等のトークンごとに生成する構造。

【Decoder RNNの処理手順】
1. Decoder RNN:Encoder RNNのfinal state（thought vector）から、各tokenの生成確率を出力していく。final stateをDecoder RNNのinitial stateとして設定し、Embeddingを入力。
2. Sampling：生成確率に基いてtokenをランダムに選ぶ。
3. Embedding：2で選ばれたtokenをEmbeddingしてDecoder RNNへの次の入力とする。
4. Detokenize：1-3を繰り返し、2で得られたtokenを文字列に直す。

• HRED
  過去n-1個の発話から次の発話を生成する。
  Seq2seqでは会話の文脈無視で応答がなされたが、HREDでは前の単語の流れに即して応答されるため、より人間らしい文章が生成される。
  Seq2seqとContext RNN（Encoderのまとめた各文章の系列をまとめて、これまでの会話コンテキスト全体を表すベクトルに変換する構造）を組み合わせた構造により、過去の発話の履歴を加味した応答を生成している。
  課題1：確率的な多様性が字面にしかなく、会話の「流れ」のような多様性がない。
  →　同じコンテキスト（発話リスト）を与えられても、答えの内容が会話の流れとしては毎回同じものしか出せない。
  課題2：短く情報量に乏しい答えをしがちである。
  →　短いよくある答えを学ぶ傾向がある。

• VHRED
  HREDにVAEの潜在変数の概念を追加することにより、HREDの課題を解決した構造。

• VAE
  潜在変数zに確率分布を仮定したもの。

Section6：Word2vec

RNNでは、単語のような可変長の文字列をNNに与えることはできない。
word2vecでは学習データからボキャブラリを作成し、「ボキャブラリ数×任意の単語ベクトルの次元数」の重み行列により、大規模データの分散表現の学習を現実的な計算速度とメモリ量で実現可能にした。
結果の数値ベクトルは、生のテキストをデータの視覚化、機械学習、および深層学習に適した数値表現に変換するために使用できる。

Section7：Attention Mechanism

「入力と出力のどの単語が関連しているか」の関連度を学習する仕組み。
Decoderが各単語を出力する際の情報として、Encoderにおける各単語の隠れ状態の重み付き平均を入力として用いるという方法で、翻訳元の文の文脈情報をより詳細に捉えることができるようになる。
Attentionの機構が導入されたことによって、ニューラル機械翻訳の精度は大きく向上し、従来の統計的機械翻訳モデルの性能を凌駕することとなった。

Python、R、Juliaそれぞれの強みがあるので組み合わせて使いたい場面がしばしばあるかと思います。コードを直接呼び出す機能もありますが、データ解析の場面ではそこまで密に結合させなくても別々のスクリプトで違うステップを担当すれば十分なことも多いでしょう。たとえばPythonでデータのスクレイピングをして、Juliaでマルチスレッドに解析を回し、Rで統計解析や可視化を行うといったケースは想像しやすいかと思います。

Why Feather?

そんなとき、Pythonでpickleで保存するともちろん他のプログラミング言語にデータを持っていくことができない、いっぽうでCSVでの保存は遅かったり読み込み時にパースし直すのが面倒、など微妙にやりづらさを覚えることがあります。今回はそんなワークフロー構築での悩みを解決してくれる Featherフォーマット とその使い方を簡単に紹介します。Featherはデータ保存用の軽量なフォーマットで、シンプルなAPIで使え、プログラミング言語間での行き来が自由で、かつ読み書きが高速なことを特徴とします。

こちらの比較記事 によるとFeatherは速度、メモリ消費の両面で優れたパフォーマンスを示すようです。実際の性能はどのようなデータを格納するかによって異なるでしょうが、ともかく簡単に使えるので試してみる価値は十分にあるのではないでしょうか。

注意

Featherフォーマットは行ラベルに対応していません。そのためpandasで行ラベルを与えている場合はあらかじめ df.reset_index() してやる必要があります。 Rはぜんぜん行ラベルは使わない気がしますし、 推奨できないという話もあります 。

Feather形式読み書きのコード

Python

python.py

import pandas as pd
import feather

# read
df = feather.read_dataframe("foobar.feather")

# write
feather.write_dataframe(df, "foobar.feather")

R

r.r

library(feather)

# read
df <- read_feather("foobar.feather")

# write
write_feather(df, "foobar.feather")

Julia

julia.jl

using DataFrames
using Feather

# read
df = Feather.read("foobar.feather")

# write
Feather.write("foobar.feather", df)

これだけです。いずれの言語でも型やヘッダなど気にすることなく読み書きができるためCSVより楽だと思います。

追記：最近 Feather V2 がリリースされました。まだJuliaでの対応パッケージが出揃っていないのでここでは扱いません。本記事の内容はFeather V1を対象としたものとなります。

結果

感想

A問題を難しいと感じて逃げてしまいました。基本に従えば難しくないので反省しています。
B問題も定数を間違えてしまって時間をロスしてしまいました。勿体ないです…。
C問題はあまりやったことのない類の問題でしたが解けました。粘れてよかったです。

A問題

①$t$の範囲は$D$に依存し$D$は極めて大きくその$t$の最小値$T$を求める必要がある②$tT$では成り立つという単調性がある、これらの二つの条件があることから$T$は二分探索により求めることがあります。また、ある$t$におけるそれぞれのスライムの合計の走った距離を$O(N)$で求められれば、計算量が$O(N \log{D})$となるので間に合うプログラムを書けると考えました。

ここで、$t$における全てのスライムの合計の進んだ距離を考えますが、二つのスライム(速さは$v_1$,$v_2$)が合体すると一つのスライムになり、残ったスライムの速さは$v_1+v_2$になります。したがって、合計の進む距離に注目すればそれぞれのスライムが合体しないものとしても変わらないので、全てのスライムが(合体しない時に)$t$まででどれだけ進むかを考えれば$O(N)$で求めることができます。

また、二分探索についてはこの記事を参照してください。さらに、この問題は$O(N \log{D})$ではなく$O(N)$でも解くことができます(本解)。

A.py

n,d=map(int,input().split())
x=list(map(int,input().split()))
v=list(map(int,input().split()))
s=sum(v)

def f(k):
    return s*k

#lは偽,rは真
l,r=0,1000000000000
while l+1<r:
    k=l+(r-l)//2
    if f(k)>=d:
        r=k
    else:
        l=k
print(r)

B問題

まず、素数判定をクエリで行うのでエラトステネスの篩によって判定する可能性のある素数を全てチェックします(これによりクエリの素数判定が$O(1)$になります)。ここで$p$が素数ではないと判定されれば$-1$を出力すれば良いので、以下では$p$を素数と仮定して議論をすすめます。

この元で$A^{P!}\ (mod \ P)$を単純に$((A^1)^2)^3…$と求めていくと$mod \ P$の条件下でも間に合いません。ここで思い出すべきはフェルマーの小定理です。$剰余$の累乗に関わる問題では常識なので覚えておくべきだと思います。この定理では、素数$p$に対して$p$の倍数でない$A$について$A^{p-1}=1 \ (mod \ p)$が成り立ちます。ここで、$A^{p!}=(A^{p-1})^x \ (x=1,2,…,p-2,p)$が成り立つので($\because \ p$は素数より$2$以上)、$A$が$p$の倍数の場合は$A^{p!}=0\ (mod \ P)$、そうでない場合は$A^{p!}=1\ (mod \ P)$となります。

以上をまとめれば、$p$が素数でない時は$-1$、$p$が素数で$A\%p=0$の時は$0$、$p$が素数で$A\%p \neq 0$の時は$1$を出力すれば良いです。

コード中のMAXRRの範囲を変え忘れてWAを出しました。辛かったです。このせいでだいぶ時間を食いました。

B.cc

//コンパイラ最適化
#pragma GCC optimize("Ofast")
//インクルードなど
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

//マクロ
//forループ
//引数は、(ループ内変数,動く範囲)か(ループ内変数,始めの数,終わりの数)、のどちらか
//Dがついてないものはループ変数は1ずつインクリメントされ、Dがついてるものはループ変数は1ずつデクリメントされる
//FORAは範囲for文(使いにくかったら消す)
#define REP(i,n) for(ll i=0;i<ll(n);i++)
#define REPD(i,n) for(ll i=n-1;i>=0;i--)
#define FOR(i,a,b) for(ll i=a;i<=ll(b);i++)
#define FORD(i,a,b) for(ll i=a;i>=ll(b);i--)
#define FORA(i,I) for(const auto& i:I)
//xにはvectorなどのコンテナ
#define ALL(x) x.begin(),x.end() 
#define SIZE(x) ll(x.size()) 
//定数
#define INF 1000000000000 //10^12:∞
#define MOD 1000000007 //10^9+7:合同式の法
#define MAXR 6000000 //10^5:配列の最大のrange
//略記
#define PB push_back //挿入
#define MP make_pair //pairのコンストラクタ
#define F first //pairの一つ目の要素
#define S second //pairの二つ目の要素

#define PN true //素数
#define NPN false //素数ではない
//非本質
#define MAXRR 3000 //√MAXR以上の数を設定する

//MAXRまでの整数は素数であると仮定する(ここから削る)
vector<bool> PN_chk(MAXR+1,PN);//0indexでi番目が整数iに対応(0~MAXR)
//素数を格納する配列を用意しておく
vector<ll> PNs;

void se(){
    //0と1は素数ではない
    PN_chk[0]=NPN;PN_chk[1]=NPN;

    FOR(i,2,MAXRR){
        //たどり着いた時に素数と仮定されているなら素数
        if(PN_chk[i]) FOR(j,i,ll(MAXR/i)){PN_chk[j*i]=NPN;}
    }
    FOR(i,2,MAXR){
        if(PN_chk[i]) PNs.PB(i);
    }
}

ll modpow(ll a,ll n,ll p){
    if(n==0)return 1;
    ll t=modpow(a,n/2,p);
    t=t*t%p;
    if(n&1)t=t*a%p;
    return t;
}

signed main(){
    //入力の高速化用のコード
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    se();
    ll t;cin>>t;
    REP(_,t){
        ll a,p;cin>>a>>p;
        if(!PN_chk[p]){
            cout<<-1<<endl;
            continue;
        }
        if(a%p==0){
            cout<<0<<endl;
            continue;
        }
        cout<<1<<endl;
    }
}

C問題

類題を知っていればあとは実装を頑張れば良いので、精進の力を感じました。

まず、以下のように$r_i,c_i$をおけば下図の四つの部分に分けることができます。

この元で、求める答えは①,②,③,④のそれぞれの長方形に含まれる数の積の積を求めれば良いことがわかります。このような長方形に含まれる部分については二次元累積で事前計算をすることによりクエリ処理を$O(1)$で行うことができます。

事前計算による表の作成において、二次元累積和の場合は以下のようになります。また、それぞれのマス目を$a[x][y]$とし、以下は全て0-indexedです。

$s[x][y] := [0, x) × [0, y)$ の長方形区間の和
$s[x+1][y+1] = s[x][y+1] + s[x+1][y] - s[x][y] + a[x][y]$

したがって、二次元累積積の場合も同様にすれば以下のようになります。

$s[x][y] := [0, x) × [0, y)$ の長方形区間の積
$s[x+1][y+1] = s[x][y+1] \times s[x+1][y] \div s[x][y] \times a[x][y]$

よって、この二次元累積積を四方向に定義すればよく残りの三つの長方形区間については以下のようになります。区間の表現は数学の定義と違うのですが、雰囲気を理解してもらえればと思います。(以下を記述せずに雰囲気でやろうとしていました。方針を記述することを徹底します。)

$s[x][y] := [w-1, x) × [0, y)$ の長方形区間の積
$s[x-1][y+1] = s[x][y+1] \times s[x-1][y] \div s[x][y] \times a[x][y]$

$s[x][y] := [0, x) × [h-1, y)$ の長方形区間の積
$s[x+1][y-1] = s[x][y-1] \times s[x+1][y] \div s[x][y] \times a[x][y]$

$s[x][y] := [w-1, x) × [h-1, y)$ の長方形区間の積
$s[x+1][y+1] = s[x][y+1] \times s[x+1][y] \div s[x][y] \times a[x][y]$

これを実装することができれば、あとは開区間であることに注意すればプログラムを書くことができます。また、$10^9+7$で割った余りを求めれば良いので、自分のmodintライブラリを用いました。

追記(8/2)

全体から塗り潰した部分を割るとすると0除算が面倒なので、複数の数の積を求める問題はできるだけ積のみで表現することを意識した方が良さそうです。

C.cc

//コンパイラ最適化
#pragma GCC optimize("Ofast")
//インクルードなど
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

//マクロ
//forループ
//引数は、(ループ内変数,動く範囲)か(ループ内変数,始めの数,終わりの数)、のどちらか
//Dがついてないものはループ変数は1ずつインクリメントされ、Dがついてるものはループ変数は1ずつデクリメントされる
//FORAは範囲for文(使いにくかったら消す)
#define REP(i,n) for(ll i=0;i<ll(n);i++)
#define REPD(i,n) for(ll i=n-1;i>=0;i--)
#define FOR(i,a,b) for(ll i=a;i<=ll(b);i++)
#define FORD(i,a,b) for(ll i=a;i>=ll(b);i--)
#define FORA(i,I) for(const auto& i:I)
//xにはvectorなどのコンテナ
#define ALL(x) x.begin(),x.end() 
#define SIZE(x) ll(x.size()) 
//定数
#define INF 1000000000000 //10^12:∞
#define MOD 1000000007 //10^9+7:合同式の法
#define MAXR 100000 //10^5:配列の最大のrange
//略記
#define PB push_back //挿入
#define MP make_pair //pairのコンストラクタ
#define F first //pairの一つ目の要素
#define S second //pairの二つ目の要素

template<ll mod> class modint{
public:
    ll val=0;
    //コンストラクタ
    modint(ll x=0){while(x<0)x+=mod;val=x%mod;}
    //コピーコンストラクタ
    modint(const modint &r){val=r.val;}
    //算術演算子
    modint operator -(){return modint(-val);} //単項
    modint operator +(const modint &r){return modint(*this)+=r;}
    modint operator -(const modint &r){return modint(*this)-=r;}
    modint operator *(const modint &r){return modint(*this)*=r;}
    modint operator /(const modint &r){return modint(*this)/=r;}
    //代入演算子
    modint &operator +=(const modint &r){
        val+=r.val;
        if(val>=mod)val-=mod;
        return *this;
    }
    modint &operator -=(const modint &r){
        if(val<r.val)val+=mod;
        val-=r.val;
        return *this;
    }
    modint &operator *=(const modint &r){
        val=val*r.val%mod;
        return *this;
    }
    modint &operator /=(const modint &r){
        ll a=r.val,b=mod,u=1,v=0;
        while(b){
            ll t=a/b;
            a-=t*b;swap(a,b);
            u-=t*v;swap(u,v);
        }
        val=val*u%mod;
        if(val<0)val+=mod;
        return *this;
    }
    //等価比較演算子
    bool operator ==(const modint& r){return this->val==r.val;}
    bool operator <(const modint& r){return this->val<r.val;}
    bool operator !=(const modint& r){return this->val!=r.val;}
};

using mint = modint<MOD>;

//入出力ストリーム
istream &operator >>(istream &is,mint& x){//xにconst付けない
    ll t;is >> t;
    x=t;
    return (is);
}
ostream &operator <<(ostream &os,const mint& x){
    return os<<x.val;
}

//累乗
mint modpow(const mint &a,ll n){
    if(n==0)return 1;
    mint t=modpow(a,n/2);
    t=t*t;
    if(n&1)t=t*a;
    return t;
}



signed main(){
    //入力の高速化用のコード
    //ios::sync_with_stdio(false);
    //cin.tie(nullptr);
    ll h,w;cin>>h>>w;
    vector<vector<ll>> A(h,vector<ll>(w,0));
    REP(i,h)REP(j,w)cin>>A[i][j];
    vector<vector<vector<mint>>> s(4,vector<vector<mint>>(h+1,vector<mint>(w+1,1)));
    REP(j,h){
        REP(k,w){
            ll x,y;
            x=j;y=k;
            s[0][x+1][y+1]=s[0][x+1][y]*s[0][x][y+1]/s[0][x][y]*A[x][y];
        }
    }
    REP(j,h){
        REP(k,w-1){
            ll x,y;
            x=j;y=w-1-k;
            s[1][x+1][y-1]=s[1][x+1][y]*s[1][x][y-1]/s[1][x][y]*A[x][y];
        }
    }
    REP(j,h-1){
        REP(k,w){
            ll x,y;
            x=h-1-j;y=k;
            s[2][x-1][y+1]=s[2][x-1][y]*s[2][x][y+1]/s[2][x][y]*A[x][y];
        }
    }
    REP(j,h-1){
        REP(k,w-1){
            ll x,y;
            x=h-1-j;y=w-1-k;
            s[3][x-1][y-1]=s[3][x-1][y]*s[3][x][y-1]/s[3][x][y]*A[x][y];
        }
    }
    #if 0
    REP(i,4){
        REP(j,h){
            REP(k,w){
                cout<<s[i][j][k]<<" ";
            }
            cout<<endl;
        }
        cout<<endl;
    }
    #endif
    ll q;cin>>q;
    REP(_,q){
        ll r,c;cin>>r>>c;
        mint ans=1;
        ll x,y;
        x=r-1;y=c-1;
        ans*=s[0][x][y];
        //x=r;y=w-c;
        ans*=s[1][x][y];
        //x=h-r;y=c;
        ans*=s[2][x][y];
        //x=h-r;y=w-c;
        ans*=s[3][x][y];
        cout<<ans<<endl;
    }
}

バージョン、パッケージ管理

brew install pyenv
brew install pipenv

設定

.bash_profile または .zshrc に追記

export PYENV_ROOT="$HOME/.pyenv"
export PATH="$PYENV_ROOT/bin:$PATH"
eval "$(pyenv init -)"

source ~/.bash_profile
source ~/.zshrc
のいずれかを実行

Python

pyenv install -l
pyenv install 3.8.5
pyenv global 3.8.5

ターミナルを終了させて、もう一度起動する

python --version で 3.8.5 になっていることを確認する

Pythonでcsvから値を抽出して、for文で回そうと思ったら抽出したデータがndarrayという、カンマ（,）がないバージョンの配列になっていたので苦戦。 => 解決したのでメモ

抽出先のcsv

sample.csv

sample_name,sample_colomn
aaaaaaaaaaa,bbbbbbbbbbbbb
ccccccccccc,ddddddddddddd

抽出するためのロジック

sample.py

import pandas as pd
import numpy as np

sample_data = pd.read_csv('sample.csv',index_col='sample_colomn')
sample_colomn = sample_data.index.values
print(sample_colomn)

# ['bbbbbbbbbbbbb' 'ddddddddddddd']

['bbbbbbbbbbbbb' 'ddddddddddddd']をfor文で回そうとしたらTypeError: only integer scalar arrays can be converted to a scalar indexと怒られました。

対策

tolist()を使用して、カンマ（,）がないバージョンの配列を、普通の配列（list）に変換する

sample.py

import pandas as pd
import numpy as np

sample_data = pd.read_csv('sample.csv',index_col='sample_colomn')
sample_colomn = sample_data.index.values
print(sample_colomn.tolist())

# ['bbbbbbbbbbbbb', 'ddddddddddddd']

以上です！

参考
https://note.nkmk.me/python-numpy-list/

日本統計学会公式認定 統計検定準1級対応 統計学実践ワークブック

が統計検定の準備だけでなく統計学の整理に役立つので、R, Pythonでの実装も試みた。

問18.1

データ

• 手打ちで転記した下記
  • 18.1.csv

L1,患者,寛解
8,2,0
10,2,0
12,3,0
14,3,0
16,3,0
18,1,1
20,3,2
22,2,1
24,1,0
26,1,1
28,1,1
32,1,0
34,1,1
38,3,2

R

df<-read.csv('19.1.csv')
head(df)

A data.frame: 6 × 3
L1  患者  寛解
<int>   <int>   <int>
1   8   2   0
2   10  2   0
3   12  3   0
4   14  3   0
5   16  3   0
6   18  1   1

• 1行1患者のデータに展開

df2<-data.frame()
for (i in (1:nrow(df))) {
  yo<-df[i, 3]
  for (j in (1:df[i, 2])) {
    rem = rem - 1
    if (rem >= 0) {
      remission = 1
    } else {
      remission = 0
    }
    tmp.df<-data.frame(L1 = df[i, 1], REM = remission)
    df2<-rbind(df2, tmp.df)
  }
}
df2$REM = as.factor(df2$REM)
head(df2)

A data.frame: 6 × 2
L1  REM
<int>   <fct>
1   8   0
2   8   0
3   10  0
4   10  0
5   12  0
6   12  0

• ロジスティック回帰

model<-glm(REM ~ L1, df2, family = binomial("logit"))
summary(model)

Call:
glm(formula = REM ~ L1, family = binomial("logit"), data = df2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9448  -0.6465  -0.4947   0.6571   1.6971  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -3.77714    1.37862  -2.740  0.00615 **
L1           0.14486    0.05934   2.441  0.01464 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 34.372  on 26  degrees of freedom
Residual deviance: 26.073  on 25  degrees of freedom
AIC: 30.073

Number of Fisher Scoring iterations: 4

Python

import pandas as pd
import statsmodels.api as sm

df = pd.read_csv('18.1.csv')
df.head()

    L1  患者  寛解
0   8   2   0
1   10  2   0
2   12  3   0
3   14  3   0
4   16  3   0

df2 = pd.DataFrame(columns = ['L1', 'REM'])
for i in range(len(df)):
  rem = df.iloc[i, 2]
  for j in range(df.iloc[i, 1]):
    rem -= 1
    if rem >= 0:
      remission = 1
    else:
      remission = 0
    df2 = df2.append({'L1': df.iloc[i, 0], 'REM': remission}, ignore_index=True)
df2.head()

L1  REM
0   8   0
1   8   0
2   10  0
3   10  0
4   12  0

df_x = pd.DataFrame(df2['L1'])
df_x = sm.add_constant(df_x)
logit_mod = sm.Logit(np.asarray(df2['REM'].astype('float')), np.asarray(df_x.astype('float')))
logit_res = logit_mod.fit(disp=0)
print(logit_res.summary())

                           Logit Regression Results                           
==============================================================================
Dep. Variable:                      y   No. Observations:                   27
Model:                          Logit   Df Residuals:                       25
Method:                           MLE   Df Model:                            1
Date:                Sun, 02 Aug 2020   Pseudo R-squ.:                  0.2414
Time:                        06:58:18   Log-Likelihood:                -13.036
converged:                       True   LL-Null:                       -17.186
Covariance Type:            nonrobust   LLR p-value:                  0.003967
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -3.7771      1.379     -2.740      0.006      -6.479      -1.075
x1             0.1449      0.059      2.441      0.015       0.029       0.261
==============================================================================

問18.2, 18.3

データ

• 手打ちで転記
  • 18.2.csv

smoking,obesity,snoring,cases,hypertension
0,0,0,60,5
1,0,0,17,2
0,1,0,8,1
1,1,0,2,0
0,0,1,187,35
1,0,1,85,13
0,1,1,51,15
1,1,1,23,8

R

df<-read.csv('18.2.csv')
head(df)

A data.frame: 6 × 5
smoking obesity snoring cases   hypertension
<int>   <int>   <int>   <int>   <int>
1   0   0   0   60  5
2   1   0   0   17  2
3   0   1   0   8   1
4   1   1   0   2   0
5   0   0   1   187 35
6   1   0   1   85  13

• 1行1患者に展開

df2<-data.frame()
for (i in (1:nrow(df))) {
  ht<-df[i, 5]
  for (j in (1:df[i, 4])) {
    ht = ht - 1
    if (ht >= 0) {
      hypertension = 1
    } else {
      hypertension = 0
    }
    tmp.df<-data.frame(df[i, 1:3], HT = hypertension)
    df2<-rbind(df2, tmp.df)
  }
}
df2$HT = as.factor(df2$HT)
head(df2)

A data.frame: 6 × 4
smoking obesity snoring HT
<int>   <int>   <int>   <fct>
1   0   0   0   1
2   0   0   0   1
3   0   0   0   1
4   0   0   0   1
5   0   0   0   1
6   0   0   0   0

• 18.2 ロジスティック回帰

model.logit<-glm(HT ~ smoking + obesity + snoring, df2, family = binomial("logit"))
summary(model.logit)

Call:
glm(formula = HT ~ smoking + obesity + snoring, family = binomial("logit"), 
    data = df2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8578  -0.6330  -0.6138  -0.4212   2.2488  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.37766    0.38010  -6.255 3.97e-10 ***
smoking     -0.06777    0.27812  -0.244   0.8075    
obesity      0.69531    0.28508   2.439   0.0147 *  
snoring      0.87194    0.39749   2.194   0.0283 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 411.42  on 432  degrees of freedom
Residual deviance: 398.92  on 429  degrees of freedom
AIC: 406.92

Number of Fisher Scoring iterations: 4

• InterceptのPrが少し違うのが?だがほかはテキスト通り

• 18.3 プロビットモデル

model.probit<-glm(HT ~ smoking + obesity + snoring, df2, family = binomial("probit"))
summary(model.probit)

Call:
glm(formula = HT ~ smoking + obesity + snoring, family = binomial("probit"), 
    data = df2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8542  -0.6337  -0.6142  -0.4211   2.2531  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.37312    0.19308  -7.112 1.15e-12 ***
smoking     -0.03865    0.15682  -0.246   0.8053    
obesity      0.39996    0.16748   2.388   0.0169 *  
snoring      0.46507    0.20464   2.273   0.0230 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 411.42  on 432  degrees of freedom
Residual deviance: 399.00  on 429  degrees of freedom
AIC: 407

Number of Fisher Scoring iterations: 4

Python

• 18.2

import pandas as pd
import statsmodels.api as sm

df = pd.read_csv('18.2.csv')
df.head()

smoking obesity snoring cases   hypertension
0   0   0   0   60  5
1   1   0   0   17  2
2   0   1   0   8   1
3   1   1   0   2   0
4   0   0   1   187 35

df2 = pd.DataFrame(columns = ['smoking', 'obesity', 'snoring', 'hypertension'])
for i in range(len(df)):
  ht = df.iloc[i, 4]
  for j in range(df.iloc[i, 3]):
    ht -= 1
    if ht >= 0:
      hypertension = 1
    else:
      hypertension = 0
    df2 = df2.append({'smoking': df.iloc[i, 0], 'obesity': df.iloc[i, 1], 'snoring': df.iloc[i, 2], 'hypertension': hypertension}, ignore_index=True)
df2.head()

smoking obesity snoring hypertension
0   0   0   0   1
1   0   0   0   1
2   0   0   0   1
3   0   0   0   1
4   0   0   0   1

• 18.2 ロジスティック回帰

df_x = pd.DataFrame(df2.iloc[:, 0:3])
df_x = sm.add_constant(df_x)
logit_mod = sm.Logit(np.asarray(df2['hypertension'].astype('float')), np.asarray(df_x.astype('float')))
logit_res = logit_mod.fit(disp=0)
print(logit_res.summary())

                           Logit Regression Results                           
==============================================================================
Dep. Variable:                      y   No. Observations:                  433
Model:                          Logit   Df Residuals:                      429
Method:                           MLE   Df Model:                            3
Date:                Sun, 02 Aug 2020   Pseudo R-squ.:                 0.03040
Time:                        07:05:53   Log-Likelihood:                -199.46
converged:                       True   LL-Null:                       -205.71
Covariance Type:            nonrobust   LLR p-value:                  0.005832
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -2.3777      0.380     -6.254      0.000      -3.123      -1.633
x1            -0.0678      0.278     -0.244      0.807      -0.613       0.477
x2             0.6953      0.285      2.439      0.015       0.137       1.254
x3             0.8719      0.398      2.193      0.028       0.093       1.651
==============================================================================

• 18.3 プロビットモデル

df_x = pd.DataFrame(df2.iloc[:, 0:3])
df_x = sm.add_constant(df_x)
probit_mod = sm.Probit(np.asarray(df2['hypertension'].astype('float')), np.asarray(df_x.astype('float')))
probit_res = probit_mod.fit(disp=0)
print(probit_res.summary())

                          Probit Regression Results                           
==============================================================================
Dep. Variable:                      y   No. Observations:                  433
Model:                         Probit   Df Residuals:                      429
Method:                           MLE   Df Model:                            3
Date:                Sun, 02 Aug 2020   Pseudo R-squ.:                 0.03019
Time:                        07:07:14   Log-Likelihood:                -199.50
converged:                       True   LL-Null:                       -205.71
Covariance Type:            nonrobust   LLR p-value:                  0.006077
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -1.3731      0.193     -7.120      0.000      -1.751      -0.995
x1            -0.0387      0.157     -0.246      0.805      -0.346       0.269
x2             0.4000      0.168      2.384      0.017       0.071       0.729
x3             0.4651      0.204      2.274      0.023       0.064       0.866
==============================================================================