This commit is contained in:
Jacob Holder 2023-05-03 09:51:45 +02:00
parent 405c25438e
commit 9152ad8e62
3 changed files with 185 additions and 53 deletions

View File

@ -1,11 +1,13 @@
from ditact_pic import plot
from spin_image import SpinImage, FFT
import sys import sys
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import glob import glob
import scipy.interpolate as ip import scipy.interpolate as ip
from lattices import VO2_Lattice
plt.style.use(["style", "colors", "one_column"]) plt.style.use(["style", "colors", "one_column"])
from spin_image import SpinImage, FFT
from ditact_pic import plot
def check_percentage(p1, p2): def check_percentage(p1, p2):
plt.figure() plt.figure()
@ -32,9 +34,9 @@ def merge(files):
out = out / summe out = out / summe
merge.append(out) 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) all = sum(merge)
summe = np.max(np.sum(all, axis=0)) summe = np.max(np.sum(all, axis=0))
@ -58,18 +60,23 @@ def debug(percentage, out):
def stacked_plot(ax, percentage, out, title=""): def stacked_plot(ax, percentage, out, title=""):
stacks = ax.stackplot(percentage, out[[0, 3, 1, 2]], colors=[ stacks = ax.stackplot(percentage, out[[0, 3, 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_hatch(hatch)
stack.set_edgecolor(color)
ax.set_xlabel("Metallic Phase (%)") ax.set_xlabel("Metallic Phase (%)")
ax.set_ylabel("normalized Intensity ") ax.set_ylabel("normalized Intensity ")
ax.set_ylim([0.4, 1]) ax.set_ylim([0.4, 1])
ax.set_xlim([0., 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, 3, 2]], colors=["None"], ec="k")
def time_scale(ax, p, o): def time_scale(ax, p, o):
@ -86,21 +93,23 @@ def time_scale(ax, p, o):
cs_rut = ip.interp1d(p[::-1], rut_perc[::-1]) cs_rut = ip.interp1d(p[::-1], rut_perc[::-1])
cs_mono = ip.interp1d(p[::-1], mono_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) time = np.linspace(0.01, 3, 1000)
phy_phase = np.exp(-time) phy_phase = 1-np.exp(-time)
rut_phase = cs_rut(phy_phase) rut_phase = cs_rut(phy_phase)
mono_phase = cs_mono(phy_phase) mono_phase = cs_mono(phy_phase)
ax.plot(time, phy_phase, "k:", label="physical") 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_xlabel("time (a.u.)")
ax.set_ylabel("Metallic Phase (%)") ax.set_ylabel("Metallic Phase (%)")
ax.set_xlim([0, 3])
ax.set_ylim([0, 1])
ax.legend() ax.legend()
@ -110,29 +119,66 @@ def read_file(file):
o = files["o"] o = files["o"]
return p, o return p, o
def intens(ax, file, p, o): def intens(ax, file, p, o):
intens = FFT() intens = FFT()
intens.load(file) intens.load(file)
plot(intens, ax) 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") #rect = plt.Rectangle((-1, -.8), 2, 1.6, facecolor="None", hatch="//")
axins.plot(p, o[2], label="mono") # ax.add_patch(rect)
axins.legend(loc='lower left', bbox_to_anchor=(1, 0.5)) lat = VO2_Lattice(20, 20)
#axins.get_yaxis().set_visible(False) reci = lat.get_spots()
axins.yaxis.tick_right() print(reci)
axins.set_yticks([0,0.5] 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[3], 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__": if __name__ == "__main__":
p, o = merge(sys.argv[2:]) p, o = merge(sys.argv[2:])
np.savez("merged.npz", p=p, o=o) np.savez("merged.npz", p=p, o=o)
# eval_data_print(f) # eval_data_print(f)
fig, axs = plt.subplots(1,3) fig, axs = plt.subplots(1, 3)
fig.set_figheight(3) fig.set_figheight(2)
stacked_plot(axs[1],p, o) stacked_plot(axs[1], p, o)
time_scale(axs[2],p, o) time_scale(axs[2], p, o)
intens(axs[0], sys.argv[1], p ,o) intens(axs[0], sys.argv[1], p, o)
plt.tight_layout() plt.tight_layout()
plt.savefig("analysis.pdf")
plt.savefig("analysis.png")
plt.show() plt.show()

View File

@ -6,7 +6,7 @@ import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
if __name__ == "__main__": if __name__ == "__main__":
plt.style.use(["style", "colors","two_column"]) plt.style.use(["style", "colors", "two_column"])
def simulate(): def simulate():
@ -52,13 +52,15 @@ def plot(fft, ax):
origin="lower" origin="lower"
) )
def norm(*intenss): def norm(*intenss):
max = 1e-10 max = 1e-10
for intens in intenss: for intens in intenss:
m = np.max(intens.intens) m = np.max(intens.intens)
max = np.maximum(max,m) max = np.maximum(max, m)
return max return max
def plot_all(intens_rutile, intens_mono, intens_mixed): def plot_all(intens_rutile, intens_mono, intens_mixed):
fig, axs = plt.subplots(3, 2) fig, axs = plt.subplots(3, 2)
fig.set_figheight(5.2) fig.set_figheight(5.2)
@ -75,19 +77,22 @@ def plot_all(intens_rutile, intens_mono, intens_mixed):
h_shift = 2*y_shift h_shift = 2*y_shift
l_shift = 0.108 l_shift = 0.108
big_shift = l_shift + h_shift 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] ax = axs[2]
cmap = plot(intens_mixed, ax) cmap = plot(intens_mixed, ax)
@ -98,12 +103,12 @@ def plot_all(intens_rutile, intens_mono, intens_mixed):
ax.set_xlim(-cut_off, cut_off) ax.set_xlim(-cut_off, cut_off)
ax.set_ylim(-cut_off, cut_off) ax.set_ylim(-cut_off, cut_off)
plt.tight_layout() 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_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,
orientation="horizontal", ticks=[1e-10, 1e-5, 1e0])
#cbar.ax.set_xticklabels(['Low', 'Medium', 'High']) #cbar.ax.set_xticklabels(['Low', 'Medium', 'High'])
fig.savefig("erklaerbaer.pdf") fig.savefig("erklaerbaer.pdf")
fig.savefig("erklaerbaer.png") fig.savefig("erklaerbaer.png")
# Plotting cuts # Plotting cuts
@ -122,10 +127,10 @@ def load():
if __name__ == "__main__": if __name__ == "__main__":
np.random.seed(1234) np.random.seed(1234)
#simulate() simulate()
# np.savez("intens.npz", r=r, mo=mo, mi=mi) # np.savez("intens.npz", r=r, mo=mo, mi=mi)
r, mo, mi = load() r, mo, mi = load()
max = norm(r,mo,mi) max = norm(r, mo, mi)
r.intens = r.intens/max r.intens = r.intens/max
mo.intens = mo.intens/max mo.intens = mo.intens/max
mi.intens = mi.intens/max mi.intens = mi.intens/max

81
clean_python/ising.py Normal file
View File

@ -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):
LEN = 120
#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_{seed}.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):
np.random.seed(1234)
print(f"runnig: {file}")
ising(file)
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])