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

removed old stuff

Browse Source
This commit is contained in:
Jacob Holder 2023-04-27 15:40:38 +02:00
parent 42e3995501
commit 405c25438e
Signed by: jacob
GPG Key ID: 2194FC747048A7FD
12 changed files with 0 additions and 1370 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

105
2d_fourie/analysis.py
View File

@ -1,105 +0,0 @@
import sys
import numpy as np
import matplotlib.pyplot as plt
import glob
def eval_data(file):
data = np.load(file)
percentage = data["percentage"]
out = data["out"]
out = np.array(out)
print(out.shape)
fig, all_axs = plt.subplots(2, out.shape[0])
axs = all_axs[0, :]
axs2 = all_axs[1, :]
for o, ax, ax2, lab in zip(out, axs, axs2, ["rutile", "mono_twin", "mono"]):
# ax.plot(percentage, o/np.max(o, axis=0))
ax.plot(percentage, o/o[0])
# ax.plot(percentage, o)
o = np.mean(o, axis=1)
o = o/o[0]
ax2.plot(percentage, o)
ax2.plot([0, 1], [o[0], o[-1]], "k:")
ax.set_title(lab)
title = ""
if "ising" in file:
title += "Ising "
else:
title += "Random "
if "rect" in file:
title += "rect "
if "voro" in file:
title += "voro "
fig.suptitle(title)
fig.savefig(f"{file}.png")
def parse_lists(out):
lists = []
for o in out:
lists.append(np.stack(o))
max = 0
for l in lists:
print(l.shape)
if max < l.shape[1]:
max = l.shape[1]
lists = [np.pad(l, ((0, 0), (0, max-l.shape[1]))) for l in lists]
for l in lists:
print(l.shape)
return np.stack(lists)
def eval_data_print(file):
data = np.load(file, allow_pickle=True)
percentage = data["percentage"]
# out = parse_lists(data["out"])
out = []
for o in ["out_1", "out_2", "out_3", "out_4"]:
out.append(np.array(data[o]))
fig, all_axs = plt.subplots(2, len(out))
axs = all_axs[0, :]
axs2 = all_axs[1, :]
for o, ax, ax2, lab in zip(out, axs, axs2, ["rutile", "monoclinic", "mono", "rest"]):
# ax.plot(percentage, o/np.max(o, axis=0))
ax.plot(percentage, o/o[0])
# ax.plot(percentage, o)
o = np.mean(o, axis=1)
# o = o/o[0]
ax2.plot(percentage, o)
ax2.plot([0, 1], [o[0], o[-1]], "k:")
ax.set_title(lab)
for ax in all_axs.flatten():
ax.set_xlabel("Rutile Phase")
ax.set_ylabel("Normalized Intensity")
title = ""
if "ising" in file:
title += "Ising "
else:
title += "Random "
if "rect" in file:
title += "rect "
if "voro" in file:
title += "voro "
fig.suptitle(title)
plt.tight_layout()
def stacked_plot(file):
data = np.load(file, allow_pickle=True)
percentage = data["percentage"]
# out = parse_lists(data["out"])
out = []
for o in ["out_1", "out_2", "out_3", "out_4"]:
out.append(np.array(data[o]))
out = np.array(out)
plt.stackplot(percentage, out)
if __name__ == "__main__":
for f in sys.argv[1:]:
#eval_data_print(f)
stacked_plot(f)
plt.show()

68
2d_fourie/cache.py
View File

@ -1,68 +0,0 @@
import numpy as np
from functools import wraps
import time
def timeit(func):
@wraps(func)
def timeit_wrapper(*args, **kwargs):
start_time = time.perf_counter()
print(f"Start Function {func.__name__}:")
result = func(*args, **kwargs)
end_time = time.perf_counter()
total_time = end_time - start_time
# first item in the args, ie `args[0]` is `self`
print(
f'Function {func.__name__} Took {total_time:.4f} seconds')
return result
return timeit_wrapper
def persist_to_file(file_name):
def decorator(original_func):
try:
file_nam = file_name
if file_name[-4:] != ".npz":
file_nam += ".npz"
file = np.load(file_nam)
cache = dict(zip((file.files), (file[k] for k in file.files)))
except (IOError, ValueError):
cache = {}
def hash_func(*param):
key = str(param)
return key
def persist_func(*param, hash=None):
if hash is None:
hash = hash_func(*param)
print("Hash: ", hash)
if cache == {} or ("hash" not in cache) or hash != cache["hash"]:
print("recalc")
data = original_func(*param)
np.savez(file_name, hash=hash, dat=data)
cache["hash"] = hash
cache["dat"] = data
print("loaded")
return cache["dat"]
return persist_func
return decorator
# test = 1234
#
#
# @persist_to_file("cache.npz", test)
# def test():
# print("calculate")
# return np.zeros((100, 100))
#
#
# if __name__ == "__main__":
# test()
# test()
# pass

213
2d_fourie/extractors.py
View File

