流れに従ひて己を失はず.

日々の研究について書きます

社会人博士課程学生周辺の環境について情報共有

社会人博士課程学生である私の周りの環境についての話.
自分の振り返りを兼ねて書きます.
この記事は社会人学生 Advent Calendar 2019の12日目の記事です.
adventar.org


・概要
大学院:情報理工学系研究科システム情報学専攻.2017年9月入学,2019年12月現在でD3.光学デバイスを用いた計測及び情報提示システムの研究に従事.
会社:2016年4月入社,2019年12月現在で4年目.職務内容は情報処理周辺の研究開発業務.

・時系列
2016年3月 修士修了
2016年4月 入社
2017年1月 博士号を意識&先輩各位から社会人博士のアドバイス
2017年3月 先生に相談
2017年6月 出願書類提出
2017年8月 入試
2017年9月 合格発表&入学

・会社のこと
職場は研究開発なので,博士号取得が推奨されています.自分の周りでも社会人博士に通っている人も何人かおり,直属の上司も通学中という状況,同期でも何人か通学してます.
社会人博士をサポートする支援制度も職場に存在しています.支援内容は,入学金・授業料・通学交通費の一部が補助されるという金銭的な面と,業務に支障が出ない範囲での平日日中の通学が認められる(大体週1日か半日というケースが多い)という時間的な面の2つです.支援対象となるにはいくつか基準があり,入社後3年半以上経過していること,入社後一定の研究成果があること,などがあります.自分の場合は入社(2016年4月)から1年半後の2017年秋入学だったため,3年半経過後の2019年10月から支援対象になり,現在は金銭面と時間面のサポートを得ました.
職場で周りの人に相談しても,ほとんどがポジティブな反応でしたが,一方で年次が若く前例が無いことを理由に反対・心配される事もありました.まあそれは2年目の終わりにある研修員発表(この苦労を経験してやっと一人前になったな!うん!というイベント)があるからだったのだろうな,と今から考えると思います.実際のところ,私見ですが大変なこともいろいろあるけれど,年次が進むと職務に加えて役職・委員の仕事が増える&ライフステージの変化があるので,行くなら若い時期から進学したほうが良いと思います.
※ここらへんの,会社の制度や待遇,そして平日何日くらい大学に行けているかの情報がもっとネットに集まってくるといいなあと考えてます.

