<jupyter_start><jupyter_text>ZigZag<jupyter_code>das_template_x = np.array([2.5*np.sqrt(2)*i for i in range(24)])
das_template_y = np.array([2.5*np.sqrt(2)*i for i in range(24)])
das_template_x2 = np.hstack([das_template_x,das_template_x[::-1],das_template_x,das_template_x[::-1]])
das_template_y2 = np.hstack([das_template_y,das_template_y+das_template_y[-1],das_template_y+2*das_template_y[-1],das_template_y+3*das_template_y[-1]])

das_x = np.hstack([das_template_x2+i*das_template_x[-1] for i in range(4)])
das_y = np.hstack([das_template_y2 for i in range(4)])

offset = (320-np.max(das_x))/2

das_x += offset
das_y += offset

azimuth_template_1 = np.array([[[45 for i in range(24)], [-45 for i in range(24)]] for i in range(2)]).flatten()
azimuth_template_2 = np.array([[[135 for i in range(24)], [215 for i in range(24)]] for i in range(2)]).flatten()
das_az = np.hstack([azimuth_template_1, azimuth_template_2, 
                    azimuth_template_1, azimuth_template_2])

raz = np.deg2rad(das_az)


cscale = 2

generate_kernels = True


L = 10 #gauge length
ll = np.linspace(-L/2, L/2, 2**5+1)
dl = ll[1]-ll[0]
p1 = das_x[:,np.newaxis]+np.sin(raz[:,np.newaxis])*ll[np.newaxis,:]
p2 = das_y[:,np.newaxis]+np.cos(raz[:,np.newaxis])*ll[np.newaxis,:]


if generate_kernels:
    os.makedirs("Kernels", exist_ok=True)
    crv = loadmat(f"../Curvelet_Basis_Construction/G_{nx}_{ny}.mat")
    G_mat = np.reshape(crv["G_mat"].T, (crv["G_mat"].shape[1], nx, ny))
    crvscales = crv["scales"].flatten()
    cvtscaler = 2.0**(cscale*crvscales)
    G1 = np.zeros((len(raz), G_mat.shape[0]))
    G2 = np.zeros((len(raz), G_mat.shape[0]))
    G3 = np.zeros((len(xs), G_mat.shape[0]))             
    for j in range(G_mat.shape[0]):
        frame = rbs(x,y,G_mat[j])
        #average derivatives of frame along gauge length
        fd1 = romb(frame.ev(p1, p2, dx=1), dl) / L
        fd2 = romb(frame.ev(p1, p2, dy=1), dl) / L
        G1[:,j] = (np.sin(raz)**2*fd1 + 
                  np.sin(2*raz)*fd2/2) / cvtscaler[j]
        G2[:,j] = (np.cos(raz)**2*fd2 + 
                  np.sin(2*raz)*fd1/2) / cvtscaler[j]
        G3[:,j] = frame.ev(xs, ys) / cvtscaler[j]

        
    G = np.hstack([G1, G2])
    Gn = np.max(np.sqrt(np.sum(G**2, axis=1)))
    G = G / Gn
#     Gn=1
    G_zigzag = G
np.linalg.slogdet(G.T@G+1e-10*np.eye(G.shape[1]))
plt.plot(np.sort(np.diag(G @ np.linalg.solve(G.T@G + 1e-10*np.eye(G.shape[1]), G.T))))
exxr = np.array([romb(s.ev(p1, p2, dx=2), dl)/L for s in shotf_itps])
eyyr = np.array([romb(s.ev(p1, p2, dy=2), dl)/L for s in shotf_itps])
exyr = np.array([romb(s.ev(p1, p2, dx=1, dy=1), dl)/L  for s in shotf_itps])
edasr = (np.sin(raz)**2*exxr+np.sin(2*raz)*exyr+np.cos(raz)**2*eyyr) 
das_wvt_data = np.array([np.hstack(pywt.wavedec(d, wvt)) for d in edasr.T])

cuxr = np.array([s.ev(xs, ys, dx=1) for s in shotf_itps])
cuyr = np.array([s.ev(xs, ys, dy=1) for s in shotf_itps])



np.save("Testing/zigzag.npy", das_wvt_data)
wvt_tmp = pywt.wavedec(edasr.T[0], wvt)
wvt_lens = [len(wc) for wc in wvt_tmp]
resi = np.load(f"Testing/zigzag_res.npy")

Gs = np.std(G)

resxi = resi[:G3.shape[1], :]
resyi = resi[G3.shape[1]:, :]

xpredi = (G3/Gn/Gs) @ resxi
ypredi = (G3/Gn/Gs) @ resyi