@ -1,213 +0,0 @@
import numpy as np
from scipy.spatial import Voronoi
import cv2
from cache import persist_to_file, timeit
class Image_Wrapper:
# def __init__(self, img, x_lower, x_res, y_lower, y_res):
# self.img = img
# self.x_lower = x_lower
# self.y_lower = y_lower
# self.x_res = x_res
# self.y_res = y_res
# self.x_upper = self.x_lower + self.img.shape[0]*self.x_res - self.x_res
# self.y_upper = self.y_lower + self.img.shape[1]*self.y_res - self.x_res
def __init__(self, img, fx, fy):
self.img = img
self.x_lower = np.min(fx)
self.y_lower = np.min(fy)
self.x_upper = np.max(fx)
self.y_upper = np.max(fy)
self.x_res = (self.x_upper - self.x_lower) / self.img.shape[0]
self.y_res = (self.y_upper - self.y_lower) / self.img.shape[0]
def val2pos(self, x, y):
x = (x - self.x_lower) / self.x_res
y = (y - self.y_lower) / self.y_res
return x, y
def check_bounds(self, xl, yl, xu, yu):
if xl > self.img.shape[0]:
print("xl lim")
return False
if yl > self.img.shape[1]:
print("yl lim")
return False
if xu < 0:
print("xu lim")
return False
if yu < 0:
print("yu lim")
return False
return True
def clean_bounds(self, xl, yl, xu, yu):
if xl < 0:
xl = 0
if yl < 0:
yl = 0
if xu > self.img.shape[0]:
xu = self.img.shape[0]
if yu > self.img.shape[1]:
yu = self.img.shape[1]
return xl, yl, xu, yu
def ext(self):
return [self.x_lower, self.x_upper, self.y_lower, self.y_upper]
class Evaluator:
def extract(self, img: Image_Wrapper):
all_val = []
all_idx = []
for ev_points in self.eval_points:
temp_val = []
temp_idx = []
for num in ev_points:
if np.sum(self.mask == num) == 0:
continue
temp_val.append(np.mean(img.img[self.mask == num]))
temp_idx.append(num)
all_val.append(temp_val)
all_idx.append(temp_idx)
all_val.append([np.mean(img.img[self.mask == -1])])
all_idx.append([-1])
return all_idx, all_val
def debug(self, img: Image_Wrapper):
count = 1
for ev_points in self.eval_points:
if count != 1 and False:
count += 1
continue
for num in ev_points:
img.img[self.mask == num] += 10000 * num
count += 1
return img.img
def get_mask(self):
return self.mask
@timeit
def generate_mask(self, img: Image_Wrapper, merge=False):
hash = str(img.img.shape)
self.mask = self.gen_mask_helper(img, hash=hash)
if merge:
self.mask = self.merge_mask_helper(hash=hash)
self.eval_points = [[a] for a in np.arange(len(self.eval_points))]
class Voronoi_Evaluator(Evaluator):
def __init__(self, list_points):
points = np.concatenate(list_points, axis=0)
self.eval_points = []
start = 0
for l in list_points:
stop = l.shape[0]
self.eval_points.append(np.arange(start, start + stop))
start += stop
self.vor = Voronoi(points)
@persist_to_file("cache_merge_voro")
def merge_mask_helper(self):
new_eval_points = np.arange(len(self.eval_points))
mask = self.mask
for nc, ev_points in zip(new_eval_points, self.eval_points):
for num in ev_points:
mask[self.mask == num] = nc
return mask
@persist_to_file("cache_voro")
def gen_mask_helper(self, img: Image_Wrapper):
mask = np.full_like(img.img, -1)
counter = -1
region_mask = self.vor.point_region
for i in np.array(self.vor.regions, dtype=list)[region_mask]:
counter += 1
if -1 in i:
continue
if len(i) == 0:
continue
pts = self.vor.vertices[i]
pts = np.stack(img.val2pos(
pts[:, 0], pts[:, 1])).astype(np.int32).T
if np.any(pts < 0):
continue
mask_2 = np.zeros_like(img.img)
cv2.fillConvexPoly(mask_2, pts, 1)
mask_2 = mask_2 > 0 # To convert to Boolean
mask[mask_2] = counter
return mask
class Rect_Evaluator(Evaluator):
# def __init__(self, points, eval_points):
# self.eval_points = eval_points
# self.points = points
# self.length = 4
def __init__(self, list_points, length=4):
self.points = np.concatenate(list_points, axis=0)
self.eval_points = []
start = 0
for l in list_points:
stop = l.shape[0]
self.eval_points.append(np.arange(start, start + stop))
start += stop
print(start)
print(self.points.shape)
self.length = length
@persist_to_file("cache_merge_rect")
def merge_mask_helper(self):
new_eval_points = np.arange(len(self.eval_points))
mask = self.mask
for nc, ev_points in zip(new_eval_points, self.eval_points):
for num in ev_points:
mask[self.mask == num] = nc
return mask
@persist_to_file("cache_rect")
def gen_mask_helper(self, img: Image_Wrapper):
mask = np.full_like(img.img, -1)
count = 0
for x, y in self.points:
x, y = img.val2pos(x, y)
x_lower = int(x - self.length)
y_lower = int(y - self.length)
x_upper = int(x + self.length + 1)
y_upper = int(y + self.length + 1)
if img.check_bounds(x_lower, y_lower, x_upper, y_upper):
x_lower, y_lower, x_upper, y_upper = img.clean_bounds(
x_lower, y_lower, x_upper, y_upper)
mask[y_lower:y_upper, x_lower:x_upper] = count
count += 1
return mask
#
# def main():
# np.random.seed(10)
# points = (np.random.rand(100, 2)-0.5) * 2
# voro = Voronoi_Evaluator(points, [[1],[2]])
# rect = Rect_Evaluator(points, [[1], [2]])
# Z = np.ones((1000, 1000))
# img = Image_Wrapper(Z, -5, .01, -5, .01)
# voro.extract(img)
# rect.extract(img)
#
# plt.scatter(points[[1], 0], points[[1], 1])
# plt.scatter(points[[2], 0], points[[2], 1])
# plt.imshow(img.img, extent=img.ext(), origin="lower")
# #plt.imshow(img.img, origin="lower")
# plt.show()
#
#
# if __name__ == "__main__":
# main()

152
2d_fourie/lattices.py
View File

