Skip to content

Commit

Permalink
copied files from dev branch to here
Browse files Browse the repository at this point in the history
  • Loading branch information
rmdocherty committed Jul 30, 2024
1 parent 7c6bfda commit a7cf40c
Show file tree
Hide file tree
Showing 11 changed files with 698 additions and 787 deletions.
104 changes: 50 additions & 54 deletions paper_figures/fig1.py
100755 → 100644
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,37 +4,34 @@
import json
import numpy as np
import torch
from representativity.old import util
from representativity import util
from torch.nn.functional import interpolate


def sample_error(img, ls, vf=0.5):
err = []
dims = len(img.shape) - 1
for l in ls:
if dims == 1:
vfs = torch.mean(img[:, : l * l], dim=(1))
elif dims == 2:
vfs = torch.mean(img[:, :l*l], dim=(1))
elif dims == 2:
vfs = torch.mean(img[:, :l, :l], dim=(1, 2))
else: # 3D
new_l = int(l ** (2 / 3))
new_l = int(l**(2/3))
vfs = torch.mean(img[:, :new_l, :new_l, :new_l], dim=(1, 2, 3))
std = torch.std(vfs.cpu())
err.append(100 * ((1.96 * std) / vf)) # percentage of error from 0.5
err.append(100*((1.96*std)/vf)) # percentage of error from 0.5
return err


errs = []
berns = []


def generate_microlib_data():
mode = "2D"
mode = '2D'
# Dataset path and list of subfolders
# with open("micro_names.json", "r") as fp: # TODO change this later
# micro_names = json.load(fp)
plotting = [f"microstructure{f}" for f in [228, 235, 205, 177]]
projects = [f"/home/amir/microlibDataset/{p}/{p}" for p in plotting]
# micro_names = json.load(fp)
plotting = [f'microstructure{f}' for f in [228, 235, 205, 177]]
projects = [f'/home/amir/microlibDataset/{p}/{p}' for p in plotting]
# Load generator network
netG = util.load_generator(projects[0])
# Edge lengths to test
Expand All @@ -46,74 +43,73 @@ def generate_microlib_data():

num_projects = len(projects)
projects = projects[:num_projects]
reps = 1000 if mode == "2D" else 400
reps = 1000 if mode=='2D' else 400
for j, proj in enumerate(projects):
# Make an image of micro
netG.load_state_dict(torch.load(proj + "_Gen.pt"))
img = util.generate_image(netG, lf=l, reps=reps, threed=mode == "3D")
img = util.generate_image(netG,lf=l,reps=reps, threed=mode=='3D')
print(img.size())
if img.any():
microstructure_name = proj.split("/")[-1]
vf = torch.mean(img).cpu().item()
img = [img]
print(f"starting tpc")
img_dims = [np.array((l,) * (len(img.shape) - 1)) for l in edge_lengths]
print(f'starting tpc')
img_dims = [np.array((l,)*(len(img.shape)-1)) for l in edge_lengths]
err_exp = util.real_image_stats(img, edge_lengths, vf)
return data


# with open("data_gen.json", "r") as fp:
# data = json.load(fp)['generated_data']
# data = json.load(fp)['generated_data']

n = "microstructure205"
img3 = tifffile.imread(f"/home/amir/microlibDataset/{n}/{n}.tif")
n = 'microstructure205'
img3 = tifffile.imread(f'/home/amir/microlibDataset/{n}/{n}.tif')
d = generate_microlib_data()


ls = torch.arange(300, 800, 20)
ls = torch.arange(300,800, 20)
img_dims = [np.array((l, l)) for l in ls]
ns = util.ns_from_dims(img_dims, 1)
img1 = torch.rand((1000, 1, ls[-1], ls[-1]), device=torch.device("cuda:0")).float()
img1[img1 <= d["vf"]] = 0
img1[img1 > d["vf"]] = 1
img1 = 1 - img1
errs.append(sample_error(img1[:, 0], ls, d["vf"]))
berns.append(util.bernouli(d["vf"], ns))

