jacob/FFT
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

a lot of py stuff

Browse Source
This commit is contained in:
Jacob Holder 2023-03-14 22:16:07 +01:00
parent ed67e6f7f5
commit 7e5eb6e8c7
Signed by: jacob
GPG Key ID: 2194FC747048A7FD
3 changed files with 98 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
2d_fourie/lattices.py
View File

@ -1,6 +1,7 @@
import numpy as np import numpy as np
from cache import timeit from cache import timeit
def deg_2_rad(winkel): def deg_2_rad(winkel):
return winkel / 180.0 * np.pi return winkel / 180.0 * np.pi
@ -42,7 +43,9 @@ class VO2_Lattice(Lattice):
base_c_m = 5.38 base_c_m = 5.38
base_c_r = 2.856 base_c_r = 2.856
base_c_r = 2.8
base_b_r = 4.554 base_b_r = 4.554
base_b_r = 4.5
base_a_r = base_b_r base_a_r = base_b_r
alpha_m = 122.64 # degree alpha_m = 122.64 # degree
@ -64,6 +67,14 @@ class VO2_Lattice(Lattice):
offset_a_r, offset_c_r = self._mono_2_rutile(offset_c_m, offset_a_m) offset_a_r, offset_c_r = self._mono_2_rutile(offset_c_m, offset_a_m)
res = 0.05
offset_a_r = res * int(offset_a_r/res)
offset_c_r = res * int(offset_c_r/res)
offset_a_r = 0.5
offset_c_r = 0.5
print(offset_a_r, offset_c_r)
self.atom_x_mono = offset_a_r + X * \ self.atom_x_mono = offset_a_r + X * \
self.base_c_r + np.mod(Y, 4) * 0.5 * self.base_c_r self.base_c_r + np.mod(Y, 4) * 0.5 * self.base_c_r
self.atom_x_mono[np.mod(X, 2) == 0] -= 2 * offset_a_r self.atom_x_mono[np.mod(X, 2) == 0] -= 2 * offset_a_r

143
2d_fourie/main.py
View File