@ -1,152 +0,0 @@
import numpy as np
from cache import timeit
def deg_2_rad(winkel):
return winkel / 180.0 * np.pi
# all units in angstrom
class Lattice:
def __init__(self, x_len, y_len):
pass
def get_from_mask(self, maske):
pass
def get_both(self):
pass
def reci(self):
pass
class SCC_Lattice(Lattice):
@timeit
def __init__(self, x_len, y_len):
x = np.arange(x_len) * 5
y = np.arange(x_len) * 5
self.X, self.Y = np.meshgrid(x, y)
def get_from_mask(self, maske):
return self.X, self.Y
def get_both(self):
return [(self.X, self.Y), (self.X, self.Y)]
def reci(self):
x = np.arange(-3, 4) * 0.2
y = np.arange(-3, 4) * 0.2
X, Y = np.meshgrid(x, y)
return [(X, Y)]
class VO2_Lattice(Lattice):
base_a_m = 5.75
base_b_m = 4.5
base_c_m = 5.38
base_c_r = 2.856
base_c_r = 2.8
base_b_r = 4.554
base_b_r = 4.5
base_a_r = base_b_r
alpha_m = 122.64 # degree
def _mono_2_rutile(self, c_m, a_m):
a_r = np.cos(deg_2_rad(self.alpha_m - 90)) * c_m * self.base_c_m
c_r = (a_m) * self.base_a_m + \
np.sin(deg_2_rad(self.alpha_m - 90)) * c_m * self.base_c_m
return a_r, c_r
def _get_rutile(self, X, Y):
self.atom_x_rut = X * self.base_c_r + \
np.mod(Y, 4) * 0.5 * self.base_c_r
self.atom_y_rut = Y * 0.5 * self.base_a_r
def _get_mono(self, X, Y):
offset_a_m = 0.25 - 0.23947
offset_c_m = 0.02646
offset_a_r, offset_c_r = self._mono_2_rutile(offset_c_m, offset_a_m)
res = 0.05
offset_a_r = res * int(offset_a_r/res)
offset_c_r = res * int(offset_c_r/res)
print(offset_a_r, offset_c_r)
self.atom_x_mono = offset_a_r + X * \
self.base_c_r + np.mod(Y, 4) * 0.5 * self.base_c_r
self.atom_x_mono[np.mod(X, 2) == 0] -= 2 * offset_a_r
self.atom_y_mono = offset_c_r + 0.5 * Y * self.base_a_r
self.atom_y_mono[np.mod(X, 2) == 0] -= 2 * offset_c_r
def _generate_vec(self, x_len: int, y_len: int):
x = np.arange(x_len)
y = np.arange(y_len)
X, Y = np.meshgrid(x, y)
X[np.mod(Y, 4) == 3] = X[np.mod(Y, 4) == 3] - 1
X[np.mod(Y, 4) == 2] = X[np.mod(Y, 4) == 2] - 1
assert np.mod(x.size, 2) == 0
assert np.mod(y.size, 2) == 0
return X, Y
@timeit
def __init__(self, x_len: int, y_len: int):
X, Y = self._generate_vec(x_len * 2, y_len * 2)
self._get_mono(X, Y)
self._get_rutile(X, Y)
def get_from_mask(self, maske: np.array):
inplace_pos_x = np.zeros_like(self.atom_x_mono)
inplace_pos_y = np.zeros_like(self.atom_x_mono)
mask = np.empty_like(self.atom_x_mono, dtype=bool)
mask[0::2, 0::2] = maske
mask[1::2, 0::2] = maske
mask[0::2, 1::2] = maske
mask[1::2, 1::2] = maske
inplace_pos_x[mask] = self.atom_x_rut[mask]
inplace_pos_y[mask] = self.atom_y_rut[mask]
mask = np.invert(mask)
inplace_pos_x[mask] = self.atom_x_mono[mask]
inplace_pos_y[mask] = self.atom_y_mono[mask]
return inplace_pos_x, inplace_pos_y
def get_both(self):
return [(self.atom_x_rut, self.atom_y_rut), (self.atom_x_mono, self.atom_y_mono)]
def reci_rutile(self):
num = 20
x = np.arange(-num, num + 1)
y = np.arange(-num, num + 1)
X, Y = np.meshgrid(x, y)
return (X * 0.22 + Y * 0.44).flatten(), (X * 0.349).flatten()
def reci_mono(self):
x, y = self.reci_rutile()
return x + 0.1083, y + 0.1719
def reci_mono_2(self):
x, y = self.reci_rutile()
return x - 0.1083, y + 0.1719
def reci(self):
cutoff = 5.
x, y = self.reci_rutile()
mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff)
p1 = (x[mask], y[mask])
x, y = self.reci_mono()
mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff)
p2 = (x[mask], y[mask])
x, y = self.reci_mono_2()
mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff)
p3 = (x[mask], y[mask])
return [p1, p2, p3]

329
2d_fourie/main.py
View File

