Skip to content

Commit a7cf40c

Browse files
author
rmdocherty
committed
copied files from dev branch to here
1 parent 7c6bfda commit a7cf40c

File tree

11 files changed

+698
-787
lines changed

11 files changed

+698
-787
lines changed

‎paper_figures/fig1.py

100755100644
+50-54
Original file line numberDiff line numberDiff line change
@@ -4,37 +4,34 @@
44
import json
55
import numpy as np
66
import torch
7-
from representativity.old import util
7+
from representativity import util
88
from torch.nn.functional import interpolate
99

10-
1110
def sample_error(img, ls, vf=0.5):
1211
err = []
1312
dims = len(img.shape) - 1
1413
for l in ls:
1514
if dims == 1:
16-
vfs = torch.mean(img[:, : l * l], dim=(1))
17-
elif dims == 2:
15+
vfs = torch.mean(img[:, :l*l], dim=(1))
16+
elif dims == 2:
1817
vfs = torch.mean(img[:, :l, :l], dim=(1, 2))
1918
else: # 3D
20-
new_l = int(l ** (2 / 3))
19+
new_l = int(l**(2/3))
2120
vfs = torch.mean(img[:, :new_l, :new_l, :new_l], dim=(1, 2, 3))
2221
std = torch.std(vfs.cpu())
23-
err.append(100 * ((1.96 * std) / vf)) # percentage of error from 0.5
22+
err.append(100*((1.96*std)/vf)) # percentage of error from 0.5
2423
return err
2524

26-
2725
errs = []
2826
berns = []
2927

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

4744
num_projects = len(projects)
4845
projects = projects[:num_projects]
49-
reps = 1000 if mode == "2D" else 400
46+
reps = 1000 if mode=='2D' else 400
5047
for j, proj in enumerate(projects):
5148
# Make an image of micro
5249
netG.load_state_dict(torch.load(proj + "_Gen.pt"))
53-
img = util.generate_image(netG, lf=l, reps=reps, threed=mode == "3D")
50+
img = util.generate_image(netG,lf=l,reps=reps, threed=mode=='3D')
5451
print(img.size())
5552
if img.any():
5653
microstructure_name = proj.split("/")[-1]
5754
vf = torch.mean(img).cpu().item()
5855
img = [img]
59-
print(f"starting tpc")
60-
img_dims = [np.array((l,) * (len(img.shape) - 1)) for l in edge_lengths]
56+
print(f'starting tpc')
57+
img_dims = [np.array((l,)*(len(img.shape)-1)) for l in edge_lengths]
6158
err_exp = util.real_image_stats(img, edge_lengths, vf)
6259
return data
6360

6461

6562
# with open("data_gen.json", "r") as fp:
66-
# data = json.load(fp)['generated_data']
63+
# data = json.load(fp)['generated_data']
6764

68-
n = "microstructure205"
69-
img3 = tifffile.imread(f"/home/amir/microlibDataset/{n}/{n}.tif")
65+
n = 'microstructure205'
66+
img3 = tifffile.imread(f'/home/amir/microlibDataset/{n}/{n}.tif')
7067
d = generate_microlib_data()
7168

7269

