ふと思い出した

DLのプログラムを書いていて、ふと大学時代にmatlabで書いていた幾何問題を解く問題（※）に似ているな、と思い、改めてpythonで簡単な問題を解いてみた。

※お車の話。最後に余談として記載。

最新MATLABハンドブック第六版

作者: 小林一行
出版社/メーカー: 秀和システム
発売日: 2018/04/05
メディア: Kindle版
この商品を含むブログを見る

二次元平面で2つの固定長の線の軌跡

一体なんなのか。

こんな線ACと線BCがあったとする。

f:id:kazuhitogo:20190211115918p:plain

線は剛体で固定長だとする。

そのとき、Bを動かしたとき、 Cはどこに行くの？

f:id:kazuhitogo:20190211115901p:plain

という問題。

ちなみに数学的にこの問題を解こうとすると、 Aを中心とする長さACの円と、 B'を中心とする長さBCの円の交点を求めればよいので、下記方程式を解けばよい…のだが、それなりにめんどくさい。

$\displaystyle ( x - A_x ) ^ 2 + ( y - A_y ) ^ 2 = (\vec{AC}) ^ 2$

$\displaystyle ( x - B ' _ x ) ^ 2 + ( y - B ' _ y ) ^ 2 = (\vec{BC}) ^ 2$

というかこのような非線形な方程式は私には解けない。で、それを現実的にどう解くか…というのが今回のエントリ。

ちなみにさきほどの方程式の解は基本的に２つある（２つの円の交点は基本２つ）が、現実世界においてはBをB'に移動したらCはCに近い方の解にしか移動しないため、そちらの座標を計算する。

どうやるのか

ある程度の誤差を許容しつつ（現実世界では0.0001mmの誤差くらいは許容、というより製造限界)、コンピュータは簡単な問題を繰り返し解くのが得意なので、その方法で解く。

ニュートンラフソン法というのがある。 ja.wikipedia.org

簡単にいうと複雑な問題（今回で言えば二元二次方程式）を、線形の連立方程式として解く。
複雑な問題を簡易に問いているので当然結果には誤差が出る。
解き終わった後誤差を評価し、また同じ方法で解く…を繰り返し、ある一定の誤差より低くなったらやめる

というものである。

1の連立方程式として解くところがミソであるのだが、 Cの座標を微小に動かしたときに、それぞれの線の長さがどう変わるか、を算出し、連立方程式化する。

で、DLやってて思い出した理由がここにあって、誤差を評価して微小（＝つまり微分）に動かして修正して、を繰り返して行くのが、まさにDLの発想だな、と思い出した次第。

改めて図示すると以下の通り。

初期位置

f:id:kazuhitogo:20190211115918p:plain

Bを移動

f:id:kazuhitogo:20190211185659p:plain

Cはどこに行く？

f:id:kazuhitogo:20190211185732p:plain

Bを動かすと赤に誤差が生じている

f:id:kazuhitogo:20190211185815p:plain

CをX方向に微小に動かして線分AC,BCの長さがどうなるかを求める

f:id:kazuhitogo:20190211185853p:plain

CをY方向に微小に動かして線分AC,BCの長さがどうなるかを求める

f:id:kazuhitogo:20190211185917p:plain

そこからは連立方程式を立てて解きましょう…。

コード

from matplotlib import pyplot as plt
import math
import numpy as np

# point A,B,C
pa = np.array([1,1])
pb =  np.array([3,-1])
pc = np.array([6,3])

# 元のデータ(描画用)
pa_org = np.array([1,1])
pb_org =  np.array([3,-1])
pc_org = np.array([6,3])


# 2点間の距離を求める
# a = np.array([x1,y1])
# b = np.array([x2,y2])
def len_2p(a,b):
    return np.linalg.norm(b-a)

# 2つの誤差の絶対値の合計を返す（誤差関数とする）
def e_2(a,b):
    return np.abs(a) + np.abs(b)

# 描画
def plt3p(a,b,c):
    x = np.array([a[0],b[0],c[0]])
    y = np.array([a[1],b[1],c[1]])
    plt.gca().set_aspect('equal', adjustable='box')
    plt.plot(x,y,"-*")

# 初期の長さ(i=init)
lab_i = len_2p(pa,pb)
lbc_i = len_2p(pb,pc)

# bをx方向に1動かす(動かしたときにcがどこに行くかを計算)
pb = pb + np.array([1,0])

# 移動後長さ(a=after)
lac_a = len_2p(pa,pc)
lbc_a = len_2p(pb,pc)

# 誤差(e=error)
lac_e = lac_i - lac_a
lbc_e = lbc_i - lbc_a

# 誤差計算
e = e_2(lac_e,lbc_e)

# 微小に動かす単位(今回はx,yともに1動かす)
unit_d = 1

# cが繰り返しの演算の中でどう動いていったかを残す配列
pc_hist = np.empty((2,0))
pc_hist = np.append(pc_hist,pc.reshape(2,-1),axis=1)

# 誤差が一定の数値以下になるまで回し続ける
# 解けない問題（先の連立方程式だと解が無いくらい離した場合など）は無限ループになるので気をつけるべし
while(e > 0.000000001): 
    # x方向に微小に動かしたpc
    pc_dx = pc + np.array([unit_d,0])
    # y方向に微小に動かしたpc
    pc_dy = pc + np.array([0,unit_d])
    
    # x方向に微小に動かしたpcの長さ評価
    lac_dx = len_2p(pa,pc_dx)
    lbc_dx = len_2p(pc_dx,pb)

    # y方向に微小に動かしたpcの長さ評価
    lac_dy = len_2p(pa,pc_dy)
    lbc_dy = len_2p(pc_dy,pb)

    # 連立方程式を解く
    # lac_dx * x + lac_dy * y  = lac_e * unit_d
    # lbc_dx * x + lbc_dy * y = lbc_e * unit_d
    # を解く
    move = np.linalg.solve(np.array([[lac_dx,lac_dy],[lbc_dx,lbc_dy]]),np.array([lac_e,lbc_e])) * unit_d
    pc = pc + move
    pc_hist = np.append(pc_hist,pc.reshape(2,-1),axis=1)
    
    # 移動後長さ(a=after)
    lac_a = len_2p(pa,pc)
    lbc_a = len_2p(pb,pc)
    
    # それぞれの線の長さの誤差(e=error)
    lac_e = lac_i - lac_a
    lbc_e = lbc_i - lbc_a
    
    # 先の長さ２つを誤差関数に入れる
    e = e_2(lac_e,lbc_e)

# 描画
plt3p(pa_org,pc_org,pb_org)
plt3p(pa,pc,pb)
plt.plot(pc_hist[0],pc_hist[1])