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(): fig, axs = plt.subplots(3, 3) LEN = 80 lat = VO2_Lattice(LEN, LEN) plot = Plotter(lat) [(x, y), (x1, y1)] = lat.get_both() si = SpinImage_Two_Phase(x, y, x1, y1) mask_misk = np.ones((LEN, LEN), dtype=bool) ind = np.arange(mask_misk.size) np.random.shuffle(ind) mask_misk[np.unravel_index(ind[:800], (LEN, LEN))] = False si.apply_mask(np.zeros((LEN, LEN), dtype=bool)) si.gaussian(20) fx, fy, intens_mono = si.fft() si.plot(axs[0, 2]) si.apply_mask(np.ones((LEN, LEN))) si.gaussian(20) fx, fy, intens_rutile = si.fft() si.plot(axs[0, 0]) si.apply_mask(mask_misk) si.gaussian(20) fx, fy, intens_mixed = si.fft() si.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) # Plotting cuts def random(seed): np.random.seed(seed) LEN = 80 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, 4) 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) axs[1].imshow(si.img) tmp = si.img maske = np.invert(maske) si.apply_mask(maske) axs[2].imshow(si.img) axs[3].imshow(si.img+tmp) 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()