<jupyter_start><jupyter_text>MACS1206 1D analysis<jupyter_code>from IPython.core.display import display, HTML
display(HTML("<style>.container { width:95% !important; }</style>"))

%matplotlib notebook
import matplotlib.pylab as plt
import seaborn as sns
sns.set(style='darkgrid')

import numpy as np
from scipy import stats

from astropy.stats import sigma_clipped_stats
from astropy.io import fits
from astropy.wcs import WCS
from astropy.convolution import convolve,Gaussian2DKernel
from reproject import reproject_exact,reproject_interp

import linmix
from sklearn.metrics import mean_absolute_error

import warnings
warnings.filterwarnings('ignore')<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code>met_map = fits.getdata('Maps/Map_metallicity.fits')
emet_map = fits.getdata('Maps/Map_metallicity_unc.fits')
header = fits.getheader('Maps/Map_metallicity.fits')
dist_map, _ = reproject_interp('../../Data/Lensing/MACS1206/simul_MACS1206_distance_kpc_source_plane.fits',header) 

# The interpolation and lensing makes the borders (where the distance is greater) to have smaller values. we're correcting for this, a little manually
bad_pix = np.where(dist_map==0)
for p,q in zip(bad_pix[0],bad_pix[1]):
    dist_map[p-2:p+2,q-2:q+2] = np.nan
    
## Convolve the distance map with the seeing
seeing_fwhm_pix = 0.78/0.12
dist_map_conv = convolve(dist_map,Gaussian2DKernel(seeing_fwhm_pix/2.355),boundary='extend')

# Plotting
fig, ax = plt.subplots(1,2)
ax[0].imshow(dist_map)
ax[1].imshow(dist_map_conv)<jupyter_output><empty_output><jupyter_text>Measure the radial profiles<jupyter_code>delta_r = 1
extract_annuli = np.arange(0,8,delta_r)
def measure_in_annuli(im,dist,stat='mean'):
    if stat == 'mean':
        profile = [ sigma_clipped_stats(im[np.where((dist >= r) & (dist < r +delta_r))])[0] for r in extract_annuli ]
    if stat == 'std':
        profile = [ sigma_clipped_stats(im[np.where((dist >= r) & (dist < r +delta_r))])[2] for r in extract_annuli ]
    return np.array(profile)


# Measure gradients
met_annuli = measure_in_annuli(met_map,dist_map,'mean') 
met_annuli_std = measure_in_annuli(met_map,dist_map,'std') 
emet_annuli = measure_in_annuli(emet_map,dist_map,'mean')
met_annuli_err = np.sqrt(emet_annuli**2 + met_annuli_std**2)

dist_annuli = np.arange(0.5,8.5,delta_r)
dist_annuli_std = measure_in_annuli(dist_map,dist_map,'std')


# Measure gradients in distance convolved map
met_conv_annuli     = measure_in_annuli(met_map,  dist_map_conv, 'mean') 
met_conv_annuli_std = measure_in_annuli(met_map,  dist_map_conv, 'std') 
emet_conv_annuli    = measure_in_annuli(emet_map, dist_map_conv, 'mean')
met_conv_annuli_err = np.sqrt(emet_conv_annuli**2 + met_conv_annuli_std**2)

dist_conv_annuli = np.arange(0.5,8.5,delta_r)
dist_conv_annuli_std = measure_in_annuli(dist_map,dist_map_conv,'std')

