少数で表される数値の出現頻度を数えてヒストグラムを作成したいAさんの話(問題)

題名の通りの問題です。問題は以下の通りです。汚文字ですいません。

問題文では少数の最小値と最大値が-2.0と2.0で、それを8つの領域に均等に分割していますが、これは一例にしかすぎません。これら(最小値、最大値、分割数)の数字が変化しても機能するように一般化した関数を求めてほしいです。

私の解答は一番下に載せます。
関数に境界値を与えると動作が予測困難な振る舞いになりますので注意が必要です。

もっと良い方法があるよなど、コメントいただけると助かります。

以上。

f:id:keimina:20181017000415p:plain

解答例 (境界値に弱いので注意):

何が言いたいのかというと`index = int((v - A) / (B - A) * N)`の1行でこの関数は実装できると言うことです。
ただし、テストデータ 2 の検証を見ていただくと分かるように境界値に弱いですので、境界値が入力として来た場合の対策は別途行う必要があります。
境界値が来たら右に少しづらすとかで対処できるかもしれません,

# -*- coding:utf-8 -*-
import numpy as np

# 分割数
N = 8

# 下限
A = -2.0

# 上限
B = 2.0

# 均等に分割した数直線
linespace = np.linspace(A, B, N + 1)

# 数直線を表示

print("================ 数直線 ================")
print("\n|\n".join(linespace.astype(str)))


################################################################
# Notice: len(linespace) have to be greater than 2.
################################################################

# テストデータ 1
print("")
print("================ テストデータ1 ================")
middle_points_of_linespace = np.array([(x2 + x1)/2 for x1, x2 in zip(linespace, linespace[1:])])
print("\n|\n".join(middle_points_of_linespace.astype(str)))

# テストデータ 2 境界値のデータ(関数が境界値でどのような振る舞いをするか調べるために最後に使用する)
print("")
print("================ テストデータ2 ================")
edge_points_of_linespace = linespace[:]
print("\n|\n".join(edge_points_of_linespace.astype(str)))


# 数直線とテストデータの関係(期待値)を表示
print("")
print("================数直線と テストデータ 1 の関係================")
print("\n|\n".join("{}\n|\n(テストデータ{})".format(i, mi) for i, mi in zip(linespace, middle_points_of_linespace)))
print("|\n{}".format(linespace[-1]))

# 変換関数の処理(面倒なので関数化していない)
# テストデータを一つづつ読み込んで、プリントする
print("")
print("================ 検証 テストデータ1 ================")
for v in middle_points_of_linespace:
    # 数値をインデックスに変換する(1行で終わり)
    index = int((v - A) / (B - A) * N)
    # 以下は数値が linespace の範囲外にあった場合の対処
    if index < 0:
        print("変換先がlinespaceのインデックスの範囲外です")
    elif index > len(linespace) - 1:
        print("変換先がlinespaceのインデックスの範囲外です")
    # テストデータの値 と linespace の値の関係を プリントする
    else:
        print("{} < {} < {}".format(linespace[index], v, linespace[index+1]))


# 境界値のテストデータでテストすると予測しにくい動きとなる
print("")
print("================ 検証 テストデータ2 ================")
for v in edge_points_of_linespace:
    # 数値をインデックスに変換する(1行で終わり)
    index = int((v - A) / (B - A) * N)
    # 以下は数値が linespace の範囲外にあった場合の対処
    if index < 0:
        print("変換先がlinespaceのインデックスの範囲外です")
    elif index > len(linespace) - 2:
        print("変換先がlinespaceのインデックスの範囲外です")
    # テストデータの値 と linespace の値の関係を プリントする
    else:
        print("{} < {} < {}".format(linespace[index], v, linespace[index+1]))



出力結果:

================ 数直線 ================
-2.0
|
-1.5
|
-1.0
|
-0.5
|
0.0
|
0.5
|
1.0
|
1.5
|
2.0

================ テストデータ1 ================
-1.75
|
-1.25
|
-0.75
|
-0.25
|
0.25
|
0.75
|
1.25
|
1.75

================ テストデータ2 ================
-2.0
|
-1.5
|
-1.0
|
-0.5
|
0.0
|
0.5
|
1.0
|
1.5
|
2.0