img2 = interpolate(
img1[:, :, : ls[-1] // 2, : ls[-1] // 2], scale_factor=(2, 2), mode="nearest"
)
errs.append(sample_error(img2[:, 0], ls, d["vf"]))
img1 = torch.rand((1000, 1, ls[-1],ls[-1]), device = torch.device('cuda:0')).float()
img1[img1<=d['vf']] = 0
img1[img1>d['vf']] = 1
img1 = 1-img1
errs.append(sample_error(img1[:,0], ls, d['vf']))
berns.append(util.bernouli(d['vf'], ns))

img2 = interpolate(img1[:,:,:ls[-1]//2, :ls[-1]//2], scale_factor=(2,2), mode='nearest')
errs.append(sample_error(img2[:,0], ls, d['vf']))
ns = util.ns_from_dims(img_dims, 2)
berns.append(util.bernouli(d["vf"], ns))
berns.append(util.bernouli(d['vf'], ns))

errs.append(d[f"err_exp_vf"][::4])
ns = util.ns_from_dims(img_dims, d["fac_vf"])
berns.append(util.bernouli(d["vf"], ns))
errs.append(d[f'err_exp_vf'][::4])
ns = util.ns_from_dims(img_dims, d['fac_vf'])
berns.append(util.bernouli(d['vf'], ns))

fig, axs = plt.subplots(2, 2)
fig.set_size_inches(8, 8)
l = 150
axs[0, 0].imshow(img1[0, 0, :l, :l].cpu(), cmap="gray")
axs[0, 1].imshow(img2[0, 0, :l, :l].cpu(), cmap="gray")
axs[1, 0].imshow(img3[0, :150, :150], cmap="gray")
fig, axs = plt.subplots(2,2)
fig.set_size_inches(8,8)
l=150
axs[0,0].imshow(img1[0,0, :l, :l].cpu(), cmap='gray')
axs[0,1].imshow(img2[0,0, :l, :l].cpu(), cmap='gray')
axs[1,0].imshow(img3[0, :150, :150], cmap='gray')

cs = ["r", "b", "g"]
labs = ["a) random 1x1", "b) random 2x2", "c) micrograph 205"]
cs = ['r', 'b', 'g']
labs = ['a) random 1x1', 'b) random 2x2', 'c) micrograph 205']
for l, err, bern, c in zip(labs, errs, berns, cs):
axs[1, 1].scatter(ls, err, c=c, s=8, marker="x")
axs[1, 1].plot(ls, bern, lw=1, c=c, label=l)
axs[1, 1].legend()
axs[0, 0].set_title("a) Random 1x1 (vf = 0.845, $a_2$ = 1)")
axs[0, 1].set_title("b) Random 2x2 (vf = 0.845, $a_2$ = 2)")
axs[1, 0].set_title("c) Micro. 205 (vf = 0.845, $a_2$ = 7.526)")
axs[1, 1].set_title("d) Experimental integral range fit")
axs[1, 1].set_xlabel("Image length size [pixels]")
axs[1, 1].set_ylabel("Volume fraction percentage error [%]")
axs[1, 1].set_ylim(bottom=0)
axs[1,1].scatter(ls, err, c=c, s=8,marker = 'x')
axs[1,1].plot(ls, bern, lw=1, c=c, label=l)
axs[1,1].legend()
axs[0,0].set_title('a) Random 1x1 (vf = 0.845, $a_2$ = 1)')
axs[0,1].set_title('b) Random 2x2 (vf = 0.845, $a_2$ = 2)')
axs[1,0].set_title('c) Micro. 205 (vf = 0.845, $a_2$ = 7.526)')
axs[1,1].set_title('d) Experimental integral range fit')
axs[1,1].set_xlabel('Image length size [pixels]')
axs[1,1].set_ylabel('Volume fraction percentage error [%]')
axs[1,1].set_ylim(bottom=0)
plt.tight_layout()
for ax in axs.ravel()[:3]:
ax.set_xticks([])
ax.set_yticks([])

plt.savefig("fig1.pdf", format="pdf")
plt.savefig('fig1.pdf', format='pdf')

207 changes: 101 additions & 106 deletions paper_figures/fig2.py
100755 → 100644
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,128 +3,118 @@
import tifffile
import json
import numpy as np

# import kneed
from scipy.optimize import curve_fit
from scipy import stats
import torch
from representativity.old import util
from representativity import util, microlib_statistics
from torch.nn.functional import interpolate
from mpl_toolkits.axes_grid1 import make_axes_locatable

with open("data_gen2.json", "r") as fp:
data = json.load(fp)["generated_data"]
# with open("data_gen2.json", "r") as fp:
# data = json.load(fp)["generated_data"]

l = len(list(data.keys()))
# l=len(list(data.keys()))
# l=3
c = [(0, 0, 0), (0.5, 0.5, 0.5)]
c=[(0,0,0), (0.5,0.5,0.5)]
# plotting = [f'microstructure{f}' for f in [235, 209,205,177]]
plotting = [f"microstructure{f}" for f in [235, 228, 205, 177]]
plotting_ims = [f'microstructure{f}' for f in [235, 228, 205,177]]

# plotting = [k for k in data.keys()]
l = len(plotting)
l = len(plotting_ims)
fig, axs = plt.subplots(l, 3)
fig.set_size_inches(12, l * 4)
preds = [[], []]
irs = [[], []]
sas = []
i = 0
for n in list(data.keys()):
if n not in plotting:
fig.set_size_inches(12, l*3.5)
preds = [[],[]]
irs = [[],[]]
colors = {"cls": 'tab:orange', "stds": 'tab:green'}

all_data, micros, netG, v_names, run_v_names = microlib_statistics.json_preprocessing()
lens_for_fit = list(range(500, 1000, 20)) # lengths of images for fitting L_characteristic
plotting_ims = [micro for micro in micros if micro.split('/')[-1] in plotting_ims]
# run the statistical analysis on the microlib dataset
for i, p in enumerate(plotting_ims):

try:
netG.load_state_dict(torch.load(p + "_Gen.pt"))
except: # if the image is greayscale it's excepting because there's only 1 channel
continue
img = tifffile.imread(f"/home/amir/microlibDataset/{n}/{n}.tif")
d = data[n]
d1 = data[n]

csets = [["black", "black"], ["gray", "gray"]]
for j, met in enumerate(["vf", "sa"]):
cs = csets[j]
img_dims = [np.array([int(im_len)] * 2) for im_len in d["ls"]]
ns = util.ns_from_dims(img_dims, d[f"ir_{met}"])
berns_vf = util.bernouli(d[f"{met}"], ns)
axs[i, 1].scatter(
d["ls"],
d[f"err_exp_{met}"],
c=cs[0],
s=8,
marker="x",
label=f"{met} errors from sampling",
)
axs[i, 1].plot(d["ls"], berns_vf, c=cs[0], label=f"{met} errors from fitted IR")
# axs[i, 1].plot(d[f'ls'], d[f'err_model_{met}'], c=cs[0], label = f'{met} errors from bernouli')
y = d[f"tpc_{met}"][
0
] # TODO write that in sa tpc, only the first direction is shown, or do something else, maybe normalise the tpc? We can do sum, because that's how we calculate the ir!
x = d[f"tpc_{met}_dist"]
y = np.array(y)

# TODO erase this afterwards:
if met == "vf":
ir = np.round(d[f"ir_vf"], 1)
axs[i, 2].plot(x, y, c=cs[1], label=f"Volume fraction 2PC")
axs[i, 2].axhline(d["vf"] ** 2, linestyle="dashed", label="$p^2$")
axs[i, 2].plot(
[0, ir],
[d["vf"] ** 2 - 0.02, d["vf"] ** 2 - 0.02],
c="green",
linewidth=3,
label=r"$\tilde{a}_2$",
)
ticks = [0, int(ir), 20, 40, 60, 80, 100]
ticklabels = map(str, ticks)
axs[i, 2].set_xticks(ticks)
axs[i, 2].set_xticklabels(ticklabels)
axs[i, 2].fill_between(
x, d["vf"] ** 2, y, alpha=0.5, label=f"Integrated part"
)
axs[i, 2].legend()

# axs[i, 2].scatter(x[knee], y_data[knee]/y.max(), c =cs[1], marker = 'x', label=f'{met} ir from tpc', s=100)
ir = d[f"ir_{met}"]
pred_ir = util.tpc_to_ir(d[f"tpc_{met}_dist"], d[f"tpc_{met}"])
pred_ir = pred_ir * 1.61
# axs[i, 2].scatter(x[round(pred_ir)], y[round(pred_ir)], c =cs[1], marker = 'x', label=f'{met} predicted tpc IR', s=100)
# axs[i, 2].scatter(x[round(ir)], y[round(ir)], facecolors='none', edgecolors = cs[1], label=f'{met} fitted IR', s=100)

irs[j].append(ir)
if i == 0:
axs[i, 1].legend()
axs[i, 2].legend()
axs[i, 1].set_xlabel("Image length size [pixels]")
axs[i, 1].set_ylabel("Volume fraction percentage error [%]")
axs[i, 2].set_ylabel("2-point correlation function")
ir = np.round(d[f"ir_vf"], 2)
sas.append(d["sa"])
im = img[0] * 255
si_size, nirs = 160, 5
sicrop = int(ir * nirs)
print(ir, sicrop)
subim = torch.tensor(im[-sicrop:, -sicrop:]).unsqueeze(0).unsqueeze(0).float()
subim = interpolate(subim, size=(si_size, si_size), mode="nearest")[0, 0]
subim = np.stack([subim] * 3, axis=-1)

subim[:5, :, :] = 125
subim[:, :5, :] = 125
# subim[5:20, 5:50, :]
subim[10:15, 10 : 10 + si_size // nirs, :] = 0
subim[10:15, 10 : 10 + si_size // nirs, 1:-1] = 125

im = np.stack([im] * 3, axis=-1)
im[-si_size:, -si_size:] = subim
axs[i, 0].imshow(im)
axs[i, 0].set_xticks([])
axs[i, 0].set_yticks([])
axs[i, 0].set_ylabel(f"M{n[1:]}")
axs[i, 0].set_xlabel(
f"Volume fraction "
+ r"$\tilde{a}_2$: "
+ f"{ir} Inset mag: x{np.round(si_size/sicrop, 2)}"
)

i += 1
imsize = 1600
lf = imsize//32 + 2 # the size of G's input
many_images = util.generate_image(netG, lf=lf, threed=False, reps=10)
pf = many_images.mean().cpu().numpy()
small_imsize = 512
img = many_images[0].detach().cpu().numpy()[:small_imsize, :small_imsize]

csets = [['black', 'black'], ['gray', 'gray']]
conf = 0.95
errs = util.real_image_stats(many_images, lens_for_fit, pf, conf=conf)
sizes_for_fit = [[lf,lf] for lf in lens_for_fit]
real_cls = util.fit_cls(errs, sizes_for_fit, pf)
stds = errs / stats.norm.interval(conf)[1] * pf / 100
std_fit = (real_cls**2/(np.array(lens_for_fit)**2)*pf*(1-pf))**0.5
# print(stds)
vars = stds**2
# from variations to L_characteristic using image size and phase fraction
clss = (np.array(lens_for_fit)**2*vars/pf/(1-pf))**0.5
print(clss)

axs_twin = axs[i, 1].twinx()
stds_scatter = axs_twin.scatter(lens_for_fit, stds, s=8, color=colors["stds"], label = f'Standard deviations')
twin_fit = axs_twin.plot(lens_for_fit, std_fit, color=colors["stds"], label = f'Standard deviations fit')
axs_twin.tick_params(axis='y', labelcolor=colors["stds"])
axs_twin.set_ylim(0, 0.025)
cls_scatter = axs[i, 1].scatter(lens_for_fit, clss, s=8, color=colors["cls"], label = f'Characteristic length scales')
cls_fit = axs[i, 1].hlines(real_cls, lens_for_fit[0], lens_for_fit[-1], color=colors["cls"], label = f'Characteristic length scales fit')
axs[i, 1].tick_params(axis='y', labelcolor=colors["cls"])
axs[i, 1].set_ylim(0, 28)

dims = len(img.shape)
# print(f'starting tpc radial')
tpc = util.tpc_radial(img, threed=False)
cls = util.tpc_to_cls(tpc, img, img.shape)

contour = axs[i, 2].contourf(tpc, cmap='plasma', levels = 200)
divider = make_axes_locatable(axs[i, 2])
cax = divider.append_axes('right', size='5%', pad=0.05)
fig.colorbar(contour, cax=cax, orientation='vertical')

if i == 0:
plots = [stds_scatter, twin_fit[0], cls_scatter, cls_fit]
label_plots = [plot.get_label() for plot in plots]
axs[i,1].legend(plots, label_plots)

# axs[i,2].legend()
# axs[i,1].set_xlabel('Image length size [pixels]')
# axs[i,1].set_ylabel('Volume fraction percentage error [%]')
# axs[i,2].set_ylabel('2-point correlation function')
# ir = np.round(d[f'ir_vf'], 2)
# im = img[0]*255
# si_size, nirs = 160, 5
# sicrop = int(ir*nirs)
# print(ir, sicrop)
# subim=torch.tensor(im[-sicrop:,-sicrop:]).unsqueeze(0).unsqueeze(0).float()
# subim = interpolate(subim, size=(si_size,si_size), mode='nearest')[0,0]
# subim = np.stack([subim]*3, axis=-1)

# subim[:5,:,:] = 125
# subim[:,:5,:] = 125
# # subim[5:20, 5:50, :]
# subim[10:15, 10:10+si_size//nirs, :] = 0
# subim[10:15, 10:10+si_size//nirs, 1:-1] = 125

# im = np.stack([im]*3, axis=-1)
# im[-si_size:,-si_size:] = subim
# axs[i, 0].imshow(im)
# axs[i, 0].set_xticks([])
# axs[i, 0].set_yticks([])
# axs[i, 0].set_ylabel(f'M{n[1:]}')
# axs[i, 0].set_xlabel(f'Volume fraction '+ r'$\tilde{a}_2$: '+ f'{ir} Inset mag: x{np.round(si_size/sicrop, 2)}')

i+=1


plt.tight_layout()
plt.savefig("fig2.pdf", format="pdf")
plt.savefig('fig2.pdf', format='pdf')

# fig, axs = plt.subplots(1,2)
# fig.set_size_inches(10,5)
Expand Down Expand Up @@ -154,3 +144,8 @@
# # res = 100*(np.mean((y-targ)**2/targ))**0.5

# print(res,res2)





Loading

0 comments on commit a7cf40c

Please sign in to comment.