5

繰返し引張試験のデータを分析しています。入力として、x 値と y 値の膨大なリストが使用されます。材料が硬化するか軟化するかを説明するには、各サイクル ループの青色の勾配を取得する必要があります。

引張試験

スロープ

坂の下りは子供祭りですが、上り坂が課題です。

チャレンジ

私はこれまでこのアプローチを行ってきました。各ループの極大値を下回るポイントがいくつかあるループをスライスし、ハードナンバーのポイント数から赤い線を作成しました。赤い線の近似は、交点を取得するために使用されますpoly1d(polyfit(x1,x2,1))。fsolveただし、ポイントの分布は常に同じではないため、常に正しく機能するとは限りません。

問題は、交差する 2 本の線 (赤)の間隔を正しく定義する方法です。上の写真では、平均勾配とともに 3 つの実験が行われています。これが最善のアプローチではないと判断して、ループごとに4つの最も近いポイントを見つけるのに数日を費やしました. 最後に、ここで stackowerflow を終了しました。

必要な出力は、交点のおおよその座標を含むリストです。再生する場合は、ここに曲線のデータがあります (0,[[xvals],[yvals]])。これらは簡単に読むことができます

import csv
import sys
csv. field_size_limit(sys.maxsize)     

csvfile = 'data.csv'
tc_data = {}
for key, val in csv.reader(open(csvfile, "r")):
    tc_data[key] = val
for key in tc_data:
  tc = eval(tc_data[key])

x = tc[0]
y = tc[1]
4

2 に答える 2

6

これは少しやり過ぎかもしれませんが、交差点を見つける適切な方法は、カーブをチャンクに分割したら、最初のチャンクのセグメントが 2 番目のチャンクのセグメントと交差するかどうかを確認することです。

簡単なデータ、長サイクロイドの一部を作成し、y 座標が増加から減少に反転する場所を見つけます。ここと同様です。

a, b = 1, 2
phi = np.linspace(3, 10, 100)
x = a*phi - b*np.sin(phi)
y = a - b*np.cos(phi)
y_growth_flips = np.where(np.diff(np.diff(y) > 0))[0] + 1

plt.plot(x, y, 'rx')
plt.plot(x[y_growth_flips], y[y_growth_flips], 'bo')
plt.axis([2, 12, -1.5, 3.5])
plt.show()

ここに画像の説明を入力

2 つの線分があり、1 つは点P0から へP1、もう 1 つは点から へ向かうQ0場合Q1、ベクトル方程式 を解くことでそれらの交点を見つけることができます。2 つの線分は、との両方が 内にあるP0 + s*(P1-P0) = Q0 + t*(Q1-Q0)場合に実際に交差します。すべてのセグメントでこれを試す:st[0, 1]

x_down = x[y_growth_flips[0]:y_growth_flips[1]+1]
y_down = y[y_growth_flips[0]:y_growth_flips[1]+1]
x_up = x[y_growth_flips[1]:y_growth_flips[2]+1]
y_up = y[y_growth_flips[1]:y_growth_flips[2]+1]

def find_intersect(x_down, y_down, x_up, y_up):
    for j in xrange(len(x_down)-1):
        p0 = np.array([x_down[j], y_down[j]])
        p1 = np.array([x_down[j+1], y_down[j+1]])
        for k in xrange(len(x_up)-1):
            q0 = np.array([x_up[k], y_up[k]])
            q1 = np.array([x_up[k+1], y_up[k+1]])
            params = np.linalg.solve(np.column_stack((p1-p0, q0-q1)),
                                     q0-p0)
            if np.all((params >= 0) & (params <= 1)):
                return p0 + params[0]*(p1 - p0)

>>> find_intersect(x_down, y_down, x_up, y_up)
array([ 6.28302264,  1.63658676])

crossing_point = find_intersect(x_down, y_down, x_up, y_up)
plt.plot(crossing_point[0], crossing_point[1], 'ro')
plt.show()

ここに画像の説明を入力

私のシステムでは、これは 1 秒あたり約 20 の交差を処理できます。これは超高速ではありませんが、時々グラフを分析するにはおそらく十分です。2x2 線形システムの解をベクトル化することで速度を上げることができる場合があります。

def find_intersect_vec(x_down, y_down, x_up, y_up):
    p = np.column_stack((x_down, y_down))
    q = np.column_stack((x_up, y_up))
    p0, p1, q0, q1 = p[:-1], p[1:], q[:-1], q[1:]
    rhs = q0 - p0[:, np.newaxis, :]
    mat = np.empty((len(p0), len(q0), 2, 2))
    mat[..., 0] = (p1 - p0)[:, np.newaxis]
    mat[..., 1] = q0 - q1
    mat_inv = -mat.copy()
    mat_inv[..., 0, 0] = mat[..., 1, 1]
    mat_inv[..., 1, 1] = mat[..., 0, 0]
    det = mat[..., 0, 0] * mat[..., 1, 1] - mat[..., 0, 1] * mat[..., 1, 0]
    mat_inv /= det[..., np.newaxis, np.newaxis]
    import numpy.core.umath_tests as ut
    params = ut.matrix_multiply(mat_inv, rhs[..., np.newaxis])
    intersection = np.all((params >= 0) & (params <= 1), axis=(-1, -2))
    p0_s = params[intersection, 0, :] * mat[intersection, :, 0]
    return p0_s + p0[np.where(intersection)[0]]

はい、面倒ですが、動作し、100 倍高速に動作します。

find_intersect(x_down, y_down, x_up, y_up)
Out[67]: array([ 6.28302264,  1.63658676])

find_intersect_vec(x_down, y_down, x_up, y_up)
Out[68]: array([[ 6.28302264,  1.63658676]])

%timeit find_intersect(x_down, y_down, x_up, y_up)
10 loops, best of 3: 66.1 ms per loop

%timeit find_intersect_vec(x_down, y_down, x_up, y_up)
1000 loops, best of 3: 375 us per loop
于 2013-07-29T18:46:35.047 に答える