In [1]:
from packages.display.core import *
%pylab inline
pylab.rcParams['figure.figsize'] = (15, 2.5)
In [ ]:
%qtconsole
https://www.iram.fr/IRAMFR/ARC/documents/cycle0/ALMA_EarlyScience_Cycle0_HighestPriority.pdf
In [61]:
file_path = '../data/2011.0.00419.S/sg_ouss_id/group_ouss_id/member_ouss_2013-03-06_id/product/IRAS16547-4247_Jet_CH3OH7-6.clean.fits'
noise_pixel = (15, 4)
train_pixels = [(133, 135),(134, 135),(133, 136),(134, 136)]
In [62]:
img = fits.open(file_path)
meta = img[0].data
hdr = img[0].header
# V axis
naxisv = hdr['NAXIS3']
onevpix = hdr['CDELT3']*0.000001
v0 = hdr['CRVAL3']*0.000001
v0pix = int(hdr['CRPIX3'])
vaxis = onevpix * (np.arange(naxisv)+1-v0pix) + v0
In [63]:
values = meta[0, :, train_pixels[0][0], train_pixels[0][1]] - np.mean(meta[0, :, train_pixels[0][0], train_pixels[0][1]])
values = values/np.max(values)
plt.plot(vaxis, values)
plt.xlim(np.min(vaxis), np.max(vaxis))
plt.ylim(-1, 1)
gca().xaxis.set_major_formatter(FormatStrFormatter('%d'))
In [64]:
noise = meta[0, :, noise_pixel[0], noise_pixel[1]] - np.mean(meta[0, :, noise_pixel[0], noise_pixel[1]])
noise = noise/np.max(noise)
plt.plot(vaxis, noise)
plt.ylim(-1, 1)
plt.xlim(np.min(vaxis), np.max(vaxis))
gca().xaxis.set_major_formatter(FormatStrFormatter('%d'))
In [10]:
cube_params = {
'freq' : vaxis[naxisv/2],
'alpha' : 0,
'delta' : 0,
'spe_bw' : naxisv*onevpix,
'spe_res' : onevpix*v0pix,
's_f' : 8,
's_a' : 0}
dictionary = gen_all_words(cube_params, True)
In [29]:
prediction = pd.DataFrame([])
for train_pixel in train_pixels:
dictionary_recal, detected_peaks = recal_words(file_path, dictionary, cube_params,
train_pixel, noise_pixel)
X = get_values_filtered_normalized(file_path, train_pixel, cube_params)
y_train = get_fortran_array(np.asmatrix(X))
dictionary_recal_fa = np.asfortranarray(dictionary_recal,
dtype= np.double)
lambda_param = 0
for idx in range(0, len(detected_peaks)):
if detected_peaks[idx] != 0:
lambda_param += 1
param = {
'lambda1' : lambda_param,
# 'L': 1,
'pos' : True,
'mode' : 0,
'ols' : True,
'numThreads' : -1}
alpha = spams.lasso(y_train, dictionary_recal_fa, **param).toarray()
total = np.inner(dictionary_recal_fa, alpha.T)
for i in range(0, len(alpha)):
iso_col = dictionary_recal.columns[i]
if(not prediction.columns.isin([iso_col]).any()):
prediction[iso_col] = alpha[i]
else:
prediction[iso_col] = prediction[iso_col]*alpha[i]
In [60]:
for p in prediction.columns:
if(prediction[p][0] != 0):
print(prediction[p])
In [76]:
latexify(8)
# Step 1: Read Cube
ax = plt.subplot(6, 1, 1)
data = get_data_from_fits(file_path)
y = data[0, :, train_pixel[0], train_pixel[1]]
plt.xticks([])
plt.plot(vaxis, y)
lines = get_lines_from_fits(file_path)
for line in lines:
# Shows lines really present
isotope_frequency = int(line[1])
isotope_name = line[0] + "-f" + str(line[1])
plt.axvline(x=isotope_frequency, ymin=0, ymax= 3, color='g')
# 2. Normalize, filter dada
ax = plt.subplot(6, 1, 2)
plt.ylim(ymin =0,ymax = 1.15)
y = get_values_filtered_normalized(file_path, train_pixel, cube_params)
plt.xticks([])
plt.plot(vaxis, y)
# 3. Possible Words
ax = plt.subplot(6, 1, 3)
plt.ylim(ymin =0,ymax = 1.15)
plt.xticks([])
plt.plot(vaxis, dictionary)
# 4. Detect Lines
ax = plt.subplot(6, 1, 4)
plt.ylim(ymin =0,ymax = 1.15)
plt.plot(vaxis, y)
plt.xticks([])
plt.ylabel("Temperature")
for idx in range(0, len(detected_peaks)):
if detected_peaks[idx] != 0:
plt.axvline(x=vaxis[idx], ymin=0, ymax= 1, color='r')
# 6. Recalibrate Dictionary
ax = plt.subplot(6, 1, 5)
plt.ylim(ymin =0,ymax = 1.15)
plt.plot(vaxis, dictionary_recal_fa)
plt.xticks([])
# 6. Recover Signal
ax = plt.subplot(6, 1, 6)
plt.ylim(ymin =0,ymax = 1.15)
plt.plot(vaxis, total)
plt.xlabel("Frequency [MHz]")
gca().xaxis.set_major_formatter(FormatStrFormatter('%d'))
In [68]:
def latexify(fig_width=None, fig_height=None, columns=1):
"""Set up matplotlib's RC params for LaTeX plotting.
Call this before plotting a figure.
Parameters
----------
fig_width : float, optional, inches
fig_height : float, optional, inches
columns : {1, 2}
"""
# code adapted from http://www.scipy.org/Cookbook/Matplotlib/LaTeX_Examples
# Width and max height in inches for IEEE journals taken from
# computer.org/cms/Computer.org/Journal%20templates/transactions_art_guide.pdf
assert(columns in [1,2])
if fig_width is None:
fig_width = 4.89 if columns==1 else 6.9 # width in inches
if fig_height is None:
golden_mean = (sqrt(5)-1.0)/2.0 # Aesthetic ratio
fig_height = fig_width*golden_mean # height in inches
MAX_HEIGHT_INCHES = 24.0
if fig_height > MAX_HEIGHT_INCHES:
print("WARNING: fig_height too large:" + fig_height +
"so will reduce to" + MAX_HEIGHT_INCHES + "inches.")
fig_height = MAX_HEIGHT_INCHES
params = {'backend': 'ps',
'text.latex.preamble': ['\usepackage{gensymb}'],
'axes.labelsize': 8, # fontsize for x and y labels (was 10)
'axes.titlesize': 8,
'text.fontsize': 8, # was 10
'legend.fontsize': 8, # was 10
'xtick.labelsize': 10,
'ytick.labelsize': 8,
'text.usetex': True,
'figure.figsize': [fig_width,fig_height],
'font.family': 'serif'
}
matplotlib.rcParams.update(params)
def format_axes(ax):
for spine in ['top', 'right']:
ax.spines[spine].set_visible(False)
for spine in ['left', 'bottom']:
ax.spines[spine].set_color(SPINE_COLOR)
ax.spines[spine].set_linewidth(0.5)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
for axis in [ax.xaxis, ax.yaxis]:
axis.set_tick_params(direction='out', color=SPINE_COLOR)
return ax
In [ ]:
for i in range(0, len((alpha > 0))):
if((alpha > 0)[i]):
print(dictionary_recal.columns[i])
print(prediction)
In [ ]:
for i in range(0, len(dictionary.index)):
print(calculate_probability(alpha, dictionary.index[i], dictionary_recal))
print(dictionary.index[i])