================数直線と テストデータ 1 の関係================
-2.0
|
(テストデータ-1.75)
|
-1.5
|
(テストデータ-1.25)
|
-1.0
|
(テストデータ-0.75)
|
-0.5
|
(テストデータ-0.25)
|
0.0
|
(テストデータ0.25)
|
0.5
|
(テストデータ0.75)
|
1.0
|
(テストデータ1.25)
|
1.5
|
(テストデータ1.75)
|
2.0

================ 検証 テストデータ1 ================
-2.0 < -1.75 < -1.5
-1.5 < -1.25 < -1.0
-1.0 < -0.75 < -0.5
-0.5 < -0.25 < 0.0
0.0 < 0.25 < 0.5
0.5 < 0.75 < 1.0
1.0 < 1.25 < 1.5
1.5 < 1.75 < 2.0

================ 検証 テストデータ2 ================
-2.0 < -2.0 < -1.5
-1.5 < -1.5 < -1.0
-1.0 < -1.0 < -0.5
-0.5 < -0.5 < 0.0
0.0 < 0.0 < 0.5
0.5 < 0.5 < 1.0
1.0 < 1.0 < 1.5
1.5 < 1.5 < 2.0
変換先がlinespaceのインデックスの範囲外です

指数移動平均(ema)について考察する

指数移動平均(Exponential Moving Average 以下 ema)の仕組みや動作がよくわからないので、計算式をPythonで実装し、関数に適用してみます。
emaの計算式は以下のWikipediaのページを参考にしました。

https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average



これによると、emaの式は以下のように表されるようです。
SとYをtの関数とし、SをYのemaとすると以下のように表せます。

・t = 1 のとき
 S(t) =Y(1)

・t > 1 のとき
 S(t) = a * Y(t) + (1 - a) * S(t-1)

a の取りうる値は 0.0 〜 1.0 で、
ema は a の値が 1.0 のとき現在の Y(t) の値そのものに、 0.0 のとき Y(0) の値そのものになります。

ここでは、 a の値を変化させたときの ema の値を見てみます。さらに、過去のデータを直近のものに限定した場合のemaの値も見てみます。

ソースコード

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider

def St_a(y):
    return y

def St_b(alpha, St_minus, Yt):
    return alpha * Yt + (1.0 - alpha) * St_minus

def S(alpha, y):
    ret = []
    for i in range(150):
        if i==0:
            ret.append(St_a(y[i]))
        else:
            ret.append(St_b(alpha, ret[i-1], y[i]))
    return ret

N = 3
alpha_0 = 0.5
x = np.arange(150)
y = 2*np.sin(N*2*np.pi/150*x) + np.random.randn(150)
s = S(alpha_0, y)

fig, (ax, ax_slider_1, ax_slider_2) = plt.subplots(3, 1, gridspec_kw = {'height_ratios':[14, 1, 1]})
ax.plot(x, y, linestyle='--')
ax.plot(x, 2*np.sin(N*2*np.pi/150*x), linestyle=':')
ax_s, = ax.plot(x, s)

sldr_1 = Slider(ax_slider_1, 'alpha', 0.0, 1.0, valinit=alpha_0)
sldr_2 = Slider(ax_slider_2, 'start_index', 0, 149, valinit=0, valfmt="%d")

def update_s(alpha_val):
    ax_s.set_xdata(x)
    ax_s.set_ydata(S(alpha_val, y))
    sldr_2.set_val(0)
    fig.canvas.draw_idle()

def update_s_start(s_start_index):
    s_start_index = int(s_start_index)
    tmp_y = y.copy()
    tmp_y[:s_start_index] = tmp_y[s_start_index]
    ax_s.set_xdata(x[s_start_index:])
    ax_s.set_ydata(S(sldr_1.val, tmp_y)[s_start_index:])
    fig.canvas.draw_idle()

sldr_1.on_changed(update_s)
sldr_2.on_changed(update_s_start)
    
plt.show()


実行結果:

f:id:keimina:20180825220145g:plain

xが横軸, yが縦軸です。

