jacob/FFT
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Some small changes

Browse Source
This commit is contained in:
Jacob Holder 2023-04-20 19:49:28 +02:00
parent a684504213
commit e1a921c2eb
Signed by: jacob
GPG Key ID: 2194FC747048A7FD
7 changed files with 278 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
clean_python/analysis.py
View File

@ -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()

136
clean_python/ditact_pic.py Normal file
View File

@ -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()

5
clean_python/extractors.py
View File

@ -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):

27
clean_python/lattices.py
View File

@ -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

59
clean_python/main.py
View File

@ -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()

6
clean_python/plotter.py
View File

@ -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,

26
clean_python/spin_image.py
View File

@ -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):