txpredi = np.real(np.array([reconstruction(w, wvt_lens, wvt) for w in xpredi]))
typredi = np.real(np.array([reconstruction(w, wvt_lens, wvt) for w in ypredi]))


res = np.sqrt(np.mean(np.hstack([np.square(txpredi-cuxr.T), np.square(typredi-cuyr.T)])))/np.std(np.hstack([cuxr, cuyr]))

plt.plot(np.std(resi, axis=1))
cax = plt.scatter(das_x, das_y,color='k', alpha=0.25)
plt.xlim(0,320)
plt.ylim(0,320)
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")

plt.scatter(xs, ys, c=  np.sqrt(np.mean(np.hstack([np.square(txpredi-cuxr.T)/np.std(cuxr, axis=0)[:,np.newaxis]**2, np.square(typredi-cuyr.T)/np.std(cuyr, axis=0)[:,np.newaxis]**2]), axis=1))
, norm=norm)
plt.colorbar()
res
plt.plot(cuxr.T[100])
plt.plot(txpredi[100])
cax = plt.scatter(das_x, das_y,c=das_az)
plt.xlim(0,320)
plt.ylim(0,320)
plt.colorbar(cax, label="Cable Azimuth")
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")

np.sqrt(np.square(das_x[1:]-das_x[:-1])+np.square(das_y[1:]-das_y[:-1]))<jupyter_output><empty_output><jupyter_text>Spiral<jupyter_code>das_theta2 = np.linspace(0,(360*4)**2, 192*2)
das_theta = np.deg2rad(np.sqrt(das_theta2))
a = 0
b = 1
das_r = b*das_theta

das_x = das_r * np.cos(das_theta)
das_y = das_r * np.sin(das_theta)
raz = np.pi/2-np.arctan2(b*np.tan(das_theta)+(a+b*das_theta), b-(a+b*das_theta)*np.tan(das_theta))
das_az = np.rad2deg(raz)

xwidth = np.max(das_x)-np.min(das_x)
das_x = das_x / xwidth * 320
das_y = das_y / xwidth * 320
x_offset = 320 - np.max(das_x)
das_x = das_x + x_offset
y_offset = np.min(das_y)
das_y = das_y - y_offset

L = 10 #gauge length
ll = np.linspace(-L/2, L/2, 2**5+1)
dl = ll[1]-ll[0]
p1 = das_x[:,np.newaxis]+np.sin(raz[:,np.newaxis])*ll[np.newaxis,:]
p2 = das_y[:,np.newaxis]+np.cos(raz[:,np.newaxis])*ll[np.newaxis,:]


if generate_kernels:
    os.makedirs("Kernels", exist_ok=True)
    crv = loadmat(f"../Curvelet_Basis_Construction/G_{nx}_{ny}.mat")
    G_mat = np.reshape(crv["G_mat"].T, (crv["G_mat"].shape[1], nx, ny))
    crvscales = crv["scales"].flatten()
    cvtscaler = 2.0**(cscale*crvscales)
    G1 = np.zeros((len(raz), G_mat.shape[0]))
    G2 = np.zeros((len(raz), G_mat.shape[0]))
    for j in range(G_mat.shape[0]):
        frame = rbs(x,y,G_mat[j])
        #average derivatives of frame along gauge length
        fd1 = romb(frame.ev(p1, p2, dx=1), dl) / L
        fd2 = romb(frame.ev(p1, p2, dy=1), dl) / L
        G1[:,j] = (np.sin(raz)**2*fd1 + 
                  np.sin(2*raz)*fd2/2) / cvtscaler[j]
        G2[:,j] = (np.cos(raz)**2*fd2 + 
                  np.sin(2*raz)*fd1/2) / cvtscaler[j]
        
    G = np.hstack([G1, G2])
    Gn = np.max(np.sqrt(np.sum(G**2, axis=1)))
    G = G / Gn
    G_spiral = G
np.linalg.slogdet(G.T@G+1e-10*np.eye(G.shape[1]))
exxr = np.array([romb(s.ev(p1, p2, dx=2), dl)/L for s in shotf_itps])
eyyr = np.array([romb(s.ev(p1, p2, dy=2), dl)/L for s in shotf_itps])
exyr = np.array([romb(s.ev(p1, p2, dx=1, dy=1), dl)/L  for s in shotf_itps])
edasr = (np.sin(raz)**2*exxr+np.sin(2*raz)*exyr+np.cos(raz)**2*eyyr) 
das_wvt_data = np.array([np.hstack(pywt.wavedec(d, wvt)) for d in edasr.T])