オレンジ色の点線はSin波を表しています。
ema を適用する前の オリジナルの関数は,このオレンジ色のSin波にランダムな値を加えたもので、青色の点線です。
それに ema を適用したのが緑色の線になります。

グラフの下にある スライダーの alpha が 数式の a を表しています。
スライダーの start_index は ema の初期値を表しています。このスライダーは過去のデータを直近のものに限定するために使用します。

まず alpha を動かして見ます。
・a = 0 のときは ema はx軸に関係なく常に初期値を維持し続けます。
・a = 1.0 のときは、青い点線と緑色の線が重なりますので、オリジナルの関数と同じになることがわかります。

次に start_index を動かしてみます。
・a = 0.0 にして start_index をプラス方向に動かすと緑色の線がy軸方向に大きく変動することがわかります。
・a = 1.0 にして start_index をプラス方向に動かしても緑色の線がy軸方向に変動することはないことがわかります。
・a = 0.4 にして start_index をプラス方向に動かすと緑色の線がy軸方向に少しだけ変動します。
・a = 0.05 にして start_index をプラス方向に動かすと緑色の線がy軸方向に大きく変動します。

このことから、 a が 1.0 に近ければ 、大雑把にいうと単純移動平均のように動作することがわかりました。
また、start_index を動かした時のグラフの変化を見ると、データ量が少ない場合はa が0に近いと初期値の影響を大きく受けてしまうため、ノイズの除去をして低周波成分を見たいときなどにはデータ量が必要だということがわかりました。そのためema を 単純移動平均と同じような場面で使うときには a を 1.0 に近付けて使うのが良いということがわかりました。emaは同一のシンプルなアルゴリズムでデータ数に関係なく過去の全ての情報を含む値を算出することができるので便利だと思いました。

以上。

TensorFlow で「y=x*(x-1)*(x-2)」の極小値を求める

TensorFlow を利用し関数の局所的に凹になってる箇所ー Local Minima ーすなわち数学的にいう極小値を求めてみようと思います。今回、対象とする関数は 「y = x*(x-1)*(x-2)」です。グラフにすると以下のようになると思います。汚くてすいません。

f:id:keimina:20180822204257j:plain

ちなみにこの関数は x = 1 ± 1/√3 (約1.577, 約0.422)のとき極小値と極大値をとります。極小値と極大値はy=f(x)のグラフの傾きが0となる箇所を求めればよいです。すなわち、関数を微分したものをf'(x)とすればf’(x)=0となるxを求めればよいことになります。試しに Sympy で確かめてみます。

from sympy import Symbol, solve
x = Symbol('x')
y = x*(x-1)*(x-2)
print(solve(y.diff(x)))
# output:
# [-sqrt(3)/3 + 1, sqrt(3)/3 + 1]

Sympy により、すでに答え(x =1.577)は出てしまったのですが、これを TensorFlow を用いて数値計算的に求めることが今回のゴールです。(上記 Sympy の例では 記号的なアルゴリズムで求めていると思われるため今回の TensorFlow で求めるやり方と内部の動きが異なります。)

ここからTensorFlow を用いて最急降下法(Gradient Descent)によって関数の極小値を見つけます。TensorFlow では、最急降下法により関数の極小値を求めるための「TensorFlow.train.GradientDescentOptimizer」クラスとその「minimize」オペレーションが用意されています。今回は、それと同様のオペレーションのある「TensorFlow.train.MomentumOptimizer」クラスを使用し、得られたxとyの値をグラフにプロットし、画像として保存し得られたxとyを見てみます。

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt

# 関数 y = x*(x-1)*(x-2) を定義
def f(x):
    return x * (x - 1.0) * (x - 2.0)

