In [1]:
%cd H:\GitHub\pytracer\
In [2]:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
%matplotlib inline
import numpy as np
from scripts.utils import nice_double_plot
import pytracer.algorithms as algorithms
In [ ]:
response_single = np.load(r'scripts\data\fission_response_single.npy')
response_double = np.load(r'scripts\data\fission_response_double.npy')
single_probs = np.load(r'scripts\data\single_probs.npy')
double_probs = np.load(r'scripts\data\double_probs.npy')
d1 = single_probs.reshape((-1))
G1 = response_single.reshape((response_single.shape[0], -1)).T
d2 = single_probs.reshape((-1))
G2 = response_single.reshape((response_double.shape[0], -1)).T
In [ ]:
alphas = np.logspace(-5, -1, 100)
norms, residuals = algorithms.trace_lcurve(d1, G1, alphas)
plt.loglog(residuals, norms)
plt.figure()
c_alphas, curve = algorithms.lcurve_curvature(alphas, norms, residuals)
plt.semilogx(c_alphas, curve)
m1_alpha = algorithms.solve_tikhonov(d1, G1, 0.007)
m1_alpha = m1_alpha.reshape((30, 50))
plt.figure()
plt.imshow(m1_alpha, interpolation='none', vmax=1.0, vmin=0)
plt.colorbar()
In [ ]:
alphas = np.logspace(-5, -1, 100)
norms, residuals = algorithms.trace_lcurve(d2, G2, alphas)
plt.loglog(residuals, norms)
plt.figure()
c_alphas, curve = algorithms.lcurve_curvature(alphas, norms, residuals)
plt.semilogx(c_alphas, curve)
m2_alpha = algorithms.solve_tikhonov(d2, G2, 0.007)
m2_alpha = m2_alpha.reshape((30, 50))
plt.figure()
plt.imshow(m2_alpha, interpolation='none', vmax=1.0, vmin=0)
plt.colorbar()
In [4]:
d = np.concatenate((d1, d2)) * 0.2
G = np.concatenate((G1, G2))
alphas = np.logspace(-6, 0, 50)
norms, residuals = algorithms.trace_lcurve(d, G, alphas)
plt.figure()
c_alphas, curve = algorithms.lcurve_curvature(alphas, norms, residuals)
plt.semilogx(c_alphas, curve)
best_alpha_arg = np.argmax(curve)
best_alpha = c_alphas[best_alpha_arg]
norms[0] *= 10
plt.figure()
plt.loglog(residuals, norms, 'k', lw=2)
plt.xlabel(r'$\log{||A X - Y||_2^2}$', size=20)
plt.ylabel(r'$\log{||X||_2^2}$', size=20)
plt.title(r'L-Curve Criterion with Optimal $\alpha$ = {:.3f}'.format(best_alpha), size=15)
plt.scatter(residuals[best_alpha_arg], norms[best_alpha_arg], c='red', s=100, edgecolors='none')
plt.tight_layout()
m_alpha = algorithms.solve_tikhonov(d, G, 0.01)
m_alpha = m_alpha.reshape((30, 50))
plt.figure(figsize=(6, 4))
ax = plt.gca()
im = plt.imshow(m_alpha, interpolation='none', vmax=0.2, vmin=0, extent=[-25. / 2, 25. / 2, -15. / 2, 15. / 2],
cmap='viridis')
plt.title('Single & Double Reconstruction', size=20)
plt.xlabel('X (cm)', size=18)
plt.ylabel('Y (cm)', size=18)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cb = plt.colorbar(im, cax=cax)
cb.set_label(r'$\mu_f / \mu_{total}$', size=20, labelpad=15)
plt.show()