diff --git a/.gitignore b/.gitignore index 1f67347c..4a387ecb 100644 --- a/.gitignore +++ b/.gitignore @@ -16,4 +16,5 @@ dist/ wiki/ doc/_build/ .vscode -env/* \ No newline at end of file +env/* +.DS_Store \ No newline at end of file diff --git a/AegeanTools/AeRes.py b/AegeanTools/AeRes.py index 395fdf7c..4641ee5e 100644 --- a/AegeanTools/AeRes.py +++ b/AegeanTools/AeRes.py @@ -12,16 +12,25 @@ from AegeanTools import catalogs, fitting, wcs_helpers from AegeanTools.logging import logger, logging +from AegeanTools.models import ComponentSource, ComponentSource3D +from AegeanTools.fits_tools import load_image_band, write_fits __author__ = "Paul Hancock" FWHM2CC = 1 / (2 * np.sqrt(2 * np.log(2))) -def load_sources(filename, - ra_col='ra', dec_col='dec', - peak_col='peak_flux', - a_col='a', b_col='b', pa_col='pa'): +def load_sources( + filename, + ra_col="ra", + dec_col="dec", + peak_col="peak_flux", + a_col="a", + b_col="b", + pa_col="pa", + alpha_col="alpha", + nu0_col="nu0", +): """ Open a file, read contents, return a list of all the sources in that file. @@ -32,13 +41,14 @@ def load_sources(filename, ra_col, dec_col, peak_col, a_col, b_col, pa_col : str The column names for each of the parameters. - Default = ['ra', 'dec', 'peak_flux', 'a', 'b', 'pa'] + Default = ['ra', 'dec', 'peak_flux', 'a', 'b', 'pa', 'alpha' ,'nu0'] Return ------ catalog : [`class:AegeanTools.models.ComponentSource`, ...] A list of source components """ + srctype = ComponentSource table = catalogs.load_table(filename) required_cols = [ra_col, dec_col, peak_col, a_col, b_col, pa_col] # required_cols = ['ra','dec','peak_flux','a','b','pa'] @@ -51,11 +61,19 @@ def load_sources(filename, logger.error("Some required columns missing or mis-labeled") return None # rename the table columns - for old, new in zip([ra_col, dec_col, peak_col, a_col, b_col, pa_col], - ['ra', 'dec', 'peak_flux', 'a', 'b', 'pa']): + for old, new in zip( + [ra_col, dec_col, peak_col, a_col, b_col, pa_col], + ["ra", "dec", "peak_flux", "a", "b", "pa"], + ): table.rename_column(old, new) - catalog = catalogs.table_to_source_list(table) + # rename the alpha column if it exists + if alpha_col in table.colnames: + table.rename_column(alpha_col, "alpha") + table.rename_column(nu0_col, "nu0") + srctype = ComponentSource3D + + catalog = catalogs.table_to_source_list(table, src_type=srctype) logger.info("read {0} sources from {1}".format(len(catalog), filename)) return catalog @@ -69,7 +87,7 @@ def make_model(sources, shape, wcshelper, mask=False, frac=None, sigma=4): sources : [`class:AegeanTools.models.ComponentSource`, ...] a list of sources - shape : [float, float] + shape : [int, int, int] the shape of the input (and output) image wcshelper : 'class:AegeanTools.wcs_helpers.WCSHelper' @@ -96,79 +114,98 @@ def make_model(sources, shape, wcshelper, mask=False, frac=None, sigma=4): factor = 5 i_count = 0 - for src in sources: - xo, yo, sx, sy, theta = wcshelper.sky2pix_ellipse([src.ra, src.dec], - src.a/3600, - src.b/3600, src.pa) - phi = np.radians(theta) - - # skip sources that have a center that is outside of the image - if not 0 < xo < shape[0]: - logger.debug("source {0} is not within image".format(src.island)) - continue - if not 0 < yo < shape[1]: - logger.debug("source {0} is not within image".format(src.island)) - continue - - # pixels over which this model is calculated - xoff = factor*(abs(sx*np.cos(phi)) + abs(sy*np.sin(phi))) - xmin = xo - xoff - xmax = xo + xoff - - yoff = factor*(abs(sx*np.sin(phi)) + abs(sy*np.cos(phi))) - ymin = yo - yoff - ymax = yo + yoff - - # clip to the image size - ymin = max(np.floor(ymin), 0) - ymax = min(np.ceil(ymax), shape[1]) - - xmin = max(np.floor(xmin), 0) - xmax = min(np.ceil(xmax), shape[0]) - - if not np.all(np.isfinite([ymin, ymax, xmin, xmax])): - continue - - if logger.isEnabledFor(logging.DEBUG): # pragma: no cover - logger.debug("Source ({0},{1})".format(src.island, src.source)) - logger.debug(" xo, yo: {0} {1}".format(xo, yo)) - logger.debug(" sx, sy: {0} {1}".format(sx, sy)) - logger.debug(" theta, phi: {0} {1}".format(theta, phi)) - logger.debug(" xoff, yoff: {0} {1}".format(xoff, yoff)) - logger.debug(" xmin, xmax, ymin, ymax: {0}:{1} {2}:{3}".format( - xmin, xmax, ymin, ymax)) - logger.debug(" flux, sx, sy: {0} {1} {2}".format( - src.peak_flux, sx, sy)) - - # positions for which we want to make the model - x, y = np.mgrid[int(xmin):int(xmax), int(ymin):int(ymax)] - x = x.ravel() - y = y.ravel() - - # TODO: understand why xo/yo -1 is needed - model = fitting.elliptical_gaussian(x, y, src.peak_flux, xo-1, yo-1, - sx*FWHM2CC, sy*FWHM2CC, theta) - - # Mask the output image if requested - if mask: - if frac is not None: - indices = np.where(model >= (frac*src.peak_flux)) + for s in range(shape[0]): + for src in sources: + xo, yo, sx, sy, theta = wcshelper.sky2pix_ellipse( + [src.ra, src.dec], src.a / 3600, src.b / 3600, src.pa + ) + phi = np.radians(theta) + + # skip sources that have a center that is outside of the image + if not 0 < xo < shape[1]: + logger.debug("source {0} is not within image".format(src.island)) + continue + if not 0 < yo < shape[2]: + logger.debug("source {0} is not within image".format(src.island)) + continue + + # pixels over which this model is calculated + xoff = factor * (abs(sx * np.cos(phi)) + abs(sy * np.sin(phi))) + xmin = xo - xoff + xmax = xo + xoff + + yoff = factor * (abs(sx * np.sin(phi)) + abs(sy * np.cos(phi))) + ymin = yo - yoff + ymax = yo + yoff + + # clip to the image size + ymin = max(np.floor(ymin), 0) + ymax = min(np.ceil(ymax), shape[2]) + + xmin = max(np.floor(xmin), 0) + xmax = min(np.ceil(xmax), shape[1]) + + if not np.all(np.isfinite([ymin, ymax, xmin, xmax])): + continue + + # Compute the peak flux at this frequency + peak_flux = src.peak_flux + if hasattr(src, "nu0") and wcshelper.has_freq: + peak_flux *= (wcshelper.pix2freq(s) / src.nu0) ** src.alpha + + if logger.isEnabledFor(logging.DEBUG): # pragma: no cover + logger.debug("Source ({0},{1})".format(src.island, src.source)) + logger.debug(" xo, yo: {0} {1}".format(xo, yo)) + logger.debug(" sx, sy: {0} {1}".format(sx, sy)) + logger.debug(" theta, phi: {0} {1}".format(theta, phi)) + logger.debug(" xoff, yoff: {0} {1}".format(xoff, yoff)) + logger.debug( + " xmin, xmax, ymin, ymax: {0}:{1} {2}:{3}".format( + xmin, xmax, ymin, ymax + ) + ) + logger.debug(" flux, sx, sy: {0} {1} {2}".format(peak_flux, sx, sy)) + + # positions for which we want to make the model + x, y = np.mgrid[int(xmin) : int(xmax), int(ymin) : int(ymax)] + x = x.ravel() + y = y.ravel() + + # TODO: understand why xo/yo -1 is needed + model = fitting.elliptical_gaussian( + x, y, src.peak_flux, xo - 1, yo - 1, sx * FWHM2CC, sy * FWHM2CC, theta + ) + + # Mask the output image if requested + if mask: + if frac is not None: + indices = np.where(model >= (frac * peak_flux)) + else: + indices = np.where(model >= (sigma * src.local_rms)) + # somehow m[x,y][indices] = np.nan doesn't assign any values + # so we have to do the more complicated + # m[x[indices],y[indices]] = np.nan + m[s, x[indices], y[indices]] = np.nan else: - indices = np.where(model >= (sigma*src.local_rms)) - # somehow m[x,y][indices] = np.nan doesn't assign any values - # so we have to do the more complicated - # m[x[indices],y[indices]] = np.nan - m[x[indices], y[indices]] = np.nan - else: - m[x, y] += model - i_count += 1 + m[s, x, y] += model + + if s == 0: + i_count += 1 logger.info("modeled {0} sources".format(i_count)) return m -def make_residual(fitsfile, catalog, rfile, mfile=None, - add=False, mask=False, frac=None, sigma=4, - colmap=None): +def make_residual( + fitsfile, + catalog, + rfile, + mfile=None, + add=False, + mask=False, + frac=None, + sigma=4, + colmap=None, +): """ Take an input image and catalogue, make a model of the catalogue, and then add/subtract or mask the input image. Saving the residual and (optionally) @@ -205,7 +242,7 @@ def make_residual(fitsfile, catalog, rfile, mfile=None, colmap : dict A mapping of column names. Default is: {'ra_col':'ra', 'dec_col':'dec', - 'peak_col':'peak_flux', 'a_col':'a', 'b_col':'b', 'pa_col':'pa} + 'peak_col':'peak_flux', 'a_col':'a', 'b_col':'b', 'pa_col':'pa, 'alpha_col':'alpha'} Return ------ @@ -219,12 +256,8 @@ def make_residual(fitsfile, catalog, rfile, mfile=None, if source_list is None: return None - # force two axes so that we dump redundant stokes/freq axes if they are - # present. - hdulist = fits.open(fitsfile, naxis=2) - # ignore dimensions of length 1 - data = np.squeeze(hdulist[0].data) - header = hdulist[0].header + + data, header = load_image_band(fitsfile, as_cube=True) wcshelper = wcs_helpers.WCSHelper.from_header(header) @@ -235,11 +268,8 @@ def make_residual(fitsfile, catalog, rfile, mfile=None, else: residual = data - model - hdulist[0].data = residual - hdulist.writeto(rfile, overwrite=True) - logger.info("wrote residual to {0}".format(rfile)) + write_fits(residual, header, rfile) if mfile is not None: - hdulist[0].data = model - hdulist.writeto(mfile, overwrite=True) - logger.info("wrote model to {0}".format(mfile)) + write_fits(model, header, mfile) + return diff --git a/AegeanTools/BANE.py b/AegeanTools/BANE.py index 74d76428..ef4324d5 100755 --- a/AegeanTools/BANE.py +++ b/AegeanTools/BANE.py @@ -104,8 +104,6 @@ def sigmaclip(arr, lo, hi, reps=10): std = np.std(clipped) mean = np.mean(clipped) prev_valid = curr_valid - # else: # disable logging so that numba can be fast - # logger.debug("No stopping criteria was reached after {0} cycles".format(count)) return mean, std @@ -130,11 +128,11 @@ def _sf2(args): except Exception as e: import traceback - logger.warn(e) + logger.warning(e) raise Exception("".join(traceback.format_exception(*sys.exc_info()))) -def sigma_filter(filename, region, step_size, box_size, shape, domask, cube_index): +def sigma_filter(filename, region, step_size, box_size, shape, domask, cube_index=None): """ Calculate the background and rms for a sub region of an image. The results are written to shared memory - irms and ibkg. @@ -154,14 +152,15 @@ def sigma_filter(filename, region, step_size, box_size, shape, domask, cube_inde box_size : (int, int) The size of the box over which the filter is applied (each step). - shape : tuple - The shape of the fits image + shape : (int, int, int) + The shape of the fits cube (1,y,x) for an image. domask : bool If true then copy the data mask to the output. - cube_index : int - The index into the 3rd dimension (if present) + cube_index : int or None + The index into the 3rd dimension (frequency?) to process. + Default = None => process all. Returns ------- @@ -178,29 +177,38 @@ def sigma_filter(filename, region, step_size, box_size, shape, domask, cube_inde # cut out the region of interest plus 1/2 the box size # and clip to the image size data_row_min = max(0, ymin - box_size[0] // 2) - data_row_max = min(shape[0], ymax + box_size[0] // 2) + data_row_max = min(shape[1], ymax + box_size[0] // 2) # Figure out how many axes are in the datafile NAXIS = fits.getheader(filename)["NAXIS"] + sz = slice(None) + sy = slice(data_row_min, data_row_max) + sx = slice(None) + slices = range(shape[0]) + # if the cube_index is not none, then only + # load/process one slice of the cube + if cube_index is not None: + sz = slice(cube_index, cube_index + 1) + slices = [0] + # For some reason we can't memmap a file with BSCALE not 1.0 # so we ignore it now and scale it later with fits.open(filename, memmap=True, do_not_scale_image_data=True) as a: if NAXIS == 2: - data = a[0].section[data_row_min:data_row_max, 0 : shape[1]] + data = a[0].section[sy, sx] + # ensure that we always end up with a 3d image + data = data[None, :, :] elif NAXIS == 3: - data = np.squeeze( - a[0].section[cube_index, data_row_min:data_row_max, 0 : shape[1]] - ) + data = a[0].section[sz, sy, sx] elif NAXIS == 4: - data = np.squeeze( - a[0].section[0, cube_index, data_row_min:data_row_max, 0 : shape[1]] - ) + data = a[0].section[0, sz, sy, sx] else: - logger.error("Too many NAXIS for me {0}".format(NAXIS)) + logger.error(f"Too many NAXIS for me {NAXIS}") logger.error("fix your file to be more sane") raise Exception("Too many NAXIS") + logger.debug(f"loaded data shape {data.shape}") # Manually scale the data if BSCALE is not 1.0 header = fits.getheader(filename) if "BSCALE" in header: @@ -209,10 +217,8 @@ def sigma_filter(filename, region, step_size, box_size, shape, domask, cube_inde # force float64 for consistency data = data.astype(np.float64) - # row_len = shape[1] - - logger.debug("data size is {0}".format(data.shape)) - logger.debug("data format is {0}".format(data.dtype)) + logger.debug(f"data size is {data.shape}") + logger.debug(f"data format is {data.dtype}") def box(r, c): """ @@ -220,30 +226,17 @@ def box(r, c): with size = box_size """ r_min = max(0, r - box_size[0] // 2) - r_max = min(data.shape[0] - 1, r + box_size[0] // 2) + r_max = min(data.shape[1] - 1, r + box_size[0] // 2) c_min = max(0, c - box_size[1] // 2) - c_max = min(data.shape[1] - 1, c + box_size[1] // 2) + c_max = min(data.shape[2] - 1, c + box_size[1] // 2) return r_min, r_max, c_min, c_max # set up a grid of rows/cols at which we will compute the bkg/rms rows = list(range(ymin - data_row_min, ymax - data_row_min, step_size[0])) rows.append(ymax - data_row_min) - cols = list(range(0, shape[1], step_size[1])) - cols.append(shape[1]) - - # store the computed bkg/rms in this smaller array - vals = np.zeros(shape=(len(rows), len(cols))) - for i, row in enumerate(rows): - for j, col in enumerate(cols): - r_min, r_max, c_min, c_max = box(row, col) - new = data[r_min:r_max, c_min:c_max] - new = np.ravel(new) - bkg, _ = sigmaclip(new, 3, 3) - vals[i, j] = bkg - - # indices of all the pixels within our region - gr, gc = np.mgrid[ymin - data_row_min : ymax - data_row_min, 0 : shape[1]] + cols = list(range(0, shape[2], step_size[1])) + cols.append(shape[2]) # Find the shared memory and create a numpy array interface ibkg_shm = SharedMemory(name=f"ibkg_{memory_id}", create=False) @@ -251,55 +244,84 @@ def box(r, c): irms_shm = SharedMemory(name=f"irms_{memory_id}", create=False) irms = np.ndarray(shape, dtype=np.float64, buffer=irms_shm.buf) - logger.debug("Interpolating bkg to sharemem") - ifunc = RegularGridInterpolator((rows, cols), vals) - interp_bkg = np.array(ifunc((gr, gc)), dtype=np.float64) - ibkg[ymin:ymax, :] = interp_bkg - del ifunc, interp_bkg - logger.debug(" ... done writing bkg") - - # wait for all to complete - i = barrier.wait() - if i == 0: - barrier.reset() - - logger.debug("background subtraction") - data[0 + ymin - data_row_min : data.shape[0] - (data_row_max - ymax), :] -= ibkg[ - ymin:ymax, : - ] - logger.debug(".. done ") - - # reset/recycle the vals array - vals[:] = 0 - - for i, row in enumerate(rows): - for j, col in enumerate(cols): - r_min, r_max, c_min, c_max = box(row, col) - new = data[r_min:r_max, c_min:c_max] - new = np.ravel(new) - _, rms = sigmaclip(new, 3, 3) - vals[i, j] = rms - - logger.debug("Interpolating rms to sharemem") - ifunc = RegularGridInterpolator((rows, cols), vals) - interp_rms = np.array(ifunc((gr, gc)), dtype=np.float64) - irms[ymin:ymax, :] = interp_rms - del ifunc, interp_rms - logger.debug(" .. done writing rms") - - if domask: + for k in slices: + logger.debug(f"Working on slice {k}") + # store the computed bkg/rms in this smaller array + vals = np.zeros(shape=(len(rows), len(cols))) + + # loop over + for i, row in enumerate(rows): + for j, col in enumerate(cols): + r_min, r_max, c_min, c_max = box(row, col) + new = data[k, r_min:r_max, c_min:c_max] + new = np.ravel(new) + bkg, _ = sigmaclip(new, 3, 3) + vals[i, j] = bkg + + # indices of all the pixels within our region + gr, gc = np.mgrid[ymin - data_row_min : ymax - data_row_min, 0 : shape[2]] + logger.debug(f"gr has shape {gr.shape}") + logger.debug("Interpolating bkg to sharemem") + ifunc = RegularGridInterpolator((rows, cols), vals) + interp_bkg = np.array(ifunc((gr, gc)), dtype=np.float64) + ibkg[k, ymin:ymax, :] = interp_bkg + del ifunc, interp_bkg + logger.debug(" ... done writing bkg") + # wait for all to complete i = barrier.wait() if i == 0: barrier.reset() - logger.debug("applying mask") - mask = ~np.isfinite( - data[0 + ymin - data_row_min : data.shape[0] - (data_row_max - ymax), :] + logger.debug("background subtraction") + logger.debug(f"data shape {data.shape}, ibkg shape {ibkg.shape}") + logger.debug( + f"k {k}, ymin {ymin}, ymax {ymax}, data_row_min {data_row_min}, data_row_max {data_row_max}" ) - ibkg[ymin:ymax, :][mask] = np.nan - irms[ymin:ymax, :][mask] = np.nan - logger.debug("... done applying mask") + logger.debug( + f"data slice is {0 + ymin - data_row_min}:{data.shape[1] - (data_row_max - ymax)}" + ) + data[ + k, 0 + ymin - data_row_min : data.shape[1] - (data_row_max - ymax), : + ] -= ibkg[k, ymin:ymax, :] + logger.debug(".. done ") + + # reset/recycle the vals array + vals[:] = 0 + + for i, row in enumerate(rows): + for j, col in enumerate(cols): + r_min, r_max, c_min, c_max = box(row, col) + new = data[k, r_min:r_max, c_min:c_max] + new = np.ravel(new) + _, rms = sigmaclip(new, 3, 3) + vals[i, j] = rms + + logger.debug("Interpolating rms to sharemem") + ifunc = RegularGridInterpolator((rows, cols), vals) + interp_rms = np.array(ifunc((gr, gc)), dtype=np.float64) + irms[k, ymin:ymax, :] = interp_rms + del ifunc, interp_rms + logger.debug(" .. done writing rms") + + if domask: + # wait for all to complete + i = barrier.wait() + if i == 0: + barrier.reset() + + logger.debug("applying mask") + mask = ~np.isfinite( + data[ + k, # TODO: CHECK HERE + 0 + ymin - data_row_min : data.shape[1] - (data_row_max - ymax), + :, + ] + ) + ibkg[k, ymin:ymax, :][mask] = np.nan + irms[k, ymin:ymax, :][mask] = np.nan + logger.debug("... done applying mask") + logger.debug( "rows {0}-{1} finished at {2}".format( ymin, ymax, strftime("%Y-%m-%d %H:%M:%S", gmtime()) @@ -309,7 +331,14 @@ def box(r, c): def filter_mc_sharemem( - filename, step_size, box_size, cores, shape, nslice=None, domask=True, cube_index=0 + filename, + step_size, + box_size, + cores, + shape, + nslice=None, + domask=True, + cube_index=None, ): """ Calculate the background and noise images corresponding to the input file. @@ -330,18 +359,20 @@ def filter_mc_sharemem( cores : int Number of cores to use. If None then use all available. + shape : (int, int) or (int, int, int) + The shape of the image or cube in the given file. + nslice : int The image will be divided into this many horizontal stripes for processing. Default = None = equal to cores - shape : (int, int) - The shape of the image in the given file. - domask : bool True(Default) = copy data mask to output. - cube_index : int - For 3d data use this index into the third dimension. Default = 0 + cube_index : int or None + For 3d data use this index into the third dimension. + Default = None -> 3d input gives 3d output. + Returns ------- @@ -354,7 +385,8 @@ def filter_mc_sharemem( if (nslice is None) or (cores == 1): nslice = cores - img_y, img_x = shape + # shape is (z,y,x) + img_y = shape[1] logger.info("using {0} cores".format(cores)) logger.info("using {0} stripes".format(nslice)) @@ -435,7 +467,6 @@ def filter_image( out_base, step_size=None, box_size=None, - twopass=False, # Deprecated cores=None, mask=True, compressed=False, @@ -461,10 +492,6 @@ def filter_image( box_size : (int,int) The size of the box in pixels - twopass : bool - Perform a second pass calculation to ensure that the noise is not - contaminated by the background. Default = False. DEPRECATED - cores : int Number of CPU corse to use. Default = all available @@ -480,29 +507,27 @@ def filter_image( Return a compressed version of the background/noise images. Default = False - cube_index : int + cube_index : int or None If the input data is 3d, then use this index for the 3rd dimension. - Default = None, use the first index. + Default = None, 3d inputs give 3d outputs. Returns ------- bkg, rms : `numpy.ndarray` The computed background and rms maps (not compressed) """ - # Use the first slice of the 3rd dimension if not specified - if cube_index is None: - cube_index = 0 header = fits.getheader(im_name) - shape = (header["NAXIS2"], header["NAXIS1"]) + shape = (1, header["NAXIS2"], header["NAXIS1"]) + if ("NAXIS3" in header) and cube_index is None: + shape = (header["NAXIS3"], header["NAXIS2"], header["NAXIS1"]) + naxis = header["NAXIS"] if naxis > 2: naxis3 = header["NAXIS3"] - if cube_index >= naxis3: + if cube_index is not None and cube_index >= naxis3: logger.error( - "3rd dimension has len {0} but index {1} was passed".format( - naxis3, cube_index - ) + f"3rd dimension has len {naxis3} but index {cube_index} was passed" ) return None @@ -517,9 +542,7 @@ def filter_image( if not step_size[0] == step_size[1]: step_size = (min(step_size), min(step_size)) logger.info( - "Changing grid to be {0} so we can compress the output".format( - step_size - ) + f"Changing grid to be {step_size} so we can compress the output" ) logger.info("using grid_size {0}, box_size {1}".format(step_size, box_size)) diff --git a/AegeanTools/CLI/AeRes.py b/AegeanTools/CLI/AeRes.py index 35c264f8..13ff7ce8 100755 --- a/AegeanTools/CLI/AeRes.py +++ b/AegeanTools/CLI/AeRes.py @@ -6,8 +6,8 @@ from AegeanTools.logging import logger, logging __author__ = "Paul Hancock" -__version__ = "v0.2.7" -__date__ = "2025-01-14" +__version__ = "v0.3" +__date__ = "2025-02-03" def main(): @@ -92,6 +92,16 @@ def main(): group3.add_argument( "--pacol", dest="pa_col", default="pa", help="Position angle column name" ) + group3.add_argument( + "--alphacol", dest="alpha_col", default="alpha", help="Alpha column name" + ) + + group3.add_argument( + "--nu0col", + dest="nu0_col", + default="nu0", + help="Reference frequency column name", + ) group4 = parser.add_argument_group("Extra options") group4.add_argument( @@ -136,6 +146,8 @@ def main(): "a_col": options.a_col, "b_col": options.b_col, "pa_col": options.pa_col, + "alpha_col": options.alpha_col, + "nu0_col": options.nu0_col, } make_residual( diff --git a/AegeanTools/CLI/BANE.py b/AegeanTools/CLI/BANE.py index 3eb0a70b..18daa301 100755 --- a/AegeanTools/CLI/BANE.py +++ b/AegeanTools/CLI/BANE.py @@ -52,7 +52,7 @@ def main(): "--slice", dest="cube_index", type=int, - default=0, + default=None, help="If the input data is a cube, then this slice " "will determine the array index of the image " "which will be processed by BANE", @@ -96,6 +96,7 @@ def main(): cores=multiprocessing.cpu_count(), stripes=multiprocessing.cpu_count() - 1, debug=False, + cube_index=None, ) options = parser.parse_args() diff --git a/AegeanTools/CLI/aegean.py b/AegeanTools/CLI/aegean.py index 12b9d976..6177006a 100755 --- a/AegeanTools/CLI/aegean.py +++ b/AegeanTools/CLI/aegean.py @@ -8,19 +8,48 @@ import sys import astropy +from astropy.table import vstack import lmfit import numpy as np import scipy from AegeanTools import __citation__, __date__, __version__ -from AegeanTools.catalogs import check_table_formats, save_catalog, show_formats +from AegeanTools.catalogs import check_table_formats, save_catalog, show_formats, load_table, write_table from AegeanTools.logging import logger, logging from AegeanTools.source_finder import get_aux_files from AegeanTools.wcs_helpers import Beam +from AegeanTools.mpi import MPI_AVAIL, MPI +from AegeanTools.exceptions import AegeanSuffixError header = """#Aegean version {0} # on dataset: {1}""" +def addSuffix(file, suffix): + """ + A function to add a specified suffix before the extension. + + parameters + ---------- + file: str + The current name of the file + + suffix: str or int + The desired suffix to be inserted before the extension + """ + if isinstance(suffix, int): + base, ext = os.path.splitext(file) + base += f"_{suffix:02d}" + fname = base + ext + elif isinstance(suffix, str): + if suffix[0] == "_": + suffix = suffix[1:] + base, ext = os.path.splitext(file) + base += f"_{suffix}" + fname = base + ext + else: + raise AegeanSuffixError(f"This suffix type is not support: {suffix}") + + return fname def check_projection(filename, options, log=logger): """ @@ -158,6 +187,13 @@ def main(): "region, and psf files using the input filename " "as a hint. [default: don't do this]", ) + group2.add_argument( + "--3d", + dest="threeD", + action="store_true", + default=False, + help="Treat the input image as a 3D cube. " + ) # Output group3 = parser.add_argument_group("Output Options") @@ -381,11 +417,11 @@ def main(): print(__citation__) return 0 - import AegeanTools - from AegeanTools.source_finder import SourceFinder + import AegeanTools #! This import is not at the top + from AegeanTools.source_finder import SourceFinder #! This import is not at the top # source finding object - sf = SourceFinder() + sf = SourceFinder() #! <--- This is the source finder object if options.table_formats: show_formats() @@ -585,28 +621,53 @@ def main(): logger.warning( "--island requested but not yet supported for " "priorized fitting" ) - sf.priorized_fit_islands( + if options.threeD: + sf.priorized_fit_islands_3D( filename, catalogue=options.input, - hdu_index=options.hdu_index, - rms=options.rms, - bkg=options.bkg, + # hdu_index=options.hdu_index, + # rms=options.rms, + # bkg=options.bkg, outfile=options.outfile, bkgin=options.backgroundimg, rmsin=options.noiseimg, - beam=options.beam, - imgpsf=options.imgpsf, - catpsf=options.catpsf, + # beam=options.beam, + # imgpsf=options.imgpsf, + # catpsf=options.catpsf, stage=options.priorized, ratio=options.ratio, outerclip=options.outerclip, cores=options.cores, doregroup=options.regroup, docov=options.docov, - cube_index=options.slice, + # cube_index=options.slice, progress=options.progress, regroup_eps=options.regroup_eps, - ) + ) + + else: + sf.priorized_fit_islands( + filename, + catalogue=options.input, + hdu_index=options.hdu_index, + rms=options.rms, + bkg=options.bkg, + outfile=options.outfile, + bkgin=options.backgroundimg, + rmsin=options.noiseimg, + beam=options.beam, + imgpsf=options.imgpsf, + catpsf=options.catpsf, + stage=options.priorized, + ratio=options.ratio, + outerclip=options.outerclip, + cores=options.cores, + doregroup=options.regroup, + docov=options.docov, + cube_index=options.slice, + progress=options.progress, + regroup_eps=options.regroup_eps, + ) sources = sf.sources @@ -618,6 +679,34 @@ def main(): "FITSFILE": filename, "RUN-AS": invocation_string, } + + # collect catalogues and clean up + if MPI_AVAIL: + logger.info("MPI is available") + catalog_list = [] for t in options.tables.split(","): + final_file_name = t + if MPI_AVAIL: + t = addSuffix(t,MPI.COMM_WORLD.Get_rank()) save_catalog(t, sources, prefix=options.column_prefix, meta=meta) + + if MPI_AVAIL: + MPI.COMM_WORLD.Barrier() + if MPI.COMM_WORLD.Get_rank() == 0: + for t in options.tables.split(","): + flist = [] + final_file_name = addSuffix(t,"comp") + for n in range(0, MPI.COMM_WORLD.Get_size()): + base = addSuffix(t,n) + base = addSuffix(base, "comp") + flist.append(base) + final_table = load_table(flist[0]) + for f in flist[1:]: + table = load_table(f) + final_table = vstack([final_table, table]) + os.remove(f) + + write_table(final_table, final_file_name) + logger.info(f"Wrote table name: {final_file_name}") + return 0 diff --git a/AegeanTools/catalogs.py b/AegeanTools/catalogs.py index 92994aef..0fc1767a 100755 --- a/AegeanTools/catalogs.py +++ b/AegeanTools/catalogs.py @@ -198,7 +198,7 @@ def load_catalog(filename): """ supported = get_table_formats() - fmt = os.path.splitext(filename)[-1][1:].lower() # extension sans '.' + fmt = os.path.splitext(filename)[-1][1:].lower() # extension sans '.' #!There is a pythonic built in function to do this, it is in the os.system library if fmt in ["csv", "tab", "tex"] and fmt in supported: logger.info("Reading file {0}".format(filename)) @@ -247,7 +247,7 @@ def load_table(filename): """ supported = get_table_formats() - fmt = os.path.splitext(filename)[-1][1:].lower() # extension sans '.' + fmt = os.path.splitext(filename)[-1][1:].lower() # extension sans '.' #!Ditto here look up ^ if fmt in ["csv", "tab", "tex"] and fmt in supported: logger.info("Reading file {0}".format(filename)) @@ -300,7 +300,8 @@ def table_to_source_list(table, src_type=ComponentSource): A single table must have consistent source types given by src_type. src_type should be one of :class:`AegeanTools.models.ComponentSource`, :class:`AegeanTools.models.SimpleSource`, or - :class:`AegeanTools.models.IslandSource`. + :class:`AegeanTools.models.IslandSource` or + :class:`AegeanTools.models.ComponentSource3D`. Parameters @@ -311,7 +312,8 @@ def table_to_source_list(table, src_type=ComponentSource): src_type : class Sources must be of type :class:`AegeanTools.models.ComponentSource`, :class:`AegeanTools.models.SimpleSource`, or - :class:`AegeanTools.models.IslandSource`. + :class:`AegeanTools.models.IslandSource`or + :class:`AegeanTools.models.ComponentSource3D`. Returns ------- @@ -357,6 +359,7 @@ def write_catalog(filename, catalog, fmt=None, meta=None, prefix=None): catalog : list A list of source objects. Sources must be of type :class:`AegeanTools.models.ComponentSource`, + :class:`AegeanTools.models.ComponentSource3D`, :class:`AegeanTools.models.SimpleSource`, or :class:`AegeanTools.models.IslandSource`. diff --git a/AegeanTools/exceptions.py b/AegeanTools/exceptions.py index 90235909..b1bc2368 100644 --- a/AegeanTools/exceptions.py +++ b/AegeanTools/exceptions.py @@ -10,3 +10,9 @@ class AegeanNaNModelError(AegeanError): be NaNs in the data. It is a little misleading. """ + +class AegeanSuffixError(AegeanError): + """ + This is a child class of AegeanError that is raised when an incorrect suffix is provided + """ + pass \ No newline at end of file diff --git a/AegeanTools/fits_tools.py b/AegeanTools/fits_tools.py index 2b760a6c..039a9b39 100755 --- a/AegeanTools/fits_tools.py +++ b/AegeanTools/fits_tools.py @@ -251,7 +251,7 @@ def write_fits(data, header, file_name): return -def load_image_band(filename, band=(0, 1), hdu_index=0, cube_index=0): +def load_image_band(filename, band=(0, 1), hdu_index=0, cube_index=0, as_cube=False): """ Load a subset of an image from a given filename. The subset is controlled using the band, which is (this band, total bands) @@ -264,6 +264,16 @@ def load_image_band(filename, band=(0, 1), hdu_index=0, cube_index=0): (this band, total bands) Default (0,1) + hdu_index : int + The index of the HDU to load from the fits file. + + cube_index : int + The index of the cube to load from the fits file. + + as_cube : boolean + This is a flag that determines whether the data to be processed is a cube + or a frequency slice, if the data is 3 dimensional and as_cube is True, then a 3-dimensional array will be returned. + returns ------- data, header : :class:`numpy.ndarray`, :class:`astropy.io.fits.header.Header` @@ -296,14 +306,22 @@ def load_image_band(filename, band=(0, 1), hdu_index=0, cube_index=0): with fits.open(filename, memmap=True, do_not_scale_image_data=True) as a: if NAXIS == 2: data = a[hdu_index].section[row_min:row_max, 0 : header["NAXIS1"]] + if as_cube: + data = np.expand_dims(data, axis=0) elif NAXIS == 3: - data = a[hdu_index].section[ - cube_index, row_min:row_max, 0 : header["NAXIS1"] - ] + if not as_cube: + data = a[hdu_index].section[ + cube_index, row_min:row_max, 0 : header["NAXIS1"] + ] + else: + data = a[hdu_index].section[:, row_min:row_max, 0 : header["NAXIS1"]] elif NAXIS == 4: - data = a[hdu_index].section[ - 0, cube_index, row_min:row_max, 0 : header["NAXIS1"] - ] + if not as_cube: + data = a[hdu_index].section[ + 0, cube_index, row_min:row_max, 0 : header["NAXIS1"] + ] + else: + data = a[hdu_index].section[0, :, row_min:row_max, 0 : header["NAXIS1"]] else: raise Exception(f"Too many NAXIS: {NAXIS}>4") if "BSCALE" in header: diff --git a/AegeanTools/fitting.py b/AegeanTools/fitting.py index 83453679..fbfae550 100644 --- a/AegeanTools/fitting.py +++ b/AegeanTools/fitting.py @@ -268,7 +268,7 @@ def emp_jacobian(pars, x, y): eps = 1e-5 matrix = [] model = ntwodgaussian_lmfit(pars)(x, y) - for i in range(pars["components"].value): + for i in range(int(pars["components"].value)): prefix = "c{0}_".format(i) for p in ["amp", "xo", "yo", "sx", "sy", "theta"]: if pars[prefix + p].vary: @@ -335,544 +335,6 @@ def lmfit_jacobian(pars, x, y, errs=None, B=None, emp=False): return matrix -def hessian(pars, x, y): - """ - Create a hessian matrix corresponding to the source model 'pars' - Only parameters that vary will contribute to the hessian. - Thus there will be a total of nvar x nvar entries, each of which is a - len(x) x len(y) array. - - Parameters - ---------- - pars : lmfit.Parameters - The model - x, y : list - locations at which to evaluate the Hessian - - Returns - ------- - h : np.array - Hessian. Shape will be (nvar, nvar, len(x), len(y)) - - See Also - -------- - :func:`AegeanTools.fitting.emp_hessian` - """ - j = 0 # keeping track of the number of variable parameters - # total number of variable parameters - ntvar = np.sum([pars[k].vary for k in pars.keys() if k != "components"]) - # construct an empty matrix of the correct size - hmat = np.zeros((ntvar, ntvar, x.shape[0], x.shape[1])) - npvar = 0 - - for i in range(pars["components"].value): - prefix = "c{0}_".format(i) - amp = pars[prefix + "amp"].value - xo = pars[prefix + "xo"].value - yo = pars[prefix + "yo"].value - sx = pars[prefix + "sx"].value - sy = pars[prefix + "sy"].value - theta = pars[prefix + "theta"].value - - amp_var = pars[prefix + "amp"].vary - xo_var = pars[prefix + "xo"].vary - yo_var = pars[prefix + "yo"].vary - sx_var = pars[prefix + "sx"].vary - sy_var = pars[prefix + "sy"].vary - theta_var = pars[prefix + "theta"].vary - - # precomputed for speed - model = elliptical_gaussian(x, y, amp, xo, yo, sx, sy, theta) - sint = np.sin(np.radians(theta)) - sin2t = np.sin(np.radians(2 * theta)) - cost = np.cos(np.radians(theta)) - cos2t = np.cos(np.radians(2 * theta)) - sx2 = sx**2 - sy2 = sy**2 - xxo = x - xo - yyo = y - yo - xcos, ycos = xxo * cost, yyo * cost - xsin, ysin = xxo * sint, yyo * sint - - if amp_var: - k = npvar # second round of keeping track of variable params - # H(amp,amp)/G = 0 - hmat[j][k] = 0 - k += 1 - - if xo_var: - # H(amp,xo)/G = 1.0*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t))/(amp*sx**2*sy**2) - hmat[j][k] = (xsin - ycos) * sint / sy2 + (xcos + ysin) * cost / sx2 - hmat[j][k] *= model - k += 1 - - if yo_var: - # H(amp,yo)/G = 1.0*(-sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*cos(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t))/(amp*sx**2*sy**2) - hmat[j][k] = -(xsin - ycos) * cost / sy2 + (xcos + ysin) * sint / sx2 - hmat[j][k] *= model / amp - k += 1 - - if sx_var: - # H(amp,sx)/G = 1.0*((x - xo)*cos(t) + (y - yo)*sin(t))**2/(amp*sx**3) - hmat[j][k] = (xcos + ysin) ** 2 - hmat[j][k] *= model / (amp * sx**3) - k += 1 - - if sy_var: - # H(amp,sy) = 1.0*((x - xo)*sin(t) + (-y + yo)*cos(t))**2/(amp*sy**3) - hmat[j][k] = (xsin - ycos) ** 2 - hmat[j][k] *= model / (amp * sy**3) - k += 1 - - if theta_var: - # H(amp,t) = (-1.0*sx**2 + sy**2)*((x - xo)*sin(t) + (-y + yo)*cos(t))*((x - xo)*cos(t) + (y - yo)*sin(t))/(amp*sx**2*sy**2) - hmat[j][k] = (xsin - ycos) * (xcos + ysin) - hmat[j][k] *= sy2 - sx2 - hmat[j][k] *= model / (amp * sx2 * sy2) - # k += 1 - j += 1 - - if xo_var: - k = npvar - if amp_var: - # H(xo,amp)/G = H(amp,xo) - hmat[j][k] = hmat[k][j] - k += 1 - - # if xo_var: - # H(xo,xo)/G = 1.0*(-sx**2*sy**2*(sx**2*sin(t)**2 + sy**2*cos(t)**2) + (sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t))**2)/(sx**4*sy**4) - hmat[j][k] = -sx2 * sy2 * (sx2 * sint**2 + sy2 * cost**2) - hmat[j][k] += (sx2 * (xsin - ycos) * sint + sy2 * (xcos + ysin) * cost) ** 2 - hmat[j][k] *= model / (sx2**2 * sy2**2) - k += 1 - - if yo_var: - # H(xo,yo)/G = 1.0*(sx**2*sy**2*(sx**2 - sy**2)*sin(2*t)/2 - (sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*cos(t) - sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t)))/(sx**4*sy**4) - hmat[j][k] = sx2 * sy2 * (sx2 - sy2) * sin2t / 2 - hmat[j][k] -= ( - sx2 * (xsin - ycos) * sint + sy2 * (xcos + ysin) * cost - ) * (sx2 * (xsin - ycos) * cost - sy2 * (xcos + ysin) * sint) - hmat[j][k] *= model / (sx**4 * sy**4) - k += 1 - - if sx_var: - # H(xo,sx) = ((x - xo)*cos(t) + (y - yo)*sin(t))*(-2.0*sx**2*sy**2*cos(t) + 1.0*((x - xo)*cos(t) + (y - yo)*sin(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t)))/(sx**5*sy**2) - hmat[j][k] = xcos + ysin - hmat[j][k] *= -2 * sx2 * sy2 * cost + (xcos + ysin) * ( - sx2 * (xsin - ycos) * sint + sy2 * (xcos + ysin) * cost - ) - hmat[j][k] *= model / (sx**5 * sy2) - k += 1 - - if sy_var: - # H(xo,sy) = ((x - xo)*sin(t) + (-y + yo)*cos(t))*(-2.0*sx**2*sy**2*sin(t) + 1.0*((x - xo)*sin(t) + (-y + yo)*cos(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t)))/(sx2*sy**5) - hmat[j][k] = xsin - ycos - hmat[j][k] *= -2 * sx2 * sy2 * sint + (xsin - ycos) * ( - sx2 * (xsin - ycos) * sint + sy2 * (xcos + ysin) * cost - ) - hmat[j][k] *= model / (sx2 * sy**5) - k += 1 - - if theta_var: - # H(xo,t) = 1.0*(sx**2*sy**2*(sx**2 - sy**2)*(x*sin(2*t) - xo*sin(2*t) - y*cos(2*t) + yo*cos(2*t)) + (-sx**2 + 1.0*sy**2)*((x - xo)*sin(t) + (-y + yo)*cos(t))*((x - xo)*cos(t) + (y - yo)*sin(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*sin(t) + sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*cos(t)))/(sx**4*sy**4) - # second part - hmat[j][k] = (sy2 - sx2) * (xsin - ycos) * (xcos + ysin) - hmat[j][k] *= sx2 * (xsin - ycos) * sint + sy2 * (xcos + ysin) * cost - # first part - hmat[j][k] += sx2 * sy2 * (sx2 - sy2) * (xxo * sin2t - yyo * cos2t) - hmat[j][k] *= model / (sx**4 * sy**4) - # k += 1 - j += 1 - - if yo_var: - k = npvar - if amp_var: - # H(yo,amp)/G = H(amp,yo) - hmat[j][k] = hmat[0][2] - k += 1 - - if xo_var: - # H(yo,xo)/G = H(xo,yo)/G - hmat[j][k] = hmat[1][2] - k += 1 - - # if yo_var: - # H(yo,yo)/G = 1.0*(-sx**2*sy**2*(sx**2*cos(t)**2 + sy**2*sin(t)**2) + (sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*cos(t) - sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t))**2)/(sx**4*sy**4) - hmat[j][k] = ( - sx2 * (xsin - ycos) * cost - sy2 * (xcos + ysin) * sint - ) ** 2 / (sx2**2 * sy2**2) - hmat[j][k] -= cost**2 / sy2 + sint**2 / sx2 - hmat[j][k] *= model - k += 1 - - if sx_var: - # H(yo,sx)/G = -((x - xo)*cos(t) + (y - yo)*sin(t))*(2.0*sx**2*sy**2*sin(t) + 1.0*((x - xo)*cos(t) + (y - yo)*sin(t))*(sx**2*((x - xo)*sin(t) - (y - yo)*cos(t))*cos(t) - sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t)))/(sx**5*sy**2) - hmat[j][k] = -1 * (xcos + ysin) - hmat[j][k] *= 2 * sx2 * sy2 * sint + (xcos + ysin) * ( - sx2 * (xsin - ycos) * cost - sy2 * (xcos + ysin) * sint - ) - hmat[j][k] *= model / (sx**5 * sy2) - k += 1 - - if sy_var: - # H(yo,sy)/G = ((x - xo)*sin(t) + (-y + yo)*cos(t))*(2.0*sx**2*sy**2*cos(t) - 1.0*((x - xo)*sin(t) + (-y + yo)*cos(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*cos(t) - sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t)))/(sx**2*sy**5) - hmat[j][k] = xsin - ycos - hmat[j][k] *= 2 * sx2 * sy2 * cost - (xsin - ycos) * ( - sx2 * (xsin - ycos) * cost - sy2 * (xcos + ysin) * sint - ) - hmat[j][k] *= model / (sx2 * sy**5) - k += 1 - - if theta_var: - # H(yo,t)/G = 1.0*(sx**2*sy**2*(sx**2*(-x*cos(2*t) + xo*cos(2*t) - y*sin(2*t) + yo*sin(2*t)) + sy**2*(x*cos(2*t) - xo*cos(2*t) + y*sin(2*t) - yo*sin(2*t))) + (1.0*sx**2 - sy**2)*((x - xo)*sin(t) + (-y + yo)*cos(t))*((x - xo)*cos(t) + (y - yo)*sin(t))*(sx**2*((x - xo)*sin(t) + (-y + yo)*cos(t))*cos(t) - sy**2*((x - xo)*cos(t) + (y - yo)*sin(t))*sin(t)))/(sx**4*sy**4) - hmat[j][k] = (sx2 - sy2) * (xsin - ycos) * (xcos + ysin) - hmat[j][k] *= sx2 * (xsin - ycos) * cost - sy2 * (xcos + ysin) * sint - hmat[j][k] += ( - sx2 - * sy2 - * (sx2 - sy2) - * (-x * cos2t + xo * cos2t - y * sin2t + yo * sin2t) - ) - hmat[j][k] *= model / (sx**4 * sy**4) - # k += 1 - j += 1 - - if sx_var: - k = npvar - if amp_var: - # H(sx,amp)/G = H(amp,sx)/G - hmat[j][k] = hmat[k][j] - k += 1 - - if xo_var: - # H(sx,xo)/G = H(xo,sx)/G - hmat[j][k] = hmat[k][j] - k += 1 - - if yo_var: - # H(sx,yo)/G = H(yo/sx)/G - hmat[j][k] = hmat[k][j] - k += 1 - - # if sx_var: - # H(sx,sx)/G = (-3.0*sx**2 + 1.0*((x - xo)*cos(t) + (y - yo)*sin(t))**2)*((x - xo)*cos(t) + (y - yo)*sin(t))**2/sx**6 - hmat[j][k] = -3 * sx2 + (xcos + ysin) ** 2 - hmat[j][k] *= (xcos + ysin) ** 2 - hmat[j][k] *= model / sx**6 - k += 1 - - if sy_var: - # H(sx,sy)/G = 1.0*((x - xo)*sin(t) + (-y + yo)*cos(t))**2*((x - xo)*cos(t) + (y - yo)*sin(t))**2/(sx**3*sy**3) - hmat[j][k] = (xsin - ycos) ** 2 * (xcos + ysin) ** 2 - hmat[j][k] *= model / (sx**3 * sy**3) - k += 1 - - if theta_var: - # H(sx,t)/G = (-2.0*sx**2*sy**2 + 1.0*(-sx**2 + sy**2)*((x - xo)*cos(t) + (y - yo)*sin(t))**2)*((x - xo)*sin(t) + (-y + yo)*cos(t))*((x - xo)*cos(t) + (y - yo)*sin(t))/(sx**5*sy**2) - hmat[j][k] = -2 * sx2 * sy2 + (sy2 - sx2) * (xcos + ysin) ** 2 - hmat[j][k] *= (xsin - ycos) * (xcos + ysin) - hmat[j][k] *= model / (sx**5 * sy**2) - # k += 1 - j += 1 - - if sy_var: - k = npvar - if amp_var: - # H(sy,amp)/G = H(amp,sy)/G - hmat[j][k] = hmat[k][j] - k += 1 - if xo_var: - # H(sy,xo)/G = H(xo,sy)/G - hmat[j][k] = hmat[k][j] - k += 1 - if yo_var: - # H(sy,yo)/G = H(yo/sy)/G - hmat[j][k] = hmat[k][j] - k += 1 - if sx_var: - # H(sy,sx)/G = H(sx,sy)/G - hmat[j][k] = hmat[k][j] - k += 1 - - # if sy_var: - # H(sy,sy)/G = (-3.0*sy**2 + 1.0*((x - xo)*sin(t) + (-y + yo)*cos(t))**2)*((x - xo)*sin(t) + (-y + yo)*cos(t))**2/sy**6 - hmat[j][k] = -3 * sy2 + (xsin - ycos) ** 2 - hmat[j][k] *= (xsin - ycos) ** 2 - hmat[j][k] *= model / sy**6 - k += 1 - - if theta_var: - # H(sy,t)/G = (2.0*sx**2*sy**2 + 1.0*(-sx**2 + sy**2)*((x - xo)*sin(t) + (-y + yo)*cos(t))**2)*((x - xo)*sin(t) + (-y + yo)*cos(t))*((x - xo)*cos(t) + (y - yo)*sin(t))/(sx**2*sy**5) - hmat[j][k] = 2 * sx2 * sy2 + (sy2 - sx2) * (xsin - ycos) ** 2 - hmat[j][k] *= (xsin - ycos) * (xcos + ysin) - hmat[j][k] *= model / (sx**2 * sy**5) - # k += 1 - j += 1 - - if theta_var: - k = npvar - if amp_var: - # H(t,amp)/G = H(amp,t)/G - hmat[j][k] = hmat[k][j] - k += 1 - if xo_var: - # H(t,xo)/G = H(xo,t)/G - hmat[j][k] = hmat[k][j] - k += 1 - if yo_var: - # H(t,yo)/G = H(yo/t)/G - hmat[j][k] = hmat[k][j] - k += 1 - if sx_var: - # H(t,sx)/G = H(sx,t)/G - hmat[j][k] = hmat[k][j] - k += 1 - if sy_var: - # H(t,sy)/G = H(sy,t)/G - hmat[j][k] = hmat[k][j] - k += 1 - # if theta_var: - # H(t,t)/G = (sx**2*sy**2*(sx**2*(((x - xo)*sin(t) + (-y + yo)*cos(t))**2 - 1.0*((x - xo)*cos(t) + (y - yo)*sin(t))**2) + sy**2*(-1.0*((x - xo)*sin(t) + (-y + yo)*cos(t))**2 + ((x - xo)*cos(t) + (y - yo)*sin(t))**2)) + (sx**2 - 1.0*sy**2)**2*((x - xo)*sin(t) + (-y + yo)*cos(t))**2*((x - xo)*cos(t) + (y - yo)*sin(t))**2)/(sx**4*sy**4) - hmat[j][k] = sx2 * sy2 - hmat[j][k] *= sx2 * ((xsin - ycos) ** 2 - (xcos + ysin) ** 2) + sy2 * ( - (xcos + ysin) ** 2 - (xsin - ycos) ** 2 - ) - hmat[j][k] += (sx2 - sy2) ** 2 * (xsin - ycos) ** 2 * (xcos + ysin) ** 2 - hmat[j][k] *= model / (sx**4 * sy**4) - # j += 1 - - # save the number of variables for the next iteration - # as we need to start our indexing at this number - npvar = k - return np.array(hmat) - - -def emp_hessian(pars, x, y): - """ - Calculate the hessian matrix empirically. - Create a hessian matrix corresponding to the source model 'pars' - Only parameters that vary will contribute to the hessian. - Thus there will be a total of nvar x nvar entries, each of which is a - len(x) x len(y) array. - - Parameters - ---------- - pars : lmfit.Parameters - The model - x, y : list - locations at which to evaluate the Hessian - - Returns - ------- - h : np.array - Hessian. Shape will be (nvar, nvar, len(x), len(y)) - - Notes - ----- - Uses :func:`AegeanTools.fitting.emp_jacobian` to calculate the first order - derivatives. - - See Also - -------- - :func:`AegeanTools.fitting.hessian` - """ - eps = 1e-5 - matrix = [] - for i in range(pars["components"].value): - model = emp_jacobian(pars, x, y) - prefix = "c{0}_".format(i) - for p in ["amp", "xo", "yo", "sx", "sy", "theta"]: - if pars[prefix + p].vary: - pars[prefix + p].value += eps - dm2didj = emp_jacobian(pars, x, y) - model - matrix.append(dm2didj / eps) - pars[prefix + p].value -= eps - matrix = np.array(matrix) - return matrix - - -def nan_acf(noise): - """ - Calculate the autocorrelation function of the noise - where the noise is a 2d array that may contain nans - - - Parameters - ---------- - noise : 2d-array - Noise image. - - Returns - ------- - acf : 2d-array - The ACF. - """ - corr = np.zeros(noise.shape) - ix, jx = noise.shape - for i in range(ix): - si_min = slice(i, None, None) - si_max = slice(None, ix - i, None) - for j in range(jx): - sj_min = slice(j, None, None) - sj_max = slice(None, jx - j, None) - if np.all(np.isnan(noise[si_min, sj_min])) or np.all( - np.isnan(noise[si_max, sj_max]) - ): - corr[i, j] = np.nan - else: - corr[i, j] = np.nansum(noise[si_min, sj_min] * noise[si_max, sj_max]) - # return the normalised acf - return corr / np.nanmax(corr) - - -def make_ita(noise, acf=None): - """ - Create the matrix ita of the noise where the noise may be a masked array - where ita(x,y) is the correlation between pixel pairs that have the same - separation as x and y. - - Parameters - ---------- - noise : 2d-array - The noise image - - acf : 2d-array - The autocorrelation matrix. (None = calculate from data). - Default = None. - - Returns - ------- - ita : 2d-array - The matrix ita - """ - if acf is None: - acf = nan_acf(noise) - # s should be the number of non-masked pixels - s = np.count_nonzero(np.isfinite(noise)) - # the indices of the non-masked pixels - xm, ym = np.where(np.isfinite(noise)) - ita = np.zeros((s, s)) - # iterate over the pixels - for i, (x1, y1) in enumerate(zip(xm, ym)): - for j, (x2, y2) in enumerate(zip(xm, ym)): - k = abs(x1 - x2) - l = abs(y1 - y2) - ita[i, j] = acf[k, l] - return ita - - -def RB_bias(data, pars, ita=None, acf=None): - """ - Calculate the expected bias on each of the parameters in the model pars. - Only parameters that are allowed to vary will have a bias. - Calculation follows the description of Refrieger & Brown 1998 (cite). - - - Parameters - ---------- - data : 2d-array - data that was fit - - pars : lmfit.Parameters - The model - - ita : 2d-array - The ita matrix (optional). - - acf : 2d-array - The acf for the data. - - Returns - ------- - bias : array - The bias on each of the parameters - """ - logger.info("data {0}".format(data.shape)) - nparams = np.sum([pars[k].vary for k in pars.keys() if k != "components"]) - # masked pixels - xm, ym = np.where(np.isfinite(data)) - # all pixels - x, y = np.indices(data.shape) - # Create the jacobian as an AxN array accounting for the masked pixels - j = np.array(np.vsplit(lmfit_jacobian(pars, xm, ym).T, nparams)).reshape( - nparams, -1 - ) - - h = hessian(pars, x, y) - # mask the hessian to be AxAxN array - h = h[:, :, xm, ym] - Hij = np.einsum("ik,jk", j, j) - Dij = np.linalg.inv(Hij) - Bijk = np.einsum("ip,jkp", j, h) - Eilkm = np.einsum("il,km", Dij, Dij) - - Cimn_1 = -1 * np.einsum("krj,ir,km,jn", Bijk, Dij, Dij, Dij) - Cimn_2 = -1.0 / 2 * np.einsum("rkj,ir,km,jn", Bijk, Dij, Dij, Dij) - Cimn = Cimn_1 + Cimn_2 - - if ita is None: - # N is the noise (data-model) - N = data - ntwodgaussian_lmfit(pars)(x, y) - if acf is None: - acf = nan_acf(N) - ita = make_ita(N, acf=acf) - logger.info("acf.shape {0}".format(acf.shape)) - logger.info("acf[0] {0}".format(acf[0])) - logger.info("ita.shape {0}".format(ita.shape)) - logger.info("ita[0] {0}".format(ita[0])) - - # Included for completeness but not required - - # now mask/ravel the noise - # N = N[np.isfinite(N)].ravel() - # Pi = np.einsum('ip,p', j, N) - # Qij = np.einsum('ijp,p', h, N) - - Vij = np.einsum("ip,jq,pq", j, j, ita) - Uijk = np.einsum("ip,jkq,pq", j, h, ita) - - bias_1 = np.einsum("imn, mn", Cimn, Vij) - bias_2 = np.einsum("ilkm, mlk", Eilkm, Uijk) - bias = bias_1 + bias_2 - logger.info("bias {0}".format(bias)) - return bias - - -def bias_correct(params, data, acf=None): - """ - Calculate and apply a bias correction to the given fit parameters - - - Parameters - ---------- - params : lmfit.Parameters - The model parameters. These will be modified. - - data : 2d-array - The data which was used in the fitting - - acf : 2d-array - ACF of the data. Default = None. - - Returns - ------- - None - - See Also - -------- - :func:`AegeanTools.fitting.RB_bias` - """ - bias = RB_bias(data, params, acf=acf) - i = 0 - for p in params: - if "theta" in p: - continue - if params[p].vary: - params[p].value -= bias[i] - i += 1 - return - - def errors(source, model, wcshelper): """ Convert pixel based errors into sky coord errors @@ -1227,7 +689,65 @@ def rfunc(x, y): return rfunc -def do_lmfit(data, params, B=None, errs=None, dojac=True): +def nthreedgaussian_lmfit(params): + """ + Convert an lmfit.Parameters object into a function which calculates the + model. + + + Parameters + ---------- + params : lmfit.Parameters + Model parameters, can have multiple components. + + Returns + ------- + model : func + A function f(x,y) that will compute the model. + """ + + def rfunc(v, x, y): # TODO: Update the Doc string + #! v is not a pixel coordinate it should be actual frequency i.e. pix2freq + """ + Compute the model given by params, at pixel coordinates x,y + + Parameters + ---------- + x, y : numpy.ndarray + The x/y pixel coordinates at which the model is being evaluated + + Returns + ------- + result : numpy.ndarray + Model + """ + result = None + + for i in range(int(params["components"].value)): + prefix = f"c{i}_" + # I hope this doesn't kill our run time + amp = np.nan_to_num(params[prefix + "amp"].value) + xo = params[prefix + "xo"].value + yo = params[prefix + "yo"].value + sx = params[prefix + "sx"].value + sy = params[prefix + "sy"].value + theta = params[prefix + "theta"].value + alpha = params[prefix + "alpha"].value + nu0 = params[prefix + "nu0"].value + if result is not None: + result += elliptical_gaussian_with_alpha( + x, y, v, amp, xo, yo, nu0, sx, sy, theta, alpha + ) # TODO: Pass the frequency into this, which is the current frequency + else: + result = elliptical_gaussian_with_alpha( + x, y, v, amp, xo, yo, nu0, sx, sy, theta, alpha + ) + return result + + return rfunc + + +def do_lmfit(data, params, B=None, errs=None, dojac=False): """ Fit the model to the data data may contain 'flagged' or 'masked' data with the value of np.nan @@ -1313,6 +833,100 @@ def residual(params, x, y, B=None, errs=None): return result, params +def do_lmfit_3D( + data, + params, + freq_mapping=None, + B=None, + errs=None, + dojac=False, +): # TODO: Go through B matrix + """ + Fit the model to the data + data may contain 'flagged' or 'masked' data with the value of np.NaN + + Parameters + ---------- + data : 2d-array + Image data + + params : lmfit.Parameters + Initial model guess. + + B : 2d-array + B matrix to be used in residual calculations. + Default = None. + + errs : 1d-array + + dojac : bool + If true then an analytic jacobian will be passed to the fitter + + Returns + ------- + result : ? + lmfit.minimize result. + + params : lmfit.Params + Fitted model. + + See Also + -------- + :func:`AegeanTools.fitting.lmfit_jacobian` + + """ + # copy the params so as not to change the initial conditions + # in case we want to use them elsewhere + params = copy.deepcopy(params) + data = np.array(data) + mask = np.where(np.isfinite(data)) + fmask = freq_mapping[mask[0]] + + def residual(params, x, y, B=None, errs=None): + """ + The residual function required by lmfit + + Parameters + ---------- + params: lmfit.Params + The parameters of the model being fit + + Returns + ------- + result : numpy.ndarray + Model - Data + """ + f = nthreedgaussian_lmfit(params) # A function describing the model + model = f(fmask, mask[1], mask[2]) # The actual model + + if np.any(~np.isfinite(model)): + logger.debug(f"The parameters are {params}") + raise AegeanNaNModelError( + "lmfit optimisation has return NaN in the parameter set. " + ) + + if B is None: + return model - data[mask] + else: + return (model - data[mask]).dot(B) + + # if dojac: + # result = lmfit.minimize( + # residual, + # params, + # kws={"x": mask[0], "y": mask[1], "B": B, "errs": errs}, + # Dfun=lmfit_jacobian, + # ) + result = lmfit.minimize( + residual, params, kws={"x": mask[0], "y": mask[1], "B": B, "errs": errs} + ) + + # Remake the residual so that it is once again (model - data) + # if B is not None: + # result.residual = result.residual.dot(inv(B)) + return result, params + + def covar_errors(params, data, errs, B, C=None): r""" Take a set of parameters that were fit with lmfit, and replace the errors diff --git a/AegeanTools/models.py b/AegeanTools/models.py index 065aba45..2b2fc70f 100644 --- a/AegeanTools/models.py +++ b/AegeanTools/models.py @@ -49,16 +49,32 @@ class SimpleSource(object): -------- :mod:`AegeanTools.flags` """ - header = "#RA DEC Flux err a b pa flags\n" + \ - "# Jy/beam Jy/beam '' '' deg ZWNCPES\n" + \ - "#========================================================================" - - formatter = "{0.ra:11.7f} {0.dec:11.7f} {0.peak_flux: 8.6f} " + \ - "{0.err_peak_flux: 8.6f} {0.a:5.2f} {0.b:5.2f} " + \ - "{0.pa:6.1f} {0.flags:07b}" - names = ['background', 'local_rms', 'ra', 'dec', 'peak_flux', - 'err_peak_flux', 'flags', 'peak_pixel', 'a', 'b', - 'pa', 'uuid'] + + header = ( + "#RA DEC Flux err a b pa flags\n" + + "# Jy/beam Jy/beam '' '' deg ZWNCPES\n" + + "#========================================================================" + ) + + formatter = ( + "{0.ra:11.7f} {0.dec:11.7f} {0.peak_flux: 8.6f} " + + "{0.err_peak_flux: 8.6f} {0.a:5.2f} {0.b:5.2f} " + + "{0.pa:6.1f} {0.flags:07b}" + ) + names = [ + "background", + "local_rms", + "ra", + "dec", + "peak_flux", + "err_peak_flux", + "flags", + "peak_pixel", + "a", + "b", + "pa", + "uuid", + ] galactic = False def __init__(self): @@ -184,18 +200,38 @@ class IslandSource(SimpleSource): :mod:`AegeanTools.flags` """ - names = ['island', 'components', 'background', 'local_rms', 'ra_str', - 'dec_str', 'ra', 'dec', 'peak_flux', 'int_flux', 'err_int_flux', - 'eta', 'x_width', 'y_width', 'max_angular_size', 'pa', - 'pixels', 'area', 'beam_area', 'flags', 'uuid'] + + names = [ + "island", + "components", + "background", + "local_rms", + "ra_str", + "dec_str", + "ra", + "dec", + "peak_flux", + "int_flux", + "err_int_flux", + "eta", + "x_width", + "y_width", + "max_angular_size", + "pa", + "pixels", + "area", + "beam_area", + "flags", + "uuid", + ] def __init__(self): SimpleSource.__init__(self) self.island = 0 # island number # background = None # local background zero point # local_rms= None #local image rms - self.ra_str = '' # str - self.dec_str = '' # str + self.ra_str = "" # str + self.dec_str = "" # str # ra = None # degrees # dec = None # degrees # peak_flux = None # Jy/beam @@ -221,19 +257,19 @@ def __str__(self): return "({0:d})".format(self.island) def __eq__(self, other): - if hasattr(other, 'island'): + if hasattr(other, "island"): return self.island == other.island else: return False def __ne__(self, other): - if hasattr(other, 'island'): + if hasattr(other, "island"): return self.island != other.island else: return True def __lt__(self, other): - if hasattr(other, 'island'): + if hasattr(other, "island"): return self.island < other.island else: return True @@ -242,7 +278,7 @@ def __le__(self, other): return self.__lt__(other) or self.__eq__(other) def __gt__(self, other): - if hasattr(other, 'island'): + if hasattr(other, "island"): return self.island > other.island else: return False @@ -251,7 +287,7 @@ def __ge__(self, other): return self.__gt__(other) or self.__eq__(other) -class ComponentSource(SimpleSource): +class ComponentSource(SimpleSource): #! <- This is the entry point """ A Gaussian component, aka a source, that was measured by Aegean. @@ -308,26 +344,57 @@ class ComponentSource(SimpleSource): :mod:`AegeanTools.flags` """ + # header for the output - header = "#isl,src bkg rms RA DEC RA err DEC err Peak err S_int err a err b err pa err flags\n" + \ - "# Jy/beam Jy/beam deg deg deg deg Jy/beam Jy/beam Jy Jy '' '' '' '' deg deg ZWNCPES\n" + \ - "#============================================================================================================================================================================================" + header = ( + "#isl,src bkg rms RA DEC RA err DEC err Peak err S_int err a err b err pa err flags\n" + + "# Jy/beam Jy/beam deg deg deg deg Jy/beam Jy/beam Jy Jy '' '' '' '' deg deg ZWNCPES\n" + + "#============================================================================================================================================================================================" + ) # formatting strings for making nice output - formatter = "({0.island:04d},{0.source:02d}) {0.background: 8.6f} " + \ - "{0.local_rms: 8.6f} {0.ra_str:12s} {0.dec_str:12s} " + \ - "{0.ra:11.7f} {0.err_ra: 9.7f} {0.dec:11.7f} " + \ - "{0.err_dec: 9.7f} {0.peak_flux: 8.6f} " + \ - "{0.err_peak_flux: 8.6f} {0.int_flux: 8.6f} " + \ - "{0.err_int_flux: 8.6f} {0.a:5.2f} {0.err_a:5.2f} " + \ - "{0.b:5.2f} {0.err_b:5.2f} {0.pa:6.1f} {0.err_pa:5.1f} " + \ - "{0.flags:07b}" - names = ['island', 'source', 'background', 'local_rms', - 'ra_str', 'dec_str', 'ra', 'err_ra', 'dec', 'err_dec', - 'peak_flux', 'err_peak_flux', 'int_flux', 'err_int_flux', - 'a', 'err_a', 'b', 'err_b', 'pa', 'err_pa', - 'flags', 'residual_mean', 'residual_std', - 'uuid', 'psf_a', 'psf_b', 'psf_pa'] + formatter = ( + "({0.island:04d},{0.source:02d}) {0.background: 8.6f} " + + "{0.local_rms: 8.6f} {0.ra_str:12s} {0.dec_str:12s} " + + "{0.ra:11.7f} {0.err_ra: 9.7f} {0.dec:11.7f} " + + "{0.err_dec: 9.7f} {0.peak_flux: 8.6f} " + + "{0.err_peak_flux: 8.6f} {0.int_flux: 8.6f} " + + "{0.err_int_flux: 8.6f} {0.a:5.2f} {0.err_a:5.2f} " + + "{0.b:5.2f} {0.err_b:5.2f} {0.pa:6.1f} {0.err_pa:5.1f} " + + "{0.flags:07b}" + ) # TODO: Update these in the new class, try to inherit it and add to it, alpha is either positive or negative single digit + # TODO : 2 decimal places e.g. 1.04 +- would be good if they have + infront of them + # TODO: formatter is +4.2f for the alpha + # TODO: For nu0 unsure if megahertz or gighertz, for now assume megahertz between 100Mhz and a few 10s of GHz, within 0.5 Mhz resolution + names = [ + "island", + "source", + "background", + "local_rms", + "ra_str", + "dec_str", + "ra", + "err_ra", + "dec", + "err_dec", + "peak_flux", + "err_peak_flux", + "int_flux", + "err_int_flux", + "a", + "err_a", + "b", + "err_b", + "pa", + "err_pa", + "flags", + "residual_mean", + "residual_std", + "uuid", + "psf_a", + "psf_b", + "psf_pa", + ] # TODO: Ditto for this, add alpha and nu0 to the names def __init__(self): SimpleSource.__init__(self) @@ -335,8 +402,8 @@ def __init__(self): self.source = 0 # source number # background = None # local background zero point # local_rms= None #local image rms - self.ra_str = '' # str - self.dec_str = '' # str + self.ra_str = "" # str + self.dec_str = "" # str # ra = None # degrees self.err_ra = np.nan # degrees # dec = None # degrees @@ -369,22 +436,22 @@ def __repr__(self): def __eq__(self, other): if self.island != other.island: return False - if not hasattr(other, 'source'): + if not hasattr(other, "source"): return False return self.source == other.source def __ne__(self, other): if self.island != other.island: return True - if not hasattr(other, 'source'): + if not hasattr(other, "source"): return True return self.source != other.source def __lt__(self, other): - if not hasattr(other, 'island'): + if not hasattr(other, "island"): return True # Islands are always less than components - if not hasattr(other, 'source'): + if not hasattr(other, "source"): return True if self.island < other.island: return True @@ -392,10 +459,10 @@ def __lt__(self, other): return self.source < other.source def __le__(self, other): - if not hasattr(other, 'island'): + if not hasattr(other, "island"): return True # Islands are always less than components - if not hasattr(other, 'source'): + if not hasattr(other, "source"): return True if self.island < other.island: return True @@ -403,10 +470,10 @@ def __le__(self, other): return self.source <= other.source def __gt__(self, other): - if not hasattr(other, 'island'): + if not hasattr(other, "island"): return False # Islands are always less than components - if not hasattr(other, 'source'): + if not hasattr(other, "source"): return False if self.island > other.island: return True @@ -414,10 +481,10 @@ def __gt__(self, other): return self.source > other.source def __ge__(self, other): - if not hasattr(other, 'island'): + if not hasattr(other, "island"): return False # Islands are always less than components - if not hasattr(other, 'source'): + if not hasattr(other, "source"): return False if self.island > other.island: return True @@ -425,6 +492,117 @@ def __ge__(self, other): return self.source >= other.source +class ComponentSource3D(ComponentSource): + """ + A 3-D Gaussian component, aka a source, that was measured by Aegean, + with additional attributes for spectral index and frequency. + + Inherits from the ComponentSource Class. + + Attributes + ---------- + alpha : float + Spectral index of the source, showing how flux varies with frequency. + Typically a small positive or negative value. + + nu0 : float + The reference frequency (in MHz) at which the source was measured. + + Other attributes are inherited from ParentClass: + + island : int + The island which this component is part of. + + source : int + The source number within the island. + + background, local_rms : float + Background and local noise level in the image at the location + of this source. + + ra, err_ra, dec, err_dec : float + Sky location of the source including uncertainties. Decimal degrees. + + ra_str, dec_str : str + Sky location in HH:MM:SS.SS +DD:MM:SS.SS format. + + galactic : bool + If true then ra,dec are interpreted as glat,glon instead. + Default = False. + This is a class attribute, not an instance attribute. + + peak_flux, err_peak_flux : float + The peak flux and associated uncertainty. + + int_flux, err_int_flux : float + Integrated flux and associated uncertainty. + + a, err_a, b, err_b, pa, err_pa: float + Shape parameters for this source and associated uncertainties. + a/b are in arcsec, pa is in degrees East of North. + + residual_mean, residual_std : float + The mean and standard deviation of the model-data for this island + of pixels. + + psf_a, psf_b, psf_pa : float + The shape parameters for the point spread function + (degrees). + + flags : int + Flags. See :mod:`AegeanTools.flags`. + + uuid : str + Unique ID for this source. This is random and not dependent on the + source properties. + + See Also + -------- + :mod:`AegeanTools.flags` + + """ + + header = ( + "#isl,src bkg rms RA DEC RA err DEC err Peak err S_int err a err b err pa err flags alpha nu0\n" + + "# Jy/beam Jy/beam deg deg deg deg Jy/beam Jy/beam Jy Jy '' '' '' '' deg deg ZWNCPES '' Mhz\n" + + "#=============================================================================================================================================================================================================" + ) + formatter = ( + ComponentSource.formatter + " {0.alpha: 4.2f} {0.nu0: 4.2f}" + ) # Add the alpha and nu0 attributes to the formatter + names = ComponentSource.names + [ + "alpha", + "nu0", + ] # Add the alpha and nu0 attributes to the names + + def __init__(self): + super().__init__() # Call the parent class constructor + self.alpha = 0 # Add the alpha attribute + self.nu0 = 0 # Add the nu0 attribute + + @staticmethod + def from_component_source(component_source): + """ + Create a new ComponentSource3D object from a ComponentSource object. + + Parameters + ---------- + component_source : ComponentSource + The ComponentSource object to convert to a ComponentSource3D object. + + Returns + ------- + ComponentSource3D + The new ComponentSource3D object. + """ + component_source_3d = ComponentSource3D() + for name in component_source.names: + # Set the attribute of the ComponentSource3D object to the + # attribute of the ComponentSource object + setattr(component_source_3d, name, getattr(component_source, name)) + return component_source_3d + + class PixelIsland(object): """ An island of pixels within an image or cube @@ -458,9 +636,10 @@ def set_mask(self, data): """ if len(data.shape) != self.dim: raise AssertionError( - ("mask shape {0} is of the wrong dimension. " + - "Expecting {1}").format(data.shape, self.dim) + ("mask shape {0} is of the wrong dimension. " + "Expecting {1}").format( + data.shape, self.dim ) + ) self.mask = data return @@ -482,8 +661,9 @@ def calc_bounding_box(self, data, offsets): if len(offsets) != self.dim: raise AssertionError( "{0} offsets were passed but {1} are required".format( - len(offsets), self.dim) - ) + len(offsets), self.dim + ) + ) # TODO: Figure out 3d boxes # set the bounding box one dimension at a time ndrow = np.any(data, axis=0) @@ -494,7 +674,7 @@ def calc_bounding_box(self, data, offsets): cmin, cmax = np.where(ndcol)[0][[0, -1]] self.bounding_box[0][0] = offsets[0] + cmin self.bounding_box[0][1] = offsets[0] + cmax + 1 - self.set_mask(data[rmin:rmax+1, cmin:cmax+1]) + self.set_mask(data[rmin : rmax + 1, cmin : cmax + 1]) return @@ -523,9 +703,9 @@ class IslandFittingData(object): If true then also measure properties of the island. """ - def __init__(self, isle_num=0, i=None, - scalars=None, offsets=(0, 0, 1, 1), - doislandflux=False): + def __init__( + self, isle_num=0, i=None, scalars=None, offsets=(0, 0, 1, 1), doislandflux=False + ): self.isle_num = isle_num self.i = i self.scalars = scalars @@ -577,7 +757,9 @@ def classify_catalog(catalog): islands = [] simples = [] for source in catalog: - if isinstance(source, ComponentSource): + if isinstance(source, ComponentSource) or isinstance( + source, ComponentSource3D + ): # TODO This needs to be split up down the line components.append(source) elif isinstance(source, IslandSource): islands.append(source) diff --git a/AegeanTools/mpi.py b/AegeanTools/mpi.py new file mode 100644 index 00000000..bdf0134e --- /dev/null +++ b/AegeanTools/mpi.py @@ -0,0 +1,12 @@ +#!/usr/bin/env python + +MPI = None +MPI_AVAIL = False + +try: + import mpi4py + from mpi4py import MPI + if MPI.COMM_WORLD.Get_size() > 1: + MPI_AVAIL = True +except ImportError: + pass \ No newline at end of file diff --git a/AegeanTools/source_finder.py b/AegeanTools/source_finder.py index c1277c8d..6501f047 100755 --- a/AegeanTools/source_finder.py +++ b/AegeanTools/source_finder.py @@ -24,16 +24,18 @@ from .fitting import ( Bmatrix, Cmatrix, - bias_correct, covar_errors, do_lmfit, + do_lmfit_3D, elliptical_gaussian, errors, ntwodgaussian_lmfit, + nthreedgaussian_lmfit, ) from .logging import logger, logging from .models import ( ComponentSource, + ComponentSource3D, DummyLM, IslandFittingData, IslandSource, @@ -45,6 +47,8 @@ from .regions import Region from .wcs_helpers import Beam, WCSHelper +from .mpi import MPI_AVAIL, MPI + __author__ = "Paul Hancock" header = """#Aegean version {0} @@ -109,7 +113,7 @@ def find_islands( l, n = label(a, structure=np.ones((3, 3))) f = find_objects(l) - logger.debug("{1} Found {0} islands total above flood limit".format(n, im.shape)) + logger.debug(f"{n} Found {im.shape} islands total above flood limit") islands = [] for i in range(n): @@ -187,9 +191,6 @@ class SourceFinder(object): docov : bool If True, then fitting will be done using the covariance matrix. - dobais: bool - Not yet implemented. Default = False - sources : [ComponentSource|IslandSource|SimpleSource] List of sources that have been found/measured. """ @@ -207,7 +208,6 @@ def __init__(self, **kwargs): self.wcshelper = None self.blank = False self.docov = True - self.dobias = False self.cube_index = 0 self.sources = [] @@ -642,23 +642,33 @@ def result_to_components(self, result, model, island_data, isflags): sources : list A list of components, and islands if requested. """ - + threeD = len(self.img.shape) == 3 # island data isle_num = island_data.isle_num idata = island_data.i xmin, xmax, ymin, ymax = island_data.offsets - box = slice(int(xmin), int(xmax)), slice(int(ymin), int(ymax)) - rms = self.rmsimg[box] - bkg = self.bkgimg[box] + if threeD: + box = slice(0, -1), slice(int(xmin), int(xmax)), slice(int(ymin), int(ymax)) + logger.debug(f"THIS IS THE 3-D BOX: {box}") + + rms = np.mean(self.rmsimg[box], axis=0) + bkg = np.mean(self.bkgimg[box], axis=0) + + else: + box = slice(int(xmin), int(xmax)), slice(int(ymin), int(ymax)) + logger.debug(f"THIS IS THE 2-D BOX: {box}") + + rms = self.rmsimg[box] + bkg = self.bkgimg[box] + residual = np.median(result.residual), np.std(result.residual) is_flag = isflags - sources = [] j = 0 for j in range(int(model["components"].value)): src_flags = is_flag - source = ComponentSource() + source = ComponentSource3D() if threeD else ComponentSource() source.island = isle_num source.source = j logger.debug(" component {0}".format(j)) @@ -671,6 +681,10 @@ def result_to_components(self, result, model, island_data, isflags): amp = model[prefix + "amp"].value src_flags |= int(model[prefix + "flags"].value) + if threeD: + source.alpha = model[prefix + "alpha"].value + source.nu0 = model[prefix + "nu0"].value + # these are goodness of fit statistics for the entire island. source.residual_mean = residual[0] source.residual_std = residual[1] @@ -887,6 +901,7 @@ def load_globals( blank=False, docov=True, cube_index=None, + as_cube=False, ): """ Populate the global_data object by loading or calculating @@ -937,74 +952,87 @@ def load_globals( cube_index : int For an image cube, which slice to use. + as_cube : bool + If the data is to be processed as a cube or not. """ # don't reload already loaded data if self.img is not None: return img, header = load_image_band( - filename, hdu_index=hdu_index, cube_index=cube_index + filename, hdu_index=hdu_index, cube_index=cube_index, as_cube=as_cube ) debug = logger.isEnabledFor(logging.DEBUG) - if mask is not None: - # allow users to supply and object instead of a filename - if isinstance(mask, Region): - self.region = mask - elif os.path.exists(mask): - logger.info("Loading mask from {0}".format(mask)) - self.region = Region.load(mask) - else: - logger.error("File {0} not found for loading".format(mask)) - self.region = None + if not as_cube: + self.img = img + self.header = header + self.dtype = type(self.img[0][0]) + self.bkgimg = np.zeros(self.img.shape, dtype=self.dtype) + self.rmsimg = np.zeros(self.img.shape, dtype=self.dtype) + + self.cube_index = cube_index + + self.wcshelper = WCSHelper.from_header(header, beam, psf_file=psf) + self.beam = self.wcshelper.beam + + if mask is not None: + # allow users to supply an object instead of a filename + if isinstance(mask, Region): + self.region = mask + elif os.path.exists(mask): + logger.info("Loading mask from {0}".format(mask)) + self.region = Region.load(mask) + else: + logger.error("File {0} not found for loading".format(mask)) + self.region = None + + if do_curve: + logger.info("Calculating curvature") + # calculate curvature but store it as -1,0,+1 + dcurve = np.zeros(self.img.shape, dtype=np.int8) + peaks = maximum_filter(self.img, size=3) + troughs = minimum_filter(self.img, size=3) + pmask = np.where(self.img == peaks) + tmask = np.where(self.img == troughs) + dcurve[pmask] = -1 + dcurve[tmask] = 1 + self.dcurve = dcurve + + # if either of rms or bkg images are not supplied + # then calculate them both + if not (rmsin and bkgin): + if verb: + logger.info("Calculating background and rms data") + self._make_bkg_rms( + filename=filename, + forced_rms=rms, + forced_bkg=bkg, + cores=cores, + ) - self.wcshelper = WCSHelper.from_header(header, beam, psf_file=psf) - self.beam = self.wcshelper.beam + else: + self.img = img + self.header = header + self.dtype = type(self.img[0][0]) - self.img = img - self.header = header - self.dtype = type(self.img[0][0]) - self.bkgimg = np.zeros(self.img.shape, dtype=self.dtype) - self.rmsimg = np.zeros(self.img.shape, dtype=self.dtype) - - self.cube_index = cube_index - - if do_curve: - logger.info("Calculating curvature") - # calculate curvature but store it as -1,0,+1 - dcurve = np.zeros(self.img.shape, dtype=np.int8) - peaks = maximum_filter(self.img, size=3) - troughs = minimum_filter(self.img, size=3) - pmask = np.where(self.img == peaks) - tmask = np.where(self.img == troughs) - dcurve[pmask] = -1 - dcurve[tmask] = 1 - self.dcurve = dcurve - - # if either of rms or bkg images are not supplied - # then calculate them both - if not (rmsin and bkgin): - if verb: - logger.info("Calculating background and rms data") - self._make_bkg_rms( - filename=filename, - forced_rms=rms, - forced_bkg=bkg, - cores=cores, - ) + self.wcshelper = WCSHelper.from_header(header, beam) + self.beam = self.wcshelper.beam # replace the calculated images with input versions, # if the user has supplied them. if bkgin: if verb: logger.info("Loading background data from file {0}".format(bkgin)) - self.bkgimg = self._load_aux_image(img, bkgin) + self.bkgimg = self._load_aux_image(img, bkgin, as_cube=as_cube) if rmsin: if verb: logger.info("Loading rms data from file {0}".format(rmsin)) - self.rmsimg = self._load_aux_image(img, rmsin) + self.rmsimg = self._load_aux_image(img, rmsin, as_cube=as_cube) + self.bkgimg = np.squeeze(self.bkgimg) + self.rmsimg = np.squeeze(self.rmsimg) # subtract the background image from the data image and save if verb and debug: logger.debug("Data max is {0}".format(np.nanmax(img))) @@ -1017,9 +1045,6 @@ def load_globals( self.blank = blank self.docov = docov - # Default to false until I can verify that this is working - self.dobias = False - # check if the WCS is galactic if "lon" in self.header["CTYPE1"].lower(): logger.info("Galactic coordinates detected and noted") @@ -1219,7 +1244,7 @@ def _make_bkg_rms(self, filename, forced_rms=None, forced_bkg=None, cores=None): return - def _load_aux_image(self, image, auxfile): + def _load_aux_image(self, image, auxfile, as_cube=False): """ Load a fits file (bkg/rms/curve) and make sure that it is the same shape as the main image. @@ -1238,7 +1263,7 @@ def _load_aux_image(self, image, auxfile): The loaded image. """ - auximg, _ = load_image_band(auxfile) + auximg, _ = load_image_band(auxfile, as_cube=as_cube) if auximg.shape != image.shape: logger.error( @@ -1286,9 +1311,7 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): for inum, isle in enumerate(group, start=istart): logger.debug("-=-") - logger.debug( - "input island = {0}, {1} components".format(isle[0].island, len(isle)) - ) + logger.debug(f"input island = {isle[0].island}, {len(isle)} components") # set up the parameters for each of the sources within the island i = 0 @@ -1303,6 +1326,7 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): # this may be a subset of all sources in the island included_sources = [] for src in isle: + # get the beam for this source pixbeam = Beam(*self.wcshelper.get_psf_sky2pix(src.ra, src.dec)) # find the right pixels from the ra/dec source_x, source_y = self.wcshelper.sky2pix([src.ra, src.dec]) @@ -1311,9 +1335,7 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): x = int(round(source_x)) y = int(round(source_y)) - logger.debug( - "pixel location ({0:5.2f},{1:5.2f})".format(source_x, source_y) - ) + logger.debug(f"pixel location ({source_x:5.2f},{source_y:5.2f})") # reject sources that are outside the image bounds, # or which have nan data/rms values if ( @@ -1324,9 +1346,7 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): or pixbeam is None ): logger.debug( - "Source ({0},{1}) not within usable region: skipping".format( - src.island, src.source - ) + f"Source ({src.island},{src.source}) not within usable region: skipping" ) continue else: @@ -1341,14 +1361,10 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): sy *= FWHM2CC logger.debug( - "Source shape [sky coords] {0:5.2f}x{1:5.2f}@{2:05.2f}".format( - src.a, src.b, src.pa - ) + f"Source shape [sky coords] {src.a:5.2f}x{src.b:5.2f}@{src.pa:05.2f}" ) logger.debug( - "Source shape [pixel coords] {0:4.2f}x{1:4.2f}@{2:05.2f}".format( - sx, sy, theta - ) + f"Source shape [pixel coords] {sx:4.2f}x{sy:4.2f}@{theta:05.2f}" ) # choose a region that is 2x the major axis of the source, @@ -1366,7 +1382,7 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): s_lims = [0.8 * min(sx, pixbeam.b * FWHM2CC), max(sy, sx) * 1.25] # Set up the parameters for the fit, including constraints - prefix = "c{0}_".format(i) + prefix = f"c{i}_" params.add(prefix + "amp", value=src.peak_flux, vary=True) # for now the xo/yo are locations within the main image, # we correct this later @@ -1400,7 +1416,6 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): ) params.add(prefix + "theta", value=theta, vary=stage >= 3) params.add(prefix + "flags", value=0, vary=False) - # this source is being refit so add it to the list included_sources.append(src) i += 1 @@ -1408,18 +1423,20 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): # the FWHM mask that is defined further on if i == 0: - logger.debug("No sources found in island {0}".format(src.island)) + logger.debug(f"No sources found in island {src.island}") continue params.add("components", value=i, vary=False) # params.components = i - logger.debug(" {0} components being fit".format(i)) + logger.debug(f" {i} components being fit") # now we correct the xo/yo positions to be # relative to the sub-image - logger.debug("xmxxymyx {0} {1} {2} {3}".format(xmin, xmax, ymin, ymax)) + logger.debug(f"xmxxymyx {xmin} {xmax} {ymin} {ymax}") for i in range(int(params["components"].value)): - prefix = "c{0}_".format(i) + prefix = f"c{i}_" + # must update limits before the value as limits are # enforced when the value is updated + # ? From my understanding this ensures the sources are within the sub-image params[prefix + "xo"].min -= xmin params[prefix + "xo"].max -= xmin params[prefix + "xo"].value -= xmin @@ -1429,12 +1446,13 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): # logger.debug(params) # don't fit if there are no sources if params["components"].value < 1: - logger.info("Island {0} has no components".format(src.island)) + logger.info(f"Island {src.island} has no components") continue # this .copy() will stop us from modifying the parent region when # we later apply our mask. idata = data[int(xmin) : int(xmax), int(ymin) : int(ymax)].copy() + # idata = data[zmin:zmax ,int(xmin) : int(xmax), int(ymin) : int(ymax)].copy() # now convert these back to indices within the idata region # island_mask = np.array([(x-xmin, y-ymin) for x,y in island_mask]) @@ -1443,7 +1461,7 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): # of the sources being fit mask_params = copy.deepcopy(params) for i in range(int(mask_params["components"].value)): - prefix = "c{0}_".format(i) + prefix = f"c{i}_" mask_params[prefix + "amp"].value = 1 mask_model = ntwodgaussian_lmfit(mask_params) mask = np.where(mask_model(allx.ravel(), ally.ravel()) <= 0.1) @@ -1468,14 +1486,14 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): # the central 3x3 pixels of it's location # If not then we don't fit that component for i in range(int(params["components"].value)): - prefix = "c{0}_".format(i) + prefix = f"c{i}_" # figure out a box around the center of this cx, cy = ( params[prefix + "xo"].value, params[prefix + "yo"].value, ) # central pixel coords - logger.debug(" comp {0}".format(i)) - logger.debug(" x0, y0 {0} {1}".format(cx, cy)) + logger.debug(f" comp {i}") + logger.debug(f" x0, y0 {cx} {cy}") xmx = int(round(np.clip(cx + 2, 0, idata.shape[0]))) xmn = int(round(np.clip(cx - 1, 0, idata.shape[0]))) ymx = int(round(np.clip(cy + 2, 0, idata.shape[1]))) @@ -1484,7 +1502,7 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): # if there are no not-nan pixels in this region # then don't vary any parameters if not np.any(np.isfinite(square)): - logger.debug(" not fitting component {0}".format(i)) + logger.debug(f" not fitting component {i}") params[prefix + "amp"].value = np.nan for p in ["amp", "xo", "yo", "sx", "sy", "theta"]: params[prefix + p].vary = False @@ -1503,9 +1521,7 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): else: if non_nan_pix < nfree: logger.debug( - "More free parameters {0} than available pixels {1}".format( - nfree, non_nan_pix - ) + f"More free parameters {nfree} than available pixels {non_nan_pix}" ) if non_nan_pix >= params["components"].value: logger.debug("Fixing all parameters except amplitudes") @@ -1519,6 +1535,7 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): # do the fit # if the pixel beam is not valid, then recalculate using # the location of the last source to have a valid psf + #! If we don't use the psf, then we don't need to recalculate the pixel beam if pixbeam is None: if src_valid_psf is not None: pixbeam = self.wcshelper.get_pixbeam( @@ -1574,6 +1591,329 @@ def _refit_islands(self, group, stage, outerclip=None, istart=0): sources.extend(new_src) return sources + def _refit_islands_3D(self, group, stage, outerclip=None, istart=0): + """ + Do island refitting (priorized fitting) on a group of islands. + + Parameters + ---------- + group : list + A list of components grouped by island. + + stage : int + Refitting stage. + + outerclip : float + Ignored, placed holder for future development. + + istart : int + The starting island number. + + Returns + ------- + sources : list + List of sources (and islands). + """ + sources = [] + + data = self.img + rmsimg = self.rmsimg + + for inum, isle in enumerate(group, start=istart): + logger.debug("-=-") + logger.debug(f"input island = {isle[0].island}, {len(isle)} components") + + # set up the parameters for each of the sources within the island + i = 0 + params = lmfit.Parameters() + shape = data.shape + _, xmin, ymin = shape + xmax = ymax = 0 + + # island_mask = [] + src_valid_psf = None + # keep track of the sources that are actually being refit + # this may be a subset of all sources in the island + included_sources = [] + for src in isle: + # get the beam for this source + pixbeam = Beam(*self.wcshelper.get_psf_sky2pix(src.ra, src.dec)) + # find the right pixels from the ra/dec + source_x, source_y = self.wcshelper.sky2pix([src.ra, src.dec]) + source_x -= 1 + source_y -= 1 + x = int(round(source_x)) + y = int(round(source_y)) + + logger.debug(f"pixel location ({source_x:5.2f},{source_y:5.2f})") + # reject sources that are outside the image bounds, + # or which have nan data/rms values + if ( + not 0 <= x < shape[1] + or not 0 <= y < shape[2] + or not np.isfinite(data[0, x, y]) + or not np.isfinite(rmsimg[0, x, y]) + or pixbeam is None + ): + logger.debug( + f"Source ({src.island},{src.source}) not within usable region: skipping" + ) + continue + else: + # Keep track of the last source to have a valid psf + # so that we can use it later on + src_valid_psf = src + # determine the shape parameters in pixel values + _, _, sx, sy, theta = self.wcshelper.sky2pix_ellipse( + [src.ra, src.dec], src.a / 3600, src.b / 3600, src.pa + ) + sx *= FWHM2CC + sy *= FWHM2CC + + logger.debug( + f"Source shape [sky coords] {src.a:5.2f}x{src.b:5.2f}@{src.pa:05.2f}" + ) + logger.debug( + f"Source shape [pixel coords] {sx:4.2f}x{sy:4.2f}@{theta:05.2f}" + ) + + # choose a region that is 2x the major axis of the source, + # 4x semimajor axis a + width = 4 * sx + ywidth = int(round(width)) + 1 + xwidth = int(round(width)) + 1 + + # adjust the size of the island to include this source + xmin = min(xmin, max(0, x - xwidth / 2)) + ymin = min(ymin, max(0, y - ywidth / 2)) + xmax = max(xmax, min(shape[1], x + xwidth / 2 + 1)) + ymax = max(ymax, min(shape[2], y + ywidth / 2 + 1)) + + s_lims = [0.8 * min(sx, pixbeam.b * FWHM2CC), max(sy, sx) * 1.25] + + # Set up the parameters for the fit, including constraints + prefix = f"c{i}_" + params.add(prefix + "amp", value=src.peak_flux, vary=True) + # for now the xo/yo are locations within the main image, + # we correct this later + params.add( + prefix + "xo", + value=source_x, + min=source_x - sx / 2.0, + max=source_x + sx / 2.0, + vary=stage >= 2, + ) + params.add( + prefix + "yo", + value=source_y, + min=source_y - sy / 2.0, + max=source_y + sy / 2.0, + vary=stage >= 2, + ) + params.add( + prefix + "sx", + value=sx, + min=s_lims[0], + max=s_lims[1], + vary=stage >= 3, + ) + params.add( + prefix + "sy", + value=sy, + min=s_lims[0], + max=s_lims[1], + vary=stage >= 3, + ) + params.add( + prefix + "alpha", + value=-1, + min=-2, + max=1, + ) + params.add(prefix + "theta", value=theta, vary=stage >= 3) + params.add(prefix + "flags", value=0, vary=False) + params.add(prefix + "nu0", value=self.wcshelper.pix2freq(0), vary=False) + included_sources.append(src) + i += 1 + + # TODO: Allow this mask to be used in conjunction with + # the FWHM mask that is defined further on + + if i == 0: + logger.debug(f"No sources found in island {src.island}") + continue + params.add("components", value=i, vary=False) + # params.components = i + logger.debug(f" {i} components being fit") + # now we correct the xo/yo positions to be + # relative to the sub-image + logger.debug(f"xmxxymyx {xmin} {xmax} {ymin} {ymax}") + for i in range(int(params["components"].value)): + prefix = f"c{i}_" + # must update limits before the value as limits are + # enforced when the value is updated + params[prefix + "xo"].min -= xmin + params[prefix + "xo"].max -= xmin + params[prefix + "xo"].value -= xmin + params[prefix + "yo"].min -= ymin + params[prefix + "yo"].max -= ymin + params[prefix + "yo"].value -= ymin + # logger.debug(params) + # don't fit if there are no sources + if params["components"].value < 1: + logger.info(f"Island {src.island} has no components") + continue + + # this .copy() will stop us from modifying the parent region when + # we later apply our mask. + idata = data[:, int(xmin) : int(xmax), int(ymin) : int(ymax)].copy() + # now convert these back to indices within the idata region + # island_mask = np.array([(x-xmin, y-ymin) for x,y in island_mask]) + + allz, allx, ally = np.indices(idata.shape) + # mask to include pixels that are withn the FWHM + # of the sources being fit + mask_params = copy.deepcopy(params) + for i in range(int(mask_params["components"].value)): + prefix = f"c{i}_" + mask_params[prefix + "amp"].value = 1 + mask_model = nthreedgaussian_lmfit(mask_params) + mask = np.where(mask_model(allz.ravel(), allx.ravel(), ally.ravel()) <= 0.1) + mask = allz.ravel()[mask], allx.ravel()[mask], ally.ravel()[mask] + del mask_params + + idata[mask] = np.nan + + _, mx, my = np.where(np.isfinite(idata)) + non_nan_pix = len(mx) + total_pix = len(allx.ravel()) + logger.debug("island extracted:") + logger.debug(" x[{0}:{1}] y[{2}:{3}]".format(xmin, xmax, ymin, ymax)) + logger.debug(f"This is the idata: {idata}") + logger.debug(" max = {0}".format(np.nanmax(idata))) + logger.debug( + " total {0}, masked {1}, not masked {2}".format( + total_pix, total_pix - non_nan_pix, non_nan_pix + ) + ) + + # Check to see that each component has some data within + # the central 3x3 pixels of it's location + # If not then we don't fit that component + for i in range(int(params["components"].value)): + prefix = f"c{i}_" + # figure out a box around the center of this + cx, cy = ( + params[prefix + "xo"].value, + params[prefix + "yo"].value, + ) # central pixel coords + logger.debug(f" comp {i}") + logger.debug(f" x0, y0 {cx} {cy}") + xmx = int(round(np.clip(cx + 2, 0, idata.shape[1]))) + xmn = int(round(np.clip(cx - 1, 0, idata.shape[1]))) + ymx = int(round(np.clip(cy + 2, 0, idata.shape[2]))) + ymn = int(round(np.clip(cy - 1, 0, idata.shape[2]))) + square = idata[xmn:xmx, ymn:ymx] + # if there are no not-nan pixels in this region + # then don't vary any parameters + if not np.any(np.isfinite(square)): #! Not sure what to do here + logger.debug(f" not fitting component {i}") + params[prefix + "amp"].value = np.nan + for p in ["amp", "xo", "yo", "sx", "sy", "theta", "alpha"]: + params[prefix + p].vary = False + params[prefix + p].stderr = np.nan + # the above results in an error of -1 later on + params[prefix + "flags"].value |= flags.NOTFIT + + # determine the number of free parameters and + # if we have enough data for a fit + nfree = np.count_nonzero([params[p].vary for p in params.keys()]) + logger.debug(params) + if nfree < 1: + logger.debug(" Island has no components to fit") + result = DummyLM() + model = params + else: + if non_nan_pix < nfree: + logger.debug( + f"More free parameters {nfree} than available pixels {non_nan_pix}" + ) + if non_nan_pix >= params["components"].value: + logger.debug("Fixing all parameters except amplitudes") + for p in params.keys(): + if "amp" not in p: + params[p].vary = False + else: + logger.debug(" no not-masked pixels, skipping") + continue + + # do the fit + # if the pixel beam is not valid, then recalculate using + # the location of the last source to have a valid psf + if pixbeam is None: + if src_valid_psf is not None: + pixbeam = self.wcshelper.get_pixbeam( + src_valid_psf.ra, src_valid_psf.dec + ) + else: + logger.critical("Cannot determine pixel beam") + fac = 1 / np.sqrt(2) + if False: # TODO: Add this back alter -> self.docov: + C = Cmatrix( + mx, + my, + pixbeam.a * FWHM2CC * fac, + pixbeam.b * FWHM2CC * fac, + pixbeam.pa, + ) + B = Bmatrix(C) + else: + C = B = None + errs = np.nanmax( + rmsimg[:, int(xmin) : int(xmax), int(ymin) : int(ymax)] + ) + mask = np.where(np.isfinite(idata)) + freq_mapping = np.array( + [self.wcshelper.pix2freq(f) for f in range(idata.shape[0])] + ) + logger.debug(f"Freq mapping: {freq_mapping}") + result, _ = do_lmfit_3D(idata, params, B=B, freq_mapping=freq_mapping) + model = covar_errors( + result.params, idata, errs=errs, B=B, C=C + ) #! Check the covar_errors function + + # convert the results to a source object + offsets = (xmin, xmax, ymin, ymax) + # TODO allow for island fluxes in the refitting. + island_data = IslandFittingData( + inum, i=idata, offsets=offsets, doislandflux=False, scalars=(4, 4, None) + ) + new_src = self.result_to_components(result, model, island_data, src.flags) + + for ns, s in zip(new_src, included_sources): + # preserve the uuid so we can do exact + # matching between catalogs + ns.uuid = s.uuid + + # flag the sources as having been priorized + ns.flags |= flags.PRIORIZED + + # if the position wasn't fit then copy the errors + # from the input catalog + if stage < 2: + ns.err_ra = s.err_ra + ns.err_dec = s.err_dec + ns.flags |= flags.FIXED2PSF + + # if the shape wasn't fit then copy the errors + # from the input catalog + if stage < 3: + ns.err_a = s.err_a + ns.err_b = s.err_b + ns.err_pa = s.err_pa + sources.extend(new_src) + return sources + def _fit_island(self, island_data): """ Take an Island, do all the parameter estimation and fitting. @@ -1775,20 +2115,6 @@ def _fit_island(self, island_data): # get the real (sky) parameter errors model = covar_errors(result.params, idata, errs=errs, B=B, C=C) - if self.dobias and self.docov: - x, y = np.indices(idata.shape) - acf = elliptical_gaussian( - x, - y, - 1, - 0, - 0, - pixbeam.a * FWHM2CC * fac, - pixbeam.b * FWHM2CC * fac, - pixbeam.pa, - ) - bias_correct(model, idata, acf=acf * errs**2) - if not result.success: is_flag |= flags.FITERR @@ -1834,6 +2160,7 @@ def find_sources_in_image( docov=True, cube_index=None, progress=True, + threeD=False, ): """ Run the Aegean source finder. @@ -1930,8 +2257,8 @@ def find_sources_in_image( blank=blank, docov=docov, cube_index=cube_index, + as_cube=threeD, ) - logger.info( "beam = {0:5.2f}'' x {1:5.2f}'' at {2:5.2f}deg".format( self.beam.a * 3600, @@ -1942,8 +2269,8 @@ def find_sources_in_image( # stop people from doing silly things. if outerclip > innerclip: outerclip = innerclip - logger.info("seedclip={0}".format(innerclip)) - logger.info("floodclip={0}".format(outerclip)) + logger.info(f"seedclip={innerclip}") + logger.info(f"floodclip={outerclip}") islands = find_islands( im=self.img, @@ -1954,7 +2281,7 @@ def find_sources_in_image( region=self.region, wcs=self.wcshelper, ) - logger.info("Found {0} islands".format(len(islands))) + logger.info(f"Found {len(islands)} islands") logger.info("Begin fitting") island_group = [] @@ -1969,7 +2296,7 @@ def find_sources_in_image( # ignore empty islands # This should now be impossible to trigger if not np.any(np.isfinite(i)): - logger.warning("Empty island detected, this should be imposisble.") + logger.warning("Empty island detected, this should be impossible.") continue isle_num += 1 scalars = (innerclip, outerclip, max_summits) @@ -1979,12 +2306,22 @@ def find_sources_in_image( # now fit all the islands sources = [] - with tqdm( - total=isle_num, desc="Fitting Islands:", disable=not progress - ) as pbar: - for i in island_group: + if MPI_AVAIL: + + size = MPI.COMM_WORLD.Get_size() + rank = MPI.COMM_WORLD.Get_rank() + + block_size = len(island_group) // size + + start = block_size * rank + end = block_size * (rank + 1) + + # highest rank will deal with the remainder (which is < size) + if rank == size: + end = None + + for i in island_group[start:end]: try: - pbar.update(1) srcs = self._fit_island(i) except AegeanNaNModelError: continue @@ -1997,6 +2334,25 @@ def find_sources_in_image( continue sources.append(src) + else: + with tqdm( + total=isle_num, desc="Fitting Islands:", disable=not progress + ) as pbar: + for i in island_group: + try: + pbar.update(1) + srcs = self._fit_island(i) + except AegeanNaNModelError: + continue + # update bar as each individual island is fit + for src in srcs: + # ignore sources that we have been told to ignore + if (src.peak_flux > 0 and nopositive) or ( + src.peak_flux < 0 and nonegative + ): + continue + sources.append(src) + # Write the output to the output file if outfile: print( @@ -2146,6 +2502,235 @@ def priorized_fit_islands( # reject sources with missing params ok = True for param in ["ra", "dec", "peak_flux", "a", "b", "pa"]: + if np.isnan(getattr(input_sources[0], param)): + logger.info(f"Source 0, is missing param '{param}'") + ok = False + if not ok: + logger.error("Missing parameters! Not fitting.") + logger.error("Maybe your table is missing or mis-labeled columns?") + return [] + del ok + + # Do the resizing + logger.info(f"{len(input_sources)} sources in catalog") + sources = cluster.resize(input_sources, ratio=ratio, wcshelper=self.wcshelper) + logger.info(f"{len(sources)} sources accepted") + + if len(sources) < 1: + logger.debug("No sources accepted for priorized fitting") + return [] + + # compute eps if it's not defined + if regroup_eps is None: + # s.a is in arcsec but we assume regroup_eps is in arcmin + regroup_eps = 4 * np.mean([s.a / 60 for s in sources]) + # convert regroup_eps into a value appropriate for a cartesian measure + regroup_eps = np.sin(np.radians(regroup_eps / 60)) + input_sources = sources + # redo the grouping if required + if doregroup: + # TODO: scale eps in some appropriate manner + groups = regroup_dbscan(input_sources, eps=regroup_eps) + else: + groups = list(island_itergen(input_sources)) + + logger.info("Begin fitting") + + island_groups = [] # will be a list of groups of islands + island_group = [] # will be a list of islands + group_size = 20 + + # ? Why split the islands into groups if we are going to append them all + # ? to the island_groups list? + for island in groups: + island_group.append(island) + # If the island group is full queue it for the subprocesses to fit + if len(island_group) >= group_size: + island_groups.append(island_group) + island_group = [] + # The last partially-filled island group also needs to + # be queued for fitting + if len(island_group) > 0: + island_groups.append(island_group) + + sources = [] + with tqdm( + total=len(island_groups), + desc="Refitting Island Groups", + disable=not progress, + ) as pbar: + for i, g in enumerate(island_groups): + srcs = self._refit_islands(g, stage, outerclip, istart=i) + # update bar as each individual island is fit + pbar.update(1) + sources.extend(srcs) + + sources = sorted(sources) + + # Write the output to the output file + if outfile: + print( + header.format("{0}-({1})".format(__version__, __date__), filename), + file=outfile, + ) + print(ComponentSource.header, file=outfile) + for source in sources: + print(str(source), file=outfile) + + logger.info("fit {0} components".format(len(sources))) + self.sources.extend(sources) + return sources + + def priorized_fit_islands_3D( + self, + filename, + catalogue, + # hdu_index=0, + outfile=None, + bkgin=None, + rmsin=None, + cores=1, + # rms=None, + # bkg=None, + # beam=None, + # imgpsf=None, + # catpsf=None, + stage=3, + ratio=None, + outerclip=3, + doregroup=True, + regroup_eps=None, + docov=False, + # cube_index=None, + progress=True, + ): + """ + Take an input catalog, image, and background/noise images + fit the flux and ra/dec for each of the given sources. #! Delete the parameters that are not used in the function + + Parameters + ---------- + filename : str or HDUList + Image filename or HDUList. + + catalogue : str or list + Input catalogue file name or list of ComponentSource objects. + + hdu_index : int + The index of the FITS HDU (extension). + + outfile : str + file for printing catalog + (NOT a table, just a text file of my own design) + + rmsin, bkgin : str or HDUList + Filename or HDUList for the noise and background images. + If either are None, then it will be calculated internally. + + cores : int + Number of CPU cores to use. None means all cores. + + rms : float + Use this rms for the entire image + (will also assume that background is 0) + + beam : (float, float, float) + (major, minor, pa) representing the synthesised beam (degrees). + Replaces whatever is given in the FITS header. + If the FITS header has no BMAJ/BMIN then this is required. + + imgpsf : str or HDUList + Filename or HDUList for a psf image. + + catpsf : str or HDUList + Filename or HDUList for the catalogue psf image. + + stage : int + Refitting stage + + ratio : float + If not None - ratio of image psf to catalog psf, + otherwise interpret from catalogue or image if possible + + outerclip : float + The flood (outer) clipping level (sigmas). + + doregroup : bool, Default=True + Relabel all the islands/groups to ensure that nearby + components are jointly fit. + + regroup_eps: float, Default=None + The linking parameter for regouping. Components that are + closer than this distance (in arcmin) will be jointly fit. + If NONE, then use 4x the average source major axis size (after + rescaling if required). + + docov : bool, Default=True + If True then include covariance matrix in the fitting process. + + cube_index : int, Default=None + For image cubes, slice determines which slice is used. + + progress : bool, Default=True + Show a progress bar when fitting island groups + + Returns + ------- + sources : list + List of sources measured. + + """ + + from AegeanTools.cluster import regroup_dbscan + + self.load_globals( + filename, + # hdu_index=hdu_index, + bkgin=bkgin, + rmsin=rmsin, + # beam=beam, + verb=True, + # rms=rms, + bkg=bkgin, + cores=cores, + do_curve=False, + # psf=imgpsf, + docov=docov, + # cube_index=cube_index, + as_cube=True, + ) + + # load the table and convert to an input source list + if isinstance(catalogue, str): + input_table = load_table(catalogue) + input_sources = np.array( + table_to_source_list(input_table, src_type=ComponentSource3D) + ) + else: + input_sources = np.array(catalogue) # Assumed to be a list of objects + + nu0 = self.wcshelper.pix2freq(0) + + for src in input_sources: + src.alpha = -1 + src.nu0 = nu0 + + if len(input_sources) < 1: + logger.debug("No input sources for priorized fitting") + return [] + + # reject sources with missing params + ok = True + for param in [ + "ra", + "dec", + "peak_flux", + "a", + "b", + "pa", + "alpha", + "nu0", + ]: #! Added freq and alpha if np.isnan(getattr(input_sources[0], param)): logger.info("Source 0, is missing param '{0}'".format(param)) ok = False @@ -2165,7 +2750,8 @@ def priorized_fit_islands( return [] # compute eps if it's not defined - if regroup_eps is None: + if regroup_eps is None: #! How do we handle this as the elliptical beam is + #! not defined in the 3D case # s.a is in arcsec but we assume regroup_eps is in arcmin regroup_eps = 4 * np.mean([s.a / 60 for s in sources]) # convert regroup_eps into a value appropriate for a cartesian measure @@ -2202,7 +2788,9 @@ def priorized_fit_islands( disable=not progress, ) as pbar: for i, g in enumerate(island_groups): - srcs = self._refit_islands(g, stage, outerclip, istart=i) + srcs = self._refit_islands_3D( + g, stage, outerclip, istart=i + ) #! <- Update it after the function is implemented # update bar as each individual island is fit pbar.update(1) sources.extend(srcs) diff --git a/AegeanTools/wcs_helpers.py b/AegeanTools/wcs_helpers.py index 87bdd645..34497297 100644 --- a/AegeanTools/wcs_helpers.py +++ b/AegeanTools/wcs_helpers.py @@ -63,7 +63,9 @@ class WCSHelper(object): The WCS object for the psf map """ - def __init__(self, wcs, beam, pixscale, refpix, psf_file=None): + def __init__( + self, wcs, beam, pixscale, refpix, psf_file=None, freq_wcs=(None, None, None) + ): """ Parameters ---------- @@ -83,6 +85,8 @@ def __init__(self, wcs, beam, pixscale, refpix, psf_file=None): Filename for a psf map """ self.wcs = wcs + self.CDELT3, self.CRVAL3, self.CRPIX3 = freq_wcs + self.has_freq = self.CDELT3 is not None self.beam = ( # the beam as per the fits header (at the reference coordiante) beam @@ -176,7 +180,11 @@ def from_header(cls, header, beam=None, psf_file=None): _, pixscale = get_pixinfo(header) refpix = (header["CRPIX1"], header["CRPIX2"]) - return cls(wcs, beam, pixscale, refpix, psf_file=psf_file) + if "CRVAL3" in header: + freq_wcs = (header["CDELT3"], header["CRVAL3"], header["CRPIX3"]) + else: + freq_wcs = (None, None, None) + return cls(wcs, beam, pixscale, refpix, psf_file=psf_file, freq_wcs=freq_wcs) @classmethod def from_file(cls, filename, beam=None, psf_file=None): @@ -243,6 +251,48 @@ def sky2pix(self, pos): # wcs and python have opposite ideas of x/y return [pixel[0][1], pixel[0][0]] + def pix2freq( + self, pos_z + ): #! Write a test case for this, afterwards create a pull request + """ + Convert z coordinates into frequency. + + Parameters + ---------- + pos_z : float + The z axis pixel coordinates. + + Returns + ------- + frequency : float + The frequency in MHz + + """ + + return (self.CDELT3 / 1_000_000) * (pos_z - self.CRPIX3) + ( + self.CRVAL3 / 1_000_000 + ) + + def freq2pix(self, freq): #! Write a test case for this + """ + Convert frequency into pixel coordinates. + + Parameters + ---------- + freq : float + The frequency in MHz. + + Returns + ------- + pixel : (float, float) + The (x,y) pixel coordinates + + """ + + return (freq - (self.CRVAL3 / 1_000_000)) / ( + self.CDELT3 / 1_000_000 + ) + self.CRPIX3 + def psf_sky2pix(self, pos): """ Convert sky coordinates into pixel coordinates. diff --git a/CHANGELOG.md b/CHANGELOG.md index 27416081..40189ed1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,24 @@ General - Supported (tested) versions of python are now 3.10-3.13, python 3.8 and 3.9 are no longer supported (but may work) - Updated the minimum version of libraries required for AegeanTools this co-incides with the change in supported versions of python +### Jan 15 2025 + +General +- All command line scripts now use `configargparse`, so that users can supply a + configuration file via `--config` instead of a long list of command line arguments. + +Aegean +- Output format '.crtf' is now supported for catalogues. +- See [Casa Docs](https://casaguides.nrao.edu/index.php/CASA_Region_Format) for details on this format. + +BANE +- BANE is now able input/ouput data cubes in their entirety. +- Default now is that if the input is a cube then the output will be a cube. +- Use `--slice` to choose a single plane from a cube and output a 2d image. + +AeRes +- Will now work with 3d image cubes and catalogues with a spectral index + ### Dec 09 2024 General diff --git a/Developer_Documentation.md b/Developer_Documentation.md new file mode 100644 index 00000000..5ba10372 --- /dev/null +++ b/Developer_Documentation.md @@ -0,0 +1,81 @@ + +# Aegean + +## Table of Contents +1. [Introduction](#introduction) +2. [Installation](#installation) +3. [Usage](#usage) +4. [Testing](#testing) + +## Introduction +Aegean is a source finding program designed to find and characterise compact continuum radio sources. + +## Installation +To install Aegean, follow the steps below. +### Prerequisites +- Python +- OpenMPI +- Git + +### Steps +1. Clone the repository. + ``` + git clone https://github.com/PaulHancock/Aegean.git + ``` +2. Navigate to the project directory. + ``` + cd AEGEAN + ``` +3. Install. + ``` + pip install . + ``` + +## Usage +Aegean can be used as a command line tool. The command line tool can be used to find sources in a FITS image and output the results in a variety of formats. The command line tool can also be used to load background and rms files that might have been previously created by BANE. +### Basic Usage +``` +aegean path/to/file +``` + +### Advanced Usage + +There are many options available for the Aegean command line tool. These include: +-- table : for specifying the format of the output +-- autoload : for automatically loading the background and rms files that might have been previously created by BANE + +### Example +``` +aegean path/to/file --table table_name.csv,table_name.fits --autoload +``` + +### MPI usage +To use the you can use the following bash script with minor modifications to suit your needs. +``` +#!/bin/bash -l +#SBATCH --account=pawsey000X +#SBATCH --partition=work +#SBATCH --time=00:40:00 +#SBATCH --nodes=16 +#SBATCH --cpus-per-task=128 +#SBATCH --ntasks-per-node=1 +#SBATCH --mem=32gb +#SBATCH --output=/scratch/pawsey000X/user//.o%A +#SBATCH --error=/scratch/pawsey000X/user//.e%A + +cd /scratch/pawsey000X// +source ../env/bin/activate +module load cray-mpich/8.1.25 +module list +echo $HOSTNAME +srun aegean --table ${SLURM_JOB_ID}.fits --autoload +``` + +## Testing +To run tests for the new addSuffix function, run the following file: +``` +pytest tests/unit/test_addSuffix.py +``` + +## References +For further user guide and documentation, please refer to the Aegean Github repository at https://github.com/PaulHancock/Aegean diff --git a/tests/integration/test_script_AeRes.py b/tests/integration/test_script_AeRes.py index fc9d5c93..3c2ec8c5 100644 --- a/tests/integration/test_script_AeRes.py +++ b/tests/integration/test_script_AeRes.py @@ -4,8 +4,12 @@ import os import sys -image_SIN = "tests/test_files/1904-66_SIN.fits" -catfile = "tests/test_files/1904_comp.fits" +image_2d = "tests/test_files/synthetic_image.fits" +image_3d = "tests/test_files/synthetic_cube.fits" +cat_with_alpha = "tests/test_files/synthetic_cat_with_alpha_comp.fits" +cat_no_alpha = "tests/test_files/synthetic_cat_no_alpha_comp.fits" +# image_SIN = "tests/test_files/1904-66_SIN.fits" +# catfile = "tests/test_files/1904_comp.fits" tempfile = "dlme" @@ -14,31 +18,66 @@ def test_help(): AeRes.main() -def test_nocat(): +def test_cat_with_alpha_2d_image(): try: + sys.argv = ["", "-c", cat_with_alpha, "-f", image_2d, "-r", tempfile] AeRes.main() + finally: + if os.path.exists(tempfile): + os.remove(tempfile) + else: + raise Exception("Output file not created") + - sys.argv = ["", "-c", catfile] +def test_cat_with_alpha_3d_image(): + try: + sys.argv = ["", "-c", cat_with_alpha, "-f", image_3d, "-r", tempfile] AeRes.main() + finally: + if os.path.exists(tempfile): + os.remove(tempfile) + else: + raise Exception("Output file not created") + - sys.argv = ["", "-c", catfile, "-f", image_SIN] +def test_cat_no_alpha_2d_image(): + try: + sys.argv = ["", "-c", cat_no_alpha, "-f", image_2d, "-r", tempfile] AeRes.main() + finally: + if os.path.exists(tempfile): + os.remove(tempfile) + else: + raise Exception("Output file not created") + - sys.argv = [ - "", - "-c", - catfile, - "-f", - image_SIN, - "-r", - tempfile, - "-m", - tempfile + "_model", - ] +def test_cat_no_alpha_3d_image(): + try: + sys.argv = ["", "-c", cat_no_alpha, "-f", image_3d, "-r", tempfile] AeRes.main() finally: - os.remove(tempfile) - os.remove(tempfile + "_model") + if os.path.exists(tempfile): + os.remove(tempfile) + else: + raise Exception("Output file not created") + + +def test_nocat(): + try: + sys.argv = [""] # no arguments + AeRes.main() + + sys.argv = ["", "-c", cat_no_alpha] + AeRes.main() + + sys.argv = ["", "-c", cat_no_alpha, "-f", image_2d] + AeRes.main() + + finally: + if os.path.exists(tempfile): + os.remove(tempfile) + if os.path.exists(tempfile + "_model"): + os.remove(tempfile + "_model") def test_configfile(): diff --git a/tests/integration/test_script_aegean.py b/tests/integration/test_script_aegean.py index cf60cc03..6958b4ee 100644 --- a/tests/integration/test_script_aegean.py +++ b/tests/integration/test_script_aegean.py @@ -67,8 +67,8 @@ def test_aux_images(): def test_find(): - sys.argv = ["", image_SIN, "--table", "test"] - aegean.main() + # sys.argv = ["", image_SIN, "--table", "test"] + # aegean.main() sys.argv = ["", image_SIN, "--out", "stdout"] aegean.main() @@ -105,6 +105,30 @@ def test_priorized(): aegean.main() +def test_priorized3D(): + sys.argv = [ + "", + "tests/test_files/synthetic_cube.fits", + "--noise", + "tests/test_files/synthetic_cube_rms.fits", + "--background", + "tests/test_files/synthetic_cube_bkg.fits", + "--table", + "output.csv", + "--3d", + "--priorized", + "1", + "--input", + "tests/test_files/synthetic_cat_no_alpha_comp.fits", + ] + aegean.main() + + if not os.path.exists("output_comp.csv"): + raise AssertionError("output file not created") + else: + os.remove("output_comp.csv") + + def test_configfile(): cfile = "aegean.ini" try: diff --git a/tests/test_files/synthetic_cat_no_alpha_comp.fits b/tests/test_files/synthetic_cat_no_alpha_comp.fits new file mode 100644 index 00000000..37fb8e95 --- /dev/null +++ b/tests/test_files/synthetic_cat_no_alpha_comp.fits @@ -0,0 +1,4 @@ +SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 149 / length of dimension 1 NAXIS2 = 102 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 27 / number of table fields TTYPE1 = 'island ' TFORM1 = 'J ' TTYPE2 = 'source ' TFORM2 = 'J ' TTYPE3 = 'background' TFORM3 = 'E ' TTYPE4 = 'local_rms' TFORM4 = 'E ' TTYPE5 = 'ra_str ' TFORM5 = '8A ' TTYPE6 = 'dec_str ' TFORM6 = '9A ' TTYPE7 = 'ra ' TFORM7 = 'E ' TTYPE8 = 'err_ra ' TFORM8 = 'E ' TTYPE9 = 'dec ' TFORM9 = 'E ' TTYPE10 = 'err_dec ' TFORM10 = 'E ' TTYPE11 = 'peak_flux' TFORM11 = 'E ' TTYPE12 = 'err_peak_flux' TFORM12 = 'E ' TTYPE13 = 'int_flux' TFORM13 = 'E ' TTYPE14 = 'err_int_flux' TFORM14 = 'E ' TTYPE15 = 'a ' TFORM15 = 'J ' TTYPE16 = 'err_a ' TFORM16 = 'E ' TTYPE17 = 'b ' TFORM17 = 'J ' TTYPE18 = 'err_b ' TFORM18 = 'E ' TTYPE19 = 'pa ' TFORM19 = 'J ' TTYPE20 = 'err_pa ' TFORM20 = 'E ' TTYPE21 = 'flags ' TFORM21 = 'J ' TTYPE22 = 'residual_mean' TFORM22 = 'E ' TTYPE23 = 'residual_std' TFORM23 = 'E ' TTYPE24 = 'uuid ' TFORM24 = '36A ' TTYPE25 = 'psf_a ' TFORM25 = 'J ' TTYPE26 = 'psf_b ' TFORM26 = 'J ' TTYPE27 = 'psf_pa ' TFORM27 = 'J ' END ?00:00:00-00:00:00A�ff��x�??�������a316b0fe-a761-452f-ba74-2d46085730a7?00:00:00-00:00:00A�ff��wI��?�?�������1f6d5a44-8df3-49be-a55c-52f3e0e8b29a?00:00:00-00:00:00A�ff��v���?�?�������2b20164c-f2f6-45f0-9c1b-259a9656f869?00:00:00-00:00:00A�ff��u���@?�������e28dafb5-e62f-4a66-a501-ecef854c8f75?00:00:00-00:00:00A�ff��u'��@ ?�������87b20128-6092-4256-89c7-3c506ea702da?00:00:00-00:00:00A�ff��tq��@@?�������81be6254-5e6f-4df6-a1a7-95ce5ad969e4?00:00:00-00:00:00A�ff��s���@`?�������eb12d597-a100-4c55-abdc-caa963f5f900?00:00:00-00:00:00A�ff��s��@�?�������77600b0b-8282-424e-b8d8-e934aa5da22e?00:00:00-00:00:00A�ff��rO��@�?�������646fc280-5c36-44a9-bd85-1a1e53355e9b ?00:00:00-00:00:00A�ff��q���@�?�������f1f469cd-dd83-48f1-b3db-d248a827f369 +?00:00:00-00:00:00A�I���x�@�?�������68a7e0a4-0eea-4b1d-b0e8-7e586ea2d248 ?00:00:00-00:00:00A�I���wI��@�?�������e01bae53-0577-4143-bfb4-655afb353a35 ?00:00:00-00:00:00A�I���v���@�?�������6b06afa2-327d-4d63-8958-7655b4e95cd0 ?00:00:00-00:00:00A�I���u���@�?�������7af01cca-8756-4b93-a919-d35f8c8c383d?00:00:00-00:00:00A�I���u'��@�?�������db52be3b-1daa-4a3c-b1ca-cdf4762e42a8?00:00:00-00:00:00A�I���tq��A?�������8b853852-f668-430b-9a25-69e5117389e0?00:00:00-00:00:00A�I���s���A?�������667ff57b-1f08-4032-8518-0fd20e30fa11?00:00:00-00:00:00A�I���s��A?�������950eb7eb-2f22-49e9-b859-c46b11684d89?00:00:00-00:00:00A�I���rO��A?�������60b27ba8-7708-431c-ba8f-9e1dfa99c800?00:00:00-00:00:00A�I���q���A ?�������0edf6eff-84db-4166-97c3-d0ed04dbd889?00:00:00-00:00:00A�-���x�A(?�������6182968a-469d-4ee7-b65a-bc645e8c1a0f?00:00:00-00:00:00A�-���wI��A0?�������57520487-510e-4295-8e6c-f8ea5d214f78?00:00:00-00:00:00A�-���v���A8?�������79ddf43f-18f2-4efa-9ebe-9800f161d546?00:00:00-00:00:00A�-���u���A@?�������34495849-1cce-4086-b86b-815d38455285?00:00:00-00:00:00A�-���u'��AH?�������58672cfd-685a-4f63-9c75-73e45634f920?00:00:00-00:00:00A�-���tq��AP?�������f89c9984-6e07-4e4b-a804-948457450d80?00:00:00-00:00:00A�-���s���AX?�������36059ee1-aca1-4547-8ca1-c581e1dd51be?00:00:00-00:00:00A�-���s��A`?�������14e00510-b4df-41aa-ae7e-1c3f0ea90238?00:00:00-00:00:00A�-���rO��Ah?�������9b5eb77f-c7f0-4791-a1ca-9ab9696202d2?00:00:00-00:00:00A�-���q���Ap?�������5c600cd6-5138-4600-b78f-9a242bdfffe0?00:00:00-00:00:00A���x�Ax?�������89c22dbb-6493-46ab-b3d2-8e2690276581?00:00:00-00:00:00A���wI��A�?�������a4ae120a-f7db-4fb1-b4ab-b2d0ed005de9 ?00:00:00-00:00:00A���v���A�?�������fab9f9bb-7371-4e43-a14e-2cf0cceac294!?00:00:00-00:00:00A���u���A�?�������780667d1-da09-4d1d-8ebb-3f5cc8cfa45a"?00:00:00-00:00:00A���u'��A�?�������d4a81f44-4cc1-4e3d-9fd3-2ce19205fd0c#?00:00:00-00:00:00A���tq��A�?�������ca40f757-611d-4a87-bce7-d7ef47270ed1$?00:00:00-00:00:00A���s���A�?�������a6aa3c28-8f36-4bda-aa0a-fa26c1ad7bf4%?00:00:00-00:00:00A���s��A�?�������5ffdcbf5-0392-466b-96de-91e216563d17&?00:00:00-00:00:00A���rO��A�?�������f0e2a560-f1c1-4b7a-b321-50a2f8ed2e4f'?00:00:00-00:00:00A���q���A�?�������e8603036-0a53-4500-b14d-c0fa40e8485c(?00:00:00-00:00:00A�����x�A�?�������5f0a3f2d-7c61-4a09-8e6e-ee99225d9cad)?00:00:00-00:00:00A�����wI��A�?�������f89946a0-36f6-4387-a842-a6a5a3edb30f*?00:00:00-00:00:00A�����v���A�?�������f3487366-bc26-420a-b244-9d0367650a2b+?00:00:00-00:00:00A�����u���A�?�������16378236-85b9-4977-b072-935e7ef39e6d,?00:00:00-00:00:00A�����u'��A�?�������02e8f848-2489-4078-a78f-1ead1ae4ea29-?00:00:00-00:00:00A�����tq��A�?�������3ef89b8b-2e29-4a5c-818b-94a424a0c809.?00:00:00-00:00:00A�����s���A�?�������28430f68-2c41-49b1-b04a-d7a35e5d3b7d/?00:00:00-00:00:00A�����s��A�?�������52133824-2182-4aec-94fc-76b69dfb59680?00:00:00-00:00:00A�����rO��A�?�������21f7635f-e7ce-4824-8b4f-9854c96b372c1?00:00:00-00:00:00A�����q���A�?�������367fac53-01ce-48aa-b95b-84656811fec02?00:00:00-00:00:00A��.��x�A�?�������73230a93-8365-4424-affd-1c3427421a2a3?00:00:00-00:00:00A��.��wI��A�?�������cebb773c-0cd9-46d1-af50-cccba9bf87e14?00:00:00-00:00:00A��.��v���A�?�������9651c0ed-e85b-453b-9ba2-c1ec2e2ea6985?00:00:00-00:00:00A��.��u���A�?�������aad5962e-e859-4fd5-9434-57ff77d35a7a6?00:00:00-00:00:00A��.��u'��A�?�������9bfc96e9-a4f3-482c-8dc8-0586d575df417?00:00:00-00:00:00A��.��tq��A�?�������caf02810-f06e-42ec-89a8-a4a5009b451a8?00:00:00-00:00:00A��.��s���A�?�������145ba62b-f743-478c-848b-4a53ef64c3159?00:00:00-00:00:00A��.��s��A�?�������bf50671e-fbd7-43dc-ab36-b647b6fc46e1:?00:00:00-00:00:00A��.��rO��A�?�������ea33081f-49e8-4098-aa72-b4ecac716a11;?00:00:00-00:00:00A��.��q���A�?�������d79633be-9ba2-4b48-88d6-0fd636ee7de5<?00:00:00-00:00:00A���x�A�?�������ad436a70-730e-49b1-8353-b3caeda973c9=?00:00:00-00:00:00A���wI��A�?�������66c233bc-2511-4ebd-8150-38611a5bda13>?00:00:00-00:00:00A���v���A�?�������8b6a8eff-d95e-4892-903d-7bf620531b7f??00:00:00-00:00:00A���u���B?�������f711c957-d04a-4c5e-8fab-9e96df4358e3@?00:00:00-00:00:00A���u'��B?�������ec83c28a-8dbb-4829-b1da-401850751c41A?00:00:00-00:00:00A���tq��B?�������744c0576-3576-4558-9f9d-784229fd941dB?00:00:00-00:00:00A���s���B?�������64dddf71-5230-4357-9a44-90d4abcba28fC?00:00:00-00:00:00A���s��B?�������96e076d6-79a1-45c8-b200-a292154c3764D?00:00:00-00:00:00A���rO��B +?�������07b0535d-52b5-4ade-8065-1269ba9cfed8E?00:00:00-00:00:00A���q���B ?�������4fed5d2c-0790-46cb-8694-49fa23d39610F?00:00:00-00:00:00A��J��x�B?�������35c72544-5b11-4d9b-aced-fc3c0e2edd18G?00:00:00-00:00:00A��J��wI��B?�������7da4e033-e46a-412a-8170-b8d0cda23a2cH?00:00:00-00:00:00A��J��v���B?�������3e533902-00b7-44d6-a4bd-a3e40e7945a4I?00:00:00-00:00:00A��J��u���B?�������4454e3cd-2100-4893-afd3-abf81ccd9342J?00:00:00-00:00:00A��J��u'��B?�������c4690174-5dc7-4b3e-95ee-e935c37f0012K?00:00:00-00:00:00A��J��tq��B?�������19fc1cdf-2d75-4de0-9fe2-e38d8c2af567L?00:00:00-00:00:00A��J��s���B?�������3369e286-1a63-4658-9588-2aa9e7cd4045M?00:00:00-00:00:00A��J��s��B?�������b0535e3c-fa34-4b8e-9fef-8c49ef5d5f78N?00:00:00-00:00:00A��J��rO��B?�������b339ca46-1583-4609-b0b0-a579f610426cO?00:00:00-00:00:00A��J��q���B ?�������555e05d0-3bf7-4c05-af92-4d83c471450aP?00:00:00-00:00:00A�����x�B"?�������661bdaa9-4359-42c0-ada4-612dcc8477caQ?00:00:00-00:00:00A�����wI��B$?�������1bb6ba57-f934-405e-86ad-4718cc3c83a2R?00:00:00-00:00:00A�����v���B&?�������fb125f3a-021c-4c05-a785-c16f9c382aebS?00:00:00-00:00:00A�����u���B(?�������eb86c23b-ccc1-4982-9dde-973e854a8da8T?00:00:00-00:00:00A�����u'��B*?�������c5cc8007-5b91-4ed4-9e95-44417e279a78U?00:00:00-00:00:00A�����tq��B,?�������bff9fbeb-79ca-46f6-8d43-49378bf214a4V?00:00:00-00:00:00A�����s���B.?�������481fbbe3-693f-4f4c-8b94-7ae4a0e2d028W?00:00:00-00:00:00A�����s��B0?�������dfa8a200-82f4-47c5-a374-cc8da78256fbX?00:00:00-00:00:00A�����rO��B2?�������ad6fde81-a22f-4b70-99a1-7511e041aea0Y?00:00:00-00:00:00A�����q���B4?�������26baedaf-c325-486a-bd5f-b61cf9552d9fZ?00:00:00-00:00:00A�ff��x�B6?�������f057383d-5d97-4f34-877a-6e4da1cf9428[?00:00:00-00:00:00A�ff��wI��B8?�������5a046e38-7d36-44ee-a8ef-9deb1f4141b0\?00:00:00-00:00:00A�ff��v���B:?�������f8561e86-a5f6-414c-bccf-e14037b72a84]?00:00:00-00:00:00A�ff��u���B<?�������3a1cc473-e29a-419c-97eb-57916e32b985^?00:00:00-00:00:00A�ff��u'��B>?�������db01f5a3-7a64-4811-a187-233375343e20_?00:00:00-00:00:00A�ff��tq��B@?�������320da9a0-a552-4127-8331-4e5647638adc`?00:00:00-00:00:00A�ff��s���BB?�������8dcea419-9aeb-434b-a236-adfdd17e6d4ba?00:00:00-00:00:00A�ff��s��BD?�������472372e2-bc75-41d1-a68a-c4b08e1350efb?00:00:00-00:00:00A�ff��rO��BF?�������e3a439fd-9419-42a9-abac-5606eb03cc02c?00:00:00-00:00:00A�ff��q���BH?�������911583a9-7ad0-4fce-9854-c1e9116fe956d��00:00:00-00:00:00A� +=��p���@@���Z�Z����f4dcfd49-0bc9-48f5-9325-c46ee4c86f44e��00:00:00-00:00:00A�=q��q���B ��������e15fcaa7-4aca-4004-bec3-9e668366f2e4 \ No newline at end of file diff --git a/tests/test_files/synthetic_cat_with_alpha_comp.fits b/tests/test_files/synthetic_cat_with_alpha_comp.fits new file mode 100644 index 00000000..c7d1b91c --- /dev/null +++ b/tests/test_files/synthetic_cat_with_alpha_comp.fits @@ -0,0 +1,5 @@ +SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 157 / length of dimension 1 NAXIS2 = 102 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 29 / number of table fields TTYPE1 = 'island ' TFORM1 = 'J ' TTYPE2 = 'source ' TFORM2 = 'J ' TTYPE3 = 'background' TFORM3 = 'E ' TTYPE4 = 'local_rms' TFORM4 = 'E ' TTYPE5 = 'ra_str ' TFORM5 = '8A ' TTYPE6 = 'dec_str ' TFORM6 = '9A ' TTYPE7 = 'ra ' TFORM7 = 'E ' TTYPE8 = 'err_ra ' TFORM8 = 'E ' TTYPE9 = 'dec ' TFORM9 = 'E ' TTYPE10 = 'err_dec ' TFORM10 = 'E ' TTYPE11 = 'peak_flux' TFORM11 = 'E ' TTYPE12 = 'err_peak_flux' TFORM12 = 'E ' TTYPE13 = 'int_flux' TFORM13 = 'E ' TTYPE14 = 'err_int_flux' TFORM14 = 'E ' TTYPE15 = 'a ' TFORM15 = 'J ' TTYPE16 = 'err_a ' TFORM16 = 'E ' TTYPE17 = 'b ' TFORM17 = 'J ' TTYPE18 = 'err_b ' TFORM18 = 'E ' TTYPE19 = 'pa ' TFORM19 = 'J ' TTYPE20 = 'err_pa ' TFORM20 = 'E ' TTYPE21 = 'flags ' TFORM21 = 'J ' TTYPE22 = 'residual_mean' TFORM22 = 'E ' TTYPE23 = 'residual_std' TFORM23 = 'E ' TTYPE24 = 'uuid ' TFORM24 = '36A ' TTYPE25 = 'psf_a ' TFORM25 = 'J ' TTYPE26 = 'psf_b ' TFORM26 = 'J ' TTYPE27 = 'psf_pa ' TFORM27 = 'J ' TTYPE28 = 'alpha ' TFORM28 = 'E ' TTYPE29 = 'nu0 ' TFORM29 = 'E ' END ?00:00:00-00:00:00A�ff��x�??�������a316b0fe-a761-452f-ba74-2d46085730a7>%qM~�?00:00:00-00:00:00A�ff��wI��?�?�������1f6d5a44-8df3-49be-a55c-52f3e0e8b29a>�x�M~�?00:00:00-00:00:00A�ff��v���?�?�������2b20164c-f2f6-45f0-9c1b-259a9656f869??�!M~�?00:00:00-00:00:00A�ff��u���@?�������e28dafb5-e62f-4a66-a501-ecef854c8f75����M~�?00:00:00-00:00:00A�ff��u'��@ ?�������87b20128-6092-4256-89c7-3c506ea702da�D��M~�?00:00:00-00:00:00A�ff��tq��@@?�������81be6254-5e6f-4df6-a1a7-95ce5ad969e4���'M~�?00:00:00-00:00:00A�ff��s���@`?�������eb12d597-a100-4c55-abdc-caa963f5f900��hM~�?00:00:00-00:00:00A�ff��s��@�?�������77600b0b-8282-424e-b8d8-e934aa5da22e?%)�M~�?00:00:00-00:00:00A�ff��rO��@�?�������646fc280-5c36-44a9-bd85-1a1e53355e9b?��M~� ?00:00:00-00:00:00A�ff��q���@�?�������f1f469cd-dd83-48f1-b3db-d248a827f369��L_M~� +?00:00:00-00:00:00A�I���x�@�?�������68a7e0a4-0eea-4b1d-b0e8-7e586ea2d248>�C�M~� ?00:00:00-00:00:00A�I���wI��@�?�������e01bae53-0577-4143-bfb4-655afb353a35���M~� ?00:00:00-00:00:00A�I���v���@�?�������6b06afa2-327d-4d63-8958-7655b4e95cd0?wQ�M~� ?00:00:00-00:00:00A�I���u���@�?�������7af01cca-8756-4b93-a919-d35f8c8c383d?�IM~�?00:00:00-00:00:00A�I���u'��@�?�������db52be3b-1daa-4a3c-b1ca-cdf4762e42a8��djM~�?00:00:00-00:00:00A�I���tq��A?�������8b853852-f668-430b-9a25-69e5117389e0� �mM~�?00:00:00-00:00:00A�I���s���A?�������667ff57b-1f08-4032-8518-0fd20e30fa11��ՑM~�?00:00:00-00:00:00A�I���s��A?�������950eb7eb-2f22-49e9-b859-c46b11684d89��c0M~�?00:00:00-00:00:00A�I���rO��A?�������60b27ba8-7708-431c-ba8f-9e1dfa99c800>BCM~�?00:00:00-00:00:00A�I���q���A ?�������0edf6eff-84db-4166-97c3-d0ed04dbd889?&U�M~�?00:00:00-00:00:00A�-���x�A(?�������6182968a-469d-4ee7-b65a-bc645e8c1a0f��%`M~�?00:00:00-00:00:00A�-���wI��A0?�������57520487-510e-4295-8e6c-f8ea5d214f78=b� M~�?00:00:00-00:00:00A�-���v���A8?�������79ddf43f-18f2-4efa-9ebe-9800f161d546�V� +M~�?00:00:00-00:00:00A�-���u���A@?�������34495849-1cce-4086-b86b-815d38455285�P�nM~�?00:00:00-00:00:00A�-���u'��AH?�������58672cfd-685a-4f63-9c75-73e45634f920�c25M~�?00:00:00-00:00:00A�-���tq��AP?�������f89c9984-6e07-4e4b-a804-948457450d80��=M~�?00:00:00-00:00:00A�-���s���AX?�������36059ee1-aca1-4547-8ca1-c581e1dd51be�7�M~�?00:00:00-00:00:00A�-���s��A`?�������14e00510-b4df-41aa-ae7e-1c3f0ea90238��}M~�?00:00:00-00:00:00A�-���rO��Ah?�������9b5eb77f-c7f0-4791-a1ca-9ab9696202d2�3M~�?00:00:00-00:00:00A�-���q���Ap?�������5c600cd6-5138-4600-b78f-9a242bdfffe0�X�M~�?00:00:00-00:00:00A���x�Ax?�������89c22dbb-6493-46ab-b3d2-8e2690276581�xSM~�?00:00:00-00:00:00A���wI��A�?�������a4ae120a-f7db-4fb1-b4ab-b2d0ed005de9���UM~� ?00:00:00-00:00:00A���v���A�?�������fab9f9bb-7371-4e43-a14e-2cf0cceac294����fM~�%?00:00:00-00:00:00A���s��A�?�������5ffdcbf5-0392-466b-96de-91e216563d17��UM~�&?00:00:00-00:00:00A���rO��A�?�������f0e2a560-f1c1-4b7a-b321-50a2f8ed2e4f���\M~�'?00:00:00-00:00:00A���q���A�?�������e8603036-0a53-4500-b14d-c0fa40e8485c�� M~�(?00:00:00-00:00:00A�����x�A�?�������5f0a3f2d-7c61-4a09-8e6e-ee99225d9cad>{�}M~�)?00:00:00-00:00:00A�����wI��A�?�������f89946a0-36f6-4387-a842-a6a5a3edb30f?x<�M~�*?00:00:00-00:00:00A�����v���A�?�������f3487366-bc26-420a-b244-9d0367650a2b>x�M~�+?00:00:00-00:00:00A�����u���A�?�������16378236-85b9-4977-b072-935e7ef39e6d����M~�,?00:00:00-00:00:00A�����u'��A�?�������02e8f848-2489-4078-a78f-1ead1ae4ea29>1�M~�-?00:00:00-00:00:00A�����tq��A�?�������3ef89b8b-2e29-4a5c-818b-94a424a0c809>�-^M~�.?00:00:00-00:00:00A�����s���A�?�������28430f68-2c41-49b1-b04a-d7a35e5d3b7d?fO�M~�/?00:00:00-00:00:00A�����s��A�?�������52133824-2182-4aec-94fc-76b69dfb5968>�"@M~�0?00:00:00-00:00:00A�����rO��A�?�������21f7635f-e7ce-4824-8b4f-9854c96b372c��P M~�1?00:00:00-00:00:00A�����q���A�?�������367fac53-01ce-48aa-b95b-84656811fec0�)�]M~�2?00:00:00-00:00:00A��.��x�A�?�������73230a93-8365-4424-affd-1c3427421a2a?Z��M~�3?00:00:00-00:00:00A��.��wI��A�?�������cebb773c-0cd9-46d1-af50-cccba9bf87e1���M~�4?00:00:00-00:00:00A��.��v���A�?�������9651c0ed-e85b-453b-9ba2-c1ec2e2ea698>Y�,M~�5?00:00:00-00:00:00A��.��u���A�?�������aad5962e-e859-4fd5-9434-57ff77d35a7a?9)RM~�6?00:00:00-00:00:00A��.��u'��A�?�������9bfc96e9-a4f3-482c-8dc8-0586d575df41�0M~�7?00:00:00-00:00:00A��.��tq��A�?�������caf02810-f06e-42ec-89a8-a4a5009b451a���M~�8?00:00:00-00:00:00A��.��s���A�?�������145ba62b-f743-478c-848b-4a53ef64c315��.�M~�9?00:00:00-00:00:00A��.��s��A�?�������bf50671e-fbd7-43dc-ab36-b647b6fc46e1��uM~�:?00:00:00-00:00:00A��.��rO��A�?�������ea33081f-49e8-4098-aa72-b4ecac716a11��Y�M~�;?00:00:00-00:00:00A��.��q���A�?�������d79633be-9ba2-4b48-88d6-0fd636ee7de5����M~�<?00:00:00-00:00:00A���x�A�?�������ad436a70-730e-49b1-8353-b3caeda973c9�Z�#M~�=?00:00:00-00:00:00A���wI��A�?�������66c233bc-2511-4ebd-8150-38611a5bda13���M~�>?00:00:00-00:00:00A���v���A�?�������8b6a8eff-d95e-4892-903d-7bf620531b7f�,�5M~�??00:00:00-00:00:00A���u���B?�������f711c957-d04a-4c5e-8fab-9e96df4358e3�g�UM~�@?00:00:00-00:00:00A���u'��B?�������ec83c28a-8dbb-4829-b1da-401850751c41��c�M~�A?00:00:00-00:00:00A���tq��B?�������744c0576-3576-4558-9f9d-784229fd941d���M~�B?00:00:00-00:00:00A���s���B?�������64dddf71-5230-4357-9a44-90d4abcba28f>�mM~�C?00:00:00-00:00:00A���s��B?�������96e076d6-79a1-45c8-b200-a292154c3764?B�'M~�D?00:00:00-00:00:00A���rO��B +?�������07b0535d-52b5-4ade-8065-1269ba9cfed8��>�M~�E?00:00:00-00:00:00A���q���B ?�������4fed5d2c-0790-46cb-8694-49fa23d39610��{[M~�F?00:00:00-00:00:00A��J��x�B?�������35c72544-5b11-4d9b-aced-fc3c0e2edd18���yM~�G?00:00:00-00:00:00A��J��wI��B?�������7da4e033-e46a-412a-8170-b8d0cda23a2c>��M~�H?00:00:00-00:00:00A��J��v���B?�������3e533902-00b7-44d6-a4bd-a3e40e7945a4�Ϥ�M~�I?00:00:00-00:00:00A��J��u���B?�������4454e3cd-2100-4893-afd3-abf81ccd9342>�M~�J?00:00:00-00:00:00A��J��u'��B?�������c4690174-5dc7-4b3e-95ee-e935c37f0012���/M~�K?00:00:00-00:00:00A��J��tq��B?�������19fc1cdf-2d75-4de0-9fe2-e38d8c2af567�T 0M~�L?00:00:00-00:00:00A��J��s���B?�������3369e286-1a63-4658-9588-2aa9e7cd4045? ��M~�M?00:00:00-00:00:00A��J��s��B?�������b0535e3c-fa34-4b8e-9fef-8c49ef5d5f78��IM~�N?00:00:00-00:00:00A��J��rO��B?�������b339ca46-1583-4609-b0b0-a579f610426c>)zM~�O?00:00:00-00:00:00A��J��q���B ?�������555e05d0-3bf7-4c05-af92-4d83c471450a?S�oM~�P?00:00:00-00:00:00A�����x�B"?�������661bdaa9-4359-42c0-ada4-612dcc8477ca���fM~�Q?00:00:00-00:00:00A�����wI��B$?�������1bb6ba57-f934-405e-86ad-4718cc3c83a2��dOM~�R?00:00:00-00:00:00A�����v���B&?�������fb125f3a-021c-4c05-a785-c16f9c382aeb���)M~�S?00:00:00-00:00:00A�����u���B(?�������eb86c23b-ccc1-4982-9dde-973e854a8da8>��M~�T?00:00:00-00:00:00A�����u'��B*?�������c5cc8007-5b91-4ed4-9e95-44417e279a78?�vM~�U?00:00:00-00:00:00A�����tq��B,?�������bff9fbeb-79ca-46f6-8d43-49378bf214a4���:M~�V?00:00:00-00:00:00A�����s���B.?�������481fbbe3-693f-4f4c-8b94-7ae4a0e2d028����M~�W?00:00:00-00:00:00A�����s��B0?�������dfa8a200-82f4-47c5-a374-cc8da78256fb�C��M~�X?00:00:00-00:00:00A�����rO��B2?�������ad6fde81-a22f-4b70-99a1-7511e041aea0��(M~�Y?00:00:00-00:00:00A�����q���B4?�������26baedaf-c325-486a-bd5f-b61cf9552d9f���5M~�Z?00:00:00-00:00:00A�ff��x�B6?�������f057383d-5d97-4f34-877a-6e4da1cf9428��� M~�[?00:00:00-00:00:00A�ff��wI��B8?�������5a046e38-7d36-44ee-a8ef-9deb1f4141b0?r��M~�\?00:00:00-00:00:00A�ff��v���B:?�������f8561e86-a5f6-414c-bccf-e14037b72a84>�,sM~�]?00:00:00-00:00:00A�ff��u���B<?�������3a1cc473-e29a-419c-97eb-57916e32b985?)�3M~�^?00:00:00-00:00:00A�ff��u'��B>?�������db01f5a3-7a64-4811-a187-233375343e20>�0M~�_?00:00:00-00:00:00A�ff��tq��B@?�������320da9a0-a552-4127-8331-4e5647638adc��<�M~�`?00:00:00-00:00:00A�ff��s���BB?�������8dcea419-9aeb-434b-a236-adfdd17e6d4b���)M~�a?00:00:00-00:00:00A�ff��s��BD?�������472372e2-bc75-41d1-a68a-c4b08e1350ef=���M~�b?00:00:00-00:00:00A�ff��rO��BF?�������e3a439fd-9419-42a9-abac-5606eb03cc02���M~�c?00:00:00-00:00:00A�ff��q���BH?�������911583a9-7ad0-4fce-9854-c1e9116fe956=#wnM~�d��00:00:00-00:00:00A� +=��p���@@���Z�Z����f4dcfd49-0bc9-48f5-9325-c46ee4c86f44>��M~�e��00:00:00-00:00:00A�=q��q���B ��������e15fcaa7-4aca-4004-bec3-9e668366f2e4=�`M~� \ No newline at end of file diff --git a/tests/test_files/synthetic_cube.fits b/tests/test_files/synthetic_cube.fits new file mode 100644 index 00000000..7edbc7b9 Binary files /dev/null and b/tests/test_files/synthetic_cube.fits differ diff --git a/tests/test_files/synthetic_cube_bkg.fits b/tests/test_files/synthetic_cube_bkg.fits new file mode 100644 index 00000000..0b99475c Binary files /dev/null and b/tests/test_files/synthetic_cube_bkg.fits differ diff --git a/tests/test_files/synthetic_cube_psf.fits b/tests/test_files/synthetic_cube_psf.fits new file mode 100644 index 00000000..f265a89c Binary files /dev/null and b/tests/test_files/synthetic_cube_psf.fits differ diff --git a/tests/test_files/synthetic_cube_rms.fits b/tests/test_files/synthetic_cube_rms.fits new file mode 100644 index 00000000..440d96f7 Binary files /dev/null and b/tests/test_files/synthetic_cube_rms.fits differ diff --git a/tests/test_files/synthetic_image.fits b/tests/test_files/synthetic_image.fits new file mode 100644 index 00000000..eb18455f Binary files /dev/null and b/tests/test_files/synthetic_image.fits differ diff --git a/tests/test_files/synthetic_image_psf.fits b/tests/test_files/synthetic_image_psf.fits new file mode 100644 index 00000000..78d1fc15 Binary files /dev/null and b/tests/test_files/synthetic_image_psf.fits differ diff --git a/tests/unit/make_test_data.py b/tests/unit/make_test_data.py index 38b3c07c..3fad77be 100644 --- a/tests/unit/make_test_data.py +++ b/tests/unit/make_test_data.py @@ -10,18 +10,20 @@ from AegeanTools.logging import logger, logging from AegeanTools.source_finder import FWHM2CC from AegeanTools.wcs_helpers import WCSHelper +from AegeanTools.catalogs import write_catalog -__author__ = 'Paul Hancock' +__author__ = "Paul Hancock" -imsize = (256, 512) # non square picks up more errors +imsize = (8, 256, 512) # non square picks up more errors noise = 0.5 # non unity noise picks up more errors psf = [30, 30, 0] # psf in arcsec,arcsec,degrees pix_per_beam = 3.5 seed = 123987554 ra, dec = (30, -15) +nu0 = 162 * 1e6 -def make_noise_image(): +def make_noise_image(as_cube=False): """ Create a noise-only image and a WCS object for use in testing. @@ -33,22 +35,36 @@ def make_noise_image(): wcs : :class:`astropy.wcs.WCS` A wcs object that is valid for this image. """ + np.random.seed(seed) image = np.random.random(size=imsize) - # force zero mean and unit variance - image -= np.mean(image) - image /= np.std(image) - # make rms = noise - image *= noise - image = gaussian_filter(image, sigma=pix_per_beam*FWHM2CC) + + for i in range(image.shape[0]): + image[i] = gaussian_filter(image[i], sigma=pix_per_beam * FWHM2CC) + # force zero mean and unit variance + image[i] -= np.mean(image[i]) + image[i] /= np.std(image[i]) + image[i] *= noise + image = np.array(image, dtype=np.float32) - wcs = WCS(naxis=2) - wcs.wcs.crpix = [imsize[0]/2, imsize[1]/2] - wcs.wcs.cdelt = np.array( - [psf[0]/3600/pix_per_beam, psf[1]/3600/pix_per_beam]) - wcs.wcs.crval = [ra, dec] - wcs.wcs.ctype = ["RA---SIN", "DEC--SIN"] + if as_cube: + wcs = WCS(naxis=3) + wcs.wcs.crpix = [imsize[1] / 2, imsize[2] / 2, 0] + wcs.wcs.cdelt = np.array( + [psf[0] / 3600 / pix_per_beam, psf[1] / 3600 / pix_per_beam, 1e6] + ) + wcs.wcs.crval = [ra, dec, nu0] + wcs.wcs.ctype = ["RA---SIN", "DEC--SIN", "FREQ"] + else: + wcs = WCS(naxis=2) + wcs.wcs.crpix = [imsize[1] / 2, imsize[2] / 2] + wcs.wcs.cdelt = np.array( + [psf[0] / 3600 / pix_per_beam, psf[1] / 3600 / pix_per_beam] + ) + wcs.wcs.crval = [ra, dec] + wcs.wcs.ctype = ["RA---SIN", "DEC--SIN"] + image = image[0] return image, wcs @@ -57,8 +73,8 @@ def make_psf_map(base): Create a small psf map for use in testing. """ - psf_data = np.ones((3, imsize[0], imsize[1])) - psf_data *= np.array(psf)[:, None, None]/3600. # acrsec -> degrees + psf_data = np.ones((3, imsize[1], imsize[2])) + psf_data *= np.array(psf)[:, None, None] / 3600.0 # acrsec -> degrees # create a few 'holes' in the map psf_data[:, 56, 34] = np.nan @@ -66,7 +82,7 @@ def make_psf_map(base): psf_data[:, :, 0] = np.nan hdu = deepcopy(base) hdu[0].data = psf_data - hdu[0].header['CTYPE3'] = ('Beam', '0=a,1=b,2=pa (degrees)') + hdu[0].header["CTYPE3"] = ("Beam", "0=a,1=b,2=pa (degrees)") return hdu @@ -76,42 +92,50 @@ def make_catalogue(): """ cat = [] # creage a grid of point sources with SNR from 1-100 - ras = np.linspace(ra-0.2, ra+0.8, 10) - decs = np.linspace(dec-0.5, dec-0.1, 10) - print(ras) - print(decs) - fluxes = [(a+1)*noise for a in range(100)] + ras = np.linspace(ra - 0.2, ra + 0.8, 10) + decs = np.linspace(dec - 0.5, dec - 0.1, 10) + fluxes = [(a + 1) * noise for a in range(100)] + for i, r in enumerate(ras): for j, d in enumerate(decs): - src = models.SimpleSource() + src = models.ComponentSource() src.ra = r src.dec = d + src.ra_str = "00:00:00" + src.dec_str = "-00:00:00" src.a, src.b, src.pa = psf - src.island = i*10 + j - src.peak_flux = fluxes[i*10 + j] + src.psf_a, src.psf_b, src.psf_pa = psf + src.island = i * 10 + j + src.peak_flux = fluxes[i * 10 + j] src.background = 0 src.local_rms = noise src.err_peak_flux = noise cat.append(src) # create an island from a diffuse source - src = models.SimpleSource() + src = models.ComponentSource() + src.island = 100 src.ra = 30.88 src.dec = -15.04 + src.ra_str = "00:00:00" + src.dec_str = "-00:00:00" src.a, src.b, src.pa = psf + src.psf_a, src.psf_b, src.psf_pa = psf src.a *= 3 src.b *= 3 - src.island = 101 - src.peak_flux = 6*noise + src.peak_flux = 6 * noise cat.append(src) # ensure that one of the islands has two sources in it - src = models.SimpleSource() + src = models.ComponentSource() src.ra = 30.78 src.dec = -15.102 src.a, src.b, src.pa = psf - src.island = 102 - src.peak_flux = 80*noise + src.psf_a, src.psf_b, src.psf_pa = psf + src.ra_str = "00:00:00" + src.dec_str = "-00:00:00" + src.island = 101 + src.peak_flux = 80 * noise cat.append(src) return cat @@ -119,22 +143,49 @@ def make_catalogue(): if __name__ == "__main__": # configure logging - logging_level = logging.DEBUG + logging_level = logging.INFO logger.setLevel(logging_level) - image, wcs = make_noise_image() + # make 3d noise image + image, wcs = make_noise_image(as_cube=True) + logger.info(f"Image shape {image.shape}") header = wcs.to_header() - header['BMAJ'] = psf[0]/3600 - header['BMIN'] = psf[1]/3600 - header['BPA'] = psf[2] + header["BMAJ"] = psf[0] / 3600 + header["BMIN"] = psf[1] / 3600 + header["BPA"] = psf[2] hdu = fits.HDUList(hdus=[fits.PrimaryHDU(header=header, data=image)]) - header = hdu[0].header - cat = make_catalogue() - wcshelper = WCSHelper.from_header(header) - hdu[0].data += AeRes.make_model(cat, image.shape, wcshelper) + # make catalogue without alpha + cat_no_alpha = make_catalogue() + + # make catalogue with alpha + # cat_with_alpha = deepcopy(cat_no_alpha) + cat_with_alpha = [ + models.ComponentSource3D.from_component_source(c) for c in cat_no_alpha + ] + for c in cat_with_alpha: + c.alpha = np.random.uniform(-2, 1) + c.nu0 = nu0 + + # save the image cube, catalogue, and psf map + wcshelper = WCSHelper.from_header(hdu[0].header) + hdu[0].data += AeRes.make_model(cat_with_alpha, image.shape, wcshelper) + hdu.writeto("synthetic_cube.fits", overwrite=True) + psf_hdu = make_psf_map(hdu) + psf_hdu.writeto("synthetic_cube_psf.fits", overwrite=True) + write_catalog("synthetic_cat_with_alpha.fits", cat_with_alpha, fmt="fits") - hdu.writeto('synthetic_test.fits', overwrite=True) + # make 2d image and save + image, wcs = make_noise_image(as_cube=False) + header = wcs.to_header() + header["BMAJ"] = psf[0] / 3600 + header["BMIN"] = psf[1] / 3600 + header["BPA"] = psf[2] + hdu = fits.HDUList(hdus=[fits.PrimaryHDU(header=header, data=image)]) + wcshelper = WCSHelper.from_header(hdu[0].header) + hdu[0].data += AeRes.make_model(cat_no_alpha, (1, *image.shape), wcshelper)[0, :, :] + hdu.writeto("synthetic_image.fits", overwrite=True) psf_hdu = make_psf_map(hdu) - psf_hdu.writeto('synthetic_test_psf.fits', overwrite=True) + psf_hdu.writeto("synthetic_image_psf.fits", overwrite=True) + write_catalog("synthetic_cat_no_alpha.fits", cat_no_alpha, fmt="fits") diff --git a/tests/unit/test_AeRes.py b/tests/unit/test_AeRes.py index 362d7346..d92592a1 100755 --- a/tests/unit/test_AeRes.py +++ b/tests/unit/test_AeRes.py @@ -9,27 +9,38 @@ from AegeanTools import catalogs, wcs_helpers from astropy.io import fits -__author__ = 'Paul Hancock' +__author__ = "Paul Hancock" def test_load_sources(): """Test load_sources""" - filename = 'tests/test_files/1904_comp.fits' + filename = "tests/test_files/1904_comp.fits" cat = ar.load_sources(filename) if cat is None: raise AssertionError("load_sources failed") return +def test_load_soruces_with_alpha(): + """Test load_sources with alpha column""" + filename = "tests/test_files/synthetic_cat_with_alpha_comp.fits" + cat = ar.load_sources(filename) + if cat is None: + raise AssertionError("load_sources failed with alpha column") + return + + def test_load_soruces_renamed_columns(): """Test load_sources with renamed columns""" - filename = 'tests/test_files/1904_comp_renamed_cols.fits' - colnames = {'ra_col': 'RAJ2000', - 'dec_col': 'DEJ2000', - 'peak_col': 'S', - 'a_col': 'bmaj', - 'b_col': 'bmin', - 'pa_col': 'bpa'} + filename = "tests/test_files/1904_comp_renamed_cols.fits" + colnames = { + "ra_col": "RAJ2000", + "dec_col": "DEJ2000", + "peak_col": "S", + "a_col": "bmaj", + "b_col": "bmin", + "pa_col": "bpa", + } cat = ar.load_sources(filename, **colnames) if cat is None: raise AssertionError("load_sources failed with renamed columns") @@ -37,13 +48,13 @@ def test_load_soruces_renamed_columns(): def test_load_sources_missing_columns(): - filename = 'tests/test_files/1904_comp.fits' + filename = "tests/test_files/1904_comp.fits" table = catalogs.load_table(filename) - table.rename_column('ra', 'RAJ2000') - table.write('dlme.fits') - cat = ar.load_sources('dlme.fits') - if os.path.exists('dlme.fits'): - os.remove('dlme.fits') + table.rename_column("ra", "RAJ2000") + table.write("dlme.fits") + cat = ar.load_sources("dlme.fits") + if os.path.exists("dlme.fits"): + os.remove("dlme.fits") if cat is not None: raise AssertionError("Missing columns should be caught, but weren't") @@ -52,72 +63,103 @@ def test_load_sources_missing_columns(): def test_make_model(): """Test make_model""" - filename = 'tests/test_files/1904_comp.fits' + filename = "tests/test_files/1904_comp.fits" + sources = ar.load_sources(filename) + hdulist = fits.open("tests/test_files/1904-66_SIN.fits") + wcs_helper = wcs_helpers.WCSHelper.from_header(header=hdulist[0].header) + # regular run + model = ar.make_model( + sources=sources, shape=(1, *hdulist[0].data.shape), wcshelper=wcs_helper + ) + if np.all(model == 0.0): + raise AssertionError("Model is empty") + + # model with *all* sources outside region + # shape (100,2) means we only sometimes reject a source based on it's x-coord + model = ar.make_model(sources=sources, shape=(1, 100, 2), wcshelper=wcs_helper) + if not np.all(model == 0.0): + raise AssertionError("Model is *not* empty") + + +def test_make_model_with_alpha(): + """Test make_model with alpha column""" + filename = "tests/test_files/synthetic_cat_with_alpha_comp.fits" sources = ar.load_sources(filename) - hdulist = fits.open('tests/test_files/1904-66_SIN.fits') + hdulist = fits.open("tests/test_files/synthetic_cube.fits") wcs_helper = wcs_helpers.WCSHelper.from_header(header=hdulist[0].header) # regular run - model = ar.make_model(sources=sources, shape=hdulist[0].data.shape, - wcshelper=wcs_helper) - if np.all(model == 0.): + model = ar.make_model( + sources=sources, shape=hdulist[0].data.shape, wcshelper=wcs_helper + ) + if np.all(model == 0.0): raise AssertionError("Model is empty") # model with *all* sources outside region # shape (100,2) means we only sometimes reject a source based on it's x-coord - model = ar.make_model(sources=sources, shape=(100, 2), - wcshelper=wcs_helper) - if not np.all(model == 0.): + model = ar.make_model(sources=sources, shape=(1, 100, 2), wcshelper=wcs_helper) + if not np.all(model == 0.0): raise AssertionError("Model is *not* empty") def test_make_masked_model(): """Test make_model when a mask is being used""" - filename = 'tests/test_files/1904_comp.fits' + filename = "tests/test_files/1904_comp.fits" sources = ar.load_sources(filename) - hdulist = fits.open('tests/test_files/1904-66_SIN.fits') + hdulist = fits.open("tests/test_files/1904-66_SIN.fits") wcs_helper = wcs_helpers.WCSHelper.from_header(header=hdulist[0].header) # test mask with sigma - model = ar.make_model(sources=sources, shape=hdulist[0].data.shape, - wcshelper=wcs_helper, mask=True) + model = ar.make_model( + sources=sources, + shape=(1, *hdulist[0].data.shape), + wcshelper=wcs_helper, + mask=True, + ) finite = np.where(np.isfinite(model)) - if np.all(model == 0.): + if np.all(model == 0.0): raise AssertionError("Model is empty") if not np.any(np.isnan(model)): raise AssertionError("Model is not masked") - if not np.all(model[finite] == 0.): + if not np.all(model[finite] == 0.0): raise AssertionError("Model has values that are not zero or nan") # test mask with frac - model = ar.make_model(sources=sources, shape=hdulist[0].data.shape, - wcshelper=wcs_helper, mask=True, frac=0.1) + model = ar.make_model( + sources=sources, + shape=(1, *hdulist[0].data.shape), + wcshelper=wcs_helper, + mask=True, + frac=0.1, + ) finite = np.where(np.isfinite(model)) - if np.all(model == 0.): + if np.all(model == 0.0): raise AssertionError("Model is empty") if not np.any(np.isnan(model)): raise AssertionError("Model is not masked") - if not np.all(model[finite] == 0.): + if not np.all(model[finite] == 0.0): raise AssertionError("Model has values that are not zero or nan") def test_make_residual(): """Test make_residual""" - fitsfile = 'tests/test_files/1904-66_SIN.fits' - catalog = 'tests/test_files/1904_comp.fits' - residual = 'tests/temp/residual.fits' - masked = 'tests/temp/masked.fits' + fitsfile = "tests/test_files/1904-66_SIN.fits" + catalog = "tests/test_files/1904_comp.fits" + residual = "tests/temp/residual.fits" + masked = "tests/temp/masked.fits" # default operation - ar.make_residual(fitsfile=fitsfile, catalog=catalog, - rfile=residual, mfile=masked, mask=False) + ar.make_residual( + fitsfile=fitsfile, catalog=catalog, rfile=residual, mfile=masked, mask=False + ) if not os.path.exists(masked): raise AssertionError("Mask file not written") os.remove(masked) # with masking - ar.make_residual(fitsfile=fitsfile, catalog=catalog, - rfile=residual, mfile=masked, mask=True) + ar.make_residual( + fitsfile=fitsfile, catalog=catalog, rfile=residual, mfile=masked, mask=True + ) if not os.path.exists(masked): raise AssertionError("Mask file not written") os.remove(masked) @@ -126,6 +168,6 @@ def test_make_residual(): if __name__ == "__main__": # introspect and run all the functions starting with 'test' for f in dir(): - if f.startswith('test'): + if f.startswith("test"): print(f) globals()[f]() diff --git a/tests/unit/test_BANE.py b/tests/unit/test_BANE.py index e9ef3422..8a474ff4 100755 --- a/tests/unit/test_BANE.py +++ b/tests/unit/test_BANE.py @@ -48,7 +48,7 @@ def test_filter_image(): raise AssertionError() os.remove(bkg) - BANE.filter_image(fname, cores=2, out_base=outbase, twopass=True, compressed=True) + BANE.filter_image(fname, cores=2, out_base=outbase, compressed=True) if not os.path.exists(rms): raise AssertionError() @@ -150,6 +150,20 @@ def test_BSCALE(): return +def test_cube_as_cube(): + """ + Ensure that running BANE on a cube delivers a cube output + """ + fname = "tests/test_files/1904-66_SIN_3d.fits" + # don't crash and die + try: + for index in [0, 1, 2]: + BANE.filter_image(fname, out_base=None, nslice=1, cube_index=None) + except Exception as e: + raise AssertionError("BANE failed to work on 3d image") + return + + if __name__ == "__main__": # introspect and run all the functions starting with 'test' for f in dir(): diff --git a/tests/unit/test_addSuffix.py b/tests/unit/test_addSuffix.py new file mode 100644 index 00000000..1b081bac --- /dev/null +++ b/tests/unit/test_addSuffix.py @@ -0,0 +1,29 @@ +#! /usr/bin/env python + +import unittest +from AegeanTools.exceptions import AegeanSuffixError +from AegeanTools.CLI.aegean import addSuffix +__author__ = "Ahmed Abdelmuniem Abdalla Mohammed" + + +class TestAddSuffix(unittest.TestCase): + + def test_add_suffix_int(self): + self.assertEqual(addSuffix("sample.txt", 1), "sample_01.txt") + self.assertEqual(addSuffix("sample_2.txt", 10), "sample_2_10.txt") + self.assertEqual(addSuffix("sample_3.txt", 100), "sample_3_100.txt") + + def test_add_suffix_str(self): + self.assertEqual(addSuffix("sample.txt", "_comp"), "sample_comp.txt") + self.assertEqual(addSuffix("sample_2.txt", "__comp"), "sample_2__comp.txt") + + def test_add_suffix_invalid_type(self): + with self.assertRaises(AegeanSuffixError): + addSuffix("sample.txt", ["invalid"]) + +if __name__ == "__main__": + # introspect and run all the functions starting with 'test' + for f in dir(): + if f.startswith("test"): + print(f) + globals()[f]() \ No newline at end of file diff --git a/tests/unit/test_fits_tools.py b/tests/unit/test_fits_tools.py index d1f61e0f..f36dba9c 100755 --- a/tests/unit/test_fits_tools.py +++ b/tests/unit/test_fits_tools.py @@ -166,6 +166,76 @@ def test_load_image_band_multi_bands(): return +def test_load_image_band_data_header_check(): + """Load an image to test if it is of type np.ndarray object + and that the header is a fits.hdu.Header object""" + try: + data, header = fits_tools.load_image_band( + "tests/test_files/synthetic_cube.fits", as_cube=True + ) + except Exception as e: + raise e + if not isinstance(data, np.ndarray): + raise AssertionError("Loaded data is not an np.ndarray object") + if not isinstance(header, fits.Header): + raise AssertionError("header is not a fits.hdu.Header object") + return + + +def test_load_image_band_cube_as_cube_true(): + """Load an image of a cube with as_cube = True""" + try: + data, header = fits_tools.load_image_band( + "tests/test_files/synthetic_cube.fits", as_cube=True + ) + except Exception as e: + raise e + if len(data.shape) < 3: + raise AssertionError("Loaded data is not 3-Dimensional") + return + + +def test_load_image_band_cube_as_cube_false(): + """Load an image of a cube with as_cube = False""" + try: + data, header = fits_tools.load_image_band( + "tests/test_files/synthetic_cube.fits", as_cube=False + ) + except Exception as e: + raise e + + if len(data.shape) < 2: + raise AssertionError("Loaded data is not 2-Dimensional") + return + + +def no_test_load_image_band_2d_as_cube_true(): + """Load an image of a band with as_cube = True""" + try: + data, header = fits_tools.load_image_band( + "tests/test_files/1904-66_AIT.fits", as_cube=True + ) + except AegeanError as e: + return + else: + raise AssertionError("Data passed as a cube but only 2 axes were provided") + + +def test_load_image_band_2d_as_cube_false(): + """Load an image using default values""" + try: + data, header = fits_tools.load_image_band( + "tests/test_files/1904-66_AIT.fits", as_cube=False + ) + except Exception as e: + raise e + if not isinstance(data, np.ndarray): + raise AssertionError("Loaded data is not an np.ndarray object") + if not isinstance(header, fits.Header): + raise AssertionError("header is not a fits.hdu.Header object") + return + + def test_load_image_band_cube_index(): return diff --git a/tests/unit/test_fitting.py b/tests/unit/test_fitting.py index c8453393..7e9f1601 100755 --- a/tests/unit/test_fitting.py +++ b/tests/unit/test_fitting.py @@ -6,31 +6,31 @@ import numpy as np from AegeanTools import fitting -__author__ = 'Paul Hancock' +__author__ = "Paul Hancock" def make_model(): """Test that we can make lmfit.Parameter models""" model = lmfit.Parameters() - model.add('c0_amp', 1, vary=True) - model.add('c0_xo', 5, vary=True) - model.add('c0_yo', 5, vary=True) - model.add('c0_sx', 2.001, vary=False) - model.add('c0_sy', 2, vary=False) - model.add('c0_theta', 0, vary=False) - model.add('components', 1, vary=False) + model.add("c0_amp", 1, vary=True) + model.add("c0_xo", 5, vary=True) + model.add("c0_yo", 5, vary=True) + model.add("c0_sx", 2.001, vary=False) + model.add("c0_sy", 2, vary=False) + model.add("c0_theta", 0, vary=False) + model.add("components", 1, vary=False) return model def test_elliptical_gaussian(): """Test our elliptical gaussian creation function""" x, y = np.indices((3, 3)) - gauss = fitting.elliptical_gaussian( - x, y, amp=1, xo=0, yo=1, sx=1, sy=1, theta=0) + gauss = fitting.elliptical_gaussian(x, y, amp=1, xo=0, yo=1, sx=1, sy=1, theta=0) if np.any(np.isnan(gauss)): raise AssertionError() gauss = fitting.elliptical_gaussian( - x, y, amp=1, xo=0, yo=1, sx=1, sy=1, theta=np.inf) + x, y, amp=1, xo=0, yo=1, sx=1, sy=1, theta=np.inf + ) if not (np.all(np.isnan(gauss))): raise AssertionError() @@ -46,30 +46,6 @@ def test_CandBmatrix(): raise AssertionError() -def test_hessian_shape(): - """Test that the hessian has the correct shape""" - # test a single component model - model = make_model() - nvar = 3 - x, y = np.indices((10, 10)) - Hij = fitting.hessian(model, x, y) - if not (Hij.shape == (nvar, nvar, 10, 10)): - raise AssertionError() - - # now add another component - model.add('c1_amp', 1, vary=True) - model.add('c1_xo', 5, vary=True) - model.add('c1_yo', 5, vary=True) - model.add('c1_sx', 2.001, vary=True) - model.add('c1_sy', 2, vary=True) - model.add('c1_theta', 0, vary=True) - nvar = 9 - model['components'].value = 2 - Hij = fitting.hessian(model, x, y) - if not (Hij.shape == (nvar, nvar, 10, 10)): - raise AssertionError() - - def test_jacobian_shape(): """ Test to see if the Jacobian matrix if of the right shape @@ -82,14 +58,14 @@ def test_jacobian_shape(): if not (Jij.shape == (nvar, 10, 10)): raise AssertionError() - model.add('c1_amp', 1, vary=True) - model.add('c1_xo', 5, vary=True) - model.add('c1_yo', 5, vary=True) - model.add('c1_sx', 2.001, vary=True) - model.add('c1_sy', 2, vary=True) - model.add('c1_theta', 0, vary=True) + model.add("c1_amp", 1, vary=True) + model.add("c1_xo", 5, vary=True) + model.add("c1_yo", 5, vary=True) + model.add("c1_sx", 2.001, vary=True) + model.add("c1_sy", 2, vary=True) + model.add("c1_theta", 0, vary=True) nvar = 9 - model['components'].value = 2 + model["components"].value = 2 Jij = fitting.jacobian(model, x, y) if not (Jij.shape == (nvar, 10, 10)): raise AssertionError() @@ -106,14 +82,14 @@ def test_emp_vs_ana_jacobian(): if not (np.max(diff) < 1e-5): raise AssertionError() - model.add('c1_amp', 1, vary=True) - model.add('c1_xo', 5, vary=True) - model.add('c1_yo', 5, vary=True) - model.add('c1_sx', 2.001, vary=True) - model.add('c1_sy', 2, vary=True) - model.add('c1_theta', 0, vary=True) + model.add("c1_amp", 1, vary=True) + model.add("c1_xo", 5, vary=True) + model.add("c1_yo", 5, vary=True) + model.add("c1_sx", 2.001, vary=True) + model.add("c1_sy", 2, vary=True) + model.add("c1_theta", 0, vary=True) - model['components'].value = 2 + model["components"].value = 2 emp_Jij = fitting.emp_jacobian(model, x, y) ana_Jij = fitting.jacobian(model, x, y) diff = np.abs(ana_Jij - emp_Jij) @@ -121,63 +97,9 @@ def test_emp_vs_ana_jacobian(): raise AssertionError() -def test_emp_vs_ana_hessian(): - """Test that the empirical and analytical Hessians agree""" - model = make_model() - - x, y = np.indices((10, 10)) - emp_Hij = fitting.emp_hessian(model, x, y) - ana_Hij = fitting.hessian(model, x, y) - diff = np.abs(ana_Hij - emp_Hij) - if not (np.max(diff) < 1e-5): - raise AssertionError() - - model.add('c1_amp', 1, vary=True) - model.add('c1_xo', 5, vary=True) - model.add('c1_yo', 5, vary=True) - model.add('c1_sx', 2.001, vary=True) - model.add('c1_sy', 2, vary=True) - model.add('c1_theta', 0, vary=True) - - model['components'].value = 2 - emp_Hij = fitting.emp_hessian(model, x, y) - ana_Hij = fitting.hessian(model, x, y) - diff = np.abs(ana_Hij - emp_Hij) - if not (np.max(diff) < 1): - raise AssertionError() - - -def test_make_ita(): - """Test make_ita""" - noise = np.random.random((10, 10)) - ita = fitting.make_ita(noise) - if not (ita.shape == (100, 100)): - raise AssertionError() - noise *= np.nan - ita = fitting.make_ita(noise) - if not (len(ita) == 0): - raise AssertionError() - - -def test_RB_bias(): - """Test RB_bias""" - data = np.random.random((4, 4)) - model = make_model() - bias = fitting.RB_bias(data, model) - if not (len(bias) == 3): - raise AssertionError() - - -def test_bias_correct(): - """test that bias_correct does things""" - data = np.random.random((4, 4)) - model = make_model() - fitting.bias_correct(model, data) - - if __name__ == "__main__": # introspect and run all the functions starting with 'test' for f in dir(): - if f.startswith('test'): + if f.startswith("test"): print(f) globals()[f]() diff --git a/tests/unit/test_mpi.py b/tests/unit/test_mpi.py new file mode 100644 index 00000000..da5497f7 --- /dev/null +++ b/tests/unit/test_mpi.py @@ -0,0 +1,15 @@ +#!/usr/bin/env python + +def test_import_mpi(): + try: + from AegeanTools.mpi import MPI_AVAIL, MPI + except Exception: + raise AssertionError("Error importing AegeanTools.mpi") + + + try: + import mpi4py + if MPI is None: + raise AssertionError("MPI should not be None when mpi4py is available") + except ImportError: + pass \ No newline at end of file diff --git a/tests/unit/test_source_finder.py b/tests/unit/test_source_finder.py index ef40c950..f64fa2af 100755 --- a/tests/unit/test_source_finder.py +++ b/tests/unit/test_source_finder.py @@ -159,6 +159,22 @@ def test_load_globals(): raise AssertionError() +def test_load_globals_cube(): + """Test load_globals""" + sfinder = sf.SourceFinder() + filename = "tests/test_files/synthetic_cube.fits" + # aux_files = sf.get_aux_files("tests/test_files/synthetic_cube.fits") + bkg = filename + rms = filename + sfinder.load_globals(filename, bkgin=bkg, rmsin=rms, as_cube=True) + if sfinder.img is None: + raise AssertionError() + if np.allclose(sfinder.img, 0): + print("True") + else: + raise AssertionError() + + def test_find_and_prior_sources(): """Test find sources and prior sources""" try: @@ -251,7 +267,8 @@ def test_find_and_prior_sources(): if not (os.path.exists("dlme")): raise AssertionError("Failed to create output file") finally: - os.remove("dlme") + if os.path.exists("dlme"): + os.remove("dlme") def dont_test_find_and_prior_parallel(): diff --git a/tests/unit/test_wcs_helpers.py b/tests/unit/test_wcs_helpers.py index b7278b3c..4f754938 100755 --- a/tests/unit/test_wcs_helpers.py +++ b/tests/unit/test_wcs_helpers.py @@ -267,6 +267,24 @@ def test_fix_aips_header(): _ = fix_aips_header(header) +def test_pix2freq(): + wcs = WCSHelper.from_file("tests/test_files/synthetic_cube.fits") + freq = wcs.pix2freq(8) + if freq != 170.0: + raise AssertionError( + f"The expected value is 223.76000000000005. However, {freq} was returned." + ) + + +def test_freq2pix(): + wcs = WCSHelper.from_file("tests/test_files/synthetic_cube.fits") + z = wcs.freq2pix(170.0) + if z != 8.0: + raise AssertionError( + f"The expected value is 7.999999999999999. However, {z} was returned." + ) + + if __name__ == "__main__": # introspect and run all the functions starting with 'test' for f in dir():