np.save("Testing/spiral.npy", das_wvt_data)
resi = np.load(f"Testing/spiral_res.npy")

Gs = np.std(G)

resxi = resi[:G3.shape[1], :]
resyi = resi[G3.shape[1]:, :]

xpredi = (G3/Gn/Gs) @ resxi
ypredi = (G3/Gn/Gs) @ resyi


txpredi = np.real(np.array([reconstruction(w, wvt_lens, wvt) for w in xpredi]))
typredi = np.real(np.array([reconstruction(w, wvt_lens, wvt) for w in ypredi]))


res = np.sqrt(np.mean(np.hstack([np.square(txpredi-cuxr.T), np.square(typredi-cuyr.T)])))/np.std(np.hstack([cuxr, cuyr]))

res
plt.plot(cuxr.T[100])
plt.plot(txpredi[100])
cax = plt.scatter(das_x, das_y,color='k', alpha=0.25)
plt.xlim(0,320)
plt.ylim(0,320)
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")

plt.scatter(xs, ys, c=  np.sqrt(np.mean(np.hstack([np.square(txpredi-cuxr.T)/np.std(cuxr, axis=0)[:,np.newaxis]**2, np.square(typredi-cuyr.T)/np.std(cuyr, axis=0)[:,np.newaxis]**2]), axis=1))
, norm=norm)
plt.colorbar()
np.hstack([cuxr, cuyr]).shape
cax = plt.scatter(das_x, das_y,c=das_az)
plt.xlim(0,320)
plt.ylim(0,320)
plt.colorbar(cax, label="Cable Azimuth")
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")

np.sqrt(np.square(das_x[1:]-das_x[:-1])+np.square(das_y[1:]-das_y[:-1]))<jupyter_output><empty_output><jupyter_text>Crossing<jupyter_code>template = np.linspace(0,320, 65)
template = (template[1:]+template[:-1])/2

das_x = np.hstack([template, template, template,[80 for i in range(len(template))], [160 for i in range(len(template))], [240 for i in range(len(template))]])
das_y = np.hstack([[80 for i in range(len(template))], [160 for i in range(len(template))], [240 for i in range(len(template))],template,template,template])
das_az = np.hstack([[90 for i in range(len(template))], [270 for i in range(len(template))], [90 for i in range(len(template))],[0 for i in range(len(template))], [180 for i in range(len(template))], [0 for i in range(len(template))]])
raz = np.deg2rad(das_az)

L = 10 #gauge length
ll = np.linspace(-L/2, L/2, 2**5+1)
dl = ll[1]-ll[0]
p1 = das_x[:,np.newaxis]+np.sin(raz[:,np.newaxis])*ll[np.newaxis,:]
p2 = das_y[:,np.newaxis]+np.cos(raz[:,np.newaxis])*ll[np.newaxis,:]


if generate_kernels:
    os.makedirs("Kernels", exist_ok=True)
    crv = loadmat(f"../Curvelet_Basis_Construction/G_{nx}_{ny}.mat")
    G_mat = np.reshape(crv["G_mat"].T, (crv["G_mat"].shape[1], nx, ny))
    crvscales = crv["scales"].flatten()
    cvtscaler = 2.0**(cscale*crvscales)
    G1 = np.zeros((len(raz), G_mat.shape[0]))
    G2 = np.zeros((len(raz), G_mat.shape[0]))
    for j in range(G_mat.shape[0]):
        frame = rbs(x,y,G_mat[j])
        #average derivatives of frame along gauge length
        fd1 = romb(frame.ev(p1, p2, dx=1), dl) / L
        fd2 = romb(frame.ev(p1, p2, dy=1), dl) / L
        G1[:,j] = (np.sin(raz)**2*fd1 + 
                  np.sin(2*raz)*fd2/2) / cvtscaler[j]
        G2[:,j] = (np.cos(raz)**2*fd2 + 
                  np.sin(2*raz)*fd1/2) / cvtscaler[j]
        
    G = np.hstack([G1, G2])
    Gn = np.max(np.sqrt(np.sum(G**2, axis=1)))
    G = G / Gn
    G_cross = G
np.linalg.slogdet(G.T@G+1e-10*np.eye(G.shape[1]))
exxr = np.array([romb(s.ev(p1, p2, dx=2), dl)/L for s in shotf_itps])
eyyr = np.array([romb(s.ev(p1, p2, dy=2), dl)/L for s in shotf_itps])
exyr = np.array([romb(s.ev(p1, p2, dx=1, dy=1), dl)/L  for s in shotf_itps])
edasr = (np.sin(raz)**2*exxr+np.sin(2*raz)*exyr+np.cos(raz)**2*eyyr) 
das_wvt_data = np.array([np.hstack(pywt.wavedec(d, wvt)) for d in edasr.T])

