From 0a023ae9e6591818f44e0c3d7efb1d52bd4a0e71 Mon Sep 17 00:00:00 2001 From: Connor Stone Date: Tue, 5 Mar 2024 21:45:23 -0500 Subject: [PATCH 1/3] Add option to set init R, swap out PSF method --- autoprof/Pipeline.py | 289 +++++++++++------- autoprof/autoprofutils/Diagnostic_Plots.py | 145 ++++----- .../pipeline_steps/Isophote_Initialize.py | 175 +++++------ autoprof/pipeline_steps/PSF.py | 73 +++-- 4 files changed, 356 insertions(+), 326 deletions(-) diff --git a/autoprof/Pipeline.py b/autoprof/Pipeline.py index bb2b5988..eb12f25b 100644 --- a/autoprof/Pipeline.py +++ b/autoprof/Pipeline.py @@ -1,6 +1,7 @@ import sys import os from .pipeline_steps import * + # from .pipeline_steps.Plotting_Steps import Plot_Galaxy_Image # from .pipeline_steps.Background import Background_Mode, Background_DilatedSources, Background_Unsharp, Background_Basic # from .pipeline_steps.PSF import PSF_IRAF, PSF_StarFind, PSF_Image, PSF_deconvolve @@ -30,11 +31,13 @@ import warnings import traceback from astropy.io.fits.verify import VerifyWarning -warnings.simplefilter('ignore', category=VerifyWarning) + +warnings.simplefilter("ignore", category=VerifyWarning) + class Isophote_Pipeline(object): - def __init__(self, loggername = None): + def __init__(self, loggername=None): """ Initialize pipeline object, user can replace functions with their own if they want, otherwise defaults are used. @@ -42,53 +45,70 @@ def __init__(self, loggername = None): """ # Functions avaiable by default for building the pipeline - self.pipeline_methods = {'background': Background_Mode, - 'background dilatedsources': Background_DilatedSources, - 'background unsharp': Background_Unsharp, - 'background basic': Background_Basic, - 'psf': PSF_StarFind, - 'psf IRAF': PSF_IRAF, - 'psf img': PSF_Image, - 'psf deconvolve': PSF_deconvolve, - 'center': Center_HillClimb, - 'center mean': Center_HillClimb_mean, - 'center forced': Center_Forced, - 'center 2DGaussian': Center_2DGaussian, - 'center 1DGaussian': Center_1DGaussian, - 'center OfMass': Center_OfMass, - 'crop': Crop, - 'isophoteinit': Isophote_Initialize, - 'isophoteinit forced': Isophote_Init_Forced, - 'isophoteinit mean': Isophote_Initialize_mean, - 'plot image': Plot_Galaxy_Image, - 'writefi': WriteFi, - 'isophotefit': Isophote_Fit_FFT_Robust, - 'isophotefit fixed': Isophote_Fit_FixedPhase, - 'isophotefit mean': Isophote_Fit_FFT_mean, - 'isophotefit forced': Isophote_Fit_Forced, - 'isophotefit photutils': Photutils_Fit, - 'mask badpixels': Bad_Pixel_Mask, - 'starmask': Star_Mask, - 'starmask IRAF': Star_Mask_IRAF, - 'mask segmentation map': Mask_Segmentation_Map, - 'isophoteextract': Isophote_Extract, - 'isophoteextract photutils': Isophote_Extract_Photutils, - 'isophoteextract forced': Isophote_Extract_Forced, - 'checkfit': Check_Fit, - 'writeprof': WriteProf, - 'ellipsemodel': EllipseModel, - 'radialprofiles': Radial_Profiles, - 'sliceprofile': Slice_Profile, - 'axialprofiles': Axial_Profiles} - + self.pipeline_methods = { + "background": Background_Mode, + "background dilatedsources": Background_DilatedSources, + "background unsharp": Background_Unsharp, + "background basic": Background_Basic, + "psf": PSF_Assumed, + "psf starfind": PSF_StarFind, + "psf IRAF": PSF_IRAF, + "psf img": PSF_Image, + "psf deconvolve": PSF_deconvolve, + "center": Center_HillClimb, + "center mean": Center_HillClimb_mean, + "center forced": Center_Forced, + "center 2DGaussian": Center_2DGaussian, + "center 1DGaussian": Center_1DGaussian, + "center OfMass": Center_OfMass, + "crop": Crop, + "isophoteinit": Isophote_Initialize, + "isophoteinit forced": Isophote_Init_Forced, + "isophoteinit mean": Isophote_Initialize_mean, + "plot image": Plot_Galaxy_Image, + "writefi": WriteFi, + "isophotefit": Isophote_Fit_FFT_Robust, + "isophotefit fixed": Isophote_Fit_FixedPhase, + "isophotefit mean": Isophote_Fit_FFT_mean, + "isophotefit forced": Isophote_Fit_Forced, + "isophotefit photutils": Photutils_Fit, + "mask badpixels": Bad_Pixel_Mask, + "starmask": Star_Mask, + "starmask IRAF": Star_Mask_IRAF, + "mask segmentation map": Mask_Segmentation_Map, + "isophoteextract": Isophote_Extract, + "isophoteextract photutils": Isophote_Extract_Photutils, + "isophoteextract forced": Isophote_Extract_Forced, + "checkfit": Check_Fit, + "writeprof": WriteProf, + "ellipsemodel": EllipseModel, + "radialprofiles": Radial_Profiles, + "sliceprofile": Slice_Profile, + "axialprofiles": Axial_Profiles, + } + # Default pipeline analysis order - self.pipeline_steps = {'head': ['background', 'psf', 'center', 'isophoteinit', - 'isophotefit', 'isophoteextract', 'checkfit', 'writeprof']} - + self.pipeline_steps = { + "head": [ + "background", + "psf", + "center", + "isophoteinit", + "isophotefit", + "isophoteextract", + "checkfit", + "writeprof", + ] + } + # Start the logger - logging.basicConfig(level=logging.INFO, filename = 'AutoProf.log' if loggername is None else loggername, filemode = 'w') + logging.basicConfig( + level=logging.INFO, + filename="AutoProf.log" if loggername is None else loggername, + filemode="w", + ) - def UpdatePipeline(self, new_pipeline_methods = None, new_pipeline_steps = None): + def UpdatePipeline(self, new_pipeline_methods=None, new_pipeline_steps=None): """ modify steps in the AutoProf pipeline. @@ -100,17 +120,19 @@ def UpdatePipeline(self, new_pipeline_methods = None, new_pipeline_steps = None) steps as values, the corresponding steps will be replaced. """ if new_pipeline_methods: - logging.info('PIPELINE updating these pipeline methods: %s' % str(new_pipeline_methods.keys())) + logging.info( + "PIPELINE updating these pipeline methods: %s" % str(new_pipeline_methods.keys()) + ) self.pipeline_methods.update(new_pipeline_methods) if new_pipeline_steps: - logging.info('PIPELINE new steps: %s' % (str(new_pipeline_steps))) + logging.info("PIPELINE new steps: %s" % (str(new_pipeline_steps))) if type(new_pipeline_steps) == list: - self.pipeline_steps['head'] = new_pipeline_steps + self.pipeline_steps["head"] = new_pipeline_steps elif type(new_pipeline_steps) == dict: - assert 'head' in new_pipeline_steps.keys() + assert "head" in new_pipeline_steps.keys() self.pipeline_steps = new_pipeline_steps - - def Process_Image(self, options = {}): + + def Process_Image(self, options={}): """ Function which runs the pipeline for a single image. Each sub-function of the pipeline is run in order and the outputs are passed along. If multiple images are given, the pipeline is @@ -123,46 +145,66 @@ def Process_Image(self, options = {}): for key in list(options.keys()): if options[key] is None: del options[key] - + # Seed the random number generator in numpy so each thread gets unique random numbers try: sleep(0.01) - np.random.seed(int(np.random.randint(10000)*current_process().pid*(time() % 1) % 2**15)) + np.random.seed( + int(np.random.randint(10000) * current_process().pid * (time() % 1) % 2**15) + ) except: pass - + # use filename if no name is given - if not ('ap_name' in options and type(options['ap_name']) == str): - base = os.path.split(options['ap_image_file'])[1] - options['ap_name'] = os.path.splitext(base)[0] + if not ("ap_name" in options and type(options["ap_name"]) == str): + base = os.path.split(options["ap_image_file"])[1] + options["ap_name"] = os.path.splitext(base)[0] # Read the primary image try: - dat = Read_Image(options['ap_image_file'], options) + dat = Read_Image(options["ap_image_file"], options) except: - logging.error('%s: could not read image %s' % (options['ap_name'], options['ap_image_file'])) + logging.error( + "%s: could not read image %s" % (options["ap_name"], options["ap_image_file"]) + ) return 1 - + # Check that image data exists and is not corrupted - if dat is None or np.all(dat[int(len(dat)/2.)-10:int(len(dat)/2.)+10, int(len(dat[0])/2.)-10:int(len(dat[0])/2.)+10] == 0): - logging.error('%s Large chunk of data missing, impossible to process image' % options['ap_name']) + if dat is None or np.all( + dat[ + int(len(dat) / 2.0) - 10 : int(len(dat) / 2.0) + 10, + int(len(dat[0]) / 2.0) - 10 : int(len(dat[0]) / 2.0) + 10, + ] + == 0 + ): + logging.error( + "%s Large chunk of data missing, impossible to process image" % options["ap_name"] + ) return 1 - + # Track time to run analysis start = time() - + # Run the Pipeline timers = {} results = {} - key = 'head' + key = "head" step = 0 while step < len(self.pipeline_steps[key]): try: - logging.info('%s: %s %s at: %.1f sec' % (options['ap_name'], key, self.pipeline_steps[key][step], time() - start)) - print('%s: %s %s at: %.1f sec' % (options['ap_name'], key, self.pipeline_steps[key][step], time() - start)) - if 'branch' in self.pipeline_steps[key][step]: - decision, newoptions = self.pipeline_methods[self.pipeline_steps[key][step]](dat, results, options) + logging.info( + "%s: %s %s at: %.1f sec" + % (options["ap_name"], key, self.pipeline_steps[key][step], time() - start) + ) + print( + "%s: %s %s at: %.1f sec" + % (options["ap_name"], key, self.pipeline_steps[key][step], time() - start) + ) + if "branch" in self.pipeline_steps[key][step]: + decision, newoptions = self.pipeline_methods[self.pipeline_steps[key][step]]( + dat, results, options + ) options.update(newoptions) if type(decision) == str: key = decision @@ -171,51 +213,63 @@ def Process_Image(self, options = {}): step += 1 else: step_start = time() - dat, res = self.pipeline_methods[self.pipeline_steps[key][step]](dat, results, options) + dat, res = self.pipeline_methods[self.pipeline_steps[key][step]]( + dat, results, options + ) results.update(res) timers[self.pipeline_steps[key][step]] = time() - step_start step += 1 except Exception as e: - logging.error('%s: on step %s got error: %s' % (options['ap_name'], self.pipeline_steps[key][step], str(e))) - logging.error('%s: with full trace: %s' % (options['ap_name'], traceback.format_exc())) + logging.error( + "%s: on step %s got error: %s" + % (options["ap_name"], self.pipeline_steps[key][step], str(e)) + ) + logging.error( + "%s: with full trace: %s" % (options["ap_name"], traceback.format_exc()) + ) return 1 - - print('%s: Processing Complete! (at %.1f sec)' % (options['ap_name'], time() - start)) - logging.info('%s: Processing Complete! (at %.1f sec)' % (options['ap_name'], time() - start)) + + print("%s: Processing Complete! (at %.1f sec)" % (options["ap_name"], time() - start)) + logging.info( + "%s: Processing Complete! (at %.1f sec)" % (options["ap_name"], time() - start) + ) return timers - + def Process_List(self, options): """ Wrapper function to run "Process_Image" in parallel for many images. """ - assert type(options['ap_image_file']) == list - + assert type(options["ap_image_file"]) == list + # Format the inputs so that they can be zipped with the images files # and passed to the Process_Image function. use_options = [] - for i in range(len(options['ap_image_file'])): + for i in range(len(options["ap_image_file"])): tmp_options = {} for k in options.keys(): - if type(options[k]) == list and not k in ['ap_new_pipeline_steps']: + if type(options[k]) == list and not k in ["ap_new_pipeline_steps"]: tmp_options[k] = options[k][i] else: tmp_options[k] = options[k] use_options.append(tmp_options) # Track how long it takes to run the analysis start = time() - + # Create a multiprocessing pool to parallelize image processing - n_procs = options['ap_n_procs'] if 'ap_n_procs' in options else 1 + n_procs = options["ap_n_procs"] if "ap_n_procs" in options else 1 if n_procs > 1: with Pool(int(n_procs)) as pool: - res = pool.map(self.Process_Image, use_options, - chunksize = 5 if len(options['ap_image_file']) > 100 else 1) + res = pool.map( + self.Process_Image, + use_options, + chunksize=5 if len(options["ap_image_file"]) > 100 else 1, + ) else: res = list(map(self.Process_Image, use_options)) - + # Report completed processing, and track time used - logging.info('All Images Finished Processing at %.1f' % (time() - start)) + logging.info("All Images Finished Processing at %.1f" % (time() - start)) timers = dict() counts = dict() for r in res: @@ -224,20 +278,20 @@ def Process_List(self, options): for s in r.keys(): if s in timers: timers[s] += r[s] - counts[s] += 1. + counts[s] += 1.0 else: timers[s] = r[s] - counts[s] = 1. + counts[s] = 1.0 if len(timers) == 0: - logging.error('All images failed to process!') - return + logging.error("All images failed to process!") + return for s in timers: timers[s] /= counts[s] - logging.info('%s took %.3f seconds on average' % (s, timers[s])) - + logging.info("%s took %.3f seconds on average" % (s, timers[s])) + # Return the success/fail indicators for every Process_Image excecution return res - + def Process_ConfigFile(self, config_file): """ Reads in a configuration file and sets parameters for the pipeline. The configuration @@ -247,44 +301,53 @@ def Process_ConfigFile(self, config_file): returns: timing of each pipeline step if successful. Else returns 1 """ - + # Import the config file regardless of where it is from - if '/' in config_file: - startat = config_file.rfind('/')+1 + if "/" in config_file: + startat = config_file.rfind("/") + 1 else: startat = 0 - if '.' in config_file: - use_config = config_file[startat:config_file.rfind('.')] + if "." in config_file: + use_config = config_file[startat : config_file.rfind(".")] else: use_config = config_file[startat:] - if '/' in config_file: - sys.path.append(config_file[:config_file.rfind('/')]) + if "/" in config_file: + sys.path.append(config_file[: config_file.rfind("/")]) try: c = importlib.import_module(use_config) except: sys.path.append(os.getcwd()) c = importlib.import_module(use_config) - if 'forced' in c.ap_process_mode: - self.UpdatePipeline(new_pipeline_steps = ['background', 'psf', 'center forced', 'isophoteinit forced', - 'isophoteextract forced', 'writeprof']) - + if "forced" in c.ap_process_mode: + self.UpdatePipeline( + new_pipeline_steps=[ + "background", + "psf", + "center forced", + "isophoteinit forced", + "isophoteextract forced", + "writeprof", + ] + ) + try: - self.UpdatePipeline(new_pipeline_methods = c.ap_new_pipeline_methods) + self.UpdatePipeline(new_pipeline_methods=c.ap_new_pipeline_methods) except: pass try: - self.UpdatePipeline(new_pipeline_steps = c.ap_new_pipeline_steps) + self.UpdatePipeline(new_pipeline_steps=c.ap_new_pipeline_steps) except: pass - + use_options = GetOptions(c) - - if c.ap_process_mode in ['image', 'forced image']: + + if c.ap_process_mode in ["image", "forced image"]: return self.Process_Image(use_options) - elif c.ap_process_mode in ['image list', 'forced image list']: + elif c.ap_process_mode in ["image list", "forced image list"]: return self.Process_List(use_options) else: - logging.error('Unrecognized process_mode! Should be in: [image, image list, forced image, forced image list]') + logging.error( + "Unrecognized process_mode! Should be in: [image, image list, forced image, forced image list]" + ) return 1 - diff --git a/autoprof/autoprofutils/Diagnostic_Plots.py b/autoprof/autoprofutils/Diagnostic_Plots.py index ded63485..db172c18 100644 --- a/autoprof/autoprofutils/Diagnostic_Plots.py +++ b/autoprof/autoprofutils/Diagnostic_Plots.py @@ -45,14 +45,11 @@ "Plot_EllipseModel", ) + def Plot_Background(values, bkgrnd, noise, results, options): hist, bins = np.histogram( - values[ - np.logical_and( - (values - bkgrnd) < 20 * noise, (values - bkgrnd) > -5 * noise - ) - ], + values[np.logical_and((values - bkgrnd) < 20 * noise, (values - bkgrnd) > -5 * noise)], bins=max(10, int(np.sqrt(len(values)) / 2)), ) plt.figure(figsize=(5, 5)) @@ -90,23 +87,23 @@ def Plot_Background(values, bkgrnd, noise, results, options): plt.close() -def Plot_PSF_Stars( - IMG, stars_x, stars_y, stars_fwhm, psf, results, options, flagstars=None -): +def Plot_PSF_Stars(IMG, stars_x, stars_y, stars_fwhm, psf, results, options, flagstars=None): LSBImage(IMG - results["background"], results["background noise"]) - AddScale(plt.gca(), IMG.shape[0]*options['ap_pixscale']) + AddScale(plt.gca(), IMG.shape[0] * options["ap_pixscale"]) for i in range(len(stars_fwhm)): plt.gca().add_patch( Ellipse( - xy = (stars_x[i], stars_y[i]), - width = 20 * psf, - height = 20 * psf, - angle = 0, + xy=(stars_x[i], stars_y[i]), + width=20 * psf, + height=20 * psf, + angle=0, fill=False, linewidth=1.5, - color=autocolours["red1"] - if not flagstars is None and flagstars[i] - else autocolours["blue1"], + color=( + autocolours["red1"] + if not flagstars is None and flagstars[i] + else autocolours["blue1"] + ), ) ) plt.tight_layout() @@ -126,15 +123,11 @@ def Plot_Isophote_Init_Ellipse(dat, circ_ellipse_radii, ellip, phase, results, o ranges = [ [ max(0, int(results["center"]["x"] - circ_ellipse_radii[-1] * 1.5)), - min( - dat.shape[1], int(results["center"]["x"] + circ_ellipse_radii[-1] * 1.5) - ), + min(dat.shape[1], int(results["center"]["x"] + circ_ellipse_radii[-1] * 1.5)), ], [ max(0, int(results["center"]["y"] - circ_ellipse_radii[-1] * 1.5)), - min( - dat.shape[0], int(results["center"]["y"] + circ_ellipse_radii[-1] * 1.5) - ), + min(dat.shape[0], int(results["center"]["y"] + circ_ellipse_radii[-1] * 1.5)), ], ] @@ -142,16 +135,16 @@ def Plot_Isophote_Init_Ellipse(dat, circ_ellipse_radii, ellip, phase, results, o dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]], results["background noise"], ) - AddScale(plt.gca(), (ranges[0][1] - ranges[0][0])*options['ap_pixscale']) + AddScale(plt.gca(), (ranges[0][1] - ranges[0][0]) * options["ap_pixscale"]) plt.gca().add_patch( Ellipse( - xy =( + xy=( results["center"]["x"] - ranges[0][0], results["center"]["y"] - ranges[1][0], ), - width = 2 * circ_ellipse_radii[-1], - height = 2 * circ_ellipse_radii[-1] * (1.0 - ellip), - angle = phase * 180 / np.pi, + width=2 * circ_ellipse_radii[-1], + height=2 * circ_ellipse_radii[-1] * (1.0 - ellip), + angle=phase * 180 / np.pi, fill=False, linewidth=1, color=autocolours["blue1"], @@ -229,11 +222,7 @@ def Plot_Isophote_Fit(dat, sample_radii, parameters, results, options): Rlim = sample_radii[-1] * ( 1.0 if parameters[-1]["m"] is None - else np.exp( - sum( - np.abs(parameters[-1]["Am"][m]) for m in range(len(parameters[-1]["m"])) - ) - ) + else np.exp(sum(np.abs(parameters[-1]["Am"][m]) for m in range(len(parameters[-1]["m"])))) ) ranges = [ [ @@ -249,7 +238,7 @@ def Plot_Isophote_Fit(dat, sample_radii, parameters, results, options): dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]], results["background noise"], ) - AddScale(plt.gca(), (ranges[0][1] - ranges[0][0])*options['ap_pixscale']) + AddScale(plt.gca(), (ranges[0][1] - ranges[0][0]) * options["ap_pixscale"]) for i in range(len(sample_radii)): N = max(15, int(0.9 * 2 * np.pi * sample_radii[i])) @@ -278,7 +267,7 @@ def Plot_Isophote_Fit(dat, sample_radii, parameters, results, options): R * X + results["center"]["x"] - ranges[0][0], R * Y + results["center"]["y"] - ranges[1][0], ) - + plt.plot( list(X) + [X[0]], list(Y) + [Y[0]], @@ -304,11 +293,7 @@ def _Plot_Isophotes(dat, R, parameters, results, options): Rlim = R[-1] * ( 1.0 if parameters[-1]["m"] is None - else np.exp( - sum( - np.abs(parameters[-1]["Am"][m]) for m in range(len(parameters[-1]["m"])) - ) - ) + else np.exp(sum(np.abs(parameters[-1]["Am"][m]) for m in range(len(parameters[-1]["m"])))) ) ranges = [ [ @@ -324,7 +309,7 @@ def _Plot_Isophotes(dat, R, parameters, results, options): dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]], results["background noise"], ) - AddScale(plt.gca(), (ranges[0][1] - ranges[0][0])*options['ap_pixscale']) + AddScale(plt.gca(), (ranges[0][1] - ranges[0][0]) * options["ap_pixscale"]) fitlim = results["fit R"][-1] if "fit R" in results else np.inf for i in range(len(R)): N = max(15, int(0.9 * 2 * np.pi * R[i])) @@ -380,14 +365,14 @@ def Plot_SB_Profile(dat, R, SB, SB_e, parameters, results, options): errscale = 1.0 if np.all(SB_e[CHOOSE] < 0.5): errscale = 1 / np.max(SB_e[CHOOSE]) - if 'ap_plot_sbprof_set_errscale' in options: - errscale = options['ap_plot_sbprof_set_errscale'] - + if "ap_plot_sbprof_set_errscale" in options: + errscale = options["ap_plot_sbprof_set_errscale"] + lnlist = [] - if errscale > 1.01 or 'ap_plot_sbprof_set_errscale' in options: - errlabel = ' (err$\\times$%.1f)' % errscale + if errscale > 1.01 or "ap_plot_sbprof_set_errscale" in options: + errlabel = " (err$\\times$%.1f)" % errscale else: - errlabel = '' + errlabel = "" lnlist.append( plt.errorbar( R[CHOOSE], @@ -413,12 +398,12 @@ def Plot_SB_Profile(dat, R, SB, SB_e, parameters, results, options): ) plt.xlabel("Semi-Major-Axis [arcsec]", fontsize=16) plt.ylabel("Surface Brightness [mag arcsec$^{-2}$]", fontsize=16) - if 'ap_plot_sbprof_xlim' in options: - plt.xlim(options['ap_plot_sbprof_xlim']) + if "ap_plot_sbprof_xlim" in options: + plt.xlim(options["ap_plot_sbprof_xlim"]) else: plt.xlim([0, None]) - if 'ap_plot_sbprof_ylim' in options: - plt.ylim(options['ap_plot_sbprof_ylim']) + if "ap_plot_sbprof_ylim" in options: + plt.ylim(options["ap_plot_sbprof_ylim"]) bkgrdnoise = ( -2.5 * np.log10(results["background noise"]) + zeropoint @@ -491,12 +476,12 @@ def Plot_I_Profile(dat, R, I, I_e, parameters, results, options): plt.xlabel("Semi-Major-Axis [arcsec]", fontsize=16) plt.ylabel("Intensity [flux arcsec$^{-2}$]", fontsize=16) plt.yscale("log") - if 'ap_plot_sbprof_xlim' in options: - plt.xlim(options['ap_plot_sbprof_xlim']) + if "ap_plot_sbprof_xlim" in options: + plt.xlim(options["ap_plot_sbprof_xlim"]) else: plt.xlim([0, None]) - if 'ap_plot_sbprof_ylim' in options: - plt.ylim(options['ap_plot_sbprof_ylim']) + if "ap_plot_sbprof_ylim" in options: + plt.ylim(options["ap_plot_sbprof_ylim"]) bkgrdnoise = results["background noise"] / (options["ap_pixscale"] ** 2) lnlist.append( plt.axhline( @@ -569,11 +554,8 @@ def Plot_Phase_Profile(R, parameters, results, options): ) plt.plot( R, - list( - p["Phim"][m] / (np.pi * parameters[0]["m"][m]) for p in parameters - ), - label="$\\phi_%i$ [rad/%i$\\pi$]" - % (parameters[0]["m"][m], parameters[0]["m"][m]), + list(p["Phim"][m] / (np.pi * parameters[0]["m"][m]) for p in parameters), + label="$\\phi_%i$ [rad/%i$\\pi$]" % (parameters[0]["m"][m], parameters[0]["m"][m]), ) plt.legend() plt.xlabel("Semi-Major-Axis [arcsec]", fontsize=16) @@ -650,9 +632,7 @@ def Plot_Meas_Fmodes(R, parameters, results, options): plt.close() -def Plot_Radial_Profiles( - dat, R, sb, sbE, pa, nwedges, wedgeangles, wedgewidth, results, options -): +def Plot_Radial_Profiles(dat, R, sb, sbE, pa, nwedges, wedgeangles, wedgewidth, results, options): zeropoint = options["ap_zeropoint"] if "ap_zeropoint" in options else 22.5 ranges = [ @@ -665,7 +645,7 @@ def Plot_Radial_Profiles( min(dat.shape[0], int(results["center"]["y"] + 1.5 * R[-1] + 2)), ], ] - + cmap = cm.get_cmap("hsv") colorind = (np.linspace(0, 1 - 1 / nwedges, nwedges) + 0.1) % 1.0 for sa_i in range(len(wedgeangles)): @@ -714,7 +694,7 @@ def Plot_Radial_Profiles( dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]], results["background noise"], ) - AddScale(plt.gca(), (ranges[0][1] - ranges[0][0])*options['ap_pixscale']) + AddScale(plt.gca(), (ranges[0][1] - ranges[0][0]) * options["ap_pixscale"]) cx, cy = ( results["center"]["x"] - ranges[0][0], @@ -787,15 +767,11 @@ def Plot_Axial_Profiles(dat, R, sb, sbE, pa, results, options): for ang in [1, -1]: key = (rd, ang) # cmap = matplotlib.cm.get_cmap('viridis_r') - norm = matplotlib.colors.Normalize( - vmin=0, vmax=R[-1] * options["ap_pixscale"] - ) + norm = matplotlib.colors.Normalize(vmin=0, vmax=R[-1] * options["ap_pixscale"]) for pi, pR in enumerate(R): if pi % 3 != 0: continue - CHOOSE = np.logical_and( - np.array(sb[key][pi]) < 99, np.array(sbE[key][pi]) < 1 - ) + CHOOSE = np.logical_and(np.array(sb[key][pi]) < 99, np.array(sbE[key][pi]) < 1) plt.errorbar( np.array(R)[CHOOSE] * options["ap_pixscale"], np.array(sb[key][pi])[CHOOSE], @@ -810,22 +786,20 @@ def Plot_Axial_Profiles(dat, R, sb, sbE, pa, results, options): "%s-axis position on line [arcsec]" % ( "Major" - if "ap_axialprof_parallel" in options - and options["ap_axialprof_parallel"] + if "ap_axialprof_parallel" in options and options["ap_axialprof_parallel"] else "Minor" ), fontsize=16, ) plt.ylabel("Surface Brightness [mag arcsec$^{-2}$]", fontsize=16) cb1 = plt.colorbar( - matplotlib.cm.ScalarMappable(norm=norm, cmap=autocmap.reversed()) + matplotlib.cm.ScalarMappable(norm=norm, cmap=autocmap.reversed()), ax=plt.gca() ) cb1.set_label( "%s-axis position of line [arcsec]" % ( "Minor" - if "ap_axialprof_parallel" in options - and options["ap_axialprof_parallel"] + if "ap_axialprof_parallel" in options and options["ap_axialprof_parallel"] else "Major" ), fontsize=16, @@ -866,9 +840,7 @@ def Plot_Axial_Profiles(dat, R, sb, sbE, pa, results, options): firstbad = np.argmax(np.logical_not(CHOOSE)) if firstbad > 3: CHOOSE[firstbad:] = False - outto = ( - np.array(results["prof data"]["R"])[CHOOSE][-1] * 1.5 / options["ap_pixscale"] - ) + outto = np.array(results["prof data"]["R"])[CHOOSE][-1] * 1.5 / options["ap_pixscale"] ranges = [ [ max(0, int(results["center"]["x"] - outto - 2)), @@ -883,11 +855,11 @@ def Plot_Axial_Profiles(dat, R, sb, sbE, pa, results, options): dat[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]], results["background noise"], ) - AddScale(plt.gca(), (ranges[0][1] - ranges[0][0])*options['ap_pixscale']) + AddScale(plt.gca(), (ranges[0][1] - ranges[0][0]) * options["ap_pixscale"]) count = 0 cmap = matplotlib.cm.get_cmap("hsv") colorind = (np.linspace(0, 1 - 1 / 4, 4) + 0.1) % 1 - colours = list(cmap(c) for c in colorind) + colours = list(cmap(c) for c in colorind) for rd in [1, -1]: for ang in [1, -1]: key = (rd, ang) @@ -912,11 +884,10 @@ def Plot_Axial_Profiles(dat, R, sb, sbE, pa, results, options): linewidth=0.5, color=colours[count], label=( - "%sR : pa%s90" - % ("+" if rd > 0 else "-", "+" if ang > 0 else "-") - ) - if pi == 0 - else None, + ("%sR : pa%s90" % ("+" if rd > 0 else "-", "+" if ang > 0 else "-")) + if pi == 0 + else None + ), ) count += 1 plt.legend() @@ -949,9 +920,7 @@ def Plot_EllipseModel(IMG, Model, R, modeltype, results, options): plt.figure(figsize=(7, 7)) autocmap.set_under("k", alpha=0) showmodel = Model[ranges[1][0] : ranges[1][1], ranges[0][0] : ranges[0][1]].copy() - showmodel[showmodel > 0] += np.max(showmodel) / (10 ** 3.5) - np.min( - showmodel[showmodel > 0] - ) + showmodel[showmodel > 0] += np.max(showmodel) / (10**3.5) - np.min(showmodel[showmodel > 0]) plt.imshow( showmodel, origin="lower", diff --git a/autoprof/pipeline_steps/Isophote_Initialize.py b/autoprof/pipeline_steps/Isophote_Initialize.py index 44bf5236..aea4547e 100644 --- a/autoprof/pipeline_steps/Isophote_Initialize.py +++ b/autoprof/pipeline_steps/Isophote_Initialize.py @@ -34,6 +34,7 @@ __all__ = ("Isophote_Init_Forced", "Isophote_Initialize", "Isophote_Initialize_mean") + def Isophote_Init_Forced(IMG, results, options): """Read global elliptical isophote to a galaxy from an aux file. @@ -79,9 +80,7 @@ def Isophote_Init_Forced(IMG, results, options): pa = ( PA_shift_convention( float( - line[ - line.find("pa:") + 3 : line.find("+-", line.find("pa:")) - ].strip() + line[line.find("pa:") + 3 : line.find("+-", line.find("pa:"))].strip() ), deg=True, ) @@ -89,30 +88,21 @@ def Isophote_Init_Forced(IMG, results, options): / 180 ) pa_err = ( - float( - line[ - line.find("+-", line.find("pa:")) + 2 : line.find("deg") - ].strip() - ) + float(line[line.find("+-", line.find("pa:")) + 2 : line.find("deg")].strip()) * np.pi / 180 ) R = float( - line[ - line.find("size:") + 5 : line.find("pix", line.find("size:")) - ].strip() + line[line.find("size:") + 5 : line.find("pix", line.find("size:"))].strip() ) break - auxmessage = ( - "global ellipticity: %.3f +- %.3f, pa: %.3f +- %.3f deg, size: %f pix" - % ( - ellip, - ellip_err, - PA_shift_convention(pa) * 180 / np.pi, - pa_err * 180 / np.pi, - R, - ) + auxmessage = "global ellipticity: %.3f +- %.3f, pa: %.3f +- %.3f deg, size: %f pix" % ( + ellip, + ellip_err, + PA_shift_convention(pa) * 180 / np.pi, + pa_err * 180 / np.pi, + R, ) return IMG, { @@ -207,34 +197,54 @@ def Isophote_Initialize(IMG, results, options): if not np.any(mask): mask = None - while circ_ellipse_radii[-1] < (len(IMG) / 2): - circ_ellipse_radii.append(circ_ellipse_radii[-1] * (1 + 0.2)) - isovals = _iso_extract( - dat, - circ_ellipse_radii[-1], - {"ellip": 0.0, "pa": 0.0}, - results["center"], - more=True, - mask=mask, - sigmaclip=True, - sclip_nsigma=3, - interp_mask=True, - ) - coefs = fft(isovals[0]) - allphase.append(coefs[2]) - # Stop when at 3 time background noise - if ( - np.quantile(isovals[0], 0.8) - < ( - (options["ap_fit_limit"] + 1 if "ap_fit_limit" in options else 3) - * results["background noise"] + if "ap_isoinit_R_set" in options: + circ_ellipse_radii = np.logspace( + np.log10(circ_ellipse_radii[0]), + np.log10(options["ap_isoinit_R_set"]), + 10, + ).tolist() + for r in circ_ellipse_radii: + isovals = _iso_extract( + dat, + r, + {"ellip": 0.0, "pa": 0.0}, + results["center"], + more=True, + mask=mask, + sigmaclip=True, + sclip_nsigma=3, + interp_mask=True, ) - and len(circ_ellipse_radii) > 4 - ): - break - logging.info( - "%s: init scale: %f pix" % (options["ap_name"], circ_ellipse_radii[-1]) - ) + coefs = fft(isovals[0]) + allphase.append(coefs[2]) + else: + while circ_ellipse_radii[-1] < (len(IMG) / 2): + circ_ellipse_radii.append(circ_ellipse_radii[-1] * (1 + 0.2)) + isovals = _iso_extract( + dat, + circ_ellipse_radii[-1], + {"ellip": 0.0, "pa": 0.0}, + results["center"], + more=True, + mask=mask, + sigmaclip=True, + sclip_nsigma=3, + interp_mask=True, + ) + coefs = fft(isovals[0]) + allphase.append(coefs[2]) + # Stop when at 3 time background noise + if ( + np.quantile(isovals[0], 0.8) + < ( + (options["ap_fit_limit"] + 1 if "ap_fit_limit" in options else 3) + * results["background noise"] + ) + and len(circ_ellipse_radii) > 4 + ): + break + + logging.info("%s: init scale: %f pix" % (options["ap_name"], circ_ellipse_radii[-1])) # Find global position angle. phase = (-Angle_Median(np.angle(allphase[-5:])) / 2) % np.pi if "ap_isoinit_pa_set" in options: @@ -336,15 +346,11 @@ def Isophote_Initialize(IMG, results, options): ) coefs = fft(isovals[0]) errallphase.append(coefs[2]) - sample_pas = ( - -np.angle(1j * np.array(errallphase) / np.mean(errallphase)) / 2 - ) % np.pi + sample_pas = (-np.angle(1j * np.array(errallphase) / np.mean(errallphase)) / 2) % np.pi pa_err = iqr(sample_pas, rng=[16, 84]) / 2 res_multi = map( lambda rrp: minimize( - lambda e, d, r, p, c, n, m: _fitEllip_loss( - _x_to_eps(e[0]), d, r, p, c, n, m - ), + lambda e, d, r, p, c, n, m: _fitEllip_loss(_x_to_eps(e[0]), d, r, p, c, n, m), x0=_inv_x_to_eps(ellip), args=( dat, @@ -369,9 +375,7 @@ def Isophote_Initialize(IMG, results, options): circ_ellipse_radii = np.array(circ_ellipse_radii) if "ap_doplot" in options and options["ap_doplot"]: - Plot_Isophote_Init_Ellipse( - dat, circ_ellipse_radii, ellip, phase, results, options - ) + Plot_Isophote_Init_Ellipse(dat, circ_ellipse_radii, ellip, phase, results, options) Plot_Isophote_Init_Optimize( circ_ellipse_radii, allphase, @@ -385,15 +389,12 @@ def Isophote_Initialize(IMG, results, options): options, ) - auxmessage = ( - "global ellipticity: %.3f +- %.3f, pa: %.3f +- %.3f deg, size: %f pix" - % ( - ellip, - ellip_err, - PA_shift_convention(phase) * 180 / np.pi, - pa_err * 180 / np.pi, - circ_ellipse_radii[-2], - ) + auxmessage = "global ellipticity: %.3f +- %.3f, pa: %.3f +- %.3f deg, size: %f pix" % ( + ellip, + ellip_err, + PA_shift_convention(phase) * 180 / np.pi, + pa_err * 180 / np.pi, + circ_ellipse_radii[-2], ) return IMG, { "init ellip": ellip, @@ -469,14 +470,9 @@ def Isophote_Initialize_mean(IMG, results, options): coefs = fft(isovals[0]) allphase.append(coefs[2]) # Stop when at 3 times background noise - if ( - np.mean(isovals[0]) < (3 * results["background noise"]) - and len(circ_ellipse_radii) > 4 - ): + if np.mean(isovals[0]) < (3 * results["background noise"]) and len(circ_ellipse_radii) > 4: break - logging.info( - "%s: init scale: %f pix" % (options["ap_name"], circ_ellipse_radii[-1]) - ) + logging.info("%s: init scale: %f pix" % (options["ap_name"], circ_ellipse_radii[-1])) # Find global position angle. phase = ( -Angle_Median(np.angle(allphase[-5:])) / 2 @@ -561,20 +557,14 @@ def Isophote_Initialize_mean(IMG, results, options): ) errallphase = [] for rr in RR: - isovals = _iso_extract( - dat, rr, {"ellip": 0.0, "pa": 0.0}, results["center"], more=True - ) + isovals = _iso_extract(dat, rr, {"ellip": 0.0, "pa": 0.0}, results["center"], more=True) coefs = fft(isovals[0]) errallphase.append(coefs[2]) - sample_pas = ( - -np.angle(1j * np.array(errallphase) / np.mean(errallphase)) / 2 - ) % np.pi + sample_pas = (-np.angle(1j * np.array(errallphase) / np.mean(errallphase)) / 2) % np.pi pa_err = np.std(sample_pas) res_multi = map( lambda rrp: minimize( - lambda e, d, r, p, c, n: _fitEllip_mean_loss( - _x_to_eps(e[0]), d, r, p, c, n - ), + lambda e, d, r, p, c, n: _fitEllip_mean_loss(_x_to_eps(e[0]), d, r, p, c, n), x0=_inv_x_to_eps(ellip), args=(dat, rrp[0], rrp[1], results["center"], results["background noise"]), method="Nelder-Mead", @@ -618,13 +608,13 @@ def Isophote_Initialize_mean(IMG, results, options): # origin = 'lower', cmap = 'Greys_r', norm = ImageNormalize(stretch=LogStretch())) plt.gca().add_patch( Ellipse( - xy = ( + xy=( results["center"]["x"] - ranges[0][0], results["center"]["y"] - ranges[1][0], ), - width = 2 * circ_ellipse_radii[-1], - height = 2 * circ_ellipse_radii[-1] * (1.0 - ellip), - angle = phase * 180 / np.pi, + width=2 * circ_ellipse_radii[-1], + height=2 * circ_ellipse_radii[-1] * (1.0 - ellip), + angle=phase * 180 / np.pi, fill=False, linewidth=1, color="y", @@ -673,15 +663,12 @@ def Isophote_Initialize_mean(IMG, results, options): ) plt.close() - auxmessage = ( - "global ellipticity: %.3f +- %.3f, pa: %.3f +- %.3f deg, size: %f pix" - % ( - ellip, - ellip_err, - PA_shift_convention(phase) * 180 / np.pi, - pa_err * 180 / np.pi, - circ_ellipse_radii[-2], - ) + auxmessage = "global ellipticity: %.3f +- %.3f, pa: %.3f +- %.3f deg, size: %f pix" % ( + ellip, + ellip_err, + PA_shift_convention(phase) * 180 / np.pi, + pa_err * 180 / np.pi, + circ_ellipse_radii[-2], ) return IMG, { "init ellip": ellip, diff --git a/autoprof/pipeline_steps/PSF.py b/autoprof/pipeline_steps/PSF.py index 5e7db083..0dc0cc62 100644 --- a/autoprof/pipeline_steps/PSF.py +++ b/autoprof/pipeline_steps/PSF.py @@ -27,7 +27,8 @@ from ..autoprofutils.Diagnostic_Plots import Plot_PSF_Stars from copy import deepcopy -__all__ = ("PSF_IRAF", "PSF_StarFind", "PSF_Image", "PSF_deconvolve") +__all__ = ("PSF_IRAF", "PSF_Assumed", "PSF_StarFind", "PSF_Image", "PSF_deconvolve") + def PSF_IRAF(IMG, results, options): """PSF routine which identifies stars and averages the FWHM. @@ -64,14 +65,11 @@ def PSF_IRAF(IMG, results, options): """ if "ap_set_psf" in options: - logging.info( - "%s: PSF set by user: %.4e" % (options["ap_name"], options["ap_set_psf"]) - ) + logging.info("%s: PSF set by user: %.4e" % (options["ap_name"], options["ap_set_psf"])) return IMG, {"psf fwhm": options["ap_set_psf"]} elif "ap_guess_psf" in options: logging.info( - "%s: PSF initialized by user: %.4e" - % (options["ap_name"], options["ap_guess_psf"]) + "%s: PSF initialized by user: %.4e" % (options["ap_name"], options["ap_guess_psf"]) ) fwhm_guess = options["ap_guess_psf"] else: @@ -120,6 +118,27 @@ def PSF_IRAF(IMG, results, options): return IMG, {"psf fwhm": psf, "auxfile psf": "psf fwhm: %.3f pix" % psf} +def PSF_Assumed(IMG, results, options): + """ + Most astronomical data is assumed to be nyquist sampled, thus we assume a + PSF scale of 4.0 pix to speed things up. Note that AutoProf just uses the + PSF FWHM to initialize some length scales on the image, accuracy is not very + important for the PSF. + """ + if "ap_set_psf" in options: + logging.info("%s: PSF set by user: %.4e pix" % (options["ap_name"], options["ap_set_psf"])) + fwhm_guess = options["ap_set_psf"] + elif "ap_guess_psf" in options: + logging.info( + "%s: PSF initialized by user: %.4e pix" % (options["ap_name"], options["ap_guess_psf"]) + ) + fwhm_guess = options["ap_guess_psf"] + else: + logging.info("%s: PSF assumed to be: %.4e pix" % (options["ap_name"], 4.0)) + fwhm_guess = 4.0 + return IMG, {"psf fwhm": fwhm_guess, "auxfile psf": "psf fwhm: %.3f pix" % fwhm_guess} + + def PSF_StarFind(IMG, results, options): """PSF routine which identifies stars and averages the FWHM. @@ -164,14 +183,11 @@ def PSF_StarFind(IMG, results, options): """ if "ap_set_psf" in options: - logging.info( - "%s: PSF set by user: %.4e" % (options["ap_name"], options["ap_set_psf"]) - ) + logging.info("%s: PSF set by user: %.4e" % (options["ap_name"], options["ap_set_psf"])) return IMG, {"psf fwhm": options["ap_set_psf"]} elif "ap_guess_psf" in options: logging.info( - "%s: PSF initialized by user: %.4e" - % (options["ap_name"], options["ap_guess_psf"]) + "%s: PSF initialized by user: %.4e" % (options["ap_name"], options["ap_guess_psf"]) ) fwhm_guess = options["ap_guess_psf"] else: @@ -185,7 +201,7 @@ def PSF_StarFind(IMG, results, options): int(IMG.shape[0] / 5.0) : int(4.0 * IMG.shape[0] / 5.0), int(IMG.shape[1] / 5.0) : int(4.0 * IMG.shape[1] / 5.0), ] = True - + stars = StarFind( IMG - results["background"], fwhm_guess, @@ -217,8 +233,7 @@ def PSF_StarFind(IMG, results, options): ) logging.info( - "%s: found psf: %f with deformity clip of: %f" - % (options["ap_name"], psf, def_clip) + "%s: found psf: %f with deformity clip of: %f" % (options["ap_name"], psf, def_clip) ) return IMG, {"psf fwhm": psf, "auxfile psf": "psf fwhm: %.3f pix" % psf} @@ -266,14 +281,11 @@ def PSF_Image(IMG, results, options): """ if "ap_set_psf" in options: - logging.info( - "%s: PSF set by user: %.4e" % (options["ap_name"], options["ap_set_psf"]) - ) + logging.info("%s: PSF set by user: %.4e" % (options["ap_name"], options["ap_set_psf"])) return IMG, {"psf fwhm": options["ap_set_psf"]} elif "ap_guess_psf" in options: logging.info( - "%s: PSF initialized by user: %.4e" - % (options["ap_name"], options["ap_guess_psf"]) + "%s: PSF initialized by user: %.4e" % (options["ap_name"], options["ap_guess_psf"]) ) fwhm_guess = options["ap_guess_psf"] else: @@ -331,9 +343,9 @@ def PSF_Image(IMG, results, options): or (dat.shape[1] - stars["y"][i]) < psf_size // 2 ): continue - flux = interpolate_Lanczos( - dat, XX + stars["x"][i], YY + stars["y"][i], 10 - ).reshape((1, psf_size, psf_size)) + flux = interpolate_Lanczos(dat, XX + stars["x"][i], YY + stars["y"][i], 10).reshape( + (1, psf_size, psf_size) + ) flux /= np.sum(flux) psf_img = flux if psf_img is None else np.concatenate((psf_img, flux)) @@ -448,9 +460,7 @@ def PSF_deconvolve(IMG, results, options): np.array(range(psf_size)) - psf_size // 2, ) psf_std = results["psf fwhm"] / np.sqrt(8 * np.log(2)) - psf_img = np.exp(-(XX ** 2 + YY ** 2) / (2 * psf_std ** 2)) / np.sqrt( - 2 * np.pi * psf_std ** 2 - ) + psf_img = np.exp(-(XX**2 + YY**2) / (2 * psf_std**2)) / np.sqrt(2 * np.pi * psf_std**2) psf_img /= np.sum(psf_img) if np.abs(np.sum(psf_img) - 1) > 1e-7: @@ -460,16 +470,18 @@ def PSF_deconvolve(IMG, results, options): dat_deconv = restoration.richardson_lucy( (IMG - dmin) / (dmax - dmin) - 0.5, psf_img, - iterations=options["ap_psf_deconvolution_iterations"] - if "ap_psf_deconvolution_iterations" in options - else 50, + iterations=( + options["ap_psf_deconvolution_iterations"] + if "ap_psf_deconvolution_iterations" in options + else 50 + ), ) dat_deconv = (dat_deconv + 0.5) * (dmax - dmin) + dmin if "ap_psf_deconvolve_save" in options and options["ap_psf_deconvolve_save"]: header = fits.Header() hdul = fits.HDUList([fits.PrimaryHDU(header=header), fits.ImageHDU(dat_deconv)]) - + hdul.writeto( os.path.join( options["ap_saveto"] if "ap_saveto" in options else "", @@ -477,6 +489,5 @@ def PSF_deconvolve(IMG, results, options): ), overwrite=True, ) - - + return dat_deconv, {} From b9d9cd5fb2c9bb98f61399f27972ac3c23fb9a95 Mon Sep 17 00:00:00 2001 From: Connor Stone Date: Tue, 5 Mar 2024 21:59:33 -0500 Subject: [PATCH 2/3] impose fixed radius in fitting too, update parameters list --- autoprof/pipeline_steps/Isophote_Fit.py | 151 +++++++----------- .../pipeline_steps/Isophote_Initialize.py | 3 + docs/defaultpipeline.rst | 4 +- docs/parameters.rst | 77 ++++++++- 4 files changed, 134 insertions(+), 101 deletions(-) diff --git a/autoprof/pipeline_steps/Isophote_Fit.py b/autoprof/pipeline_steps/Isophote_Fit.py index 69e7f4d9..7af963d9 100644 --- a/autoprof/pipeline_steps/Isophote_Fit.py +++ b/autoprof/pipeline_steps/Isophote_Fit.py @@ -33,7 +33,14 @@ ) from ..autoprofutils.Diagnostic_Plots import Plot_Isophote_Fit -__all__ = ("Photutils_Fit", "Isophote_Fit_FixedPhase", "Isophote_Fit_FFT_Robust", "Isophote_Fit_Forced", "Isophote_Fit_FFT_mean") +__all__ = ( + "Photutils_Fit", + "Isophote_Fit_FixedPhase", + "Isophote_Fit_FFT_Robust", + "Isophote_Fit_Forced", + "Isophote_Fit_FFT_mean", +) + def Photutils_Fit(IMG, results, options): """Photutils elliptical isophote wrapper. @@ -228,7 +235,7 @@ def _pa_smooth(R, PA, deg): def _FFT_Robust_loss( - dat, R, PARAMS, i, C, noise, mask=None, reg_scale=1.0, robust_clip=0.15, fit_coefs=None, name="" + dat, R, PARAMS, i, C, noise, mask=None, reg_scale=1.0, robust_clip=0.15, fit_coefs=None, name="" ): isovals = _iso_extract( @@ -242,11 +249,11 @@ def _FFT_Robust_loss( ) try: - coefs = fft(np.clip(isovals, a_max=np.quantile(isovals, 1. - robust_clip), a_min=None)) + coefs = fft(np.clip(isovals, a_max=np.quantile(isovals, 1.0 - robust_clip), a_min=None)) except: coefs = np.zeros(100) isovals = np.zeros(100) - + if fit_coefs is None: f2_loss = np.abs(coefs[2]) / ( len(isovals) * (max(0, np.median(isovals)) + noise / np.sqrt(len(isovals))) @@ -265,14 +272,10 @@ def _FFT_Robust_loss( reg_loss += abs( (PARAMS[i]["ellip"] - PARAMS[i + 1]["ellip"]) / (1 - PARAMS[i + 1]["ellip"]) ) - reg_loss += abs( - Angle_TwoAngles_sin(PARAMS[i]["pa"], PARAMS[i + 1]["pa"]) / (0.2) - ) + reg_loss += abs(Angle_TwoAngles_sin(PARAMS[i]["pa"], PARAMS[i + 1]["pa"]) / (0.2)) if not PARAMS[i]["m"] is None: for m in range(len(PARAMS[i]["m"])): - reg_loss += fmode_scale * abs( - (PARAMS[i]["Am"][m] - PARAMS[i + 1]["Am"][m]) / 0.2 - ) + reg_loss += fmode_scale * abs((PARAMS[i]["Am"][m] - PARAMS[i + 1]["Am"][m]) / 0.2) reg_loss += fmode_scale * abs( Angle_TwoAngles_cos( PARAMS[i]["m"][m] * PARAMS[i]["Phim"][m], @@ -286,14 +289,10 @@ def _FFT_Robust_loss( reg_loss += abs( (PARAMS[i]["ellip"] - PARAMS[i - 1]["ellip"]) / (1 - PARAMS[i - 1]["ellip"]) ) - reg_loss += abs( - Angle_TwoAngles_sin(PARAMS[i]["pa"], PARAMS[i - 1]["pa"]) / (0.2) - ) + reg_loss += abs(Angle_TwoAngles_sin(PARAMS[i]["pa"], PARAMS[i - 1]["pa"]) / (0.2)) if not PARAMS[i]["m"] is None: for m in range(len(PARAMS[i]["m"])): - reg_loss += fmode_scale * abs( - (PARAMS[i]["Am"][m] - PARAMS[i - 1]["Am"][m]) / 0.2 - ) + reg_loss += fmode_scale * abs((PARAMS[i]["Am"][m] - PARAMS[i - 1]["Am"][m]) / 0.2) reg_loss += fmode_scale * abs( Angle_TwoAngles_cos( PARAMS[i]["m"][m] * PARAMS[i]["Phim"][m], @@ -308,7 +307,7 @@ def _FFT_Robust_loss( def _FFT_Robust_Errors( - dat, R, PARAMS, C, noise, mask=None, reg_scale=1.0, robust_clip=0.15, fit_coefs=None, name="" + dat, R, PARAMS, C, noise, mask=None, reg_scale=1.0, robust_clip=0.15, fit_coefs=None, name="" ): PA_err = np.zeros(len(R)) @@ -339,7 +338,7 @@ def _FFT_Robust_Errors( noise, mask=mask, reg_scale=reg_scale, - robust_clip = robust_clip, + robust_clip=robust_clip, fit_coefs=fit_coefs, name=name, ), @@ -512,6 +511,12 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): if not np.any(mask): mask = None + if "ap_isoinit_R_set" in options: + minR = options["ap_isoinit_R_set"] + elif "init R" in results: + minR = results["init R"] * 0.8 + else: + minR = 0.0 # Determine sampling radii ###################################################################### shrink = 0 @@ -527,7 +532,8 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): mask=mask, ) if ( - np.median(isovals) + sample_radii[-1] > minR + and np.median(isovals) < (options["ap_fit_limit"] if "ap_fit_limit" in options else 2) * results["background noise"] ): @@ -547,23 +553,13 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): # Fit isophotes ###################################################################### perturb_scale = 0.03 - regularize_scale = ( - options["ap_regularize_scale"] if "ap_regularize_scale" in options else 1.0 - ) - robust_clip = ( - options["ap_isofit_robustclip"] if "ap_isofit_robustclip" in options else 0.15 - ) + regularize_scale = options["ap_regularize_scale"] if "ap_regularize_scale" in options else 1.0 + robust_clip = options["ap_isofit_robustclip"] if "ap_isofit_robustclip" in options else 0.15 N_perturb = 5 - fit_coefs = ( - options["ap_isofit_losscoefs"] if "ap_isofit_losscoefs" in options else None - ) - fit_params = ( - options["ap_isofit_fitcoefs"] if "ap_isofit_fitcoefs" in options else None - ) + fit_coefs = options["ap_isofit_losscoefs"] if "ap_isofit_losscoefs" in options else None + fit_params = options["ap_isofit_fitcoefs"] if "ap_isofit_fitcoefs" in options else None fit_superellipse = ( - options["ap_isofit_superellipse"] - if "ap_isofit_superellipse" in options - else False + options["ap_isofit_superellipse"] if "ap_isofit_superellipse" in options else False ) parameters = list( { @@ -580,17 +576,11 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): count = 0 iterlimitmax = ( - options["ap_isofit_iterlimitmax"] - if "ap_isofit_iterlimitmax" in options - else 1000 - ) - iterlimitmin = ( - options["ap_isofit_iterlimitmin"] if "ap_isofit_iterlimitmin" in options else 0 + options["ap_isofit_iterlimitmax"] if "ap_isofit_iterlimitmax" in options else 1000 ) + iterlimitmin = options["ap_isofit_iterlimitmin"] if "ap_isofit_iterlimitmin" in options else 0 iterstopnochange = ( - options["ap_isofit_iterstopnochange"] - if "ap_isofit_iterstopnochange" in options - else 3 + options["ap_isofit_iterstopnochange"] if "ap_isofit_iterstopnochange" in options else 3 ) count_nochange = 0 use_center = copy(results["center"]) @@ -618,7 +608,7 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): results["background noise"], mask=mask, reg_scale=regularize_scale if count > 4 else 0, - robust_clip = robust_clip, + robust_clip=robust_clip, fit_coefs=fit_coefs, name=options["ap_name"], ) @@ -644,23 +634,17 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): (count % param_cycle) - base_params ] += np.random.normal(loc=0, scale=perturb_scale) elif count % param_cycle < (base_params + 2 * len(parameters[i]["m"])): - phim_index = ( - (count % param_cycle) - base_params - len(parameters[i]["m"]) - ) + phim_index = (count % param_cycle) - base_params - len(parameters[i]["m"]) perturbations[-1][i]["Phim"][phim_index] = ( perturbations[-1][i]["Phim"][phim_index] + np.random.normal( loc=0, - scale=2 - * np.pi - * perturb_scale - / parameters[i]["m"][phim_index], + scale=2 * np.pi * perturb_scale / parameters[i]["m"][phim_index], ) ) % (2 * np.pi / parameters[i]["m"][phim_index]) else: raise Exception( - "Unrecognized optimization parameter id: %i" - % (count % param_cycle) + "Unrecognized optimization parameter id: %i" % (count % param_cycle) ) perturbations[-1][i]["loss"] = _FFT_Robust_loss( dat, @@ -671,7 +655,7 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): results["background noise"], mask=mask, reg_scale=regularize_scale if count > 4 else 0, - robust_clip = robust_clip, + robust_clip=robust_clip, fit_coefs=fit_coefs, name=options["ap_name"], ) @@ -687,14 +671,11 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): count_nochange < (iterstopnochange * (len(sample_radii) - 1)) or count < iterlimitmin ): - if param_cycle > 2 or ( - parameters[i]["m"] is None and not fit_superellipse - ): + if param_cycle > 2 or (parameters[i]["m"] is None and not fit_superellipse): break elif parameters[i]["m"] is None and fit_superellipse: logging.info( - "%s: Started C fitting at iteration %i" - % (options["ap_name"], count) + "%s: Started C fitting at iteration %i" % (options["ap_name"], count) ) param_cycle = 3 iterstopnochange = max(iterstopnochange, param_cycle) @@ -704,13 +685,11 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): fit_coefs = (2, 4) else: logging.info( - "%s: Started Fmode fitting at iteration %i" - % (options["ap_name"], count) + "%s: Started Fmode fitting at iteration %i" % (options["ap_name"], count) ) if fit_superellipse: logging.info( - "%s: Started C fitting at iteration %i" - % (options["ap_name"], count) + "%s: Started C fitting at iteration %i" % (options["ap_name"], count) ) param_cycle = base_params + 2 * len(parameters[i]["m"]) iterstopnochange = max(iterstopnochange, param_cycle) @@ -765,14 +744,11 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): ) % (2 * np.pi) if not ( - count_nochange < (iterstopnochange * (len(sample_radii) - 1)) - or count < iterlimitmin + count_nochange < (iterstopnochange * (len(sample_radii) - 1)) or count < iterlimitmin ): break - logging.info( - "%s: Completed isohpote fit in %i itterations" % (options["ap_name"], count) - ) + logging.info("%s: Completed isohpote fit in %i itterations" % (options["ap_name"], count)) # Compute errors ###################################################################### ellip_err, pa_err = _FFT_Robust_Errors( @@ -783,7 +759,7 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): results["background noise"], mask=mask, reg_scale=regularize_scale, - robust_clip = robust_clip, + robust_clip=robust_clip, fit_coefs=fit_coefs, name=options["ap_name"], ) @@ -796,9 +772,7 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): Plot_Isophote_Fit(dat, sample_radii, parameters, results, options) res = { - "fit ellip": np.array( - list(parameters[i]["ellip"] for i in range(len(parameters))) - ), + "fit ellip": np.array(list(parameters[i]["ellip"] for i in range(len(parameters)))), "fit pa": np.array(list(parameters[i]["pa"] for i in range(len(parameters)))), "fit R": sample_radii, "fit ellip_err": ellip_err, @@ -821,13 +795,7 @@ def Isophote_Fit_FFT_Robust(IMG, results, options): } ) if fit_superellipse: - res.update( - { - "fit C": np.array( - list(parameters[i]["C"] for i in range(len(parameters))) - ) - } - ) + res.update({"fit C": np.array(list(parameters[i]["C"] for i in range(len(parameters))))}) return IMG, res @@ -884,7 +852,7 @@ def Isophote_Fit_Forced(IMG, results, options): if "ap_doplot" in options and options["ap_doplot"]: parameters = [] for i in range(len(force["R"])): - parameters.append({'ellip': force["ellip"][i], 'pa': force["pa"][i], 'C': None}) + parameters.append({"ellip": force["ellip"][i], "pa": force["pa"][i], "C": None}) Plot_Isophote_Fit( IMG - results["background"], np.array(force["R"]), @@ -1036,9 +1004,7 @@ def Isophote_Fit_FFT_mean(IMG, results, options): # Fit isophotes ###################################################################### perturb_scale = np.array([0.03, 0.06]) - regularize_scale = ( - options["ap_regularize_scale"] if "ap_regularize_scale" in options else 1.0 - ) + regularize_scale = options["ap_regularize_scale"] if "ap_regularize_scale" in options else 1.0 N_perturb = 5 count = 0 @@ -1077,8 +1043,7 @@ def Isophote_Fit_FFT_mean(IMG, results, options): ) if count % 3 in [1, 2]: perturbations[-1]["pa"][i] = ( - perturbations[-1]["pa"][i] - + np.random.normal(loc=0, scale=perturb_scale[1]) + perturbations[-1]["pa"][i] + np.random.normal(loc=0, scale=perturb_scale[1]) ) % np.pi perturbations[-1]["loss"] = _FFT_mean_loss( dat, @@ -1101,9 +1066,7 @@ def Isophote_Fit_FFT_mean(IMG, results, options): else: count_nochange += 1 - logging.info( - "%s: Completed isohpote fit in %i itterations" % (options["ap_name"], count) - ) + logging.info("%s: Completed isohpote fit in %i itterations" % (options["ap_name"], count)) # detect collapsed center ###################################################################### for i in range(5): @@ -1139,10 +1102,10 @@ def Isophote_Fit_FFT_mean(IMG, results, options): for i in range(len(sample_radii)): plt.gca().add_patch( Ellipse( - xy = (use_center["x"] - ranges[0][0], use_center["y"] - ranges[1][0]), - width = 2 * sample_radii[i], - height = 2 * sample_radii[i] * (1.0 - ellip[i]), - angle = pa[i] * 180 / np.pi, + xy=(use_center["x"] - ranges[0][0], use_center["y"] - ranges[1][0]), + width=2 * sample_radii[i], + height=2 * sample_radii[i] * (1.0 - ellip[i]), + angle=pa[i] * 180 / np.pi, fill=False, linewidth=((i + 1) / len(sample_radii)) ** 2, color="r", @@ -1198,9 +1161,7 @@ def Isophote_Fit_FFT_mean(IMG, results, options): ellip_err[i] = np.sqrt( np.sum((ellip[i - 2 : i + 2] - smooth_ellip[i - 2 : i + 2]) ** 2) / 4 ) - pa_err[i] = np.sqrt( - np.sum((pa[i - 2 : i + 2] - smooth_pa[i - 2 : i + 2]) ** 2) / 4 - ) + pa_err[i] = np.sqrt(np.sum((pa[i - 2 : i + 2] - smooth_pa[i - 2 : i + 2]) ** 2) / 4) res = { "fit ellip": ellip, diff --git a/autoprof/pipeline_steps/Isophote_Initialize.py b/autoprof/pipeline_steps/Isophote_Initialize.py index aea4547e..2f2e5256 100644 --- a/autoprof/pipeline_steps/Isophote_Initialize.py +++ b/autoprof/pipeline_steps/Isophote_Initialize.py @@ -161,6 +161,9 @@ def Isophote_Initialize(IMG, results, options): ap_isoinit_ellip_set : float, default None User set initial ellipticity (1 - b/a), will override the calculation. + ap_isoinit_R_set : float, default None + User set initial semi-major axis length, will override the calculation. + Notes ---------- :References: diff --git a/docs/defaultpipeline.rst b/docs/defaultpipeline.rst index 4af30676..41979b3c 100644 --- a/docs/defaultpipeline.rst +++ b/docs/defaultpipeline.rst @@ -35,7 +35,7 @@ a similarly size or larger object. Put plainly, the default AutoProf pipeline is as follows: 1. Background: :func:`~autoprof.pipeline_steps.Background.Background_Mode` -#. PSF: :func:`~autoprof.pipeline_steps.PSF.PSF_StarFind` +#. PSF: :func:`~autoprof.pipeline_steps.PSF.PSF_Assumed` #. Center: :func:`~autoprof.pipeline_steps.Center.Center_HillClimb` #. Initial Isophote: :func:`~autoprof.pipeline_steps.Isophote_Initialize.Isophote_Initialize` #. Fit Isophotes: :func:`~autoprof.pipeline_steps.Isophote_Fit.Isophote_Fit_FFT_Robust` @@ -78,7 +78,7 @@ onto another image. The default forced photometry pipeline works as follows: 1. Background: :func:`~autoprof.pipeline_steps.Background.Background_Mode` -#. PSF: :func:`~autoprof.pipeline_steps.PSF.PSF_StarFind` +#. PSF: :func:`~autoprof.pipeline_steps.PSF.PSF_Assumed` #. Center: :func:`~autoprof.pipeline_steps.Center.Center_Forced` #. Initial Isophote: :func:`~autoprof.pipeline_steps.Isophote_Initialize.Isophote_Init_Forced` #. Fit Isophotes: :func:`~autoprof.pipeline_steps.Isophote_Fit.Isophote_Fit_Forced` diff --git a/docs/parameters.rst b/docs/parameters.rst index 2ca1a6b0..4e46a5aa 100644 --- a/docs/parameters.rst +++ b/docs/parameters.rst @@ -157,14 +157,28 @@ ap_fit_limit (float, default 2) noise level out to which to extend the fit in units of pixel background noise level. Default is 2, smaller values will end fitting further out in the galaxy image. +ap_fluxunits (str, default "mag") +---------------------------------------------------------------------- + +**Referencing Methods** + +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Forced` +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract` +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Photutils` + +**Description** + +units for outputted photometry. Can either be "mag" for log +units, or "intensity" for linear units. + ap_forcing_profile (string, default None) ---------------------------------------------------------------------- **Referencing Methods** - :func:`~autoprof.pipeline_steps.Isophote_Initialize.Isophote_Init_Forced` -- :func:`~autoprof.pipeline_steps.Center.Center_Forced` - :func:`~autoprof.pipeline_steps.Isophote_Fit.Isophote_Fit_Forced` +- :func:`~autoprof.pipeline_steps.Center.Center_Forced` **Description** @@ -275,9 +289,9 @@ ap_isoaverage_method (string, default 'median') **Referencing Methods** +- :func:`~autoprof.pipeline_steps.Slice_Profiles.Slice_Profile` - :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Forced` - :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract` -- :func:`~autoprof.pipeline_steps.Slice_Profiles.Slice_Profile` - :func:`~autoprof.pipeline_steps.Axial_Profiles.Axial_Profiles` **Description** @@ -525,6 +539,17 @@ from a rectangle to an ellipse to a four-armed-star like shape is C. A value of C = 2 represents an ellipse and is the starting point of the optimization. +ap_isoinit_R_set (float, default None) +---------------------------------------------------------------------- + +**Referencing Methods** + +- :func:`~autoprof.pipeline_steps.Isophote_Initialize.Isophote_Initialize` + +**Description** + +User set initial semi-major axis length, will override the calculation. + ap_isoinit_ellip_set (float, default None) ---------------------------------------------------------------------- @@ -569,6 +594,50 @@ ap_name (string, default None) Name of the current galaxy, used for making filenames. +ap_plot_sbprof_set_errscale (float, default None) +---------------------------------------------------------------------- + +**Referencing Methods** + +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Forced` +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract` +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Photutils` + +**Description** + +Float value by which to scale errorbars on the SB profile +this makes them more visible in cases where the statistical +errors are very small. + +ap_plot_sbprof_xlim (tuple, default None) +---------------------------------------------------------------------- + +**Referencing Methods** + +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Forced` +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract` +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Photutils` + +**Description** + +Tuple with axes limits for the x-axis in the SB profile +diagnostic plot. + +ap_plot_sbprof_ylim (tuple, default None) +---------------------------------------------------------------------- + +**Referencing Methods** + +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Forced` +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract` +- :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Photutils` + +**Description** + +Tuple with axes limits for the y-axis in the SB profile +diagnostic plot. Be careful when using intensity units +since this will change the ideal axis limits. + ap_psf_deconvolution_iterations (int, default 50) ---------------------------------------------------------------------- @@ -903,11 +972,11 @@ ap_zeropoint (float, default 22.5) **Referencing Methods** +- :func:`~autoprof.pipeline_steps.Slice_Profiles.Slice_Profile` +- :func:`~autoprof.pipeline_steps.Ellipse_Model.EllipseModel` - :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Forced` - :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract` - :func:`~autoprof.pipeline_steps.Isophote_Extract.Isophote_Extract_Photutils` -- :func:`~autoprof.pipeline_steps.Ellipse_Model.EllipseModel` -- :func:`~autoprof.pipeline_steps.Slice_Profiles.Slice_Profile` - :func:`~autoprof.pipeline_steps.Axial_Profiles.Axial_Profiles` **Description** From 738d763cacee08dea310f1df58e6a813b1663cde Mon Sep 17 00:00:00 2001 From: Connor Stone Date: Tue, 5 Mar 2024 22:00:44 -0500 Subject: [PATCH 3/3] increment version --- autoprof/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/autoprof/__init__.py b/autoprof/__init__.py index d47daa7b..e85bdb88 100644 --- a/autoprof/__init__.py +++ b/autoprof/__init__.py @@ -3,10 +3,11 @@ from . import autoprofutils, Pipeline, pipeline_steps -__version__ = "1.0.2" +__version__ = "1.1.0" __author__ = "Connor Stone" __email__ = "connorstone628@gmail.com" + def run_from_terminal(): assert ( @@ -26,4 +27,3 @@ def run_from_terminal(): PIPELINE = Pipeline.Isophote_Pipeline(loggername=logfile) PIPELINE.Process_ConfigFile(config_file) -