From 405c25438eb38446a48d86587268fb86cfe666d7 Mon Sep 17 00:00:00 2001 From: Jacob Date: Thu, 27 Apr 2023 15:40:38 +0200 Subject: [PATCH] removed old stuff --- 2d_fourie/analysis.py | 105 ------------ 2d_fourie/cache.py | 68 -------- 2d_fourie/extractors.py | 213 ------------------------ 2d_fourie/lattices.py | 152 ------------------ 2d_fourie/main.py | 329 -------------------------------------- 2d_fourie/plot.py | 5 - 2d_fourie/spin_image.py | 270 ------------------------------- 2d_fourie/test_voronoi.py | 46 ------ rust/Cargo.toml | 11 -- rust/plot.py | 9 -- rust/src/main.rs | 19 --- rust/src/vo2.rs | 143 ----------------- 12 files changed, 1370 deletions(-) delete mode 100644 2d_fourie/analysis.py delete mode 100644 2d_fourie/cache.py delete mode 100644 2d_fourie/extractors.py delete mode 100644 2d_fourie/lattices.py delete mode 100644 2d_fourie/main.py delete mode 100644 2d_fourie/plot.py delete mode 100644 2d_fourie/spin_image.py delete mode 100644 2d_fourie/test_voronoi.py delete mode 100644 rust/Cargo.toml delete mode 100644 rust/plot.py delete mode 100644 rust/src/main.rs delete mode 100644 rust/src/vo2.rs diff --git a/2d_fourie/analysis.py b/2d_fourie/analysis.py deleted file mode 100644 index 18f071a..0000000 --- a/2d_fourie/analysis.py +++ /dev/null @@ -1,105 +0,0 @@ -import sys -import numpy as np -import matplotlib.pyplot as plt -import glob - - -def eval_data(file): - data = np.load(file) - percentage = data["percentage"] - out = data["out"] - out = np.array(out) - print(out.shape) - fig, all_axs = plt.subplots(2, out.shape[0]) - axs = all_axs[0, :] - axs2 = all_axs[1, :] - for o, ax, ax2, lab in zip(out, axs, axs2, ["rutile", "mono_twin", "mono"]): - # ax.plot(percentage, o/np.max(o, axis=0)) - ax.plot(percentage, o/o[0]) - # ax.plot(percentage, o) - o = np.mean(o, axis=1) - o = o/o[0] - ax2.plot(percentage, o) - ax2.plot([0, 1], [o[0], o[-1]], "k:") - ax.set_title(lab) - - title = "" - if "ising" in file: - title += "Ising " - else: - title += "Random " - if "rect" in file: - title += "rect " - if "voro" in file: - title += "voro " - fig.suptitle(title) - fig.savefig(f"{file}.png") - - -def parse_lists(out): - lists = [] - for o in out: - lists.append(np.stack(o)) - - max = 0 - for l in lists: - print(l.shape) - if max < l.shape[1]: - max = l.shape[1] - lists = [np.pad(l, ((0, 0), (0, max-l.shape[1]))) for l in lists] - for l in lists: - print(l.shape) - return np.stack(lists) - - -def eval_data_print(file): - data = np.load(file, allow_pickle=True) - percentage = data["percentage"] - # out = parse_lists(data["out"]) - out = [] - for o in ["out_1", "out_2", "out_3", "out_4"]: - out.append(np.array(data[o])) - fig, all_axs = plt.subplots(2, len(out)) - axs = all_axs[0, :] - axs2 = all_axs[1, :] - for o, ax, ax2, lab in zip(out, axs, axs2, ["rutile", "monoclinic", "mono", "rest"]): - # ax.plot(percentage, o/np.max(o, axis=0)) - ax.plot(percentage, o/o[0]) - # ax.plot(percentage, o) - o = np.mean(o, axis=1) - # o = o/o[0] - ax2.plot(percentage, o) - ax2.plot([0, 1], [o[0], o[-1]], "k:") - ax.set_title(lab) - for ax in all_axs.flatten(): - ax.set_xlabel("Rutile Phase") - ax.set_ylabel("Normalized Intensity") - - title = "" - if "ising" in file: - title += "Ising " - else: - title += "Random " - if "rect" in file: - title += "rect " - if "voro" in file: - title += "voro " - fig.suptitle(title) - plt.tight_layout() - - -def stacked_plot(file): - data = np.load(file, allow_pickle=True) - percentage = data["percentage"] - # out = parse_lists(data["out"]) - out = [] - for o in ["out_1", "out_2", "out_3", "out_4"]: - out.append(np.array(data[o])) - out = np.array(out) - plt.stackplot(percentage, out) - -if __name__ == "__main__": - for f in sys.argv[1:]: - #eval_data_print(f) - stacked_plot(f) - plt.show() diff --git a/2d_fourie/cache.py b/2d_fourie/cache.py deleted file mode 100644 index be741e3..0000000 --- a/2d_fourie/cache.py +++ /dev/null @@ -1,68 +0,0 @@ -import numpy as np -from functools import wraps -import time - - -def timeit(func): - @wraps(func) - def timeit_wrapper(*args, **kwargs): - start_time = time.perf_counter() - print(f"Start Function {func.__name__}:") - result = func(*args, **kwargs) - end_time = time.perf_counter() - total_time = end_time - start_time - # first item in the args, ie `args[0]` is `self` - print( - f'Function {func.__name__} Took {total_time:.4f} seconds') - return result - return timeit_wrapper - - -def persist_to_file(file_name): - - def decorator(original_func): - try: - file_nam = file_name - if file_name[-4:] != ".npz": - file_nam += ".npz" - file = np.load(file_nam) - cache = dict(zip((file.files), (file[k] for k in file.files))) - except (IOError, ValueError): - cache = {} - - def hash_func(*param): - key = str(param) - return key - - def persist_func(*param, hash=None): - if hash is None: - hash = hash_func(*param) - - print("Hash: ", hash) - if cache == {} or ("hash" not in cache) or hash != cache["hash"]: - print("recalc") - data = original_func(*param) - np.savez(file_name, hash=hash, dat=data) - cache["hash"] = hash - cache["dat"] = data - print("loaded") - return cache["dat"] - - return persist_func - - return decorator - - -# test = 1234 -# -# -# @persist_to_file("cache.npz", test) -# def test(): -# print("calculate") -# return np.zeros((100, 100)) -# -# -# if __name__ == "__main__": -# test() -# test() -# pass diff --git a/2d_fourie/extractors.py b/2d_fourie/extractors.py deleted file mode 100644 index 38e78ee..0000000 --- a/2d_fourie/extractors.py +++ /dev/null @@ -1,213 +0,0 @@ -import numpy as np -from scipy.spatial import Voronoi -import cv2 -from cache import persist_to_file, timeit - - -class Image_Wrapper: - # def __init__(self, img, x_lower, x_res, y_lower, y_res): - # self.img = img - # self.x_lower = x_lower - # self.y_lower = y_lower - # self.x_res = x_res - # self.y_res = y_res - - # self.x_upper = self.x_lower + self.img.shape[0]*self.x_res - self.x_res - # self.y_upper = self.y_lower + self.img.shape[1]*self.y_res - self.x_res - - def __init__(self, img, fx, fy): - self.img = img - self.x_lower = np.min(fx) - self.y_lower = np.min(fy) - self.x_upper = np.max(fx) - self.y_upper = np.max(fy) - - self.x_res = (self.x_upper - self.x_lower) / self.img.shape[0] - self.y_res = (self.y_upper - self.y_lower) / self.img.shape[0] - - def val2pos(self, x, y): - x = (x - self.x_lower) / self.x_res - y = (y - self.y_lower) / self.y_res - return x, y - - def check_bounds(self, xl, yl, xu, yu): - if xl > self.img.shape[0]: - print("xl lim") - return False - if yl > self.img.shape[1]: - print("yl lim") - return False - if xu < 0: - print("xu lim") - return False - if yu < 0: - print("yu lim") - return False - return True - - def clean_bounds(self, xl, yl, xu, yu): - if xl < 0: - xl = 0 - if yl < 0: - yl = 0 - - if xu > self.img.shape[0]: - xu = self.img.shape[0] - if yu > self.img.shape[1]: - yu = self.img.shape[1] - - return xl, yl, xu, yu - - def ext(self): - return [self.x_lower, self.x_upper, self.y_lower, self.y_upper] - - -class Evaluator: - def extract(self, img: Image_Wrapper): - all_val = [] - all_idx = [] - for ev_points in self.eval_points: - temp_val = [] - temp_idx = [] - for num in ev_points: - if np.sum(self.mask == num) == 0: - continue - temp_val.append(np.mean(img.img[self.mask == num])) - temp_idx.append(num) - all_val.append(temp_val) - all_idx.append(temp_idx) - all_val.append([np.mean(img.img[self.mask == -1])]) - all_idx.append([-1]) - return all_idx, all_val - - def debug(self, img: Image_Wrapper): - count = 1 - for ev_points in self.eval_points: - if count != 1 and False: - count += 1 - continue - for num in ev_points: - img.img[self.mask == num] += 10000 * num - count += 1 - return img.img - - def get_mask(self): - return self.mask - - @timeit - def generate_mask(self, img: Image_Wrapper, merge=False): - hash = str(img.img.shape) - self.mask = self.gen_mask_helper(img, hash=hash) - if merge: - self.mask = self.merge_mask_helper(hash=hash) - self.eval_points = [[a] for a in np.arange(len(self.eval_points))] - - -class Voronoi_Evaluator(Evaluator): - def __init__(self, list_points): - points = np.concatenate(list_points, axis=0) - self.eval_points = [] - start = 0 - for l in list_points: - stop = l.shape[0] - self.eval_points.append(np.arange(start, start + stop)) - start += stop - self.vor = Voronoi(points) - - @persist_to_file("cache_merge_voro") - def merge_mask_helper(self): - new_eval_points = np.arange(len(self.eval_points)) - mask = self.mask - for nc, ev_points in zip(new_eval_points, self.eval_points): - for num in ev_points: - mask[self.mask == num] = nc - return mask - - @persist_to_file("cache_voro") - def gen_mask_helper(self, img: Image_Wrapper): - mask = np.full_like(img.img, -1) - - counter = -1 - region_mask = self.vor.point_region - for i in np.array(self.vor.regions, dtype=list)[region_mask]: - counter += 1 - if -1 in i: - continue - if len(i) == 0: - continue - pts = self.vor.vertices[i] - pts = np.stack(img.val2pos( - pts[:, 0], pts[:, 1])).astype(np.int32).T - if np.any(pts < 0): - continue - mask_2 = np.zeros_like(img.img) - cv2.fillConvexPoly(mask_2, pts, 1) - mask_2 = mask_2 > 0 # To convert to Boolean - mask[mask_2] = counter - return mask - - -class Rect_Evaluator(Evaluator): - # def __init__(self, points, eval_points): - # self.eval_points = eval_points - # self.points = points - # self.length = 4 - - def __init__(self, list_points, length=4): - self.points = np.concatenate(list_points, axis=0) - self.eval_points = [] - start = 0 - for l in list_points: - stop = l.shape[0] - self.eval_points.append(np.arange(start, start + stop)) - start += stop - print(start) - print(self.points.shape) - self.length = length - - @persist_to_file("cache_merge_rect") - def merge_mask_helper(self): - new_eval_points = np.arange(len(self.eval_points)) - mask = self.mask - for nc, ev_points in zip(new_eval_points, self.eval_points): - for num in ev_points: - mask[self.mask == num] = nc - return mask - - @persist_to_file("cache_rect") - def gen_mask_helper(self, img: Image_Wrapper): - mask = np.full_like(img.img, -1) - count = 0 - for x, y in self.points: - x, y = img.val2pos(x, y) - x_lower = int(x - self.length) - y_lower = int(y - self.length) - x_upper = int(x + self.length + 1) - y_upper = int(y + self.length + 1) - if img.check_bounds(x_lower, y_lower, x_upper, y_upper): - x_lower, y_lower, x_upper, y_upper = img.clean_bounds( - x_lower, y_lower, x_upper, y_upper) - mask[y_lower:y_upper, x_lower:x_upper] = count - count += 1 - return mask - -# -# def main(): -# np.random.seed(10) -# points = (np.random.rand(100, 2)-0.5) * 2 -# voro = Voronoi_Evaluator(points, [[1],[2]]) -# rect = Rect_Evaluator(points, [[1], [2]]) -# Z = np.ones((1000, 1000)) -# img = Image_Wrapper(Z, -5, .01, -5, .01) -# voro.extract(img) -# rect.extract(img) -# -# plt.scatter(points[[1], 0], points[[1], 1]) -# plt.scatter(points[[2], 0], points[[2], 1]) -# plt.imshow(img.img, extent=img.ext(), origin="lower") -# #plt.imshow(img.img, origin="lower") -# plt.show() -# -# -# if __name__ == "__main__": -# main() diff --git a/2d_fourie/lattices.py b/2d_fourie/lattices.py deleted file mode 100644 index d25e21b..0000000 --- a/2d_fourie/lattices.py +++ /dev/null @@ -1,152 +0,0 @@ -import numpy as np -from cache import timeit - - -def deg_2_rad(winkel): - return winkel / 180.0 * np.pi - - -# all units in angstrom - - -class Lattice: - def __init__(self, x_len, y_len): - pass - - def get_from_mask(self, maske): - pass - - def get_both(self): - pass - - def reci(self): - pass - - -class SCC_Lattice(Lattice): - @timeit - def __init__(self, x_len, y_len): - x = np.arange(x_len) * 5 - y = np.arange(x_len) * 5 - self.X, self.Y = np.meshgrid(x, y) - - def get_from_mask(self, maske): - return self.X, self.Y - - def get_both(self): - return [(self.X, self.Y), (self.X, self.Y)] - - def reci(self): - x = np.arange(-3, 4) * 0.2 - y = np.arange(-3, 4) * 0.2 - X, Y = np.meshgrid(x, y) - return [(X, Y)] - - -class VO2_Lattice(Lattice): - base_a_m = 5.75 - base_b_m = 4.5 - base_c_m = 5.38 - - base_c_r = 2.856 - base_c_r = 2.8 - base_b_r = 4.554 - base_b_r = 4.5 - base_a_r = base_b_r - - alpha_m = 122.64 # degree - - def _mono_2_rutile(self, c_m, a_m): - a_r = np.cos(deg_2_rad(self.alpha_m - 90)) * c_m * self.base_c_m - c_r = (a_m) * self.base_a_m + \ - np.sin(deg_2_rad(self.alpha_m - 90)) * c_m * self.base_c_m - return a_r, c_r - - def _get_rutile(self, X, Y): - self.atom_x_rut = X * self.base_c_r + \ - np.mod(Y, 4) * 0.5 * self.base_c_r - self.atom_y_rut = Y * 0.5 * self.base_a_r - - def _get_mono(self, X, Y): - offset_a_m = 0.25 - 0.23947 - offset_c_m = 0.02646 - - 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) - - print(offset_a_r, offset_c_r) - - self.atom_x_mono = offset_a_r + X * \ - 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_y_mono = offset_c_r + 0.5 * Y * self.base_a_r - self.atom_y_mono[np.mod(X, 2) == 0] -= 2 * offset_c_r - - def _generate_vec(self, x_len: int, y_len: int): - x = np.arange(x_len) - y = np.arange(y_len) - X, Y = np.meshgrid(x, y) - X[np.mod(Y, 4) == 3] = X[np.mod(Y, 4) == 3] - 1 - X[np.mod(Y, 4) == 2] = X[np.mod(Y, 4) == 2] - 1 - assert np.mod(x.size, 2) == 0 - assert np.mod(y.size, 2) == 0 - - return X, Y - - @timeit - def __init__(self, x_len: int, y_len: int): - X, Y = self._generate_vec(x_len * 2, y_len * 2) - self._get_mono(X, Y) - self._get_rutile(X, Y) - - def get_from_mask(self, maske: np.array): - inplace_pos_x = np.zeros_like(self.atom_x_mono) - inplace_pos_y = np.zeros_like(self.atom_x_mono) - - mask = np.empty_like(self.atom_x_mono, dtype=bool) - mask[0::2, 0::2] = maske - mask[1::2, 0::2] = maske - mask[0::2, 1::2] = maske - mask[1::2, 1::2] = maske - - inplace_pos_x[mask] = self.atom_x_rut[mask] - inplace_pos_y[mask] = self.atom_y_rut[mask] - mask = np.invert(mask) - inplace_pos_x[mask] = self.atom_x_mono[mask] - inplace_pos_y[mask] = self.atom_y_mono[mask] - return inplace_pos_x, inplace_pos_y - - def get_both(self): - return [(self.atom_x_rut, self.atom_y_rut), (self.atom_x_mono, self.atom_y_mono)] - - def reci_rutile(self): - num = 20 - x = np.arange(-num, num + 1) - y = np.arange(-num, num + 1) - X, Y = np.meshgrid(x, y) - return (X * 0.22 + Y * 0.44).flatten(), (X * 0.349).flatten() - - def reci_mono(self): - x, y = self.reci_rutile() - return x + 0.1083, y + 0.1719 - - def reci_mono_2(self): - x, y = self.reci_rutile() - return x - 0.1083, y + 0.1719 - - def reci(self): - cutoff = 5. - x, y = self.reci_rutile() - mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff) - p1 = (x[mask], y[mask]) - x, y = self.reci_mono() - mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff) - p2 = (x[mask], y[mask]) - x, y = self.reci_mono_2() - mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff) - p3 = (x[mask], y[mask]) - return [p1, p2, p3] diff --git a/2d_fourie/main.py b/2d_fourie/main.py deleted file mode 100644 index ceb3d6a..0000000 --- a/2d_fourie/main.py +++ /dev/null @@ -1,329 +0,0 @@ -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 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, evaluator=None): - if evaluator is not None: - img = Image_Wrapper(freqx,freqy, intens) - intens = evaluator.debug(img) - 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" - ) - self.add_peaks(ax_log) - 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 - - - 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]) - - 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]) - - img = - rect = Rect_Evaluator([all_rutile, all_mono, all_mono2]) - rect.generate_mask(img) - plot.plot(freqx=fx, freqy=fy, intens=intens_rutile, - ax_log=axs[1, 0], ax_lin=axs[2, 0], evaluator=rect) - plot.plot(freqx=fx, freqy=fy, intens=intens_mono, - ax_log=axs[1, 2], ax_lin=axs[2, 2], evaluator=rect) - plot.plot(freqx=fx, freqy=fy, intens=intens_mixed, - ax_log=axs[1, 1], ax_lin=axs[2, 1], evaluator=rect) - - # 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 - plot = Plotter(lat) - [(x, y), (x1, y1)] = lat.get_both() - si = SpinImage_Two_Phase(x, y, x1, y1) - 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 - - si.apply_mask(maske) - si.gaussian(20) - - fx, fy, intens = si.fft() - img = Image_Wrapper(intens, fx, fy) - #if not already_inited: - # print("start_init") - # 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) - 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__": - np.random.seed(1234) - # test_me() - # test_square() - test_mixed() - plt.show() - #random(1234) - # for i in np.random.randint(0, 10000, 1): - # random(i) - # ising(i) - - # plt.show() diff --git a/2d_fourie/plot.py b/2d_fourie/plot.py deleted file mode 100644 index 5b685d0..0000000 --- a/2d_fourie/plot.py +++ /dev/null @@ -1,5 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np - -if __name__ == "__main__": - pass diff --git a/2d_fourie/spin_image.py b/2d_fourie/spin_image.py deleted file mode 100644 index d2d26b4..0000000 --- a/2d_fourie/spin_image.py +++ /dev/null @@ -1,270 +0,0 @@ -import numpy as np -import scipy -import scipy.fftpack as sfft -import scipy.signal -import tqdm -from cache import timeit - - -class SpinImage_Two_Phase: - resolution = 0.05 - offset = 40 - - def make_list(self, x_pos, y_pos, x_inds, y_inds): - sigma = .1 - x_ind = np.arange(0, self.length_x, self.resolution) - y_ind = np.arange(0, self.length_y, self.resolution) - X, Y = np.meshgrid(x_ind, y_ind, indexing="ij") - out_list = [] - for x, y, x_ind, y_ind in zip( - x_pos.flatten(), - y_pos.flatten(), - x_inds.flatten(), - y_inds.flatten(), - ): - xl, yl, xu, yu = self.clean_bounds( - x_ind - self.offset, - y_ind - self.offset, - x_ind + self.offset, - y_ind + self.offset, - X.shape, - ) - out_list.append(np.exp(-0.5 * ((X[xl:xu, yl:yu] - x) ** 2 + - (Y[xl:xu, yl:yu] - y) ** 2) / sigma**2)) - #out_list.append(np.ones_like(X[xl:xu, yl:yu])) - out_list = np.array(out_list) - return out_list - - @timeit - def __init__(self, x_pos_low, y_pos_low, x_pos_high, y_pos_high): - assert x_pos_low.shape == y_pos_low.shape - assert x_pos_high.shape == y_pos_high.shape - assert x_pos_low.shape == x_pos_high.shape - - offset_shift = self.offset * self.resolution - - zero_shift_x = np.minimum(np.min(x_pos_low), np.min(x_pos_high)) - zero_shift_y = np.minimum(np.min(y_pos_low), np.min(y_pos_high)) - - x_pos_low = x_pos_low - zero_shift_x + offset_shift - y_pos_low = y_pos_low - zero_shift_y + offset_shift - - x_pos_high = x_pos_high - zero_shift_x + offset_shift - y_pos_high = y_pos_high - zero_shift_y + offset_shift - - max_len_x = np.maximum(np.max(x_pos_low), np.max(x_pos_high)) - max_len_y = np.maximum(np.max(y_pos_low), np.max(y_pos_high)) - - self.length_x = max_len_x + (self.offset + 1) * self.resolution - self.length_y = max_len_y + (self.offset + 1) * self.resolution - self.x_index_low, self.y_index_low = self._stuff(x_pos_low, y_pos_low) - self.x_index_high, self.y_index_high = self._stuff( - x_pos_high, y_pos_high) - - self.buffer_low = self.make_list( - x_pos_low, y_pos_low, self.x_index_low, self.y_index_low) - self.buffer_high = self.make_list( - x_pos_high, y_pos_high, self.x_index_high, self.y_index_high) - - def clean_bounds(self, xl, yl, xu, yu, shape): - if xl < 0: - xl = 0 - if yl < 0: - yl = 0 - - if xu > shape[0]: - xu = shape[0] - if yu > shape[1]: - yu = shape[1] - - return xl, yl, xu, yu - - def _stuff(self, pos_x, pos_y): - x_ind = np.arange(0, self.length_x, self.resolution) - y_ind = np.arange(0, self.length_y, self.resolution) - xind = np.searchsorted(x_ind, pos_x).astype(int) - yind = np.searchsorted(y_ind, pos_y).astype(int) - return xind, yind - - def _apply_mask(self, x, y, buffer, mask): - for x, y, dat in zip( - x.flatten()[mask], - y.flatten()[mask], - buffer[mask], - ): - xl, yl, xu, yu = self.clean_bounds( - x - self.offset, - y - self.offset, - x + self.offset, - y + self.offset, - self.img.shape, - ) - self.img[x - self.offset: x + self.offset, - y - self.offset: y + self.offset] += dat - - @timeit - def apply_mask(self, maske): - mask = np.empty_like(self.x_index_high, dtype=bool) - print(maske) - mask[0::2, 0::2] = maske - mask[1::2, 0::2] = maske - mask[0::2, 1::2] = maske - mask[1::2, 1::2] = maske - mask = mask.flatten() - self.img = np.zeros( - (int(self.length_x/self.resolution), int(self.length_y/self.resolution))) - self._apply_mask(self.x_index_low, self.y_index_low, - self.buffer_low, mask) - mask = np.invert(mask) - self._apply_mask(self.x_index_high, self.y_index_high, - self.buffer_high, mask) - - @timeit - def fft(self): - Z_fft = sfft.fft2(self.img) - Z_shift = sfft.fftshift(Z_fft) - fft_freqx = sfft.fftfreq(self.img.shape[0], self.resolution) - fft_freqy = sfft.fftfreq(self.img.shape[1], self.resolution) - fft_freqx_clean = sfft.fftshift(fft_freqx) - fft_freqy_clean = sfft.fftshift(fft_freqy) - return fft_freqx_clean, fft_freqy_clean, np.abs(Z_shift) ** 2 - - @timeit - def plot(self, ax, scale=None): - if scale is None: - ax.imshow(self.img) - else: - quad = np.ones((int(scale / self.resolution), - int(scale / self.resolution))) - img = scipy.signal.convolve2d(self.img, quad) - ax.imshow(img) - - @timeit - def pad_it_square(self, additional_pad=0, size=None): - h = self.img.shape[0] - w = self.img.shape[1] - xx = np.maximum(h, w) + 2 * additional_pad + 1 - if size is not None: - xx = np.maximum(xx, size) - yy = xx - - 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 gaussian(self, sigma): - x = np.linspace(0, self.length_x, self.img.shape[0]) - y = np.linspace(0, self.length_y, self.img.shape[1]) - X, Y = np.meshgrid(x, y) - X = X - (self.length_x / 2.) - Y = Y - (self.length_y / 2.) - z = 1 / (2 * np.pi * sigma * sigma) * \ - np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2))) - self.img = np.multiply(self.img, z.T) - - -class SpinImage: - resolution = 0.05 - - def __init__(self, x_pos, y_pos): - x_pos = x_pos - np.min(x_pos) - y_pos = y_pos - np.min(y_pos) - self.length_x = np.max(x_pos) + self.resolution - self.length_y = np.max(y_pos) + self.resolution - self.img = self.image_from_pos(x_pos, y_pos) - - def image_from_pos(self, pos_x, pos_y): - x_ind = np.arange(0, self.length_x, self.resolution) # angstrom - y_ind = np.arange(0, self.length_y, self.resolution) # angstrom - img = np.zeros((x_ind.size, y_ind.size)) - xind = np.searchsorted(x_ind, pos_x) - yind = np.searchsorted(y_ind, pos_y) - img[xind, yind] = 1 - return img - - def image_from_pos_with_gauss(self, pos_x, pos_y, sigma=1): - x_ind = np.arange(0, self.length_x, self.resolution) # angstrom - y_ind = np.arange(0, self.length_y, self.resolution) # angstrom - img = np.zeros((x_ind.size, y_ind.size)) - X, Y = np.meshgrid(y_ind, x_ind) - for px, py in tqdm.tqdm(zip(pos_x.flatten(), pos_y.flatten())): - img += np.exp(-0.5 * ((X - px) ** 2 + (Y - py) ** 2) / sigma**2) - return img - - def fft(self): - Z_fft = sfft.fft2(self.img) - Z_shift = sfft.fftshift(Z_fft) - fft_freqx = sfft.fftfreq(self.img.shape[0], self.resolution) - fft_freqy = sfft.fftfreq(self.img.shape[1], self.resolution) - fft_freqx_clean = sfft.fftshift(fft_freqx) - fft_freqy_clean = sfft.fftshift(fft_freqy) - return fft_freqx_clean, fft_freqy_clean, np.abs(Z_shift) ** 2 - - def pad_it_square(self, additional_pad=0, size=None): - h = self.img.shape[0] - w = self.img.shape[1] - xx = np.maximum(h, w) + 2 * additional_pad + 1 - if size is not None: - xx = np.maximum(xx, size) - yy = xx - 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 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): - x = np.linspace(-self.length_x / 2, self.length_x / 2, mask.shape[0]) - y = np.linspace(-self.length_y / 2, self.length_y / 2, mask.shape[1]) - X, Y = np.meshgrid(x, y) - z = 1 / (2 * np.pi * sigma * sigma) * \ - np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2))) - return np.multiply(mask, z.T) - - def gaussian(self, sigma): - x = np.arange(-self.length_x / 2, self.length_x / 2, self.resolution) - y = np.arange(-self.length_y / 2, self.length_y / 2, self.resolution) - X, Y = np.meshgrid(x, y) - z = 1 / (2 * np.pi * sigma * sigma) * \ - np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2))) - self.img = np.multiply(self.img, z.T) - - def plot(self, ax, scale=None): - if scale is None: - ax.imshow(self.img) - else: - quad = np.ones((int(scale / self.resolution), - int(scale / self.resolution))) - img = scipy.signal.convolve2d(self.img, quad) - ax.imshow(img) - - def blur(self, sigma): - self.img = scipy.ndimage.gaussian_filter(self.img, sigma) diff --git a/2d_fourie/test_voronoi.py b/2d_fourie/test_voronoi.py deleted file mode 100644 index 9ba2806..0000000 --- a/2d_fourie/test_voronoi.py +++ /dev/null @@ -1,46 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -from scipy.spatial import Voronoi, voronoi_plot_2d -import cv2 - - - -class Voronoi_Evaulator: - def __init__(self, points, eval_points): - self.eval_points = eval_points - self.vor = Voronoi(points) - - def extract(self, Z): - for debug_num, ev_points in zip([-10, -5], self.eval_points): - region_mask = self.vor.point_region[ev_points] - print(region_mask) - for i in np.array(self.vor.regions)[region_mask]: - if -1 in i: - print("Containse outside points") - continue - if len(i) == 0: - print("Containse outside points") - continue - print(i) - pts = self.vor.vertices[i] - pts = (pts * 100).astype(np.int32) - print(pts) - mask = np.zeros((Z.shape[0], Z.shape[1])) - cv2.fillConvexPoly(mask, pts, 1) - mask = mask > 0 # To convert to Boolean - Z[mask] = debug_num - return Z - -if __name__ == "__main__": - np.random.seed(20) - points = (np.random.rand(100, 2)-0.1) * 2 - voro = Voronoi_Evaulator(points, [[1, 4, 5], [2, 3, 6]]) - - x = np.linspace(0, 1, 100) - y = np.linspace(0, 1, 200) - X, Y = np.meshgrid(x, y) - Z = X*2 + Y - - Z = voro.extract(Z) - plt.imshow(Z) - plt.show() diff --git a/rust/Cargo.toml b/rust/Cargo.toml deleted file mode 100644 index 97daa1e..0000000 --- a/rust/Cargo.toml +++ /dev/null @@ -1,11 +0,0 @@ -[package] -name = "rust" -version = "0.1.0" -edition = "2021" - -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html - -[dependencies] -ndarray = "0.15.6" -ndarray-npy = "0.8.1" -plotters = "0.3.4" diff --git a/rust/plot.py b/rust/plot.py deleted file mode 100644 index ee88bbf..0000000 --- a/rust/plot.py +++ /dev/null @@ -1,9 +0,0 @@ -import numpy as np - -import matplotlib.pyplot as plt - -file = np.load("./rutile.npz") -plt.scatter(file["x"], file["y"]) -file = np.load("./mono.npz") -plt.scatter(file["x"], file["y"]) -plt.show() diff --git a/rust/src/main.rs b/rust/src/main.rs deleted file mode 100644 index 1b4cbd1..0000000 --- a/rust/src/main.rs +++ /dev/null @@ -1,19 +0,0 @@ -mod vo2; -use vo2::get_mono; - -use crate::vo2::get_rutile; -use ndarray_npy::NpzWriter; -use std::fs::File; -fn main() { - let (mono_x, mono_y) = get_mono(20, 20); - let (rutile_x, rutile_y) = get_rutile(20, 20); - let mut npz = NpzWriter::new(File::create("mono.npz").unwrap()); - npz.add_array("x", &mono_x); - npz.add_array("y", &mono_y); - npz.finish().unwrap(); - - let mut npz = NpzWriter::new(File::create("rutile.npz").unwrap()); - npz.add_array("x", &rutile_x); - npz.add_array("y", &rutile_y); - npz.finish().unwrap(); -} diff --git a/rust/src/vo2.rs b/rust/src/vo2.rs deleted file mode 100644 index 743418d..0000000 --- a/rust/src/vo2.rs +++ /dev/null @@ -1,143 +0,0 @@ -use ndarray::{Array, Array2}; - -const BASE_A_M: f64 = 5.75; -const BASE_B_M: f64 = 4.5; -const BASE_C_M: f64 = 5.38; - -const BASE_C_R: f64 = 2.856; -const BASE_B_R: f64 = 4.554; -const BASE_A_R: f64 = BASE_B_R; - -const ALPHA_M: f64 = 122.64; // degree - -fn deg_2_rad(angle_degree: f64) -> f64 { - return std::f64::consts::PI / 180. * angle_degree; -} - -fn mono_2_rutile(c_m: f64, a_m: f64) -> (f64, f64) { - let a_r = deg_2_rad(ALPHA_M - 90.).cos() * c_m * BASE_C_M; - let c_r = (a_m) * BASE_A_M + deg_2_rad(ALPHA_M - 90.).sin() * c_m * BASE_C_M; - return (a_r, c_r); -} -pub fn get_rutile(len_x: usize, len_y: usize) -> (Array2, Array2) { - let mut rutile_x_pos: Array2 = Array::zeros((len_x * 2, len_y * 2)); - let mut rutile_y_pos: Array2 = Array::zeros((len_x * 2, len_y * 2)); - for x in 0..2 * len_x { - for y in 0..2 * len_y { - let tmp_x: f64 = if y % 4 >= 2 { x as f64 - 1. } else { x as f64 }; - rutile_x_pos[[x, y]] = tmp_x as f64 * BASE_C_R + (y % 4) as f64 * 0.5 * BASE_C_R; - rutile_y_pos[[x, y]] = y as f64 * 0.5 * BASE_A_R; - } - } - return (rutile_x_pos, rutile_y_pos); -} - -pub fn get_mono(len_x: usize, len_y: usize) -> (Array2, Array2) { - const OFFSET_A_M: f64 = 0.25 - 0.23947; - const OFFSET_C_M: f64 = 0.02646; - let (offset_a_r, offset_c_r) = mono_2_rutile(OFFSET_C_M, OFFSET_A_M); - let mut mono_x_pos: Array2 = Array::zeros((len_x * 2, len_y * 2)); - let mut mono_y_pos: Array2 = Array::zeros((len_x * 2, len_y * 2)); - for x in 0..2 * len_x { - for y in 0..2 * len_y { - let tmp_x: i64 = if y % 4 >= 2 { x as i64 - 1 } else { x as i64 }; - let offset_sign = if tmp_x % 2 == 0 { -1. } else { 1. }; - mono_x_pos[[x, y]] = offset_sign * 2. * offset_a_r - + tmp_x as f64 * BASE_C_R - + (y % 4) as f64 * 0.5 * BASE_C_R; - mono_y_pos[[x, y]] = offset_sign * 2. * offset_c_r + 0.5 * y as f64 * BASE_A_R; - } - } - return (mono_x_pos, mono_y_pos); -} /* - def _mono_2_rutile(self, c_m, a_m): - return a_r, c_r - - def _get_rutile(self, X, Y): - self.atom_x_rut = X * self.base_c_r + \ - np.mod(Y, 4) * 0.5 * self.base_c_r - self.atom_y_rut = Y * 0.5 * self.base_a_r - - def _get_mono(self, X, Y): - offset_a_m = 0.25 - 0.23947 - offset_ndc_m = 0.02646 - - 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) - - print(offset_a_r, offset_c_r) - - self.atom_x_mono = offset_a_r + X * \ - 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_y_mono = offset_c_r + 0.5 * Y * self.base_a_r - self.atom_y_mono[np.mod(X, 2) == 0] -= 2 * offset_c_r - - def _generate_vec(self, x_len: int, y_len: int): - x = np.arange(x_len) - y = np.arange(y_len) - X, Y = np.meshgrid(x, y) - X[np.mod(Y, 4) == 3] = X[np.mod(Y, 4) == 3] - 1 - X[np.mod(Y, 4) == 2] = X[np.mod(Y, 4) == 2] - 1 - assert np.mod(x.size, 2) == 0 - assert np.mod(y.size, 2) == 0 - - return X, Y - - @timeit - def __init__(self, x_len: int, y_len: int): - X, Y = self._generate_vec(x_len * 2, y_len * 2) - self._get_mono(X, Y) - self._get_rutile(X, Y) - - def get_from_mask(self, maske: np.array): - inplace_pos_x = np.zeros_like(self.atom_x_mono) - inplace_pos_y = np.zeros_like(self.atom_x_mono) - - mask = np.empty_like(self.atom_x_mono, dtype=bool) - mask[0::2, 0::2] = maske - mask[1::2, 0::2] = maske - mask[0::2, 1::2] = maske - mask[1::2, 1::2] = maske - - inplace_pos_x[mask] = self.atom_x_rut[mask] - inplace_pos_y[mask] = self.atom_y_rut[mask] - mask = np.invert(mask) - inplace_pos_x[mask] = self.atom_x_mono[mask] - inplace_pos_y[mask] = self.atom_y_mono[mask] - return inplace_pos_x, inplace_pos_y - - def reci_rutile(self): - num = 20 - #num = 2 - x = np.arange(-num, num + 1) - y = np.arange(-num, num + 1) - X, Y = np.meshgrid(x, y) - return (X * 0.22 + Y * 0.44).flatten(), (X * 0.349).flatten() - - def reci_mono(self): - x, y = self.reci_rutile() - return x + 0.1083, y + 0.1719 - - def reci_mono_2(self): - x, y = self.reci_rutile() - return x - 0.1083, y + 0.1719 - - def reci(self): - cutoff = 5. - x, y = self.reci_rutile() - mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff) - p1 = (x[mask], y[mask]) - x, y = self.reci_mono() - mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff) - p2 = (x[mask], y[mask]) - x, y = self.reci_mono_2() - mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff) - p3 = (x[mask], y[mask]) - return [p1, p2, p3] - - */