@ -8,6 +8,7 @@ import tqdm
from extractors import Rect_Evaluator, Voronoi_Evaluator, Image_Wrapper from extractors import Rect_Evaluator, Voronoi_Evaluator, Image_Wrapper
from cache import timeit from cache import timeit
class Plotter: class Plotter:
def __init__(self, lat): def __init__(self, lat):
self.lattice = lat self.lattice = lat
@ -61,7 +62,7 @@ def rotate(x, y, angle):
def test_square(): def test_square():
LEN = 10 LEN = 40
lat = SCC_Lattice(LEN, LEN) lat = SCC_Lattice(LEN, LEN)
# lat = VO2_Lattice(LEN, LEN) # lat = VO2_Lattice(LEN, LEN)
plot = Plotter(lat) plot = Plotter(lat)
@ -69,16 +70,26 @@ def test_square():
# pos_x, pos_y = rotate(pos_x, pos_y, 30) # pos_x, pos_y = rotate(pos_x, pos_y, 30)
si = SpinImage(pos_x, pos_y) si = SpinImage(pos_x, pos_y)
fig, axs = plt.subplots(2, 2) fig, axs = plt.subplots(2, 2)
si.plot(axs[0, 0],1) si.plot(axs[0, 0], 1)
# si.gaussian(LEN) # si.gaussian(LEN)
# si.blur(3) # si.blur(3)
si.pad_it_square(int((5./2.)/0.1)) si.pad_it_square(size=int((5*LEN)/si.resolution))
si.plot(axs[0, 1],1) si.plot(axs[0, 1], 1)
plt.pause(0.1) plt.pause(0.1)
fx, fy, intens = si.fft() fx, fy, intens = si.fft()
plot.plot(fx, fy, intens, axs[1, 0], axs[1, 1]) plot.plot(fx, fy, intens, axs[1, 0], axs[1, 1])
plt.tight_layout()
plt.savefig("test.pdf") plt.savefig("test.pdf")
plt.figure()
index = np.abs(fy).argmin()
plt.plot(fx, intens[index, :], label="fy=0")
plt.plot(fx, intens[index+1, :], label="fy=df")
plt.plot(fx, intens[index+2, :], label="fy=2df")
plt.legend()
plt.yscale("log")
plt.show() plt.show()
@ -97,87 +108,67 @@ def helper(intens, fx, fy, voro, rect):
def test_mixed(): def test_mixed():
shift_x = -(5 + 3*28)
shift_y = -43
LEN = 40 LEN = 40
lat = VO2_Lattice(LEN, LEN) lat = VO2_Lattice(LEN, LEN)
plot = Plotter(lat) plot = Plotter(lat)
all_rutile = np.stack(lat.reci()[0]).T
all_mono = np.stack(lat.reci()[1]).T
all_mono2 = np.stack(lat.reci()[2]).T
voro = Voronoi_Evaluator([all_rutile, all_mono, all_mono2])
rect = Rect_Evaluator([all_rutile, all_mono, all_mono2])
#voro = Voronoi_Evaluator([all_mono])
#rect = Rect_Evaluator([all_mono])
pos_x, pos_y = lat.get_from_mask(np.zeros((40, 40))) pos_x, pos_y = lat.get_from_mask(np.zeros((LEN, LEN)))
si = SpinImage(pos_x, pos_y) si_mono = SpinImage(pos_x, pos_y)
si.pad_it_square(10, size=2300) si_mono.pad_it(int((LEN * 5.712) / si_mono.resolution)+shift_x,
fx, fy, intens_mono = si.fft() int((LEN * 4.554) / si_mono.resolution)+shift_y)
fx, fy, intens_mono = si_mono.fft()
pos_x, pos_y = lat.get_from_mask(np.ones((40, 40))) pos_x, pos_y = lat.get_from_mask(np.ones((LEN, LEN)))
si = SpinImage(pos_x, pos_y) si_rutile = SpinImage(pos_x, pos_y)
si.pad_it_square(10, size=2300) print(int((LEN * 5.712) / si_rutile.resolution)+shift_x,
fx, fy, intens_rutile = si.fft() int((LEN * 4.554) / si_rutile.resolution)+shift_y,
si_rutile.img.shape)
si_rutile.pad_it(int((LEN * 5.712) / si_rutile.resolution)+shift_x,
int((LEN * 4.554) / si_rutile.resolution)+shift_y)
fx, fy, intens_rutile = si_rutile.fft()
mask_misk = np.ones((40, 40)) mask_misk = np.ones((LEN, LEN))
ind = np.arange(mask_misk.size) ind = np.arange(mask_misk.size)
np.random.shuffle(ind) np.random.shuffle(ind)
mask_misk[np.unravel_index(ind[:800], (40, 40))] = False mask_misk[np.unravel_index(ind[:800], (LEN, LEN))] = False
pos_x, pos_y = lat.get_from_mask(mask_misk) pos_x, pos_y = lat.get_from_mask(mask_misk)
si = SpinImage(pos_x, pos_y) si_mixed = SpinImage(pos_x, pos_y)
si.pad_it_square(10, size=2300) si_mixed.pad_it(int((LEN * 5.712) / si_mixed.resolution)+shift_x,
fx, fy, intens_mixed = si.fft() int((LEN * 4.554) / si_mixed.resolution)+shift_y)
img = Image_Wrapper(intens_mono, fx, fy) fx, fy, intens_mixed = si_mixed.fft()
voro.generate_mask(img)
rect.generate_mask(Image_Wrapper(intens_mono, fx, fy))
fig, axs = plt.subplots(1, 2) fig, axs = plt.subplots(3, 3)
axs[0].imshow(rect.get_mask(), extent=img.ext(), origin="lower") si_rutile.plot(axs[0, 0], 1)
axs[0].plot(all_rutile[:, 0], all_rutile[:, 1], ".") si_mono.plot(axs[0, 2], 1)
axs[1].plot(all_rutile[:, 0], all_rutile[:, 1], ".") si_mixed.plot(axs[0, 1], 1)
axs[0].plot(all_mono[:, 0], all_mono[:, 1], ".") plot.plot(freqx=fx, freqy=fy, intens=intens_rutile,
axs[1].plot(all_mono[:, 0], all_mono[:, 1], ".") ax_log=axs[1, 0], ax_lin=axs[2, 0], vmax=10e7)
axs[0].plot(all_mono2[:, 0], all_mono2[:, 1], ".") plot.plot(freqx=fx, freqy=fy, intens=intens_mono,
axs[1].plot(all_mono2[:, 0], all_mono2[:, 1], ".") ax_log=axs[1, 2], ax_lin=axs[2, 2], vmax=10e7)
axs[1].imshow(voro.get_mask(), extent=img.ext(), origin="lower") plot.plot(freqx=fx, freqy=fy, intens=intens_mixed,
ax_log=axs[1, 1], ax_lin=axs[2, 1], vmax=10e7)
print("mono") index = np.abs(fx).argmin()
helper(intens=intens_mono, fx=fx, fy=fy, voro=voro, rect=rect) print(index, "/", fx.shape)
print("mixed") fig, axs = plt.subplots(1, 3)
helper(intens=intens_mixed, fx=fx, fy=fy, voro=voro, rect=rect) axs[0].plot(intens_rutile[index, :])
print("rutile") axs[0].plot(intens_rutile[index-1, :])
helper(intens=intens_rutile, fx=fx, fy=fy, voro=voro, rect=rect) axs[0].plot(intens_rutile[index+1, :])
axs[1].plot(intens_mixed[index, :])
axs[2].plot(intens_mono[index, :])
for a in axs:
a.set_yscale("log")
new_intens_mono = rect.debug(Image_Wrapper(intens_mono, fx, fy)) index = np.abs(fy).argmin()
new_intens_mixed = rect.debug(Image_Wrapper(intens_mixed, fx, fy)) print(index, "/", fy.shape)
new_intens_rutile = rect.debug(Image_Wrapper(intens_rutile, fx, fy)) fig, axs = plt.subplots(1, 3)
axs[0].plot(intens_rutile[:, index])
fig, axs = plt.subplots(2, 3) axs[1].plot(intens_mixed[:, index])
plot.plot(freqx=fx, freqy=fy, intens=new_intens_rutile, axs[2].plot(intens_mono[:, index])
ax_log=axs[0, 0], ax_lin=axs[1, 0], vmax=10e7) for a in axs:
plot.plot(freqx=fx, freqy=fy, intens=new_intens_mono, a.set_yscale("log")
ax_log=axs[0, 2], ax_lin=axs[1, 2], vmax=10e7)
plot.plot(freqx=fx, freqy=fy, intens=new_intens_mixed,
ax_log=axs[0, 1], ax_lin=axs[1, 1], vmax=10e7)
for ax in axs.flatten():
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
new_intens_mono = voro.debug(Image_Wrapper(new_intens_mono, fx, fy))
new_intens_mixed = voro.debug(Image_Wrapper(intens_mixed, fx, fy))
new_intens_rutile = voro.debug(Image_Wrapper(intens_rutile, fx, fy))
fig, axs = plt.subplots(2, 3)
plot.plot(freqx=fx, freqy=fy, intens=new_intens_rutile,
ax_log=axs[0, 0], ax_lin=axs[1, 0], vmax=10e7)
plot.plot(freqx=fx, freqy=fy, intens=new_intens_mono,
ax_log=axs[0, 2], ax_lin=axs[1, 2], vmax=10e7)
plot.plot(freqx=fx, freqy=fy, intens=new_intens_mixed,
ax_log=axs[0, 1], ax_lin=axs[1, 1], vmax=10e7)
for ax in axs.flatten():
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
def random(seed): def random(seed):
@ -199,7 +190,6 @@ def random(seed):
counter = 0 counter = 0
already_inited = False already_inited = False
for i in tqdm.tqdm(ind): for i in tqdm.tqdm(ind):
maske[np.unravel_index(i, (LEN, LEN))] = True maske[np.unravel_index(i, (LEN, LEN))] = True
counter += 1 counter += 1
@ -241,6 +231,7 @@ def sample_index(p):
i = np.random.choice(np.arange(p.size), p=p.ravel()) i = np.random.choice(np.arange(p.size), p=p.ravel())
return np.unravel_index(i, p.shape) return np.unravel_index(i, p.shape)
@timeit @timeit
def ising(seed): def ising(seed):
np.random.seed(seed) np.random.seed(seed)
@ -305,8 +296,8 @@ def ising(seed):
if __name__ == "__main__": if __name__ == "__main__":
test_square() # test_square()
#test_mixed() test_mixed()
plt.show() plt.show()
# random() # random()
# np.random.seed(1234) # np.random.seed(1234)