def save_graph(solved_x, filename):
    # グラフ描画して保存
    # x軸とy軸を描画
    plt.axis('off')
    plt.axhline(0, color='black', linewidth=1)
    plt.axvline(0, color='black', linewidth=1)
    # (x, y) = (0, 0) に点を描画
    plt.plot(0, 0, marker='o', color='black')
    plt.text(0, 0, "(0, 0)")
    # (x, y) = (1, 0) に点を描画
    plt.plot(1, 0, marker='o', color='black')
    plt.text(1, 0, "(1, 0)")
    # (x, y) = (2, 0) に点を描画
    plt.plot(2, 0, marker='o', color='black')
    plt.text(2, 0, "(2, 0)")
    # x軸にラベルを追加
    plt.text(2.3 + 0.2, 0, "x", verticalalignment='center')
    # y軸にラベルを追加
    plt.text(0, f(2.3)+0.2, "y", horizontalalignment='center')
    # 関数 y = f(x) を x = [-0.3, 2.3] の範囲で描画
    a = np.linspace(-0.3, 2.3, 100)
    b = f(a)
    plt.plot(a, b, label='y = x*(x-1)*(x-2)')
    plt.legend(loc=1)
    # 極小値の描画
    plt.plot(solved_x, f(solved_x), marker='o', color='orange')
    # f(x)=0 の点をわかりやすくするための線を描画
    plt.axvline(solved_x, linestyle='--', color='orange')
    plt.axhline(f(solved_x), linestyle='--', color='orange')
    plt.text(solved_x, f(solved_x)-0.1, "(x, y)=(%f, %f)"%(solved_x, f(solved_x)), horizontalalignment='center')
    plt.savefig(filename)
    plt.close()


################################################################
# メイン
################################################################
x_kyokusyou = 1.0 - 1.0/np.sqrt(3) # 極大値(local maxima)をとるxの値
#x_start = x_kyokusyou - 0.01 # パターン1
#x_start = x_kyokusyou + 0.01 # パターン2
x_start = 2.2 # パターン3
x = tf.Variable(x_start)

y = x * (x - 1.0) * (x - 2.0)
#grad_opt = tf.train.GradientDescentOptimizer(0.03)
grad_opt = tf.train.MomentumOptimizer(0.01, 0.95)
grad_opt_minimize = grad_opt.minimize(y)
init = tf.global_variables_initializer()

with tf.Session() as sess:
    sess.run(init)
    for i in range(200):
        sess.run(grad_opt_minimize)
        solved_x = sess.run(x)
        if solved_x < -0.3: # -infに落ちるので途中でやめる
            break
        # グラフ描画して保存
        save_graph(solved_x, "/Users/kei/Desktop/pic/%03d.png"%i)

初期値を変えて3パターン(パターン1、パターン2、パターン3)で試し、出力されたグラフを動画にしてみました。

・パターン1(初期値 x=1-1/√3 - 0.01)
f:id:keimina:20180822205755g:plain

・パターン2(初期値 x=1-1/√3 + 0.01)
f:id:keimina:20180822205848g:plain

・パターン3(初期値 x=2.2)
f:id:keimina:20180822210005g:plain

パターン1では求めた点がx軸のマイナス方向へ向かって凹があるまで移動していきます。この関数は x < 1 - 1/√3 では凹となる箇所がないため点が無限に移動し続けます。パターン2、パターン3は振り子のように減衰しながら極小値に収束していくように移動します。

最後に、パターン3で求めたxを見て見ます。

print(solved_x)
# output:
# 1.5798204

多少の誤差はありますが、現在の繰り返し回数を200回から増やせば誤差はなくなっていくものと思われます。

以上より、最急降下法によって関数の極小値を数値計算的に求められることがわかりました。
今回の例ではパラメータが一つなため計算量が少なくCPUでもそんなに時間もかからず求まりますが、、ニューラルネットワークの場合は計算量が層の深さの累乗のオーダーで大きくなっていきます。ただ、計算自体は1度、計算式を組み立ててしまえば、その計算式にひたすら数値を代入して演算を繰り返してニューラルネットワークを更新するだけです。そのためGPUのように安く大量の数値計算の並列処理ができるデバイスとこの方法を掛け合わせることによりディープラーニングが実用的なものになり、それが現在の第3次人口知能ブームのきっかけなのではないかと思っています。

以上です。おやすみなさい。

ニューラルネットワークで"Z = x*x + y*y"を学習させてみた

いきなりですが、関数 "Z=x*x + y*y" の様な曲面をニューラルネットワークで学習させてみます。ディープラーニングの話です。

