354 lines
11 KiB
Python
354 lines
11 KiB
Python
from lattices import SCC_Lattice, VO2_Lattice
|
|
from spin_image import SpinImage, SpinImage_Two_Phase
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
import matplotlib.patches as patches
|
|
import matplotlib
|
|
import tqdm
|
|
from extractors import Rect_Evaluator, Voronoi_Evaluator, Image_Wrapper
|
|
from cache import timeit
|
|
from scipy import signal
|
|
|
|
|
|
class Plotter:
|
|
def __init__(self, lat):
|
|
self.lattice = lat
|
|
self.length_2 = 0.05
|
|
|
|
def reduce(self, arr):
|
|
arr = np.array(arr)
|
|
arr = arr.flatten()
|
|
return np.mean(arr)
|
|
# return np.sum(arr[np.argpartition(arr, -8)[-8:]])
|
|
|
|
def rect_at_point(self, x, y, color):
|
|
length_2 = self.length_2
|
|
rect = patches.Rectangle(
|
|
(x - length_2, y - length_2),
|
|
2 * length_2,
|
|
2 * length_2,
|
|
linewidth=1,
|
|
edgecolor=color,
|
|
facecolor="none",
|
|
)
|
|
return rect
|
|
|
|
def plot(self, freqx, freqy, intens, ax_log=None, ax_lin=None, vmax=None):
|
|
if ax_log:
|
|
t = ax_log.imshow(
|
|
intens,
|
|
extent=(np.min(freqx), np.max(freqx),
|
|
np.min(freqy), np.max(freqy)),
|
|
norm=matplotlib.colors.LogNorm(vmin=1e-12),
|
|
cmap="viridis",
|
|
origin="lower"
|
|
)
|
|
plt.colorbar(t, ax=ax_log)
|
|
if ax_lin:
|
|
t = ax_lin.imshow(
|
|
intens,
|
|
extent=(np.min(freqx), np.max(freqx),
|
|
np.min(freqy), np.max(freqy)),
|
|
# vmax=vmax,
|
|
cmap="viridis",
|
|
origin="lower"
|
|
)
|
|
plt.colorbar(t, ax=ax_lin)
|
|
|
|
|
|
def plot_periodic(array):
|
|
plt.figure()
|
|
kernel = np.ones((10, 10))
|
|
#array = signal.convolve2d(array, kernel, boundary='symm', mode='same')
|
|
tiled = np.tile(array, (2, 2))
|
|
plt.imshow(tiled)
|
|
|
|
|
|
def rotate(x, y, angle):
|
|
radian = angle / 180 * 2 * np.pi
|
|
return np.cos(radian) * x - np.sin(radian) * y,\
|
|
np.sin(radian) * x + np.cos(radian) * y
|
|
|
|
|
|
def test_square():
|
|
LEN = 40
|
|
lat = SCC_Lattice(LEN, LEN)
|
|
# lat = VO2_Lattice(LEN, LEN)
|
|
plot = Plotter(lat)
|
|
pos_x, pos_y = lat.get_from_mask(np.zeros((40, 40)))
|
|
[(x, y), (x1, y1)] = lat.get_both()
|
|
# pos_x, pos_y = rotate(pos_x, pos_y, 30)
|
|
# si = SpinImage(pos_x, pos_y)
|
|
si = SpinImage_Two_Phase(x, y, x1, y1)
|
|
si.apply_mask(np.zeros((20, 20)))
|
|
fig, axs = plt.subplots(2, 2)
|
|
si.plot(axs[0, 0])
|
|
# si.gaussian(LEN)
|
|
# si.blur(3)
|
|
# si.pad_it_square(additional_pad=1000)
|
|
si.gaussian(20)
|
|
si.plot(axs[0, 1])
|
|
|
|
plt.pause(0.1)
|
|
|
|
fx, fy, intens = si.fft()
|
|
plot.plot(fx, fy, intens, axs[1, 0], axs[1, 1])
|
|
plt.tight_layout()
|
|
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()
|
|
|
|
|
|
def helper(intens, fx, fy, voro, rect):
|
|
print(np.mean(intens))
|
|
v_idx, v_val = voro.extract(Image_Wrapper(intens, fx, fy))
|
|
r_idx, r_val = rect.extract(Image_Wrapper(intens, fx, fy))
|
|
|
|
for iv, av, ir, ar in zip(v_idx, v_val, r_idx, r_val):
|
|
av = np.array(av)
|
|
ar = np.array(ar)
|
|
mask = np.isin(ir, iv)
|
|
print("Test: ", np.all(np.isclose(
|
|
ar[mask], av)), np.isclose(ar[mask], av))
|
|
print(iv, "\n", ar[mask], "\n", av, "\n\n\n")
|
|
|
|
|
|
def test_mixed():
|
|
shift_x = -(5 + 3*28)
|
|
shift_y = -43
|
|
LEN = 40
|
|
lat = VO2_Lattice(LEN, LEN)
|
|
plot = Plotter(lat)
|
|
|
|
pos_x, pos_y = lat.get_from_mask(np.zeros((LEN, LEN)))
|
|
si_mono = SpinImage(pos_x, pos_y)
|
|
si_mono.pad_it(int((LEN * 5.712) / si_mono.resolution)+shift_x,
|
|
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((LEN, LEN)))
|
|
si_rutile = SpinImage(pos_x, pos_y)
|
|
print(int((LEN * 5.712) / si_rutile.resolution)+shift_x,
|
|
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((LEN, LEN))
|
|
ind = np.arange(mask_misk.size)
|
|
np.random.shuffle(ind)
|
|
mask_misk[np.unravel_index(ind[:800], (LEN, LEN))] = False
|
|
pos_x, pos_y = lat.get_from_mask(mask_misk)
|
|
si_mixed = SpinImage(pos_x, pos_y)
|
|
si_mixed.pad_it(int((LEN * 5.712) / si_mixed.resolution)+shift_x,
|
|
int((LEN * 4.554) / si_mixed.resolution)+shift_y)
|
|
fx, fy, intens_mixed = si_mixed.fft()
|
|
|
|
plot_periodic(si_rutile.img)
|
|
plot_periodic(si_mono.img)
|
|
plot_periodic(si_mixed.img)
|
|
|
|
fig, axs = plt.subplots(3, 3)
|
|
si_rutile.plot(axs[0, 0])
|
|
si_mono.plot(axs[0, 2])
|
|
si_mixed.plot(axs[0, 1])
|
|
plot.plot(freqx=fx, freqy=fy, intens=intens_rutile,
|
|
ax_log=axs[1, 0], ax_lin=axs[2, 0], vmax=10e7)
|
|
plot.plot(freqx=fx, freqy=fy, intens=intens_mono,
|
|
ax_log=axs[1, 2], ax_lin=axs[2, 2], vmax=10e7)
|
|
plot.plot(freqx=fx, freqy=fy, intens=intens_mixed,
|
|
ax_log=axs[1, 1], ax_lin=axs[2, 1], vmax=10e7)
|
|
|
|
index = np.abs(fx).argmin()
|
|
print(index, "/", fx.shape)
|
|
fig, axs = plt.subplots(1, 3)
|
|
axs[0].plot(intens_rutile[index, :])
|
|
axs[0].plot(intens_rutile[index-1, :])
|
|
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")
|
|
|
|
index = np.abs(fy).argmin()
|
|
print(index, "/", fy.shape)
|
|
fig, axs = plt.subplots(1, 3)
|
|
axs[0].plot(intens_rutile[:, index])
|
|
axs[1].plot(intens_mixed[:, index])
|
|
axs[2].plot(intens_mono[:, index])
|
|
for a in axs:
|
|
a.set_yscale("log")
|
|
|
|
|
|
def random(seed):
|
|
np.random.seed(seed)
|
|
LEN = 40
|
|
lat = VO2_Lattice(LEN, LEN)
|
|
maske = np.zeros((LEN, LEN))
|
|
ind = np.arange(LEN * LEN)
|
|
np.random.shuffle(ind)
|
|
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])
|
|
|
|
out_rect = [[] for x in range(len(lat.reci())+1)]
|
|
out_voro = [[] for x in range(len(lat.reci())+1)]
|
|
percentage = []
|
|
counter = 0
|
|
already_inited = False
|
|
|
|
for i in tqdm.tqdm(ind):
|
|
maske[np.unravel_index(i, (LEN, LEN))] = True
|
|
counter += 1
|
|
if np.mod(counter, 100) != 0:
|
|
continue
|
|
|
|
pos_x, pos_y = lat.get_from_mask(maske)
|
|
si = SpinImage(pos_x, pos_y)
|
|
si.pad_it_square(10, size=2300)
|
|
fx, fy, intens = si.fft()
|
|
img = Image_Wrapper(intens, fx, fy)
|
|
if not already_inited:
|
|
print("start_init")
|
|
voro.generate_mask(img, merge=True)
|
|
print("stop_init")
|
|
rect.generate_mask(img, merge=True)
|
|
already_inited = True
|
|
|
|
iv, vv = voro.extract(img)
|
|
ir, vr = rect.extract(img)
|
|
for lis, val in zip(out_rect, vr):
|
|
lis.append(val)
|
|
for lis, val in zip(out_voro, vv):
|
|
lis.append(val)
|
|
percentage.append(np.sum(maske))
|
|
|
|
percentage = np.array(percentage)
|
|
percentage /= np.max(percentage)
|
|
|
|
np.savez(f"random_rect_{seed}.npz",
|
|
percentage=percentage, out_1=out_rect[0],
|
|
out_2=out_rect[1], out_3=out_rect[2], out_4=out_rect[3])
|
|
np.savez(f"random_voro_{seed}.npz",
|
|
percentage=percentage, out_1=out_voro[0],
|
|
out_2=out_voro[1], out_3=out_voro[2], out_4=out_voro[3])
|
|
|
|
|
|
def sample_index(p):
|
|
i = np.random.choice(np.arange(p.size), p=p.ravel())
|
|
return np.unravel_index(i, p.shape)
|
|
|
|
|
|
@timeit
|
|
def ising(seed):
|
|
np.random.seed(seed)
|
|
LEN = 40
|
|
temp = 0.1
|
|
maske = np.zeros((LEN, LEN), dtype=bool)
|
|
|
|
lat = VO2_Lattice(LEN, LEN)
|
|
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])
|
|
|
|
out_rect = [[] for x in range(len(lat.reci())+1)]
|
|
out_voro = [[] for x in range(len(lat.reci())+1)]
|
|
percentage = []
|
|
counter = 0
|
|
already_inited = False
|
|
for i in tqdm.tqdm(range(LEN*LEN)):
|
|
probability = np.roll(maske, 1, axis=0).astype(float)
|
|
probability += np.roll(maske, -1, axis=0).astype(float)
|
|
probability += np.roll(maske, 1, axis=1).astype(float)
|
|
probability += np.roll(maske, -1, axis=1).astype(float)
|
|
|
|
probability = np.exp(probability/temp)
|
|
probability[maske] = 0
|
|
probability /= np.sum(probability)
|
|
maske[sample_index(probability)] = True
|
|
|
|
counter += 1
|
|
if np.mod(counter, 100) != 0:
|
|
continue
|
|
|
|
pos_x, pos_y = lat.get_from_mask(maske)
|
|
si = SpinImage(pos_x, pos_y)
|
|
si.pad_it_square(10, size=2300)
|
|
fx, fy, intens = si.fft()
|
|
img = Image_Wrapper(intens, fx, fy)
|
|
if not already_inited:
|
|
voro.generate_mask(img, merge=True)
|
|
rect.generate_mask(img, merge=True)
|
|
already_inited = True
|
|
|
|
iv, vv = voro.extract(img)
|
|
ir, vr = rect.extract(img)
|
|
for lis, val in zip(out_rect, vr):
|
|
lis.append(val)
|
|
for lis, val in zip(out_voro, vv):
|
|
lis.append(val)
|
|
percentage.append(np.sum(maske))
|
|
|
|
percentage = np.array(percentage, dtype=np.float64)
|
|
percentage /= np.max(percentage)
|
|
|
|
np.savez(f"ising_rect_{seed}.npz",
|
|
percentage=percentage, out_1=out_rect[0],
|
|
out_2=out_rect[1], out_3=out_rect[2], out_4=out_rect[3])
|
|
np.savez(f"ising_voro_{seed}.npz",
|
|
percentage=percentage, out_1=out_voro[0],
|
|
out_2=out_voro[1], out_3=out_voro[2], out_4=out_voro[3])
|
|
|
|
|
|
def test_me():
|
|
LEN = 20
|
|
lat = VO2_Lattice(LEN, LEN)
|
|
maske = np.zeros((LEN, LEN), dtype=bool)
|
|
pos_x, pos_y = lat.get_from_mask(maske)
|
|
si = SpinImage(pos_x, pos_y)
|
|
|
|
fig, axs = plt.subplots(1, 3)
|
|
kernel = np.ones((20, 20))
|
|
array = signal.convolve2d(si.img, kernel, boundary='symm', mode='same')
|
|
axs[0].imshow(array)
|
|
|
|
[(x, y), (x1, y1)] = lat.get_both()
|
|
si = SpinImage_Two_Phase(x, y, x1, y1)
|
|
si.apply_mask(maske)
|
|
|
|
array = signal.convolve2d(si.img, kernel, boundary='symm', mode='same')
|
|
# axs[1].imshow(array)
|
|
axs[1].imshow(si.img)
|
|
|
|
np.invert(maske)
|
|
si.apply_mask(maske)
|
|
array = signal.convolve2d(si.img, kernel, boundary='symm', mode='same')
|
|
# axs[1].imshow(array)
|
|
axs[2].imshow(si.img)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
# test_me()
|
|
test_square()
|
|
# test_mixed()
|
|
plt.show()
|
|
# random()
|
|
# np.random.seed(1234)
|
|
# for i in np.random.randint(0, 10000, 1):
|
|
# random(i)
|
|
# ising(i)
|
|
|
|
# plt.show()
|