<jupyter_start><jupyter_text>Demo of MuSCADeT on the Refsdal galaxy<jupyter_code>#Import libraries
import astropy.io.fits as pf
import numpy as np
import matplotlib.pyplot as plt
from MuSCADeT import MCA
from MuSCADeT  import colour_subtraction as cs<jupyter_output><empty_output><jupyter_text>At this point, the data are loaded.<jupyter_code>## Openning data cube
cube = pf.open('./Simu_Refsdal_big/Cube.fits')[0].data
nb,n1,n2 = np.shape(cube)
#Here A was estimated from colours of other galaxies in the field
Aprior =pf.open('Simu_Refsdal_big/Estimated_A_PCA.fits')[0].data<jupyter_output><empty_output><jupyter_text>Setting up the parameters for MuSCADeT and running.<jupyter_code>##param
pca = 'noPCA'   #Estimation of the mixing coefficients from PCA. If different from PCA it will use the array provided in Aprior
n = 2000        #Number of iterations
nsig = 5        #Threshold in units of noise standard deviation
ns = 2          #Number of sources
angle = 50      #Resolution angle for the PCA colour estimation (start with 15 then adjust empirically)


## Running MuSCADeT
S,A = MCA.mMCA(cube, Aprior, nsig,n, PCA=[ns,angle], mode=pca, harder = 1, lvl=6)<jupyter_output>mode noPCA<jupyter_text>Display of the results. The first plot shows the model and a comparison with the original image. The model is simply `AS` with `A` and `S` as estimated by MuSCADeT.<jupyter_code># Models as extracted by MuSCADeT for display
model = np.dot(A,S.reshape([A.shape[1], n1*n2])).reshape(cube.shape)

normodel = cs.asinh_norm(model, Q=100, range = 1)
normcube = cs.asinh_norm((cube), Q = 100, range = 1)
normres = cs.asinh_norm(cube-model, Q = 10, range = 1)
plt.figure(figsize = (15, 5))
plt.subplot(131)
plt.title('model')
plt.imshow(normodel)
plt.subplot(132)
plt.title('data')
plt.imshow(normcube)
plt.subplot(133)
plt.title('Residuals')
plt.imshow(normres)
plt.show()<jupyter_output><empty_output><jupyter_text>This sequence of plots shows each component, with what the data looks like once the components are removed, i.e. `Y-AiSi`.<jupyter_code>for i in range(A.shape[1]):
    C = A[:,i, np.newaxis, np.newaxis]*S[np.newaxis,i,:,:]
    normC = cs.asinh_norm(C, Q = 100, range = 1)
    normCres = cs.asinh_norm((cube-C), Q = 100, range = 1)
    plt.figure(figsize = (15, 5))
    plt.subplot(131)
    plt.title('data')
    plt.imshow(normcube)
    plt.subplot(132)
    plt.title('component ' + str(i))
    plt.imshow(normC)
    plt.subplot(133)
    plt.title('data - component ' + str(i))
    plt.imshow(normCres)
    plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12417024599765716, 'fraction_numerical': 0.024599765716516987, 'mean_word_length': 4.809523809523809, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 20, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '133416', 'license': 'MIT', 'n_tokens_mistral': 1034, 'n_tokens_neox': 956, 'n_words': 254, 'path': 'Examples/Example_Refsdal.ipynb', 'repo_name': 'herjy/M2CAD', 'size': '2561'}