@ -1,329 +0,0 @@
from lattices import SCC_Lattice, VO2_Lattice
from spin_image import SpinImage, SpinImage_Two_Phase
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib
import tqdm
from extractors import Rect_Evaluator, Voronoi_Evaluator, Image_Wrapper
from cache import timeit
from scipy import signal
class Plotter:
def __init__(self, lat):
self.lattice = lat
self.length_2 = 0.05
def rect_at_point(self, x, y, color):
length_2 = self.length_2
rect = patches.Rectangle(
(x - length_2, y - length_2),
2 * length_2,
2 * length_2,
linewidth=1,
edgecolor=color,
facecolor="none",
)
return rect
def plot(self, freqx, freqy, intens, ax_log=None, ax_lin=None, vmax=None, evaluator=None):
if evaluator is not None:
img = Image_Wrapper(freqx,freqy, intens)
intens = evaluator.debug(img)
if ax_log:
t = ax_log.imshow(
intens,
extent=(np.min(freqx), np.max(freqx),
np.min(freqy), np.max(freqy)),
norm=matplotlib.colors.LogNorm(vmin=1e-12),
cmap="viridis",
origin="lower"
)
self.add_peaks(ax_log)
plt.colorbar(t, ax=ax_log)
if ax_lin:
t = ax_lin.imshow(
intens,
extent=(np.min(freqx), np.max(freqx),
np.min(freqy), np.max(freqy)),
# vmax=vmax,
cmap="viridis",
origin="lower"
)
plt.colorbar(t, ax=ax_lin)
def plot_periodic(array):
plt.figure()
kernel = np.ones((10, 10))
#array = signal.convolve2d(array, kernel, boundary='symm', mode='same')
tiled = np.tile(array, (2, 2))
plt.imshow(tiled)
def rotate(x, y, angle):
radian = angle / 180 * 2 * np.pi
return np.cos(radian) * x - np.sin(radian) * y,\
np.sin(radian) * x + np.cos(radian) * y
def test_square():
LEN = 40
lat = SCC_Lattice(LEN, LEN)
# lat = VO2_Lattice(LEN, LEN)
plot = Plotter(lat)
pos_x, pos_y = lat.get_from_mask(np.zeros((40, 40)))
[(x, y), (x1, y1)] = lat.get_both()
# pos_x, pos_y = rotate(pos_x, pos_y, 30)
# si = SpinImage(pos_x, pos_y)
si = SpinImage_Two_Phase(x, y, x1, y1)
si.apply_mask(np.zeros((20, 20)))
fig, axs = plt.subplots(2, 2)
si.plot(axs[0, 0])
# si.gaussian(LEN)
# si.blur(3)
# si.pad_it_square(additional_pad=1000)
si.gaussian(20)
si.plot(axs[0, 1])
plt.pause(0.1)
fx, fy, intens = si.fft()
plot.plot(fx, fy, intens, axs[1, 0], axs[1, 1])
plt.tight_layout()
plt.savefig("test.pdf")
plt.figure()
index = np.abs(fy).argmin()
plt.plot(fx, intens[index, :], label="fy=0")
plt.plot(fx, intens[index+1, :], label="fy=df")
plt.plot(fx, intens[index+2, :], label="fy=2df")
plt.legend()
plt.yscale("log")
plt.show()
def helper(intens, fx, fy, voro, rect):
print(np.mean(intens))
v_idx, v_val = voro.extract(Image_Wrapper(intens, fx, fy))
r_idx, r_val = rect.extract(Image_Wrapper(intens, fx, fy))
for iv, av, ir, ar in zip(v_idx, v_val, r_idx, r_val):
av = np.array(av)
ar = np.array(ar)
mask = np.isin(ir, iv)
print("Test: ", np.all(np.isclose(
ar[mask], av)), np.isclose(ar[mask], av))
print(iv, "\n", ar[mask], "\n", av, "\n\n\n")
def test_mixed():
fig, axs = plt.subplots(3, 3)
LEN = 80
lat = VO2_Lattice(LEN, LEN)
plot = Plotter(lat)
[(x, y), (x1, y1)] = lat.get_both()
si = SpinImage_Two_Phase(x, y, x1, y1)
mask_misk = np.ones((LEN, LEN), dtype=bool)
ind = np.arange(mask_misk.size)
np.random.shuffle(ind)
mask_misk[np.unravel_index(ind[:800], (LEN, LEN))] = False
all_rutile = np.stack(lat.reci()[0]).T
all_mono = np.stack(lat.reci()[1]).T
all_mono2 = np.stack(lat.reci()[2]).T
#voro = Voronoi_Evaluator([all_rutile, all_mono, all_mono2])
si.apply_mask(np.zeros((LEN, LEN), dtype=bool))
si.gaussian(20)
fx, fy, intens_mono = si.fft()
si.plot(axs[0, 2])
si.apply_mask(np.ones((LEN, LEN)))
si.gaussian(20)
fx, fy, intens_rutile = si.fft()
si.plot(axs[0, 0])
si.apply_mask(mask_misk)
si.gaussian(20)
fx, fy, intens_mixed = si.fft()
si.plot(axs[0, 1])
img =
rect = Rect_Evaluator([all_rutile, all_mono, all_mono2])
rect.generate_mask(img)
plot.plot(freqx=fx, freqy=fy, intens=intens_rutile,
ax_log=axs[1, 0], ax_lin=axs[2, 0], evaluator=rect)
plot.plot(freqx=fx, freqy=fy, intens=intens_mono,
ax_log=axs[1, 2], ax_lin=axs[2, 2], evaluator=rect)
plot.plot(freqx=fx, freqy=fy, intens=intens_mixed,
ax_log=axs[1, 1], ax_lin=axs[2, 1], evaluator=rect)
# Plotting cuts
def random(seed):
np.random.seed(seed)
LEN = 80
lat = VO2_Lattice(LEN, LEN)
maske = np.zeros((LEN, LEN))
ind = np.arange(LEN * LEN)
np.random.shuffle(ind)
all_rutile = np.stack(lat.reci()[0]).T
all_mono = np.stack(lat.reci()[1]).T
all_mono2 = np.stack(lat.reci()[2]).T
voro = Voronoi_Evaluator([all_rutile, all_mono, all_mono2])
rect = Rect_Evaluator([all_rutile, all_mono, all_mono2])
out_rect = [[] for x in range(len(lat.reci())+1)]
out_voro = [[] for x in range(len(lat.reci())+1)]
percentage = []
counter = 0
plot = Plotter(lat)
[(x, y), (x1, y1)] = lat.get_both()
si = SpinImage_Two_Phase(x, y, x1, y1)
already_inited = False
for i in tqdm.tqdm(ind):
maske[np.unravel_index(i, (LEN, LEN))] = True
counter += 1
if np.mod(counter, 100) != 0:
continue
si.apply_mask(maske)
si.gaussian(20)
fx, fy, intens = si.fft()
img = Image_Wrapper(intens, fx, fy)
#if not already_inited:
# print("start_init")
# voro.generate_mask(img, merge=True)
# rect.generate_mask(img, merge=True)
# already_inited = True
#iv, vv = voro.extract(img)
#ir, vr = rect.extract(img)
#for lis, val in zip(out_rect, vr):
# lis.append(val)
#for lis, val in zip(out_voro, vv):
# lis.append(val)
#percentage.append(np.sum(maske))
percentage = np.array(percentage)
percentage /= np.max(percentage)
np.savez(f"random_rect_{seed}.npz",
percentage=percentage, out_1=out_rect[0],
out_2=out_rect[1], out_3=out_rect[2], out_4=out_rect[3])
np.savez(f"random_voro_{seed}.npz",
percentage=percentage, out_1=out_voro[0],
out_2=out_voro[1], out_3=out_voro[2], out_4=out_voro[3])
def sample_index(p):
i = np.random.choice(np.arange(p.size), p=p.ravel())
return np.unravel_index(i, p.shape)
@timeit
def ising(seed):
np.random.seed(seed)
LEN = 40
temp = 0.1
maske = np.zeros((LEN, LEN), dtype=bool)
lat = VO2_Lattice(LEN, LEN)
all_rutile = np.stack(lat.reci()[0]).T
all_mono = np.stack(lat.reci()[1]).T
all_mono2 = np.stack(lat.reci()[2]).T
voro = Voronoi_Evaluator([all_rutile, all_mono, all_mono2])
rect = Rect_Evaluator([all_rutile, all_mono, all_mono2])
out_rect = [[] for x in range(len(lat.reci())+1)]
out_voro = [[] for x in range(len(lat.reci())+1)]
percentage = []
counter = 0
already_inited = False
for i in tqdm.tqdm(range(LEN*LEN)):
probability = np.roll(maske, 1, axis=0).astype(float)
probability += np.roll(maske, -1, axis=0).astype(float)
probability += np.roll(maske, 1, axis=1).astype(float)
probability += np.roll(maske, -1, axis=1).astype(float)
probability = np.exp(probability/temp)
probability[maske] = 0
probability /= np.sum(probability)
maske[sample_index(probability)] = True
counter += 1
if np.mod(counter, 100) != 0:
continue
pos_x, pos_y = lat.get_from_mask(maske)
si = SpinImage(pos_x, pos_y)
si.pad_it_square(10, size=2300)
fx, fy, intens = si.fft()
img = Image_Wrapper(intens, fx, fy)
if not already_inited:
voro.generate_mask(img, merge=True)
rect.generate_mask(img, merge=True)
already_inited = True
iv, vv = voro.extract(img)
ir, vr = rect.extract(img)
for lis, val in zip(out_rect, vr):
lis.append(val)
for lis, val in zip(out_voro, vv):
lis.append(val)
percentage.append(np.sum(maske))
percentage = np.array(percentage, dtype=np.float64)
percentage /= np.max(percentage)
np.savez(f"ising_rect_{seed}.npz",
percentage=percentage, out_1=out_rect[0],
out_2=out_rect[1], out_3=out_rect[2], out_4=out_rect[3])
np.savez(f"ising_voro_{seed}.npz",
percentage=percentage, out_1=out_voro[0],
out_2=out_voro[1], out_3=out_voro[2], out_4=out_voro[3])
def test_me():
LEN = 20
lat = VO2_Lattice(LEN, LEN)
maske = np.zeros((LEN, LEN), dtype=bool)
pos_x, pos_y = lat.get_from_mask(maske)
si = SpinImage(pos_x, pos_y)
fig, axs = plt.subplots(1, 4)
kernel = np.ones((20, 20))
array = signal.convolve2d(si.img, kernel, boundary='symm', mode='same')
axs[0].imshow(array)
[(x, y), (x1, y1)] = lat.get_both()
si = SpinImage_Two_Phase(x, y, x1, y1)
si.apply_mask(maske)
axs[1].imshow(si.img)
tmp = si.img
maske = np.invert(maske)
si.apply_mask(maske)
axs[2].imshow(si.img)
axs[3].imshow(si.img+tmp)
if __name__ == "__main__":
np.random.seed(1234)
# test_me()
# test_square()
test_mixed()
plt.show()
#random(1234)
# for i in np.random.randint(0, 10000, 1):
# random(i)
# ising(i)
# plt.show()

