rebased

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(w_percentage, out[0, :], "r.")
-        plt.plot(out[3, :], "b")
+        plt.plot(w_percentage, out[3, :], "b.")
-        plt.plot(out[2, :], "g")
+        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[0, :], "k")
-    plt.plot(all[3, :], "k")
+    # plt.plot(all[3, :], "k")
-    plt.plot(all[2, :], "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=[
+    stacks = ax.stackplot(percentage, out[[0, 1, 2]], colors=[
-                           "w"], ls="solid", ec="k")
+        "w"], ls=(0, (0, 1)), ec="w")
-    hatches = ["/", "", "\\", "\\"]
+    hatches = ["//", "|", "\\\\"]
-    for stack, hatch in zip(stacks, 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.1, 0.9, "monoclinic", backgroundcolor="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.75, "diffusive", backgroundcolor="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()
+    # plt.figure()
-    #ph = np.linspace(0.01, 0.99, 100)
+    # ph = np.linspace(0.01, 0.99, 100)
-    #plt.plot(ph, cs_rut(ph))
+    # plt.plot(ph, cs_rut(ph))
-    #plt.plot(ph, cs_mono(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, rut_phase, label="rutile", color="C1")
-    ax.plot(time, mono_phase, label="monoclinic")
+    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_xlim([-1, 1])
-    ax.set_ylim([-.8,1.6])
+    ax.set_ylim([-.9, .9])
-    axins = ax.inset_axes([0.0, 0.5, 0.47, 0.5])
+    ax.axis("off")
-    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)
 
-    fig, axs = plt.subplots(1,3)
+    # rect = plt.Rectangle((-1, -.8), 2, 1.6, facecolor="None", hatch="//")
-    fig.set_figheight(3)
+    # ax.add_patch(rect)
-    stacked_plot(axs[1],p, o)
+    lat = VO2_Lattice(20, 20)
-    time_scale(axs[2],p, o)
+    reci = lat.get_spots()
-    intens(axs[0], sys.argv[1], p ,o)
+    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()

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,
+    # c = plt.Circle((-l_shift, y_shift), radius=0.07,
-                   label='patch', fill=False, ec="w", ls=":")
+    #               label='patch', fill=False, ec="w", ls=":")
-    ax.add_patch(c)
+    # ax.add_patch(c)
-    c = plt.Circle((l_shift, -y_shift), radius=0.07,
+    # c = plt.Circle((l_shift, -y_shift), radius=0.07,
-                   label='patch', fill=False, ec="w", ls=":")
+    #               label='patch', fill=False, ec="w", ls=":")
-    ax.add_patch(c)
+    # ax.add_patch(c)
 
-    c = plt.Circle((-h_shift, -y_shift), radius=0.07,
+    # c = plt.Circle((-h_shift, -y_shift), radius=0.07,
-                   label='patch', fill=False, ec="w", ls=":")
+    #               label='patch', fill=False, ec="w", ls=":")
-    ax.add_patch(c)
+    # ax.add_patch(c)
-    c = plt.Circle((h_shift, y_shift), radius=0.07,
+    # c = plt.Circle((h_shift, y_shift), radius=0.07,
-                   label='patch', fill=False, ec="w", ls=":")
+    #               label='patch', fill=False, ec="w", ls=":")
-    ax.add_patch(c)
+    # 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 = fig.colorbar(cmap, cax=cbar_ax,
-    #cbar.ax.set_xticklabels(['Low', 'Medium', 'High'])
+                        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

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)