BARNACLE’s main script

This code determines the astrophysical null depth either by the NSC method (https://ui.adsabs.harvard.edu/abs/2011ApJ…729..110H/abstract) which does not necessarily requires a calibrator, either by a classic measurement (https://ui.adsabs.harvard.edu/abs/2011ApJ…729..110H/abstract) which must be calibrated.

The fit is done via the fitting algorithm Trust Region Reflective used in the scipy function least_squares (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html). In order to be more resilient to local minima, basin hopping is practiced 10 times by slighly changing the initial conditions (see the documentation of the function basin_hoppin_values below).

NOTE: the code is not fully stabilised yet so some discarded functions exist but are not used anymore. Similarly, the execution is highly customisable due to the different explorations done.

The data to load and the settings about the space parameters are handled by the script glint_fitting_config6.py. This script handles the way the null depth is determined and on which baseline. The library is glint_fitting_function6.py. However, some core functions are coded here:

  • gaussian: to generate Gaussian curves;
  • MCfunction: the core function to create and fit the histogram of synthetic null depths to the real histograms;
  • Logger: class logging all the console outputs into a log file;
  • basin_hoppin_values: function changing initial condition to do basin hopping on the fitting algorithm.

Fitted values are stored into npy file (one per basin iteration in case of breakdown of the program). Pickle files stored the results from all iterations. Plots of the fit of the histogram and plots of the intermediate quantities (histogram of (anti-)null outputs, synthetic or real, photometries, injection) are also saved.

The settings of the script are defined by:

  • linear_model: bool, use the linear model \(N = Na + Ninstr\);
  • activate_use_antinull: bool. If True, it use the antinull output;
  • activate_select_frames: bool. If True, it keeps frames for which the null depths are between boundaries for every spectral channels;
  • activate_random_init_guesses: bool. If True, it activates the basin hopping. It is automatically set to True after the first iteration;
  • activate_spectral_binning: bool. If True, it bins the spectral channels;
  • activate_use_photometry: bool. If True, it uses the photometric outputs instead of simulating sequences of photometries for the fit. It is advised to keep it at False as its use gives non-convincing results and more investigation are required;
  • activate_dark_correction: bool. If True, it corrects the distribution of the injection from the contribution of the detector noise assuming they all follow normal distribution. It is adviced to keep it at False;
  • activate_save_basic_esti: bool. If True, it gives a broadband null depth to calibrate in the old-school way;
  • activate_frame_sorting: bool. If True, it keeps the frame for which the OPD follows a normal distribution (anti-LWE filter). It is strongly recommended to keep it True.
  • activate_preview_only: bool. If True, the script stops after sorting the frame according to the phase distribution and plot the sequences of frames with highlighting the kept ones to tune the parameters of the filter;
  • nb_frames_sorting_binning: integer, number of frames to bin to set the frame sorting parameters. Typical value is 100;
  • activate_time_binning_photometry: bool. If True, it bins the photometric outputs temporally. It is strongly recommended to keep it True to mitigate the contribution of the detector noise;
  • nb_frames_binning_photometry: integer, number of frames to bin the photometric outputs. Typical value is 100;
  • global_binning: integer, frames to bin before anything starts in order to increase SNR;
  • wl_min: float, lower bound of the bandwidth to process;
  • wl_max: float, upper bound of the bandwidth to process;
  • n_samp_total: int, number of total Monte-Carlo values to generate. Typical value is 1e+8;
  • n_samp_per_loop: int, number of Monte-Carlo values to generate per loop (to not saturate the memory). Typical value is 1e+7;
  • phase_bias_switch: bool. If True, it implements a non-null achromatic phase in the null model;
  • opd_bias_switch: bool. If True, it implements an offset OPD in the null model;
  • zeta_switch: bool. If True, it uses the measured zeta coeff. If False, value are set to 0.5;
  • oversampling_switch: bool. If True, it includes the loss of coherence due to non-monochromatic spectral channels;
  • skip_fit: bool. If True, the histograms of data and model given the initial guesses are plotted without doing the fit. Useful to just display the histograms and tune the settings;
  • chi2_map_switch: bool. If True, it grids the parameter space in order to set the boundaries for the fitting algorithm;
  • nb_files_data: 2-int tuple. Load the data files comprised between the values in the tuple;
  • nb_files_dark: 2-int tuple. Load the dark files comprised between the values in the tuple;
  • basin_hopping_nloop: 2-int tuple. Lower and upper bound of the iteration loop for basin hopping method;
  • which_nulls: list. List of baselines to process;
  • map_na_sz: int. Number of values of astrophysical null depth in the grid of the parameter space;
  • map_mu_sz: int. Number of values of \(\mu\) in the grid of the parameter space;
  • map_sig_sz: int. Number of values of \(\sig\) in the grid of the parameter space.
class glint_fitting_gpu6.Logger(log_path)

Class allowing to save the content of the console inside a txt file.

glint_fitting_gpu6.MCfunction(bins0, na, mu_opd, sig_opd, *args)

Monte-Carlo function ofr NSC method. It fits polychromatic data to give one achromatic null depth. This functions runs on the GPU with cupy (https://cupy.dev/).

The arguments of the function are the parameters to fit. The other values required to calculated the synthetic sequences of null depth are imported via global variables listed at the beginning of the function.

The verbose of this function displays the input parameters to track the convergence of the fitting algorithm.

Global variables include quantities necessary to fit the histograms and intermediate products to control the efficiency of the fit.

The functions loops over the spectral channels to create independent histograms for each of them and loops over the Monte-Carlo samples to avoid the saturation of the memory. For the latter, the histograms are added together then averaged to get the normalization to 1.

Parameters:
bins0: array

Bins of the histogram of the null depth to fit

na: float

Astrophysical null depth, quantity of interest.

mu_opd: float

Instrumental OPD of the considered baseline.

sig_opd: float

Sigma of the normal distribution describing the fluctuations of the OPD.

args: optional

Array, use to fit the parameters mu and sig of the normal distribution of Ir (https://ui.adsabs.harvard.edu/abs/2011ApJ…729..110H/abstract) if the antinull output is not used. arg can be either a 2-element or a \(2 imes number of spectral channels\)-element array, depending on the fit of a unique Ir or one per spectral channel.

Returns:
accum: array

Histogram of the synthetic sequences of null depth.

TODO:

Implement the input of args.

glint_fitting_gpu6.basin_hoppin_values(mu_opd0, sig_opd0, na0, bounds_mu, bounds_sig, bounds_na)

Create as many as initial guess as there are basin hopping iterations to do. The generation of these values are done with a normal distribution.

Parameters:
mu_opd0: float

Instrumental OPD around which random initial guesses are created.

sig_opd0: float

Instrumental OPD around which random initial guesses are created.

na0: float

Instrumental OPD around which random initial guesses are created.

bounds_mu: 2-tuple

Lower and upper bounds between which the random values of mu_opd must be.

bounds_sig: 2-tuple

Lower and upper bounds between which the random values of sig_opd must be.

bounds_na: 2-tuple

Lower and upper bounds between which the random values of na must be.

Returns:
mu_opd: float

New initial guess for mu_opd.

sig_opd: float

New initial guess for sig_opd.

na: float

New initial guess for na.

glint_fitting_gpu6.gaussian(x, a, x0, sig)

Computes a Gaussian curve.

Parameters:

x: values where the curve is estimated.

a: amplitude of the Gaussian.

x0: location of the Gaussian.

sig: scale of the Gaussian.

Returns:

Gaussian curve.