5
2d_fourie/plot.py
View File

@ -1,5 +0,0 @@
import matplotlib.pyplot as plt
import numpy as np
if __name__ == "__main__":
pass

270
2d_fourie/spin_image.py
View File

@ -1,270 +0,0 @@
import numpy as np
import scipy
import scipy.fftpack as sfft
import scipy.signal
import tqdm
from cache import timeit
class SpinImage_Two_Phase:
resolution = 0.05
offset = 40
def make_list(self, x_pos, y_pos, x_inds, y_inds):
sigma = .1
x_ind = np.arange(0, self.length_x, self.resolution)
y_ind = np.arange(0, self.length_y, self.resolution)
X, Y = np.meshgrid(x_ind, y_ind, indexing="ij")
out_list = []
for x, y, x_ind, y_ind in zip(
x_pos.flatten(),
y_pos.flatten(),
x_inds.flatten(),
y_inds.flatten(),
):
xl, yl, xu, yu = self.clean_bounds(
x_ind - self.offset,
y_ind - self.offset,
x_ind + self.offset,
y_ind + self.offset,
X.shape,
)
out_list.append(np.exp(-0.5 * ((X[xl:xu, yl:yu] - x) ** 2 +
(Y[xl:xu, yl:yu] - y) ** 2) / sigma**2))
#out_list.append(np.ones_like(X[xl:xu, yl:yu]))
out_list = np.array(out_list)
return out_list
@timeit
def __init__(self, x_pos_low, y_pos_low, x_pos_high, y_pos_high):
assert x_pos_low.shape == y_pos_low.shape
assert x_pos_high.shape == y_pos_high.shape
assert x_pos_low.shape == x_pos_high.shape
offset_shift = self.offset * self.resolution
zero_shift_x = np.minimum(np.min(x_pos_low), np.min(x_pos_high))
zero_shift_y = np.minimum(np.min(y_pos_low), np.min(y_pos_high))
x_pos_low = x_pos_low - zero_shift_x + offset_shift
y_pos_low = y_pos_low - zero_shift_y + offset_shift
x_pos_high = x_pos_high - zero_shift_x + offset_shift
y_pos_high = y_pos_high - zero_shift_y + offset_shift
max_len_x = np.maximum(np.max(x_pos_low), np.max(x_pos_high))
max_len_y = np.maximum(np.max(y_pos_low), np.max(y_pos_high))
self.length_x = max_len_x + (self.offset + 1) * self.resolution
self.length_y = max_len_y + (self.offset + 1) * self.resolution
self.x_index_low, self.y_index_low = self._stuff(x_pos_low, y_pos_low)
self.x_index_high, self.y_index_high = self._stuff(
x_pos_high, y_pos_high)
self.buffer_low = self.make_list(
x_pos_low, y_pos_low, self.x_index_low, self.y_index_low)
self.buffer_high = self.make_list(
x_pos_high, y_pos_high, self.x_index_high, self.y_index_high)
def clean_bounds(self, xl, yl, xu, yu, shape):
if xl < 0:
xl = 0
if yl < 0:
yl = 0
if xu > shape[0]:
xu = shape[0]
if yu > shape[1]:
yu = shape[1]
return xl, yl, xu, yu
def _stuff(self, pos_x, pos_y):
x_ind = np.arange(0, self.length_x, self.resolution)
y_ind = np.arange(0, self.length_y, self.resolution)
xind = np.searchsorted(x_ind, pos_x).astype(int)
yind = np.searchsorted(y_ind, pos_y).astype(int)
return xind, yind
def _apply_mask(self, x, y, buffer, mask):
for x, y, dat in zip(
x.flatten()[mask],
y.flatten()[mask],
buffer[mask],
):
xl, yl, xu, yu = self.clean_bounds(
x - self.offset,
y - self.offset,
x + self.offset,
y + self.offset,
self.img.shape,
)
self.img[x - self.offset: x + self.offset,
y - self.offset: y + self.offset] += dat
@timeit
def apply_mask(self, maske):
mask = np.empty_like(self.x_index_high, dtype=bool)
print(maske)
mask[0::2, 0::2] = maske
mask[1::2, 0::2] = maske
mask[0::2, 1::2] = maske
mask[1::2, 1::2] = maske
mask = mask.flatten()
self.img = np.zeros(
(int(self.length_x/self.resolution), int(self.length_y/self.resolution)))
self._apply_mask(self.x_index_low, self.y_index_low,
self.buffer_low, mask)
mask = np.invert(mask)
self._apply_mask(self.x_index_high, self.y_index_high,
self.buffer_high, mask)
@timeit
def fft(self):
Z_fft = sfft.fft2(self.img)
Z_shift = sfft.fftshift(Z_fft)
fft_freqx = sfft.fftfreq(self.img.shape[0], self.resolution)
fft_freqy = sfft.fftfreq(self.img.shape[1], self.resolution)
fft_freqx_clean = sfft.fftshift(fft_freqx)
fft_freqy_clean = sfft.fftshift(fft_freqy)
return fft_freqx_clean, fft_freqy_clean, np.abs(Z_shift) ** 2
@timeit
def plot(self, ax, scale=None):
if scale is None:
ax.imshow(self.img)
else:
quad = np.ones((int(scale / self.resolution),
int(scale / self.resolution)))
img = scipy.signal.convolve2d(self.img, quad)
ax.imshow(img)
@timeit
def pad_it_square(self, additional_pad=0, size=None):
h = self.img.shape[0]
w = self.img.shape[1]
xx = np.maximum(h, w) + 2 * additional_pad + 1
if size is not None:
xx = np.maximum(xx, size)
yy = xx
a = (xx - h) // 2
aa = xx - a - h
b = (yy - w) // 2
bb = yy - b - w
self.img = np.pad(self.img, pad_width=(
(a, aa), (b, bb)), mode="constant")
def gaussian(self, sigma):
x = np.linspace(0, self.length_x, self.img.shape[0])
y = np.linspace(0, self.length_y, self.img.shape[1])
X, Y = np.meshgrid(x, y)
X = X - (self.length_x / 2.)
Y = Y - (self.length_y / 2.)
z = 1 / (2 * np.pi * sigma * sigma) * \
np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2)))
self.img = np.multiply(self.img, z.T)
class SpinImage:
resolution = 0.05
def __init__(self, x_pos, y_pos):
x_pos = x_pos - np.min(x_pos)
y_pos = y_pos - np.min(y_pos)
self.length_x = np.max(x_pos) + self.resolution
self.length_y = np.max(y_pos) + self.resolution
self.img = self.image_from_pos(x_pos, y_pos)
def image_from_pos(self, pos_x, pos_y):
x_ind = np.arange(0, self.length_x, self.resolution) # angstrom
y_ind = np.arange(0, self.length_y, self.resolution) # angstrom
img = np.zeros((x_ind.size, y_ind.size))
xind = np.searchsorted(x_ind, pos_x)
yind = np.searchsorted(y_ind, pos_y)
img[xind, yind] = 1
return img
def image_from_pos_with_gauss(self, pos_x, pos_y, sigma=1):
x_ind = np.arange(0, self.length_x, self.resolution) # angstrom
y_ind = np.arange(0, self.length_y, self.resolution) # angstrom
img = np.zeros((x_ind.size, y_ind.size))
X, Y = np.meshgrid(y_ind, x_ind)
for px, py in tqdm.tqdm(zip(pos_x.flatten(), pos_y.flatten())):
img += np.exp(-0.5 * ((X - px) ** 2 + (Y - py) ** 2) / sigma**2)
return img
def fft(self):
Z_fft = sfft.fft2(self.img)
Z_shift = sfft.fftshift(Z_fft)
fft_freqx = sfft.fftfreq(self.img.shape[0], self.resolution)
fft_freqy = sfft.fftfreq(self.img.shape[1], self.resolution)
fft_freqx_clean = sfft.fftshift(fft_freqx)
fft_freqy_clean = sfft.fftshift(fft_freqy)
return fft_freqx_clean, fft_freqy_clean, np.abs(Z_shift) ** 2
def pad_it_square(self, additional_pad=0, size=None):
h = self.img.shape[0]
w = self.img.shape[1]
xx = np.maximum(h, w) + 2 * additional_pad + 1
if size is not None:
xx = np.maximum(xx, size)
yy = xx
self.length_x = xx * self.resolution
self.length_y = yy * self.resolution
a = (xx - h) // 2
aa = xx - a - h
b = (yy - w) // 2
bb = yy - b - w
self.img = np.pad(self.img, pad_width=(
(a, aa), (b, bb)), mode="constant")
def pad_it(self, x_size, y_size):
h = self.img.shape[0]
w = self.img.shape[1]
xx = x_size
yy = y_size
self.length_x = xx * self.resolution
self.length_y = yy * self.resolution
a = (xx - h) // 2
aa = xx - a - h
b = (yy - w) // 2
bb = yy - b - w
self.img = np.pad(self.img, pad_width=(
(a, aa), (b, bb)), mode="constant")
def percentage_gaussian(self, mask, sigma):
x = np.linspace(-self.length_x / 2, self.length_x / 2, mask.shape[0])
y = np.linspace(-self.length_y / 2, self.length_y / 2, mask.shape[1])
X, Y = np.meshgrid(x, y)
z = 1 / (2 * np.pi * sigma * sigma) * \
np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2)))
return np.multiply(mask, z.T)
def gaussian(self, sigma):
x = np.arange(-self.length_x / 2, self.length_x / 2, self.resolution)
y = np.arange(-self.length_y / 2, self.length_y / 2, self.resolution)
X, Y = np.meshgrid(x, y)
z = 1 / (2 * np.pi * sigma * sigma) * \
np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2)))
self.img = np.multiply(self.img, z.T)
def plot(self, ax, scale=None):
if scale is None:
ax.imshow(self.img)
else:
quad = np.ones((int(scale / self.resolution),
int(scale / self.resolution)))
img = scipy.signal.convolve2d(self.img, quad)
ax.imshow(img)
def blur(self, sigma):
self.img = scipy.ndimage.gaussian_filter(self.img, sigma)

