From 458998d1f7f0235410995b6804bb06e9b2e33e58 Mon Sep 17 00:00:00 2001 From: Jacob Date: Wed, 19 Apr 2023 16:54:25 +0200 Subject: [PATCH 1/3] fixed bug --- clean_python/extractors.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/clean_python/extractors.py b/clean_python/extractors.py index 234b7af..82ee0cf 100644 --- a/clean_python/extractors.py +++ b/clean_python/extractors.py @@ -1,12 +1,13 @@ -import numpy as np +import logging +from abc import abstractmethod, ABC +from tools import clean_bounds_offset +from spin_image import FFT +from cache import persist_to_file, timeit import tqdm +import numpy as np +import matplotlib.pyplot as plt # from scipy.spatial import Voronoi # import cv2 -from cache import persist_to_file, timeit -from spin_image import FFT -from tools import clean_bounds_offset -from abc import abstractmethod, ABC -import logging logger = logging.getLogger("fft") @@ -76,11 +77,17 @@ class Rect_Evaluator(Evaluator): def merge_mask_helper(self): new_eval_points = np.arange(len(self.eval_points)) - mask = self.mask + mask = self.mask.copy() for nc, ev_points in zip(new_eval_points, self.eval_points): - maske_low = np.min(ev_points) >= self.mask - maske_high = np.max(ev_points) <= self.mask - mask[np.logical_and(maske_high, maske_low)] = nc + maske_low = np.min(ev_points) <= self.mask + maske_high = np.max(ev_points) >= self.mask + maske = np.logical_and(maske_high, maske_low) + print("Debug: ", np.sum(maske)) + mask[maske] = nc + plt.figure() + plt.imshow(mask) + plt.show() + return mask def gen_mask_helper(self, img: FFT): From a68450421310af2acbd480ad25ccd36cbcde2f39 Mon Sep 17 00:00:00 2001 From: Jacob Date: Wed, 19 Apr 2023 17:34:32 +0200 Subject: [PATCH 2/3] Fixed counter --- clean_python/main.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/clean_python/main.py b/clean_python/main.py index c78e84d..07f06e4 100644 --- a/clean_python/main.py +++ b/clean_python/main.py @@ -122,9 +122,11 @@ def random(seed): for i in tqdm.tqdm(ind): maske[np.unravel_index(i, (LEN, LEN))] = True counter += 1 - if np.mod(counter, 100) != 0 and not i == ind[-1]: + if np.mod(counter, 100) != 0 and i != ind[-1] and i != ind[0]: continue + print(counter, i) + si.apply_mask(lat.parse_mask(maske)) si.gaussian(20) From e1a921c2eb7fb8f51d860e28b81ff3a41af21abc Mon Sep 17 00:00:00 2001 From: Jacob Date: Thu, 20 Apr 2023 19:49:28 +0200 Subject: [PATCH 3/3] Some small changes --- clean_python/analysis.py | 68 +++++++++++++------ clean_python/ditact_pic.py | 136 +++++++++++++++++++++++++++++++++++++ clean_python/extractors.py | 5 -- clean_python/lattices.py | 27 ++++++++ clean_python/main.py | 59 ++++++++++++---- clean_python/plotter.py | 6 +- clean_python/spin_image.py | 26 +++++-- 7 files changed, 278 insertions(+), 49 deletions(-) create mode 100644 clean_python/ditact_pic.py diff --git a/clean_python/analysis.py b/clean_python/analysis.py index 1389e48..29ab3b9 100644 --- a/clean_python/analysis.py +++ b/clean_python/analysis.py @@ -13,25 +13,36 @@ def check_percentage(p1, p2): def merge(files): merge = [] + plt.figure() for file in files: data = np.load(file, allow_pickle=True) old_percentage = data["percentage"] w_percentage = data["w_percentage"] # check_percentage(old_percentage, w_percentage) - percentage = w_percentage + percentage = old_percentage out = [] for o in ["out_1", "out_2", "out_3", "out_4"]: out.append(np.array(data[o])) out = np.array(out)[:, :, 0] + summe = np.max(np.sum(out, axis=0)) out = out / summe merge.append(out) - print(merge) - merge = sum(merge) - summe = np.max(np.sum(merge, axis=0)) - merge = merge / summe - print(merge) - return percentage, merge + + plt.plot(out[0, :], "r") + plt.plot(out[3, :], "b") + plt.plot(out[2, :], "g") + + all = sum(merge) + + summe = np.max(np.sum(all, axis=0)) + all = all / summe + + plt.plot(all[0, :], "k") + plt.plot(all[3, :], "k") + plt.plot(all[2, :], "k") + percentage = 1-percentage + return percentage, all def debug(percentage, out): @@ -39,15 +50,19 @@ def debug(percentage, out): for o in out: plt.plot(percentage, o) + plt.plot(percentage, out[0, :], "k") + plt.plot(percentage, out[3, :], "k") + plt.plot(percentage, out[2, :], "k") + def stacked_plot(percentage, out, title=""): plt.figure() - stacks = plt.stackplot(percentage, out[[0, 3, 1, 2], :], colors=[ + stacks = plt.stackplot(percentage, out[[0, 3, 1, 2]], colors=[ "w"], ls="solid", ec="k") hatches = ["/", "", "\\", "\\"] for stack, hatch in zip(stacks, hatches): stack.set_hatch(hatch) - plt.xlabel("Insulating Phase (%)") + plt.xlabel("Metallic Phase (%)") plt.ylabel("normalized Intensity ") plt.ylim([0.4, 1]) plt.xlim([0., 1]) @@ -56,6 +71,8 @@ def stacked_plot(percentage, out, title=""): plt.text(0.6, 0.5, "rutile", backgroundcolor="w") plt.text(0.35, 0.75, "diffusive", backgroundcolor="w") plt.title(title) + plt.savefig("intens.png") + plt.savefig("intens.pdf") def time_scale(p, o): @@ -66,30 +83,43 @@ def time_scale(p, o): mono_perc = -o[2] mono_perc = mono_perc - np.min(mono_perc) mono_perc /= np.max(mono_perc) - - cs_rut = ip.CubicSpline(p[::-1], rut_perc[::-1]) - cs_mono = ip.CubicSpline(p[::-1], mono_perc[::-1]) + + # cs_rut = ip.CubicSpline(p[::-1], rut_perc[::-1]) + # cs_mono = ip.CubicSpline(p[::-1], mono_perc[::-1]) + cs_rut = ip.interp1d(p[::-1], rut_perc[::-1]) + cs_mono = ip.interp1d(p[::-1], mono_perc[::-1]) plt.figure() - ph = np.linspace(0, 1, 100) + ph = np.linspace(0.01, 0.99, 100) plt.plot(ph, cs_rut(ph)) plt.plot(ph, cs_mono(ph)) - time = np.linspace(0, 3, 1000) + time = np.linspace(0.01, 3, 1000) phy_phase = np.exp(-time) rut_phase = cs_rut(phy_phase) mono_phase = cs_mono(phy_phase) plt.figure() - plt.plot(time, phy_phase) - plt.plot(time, rut_phase) - plt.plot(time, mono_phase) + plt.plot(time, phy_phase, "k:", label="corr.") + plt.plot(time, rut_phase, label="rut.") + plt.plot(time, mono_phase, label="mono") + plt.xlabel("time (a.u.)") + plt.ylabel("Metallic Phase (%)") + plt.legend() + plt.tight_layout() + plt.savefig("timescale.png") + plt.savefig("timescale.pdf") + +def read_file(file): + files = np.load("./merged.npz") + p = files["p"] + o = files["o"] + return p, o if __name__ == "__main__": p, o = merge(sys.argv[1:]) - # eval_data_print(f) - stacked_plot(p, o) # debug(p, o) + stacked_plot(p, o) time_scale(p, o) plt.show() diff --git a/clean_python/ditact_pic.py b/clean_python/ditact_pic.py new file mode 100644 index 0000000..c523354 --- /dev/null +++ b/clean_python/ditact_pic.py @@ -0,0 +1,136 @@ +from plotter import Plotter +import matplotlib +from lattices import VO2_New +from spin_image import SpinImage, FFT +import numpy as np +import matplotlib.pyplot as plt + + +plt.style.use("two_column") + + +def simulate(): + LEN = 50 + lat = VO2_New(LEN, LEN) + plot = Plotter(lat) + si = SpinImage(lat.get_phases()) + mask_misk = np.ones((LEN, 2*LEN)) + ind = np.arange(mask_misk.size) + np.random.shuffle(ind) + mask_misk[np.unravel_index( + ind[:int(mask_misk.size/2)], mask_misk.shape)] = 0 + print(mask_misk.shape) + + si.apply_mask(lat.parse_mask(np.zeros((LEN, 2*LEN)))) + si.gaussian(20) + intens_mono = si.fft() + intens_mono.clean() + + si.apply_mask(lat.parse_mask(np.ones((LEN, 2*LEN)))) + si.gaussian(20) + intens_rutile = si.fft() + intens_rutile.clean() + + si.apply_mask(lat.parse_mask(mask_misk)) + # si.apply_mask(mask_misk) + si.gaussian(20) + intens_mixed = si.fft() + intens_mixed.clean() + + intens_rutile.save("intens_rut.npz") + intens_mono.save("intens_mono.npz") + intens_mixed.save("intens_mixed.npz") + + +def plot(fft, ax): + ax.imshow( + fft.intens, + extent=fft.extents(), + norm=matplotlib.colors.LogNorm(), + cmap="magma", + origin="lower" + ) + + +def plot_all(intens_rutile, intens_mono, intens_mixed): + fig, axs = plt.subplots(4, 2) + fig.set_figheight(6) + for ax in axs.flatten(): + ax.axis("off") + axs = axs[:, 1] + + ax = axs[0] + plot(intens_rutile, ax) + + ax = axs[1] + plot(intens_mono, ax) + y_shift = 0.175 + h_shift = 2*y_shift + l_shift = 0.108 + big_shift = l_shift + h_shift + c = plt.Circle((-l_shift, y_shift), radius=0.07, + label='patch', fill=False, ec="w", ls=":") + ax.add_patch(c) + c = plt.Circle((l_shift, -y_shift), radius=0.07, + label='patch', fill=False, ec="w", ls=":") + ax.add_patch(c) + + c = plt.Circle((-h_shift, -y_shift), radius=0.07, + label='patch', fill=False, ec="w", ls=":") + ax.add_patch(c) + c = plt.Circle((h_shift, y_shift), radius=0.07, + label='patch', fill=False, ec="w", ls=":") + ax.add_patch(c) + + ax = axs[2] + plot(intens_mixed, ax) + # c = matplotlib.patches.Ellipse((0., 0.), width=0.2, height=1.4, angle=45, + # label='patch', fill=False, ec="w", ls=":") + # height = 1.4 + # width = 0.0 + # c = matplotlib.patches.Rectangle((-width/2, -height/2), width=width, height=height, + # angle=45, rotation_point="center", + # label='patch', fill=False, ec="w", ls=":") + # ax.add_patch(c) + ax.plot([-1, 1], [1, -1], "w", linestyle=(0, (2, 6))) + + ax = axs[3] + plot(intens_mixed, ax) + ax.plot([-1, 1], [0, 0], "w", linestyle=(0, (2, 6))) + ax.plot([-1, 1], [2 * y_shift, 2 * y_shift], "w", linestyle=(0, (2, 6))) + ax.plot([-1, 1], [-2*y_shift, -2*y_shift], "w", linestyle=(0, (2, 6))) + + ax.plot([-100*l_shift + big_shift*.5, 100*l_shift + big_shift*.5], + [y_shift*100, -y_shift*100], "w", linestyle=(0, (2, 6))) + ax.plot([-100*l_shift - big_shift*.5, 100*l_shift - big_shift*.5], + [y_shift*100, -y_shift*100], "w", linestyle=(0, (2, 6))) + + for ax in axs: + ax.axis("off") + ax.set_xlim(-0.5, 0.5) + ax.set_ylim(-0.5, 0.5) + plt.tight_layout() + fig.savefig("erklaerbaer.pdf") + fig.savefig("erklaerbaer.png") + # Plotting cuts + + +def load(): + r = FFT() + r.load("intens_rut.npz") + mo = FFT() + mo.load("intens_mono.npz") + mi = FFT() + mi.load("intens_mixed.npz") + + return r, mo, mi + + +if __name__ == "__main__": + np.random.seed(1234) + simulate() + # np.savez("intens.npz", r=r, mo=mo, mi=mi) + r, mo, mi = load() + plot_all(r, mo, mi) + + plt.show() diff --git a/clean_python/extractors.py b/clean_python/extractors.py index 82ee0cf..712c583 100644 --- a/clean_python/extractors.py +++ b/clean_python/extractors.py @@ -82,12 +82,7 @@ class Rect_Evaluator(Evaluator): maske_low = np.min(ev_points) <= self.mask maske_high = np.max(ev_points) >= self.mask maske = np.logical_and(maske_high, maske_low) - print("Debug: ", np.sum(maske)) mask[maske] = nc - plt.figure() - plt.imshow(mask) - plt.show() - return mask def gen_mask_helper(self, img: FFT): diff --git a/clean_python/lattices.py b/clean_python/lattices.py index b4e3d0c..d013856 100644 --- a/clean_python/lattices.py +++ b/clean_python/lattices.py @@ -1,3 +1,4 @@ +import matplotlib.pyplot as plt import numpy as np from cache import timeit from abc import ABC, abstractmethod @@ -106,7 +107,12 @@ class VO2_Lattice(Lattice): offset_a_m = 0.25 - 0.23947 offset_c_m = 0.02646 + # offset_c_m = -offset_c_m + # offset_a_m = -offset_a_m + offset_a_r, offset_c_r = self._mono_2_rutile(offset_c_m, offset_a_m) + # offset_a_r = -offset_a_r + # offset_c_r = -offset_c_r res = 0.05 offset_a_r = res * int(offset_a_r/res) @@ -148,3 +154,24 @@ class VO2_Lattice(Lattice): def _reci_mono_2(self): x, y = self._reci_rutile() return x - 0.5 * 1. / self.base_c_r, y + 0.5 * 1./self.base_a_r + + +class VO2_New(VO2_Lattice): + # def parse_mask(self, mask: np.ndarray) -> np.ndarray: + # maske = np.empty((mask.shape[0]*2, mask.shape[1]*2)) + # maske[0::2, 0::2] = mask + # maske[1::2, 0::2] = mask + # maske[0::2, 1::2] = mask + # maske[1::2, 1::2] = mask + # maske[0::4, :] = np.roll(maske[0::4, :], axis=1, shift=1) + # maske[1::4, :] = np.roll(maske[1::4, :], axis=1, shift=1) + # return maske + #def parse_mask(self, mask: np.ndarray) -> np.ndarray: + # print(mask.shape) + # maske = np.empty((mask.shape[0]*2, mask.shape[1])) + # maske[0::2, :] = mask + # maske[1::2, :] = mask + # print(maske.shape) + # return maske + def parse_maske(self, mask: np.ndarray): + return mask diff --git a/clean_python/main.py b/clean_python/main.py index 07f06e4..9a9ec9e 100644 --- a/clean_python/main.py +++ b/clean_python/main.py @@ -1,15 +1,16 @@ +import logging +import scipy.fftpack as sfft +from plotter import Plotter +from scipy import signal +from cache import timeit +from extractors import Rect_Evaluator +import tqdm from lattices import SCC_Lattice, VO2_Lattice import sys from spin_image import SpinImage import numpy as np import matplotlib.pyplot as plt -import tqdm -from extractors import Rect_Evaluator -from cache import timeit -from scipy import signal -from plotter import Plotter -import scipy.fftpack as sfft -import logging +plt.style.use(["style", "colors", "two_column"]) logger = logging.getLogger('fft') # logger.setLevel(logging.DEBUG) ch = logging.StreamHandler() @@ -21,8 +22,9 @@ logger.addHandler(ch) def test_mixed(): + plt.style.use("one_column") fig, axs = plt.subplots(3, 3) - LEN = 40 + LEN = 50 lat = VO2_Lattice(LEN, LEN) plot = Plotter(lat) si = SpinImage(lat.get_phases()) @@ -62,9 +64,23 @@ def test_mixed(): ax_log=axs[1, 2], ax_lin=axs[2, 2]) plot.plot_fft(intens_mixed, ax_log=axs[1, 1], ax_lin=axs[2, 1]) + plt.figure() + fig, axs = plt.subplots(1,3) + fig.set_figheight(1.7) plot.plot_fft(intens_mixed, - ax_log=plt.gca()) + ax_log=axs[1]) + plot.plot_fft(intens_mono, + ax_log=axs[0]) + plot.plot_fft(intens_rutile, + ax_log=axs[2]) + for ax, t in zip(axs,["monoclinic", "mixed", "rutile"]): + ax.set_title(t) + ax.set_xlim(-1,1) + ax.set_ylim(-1,1) + plt.tight_layout() + fig.savefig("diff_pattern.pdf") + fig.savefig("diff_pattern.png") # Plotting cuts @@ -98,10 +114,26 @@ def test_pdf(): rect.purge(intens) plt.figure() plot.plot_fft(intens, ax_log=plt.gca()) + plt.xlim(-1, 1) + plt.ylim(-1, 1) + plt.savefig("diff.png") + plt.savefig("diff.pdf") + pdf = sfft.fft2(intens.intens) pdf = sfft.fftshift(pdf) + plt.figure() - plt.imshow(np.abs(pdf)) + plt.imshow(np.abs(pdf), vmax=100) + plt.xlabel("Pos") + plt.ylabel("Pos") + x = pdf.shape[1] / 2. + y = pdf.shape[0] / 2. + off = 100 + plt.xlim(x-off, x+off) + plt.ylim(y-off, y+off) + plt.tight_layout() + plt.savefig("pdf.pdf") + plt.savefig("pdf.png") def random(seed): @@ -221,11 +253,12 @@ def runner(): if __name__ == "__main__": np.random.seed(1234) - runner() + # runner() # test_me() # test_square() - # test_mixed() - # plt.show() + test_mixed() + # test_pdf() + plt.show() # random(1234) # ising(1234) # test_pdf() diff --git a/clean_python/plotter.py b/clean_python/plotter.py index 4a7c29e..660576a 100644 --- a/clean_python/plotter.py +++ b/clean_python/plotter.py @@ -39,11 +39,7 @@ class Plotter: origin="lower" ) plt.colorbar(t, ax=ax_log, extend="min") - - ax_log.set_xlim(-2, 2) - ax_log.set_ylim(-2, 2) - ax_log.set_xlim(-8, 8) - ax_log.set_ylim(-8, 8) + if ax_lin: t = ax_lin.imshow( fft.intens, diff --git a/clean_python/spin_image.py b/clean_python/spin_image.py index 81c04fe..5c0d117 100644 --- a/clean_python/spin_image.py +++ b/clean_python/spin_image.py @@ -11,24 +11,35 @@ logger = logging.getLogger("fft") @dataclass class FFT: - freqx: np.ndarray - freqy: np.ndarray - intens: np.ndarray + freqx: np.ndarray = np.array([]) + freqy: np.ndarray = np.array([]) + intens: np.ndarray = np.array([]) def extents(self): return (np.min(self.freqx), np.max(self.freqx), np.min(self.freqy), np.max(self.freqy)) def clean(self): + cutoff = 1e-10 total = np.sum(self.intens) - low = np.sum(self.intens[self.intens < 1e-12]) - print(f"{low}/{total}") - self.intens[self.intens < 1e-12] = 1e-12 + low = np.sum(self.intens[self.intens < cutoff]) + print(f"{low}/{total} = {low/total}") + self.intens[self.intens < cutoff] = cutoff def val2pos(self, x, y): xind = np.searchsorted(self.freqx, x).astype(int) yind = np.searchsorted(self.freqy, y).astype(int) return xind, yind + def save(self, filename): + np.savez(filename, intens=self.intens, + freqx=self.freqx, freqy=self.freqy) + + def load(self, filename): + files = np.load(filename) + self.freqx = files["freqx"] + self.freqy = files["freqy"] + self.intens = files["intens"] + class SpinImage: resolution = 0.05 @@ -144,7 +155,8 @@ class SpinImage: for idx, map in zip(self.true_pos, self.intens_map): (X, Y) = idx z = 1 / (2 * np.pi * sigma * sigma) * \ - np.exp(-((X - mu_x)**2 / (2 * sigma**2) + (Y-mu_y)**2 / (2 * sigma**2))) + np.exp(-((X - mu_x)**2 / (2 * sigma**2) + + (Y-mu_y)**2 / (2 * sigma**2))) map *= z def get_intens(self, mask):