73-
ls = torch.arange(300, 800, 20)
70+
ls = torch.arange(300,800, 20)
7471
img_dims = [np.array((l, l)) for l in ls]
7572
ns = util.ns_from_dims(img_dims, 1)
76-
img1 = torch.rand((1000, 1, ls[-1], ls[-1]), device=torch.device("cuda:0")).float()
77-
img1[img1 <= d["vf"]] = 0
78-
img1[img1 > d["vf"]] = 1
79-
img1 = 1 - img1
80-
errs.append(sample_error(img1[:, 0], ls, d["vf"]))
81-
berns.append(util.bernouli(d["vf"], ns))
82-
83-
img2 = interpolate(
84-
img1[:, :, : ls[-1] // 2, : ls[-1] // 2], scale_factor=(2, 2), mode="nearest"
85-
)
86-
errs.append(sample_error(img2[:, 0], ls, d["vf"]))
73+
img1 = torch.rand((1000, 1, ls[-1],ls[-1]), device = torch.device('cuda:0')).float()
74+
img1[img1<=d['vf']] = 0
75+
img1[img1>d['vf']] = 1
76+
img1 = 1-img1
77+
errs.append(sample_error(img1[:,0], ls, d['vf']))
78+
berns.append(util.bernouli(d['vf'], ns))
79+
80+
img2 = interpolate(img1[:,:,:ls[-1]//2, :ls[-1]//2], scale_factor=(2,2), mode='nearest')
81+
errs.append(sample_error(img2[:,0], ls, d['vf']))
8782
ns = util.ns_from_dims(img_dims, 2)
88-
berns.append(util.bernouli(d["vf"], ns))
83+
berns.append(util.bernouli(d['vf'], ns))
8984

90-
errs.append(d[f"err_exp_vf"][::4])
91-
ns = util.ns_from_dims(img_dims, d["fac_vf"])
92-
berns.append(util.bernouli(d["vf"], ns))
85+
errs.append(d[f'err_exp_vf'][::4])
86+
ns = util.ns_from_dims(img_dims, d['fac_vf'])
87+
berns.append(util.bernouli(d['vf'], ns))
9388

94-
fig, axs = plt.subplots(2, 2)
95-
fig.set_size_inches(8, 8)
96-
l = 150
97-
axs[0, 0].imshow(img1[0, 0, :l, :l].cpu(), cmap="gray")
98-
axs[0, 1].imshow(img2[0, 0, :l, :l].cpu(), cmap="gray")
99-
axs[1, 0].imshow(img3[0, :150, :150], cmap="gray")
89+
fig, axs = plt.subplots(2,2)
90+
fig.set_size_inches(8,8)
91+
l=150
92+
axs[0,0].imshow(img1[0,0, :l, :l].cpu(), cmap='gray')
93+
axs[0,1].imshow(img2[0,0, :l, :l].cpu(), cmap='gray')
94+
axs[1,0].imshow(img3[0, :150, :150], cmap='gray')
10095

101-
cs = ["r", "b", "g"]
102-
labs = ["a) random 1x1", "b) random 2x2", "c) micrograph 205"]
96+
cs = ['r', 'b', 'g']
97+
labs = ['a) random 1x1', 'b) random 2x2', 'c) micrograph 205']
10398
for l, err, bern, c in zip(labs, errs, berns, cs):
104-
axs[1, 1].scatter(ls, err, c=c, s=8, marker="x")
105-
axs[1, 1].plot(ls, bern, lw=1, c=c, label=l)
106-
axs[1, 1].legend()
107-
axs[0, 0].set_title("a) Random 1x1 (vf = 0.845, $a_2$ = 1)")
108-
axs[0, 1].set_title("b) Random 2x2 (vf = 0.845, $a_2$ = 2)")
109-
axs[1, 0].set_title("c) Micro. 205 (vf = 0.845, $a_2$ = 7.526)")
110-
axs[1, 1].set_title("d) Experimental integral range fit")
111-
axs[1, 1].set_xlabel("Image length size [pixels]")
112-
axs[1, 1].set_ylabel("Volume fraction percentage error [%]")
113-
axs[1, 1].set_ylim(bottom=0)
99+
axs[1,1].scatter(ls, err, c=c, s=8,marker = 'x')
100+
axs[1,1].plot(ls, bern, lw=1, c=c, label=l)
101+
axs[1,1].legend()
102+
axs[0,0].set_title('a) Random 1x1 (vf = 0.845, $a_2$ = 1)')
103+
axs[0,1].set_title('b) Random 2x2 (vf = 0.845, $a_2$ = 2)')
104+
axs[1,0].set_title('c) Micro. 205 (vf = 0.845, $a_2$ = 7.526)')
105+
axs[1,1].set_title('d) Experimental integral range fit')
106+
axs[1,1].set_xlabel('Image length size [pixels]')
107+
axs[1,1].set_ylabel('Volume fraction percentage error [%]')
108+
axs[1,1].set_ylim(bottom=0)
114109
plt.tight_layout()
115110
for ax in axs.ravel()[:3]:
116111
ax.set_xticks([])
117112
ax.set_yticks([])
118113

119-
plt.savefig("fig1.pdf", format="pdf")
114+
plt.savefig('fig1.pdf', format='pdf')
115+

‎paper_figures/fig2.py

100755100644
+101-106
Original file line numberDiff line numberDiff line change
@@ -3,128 +3,118 @@
33
import tifffile
44
import json
55
import numpy as np
6-
76
# import kneed
87
from scipy.optimize import curve_fit
98
from scipy import stats
109
import torch
11-
from representativity.old import util
10+
from representativity import util, microlib_statistics
1211
from torch.nn.functional import interpolate
12+
from mpl_toolkits.axes_grid1 import make_axes_locatable
1313

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

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

2323
# plotting = [k for k in data.keys()]
24-
l = len(plotting)
24+
l = len(plotting_ims)
2525
fig, axs = plt.subplots(l, 3)
26-
fig.set_size_inches(12, l * 4)
27-
preds = [[], []]
28-
irs = [[], []]
29-
sas = []
30-
i = 0
31-
for n in list(data.keys()):
32-
if n not in plotting:
26+
fig.set_size_inches(12, l*3.5)
27+
preds = [[],[]]
28+
irs = [[],[]]
29+
colors = {"cls": 'tab:orange', "stds": 'tab:green'}
30+
31+
all_data, micros, netG, v_names, run_v_names = microlib_statistics.json_preprocessing()
32+
lens_for_fit = list(range(500, 1000, 20)) # lengths of images for fitting L_characteristic
33+
plotting_ims = [micro for micro in micros if micro.split('/')[-1] in plotting_ims]
34+
# run the statistical analysis on the microlib dataset
35+
for i, p in enumerate(plotting_ims):
36+
37+
try:
38+
netG.load_state_dict(torch.load(p + "_Gen.pt"))
39+
except: # if the image is greayscale it's excepting because there's only 1 channel
3340
continue
34-
img = tifffile.imread(f"/home/amir/microlibDataset/{n}/{n}.tif")
35-
d = data[n]
36-
d1 = data[n]
37-
38-
csets = [["black", "black"], ["gray", "gray"]]
39-
for j, met in enumerate(["vf", "sa"]):
40-
cs = csets[j]
41-
img_dims = [np.array([int(im_len)] * 2) for im_len in d["ls"]]
42-
ns = util.ns_from_dims(img_dims, d[f"ir_{met}"])
43-
berns_vf = util.bernouli(d[f"{met}"], ns)
44-
axs[i, 1].scatter(
45-
d["ls"],
46-
d[f"err_exp_{met}"],
47-
c=cs[0],
48-
s=8,
49-
marker="x",
50-
label=f"{met} errors from sampling",
51-
)
52-
axs[i, 1].plot(d["ls"], berns_vf, c=cs[0], label=f"{met} errors from fitted IR")
53-
# axs[i, 1].plot(d[f'ls'], d[f'err_model_{met}'], c=cs[0], label = f'{met} errors from bernouli')
54-
y = d[f"tpc_{met}"][
55-
0
56-
] # 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!
57-
x = d[f"tpc_{met}_dist"]
58-
y = np.array(y)
59-
60-
# TODO erase this afterwards:
61-
if met == "vf":
62-
ir = np.round(d[f"ir_vf"], 1)
63-
axs[i, 2].plot(x, y, c=cs[1], label=f"Volume fraction 2PC")
64-
axs[i, 2].axhline(d["vf"] ** 2, linestyle="dashed", label="$p^2$")
65-
axs[i, 2].plot(
66-
[0, ir],
67-
[d["vf"] ** 2 - 0.02, d["vf"] ** 2 - 0.02],
68-
c="green",
69-
linewidth=3,
70-
label=r"$\tilde{a}_2$",
71-
)
72-
ticks = [0, int(ir), 20, 40, 60, 80, 100]
73-
ticklabels = map(str, ticks)
74-
axs[i, 2].set_xticks(ticks)
75-
axs[i, 2].set_xticklabels(ticklabels)
76-
axs[i, 2].fill_between(
77-
x, d["vf"] ** 2, y, alpha=0.5, label=f"Integrated part"
78-
)
79-
axs[i, 2].legend()
80-
81-
# axs[i, 2].scatter(x[knee], y_data[knee]/y.max(), c =cs[1], marker = 'x', label=f'{met} ir from tpc', s=100)
82-
ir = d[f"ir_{met}"]
83-
pred_ir = util.tpc_to_ir(d[f"tpc_{met}_dist"], d[f"tpc_{met}"])
84-
pred_ir = pred_ir * 1.61
85-
# 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)
86-
# axs[i, 2].scatter(x[round(ir)], y[round(ir)], facecolors='none', edgecolors = cs[1], label=f'{met} fitted IR', s=100)
87-
88-
irs[j].append(ir)
89-
if i == 0:
90-
axs[i, 1].legend()
91-
axs[i, 2].legend()
92-
axs[i, 1].set_xlabel("Image length size [pixels]")
93-
axs[i, 1].set_ylabel("Volume fraction percentage error [%]")
94-
axs[i, 2].set_ylabel("2-point correlation function")
95-
ir = np.round(d[f"ir_vf"], 2)
96-
sas.append(d["sa"])
97-
im = img[0] * 255
98-
si_size, nirs = 160, 5
99-
sicrop = int(ir * nirs)
100-
print(ir, sicrop)
101-
subim = torch.tensor(im[-sicrop:, -sicrop:]).unsqueeze(0).unsqueeze(0).float()
102-
subim = interpolate(subim, size=(si_size, si_size), mode="nearest")[0, 0]
103-
subim = np.stack([subim] * 3, axis=-1)
104-
105-
subim[:5, :, :] = 125
106-
subim[:, :5, :] = 125
107-
# subim[5:20, 5:50, :]
108-
subim[10:15, 10 : 10 + si_size // nirs, :] = 0
109-
subim[10:15, 10 : 10 + si_size // nirs, 1:-1] = 125
110-
111-
im = np.stack([im] * 3, axis=-1)
112-
im[-si_size:, -si_size:] = subim
113-
axs[i, 0].imshow(im)
114-
axs[i, 0].set_xticks([])
115-
axs[i, 0].set_yticks([])
116-
axs[i, 0].set_ylabel(f"M{n[1:]}")
117-
axs[i, 0].set_xlabel(
118-
f"Volume fraction "
119-
+ r"$\tilde{a}_2$: "
120-
+ f"{ir} Inset mag: x{np.round(si_size/sicrop, 2)}"
121-
)
122-
123-
i += 1
41+
imsize = 1600
42+
lf = imsize//32 + 2 # the size of G's input
43+
many_images = util.generate_image(netG, lf=lf, threed=False, reps=10)
44+
pf = many_images.mean().cpu().numpy()
45+
small_imsize = 512
46+
img = many_images[0].detach().cpu().numpy()[:small_imsize, :small_imsize]
47+
48+
csets = [['black', 'black'], ['gray', 'gray']]
49+
conf = 0.95
50+
errs = util.real_image_stats(many_images, lens_for_fit, pf, conf=conf)
51+
sizes_for_fit = [[lf,lf] for lf in lens_for_fit]
52+
real_cls = util.fit_cls(errs, sizes_for_fit, pf)
53+
stds = errs / stats.norm.interval(conf)[1] * pf / 100
54+
std_fit = (real_cls**2/(np.array(lens_for_fit)**2)*pf*(1-pf))**0.5
55+
# print(stds)
56+
vars = stds**2
57+
# from variations to L_characteristic using image size and phase fraction
58+
clss = (np.array(lens_for_fit)**2*vars/pf/(1-pf))**0.5
59+
print(clss)
60+
61+
axs_twin = axs[i, 1].twinx()
62+
stds_scatter = axs_twin.scatter(lens_for_fit, stds, s=8, color=colors["stds"], label = f'Standard deviations')
63+
twin_fit = axs_twin.plot(lens_for_fit, std_fit, color=colors["stds"], label = f'Standard deviations fit')
64+
axs_twin.tick_params(axis='y', labelcolor=colors["stds"])
65+
axs_twin.set_ylim(0, 0.025)
66+
cls_scatter = axs[i, 1].scatter(lens_for_fit, clss, s=8, color=colors["cls"], label = f'Characteristic length scales')
67+
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')
68+
axs[i, 1].tick_params(axis='y', labelcolor=colors["cls"])
69+
axs[i, 1].set_ylim(0, 28)
70+
71+
dims = len(img.shape)
72+
# print(f'starting tpc radial')
73+
tpc = util.tpc_radial(img, threed=False)
74+
cls = util.tpc_to_cls(tpc, img, img.shape)
75+
76+
contour = axs[i, 2].contourf(tpc, cmap='plasma', levels = 200)
77+
divider = make_axes_locatable(axs[i, 2])
78+
cax = divider.append_axes('right', size='5%', pad=0.05)
79+
fig.colorbar(contour, cax=cax, orientation='vertical')
80+
81+
if i == 0:
82+
plots = [stds_scatter, twin_fit[0], cls_scatter, cls_fit]
83+
label_plots = [plot.get_label() for plot in plots]
84+
axs[i,1].legend(plots, label_plots)
85+
86+
# axs[i,2].legend()
87+
# axs[i,1].set_xlabel('Image length size [pixels]')
88+
# axs[i,1].set_ylabel('Volume fraction percentage error [%]')
89+
# axs[i,2].set_ylabel('2-point correlation function')
90+
# ir = np.round(d[f'ir_vf'], 2)
91+
# im = img[0]*255
92+
# si_size, nirs = 160, 5
93+
# sicrop = int(ir*nirs)
94+
# print(ir, sicrop)
95+
# subim=torch.tensor(im[-sicrop:,-sicrop:]).unsqueeze(0).unsqueeze(0).float()
96+
# subim = interpolate(subim, size=(si_size,si_size), mode='nearest')[0,0]
97+
# subim = np.stack([subim]*3, axis=-1)
98+
99+
# subim[:5,:,:] = 125
100+
# subim[:,:5,:] = 125
101+
# # subim[5:20, 5:50, :]
102+
# subim[10:15, 10:10+si_size//nirs, :] = 0
103+
# subim[10:15, 10:10+si_size//nirs, 1:-1] = 125
104+
105+
# im = np.stack([im]*3, axis=-1)
106+
# im[-si_size:,-si_size:] = subim
107+
# axs[i, 0].imshow(im)
108+
# axs[i, 0].set_xticks([])
109+
# axs[i, 0].set_yticks([])
110+
# axs[i, 0].set_ylabel(f'M{n[1:]}')
111+
# axs[i, 0].set_xlabel(f'Volume fraction '+ r'$\tilde{a}_2$: '+ f'{ir} Inset mag: x{np.round(si_size/sicrop, 2)}')
112+
113+
i+=1
124114

125115

126116
plt.tight_layout()
127-
plt.savefig("fig2.pdf", format="pdf")
117+
plt.savefig('fig2.pdf', format='pdf')
128118

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

156146
# print(res,res2)
147+
148+
149+
150+
151+

0 commit comments

Comments
 (0)