46
2d_fourie/test_voronoi.py
View File

@ -1,46 +0,0 @@
import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import Voronoi, voronoi_plot_2d
import cv2
class Voronoi_Evaulator:
def __init__(self, points, eval_points):
self.eval_points = eval_points
self.vor = Voronoi(points)
def extract(self, Z):
for debug_num, ev_points in zip([-10, -5], self.eval_points):
region_mask = self.vor.point_region[ev_points]
print(region_mask)
for i in np.array(self.vor.regions)[region_mask]:
if -1 in i:
print("Containse outside points")
continue
if len(i) == 0:
print("Containse outside points")
continue
print(i)
pts = self.vor.vertices[i]
pts = (pts * 100).astype(np.int32)
print(pts)
mask = np.zeros((Z.shape[0], Z.shape[1]))
cv2.fillConvexPoly(mask, pts, 1)
mask = mask > 0 # To convert to Boolean
Z[mask] = debug_num
return Z
if __name__ == "__main__":
np.random.seed(20)
points = (np.random.rand(100, 2)-0.1) * 2
voro = Voronoi_Evaulator(points, [[1, 4, 5], [2, 3, 6]])
x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 200)
X, Y = np.meshgrid(x, y)
Z = X*2 + Y
Z = voro.extract(Z)
plt.imshow(Z)
plt.show()