# Cleaning nans
dist_conv_annuli = dist_conv_annuli[2:]
met_conv_annuli = met_conv_annuli[2:] 
dist_conv_annuli_std = dist_conv_annuli_std[2:]
met_conv_annuli_err = met_conv_annuli_err[2:]<jupyter_output><empty_output><jupyter_text>Find the radius for which we get the lowest residuals<jupyter_code>def fit_1d_gradient(x,y,xsig,ysig,min_r,max_r,ax=None,color=None,plot=True):
    
    # Crop data
    x_full, y_full, xsig_full, ysig_full  = x , y, xsig, ysig    
    x = x[min_r:max_r]
    y = y[min_r:max_r]
    xsig = xsig[min_r:max_r]
    ysig = ysig[min_r:max_r]

    # Fit
    lm = linmix.LinMix(x, y, xsig, ysig=ysig, K=2)
    lm.run_mcmc(silent=True)
    start = int(len(lm.chain)*0.5)
    
    mean_alpha, std_alpha = np.mean(lm.chain['alpha'][start:]), np.std(lm.chain['alpha'][start:])
    mean_beta, std_beta = np.mean(lm.chain['beta'][start:]), np.std(lm.chain['beta'][start:])
    y_flat = np.ones_like(x_full[max_r-1:])*np.mean(y_full[max_r-1:])
    mae_grad = mean_absolute_error(y, mean_alpha + mean_beta*x)
    mae_lin = mean_absolute_error(y_full[max_r-1:], y_flat)
    mae = np.mean((mae_grad,mae_lin))
    chi2 = (np.nansum((y-(mean_alpha + mean_beta*x))**2/ysig**2))/(len(y)-2)    

    if plot:
        # model
        for i in range(start, len(lm.chain), 20):
            ys = lm.chain[i]['alpha'] +  lm.chain[i]['beta']* x 
            ax.plot(x, ys, color=color, alpha=0.02,zorder=1)
        ax.plot(x,mean_alpha + mean_beta*x, color=color,zorder=2) 
        ax.plot(x_full[max_r-1:],y_flat, color=color,zorder=3)  

        ax.annotate('m = %0.3f$\pm$%0.3f\ny${_0}$ = %0.2f$\pm$%0.2f\n$\chi^2/dof$ = %0.2f'
                %(mean_beta,std_beta,mean_alpha,std_alpha,chi2),
                xy=(0.50, 0.70),xycoords='axes fraction')
        #data
        ax.errorbar(x_full, y_full, xerr = xsig_full, yerr = ysig_full,
               marker='o',mfc='0.4',mec='k',mew=0.5,linestyle='',
               ecolor='0.1',elinewidth=0.7,capsize=1)
        ax.plot(x,y,marker='o',color='0.2',linestyle='',zorder=10,alpha=0.9)

   
    return mae<jupyter_output><empty_output><jupyter_text>Test the flatenning hypothesis<jupyter_code>fig, ax = plt.subplots(1,4,figsize=(14,5))
fig.subplots_adjust(left=0.05,right=0.99)
ax = ax.ravel()
dummy = [x.set_xlim(0,10) for x in ax]
dummy = [x.set_ylim(8.0,10) for x in ax]

mae_array = [fit_1d_gradient(dist_annuli,met_annuli,dist_annuli_std,met_annuli_err,min_r=1,max_r=r,plot=True,ax=ax[i]) for  i,r in enumerate(np.arange(8,4,-1))]    
plt.legend()
fig, ax = plt.subplots(1,1)
ax.set_xlim(0,10)
ax.set_ylim(8.0,10) 

fit_1d_gradient(dist_conv_annuli,met_conv_annuli,dist_conv_annuli_std,met_conv_annuli_err,min_r=1,max_r=6,plot=True,ax=ax)<jupyter_output><empty_output><jupyter_text>Plotting<jupyter_code># Figure
fig, ax = plt.subplots(1,1,figsize=(5,4))
fig.subplots_adjust(left=0.2,top=0.90,bottom=0.2)

fig.suptitle('M1206-sys1')
ax.set_xlabel('Radius (kpc)')
ax.set_ylabel('12 + log(O/H)')

ax.set_ylim(8.5,9.5)
ax.set_xlim(-0.5,9)

metfit = fit_1d_gradient(dist_annuli,met_annuli,dist_annuli_std,met_annuli_err,0,8,ax,'#960056')   


fig.savefig('../../Plots/MACS1206_anulli.pdf')
# Figure
fig, ax = plt.subplots(1,1,figsize=(5,3))
fig.subplots_adjust(left=0.2,top=0.90,bottom=0.2)

fig.suptitle('M1206-sys1')
ax.set_xlabel('Radius (kpc))')
ax.set_ylabel('12 + log(O/H)')

ax.set_ylim(8.5,9.5)
ax.set_xlim(-0.5,9)

metfit = fit_1d_gradient(dist_conv_annuli, met_conv_annuli,dist_conv_annuli_std,met_conv_annuli_err,0,6,ax,'#cb416b')   

#fig.savefig('../../Plots/MACS1206_anulli.pdf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15406206143770465, 'fraction_numerical': 0.03446125058474973, 'mean_word_length': 5.544897959183674, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '135268', 'license': 'MIT', 'n_tokens_mistral': 2847, 'n_tokens_neox': 2746, 'n_words': 367, 'path': 'Metallicity/MACS1206/macs1206_1D_analysis.ipynb', 'repo_name': 'VeraPatricio/Resolved_Metallicity', 'size': '6413'}