2013-07-29 30 views
5

Döngüsel çekme testlerinden veri analiz ediyorum. Giriş olarak büyük x ve y değerlerinin listesi kullanılır. Malzemenin sertleşip yumuşamadığını açıklamak için, her devir döngüsünün mavi eğimini almam gerekir. Xy veri noktası grafiğinin tüm kesişim noktalarını uyuşturabilir misiniz?

tensile_test

yamacın alt noktası alınıyor

slope

çocuk festival, fakat üst bir, o meydan okumadır. Her döngünün yerel maksimum aşağıda birkaç puanla döngüler dışarı dilimleme ve kırmızı çizgiler yaparak bugüne kadar bu yaklaşımı yaptık

the_challenge

puan sayımı hardnumbered. Kırmızı çizgilerin tahmini, poly1d(polyfit(x1,x2,1)) ve fsolve ile kesişme noktasını elde etmek için kullanılır. Ancak, her zaman doğru şekilde çalışmaz, çünkü noktaların dağılımı her zaman aynı değildir.

sorun doğru kesişen çizgiler aralığını iki (kırmızı) tanımlar nasıl olup. Yukarıdaki resimde, ortalama eğim ile birlikte 3 deney bulunmaktadır. Her döngü için en yakın 4 noktayı bulmaya çalışarak birkaç gün geçirdim. Bu, en iyi yaklaşım değil. Ve son olarak, burada stackowerflow'ta bitti.

İstenen çıktı, kesişme noktalarının yaklaşık koordinatlarına sahip listedir - eğer oynamak istiyorsanız, here eğri için verilerdir (0, [[xvals], [yvals]]). Sayılan bu kolayca aynı x değerlerinde tüm hatların y-değerlerini yeniden örneklemek için SciPy gelen interp1d işlevini kullanarak çok basit bu yapabileceğini

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] 
+0

benim için çalışıyor. – gggg

+1

Matplotlib –

cevap

6

Bu biraz overkill ancak kesişim noktasını bulmak için uygun bir yol olabilir, ilk yığın herhangi segmenti herhangi segmentte kesiştiğini görmek için ikinci yığıntan.

kendime biraz kolay veri bir prolate cycloid bir parça yapacağım ve here benzer azalan ulaştığı artış döndürür y koordinatı nerede yerleri bulmak için gidiyorum:

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() 

enter image description here

P0 noktasından P1 noktasına ve Q0Q1 noktasına giden iki segmentiniz varsa, vektör denkleminiçözerek kesişim noktalarını bulabilirsiniz.ve iki bölüm aslında hem s hem de t[0, 1] içeriyorsa kesişir.Tüm segmentler için bu dışarı çalışıyorum: Benim sistemde

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() 

enter image description here

, bu muhtemelen yeterli yaklaşık 20 süper değil saniyede kesişimleri, ancak işleyebilir her şimdi ve sonra grafikleri analiz etmek. Sen 2x2 lineer sistemlerin çözümü Vektörizasyonu ortalığı spped mümkün olabilir:

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]] 

Evet, dağınık, ama çalışır ve böylece x100 kat daha hızlı yapar: Bağlantılarınızın ne

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 
+0

Jaime ile yakınlaştırma efektini nasıl elde ettiğinizi merak ediyorum, size bira sandığını borçluyum. Bana adresini yolla ve hemen göndereceğim! – ptaeck

+0

TypeError alıyorum: int satır içi 'intersection = np.all ((params> = 0) & (params <= 1), eksen = (- 1, -2)) '' find_intersect_vec' için gerekli sikloid örnek veriler. – ptaeck

+0

Hangi sürümleri kullanıyorsunuz? Eğer daha önceki bir sürümü kullanıyorsanız, 'np.all' için iki çağrı yapmak zorunda kalabilirsiniz, ben kesişim = np.all (np.all ((params> = 0) & (params) <= 1), eksen = -1), eksen = -1) '1.6'da aynı şeyi yapmalıdır. – Jaime