Skip to content

DDT Yorick layout

Jump to bottom Edit New page
cmsaunders edited this page Jul 11, 2014 · 19 revisions

The main entry point is the Yorick script run_ddt.i

It takes as an argument the filename of another Yorick script. This script contains all the inputs and output filenames and optional parameters. run_ddt.i does the following:

  • includes ddt_startup.i. This includes all the necessary modules (Yorick and DDT-specific), like a lot of imports in Python.

  • For job=0, it includes run_ddt_galaxy_subtraction_all_regul_variable_all_epoch.i. For other job numbers, it includes different scripts.

The following details the script run_ddt_galaxy_subtraction_all_regul_variable_all_epoch.i

1. Set-up

ddt = ddt_setup_ddt(option_include_file) [defined in ddt_setup.i]

does a few scalar parameter checks

cube I/O

ddt_data = ddt_read_dataset(in_cube, ... ) [defined in ddt_io.i]

in_cube is a list of cube filenames

Calls ddt_read_datacube() repeatedly. Each call:

  • calls ddt_get_axis() - kinda gets WCS for single axis
  • calls ddt_read_array() - simple read with BSCALE/BZERO FITS handling
  • returns hash table with members x, y, lambda, data, weight. x, y apparently not used.

Header keys

Reads select keywords from a single FITS header. Uses the final ref if there is one, otherwise, just uses the first FITS file in the in_cube list.

PSF creation

Does different things depending on PARAM_PSF_TYPE variable. For "ES-PSF":

  • calls ddt_psf_from_GS(PARAM_PSF_ES, ...) (input is 2-d array from option file) [defined in ddt_psf_toolbox.i]
    • calls ddt_ES_param_from_GS() to calculate ellipticity and alpha.
    • calls ddt_gaussian_moffat_psf() with inputs ellipticity and alpha.
      • calls gaussian_moffat_psf_2d()
        • calls gaussian_moffat_psf()

In the above functions, ellipticity and alpha are both 2-d arrays where the first axis corresponds to wavelength and the second axis corresponds to time.

For "G-PSF" it loads the PSF from FITS files... which is strange because in the example input files, there YAML filenames where the FITS files should be.

Loading the model

  • calls ddt_setup_model() with different parameters depending on whether galaxy_file is defined in the input setup script. Basically, this just initializes some empty arrays of the right size.

ADR

  • calls ddt_calculate_adr() [defined in ddt_adr.i]
  • In turn, calls other functions in ddt_adr.i.
  • Key return values are paralactic angle of each observation and differential refraction as a function of wavelength for each observation.

Other stuff

  • calls ddt_setup_fft, ddt : Initializes an FFT (not clear to what extent: arrays allocated at this point?)
  • calls ddt_setup_R, ddt : sets up supernova offsets, "ranges" and an operator that returns a subset of an input array (as a copy) based on ranges.
  • ... more to document ...

2. "Prepare"

Run ddt_prepare_ddt() [defined in ddt_run.i]

  • calls ddt_precondition : divides data array (and other arrays) by the mean
  • calls ddt_sky_guess_all(ddt.ddt_data) [defined in ddt_data_toolbox.i]: guesses the sky based on all the data.
    • for each time slice, calls ddt_sky_guess
    • returns 2-d array of shape (n_lambda, n_t): spatially constant sky estimation.
  • sets ddt.ddt_data.guess_sky to the above result.
  • sets ddt.ddt_model.final_ref_sky to the slice of the above result that is for the final ref.
  • calculate the difference between the data and sky guess for the final ref
  • set ddt.galaxy_min_limit to the global minimum of that difference (scalar)
  • set ddt.ddt_data.guess_galaxy to the spatial average of that difference (1-d array in wavelength).
  • set ddt.data_ref to that difference.
  • set ddt.weight_ref to the final ref's 3-d slice of the main 4-d weight array.

3. Set up regularization

  • calls ddt_setup_regularization_galaxy_xy_normalized_1 (spatial component)
    • turned this into an object 'regul_galaxy_xy'
    • sets regularization function
    • sets regul_galaxy_xy.q to inverse spectrum
    • sets regul_galaxy_xy.mu_estim to 1 (originally was number of data points/estimate of data variation)
    • sets several other variables to 1 or -1.
    • defines call function that uses the class attributes and the set regularization fn to get a regularization value. Also changes some 'grd' argument. and then
  • calls ddt_setup_regularization_galaxy_lambda_normalized_3 (wavelength component)
    • Also turned this into an object
    • Similar to above, but doesn't set regularization fn, no extra variables set to -1 or 1.
    • Also just sets mu_estim to 1. Comment in DDT says calculation not implemented yet.
    • Also has call function, sets regularization value, changes some 'grd' argument.

These take as inputs:

  • The data and weight cubes of the final ref
  • A prior (mu_galaxy_xy_prior or mu_galaxy_lambda_prior)
  • The "sky guess" for the final ref (ddt.ddt_data.guess_sky)

and set regul_galaxy_xy and regul_galaxy_lambda as attributes of the model (ddt.ddt_model).

  • Finally, initializes ddt_model.galaxy to an 3-d array of zeros (3-d = spatial and wavelength components). This is actually already done in ddt_setup_model (in the case no galaxy file is provided).

4. Fit

  • set ddt.i_fit to the index of the final ref.

  • set ddt.fit_final_ref_sky to false.

  • run ddt_fit_model_all_epoch() [defined in ddt_optim_toolbox.i]

    This operates on the main ddt object and takes a few options, including the "minimum limit" for the galaxy: ddt.galaxy_min_limit.

    • run op_mnb(penalty, x, ...) from OptimPack
      • penalty: a function, in this case ddt_penalty_g_all_epoch
      • x : starting point (set to ddt.ddt_model.galaxy)
      • xmin is set to ddt.galaxy_min_limit
    • set ddt.ddt_model.galaxy_old to the galaxy before the fit.
    • set ddt.ddt_model.galaxy to the result of the optimization.

  • h_set, ddt.ddt_model, galaxy = x_new;

  • ddt_setup_G, ddt : set the ddt.G to a list of functions (1 per time slice)

  • calls sn_sky = ddt_extract_eta_sn_sky_all(ddt, ...): calculates the optimal SN and Sky in the chi^2 sense for given PSF and galaxy, including a possible offset of the supernova

  • Registration part (maybe this should be separate section):

    • Sets some parameters, sets the weighting for the alignment refinement

      • includes ddt_calc_residual (only uses attribute 'm')
      • the relevant part of this involves making a model cube out of the (fits of) the galaxy, sn, and sky, then doing R(H(model cube)), which I think convolves the model with the psf and then corrects the position.

    • Register the galaxy in the other final refs:

      • calls galaxy_position_fit = ddt_sn_galaxy_registration(ddt, ...)
      • the only part used of what this function returns is the galaxy offset a, which is x,y coordinates.
      • offset is found by using op_nllsq to I think minimize (data-model)*weight**2 as a function of a.
      • op_nllsq is a non-linear least squares fit in OptimPack-3.0/yorick/OptimPack3.i
      • recalculate is set to 1, so this step also redoes ddt_extract_eta_sn_sky and updates sn and sky, but it looks like the option is turned off that updates the whole DDT object.
      • if the offset is greater than n_spaxel_max_galaxy=3, the final ref is thrown out.

    • calls ddr_adr_sn = ddt_setup_new_pointing_all

    • calls ddt_extract_eta_sn_sky_all which runs ddt_extract_eta_sn_sky for all final refs and updates the DDT object.

Add a custom footer
Add a custom sidebar
Clone this wiki locally