目的:多層ニューラルネットワークが色々な関数を近似できるとは聞いていたが、それが本当かどうか確かめる。多層ニューラルネットワークを学習させた関数のグラフが関数"Z=x*x + y*y"と同様な曲面になるか確かめる。

注意:細かい話ですが、実際の関数は"Z=k(x*x + y*y) + c"(kとcは定数)です。xとyの範囲は[-1,1]です。

・モデルの構築
ニューラルネットワークのモデルは以下の条件1、2を満たす様に構築しました。

(条件1)「3層以上であること(全て全結合層)」
(条件2)「1層のニューロンの数が2以上であること」

(条件1の理由):3層以上でないと関数の非線形性を表現しきれないから
(条件2の理由):独立変数x,yを一つの変数で現すことは不可能、よって1層のニューロンの数は入力の独立変数の数より多くなければならないから

f:id:keimina:20180324025738p:plain

構築したネットワークを連立方程式で書くと以下の様になります。
今回は活性化関数としてtanh(ハイパボリックタンジェント)を用いることとしました。

# layer 1
x11 = tanh(a*x + b*y + c)
x21 = tanh(d*x + e*y + f)

# layer 2
x13 = tanh(g*x11 + h*x21 + i)
x23 = tanh(j*x11 + k*x21 + l)

# layer 3(output)
z= m*x13 + n*x23 + o

範囲「-1<x<1, -1<y<1」において、目的(Z=x*x + y*y)の関数と出力が同じになる様にパラメータを学習させます。パラメータの最適化には勾配下降法を用いています。

・実装
参考として、TensorFlow のコードを載せておきます。

import numpy as np
import tensorflow as tf
import random
from itertools import product
import time
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import math

def Mokuteki(x1, x2):
    return x1**2 + x2**2

k1 = math.sqrt(15)

I = tf.placeholder(tf.float32, (1, 2))
W1 = tf.Variable(k1*(1-2*tf.random_uniform((2, 2))))
B1 = tf.Variable(k1*(1-2*tf.random_uniform((1, 2))))

I2 = tf.nn.tanh(tf.matmul(I, W1) + B1)

W2 = tf.Variable(k1*(1-2*tf.random_uniform((2, 2))))
B2 = tf.Variable(k1*(1-2*tf.random_uniform((1, 2))))
I3 = tf.nn.tanh(tf.matmul(I2, W2) +  B2)

W3 = tf.Variable(k1*(1-2*tf.random_uniform((2, 1))))
B3 = tf.Variable(k1*(1-2*tf.random_uniform((1,))))
Y = tf.matmul(I3, W3) + B3

T = tf.placeholder(tf.float32)
L = tf.square(T - Y)
train_step = tf.train.GradientDescentOptimizer(0.03).minimize(L)

#training data
N = 16
x1_range = np.linspace(-1,1,N, dtype=np.float32)
x2_range = np.linspace(-1,1,N, dtype=np.float32)
train_x1, train_x2 = np.meshgrid(x1_range, x2_range)
train_y = np.array([Mokuteki(val1, val2) for val1, val2 in zip(train_x1.flatten(), train_x2.flatten())], dtype=float).reshape(N,N)
train_y -= np.mean(train_y)
train_y /= np.std(train_y)

init = tf.global_variables_initializer()

sess = tf.Session()
sess.run(init)

start_time = time.time()
for n_epoch in range(100):
    sample_all = random.sample(list(zip(train_x1.flatten(), train_x2.flatten(), train_y.flatten())), N*N)
    for x1v, x2v, tv in sample_all:
        sess.run(train_step, feed_dict={I:[[x1v,x2v]], T:tv})

    total_err = 0.0
    for x1v, x2v, tv in sample_all:
        total_err +=  sess.run(L, feed_dict={I:[[x1v, x2v]], T:tv})
    print("%d epoch:"%n_epoch, "loss:%f"%(total_err/N/N))

print("%d sec"%(time.time()-start_time))
print("finished")

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(train_x1, train_x2, np.array([sess.run(Y, feed_dict={I:[[x1v, x2v]]}) for x1v, x2v in zip(train_x1.flatten(), train_x2.flatten())], dtype=float).reshape(N,N))
plt.show()