22
2d_fourie/spin_image.py
View File

@ -3,8 +3,9 @@ import scipy
import scipy.fftpack as sfft import scipy.fftpack as sfft
import scipy.signal import scipy.signal
class SpinImage: class SpinImage:
resolution = 0.1 resolution = 0.05
def __init__(self, x_pos, y_pos): def __init__(self, x_pos, y_pos):
x_pos = x_pos - np.min(x_pos) x_pos = x_pos - np.min(x_pos)
@ -34,7 +35,7 @@ class SpinImage:
def pad_it_square(self, additional_pad=0, size=None): def pad_it_square(self, additional_pad=0, size=None):
h = self.img.shape[0] h = self.img.shape[0]
w = self.img.shape[1] w = self.img.shape[1]
xx = np.maximum(h, w) + 2 * additional_pad xx = np.maximum(h, w) + 2 * additional_pad + 1
if size is not None: if size is not None:
xx = np.maximum(xx, size) xx = np.maximum(xx, size)
yy = xx yy = xx
@ -50,6 +51,23 @@ class SpinImage:
self.img = np.pad(self.img, pad_width=( self.img = np.pad(self.img, pad_width=(
(a, aa), (b, bb)), mode="constant") (a, aa), (b, bb)), mode="constant")
def pad_it(self, x_size, y_size):
h = self.img.shape[0]
w = self.img.shape[1]
xx = x_size
yy = y_size
self.length_x = xx * self.resolution
self.length_y = yy * self.resolution
a = (xx - h) // 2
aa = xx - a - h
b = (yy - w) // 2
bb = yy - b - w
self.img = np.pad(self.img, pad_width=(
(a, aa), (b, bb)), mode="constant")
def percentage_gaussian(self, mask, sigma): def percentage_gaussian(self, mask, sigma):
x = np.linspace(-self.length_x / 2, x = np.linspace(-self.length_x / 2,
self.length_x / 2, mask.shape[0]) self.length_x / 2, mask.shape[0])