Individuare tutte le intersezioni del grafico del punto di dati xy con numpy?

Question

Sto analizzando i dati da prove di trazione ciclica. Come input vengono utilizzati enormi elenchi di valori xey. Per descrivere se il materiale si indurisce o si ammorbidisce, ho bisogno di ottenere la pendenza blu di ciascun ciclo.

Ottenere il punto più basso della pista è il festival del bambino, ma quello superiore, che è la sfida.

Ho fatto questo approccio finora, affettare le anse con pochi punti indicati massimo locale di ogni ciclo e rendendo le linee rosse da hardnumbered conteggio di punti. La approssimazione delle linee rosse viene effettuata da poly1d(polyfit(x1,x2,1)) e viene utilizzato fsolve per ottenere il punto di intersezione. Tuttavia non funziona sempre correttamente, perché la distribuzione dei punti non è sempre la stessa.

Il problema è come correttamente definire l'intervallo di due (rosso) linee che si intersecano. Nell'immagine sopra sono 3 esperimenti insieme a una pendenza media. Ho trascorso alcuni giorni cercando di trovare 4 punti più vicini per ogni ciclo, decidendo che questo non è l'approccio migliore. E alla fine, sono finito qui su Stackowerflow.

L'output desiderato è un elenco con le coordinate approssimative dei punti di intersezione: se si desidera riprodurre, here sono dati per la curva (0, [[xvals], [yvals]]). Theese può essere facilmente letta con

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] 

Answer 1

Questo può essere un po 'eccessivo, ma il modo corretto per trovare il tuo punto di intersezione, una volta che avete dividere la curva in blocchi, è quello di vedere se qualsiasi segmento dal primo pezzo si interseca con qualsiasi segmento dal secondo pezzo.

ho intenzione di farmi alcuni dati semplici, un pezzo di un prolate cycloid, e sto andando a trovare luoghi in cui la coordinata y lanci di aumentare a diminuire in modo simile a here:

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

Se si dispone di due segmenti, uno che va dal punto a P1 e un altro che va dal punto Q0 a Q1, è possibile trovare il loro punto di intersezione risolvendo l'equazione vettoriale P0 + s*(P1-P0) = Q0 + t*(Q1-Q0) e i due segmenti si intersecano effettivamente se entrambi i numeri s e t sono in [0, 1].Provare questo per tutti i segmenti:

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

Sul mio sistema, questo in grado di gestire circa 20 incroci al secondo, che non è superveloce, ma probabilmente sufficienti per analizzare i grafici di tanto in tanto. Si può essere in grado di spped le cose da vettorizzazione soluzione dei sistemi 2x2 lineari:

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

Sì, è disordinato, ma funziona, e lo fa volte in modo X100 più veloce:

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 

Answer 2

Si potrebbe fare questo in modo molto semplice utilizzando la funzione interp1d da SciPy per ricampionare i valori y di tutte le linee agli stessi valori x.

http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html