11
rust/Cargo.toml
View File

@ -1,11 +0,0 @@
[package]
name = "rust"
version = "0.1.0"
edition = "2021"
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
[dependencies]
ndarray = "0.15.6"
ndarray-npy = "0.8.1"
plotters = "0.3.4"

9
rust/plot.py
View File

@ -1,9 +0,0 @@
import numpy as np
import matplotlib.pyplot as plt
file = np.load("./rutile.npz")
plt.scatter(file["x"], file["y"])
file = np.load("./mono.npz")
plt.scatter(file["x"], file["y"])
plt.show()

19
rust/src/main.rs
View File

@ -1,19 +0,0 @@
mod vo2;
use vo2::get_mono;
use crate::vo2::get_rutile;
use ndarray_npy::NpzWriter;
use std::fs::File;
fn main() {
let (mono_x, mono_y) = get_mono(20, 20);
let (rutile_x, rutile_y) = get_rutile(20, 20);
let mut npz = NpzWriter::new(File::create("mono.npz").unwrap());
npz.add_array("x", &mono_x);
npz.add_array("y", &mono_y);
npz.finish().unwrap();
let mut npz = NpzWriter::new(File::create("rutile.npz").unwrap());
npz.add_array("x", &rutile_x);
npz.add_array("y", &rutile_y);
npz.finish().unwrap();
}

143
rust/src/vo2.rs
View File