sess.close()


やっていることは、まず、範囲「-1<x<1, -1<y<1」における Z = x*x + x*x の集合を訓練データとして生成し、あらかじめシャッフルしておく。ニューラルネットワークの出力を表す関数と、訓練データの2乗誤差をとる新たな関数Lを定義。勾配下降法により、その関数Lを最小にする様にパラメータを調整。
注意:表示される"loss:"は2乗誤差の値の平均値が記載されます。

最後に、6000エポックまでの学習中の関数のグラフの変化を動画にしました。
関数のグラフが、平衡状態から崩れて別の形になったりと興味深いです。

youtu.be

パラメータの初期値と学習係数を変えて色々と試してみましたが、ソースコードにあるもので落ち着いています。一番困ったのが、局所解で平衡状態になって抜け出せなくなってしまう場合です。局所解におちいらない様な設定値の作り方やアルゴリズムを知っていたら教えてください。

以上です。

numpy.sumのaxisを指定したときの挙動が覚えられないのでイメージ化してみた

numpyの関数にはnumpy.sumやnumpy.anyなど引数にaxisを指定できる場合があるのだけど、いまいちイメージがつかめなかったのでnp.sumのaxis指定したときの挙動のイメージをお絵描きしました。2次元の行列に限定した話です。以下の図のように「axis=0のときは踏みつぶす」、「axis=1のときは握り潰す」とイメージすることで解決です。あれ、踏み潰すのはどっちだっけ。それでは、おやすみなさい。

 

f:id:keimina:20180302222915j:image

エクセルのアルファベットの列番号を数値に変換するPythonプログラム

昼休みの暇つぶしに書いたコードを載せます。
エクセルのアルファベットの列番号を数値に変換する簡単なPythonプログラムになります。

例:
"A"は1
"Z"は26
"AA"は27

アルファベットを数値に変換し、各桁の重みをかけて右から左へ順に足しこんで行くプログラムです。
以下にアルファベットを数値に変換する関数conv(alphabet)とその関数のテストプログラムを示します。

inps = ['A', 'Z', 'AA', 'AZ', 'ZZ']
out_expects = [1, 26, 26+1, 26+26, 26*26+26]

def conv(alphabet):
	total = 0
	for n, i in enumerate(alphabet[::-1]):
		chr_num = ord(i.lower())
		chr_A_num = ord('a')
		total +=  (1 + chr_num - chr_A_num)* (26**n)
	return total
		
def test(inps, out_expects):
	for alphabet, out_expect in zip(inps, out_expects):
		out = conv(alphabet)
		print('input:%s, output:%s, out expect:%s'%(alphabet, str(out), str(out_expect)))
		print('%s' %('OK' if out==out_expect else 'NG'))
		
test(inps, out_expects)

テスト結果

input:A, output:1, out expect:1
OK
input:Z, output:26, out expect:26
OK
input:AA, output:27, out expect:27
OK
input:AZ, output:52, out expect:52
OK
input:ZZ, output:702, out expect:702
OK

昆布関数conv(alphabet)がアルファベットを数値に変換します。
出力結果の期待値については以下のように考えるとわかりやすいと思います。

'A'は1
'AA'は26+1
'AZ'は26+26
'BA'は26*2 + 1
'BZ'は26*2 + 26
、、、
'ZZ'は26*26+26

もちろん「conv(alphabet) - 1」とすれば0から始まるインデックスとして使うこともできます。
2進数変換はよく聞くけどこの形の変換は聞いたことがないです。
この形式の数はなんというのでしょうかね?16進数のように数値とアルファベットが混ざっているわけではないのでN進数という言葉ではない気がします。知っている人が入ればコメントください。


以上

2つの線の交点を表す数式を作ってpython matplotlibで検証

XY座標に点A1(Xa1, Ya1), A2(Xa2, Ya2), B1(Xb1, Yb1), B2(Xb2, Yb2)があったとして、

線A1A2 と 線B1B2の交点を数学的に求め、数式があってるか確かめようという話。

以下のように数学的に求めます。

  • まず交点の式をだします。