np.save("Testing/crossing.npy", das_wvt_data )
resi = np.load(f"Testing/crossing_res.npy")

Gs = np.std(G)

resxi = resi[:G3.shape[1], :]
resyi = resi[G3.shape[1]:, :]

xpredi = (G3/Gn/Gs) @ resxi
ypredi = (G3/Gn/Gs) @ resyi


txpredi = np.real(np.array([reconstruction(w, wvt_lens, wvt) for w in xpredi]))
typredi = np.real(np.array([reconstruction(w, wvt_lens, wvt) for w in ypredi]))


res = np.sqrt(np.mean(np.hstack([np.square(txpredi-cuxr.T), np.square(typredi-cuyr.T)])))/np.std(np.hstack([cuxr, cuyr]))

res
plt.plot(cuxr.T[150])
plt.plot(txpredi[150])
cax = plt.scatter(das_x, das_y,color='k', alpha=0.25)
plt.xlim(0,320)
plt.ylim(0,320)
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")

plt.scatter(xs, ys, c=  np.sqrt(np.mean(np.hstack([np.square(txpredi-cuxr.T)/np.std(cuxr, axis=0)[:,np.newaxis]**2, np.square(typredi-cuyr.T)/np.std(cuyr, axis=0)[:,np.newaxis]**2]), axis=1))
, norm=norm)
plt.scatter(xs[150], ys[150], color='r')
plt.colorbar()
cax = plt.scatter(das_x, das_y,c=das_az)
plt.xlim(0,320)
plt.ylim(0,320)
plt.colorbar(cax, label="Cable Azimuth")
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")<jupyter_output><empty_output><jupyter_text>Random<jupyter_code>template = np.linspace(0,320, 65)
template = (template[1:]+template[:-1])/2

np.random.seed(94899109)
das_x = np.random.uniform(5,315,384)
das_y = np.random.uniform(5,315,384)
das_az = np.random.uniform(0,360,384)

raz = np.deg2rad(das_az)

L = 10 #gauge length
ll = np.linspace(-L/2, L/2, 2**5+1)
dl = ll[1]-ll[0]
p1 = das_x[:,np.newaxis]+np.sin(raz[:,np.newaxis])*ll[np.newaxis,:]
p2 = das_y[:,np.newaxis]+np.cos(raz[:,np.newaxis])*ll[np.newaxis,:]


if generate_kernels:
    os.makedirs("Kernels", exist_ok=True)
    crv = loadmat(f"../Curvelet_Basis_Construction/G_{nx}_{ny}.mat")
    G_mat = np.reshape(crv["G_mat"].T, (crv["G_mat"].shape[1], nx, ny))
    crvscales = crv["scales"].flatten()
    cvtscaler = 2.0**(cscale*crvscales)
    G1 = np.zeros((len(raz), G_mat.shape[0]))
    G2 = np.zeros((len(raz), G_mat.shape[0]))
    for j in range(G_mat.shape[0]):
        frame = rbs(x,y,G_mat[j])
        #average derivatives of frame along gauge length
        fd1 = romb(frame.ev(p1, p2, dx=1), dl) / L
        fd2 = romb(frame.ev(p1, p2, dy=1), dl) / L
        G1[:,j] = (np.sin(raz)**2*fd1 + 
                  np.sin(2*raz)*fd2/2) / cvtscaler[j]
        G2[:,j] = (np.cos(raz)**2*fd2 + 
                  np.sin(2*raz)*fd1/2) / cvtscaler[j]
        
    G = np.hstack([G1, G2])
    Gn = np.max(np.sqrt(np.sum(G**2, axis=1)))
    G = G / Gn
    G_random = G
np.linalg.slogdet(G.T@G+1e-10*np.eye(G.shape[1]))
exxr = np.array([romb(s.ev(p1, p2, dx=2), dl)/L for s in shotf_itps])
eyyr = np.array([romb(s.ev(p1, p2, dy=2), dl)/L for s in shotf_itps])
exyr = np.array([romb(s.ev(p1, p2, dx=1, dy=1), dl)/L  for s in shotf_itps])
edasr = (np.sin(raz)**2*exxr+np.sin(2*raz)*exyr+np.cos(raz)**2*eyyr) 
das_wvt_data = np.array([np.hstack(pywt.wavedec(d, wvt)) for d in edasr.T])

np.save("Testing/random.npy", das_wvt_data)
resi = np.load(f"Testing/random_res.npy")