@ -1,143 +0,0 @@
use ndarray::{Array, Array2};
const BASE_A_M: f64 = 5.75;
const BASE_B_M: f64 = 4.5;
const BASE_C_M: f64 = 5.38;
const BASE_C_R: f64 = 2.856;
const BASE_B_R: f64 = 4.554;
const BASE_A_R: f64 = BASE_B_R;
const ALPHA_M: f64 = 122.64; // degree
fn deg_2_rad(angle_degree: f64) -> f64 {
return std::f64::consts::PI / 180. * angle_degree;
}
fn mono_2_rutile(c_m: f64, a_m: f64) -> (f64, f64) {
let a_r = deg_2_rad(ALPHA_M - 90.).cos() * c_m * BASE_C_M;
let c_r = (a_m) * BASE_A_M + deg_2_rad(ALPHA_M - 90.).sin() * c_m * BASE_C_M;
return (a_r, c_r);
}
pub fn get_rutile(len_x: usize, len_y: usize) -> (Array2<f64>, Array2<f64>) {
let mut rutile_x_pos: Array2<f64> = Array::zeros((len_x * 2, len_y * 2));
let mut rutile_y_pos: Array2<f64> = Array::zeros((len_x * 2, len_y * 2));
for x in 0..2 * len_x {
for y in 0..2 * len_y {
let tmp_x: f64 = if y % 4 >= 2 { x as f64 - 1. } else { x as f64 };
rutile_x_pos[[x, y]] = tmp_x as f64 * BASE_C_R + (y % 4) as f64 * 0.5 * BASE_C_R;
rutile_y_pos[[x, y]] = y as f64 * 0.5 * BASE_A_R;
}
}
return (rutile_x_pos, rutile_y_pos);
}
pub fn get_mono(len_x: usize, len_y: usize) -> (Array2<f64>, Array2<f64>) {
const OFFSET_A_M: f64 = 0.25 - 0.23947;
const OFFSET_C_M: f64 = 0.02646;
let (offset_a_r, offset_c_r) = mono_2_rutile(OFFSET_C_M, OFFSET_A_M);
let mut mono_x_pos: Array2<f64> = Array::zeros((len_x * 2, len_y * 2));
let mut mono_y_pos: Array2<f64> = Array::zeros((len_x * 2, len_y * 2));
for x in 0..2 * len_x {
for y in 0..2 * len_y {
let tmp_x: i64 = if y % 4 >= 2 { x as i64 - 1 } else { x as i64 };
let offset_sign = if tmp_x % 2 == 0 { -1. } else { 1. };
mono_x_pos[[x, y]] = offset_sign * 2. * offset_a_r
+ tmp_x as f64 * BASE_C_R
+ (y % 4) as f64 * 0.5 * BASE_C_R;
mono_y_pos[[x, y]] = offset_sign * 2. * offset_c_r + 0.5 * y as f64 * BASE_A_R;
}
}
return (mono_x_pos, mono_y_pos);
} /*
def _mono_2_rutile(self, c_m, a_m):
return a_r, c_r
def _get_rutile(self, X, Y):
self.atom_x_rut = X * self.base_c_r + \
np.mod(Y, 4) * 0.5 * self.base_c_r
self.atom_y_rut = Y * 0.5 * self.base_a_r
def _get_mono(self, X, Y):
offset_a_m = 0.25 - 0.23947
offset_ndc_m = 0.02646
offset_a_r, offset_c_r = self._mono_2_rutile(offset_c_m, offset_a_m)
res = 0.05
offset_a_r = res * int(offset_a_r/res)
offset_c_r = res * int(offset_c_r/res)
print(offset_a_r, offset_c_r)
self.atom_x_mono = offset_a_r + X * \
self.base_c_r + np.mod(Y, 4) * 0.5 * self.base_c_r
self.atom_x_mono[np.mod(X, 2) == 0] -= 2 * offset_a_r
self.atom_y_mono = offset_c_r + 0.5 * Y * self.base_a_r
self.atom_y_mono[np.mod(X, 2) == 0] -= 2 * offset_c_r
def _generate_vec(self, x_len: int, y_len: int):
x = np.arange(x_len)
y = np.arange(y_len)
X, Y = np.meshgrid(x, y)
X[np.mod(Y, 4) == 3] = X[np.mod(Y, 4) == 3] - 1
X[np.mod(Y, 4) == 2] = X[np.mod(Y, 4) == 2] - 1
assert np.mod(x.size, 2) == 0
assert np.mod(y.size, 2) == 0
return X, Y
@timeit
def __init__(self, x_len: int, y_len: int):
X, Y = self._generate_vec(x_len * 2, y_len * 2)
self._get_mono(X, Y)
self._get_rutile(X, Y)
def get_from_mask(self, maske: np.array):
inplace_pos_x = np.zeros_like(self.atom_x_mono)
inplace_pos_y = np.zeros_like(self.atom_x_mono)
mask = np.empty_like(self.atom_x_mono, dtype=bool)
mask[0::2, 0::2] = maske
mask[1::2, 0::2] = maske
mask[0::2, 1::2] = maske
mask[1::2, 1::2] = maske
inplace_pos_x[mask] = self.atom_x_rut[mask]
inplace_pos_y[mask] = self.atom_y_rut[mask]
mask = np.invert(mask)
inplace_pos_x[mask] = self.atom_x_mono[mask]
inplace_pos_y[mask] = self.atom_y_mono[mask]
return inplace_pos_x, inplace_pos_y
def reci_rutile(self):
num = 20
#num = 2
x = np.arange(-num, num + 1)
y = np.arange(-num, num + 1)
X, Y = np.meshgrid(x, y)
return (X * 0.22 + Y * 0.44).flatten(), (X * 0.349).flatten()
def reci_mono(self):
x, y = self.reci_rutile()
return x + 0.1083, y + 0.1719
def reci_mono_2(self):
x, y = self.reci_rutile()
return x - 0.1083, y + 0.1719
def reci(self):
cutoff = 5.
x, y = self.reci_rutile()
mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff)
p1 = (x[mask], y[mask])
x, y = self.reci_mono()
mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff)
p2 = (x[mask], y[mask])
x, y = self.reci_mono_2()
mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff)
p3 = (x[mask], y[mask])
return [p1, p2, p3]
*/