f:id:keimina:20171014081940j:plain

  • 次に上画像の一番下に書いてある方程式を解きます。

f:id:keimina:20171014081906j:plain

  • 逆行列をどうやって求めるかは、numpyライブラリを使わずに数学公式を使います。

行列Mとその逆行列invMの関係は以下になります。

M = ((a,b),(c,d))
invM = 1/(a*d - b*c) * ((d, -b),(-c, a))

参考:
2x2行列と3x3行列と4x4行列の逆行列の公式

  • 準備が整ったのでpythonで実装します。

交点を求める自作関数get_line_intersectionで上記で求めた数式を使っています。

from numpy import sin, cos, pi
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider


class point:
    def __init__(self, x, y):
        self.x = float(x)
        self.y = float(y)

def get_line_intersection(A1, A2, B1, B2):
    a = A2.y - A1.y
    b = -B2.y + B1.y
    c = A2.x - A1.x
    d = -B2.x + B1.x
    C1 = B1.y - A1.y
    C2 = B1.x - A1.x
    
    tmp = a * d - b * c
    if tmp:
        invMa = d  / tmp
        invMb = -b / tmp
        invMc = -c / tmp
        invMd = a  / tmp
        
        m = invMa * C1 + invMb * C2
        n = invMc * C1 + invMd * C2

        if 0<=m<=1 and 0<=n<=1:#両方の線上に点がある場合は交点の座標を返す
            return point(A1.x + m * (A2.x - A1.x), A1.y + m * (A2.y - A1.y))
        else:
            return point(0, 0)#辺の延長線上に交点がある場合は(0,0)を返すこととする
    else:
        return point(0, 0)#逆行列が存在しない時は強制的に(0,0)を返すことにする

################################################################
# 描画
################################################################
fig, ax = plt.subplots()
ax.set_xlim(-11,29)
ax.set_ylim(-11,29)
ax.set_aspect('equal')

plt.subplots_adjust(left=0.25, bottom=0.35)
l1, = plt.plot([6, 6], [6 + 5 * cos(pi/2), 6 + 5 * sin(pi/2)])
l2, = plt.plot([12, 12], [12 + 5 * cos(pi/2), 12 + 5 * sin(pi/2)])
l3, = plt.plot([0],[0], marker="o")


axcolor = 'lightgoldenrodyellow'
axRA = plt.axes([0.25, 0.1, 0.65, 0.03], facecolor=axcolor)
axTA = plt.axes([0.25, 0.15, 0.65, 0.03], facecolor=axcolor)

axRB = plt.axes([0.25, 0.2, 0.65, 0.03], facecolor=axcolor)
axTB = plt.axes([0.25, 0.25, 0.65, 0.03], facecolor=axcolor)


RA = Slider(axRA, 'RadiusA', 0.0, 20.0, valinit=5.0)
TA = Slider(axTA, 'ThetaA', 0.0, 4 * pi, valinit=pi/2)

RB = Slider(axRB, 'RadiusB', 0.0, 20.0, valinit=5.0)
TB = Slider(axTB, 'ThetaB', 0.0, 4 * pi, valinit=pi/2)

def update(val):
    l1.set_xdata([6, 6 + RA.val * cos(TA.val)])
    l1.set_ydata([6, 6 + RA.val * sin(TA.val)])
    l2.set_xdata([12, 12 + RB.val * cos(TB.val)])
    l2.set_ydata([12, 12 + RB.val * sin(TB.val)])
    p = get_line_intersection(point(6,6), point(6+RA.val*cos(TA.val), 6+RA.val*sin(TA.val)), point(12,12), point(12+RB.val*cos(TB.val), 12+RB.val*sin(TB.val)))
    l3.set_xdata([p.x])
    l3.set_ydata([p.y])
    fig.canvas.draw_idle()

RA.on_changed(update)
TA.on_changed(update)
RB.on_changed(update)
TB.on_changed(update)

plt.show()
  • 実行結果

youtu.be

  • 計算時間

CPU Core i7 2.0Ghzで計算したが計算時間は100万回やっても5秒くらいらしい

100000 loops, best of 3: 4.76 µs per loop