diff --git a/clean_python/analysis.py b/clean_python/analysis.py index 4bb1b08..3c4bda2 100644 --- a/clean_python/analysis.py +++ b/clean_python/analysis.py @@ -1,17 +1,87 @@ +from ditact_pic import plot +from spin_image import SpinImage, FFT import sys import numpy as np import matplotlib.pyplot as plt import glob import scipy.interpolate as ip -plt.style.use(["style", "colors", "one_column"]) from spin_image import SpinImage, FFT from ditact_pic import plot +from lattices import VO2_Lattice +plt.style.use(["style", "colors", "one_column"]) + def check_percentage(p1, p2): plt.figure() plt.plot(p1, p2) +def average_mean(arr, window_size=20): + arr_sum = np.cumsum(arr) + arr = (arr_sum[window_size:] - arr_sum[:-window_size]) / window_size + return arr + + +def new_merge(files): + wp = [] + op = [] + spot_1 = [] + spot_2 = [] + spot_3 = [] + plt.figure() + for file in files: + print(file) + data = np.load(file, allow_pickle=True) + old_percentage = data["percentage"] + w_percentage = data["w_percentage"] + wp.append(w_percentage) + op.append(old_percentage) + # check_percentage(old_percentage, w_percentage) + out = [] + for o in ["out_1", "out_2", "out_3", "out_4"]: + out.append(np.array(data[o])) + print(out) + out = np.array(out)[:, :, 0] + + spot_1.append(out[0, :]) + spot_2.append(out[3, :]) + spot_3.append(out[2, :]) + wp = np.concatenate(wp, axis=0) + op = np.concatenate(op, axis=0) + spot_1 = np.concatenate(spot_1, axis=0) + spot_2 = np.concatenate(spot_2, axis=0) + spot_3 = np.concatenate(spot_3, axis=0) + + arg_sort = np.argsort(op) + wp = wp[arg_sort] + op = op[arg_sort] + spot_1 = spot_1[arg_sort] + spot_2 = spot_2[arg_sort] + spot_3 = spot_3[arg_sort] + + win = 100 + wp = average_mean(wp, win) + op = average_mean(op, win) + spot_1 = average_mean(spot_1, win) + spot_2 = average_mean(spot_2, win) + spot_3 = average_mean(spot_3, win) + + x = op + plt.plot(x, spot_1, "r.") + plt.plot(x, spot_2, "g.") + plt.plot(x, spot_3, "b.") + + ma = np.max(spot_1+spot_2+spot_3) + spot_1 /= ma + spot_2 /= ma + spot_3 /= ma + + print("debug....") + print(wp.shape) + plt.savefig("debug.png") + return op, np.stack([spot_2, spot_1, spot_3]) + + def merge(files): merge = [] plt.figure() @@ -32,17 +102,18 @@ def merge(files): out = out / summe merge.append(out) - plt.plot(out[0, :], "r") - plt.plot(out[3, :], "b") - plt.plot(out[2, :], "g") + plt.plot(w_percentage, out[0, :], "r.") + plt.plot(w_percentage, out[3, :], "b.") + plt.plot(w_percentage, 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") + # plt.plot(all[0, :], "k") + # plt.plot(all[3, :], "k") + # plt.plot(all[2, :], "k") + plt.savefig("debug.png") percentage = 1-percentage return percentage, all @@ -58,18 +129,23 @@ def debug(percentage, out): def stacked_plot(ax, percentage, out, title=""): - stacks = ax.stackplot(percentage, out[[0, 3, 1, 2]], colors=[ - "w"], ls="solid", ec="k") - hatches = ["/", "", "\\", "\\"] - for stack, hatch in zip(stacks, hatches): + stacks = ax.stackplot(percentage, out[[0, 1, 2]], colors=[ + "w"], ls=(0, (0, 1)), ec="w") + hatches = ["//", "|", "\\\\"] + for stack, hatch, color in zip(stacks, hatches, ["C1", "C0", "C2"]): stack.set_hatch(hatch) + stack.set_edgecolor(color) ax.set_xlabel("Metallic Phase (%)") ax.set_ylabel("normalized Intensity ") ax.set_ylim([0.4, 1]) ax.set_xlim([0., 1]) - ax.text(0.1, 0.9, "monoclinic", backgroundcolor="w") - ax.text(0.6, 0.5, "rutile", backgroundcolor="w") - ax.text(0.35, 0.75, "diffusive", backgroundcolor="w") + ax.text(0.1, 0.9, "monoclinic", backgroundcolor="w", + bbox=dict(boxstyle='square,pad=0.0', ec="None", fc="w")) + ax.text(0.6, 0.5, "rutile", backgroundcolor="w", + bbox=dict(boxstyle='square,pad=0.0', ec="None", fc="w")) + ax.text(0.35, 0.73, "diffusive", backgroundcolor="w", + bbox=dict(boxstyle='square,pad=0.0', ec="None", fc="w")) + ax.stackplot(percentage, out[[0, 1, 2]], colors=["None"], ec="k") def time_scale(ax, p, o): @@ -86,21 +162,23 @@ def time_scale(ax, p, o): cs_rut = ip.interp1d(p[::-1], rut_perc[::-1]) cs_mono = ip.interp1d(p[::-1], mono_perc[::-1]) - #plt.figure() - #ph = np.linspace(0.01, 0.99, 100) - #plt.plot(ph, cs_rut(ph)) - #plt.plot(ph, cs_mono(ph)) + # plt.figure() + # ph = np.linspace(0.01, 0.99, 100) + # plt.plot(ph, cs_rut(ph)) + # plt.plot(ph, cs_mono(ph)) time = np.linspace(0.01, 3, 1000) - phy_phase = np.exp(-time) + phy_phase = 1-np.exp(-time) rut_phase = cs_rut(phy_phase) mono_phase = cs_mono(phy_phase) ax.plot(time, phy_phase, "k:", label="physical") - ax.plot(time, rut_phase, label="rutile") - ax.plot(time, mono_phase, label="monoclinic") + ax.plot(time, rut_phase, label="rutile", color="C1") + ax.plot(time, mono_phase, label="monoclinic", color="C2") ax.set_xlabel("time (a.u.)") ax.set_ylabel("Metallic Phase (%)") + ax.set_xlim([0, 3]) + ax.set_ylim([0, 1]) ax.legend() @@ -114,25 +192,61 @@ def intens(ax, file, p, o): intens = FFT() intens.load(file) plot(intens, ax) - ax.set_xlim([-.8,0.8]) - ax.set_ylim([-.8,1.6]) - axins = ax.inset_axes([0.0, 0.5, 0.47, 0.5]) - axins.plot(p, o[0], label="rut") - axins.plot(p, o[3], label="diff") - axins.plot(p, o[2], label="mono") - axins.legend(loc='lower left', bbox_to_anchor=(1, 0.5)) - #axins.get_yaxis().set_visible(False) - axins.yaxis.tick_right() - axins.set_yticks([0,0.5] -if __name__ == "__main__": - p, o = merge(sys.argv[2:]) - np.savez("merged.npz", p=p, o=o) - # eval_data_print(f) + ax.set_xlim([-1, 1]) + ax.set_ylim([-.9, .9]) + ax.axis("off") - fig, axs = plt.subplots(1,3) - fig.set_figheight(3) - stacked_plot(axs[1],p, o) - time_scale(axs[2],p, o) - intens(axs[0], sys.argv[1], p ,o) + # rect = plt.Rectangle((-1, -.8), 2, 1.6, facecolor="None", hatch="//") + # ax.add_patch(rect) + lat = VO2_Lattice(20, 20) + reci = lat.get_spots() + print(reci) + size = (intens.freqx[1] - intens.freqx[0]) * 20 + size2 = size/2 + # big_rect = plt.Rectangle((-10, -10), 20, 20, fc="None", ec="k", hatch="//") + # ax.add_patch(big_rect) + for x, y in zip(reci[0][0], reci[0][1]): + if x < 1 and x > -1: + if y < 1 and y > -1: + print(x, y) + rect = plt.Rectangle((-y-size2, x-size2), + size, size, fc="C1", ec="k", alpha=0.5) + # big_rect.set_clip_path(rect) + ax.add_patch(rect) + for x, y in zip(reci[1][0], reci[1][1]): + if x < 1 and x > -1: + if y < 1 and y > -1: + print(x, y) + rect = plt.Rectangle((-y-size2, x-size2), + size, size, fc="C2", ec="k", alpha=0.5) + ax.add_patch(rect) + + axins = ax.inset_axes([0.0, 0.0, 0.5, 0.5]) + axins.plot(p, o[0], label="rut.", color="C1") + axins.plot(p, o[2], label="mono.", color="C2") + axins.plot(p, o[1], label="diff.", color="C0") + axins.legend(loc='center left', bbox_to_anchor=(1, 0.5)) + axins.set_xlim([0, 1]) + axins.set_ylim([0, 1]) + axins.set_xlabel("phase (%)") + axins.set_ylabel("signal", labelpad=-5) + + # axins.get_yaxis().set_visible(False) + # axins.yaxis.tick_right() + axins.set_yticks([0, 1]) + + +if __name__ == "__main__": + p, o = new_merge(sys.argv[2:]) + np.savez("merged.npz", p=p, o=o) + + fig, axs = plt.subplots(1, 3) + fig.set_figheight(2) + stacked_plot(axs[1], p, o) + time_scale(axs[2], p, o) + if "intens" in sys.argv[1]: + intens(axs[0], sys.argv[1], p, o) plt.tight_layout() + plt.savefig("analysis.pdf") + plt.savefig("analysis.png") plt.show() diff --git a/clean_python/ditact_pic.py b/clean_python/ditact_pic.py index 92754b2..8425a83 100644 --- a/clean_python/ditact_pic.py +++ b/clean_python/ditact_pic.py @@ -6,7 +6,7 @@ import numpy as np import matplotlib.pyplot as plt if __name__ == "__main__": - plt.style.use(["style", "colors","two_column"]) + plt.style.use(["style", "colors", "two_column"]) def simulate(): @@ -47,18 +47,28 @@ def plot(fft, ax): fft.intens, extent=fft.extents(), norm=matplotlib.colors.LogNorm(vmin=1e-10, vmax=1), - #norm=matplotlib.colors.Normalize(vmax=1, vmin=1e-10), + # norm=matplotlib.colors.Normalize(vmax=1, vmin=1e-10), cmap="magma", origin="lower" ) + def norm(*intenss): max = 1e-10 for intens in intenss: m = np.max(intens.intens) - max = np.maximum(max,m) + max = np.maximum(max, m) return max + +def norm(*intenss): + max = 1e-10 + for intens in intenss: + m = np.max(intens.intens) + max = np.maximum(max, m) + return max + + def plot_all(intens_rutile, intens_mono, intens_mixed): fig, axs = plt.subplots(3, 2) fig.set_figheight(5.2) @@ -75,19 +85,22 @@ def plot_all(intens_rutile, intens_mono, intens_mixed): 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((-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) + # 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.annotate("", (-l_shift, y_shift), (0, 0), + arrowprops=dict(ec='w', facecolor='white', width=.8, headwidth=3, headlength=5)) ax = axs[2] cmap = plot(intens_mixed, ax) @@ -98,11 +111,11 @@ def plot_all(intens_rutile, intens_mono, intens_mixed): ax.set_xlim(-cut_off, cut_off) ax.set_ylim(-cut_off, cut_off) plt.tight_layout() - fig.subplots_adjust(bottom=0.1,right=0.95,left=0.15,wspace=0.) + fig.subplots_adjust(bottom=0.1, right=0.95, left=0.15, wspace=0.) cbar_ax = fig.add_axes([0.55, 0.07, 0.4, 0.015]) - cbar = fig.colorbar(cmap, cax=cbar_ax, orientation="horizontal", ticks=[1e-10, 1e-5, 1e0]) - #cbar.ax.set_xticklabels(['Low', 'Medium', 'High']) - + cbar = fig.colorbar(cmap, cax=cbar_ax, + orientation="horizontal", ticks=[1e-10, 1e-5, 1e0]) + # cbar.ax.set_xticklabels(['Low', 'Medium', 'High']) fig.savefig("erklaerbaer.pdf") fig.savefig("erklaerbaer.png") @@ -122,10 +135,10 @@ def load(): if __name__ == "__main__": np.random.seed(1234) - #simulate() + # simulate() # np.savez("intens.npz", r=r, mo=mo, mi=mi) r, mo, mi = load() - max = norm(r,mo,mi) + max = norm(r, mo, mi) r.intens = r.intens/max mo.intens = mo.intens/max mi.intens = mi.intens/max diff --git a/clean_python/ising.py b/clean_python/ising.py new file mode 100644 index 0000000..5b62dbb --- /dev/null +++ b/clean_python/ising.py @@ -0,0 +1,81 @@ +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, VO2_New +import sys +from spin_image import SpinImage +import numpy as np +import matplotlib.pyplot as plt +plt.style.use(["style", "colors", "two_column"]) +logger = logging.getLogger('fft') +# logger.setLevel(logging.DEBUG) +ch = logging.StreamHandler() +ch.setLevel(logging.DEBUG) +formatter = logging.Formatter( + '%(asctime)s - %(name)s - %(levelname)s - %(message)s') +ch.setFormatter(formatter) +logger.addHandler(ch) + + +def ising(file, num): + LEN = 60 + #lat = VO2_New(LEN, LEN) + lat = VO2_New(LEN, LEN) + rect = Rect_Evaluator(lat.get_spots()) + + out_rect = [[] for x in range(4)] + percentage = [] + weighted_percentage = [] + si = SpinImage(lat.get_phases()) + already_inited = False + + spins = np.load(file)["s1"] + spins[spins==-1] = 0 + + for i in tqdm.tqdm(range(1000)): + + (ix,iy) = np.random.randint(2000-LEN-LEN,size=2) + maske = spins[ix:ix+2*LEN, iy:iy+2*LEN] + + si.apply_mask(lat.parse_mask(maske)) + si.gaussian(20) + + intens = si.fft() + if not already_inited: + rect.generate_mask(intens, merge=True) + already_inited = True + + ir, vr = rect.extract(intens) + for lis, val in zip(out_rect, vr): + lis.append(val) + percentage.append(np.sum(maske)) + [p1, p2] = si.get_intens(lat.parse_mask(maske)) + weighted_percentage.append(p1/(p1+p2)) + + percentage = np.array(percentage) + weighted_percentage = np.array(weighted_percentage) + percentage /= np.max(percentage) + + np.savez(f"ising_rect_{num}.npz", + w_percentage=weighted_percentage, percentage=percentage, out_1=out_rect[0], + out_2=out_rect[1], out_3=out_rect[2], out_4=out_rect[3]) + +def runner(file, idx): + np.random.seed(1234) + print(f"runnig: {file}") + ising(file,idx) + + +if __name__ == "__main__": + files = sys.argv[2:] + idx = int(sys.argv[1]) + print(f"{idx}/{len(files)}") + if idx > len(files): + exit() + if idx < 1: + exit() + runner(files[idx-1], idx)