How to use MUSTARD estimator
An script version of this notebook is available here
Build the estimator
MUSTARD is object-oriented. The estimator object need the ADI cube centred and normalized and the angle to be build.
First, we import the package we will need (don”t forget to install mustard)
[1]:
# Mustard estimator
from mustard import mustard_estimator
# Also, this will help you build the mask fo regularization
from mustard.utils import gaussian
# Misc
from vip_hci.fits import open_fits
from vip_hci.preproc import cube_crop_frames
Fits HDU-0 data successfully loaded. Data shape: (66, 256, 256)
Fits HDU-0 data successfully loaded. Data shape: (66,)
Then, we load the data and build the estimator
[ ]:
# Load the ADI cube (1 channel) and the associated angles
datadir = "./example-data/"
science_data = open_fits(datadir+"cube")
angles = open_fits(datadir+"angles")
# Note : Don't hestiate to crop the cubes it is time/ressources-consuming
science_data = cube_crop_frames(science_data, 256)
estimator = mustard_estimator(science_data, angles=angles, coro=10, pupil="edge", Badframes=None, hid_mask=None)
This operation builed the object that represent the estimator. At this point, it only stores information about the dataset.
Regularization
Circular ambiguities
State-of-the-art algorithms suffer from a lack of a circularly invariant component in the disk estimate due to their incapacity of distinguishing the disk’s flux invariant to the rotation from a static component (i.e., speckle field). In a dataset, both disk and speckle field morphologies contain flux invariant to the rotation. Depending on how aggressively the algorithm will remove the quasi-static contribution from the dataset, more deformation will appear.
We aim to correct for the geometrical biases of standard PSF-subtraction techniques, which generally assign any rotation-invariant flux fully to the reconstructed PSF, by making an assumption on the morphology of the speckle field.
For that purpose, we will defin a Gaussians with a width set to match the spread of the stellar halo.
While fine-tuning of the mask from a dataset to another is expected to improve our results, we rather considered a realistic case scenario of not knowing beforehand the stellar halo profile due to the potential presence of bright disk signals. \texttt{MUSTARD} is initialized with the mean derotated ADI sequence (as in noADI) in order to assign all ambiguous flux to the disk as first estimate. Hence, this regularization will push out the ambiguous flux that does not belong to the disk accordingly to the shape of the mask.
In addition, we define a smoothness regularization term using the spatial gradient of the estimated circumstellar and stellar signals \(d\) and \(s\), respectively. It is a common regularization to compensate for noise.
[3]:
# Create a mask
shape = estimator.model.frame_shape # Get the frame shape here.
M = circle(shape, shape[0]//2) + 10*circle(shape, 13)
# Configure R2
estimator.configR2(Msk=None, mode="l1", penaliz="X", invert=True)
# Configure R1
estimator.configR1(mode="smooth", p_L=0.5)
Parameters
Finally, MUSTARD required to set few parameters (you can keep defaults) :
[2]:
param = {'w_r' : 0.15, # Proportion of Regul over J
'w_r2' : 0.03, # Proportion of Regul2 over J
'w_r3' : 0.001, # Proportion of Regul2 over J
'w_way' : (1, 0), # You‡ can either work with derotated_cube or rotated cube. Or both
'gtol' : 1e-100, # Gradient tolerence. Stop the estimation when the mean of gradient will hit the value
'kactiv' : 0, # Iter before activate regul (i.e when to compute true weight base on w_r proportion)
'estimI' : "None", # Estimate frames flux is highly recommended ! possible value : {"Frame","L","Both"}
'weighted_rot' : False, # Compute weight for each frame according to PA angle separations.
'suffix' : "test", # # Name of your simulation (this is optional)
'maxiter': 13,
'mask_L': (0 , 100),
'init_maxL': True}
-w_r and w_r2, define respectivly the weight of the regularization terms to smooth and to sort circular ambiguity (the last one will be explain later). A weight between 10%-5% is recommanded in order to let the data attachment terme drive the estimation to avoid deformations. R2 can be higher if you purposly want to get rid of all circular ambiguities from the disk estimations.
-w_r define the model definition. If your signal is corrupted by Gibbs artifact, you can swhitch to reverse mode. In reverse mode, the speckels map is rotating and it is less likely to cause such problem. You can also try both way at the same time but it takes more time and will not necessarly give better results.
Start the minimization
The demo from “exemple-data” usually takes ~4mins.
TIP : Click on the ◼︎ icon to stop prematurely the minimization (or ctrl+C if you are using a console). It will quit proprely the iterative loop. Results will be stored/saved/return.
[ ]:
L_est, X_est = estimator.estimate(**param, save=datadir, gif=False, verbose=True)
__________________________________________________
Resolving IP-ADI optimization problem - name :
Outputs will be saved ./example-data/
Regul R1 : 'smooth' and R2 : 'l1 on X inverted'
No deconvolution and with frame weighted based on rotations
Relative amplitude of BothX and L will be estimated
Regul weight are set to w_r=5.00e-02 and w_r2=1.00e-01, maxiter=60
|it | loss | R1 | R2 | Rpos | total |
| 0 |6.503639e-05 (100%)|0.000000e+00 (0 %)|0.000000e+00 (0 %)|0.000000e+00 (0 %)|6.503639112875e-05|
INFO : The iteration k+1/2 is the step where regularization weights are computed.
Results
The estimator store all the results. It does also provide some method to get relevant information.
[ ]:
import matplotlib.pyplot as plt
from numpy import percentile
plt.figure("Results",figsize=(16,9))
plt.subplot(121), plt.imshow(L_est,cmap='jet',vmax=percentile(L_est,99.5))
plt.subplot(122), plt.imshow(X_est,cmap='jet',vmax=percentile(L_est,99))
Evolution of the criteria
This method return the array of the values of the criteria at each iteration. If you set show=True or save=$path$, it will also plot the convergence curve.
[ ]:
evo = estimator.get_evo_convergence()
Frame weight
The frame weight is computed based on the PA angle vector of each frame. Indeed for the sake of the estimation, if the angle between two frame is two small, there is chance that the signal will overllap and it can biais the results.
The fuction below return the vector of frame weights. If you set show=True or save=$path$, it will also plot the weight bars and the PA angle cuve.
[ ]:
wr = estimator.get_rot_weight()
Residual
This is the noise. If you see some structures that is alright. This map containt every kind of noise - not only white. MUSTARD can distangle noise that have a life time between ( exposure time < lifetime < aquisition-time/2 ).
Look carefully, because if you see very clear structes that look like the disk and/or that appear on more that half of the frames it might be a sign that something went wrong.
[ ]:
from numpy import percentile
from hciplot import plot_cubes
residual = estimator.get_residual()
lim = percentile(abs(residual),99) # TIP : try without percentile.
plot_cubes(residual,cmap="jet",vmax=lim,vmin=-lim)
Reconstruction
This is the reconstruction. It should look like the initial ADI cube.
[ ]:
reconstruction = estimator.get_reconstruction()
plot_cubes(reconstruction, cmap="jet", vmax=percentile(reconstruction,99))
BONUS : You can generate the mustard gif of your simulation (the one in introcution)
[ ]:
estimator.mustard_results()
The result will be print saved in your datadir.