・大学のこと
私は出身研究室に再び入れていただくこととなりました.先生にご相談させていただいてから入学までの半年間は,入試の準備と入学後の研究の下準備を進めていました.6月提出の出願書類にすでに研究計画書を入れなければならず,ビジョンが固まっていない状態でしたので作成に苦労しました.あとは,社会人博士なので「職場としてこの人を通学させていいと考えてる」ことを証明する書類を会社に作ってもらって出願時に大学に提出しました.
入試は専門科目と英語,研究内容のプレゼンがありますが,自分の場合修士の出身研究科ということで専門科目と英語は免除でした.研究内容のプレゼンの内容は修士の研究紹介と博士課程での研究計画です.出願書類でも大変でしたが研究計画プレゼンの準備が結構大変でした.あんまりできるビジョンが見えていすぎても駄目だし,全くできる気がしなくても駄目なのでそのへんの塩梅が難しかったです.先輩方からは「大丈夫だよ〜」と言われてましたが,プレゼン前は結構緊張しました.落ちる人は落ちます.
8月23日にプレゼンがあって,9月4日に合格発表でした.プレゼン後の質疑の雰囲気は良かったのでそこまで心配はしていなかったんですが,掲示をみてやっと一安心という感じでした.
社会人博士が受けられる金銭面のサポートとして,博士課程研究遂行協力制度(参考:https://giju.t.u-tokyo.ac.jp/eejim/wp-content/uploads/2018/05/web_guide_suiko_20180531.pdf)というものがあります.半年間月5万円支給なので学費に充てることができ,とても助かりました.
授業は,博士課程研究以外に通算4コマ8単位を取らなければなりません.修士のときに余分に単位をとっていれば博士課程での修了要件の単位数に加えることができたんですが,私はとっていなかったので講義をとる必要があります.大学を出た後にとっておけばよかったなと思っていた,他分野の講義をとるようにしました.
本来は博士課程学生がやらなかればならない研究室の仕事を免除していただいており,研究室の関係者各位には頭が上がりません.

・普段の生活
土日休み,平日4日会社,1日大学です.
夜は自宅で毎日だいたい21時〜24時で作業します.大学に行く日には装置を組んだりこれを使った実験をしてます.シミュレーションや調べ物,コーディングは自宅からやってます.基本的にゲームやったり工作したりなどの時間は無いですが,週2回は身体を動かす時間を確保できている状況です.

ベクトル正規化(3次元)の速度比較

np.linalg.norm()を使うと遅い.
二乗和はぐりぐり計算して,平方根はmath.sqrt()を使うのが良い.
という話.


np.linalg.norm()を使うのと,np.dot()を使うのと,math.sqrt()を使うのと,愚直に計算するのではどれが一番速いか比べた.

ライブラリインポートとデータ作成

import numpy as np
import math

vec = np.random.rand(1000,3)


np.linalg.norm()を使った場合:9.32 ms

def vec_3d_length_np(v):
    return v/np.linalg.norm(v)

%%time
for _ in range(1000):
    vn = vec_3d_length_np(vec[_,:])

CPU times: user 8.94 ms, sys: 3.39 ms, total: 12.3 ms
Wall time: 9.32 ms


np.dot()を使った場合:7.08 ms

def vec_3d_length_npdot(v):
    return v/(np.dot(v,v))**(0.5)

%%time
for _ in range(1000):
    vn = vec_3d_length_npdot(vec[_,:])

CPU times: user 6.08 ms, sys: 4.31 ms, total: 10.4 ms
Wall time: 7.08 ms


math.sqrt()を使った場合:2.35 ms

def vec_3d_length_math(v):
    return v/math.sqrt(v[0]*v[0] + v[1]*v[1]+v[2]*v[2])

%%time
for _ in range(1000):
    vn = vec_3d_length_math(vec[_,:])

CPU times: user 2.33 ms, sys: 6 µs, total: 2.33 ms
Wall time: 2.35 ms


愚直に計算した場合:2.73 ms

def vec_3d_length_(v):
    return v/(v[0]**2 + v[1]**2+v[2]**2)**(0.5)

%%time
for _ in range(1000):
    vn = vec_3d_length_(vec[_,:])

CPU times: user 2.72 ms, sys: 3 µs, total: 2.72 ms
Wall time: 2.73 ms

スネルの法則を使って屈折方向の計算

スネルの法則で3次元空間上の屈折方向ベクトルを求めるのをPythonで実装.


ベクトルで解く

各要素の関係は以下.ただし各ベクトルの大きさを1とする.
入射ベクトル:v_{inter}
法線ベクトル:v_{norm}
出射ベクトル(これを求める):v_{exit}
入射角:\theta_{inter}
出射角(これもわかってない):\theta_{exit}
入射側屈折率:n_{inter}
入射側屈折率:n_{exit}

f:id:Yasutchi:20190914234142j:plain


スネルの法則
n_{exit} sin\theta_{exit} = n_{inter} sin\theta_{inter}. \tag{0}
つまり,\theta_{exit}はぐりぐりやればわかる.


ここで,入射ベクトル,法線ベクトル,出射ベクトルは同一平面上にあるから, \alpha, \beta(ともにスカラ)を使って以下のように表現できる.

v_{exit} = \alpha\ v_{inter} + \beta\ v_{norm}. \tag{1}

まず \alphaから求める.
式(1)の両辺について,法線ベクトルとの外積をとって,更にその大きさをとると,各ベクトルの大きさが1であること,最後は式(0)を用いて\alphaが決まる.


                \begin{eqnarray}

                |v_{norm} \times v_{exit}| &=& |v_{norm} \times (\alpha\ v_{inter} + \beta\ v_{norm})|,\\|v_{norm} \times v_{exit}| &=& |v_{norm} \times \alpha\ v_{inter} + v_{norm} \times \beta\ v_{norm}|,\\|v_{norm}||v_{exit}|sin(\pi - \theta_{exit}) &=& \alpha|v_{norm}||v_{inter}|sin(\pi - \theta_{inter}) + 0,\\sin\theta_{exit} &=& \alpha sin\theta_{inter},\\ \alpha &=& \dfrac{sin\theta_{exit}}{sin\theta_{inter}} = \dfrac{n_{inter}}{n_{exit}}.

                \end{eqnarray}



次に \betaを求める.
式(1)の両辺について,法線ベクトルとの内積をとると,先程と同様に各ベクトルの大きさが1であることを用いて\betaが決まる.\theta_{exit}が残っているが,式(0)から求めることができる.

                \begin{eqnarray}

                v_{norm} \cdot v_{exit} &=& v_{norm} \cdot (\alpha\ v_{inter} + \beta\ v_{norm}),\\v_{norm} \cdot v_{exit} &=& \alpha\ v_{norm} \cdot v_{inter} + \beta\ v_{norm} \cdot v_{norm},\\ |v_{norm}||v_{exit}|cos(\pi - \theta_{exit}) &=& \alpha|v_{norm}||v_{inter}|cos(\pi - \theta_{inter}) + \beta,\\-cos\theta_{exit} &=& -\alpha cos\theta_{inter} + \beta,\\ \beta &=& \alpha cos\theta_{inter} -cos\theta_{exit}.

                \end{eqnarray}



Python実装

import numpy as np 

def snell(v_inter,v_norm,n_inter,n_exit):
    #3次元スネルの公式の適用結果を返す関数
    #入力:入射光線の方向ベクトル,入射点での法線ベクトル,入射元の屈折率,入射先の屈折率
    
    #入力されたベクトルの正規化
    v_inter = v_inter / np.linalg.norm(v_inter)
    v_norm = v_norm / np.linalg.norm(v_norm)
    
    #係数a
    a =  n_inter/n_exit
    
    cos_t_inter = -1.0*np.dot(v_inter,v_norm)
    #量子化誤差対策
    if cos_t_inter<-1.:
        cos_t_inter = -1.
    elif cos_t_inter>1.:
        cos_t_inter = 1.
        
    sin_t_inter = (1.0 - cos_t_inter**2)**(1/2)
    sin_t_exit = sin_t_inter*a
    
    if sin_t_exit>1.0:
        #全反射する場合
        return np.zeros(3), False
    
    cos_t_exit = (1 - sin_t_exit**2)**(1/2)
    
    #係数b
    b = -1.0*cos_t_exit + a*cos_t_inter
    
    #出射光線の方向ベクトル
    v_exit = a*v_inter + b*v_norm
    #正規化
    v_exit = v_exit / np.linalg.norm(v_exit)
    
    return v_exit, True

vi = np.array([0.08584,0.17301,0.9811726])
vn = np.array([0.050, 0.060, -0.9969453])
ni = 1.0
ne = 1.5

snell(vi,vn,ni,ne)#(array([0.0401461 , 0.0948433 , 0.99468238]), True)

#方向ベクトルなどの実値参考:
#https://www.cybernet.co.jp/optical/course/optics/opt02/opt06.html


参考:http://yees.g.dgdg.jp/RayTrace/RayTraceVersion001b.pdf

SymPyでの幾何計算速度

pythonのライブラリであるSymPyを使って幾何計算をやるととても遅いという話.

以下は2直線の交点導出を100セット解いた例.
SymPyを使った場合とベクトルで解いた場合の比較.10^3オーダーで変わってくる.

SymPyを使った場合

import numpy as np
import sympy.geometry as sg
import time

#点[listA]を通り方向ベクトル[listB]を持つ直線と
#点[listC]を通り方向ベクトル[listD]を持つ直線との交点を求める
listA = np.random.rand(100,2)
listB = np.random.rand(100,2)
listC = np.random.rand(100,2)
listD = np.random.rand(100,2)

start = time.time()

for i in range(0,100):
    pointA = listA[i,:]
    direcB = listB[i,:]
    pointC = listC[i,:]
    direcD = listD[i,:]
    line1 = sg.Line( sg.Point(pointA), 
                                        sg.Point(pointA+direcB))
    line2 = sg.Line( sg.Point(pointC), 
                                        sg.Point(pointC+direcD))
    result = sg.intersection(line1, line2)

elapsed_time = time.time() - start
print ("sympy calc_time:{0}".format(round(elapsed_time,3)) + "[sec]")

sympy calc_time:15.911[sec]


ベクトルで解いた場合

import numpy as np
import sympy.geometry as sg
import time

#点[listA]を通り方向ベクトル[listB]を持つ直線と
#点[listC]を通り方向ベクトル[listD]を持つ直線との交点を求める
listA = np.random.rand(100,2)
listB = np.random.rand(100,2)
listC = np.random.rand(100,2)
listD = np.random.rand(100,2)

#90度の回転行列
rot90 = np.array([[np.cos(np.pi/2), -np.sin(np.pi/2)],
                              [np.sin(np.pi/2), np.cos(np.pi/2)]])

def calcintersection(point1,direc1,point2,direc2):
    # 参考
    # https://qiita.com/tmakimoto/items/2da05225633272ef935c
    
    #2直線が平行なら[]を返す
    if direc1[0]*direc2[1]==direc1[1]*direc2[0]:
        return []
    
    #方向ベクトルの正規化
    direc1_n = direc1 / np.linalg.norm(direc1) 
    direc2_n = direc2 / np.linalg.norm(direc2) 
    
    #direc2_nに垂直な単位ベクトル
    direc3_n = np.dot(direc2_n, rot90)
    
    #直線1における交点のパラメータloを求める
    lo = np.dot(direc3_n,(point2-point1)) / np.dot(direc3_n,direc1_n)
    
    return point1+direc1_n*lo


start = time.time()

for i in range(0,100):
    pointA = listA[i,:]
    direcB = listB[i,:]
    pointC = listC[i,:]
    direcD = listD[i,:]
    
    result = calcintersection(pointA,direcB,pointC,direcD)

elapsed_time = time.time() - start
print ("vector calc_time:{0}".format(round(elapsed_time,3)) + "[sec]")

vector calc_time:0.004[sec]

pythonのmultiprocessingで複数引数複数戻値を並列処理

毎回忘れるのでメモ.

pythonで複数引数をとって,複数戻値しかも配列を返す処理を並列で行う.

参考: https://qiita.com/wikipediia/items/2919362de582a7d8de9e

import time
import numpy as np
import os
from multiprocessing import Pool

#総スレッド数の取得
cpu_count = os.cpu_count()

#最終的に作りたい出力結果の行数
df_size = 1234

#作りたいデータその1
#df_size*2の2次元データ
#n行目のデータは[n,2*n+1]にしたい
data_frame = np.ones([df_size,2],dtype=int)

#作りたいデータその2
#df_size*2の2次元データ
#n行目のデータは[n,n^2+1]にしたい
data_frame_s = np.ones([df_size,2],dtype=int)

#それぞれのスレッドに割振る行数の開始番号と終了番号の組
sta_gol = np.ones([cpu_count,2],dtype=int)
for i in range(0,cpu_count):
    sta = (df_size//cpu_count)*i
    gol = (df_size//cpu_count)*(i+1)-1
    if i == cpu_count-1:
        gol = df_size-1
    sta_gol[i,:] = [sta,gol]
print(sta_gol)

# ここを並列処理する
def func(i,arg1, arg2):
    sta = int(arg1)
    gol = int(arg2)
    print("Thread for sta:",sta," gol:", gol)
    
    #作りたいデータの割り当てられた各部分を作る処理
    data_frame_sg = np.ones([gol+1-sta,2],dtype=int)
    data_frame_sg_s = np.ones([gol+1-sta,2],dtype=int)
    for r in range(sta,gol+1):
        data_frame_sg[r-sta,:] = [r,2*r+1]
        data_frame_sg_s[r-sta,:] = [r,r*r+1]
    
    #ダミーの処理
    time.sleep(5/(i+1))
    
    return i,data_frame_sg,data_frame_sg_s
    
def wrapper(args):
    # argsは(i,sta, gol)となっている
    return func(*args)

def multi_process(sampleList):
    # プロセス数:cpu_count
    p = Pool(cpu_count)
    #20220718追記
    # M1Macだとp = get_context("fork").Pool(cpu_count)
    # 参考:https://stackoverflow.com/questions/67999589/
    
    #ラッパにわたす
    output = p.map(wrapper, sampleList)
    
    # プロセスの終了
    p.close()
    return output

def mproc():
    # (i,sta, gol)が引数になる
    sampleList = [(i,sta_gol[i,0],sta_gol[i,1]) for i in range(cpu_count)]

    return multi_process(sampleList)

start = time.time()

#実行するところ
mproc_output = mproc()

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(round(elapsed_time,2)) + "[sec]")

#並列処理の結果をもとのスキーマに直すところ
for i in range(0,len(mproc_output)):
    ind = mproc_output[i][0]
    data_frame[sta_gol[ind,0]:sta_gol[ind,1]+1,:] = mproc_output[i][1]
    data_frame_s[sta_gol[ind,0]:sta_gol[ind,1]+1,:] = mproc_output[i][2]
print(data_frame)
print(data_frame_s)


深夜の馬鹿力で出題された問題を解く

TBSラジオの超人気番組,月曜JUNK伊集院光深夜の馬鹿力の「少しおもしろい」のコーナーで気になる投稿があった.

【問題】
都道府県の最初の漢字だけで表すことのできる(多分)唯一の芸能人は長島三奈
(月曜JUNK伊集院光深夜の馬鹿力 20181224深夜放送分より)


www.tbsradio.jp



伊集院さんも他にないのかなー??と唸っていた.
芸能人のお名前を片っ端から調べてみる.
芸能人のリストとしてgooのタレント一覧のページ(https://news.goo.ne.jp/entertainment/talent)を使う.

news.goo.ne.jp



残念ながらこちらには,放送内で先に正解例として挙がっている長島三奈さんは含まれていない.他に良さそうな情報源が発見できなかったのでこちらを使用.

目視で調べるのもあれなのでスクレイピングして芸能人の名前リストを取ってくる.
スクレイピングについての参考と注意点は以下.

https://qiita.com/Azunyan1111/items/9b3d16428d2bcc7c9406
https://docs.pyq.jp/column/crawler.html



gooのタレント一覧のページのソースをみると,<p class="list-title-news margin-bottom5">のあとにタレント名が記載されているようなので,
この要素だけ見るようにする.

【https://news.goo.ne.jp/entertainment/talent/j-a/ より】
(中略)
<ul class="gn-thumbs-list thumbs-list-large">
<li><a href="/entertainment/talent/W17-0923.html">
<div class="gn-news-item news-list-item">
<div class="news-item-thumbs"><img src="https://img.news.goo.ne.jp/talent/MW-W17-0923.jpg" alt="" class="thumbs-xx-small"></div>
</div>
<div class="image-news-block">
<p class="list-title-news margin-bottom5">あぁ〜しらき</p>
<p class="list-news-source">アァ〜シラキ</p>
</div>
</a></li>
<li><a href="/entertainment/talent/M10-0879.html">
<div class="gn-news-item news-list-item">
<div class="news-item-thumbs"><img src="https://img.news.goo.ne.jp/talent/MM-M10-0879.jpg" alt="" class="thumbs-xx-small"></div>
</div>
<div class="image-news-block">
<p class="list-title-news margin-bottom5">Asu</p>
<p class="list-news-source">アアス</p>
</div>
(以下略)


https://news.goo.ne.jp/entertainment/talent/j-以下がa-1,a-2,a-3...となっていくので総当たりで調べていく.
ページのソースが回収できたら,先程の<p class="list-title-news margin-bottom5">のあとに続く文字列が,都道府県の最初の漢字だけかを判定.

【name.py】
import urllib.request
from bs4 import BeautifulSoup
import time

#都道府県名最初の一文字と全角スペースのリスト
list_todohuken=[" ","北","青","岩","秋","宮","山","福","茨","栃","群","埼","東","長"
               ,"千","神","静","新","富","石","岐","愛","滋","三","京","奈","和","大","兵","鳥"
               ,"島","岡","広","香","徳","高","佐","熊","鹿","沖"]

#与えられた人名が,すべて都道府県名最初の一文字かを検証する関数
def check_name_todohuken(name):
    for c in name:
        if c in list_todohuken:
            continue
        else:
            return

    #与えられた人名が,すべて都道府県名最初の一文字なら出力
    print(name)

def check_list_url(url):
    # html取得
    html = urllib.request.urlopen(url)

    # htmlをBeautifulSoupで扱う
    soup = BeautifulSoup(html, "html.parser")

    # p要素全てを摘出する
    p = soup.find_all("p")

    # class=""list-title-news margin-bottom5"となっているところを探す
    for tag in p:
        try:
            string_ = tag.get("class").pop(0)

            if string_ in "list-title-news margin-bottom5":
                check_name_todohuken(tag.string)
        except:
            pass

shiin = ["","k","s","t","n","h","m","r"]
boin = ["a","i","u","e","o"]
other = ["ya","yu","yo","wa"]
num = ["1","2","3","4","5"]

#子音と母音のセットで各ページをチェック
for s in shiin:
    for b in boin:
        for n in num:
            check_list_url("https://news.goo.ne.jp/entertainment/talent/j-"+s+b+"-"+n+"/")
            #サーバアクセスの間隔を1秒以上空ける
            time.sleep(1)
            
#「や」「ゆ」「よ」「わ」は個別にチェック
for o in other:
        for n in num:
            check_list_url("https://news.goo.ne.jp/entertainment/talent/j-"+o+"-"+n+"/")
            #サーバアクセスの間隔を1秒以上空ける
            time.sleep(1)

【結果】
秋山 静香
東 千秋
大山 秋

神山 大和

島 和三
千石 愛
千秋
新奈
新山 大
東 大
大和
大和
和香
和香奈

結構いらっしゃいました.
あまりタレントさんには明るくないので知らない人ばかり.
ああ!確かに!となるのは千秋さん.
データセットを変えれば他にも見つかると思われる.

WaybackMachine

WaybackMachine(https://archive.org/web/)というWebアーカイビングサービス.
Web上のあらゆるサイトのデータを残しておいてくれている.

例えば昔書いていたブログ(http://blog.masahikoyasui.jp/)はOCNのブログサービス閉鎖に伴って現在は見れないけれど,
ここだと残っている.
web.archive.org


もともとは小学校のころ(2000年)の情報の時間に作ったクラスのホームページを探していたが,
そちらはアクセス数が少なかったからか一部しか残っておらず,目的のページは発見できなかった.
友達の書いた記事はかろうじて残っていた.
(30年くらい時代のズレた20世紀少年にありそうなネタ).

ポケモンだいすきクラブとか休み時間にコンピュータ室でみてたなぁ.
web.archive.org