Gs = np.std(G)

resxi = resi[:G3.shape[1], :]
resyi = resi[G3.shape[1]:, :]

xpredi = (G3/Gn/Gs) @ resxi
ypredi = (G3/Gn/Gs) @ resyi


txpredi = np.real(np.array([reconstruction(w, wvt_lens, wvt) for w in xpredi]))
typredi = np.real(np.array([reconstruction(w, wvt_lens, wvt) for w in ypredi]))


res = np.sqrt(np.mean(np.hstack([np.square(txpredi-cuxr.T), np.square(typredi-cuyr.T)])))/np.std(np.hstack([cuxr, cuyr]))

plt.plot(np.std(resi, axis=1))
res
plt.plot(cuxr.T[100])
plt.plot(txpredi[100])
cax = plt.scatter(das_x, das_y,color='k', alpha=0.25)
plt.xlim(0,320)
plt.ylim(0,320)
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")

plt.scatter(xs, ys, c=  np.sqrt(np.mean(np.hstack([np.square(txpredi-cuxr.T)/np.std(cuxr, axis=0)[:,np.newaxis]**2, np.square(typredi-cuyr.T)/np.std(cuyr, axis=0)[:,np.newaxis]**2]), axis=1))
, norm=norm)
plt.colorbar()
cax = plt.scatter(das_x, das_y,c=das_az)
plt.xlim(0,320)
plt.ylim(0,320)
plt.colorbar(cax, label="Cable Azimuth")
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")
np.save("Kernels/G_zigzag.npy", G_zigzag)
np.save("Kernels/G_spiral.npy", G_spiral)
np.save("Kernels/G_cross.npy", G_cross)
np.save("Kernels/G_random.npy", G_random)<jupyter_output><empty_output><jupyter_text>Eigenvalue Spectrum<jupyter_code>G_full = np.vstack([np.hstack([G3, np.zeros(G3.shape)]), np.hstack([np.zeros(G3.shape), G3])])

idet = 1e-10*np.eye(G_zigzag.shape[1])
ezig = np.sort(np.linalg.eigvals(G_zigzag.T@G_zigzag+idet))[::-1]
espi = np.sort(np.linalg.eigvals(G_spiral.T@G_spiral+idet))[::-1]
ecro = np.sort(np.linalg.eigvals(G_cross.T@G_cross+idet))[::-1]
eran = np.sort(np.linalg.eigvals(G_random.T@G_random+idet))[::-1]
# efull= np.sort(np.linalg.eigvals(G_full.T@G_full+idet))[::-1]

ezign = ezig / ezig[0]
espin = espi / espi[0]
ecron = ecro / ecro[0]
erann = eran / eran[0]
# efulln = efull / efull[0]

plt.plot(np.sort(np.diag(G_zigzag @ np.linalg.solve(G_zigzag.T@G_zigzag + 1e-0*np.eye(G_zigzag.shape[1]), G_zigzag.T))), label="ZigZag")
plt.plot(np.sort(np.diag(G_spiral @ np.linalg.solve(G_spiral.T@G_spiral + 1e-0*np.eye(G_spiral.shape[1]), G_spiral.T))), label="Spiral")
plt.plot(np.sort(np.diag(G_cross @ np.linalg.solve(G_cross.T@G_cross + 1e-0*np.eye(G_cross.shape[1]), G_cross.T))), label="Crossing")
plt.plot(np.sort(np.diag(G_random @ np.linalg.solve(G_random.T@G_random + 1e-0*np.eye(G_random.shape[1]), G_random.T))), label="Random")

plt.xlim(0,384)
plt.ylabel("Coherence")
plt.xlabel("Sorted Diagonal")
plt.legend(loc="lower right")
plt.plot(np.log10(np.real(ezig)), label="ZigZag")
plt.plot(np.log10(np.real(espi)), label="Spiral")
plt.plot(np.log10(np.real(ecro)), label="Crossing")
plt.plot(np.log10(np.real(eran)), label="Random")

plt.xlim(0,384)
plt.ylabel("Log10 Normalized Eigenvalues")
plt.xlabel("Eigenvalue Index")
plt.legend(loc="upper right")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.22058545797922569, 'fraction_numerical': 0.03896757947749449, 'mean_word_length': 4.605504587155964, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '2814', 'license': 'MIT', 'n_tokens_mistral': 8005, 'n_tokens_neox': 7632, 'n_words': 1090, 'path': 'OptimalDesign/DAS_Designs.ipynb', 'repo_name': 'jbmuir/DAS-Reconstruction', 'size': '15885'}