From 21207d61bbf07fd8733b1d1df4b0d1de90bfd31d Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Thu, 4 Apr 2024 15:33:46 -0700 Subject: [PATCH 01/11] Fetch inputs from EHT repo --- data_validation/scripts/unpack_data.sh | 3 +++ 1 file changed, 3 insertions(+) diff --git a/data_validation/scripts/unpack_data.sh b/data_validation/scripts/unpack_data.sh index cdb16be..3f8ec8b 100755 --- a/data_validation/scripts/unpack_data.sh +++ b/data_validation/scripts/unpack_data.sh @@ -1,5 +1,8 @@ #!/bin/bash mkdir -p images mkdir -p data +wget https://github.com/eventhorizontelescope/2019-D01-01/raw/master/EHTC_FirstM87Results_Apr2019_uvfits.tgz +wget https://github.com/eventhorizontelescope/2019-D01-01/raw/master/EHTC_FirstM87Results_Apr2019_csv.tgz +tar -xvzf EHTC_FirstM87Results_Apr2019_uvfits.tgz -C data --strip-components=1 tar -xvzf EHTC_FirstM87Results_Apr2019_csv.tgz -C data --strip-components=1 mkdir -p data/csv/converted From 2acca27e59064be3cf6775c97e727ecdb5722b27 Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Thu, 4 Apr 2024 15:34:56 -0700 Subject: [PATCH 02/11] Create output dir before commands and fix path to data dir --- src/difmap/run-postprocessing.sh | 6 ++++-- src/difmap/run.sh | 4 +++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/difmap/run-postprocessing.sh b/src/difmap/run-postprocessing.sh index a3ccaf9..9573036 100755 --- a/src/difmap/run-postprocessing.sh +++ b/src/difmap/run-postprocessing.sh @@ -1,14 +1,16 @@ #!/usr/bin/env bash + +mkdir -p difmap-pdfs difmap-imgsums cp difmap-output/*.fits . -for f in *.fits; do +for f in *.fits; do python difmap-postprocessing.py -i $f --all done for d in 095 096 100 101; do python difmap-imgsum.py \ -i SR1_M87_2017_${d}_lo_hops_netcal_StokesI.CircMask_r30_x-0.002_y0.022.RT-10.CF0.5.ALMA0.1.UVW2_-1.noresiduals.fits \ - -o ../data/uvfits/SR1_M87_2017_${d}_lo_hops_netcal_StokesI.uvfits + -o ../../data_validation/data/uvfits/SR1_M87_2017_${d}_lo_hops_netcal_StokesI.uvfits done rm *.fits diff --git a/src/difmap/run.sh b/src/difmap/run.sh index 2ac5d58..5e217bc 100755 --- a/src/difmap/run.sh +++ b/src/difmap/run.sh @@ -15,7 +15,9 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -cp ../data/uvfits/*.uvfits . + +mkdir -p difmap-output +cp ../../data_validation/data/uvfits/*.uvfits . for f in *.uvfits; do ./difmap.sh $f done From f70a5aa1ede3edce4fc3e6b596945eefe534a2ef Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Thu, 4 Apr 2024 15:35:55 -0700 Subject: [PATCH 03/11] Create output dir before commands and fix path to data dir --- src/eht-imaging/run-pipeline.sh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/eht-imaging/run-pipeline.sh b/src/eht-imaging/run-pipeline.sh index 4e53fcc..78704d8 100755 --- a/src/eht-imaging/run-pipeline.sh +++ b/src/eht-imaging/run-pipeline.sh @@ -9,12 +9,13 @@ # To save .pdf of image summary statistics, uncomment --imsum (This doesn't work currently) # +mkdir -p ../output + echo "=============== Beginning EHT-Imaging Pipeline Execution =========================" -cd scripts for d in 095 096 100 101; do python eht-imaging_pipeline.py \ - -i ../data/uvfits/SR1_M87_2017_${d}_lo_hops_netcal_StokesI.uvfits \ - -i2 ../data/uvfits/SR1_M87_2017_${d}_hi_hops_netcal_StokesI.uvfits \ + -i ../../data_validation/data/uvfits/SR1_M87_2017_${d}_lo_hops_netcal_StokesI.uvfits \ + -i2 ../../data_validation/data/uvfits/SR1_M87_2017_${d}_hi_hops_netcal_StokesI.uvfits \ -o ../output/SR1_M87_2017_${d}.fits \ --savepdf \ --imgsum @@ -24,4 +25,3 @@ echo " echo "================ Beginning EHT-Imaging Post-processing ===========================" bash run-postprocessing.sh echo "================ Completed EHT-Imaging Post-processing ===========================" -cd ~ \ No newline at end of file From 1c4d0da1fec30455dcffbd97927073ac039eb666 Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Thu, 4 Apr 2024 15:36:32 -0700 Subject: [PATCH 04/11] Add missing python file --- src/eht-imaging/eht-imaging_pipeline.py | 370 ++++++++++++++++++++++++ 1 file changed, 370 insertions(+) create mode 100644 src/eht-imaging/eht-imaging_pipeline.py diff --git a/src/eht-imaging/eht-imaging_pipeline.py b/src/eht-imaging/eht-imaging_pipeline.py new file mode 100644 index 0000000..b7b82db --- /dev/null +++ b/src/eht-imaging/eht-imaging_pipeline.py @@ -0,0 +1,370 @@ +"""eht-imaging M87 Stokes I Imaging Pipeline for EHT observations in April 2017 + +Authors: The Event Horizon Telescope Collaboration et al. +Date: April 10, 2019 +Primary Reference: The Event Horizon Telescope Collaboration, et al. 2019d, ApJL, 875, L4 (M87 Paper IV) +Data Product Code: 2019-D01-02 + +Brief Description: +The pipeline reconstructs an image from uvfits files simultaneously +released in the EHT website (data release ID: 2019-D01-01) using a +python-interfaced imaging package eht-imaging (Chael et al. 2016,2018). + +To run the pipeline, specify the input uvfits data file. Multiple bands +can be specified separately. Additional flags control the output, which +is only the reconstructed image as a FITS image by default. + +Example call: +python eht-imaging_pipeline.py -i ../SR1_M87_2017_101_lo_hops_netcal_StokesI.uvfits -i2 ../SR1_M87_2017_101_hi_hops_netcal_StokesI.uvfits --savepdf + +For additional details, please read the help document associated in the +imaging script: "python eht-imaging_pipeline.py --help". + +Additional References: + - EHT Collaboration Data Portal Website: + https://eventhorizontelescope.org/for-astronomers/data + - The Event Horizon Telescope Collaboration, et al. 2019a, ApJL, 875, L1 (M87 Paper I) + - The Event Horizon Telescope Collaboration, et al. 2019b, ApJL, 875, L2 (M87 Paper II) + - The Event Horizon Telescope Collaboration, et al. 2019c, ApJL, 875, L3 (M87 Paper III) + - The Event Horizon Telescope Collaboration, et al. 2019d, ApJL, 875, L4 (M87 Paper IV) + - The Event Horizon Telescope Collaboration, et al. 2019e, ApJL, 875, L5 (M87 Paper V) + - The Event Horizon Telescope Collaboration, et al. 2019f, ApJL, 875, L6 (M87 Paper VI) + - Chael, A., Johnson, M., Narayan, R., et al. 2016, ApJ, 829, 11C + - Chael, A., Johnson, M., Bouman, K., et al. 2018, ApJ, 857, 23C + - Chael, A., Bouman, K., Johnson, M. et al., Zenodo (eht-imaging version 1.1.0) + - eht-imaging: https://github.com/achael/eht-imaging +""" + +#------------------------------------------------------------------------------- +# Author Information +#------------------------------------------------------------------------------- +__author__ = "The Event Horizon Telescope Collaboration et al." +__copyright__ = "Copyright 2019, the Event Horizon Telescope Collaboration et al." +__license__ = "GPL version 3" +__version__ = "1.0" +__date__ = "April 10 2019" + +#------------------------------------------------------------------------------- +# Modules +#------------------------------------------------------------------------------- +import matplotlib +matplotlib.use('Agg') + +import os +import argparse +import ehtim as eh +import numpy as np + +#------------------------------------------------------------------------------- +# Load command-line arguments +#------------------------------------------------------------------------------- +parser = argparse.ArgumentParser(description="Fiducial eht-imaging script for M87") +parser.add_argument('-i', '--infile', default="obs.uvfits",help="input UVFITS file") +parser.add_argument('-i2', '--infile2', default="", help="optional 2nd input file (different band) for imaging") +parser.add_argument('-o', '--outfile', default='out.fits', help='output FITS image') +parser.add_argument('--savepdf', default=False, help='saves image pdf (True or False)?',action='store_true') +parser.add_argument('--imgsum', default=False, help='generate image summary pdf',action='store_true') +args = parser.parse_args() + +#------------------------------------------------------------------------------- +# Fiducial imaging parameters obtained from the eht-imaging parameter survey +#------------------------------------------------------------------------------- +zbl = 0.60 # Total compact flux density (Jy) +prior_fwhm = 40.0*eh.RADPERUAS # Gaussian prior FWHM (radians) +sys_noise = 0.02 # fractional systematic noise + # added to complex visibilities + +# constant regularization weights +reg_term = {'simple' : 100, # Maximum-Entropy + 'tv' : 1.0, # Total Variation + 'tv2' : 1.0, # Total Squared Variation + 'l1' : 0.0, # L1 sparsity prior + 'flux' : 1e4} # compact flux constraint + +# initial data weights - these are updated throughout the imaging pipeline +data_term = {'amp' : 0.2, # visibility amplitudes + 'cphase' : 1.0, # closure phases + 'logcamp': 1.0} # log closure amplitudes + + +#------------------------------------------------------------------------------- +# Fixed imaging parameters +#------------------------------------------------------------------------------- +obsfile = args.infile # Pre-processed observation file +ttype = 'nfft' # Type of Fourier transform ('direct', 'nfft', or 'fast') +npix = 64 # Number of pixels across the reconstructed image +fov = 128*eh.RADPERUAS # Field of view of the reconstructed image +maxit = 100 # Maximum number of convergence iterations for imaging +stop = 1e-4 # Imager stopping criterion +gain_tol = [0.02,0.2] # Asymmetric gain tolerance for self-cal; we expect larger values + # for unaccounted sensitivity loss + # than for unaccounted sensitivity improvement +uv_zblcut = 0.1e9 # uv-distance that separates the inter-site "zero"-baselines + # from intra-site baselines +reverse_taper_uas = 5.0 # Finest resolution of reconstructed features + +# Specify the SEFD error budget +# (reported in First M87 Event Horizon Telescope Results III: Data Processing and Calibration) +SEFD_error_budget = {'AA':0.10, + 'AP':0.11, + 'AZ':0.07, + 'LM':0.22, + 'PV':0.10, + 'SM':0.15, + 'JC':0.14, + 'SP':0.07} + +# Add systematic noise tolerance for amplitude a-priori calibration errors +# Start with the SEFD noise (but need sqrt) +# then rescale to ensure that final results respect the stated error budget +systematic_noise = SEFD_error_budget.copy() +for key in systematic_noise.keys(): + systematic_noise[key] = ((1.0+systematic_noise[key])**0.5 - 1.0) * 0.25 + +# Extra noise added for the LMT, which has much more variability than the a-priori error budget +systematic_noise['LM'] += 0.15 + +#------------------------------------------------------------------------------- +# Define helper functions +#------------------------------------------------------------------------------- + +# Rescale short baselines to excise contributions from extended flux. +# setting zbl < zbl_tot assumes there is an extended constant flux component of zbl_tot-zbl Jy +def rescale_zerobaseline(obs, totflux, orig_totflux, uv_max): + multiplier = zbl / zbl_tot + for j in range(len(obs.data)): + if (obs.data['u'][j]**2 + obs.data['v'][j]**2)**0.5 >= uv_max: continue + for field in ['vis','qvis','uvis','vvis','sigma','qsigma','usigma','vsigma']: + obs.data[field][j] *= multiplier + +# repeat imaging with blurring to assure good convergence +def converge(major=3, blur_frac=1.0): + for repeat in range(major): + init = imgr.out_last().blur_circ(blur_frac*res) + imgr.init_next = init + imgr.make_image_I(show_updates=False) + +#------------------------------------------------------------------------------- +# Prepare the data +#------------------------------------------------------------------------------- + +# Load a single uvfits file +if args.infile2 == '': + # load the uvfits file + obs = eh.obsdata.load_uvfits(obsfile) + + # scan-average the data + # identify the scans (times of continous observation) in the data + obs.add_scans() + + # coherently average the scans, which can be averaged due to ad-hoc phasing + obs = obs.avg_coherent(0.,scan_avg=True) + +# If two uvfits files are passed as input (e.g., high and low band) then use both datasets, +# but do not form closure quantities between the two datasets +else: + # load the two uvfits files + obs1 = eh.obsdata.load_uvfits(obsfile) + obs2 = eh.obsdata.load_uvfits(args.infile2) + + # Average data based on individual scan lengths + obs1.add_scans() + obs2.add_scans() + obs1 = obs1.avg_coherent(0.,scan_avg=True) + obs2 = obs2.avg_coherent(0.,scan_avg=True) + + # Add a slight offset to avoid mixed closure products + obs2.data['time'] += 0.00001 + + # concatenate the observations into a single observation object + obs = obs1.copy() + obs.data = np.concatenate([obs1.data,obs2.data]) + +# Estimate the total flux density from the ALMA(AA) -- APEX(AP) zero baseline +zbl_tot = np.median(obs.unpack_bl('AA','AP','amp')['amp']) +if zbl > zbl_tot: + print('Warning: Specified total compact flux density ' + + 'exceeds total flux density measured on AA-AP!') + +# Flag out sites in the obs.tarr table with no measurements +allsites = set(obs.unpack(['t1'])['t1'])|set(obs.unpack(['t2'])['t2']) +obs.tarr = obs.tarr[[o in allsites for o in obs.tarr['site']]] +obs = eh.obsdata.Obsdata(obs.ra, obs.dec, obs.rf, obs.bw, obs.data, obs.tarr, + source=obs.source, mjd=obs.mjd, + ampcal=obs.ampcal, phasecal=obs.phasecal) + +obs_orig = obs.copy() # save obs before any further modifications + +# Rescale short baselines to excize contributions from extended flux +if zbl != zbl_tot: + rescale_zerobaseline(obs, zbl, zbl_tot, uv_zblcut) + +# Order the stations by SNR. +# This will create a minimal set of closure quantities +# with the highest snr and smallest covariance. +obs.reorder_tarr_snr() + +#------------------------------------------------------------------------------- +# Pre-calibrate the data +#------------------------------------------------------------------------------- + +obs_sc = obs.copy() # From here on out, don't change obs. Use obs_sc to track gain changes +res = obs.res() # The nominal array resolution: 1/(longest baseline) + +# Make a Gaussian prior image for maximum entropy regularization +# This Gaussian is also the initial image +gaussprior = eh.image.make_square(obs_sc, npix, fov) +gaussprior = gaussprior.add_gauss(zbl, (prior_fwhm, prior_fwhm, 0, 0, 0)) + +# To avoid gradient singularities in the first step, add an additional small Gaussians +gaussprior = gaussprior.add_gauss(zbl*1e-3, (prior_fwhm, prior_fwhm, 0, prior_fwhm, prior_fwhm)) + +# Reverse taper the observation: this enforces a maximum resolution on reconstructed features +if reverse_taper_uas > 0: + obs_sc = obs_sc.reverse_taper(reverse_taper_uas*eh.RADPERUAS) + +# Add non-closing systematic noise to the observation +obs_sc = obs_sc.add_fractional_noise(sys_noise) + +# Make a copy of the initial data (before any self-calibration but after the taper) +obs_sc_init = obs_sc.copy() + +# Self-calibrate the LMT to a Gaussian model +# (Refer to Section 4's "Pre-Imaging Considerations") +print("Self-calibrating the LMT to a Gaussian model for LMT-SMT...") + +obs_LMT = obs_sc_init.flag_uvdist(uv_max=2e9) # only consider the short baselines (LMT-SMT) +if reverse_taper_uas > 0: + # start with original data that had no reverse taper applied. + # Re-taper, if necessary + obs_LMT = obs_LMT.taper(reverse_taper_uas*eh.RADPERUAS) + +# Make a Gaussian image that would result in the LMT-SMT baseline visibility amplitude +# as estimated in Section 4's "Pre-Imaging Considerations". +# This is achieved with a Gaussian of size 60 microarcseconds and total flux of 0.6 Jy +gausspriorLMT = eh.image.make_square(obs, npix, fov) +gausspriorLMT = gausspriorLMT.add_gauss(0.6, (60.0*eh.RADPERUAS, 60.0*eh.RADPERUAS, 0, 0, 0)) + +# Self-calibrate the LMT visibilities to the gausspriorLMT image +# to enforce the estimated LMT-SMT visibility amplitude +caltab = eh.selfcal(obs_LMT, gausspriorLMT, sites=['LM'], gain_tol=1.0, + method='both', ttype=ttype, caltable=True) + +# Spply the calibration solution to the full (and potentially tapered) dataset +obs_sc = caltab.applycal(obs_sc, interp='nearest', extrapolate=True) + +#------------------------------------------------------------------------------- +# Reconstruct an image +#------------------------------------------------------------------------------- + + +# First Round of Imaging +#------------------------- +print("Round 1: Imaging with visibility amplitudes and closure quantities...") + +# Initialize imaging with a Gaussian image +imgr = eh.imager.Imager(obs_sc, gaussprior, prior_im=gaussprior, + flux=zbl, data_term=data_term, maxit=maxit, + norm_reg=True, systematic_noise=systematic_noise, + reg_term=reg_term, ttype=ttype, cp_uv_min=uv_zblcut, stop=stop) + +# Imaging +imgr.make_image_I(show_updates=False) +converge() + +# Self-calibrate to the previous model (phase-only); +# The solution_interval is 0 to align phases from high and low bands if needed +obs_sc = eh.selfcal(obs_sc, imgr.out_last(), method='phase', ttype=ttype, solution_interval=0.0) + + +# Second Round of Imaging +#------------------------- +print("Round 2: Imaging with visibilities and closure quantities...") + +# Blur the previous reconstruction to the intrinsic resolution of ~25 uas +init = imgr.out_last().blur_circ(res) + +# Increase the weights on the data terms and reinitialize imaging +data_term_intermediate = {'vis':imgr.dat_terms_last()['amp']*10, + 'cphase':imgr.dat_terms_last()['cphase']*10, + 'logcamp':imgr.dat_terms_last()['logcamp']*10} + +imgr = eh.imager.Imager(obs_sc, init, prior_im=gaussprior, flux=zbl, + data_term=data_term_intermediate, maxit=maxit, norm_reg=True, + systematic_noise=systematic_noise, reg_term = reg_term, ttype=ttype, + cp_uv_min=uv_zblcut, stop=stop) +# Imaging +imgr.make_image_I(show_updates=False) +converge() + +# Self-calibrate to the previous model starting from scratch +# phase for all sites; amplitudes for LMT +obs_sc = eh.selfcal(obs_sc_init, imgr.out_last(), method='phase', ttype=ttype) +caltab = eh.selfcal(obs_sc, imgr.out_last(), sites=['LM'], method='both', gain_tol=gain_tol, + ttype=ttype, caltable=True) +obs_sc = caltab.applycal(obs_sc, interp='nearest',extrapolate=True) + + +# Third and Fourth Rounds of Imaging +#----------------------------------- +print("Rounds 3+4: Imaging with visibilities and closure quantities...") + +# Increase the data weights before imaging again +data_term_final = {'vis':imgr.dat_terms_last()['vis']*5, + 'cphase':imgr.dat_terms_last()['cphase']*2, + 'logcamp':imgr.dat_terms_last()['logcamp']*2} + +# Repeat imaging twice +for repeat_selfcal in range(2): + # Blur the previous reconstruction to the intrinsic resolution of ~25 uas + init = imgr.out_last().blur_circ(res) + + # Reinitialize imaging now using complex visibilities; common systematic noise + imgr = eh.imager.Imager(obs_sc, init, prior_im=gaussprior, flux=zbl, + data_term=data_term_final, maxit=maxit, norm_reg=True, + systematic_noise=0.01, reg_term=reg_term, ttype=ttype, + cp_uv_min=uv_zblcut, stop=stop) + # Imaging + imgr.make_image_I(show_updates=False) + converge() + + # Self-calibrate + caltab = eh.selfcal(obs_sc_init, imgr.out_last(), method='both', + sites=['SM','JC','AA','AP','LM','SP','AZ'], + ttype=ttype, gain_tol=gain_tol, caltable=True) + obs_sc = caltab.applycal(obs_sc_init, interp='nearest',extrapolate=True) + +#------------------------------------------------------------------------------- +# Output the results +#------------------------------------------------------------------------------- + +# Save the reconstructed image +im_out = imgr.out_last().copy() + +# If an inverse taper was used, restore the final image +# to be consistent with the original data +if reverse_taper_uas > 0.0: + im_out = im_out.blur_circ(reverse_taper_uas*eh.RADPERUAS) + +# Save the final image +im_out.save_fits(args.outfile) + +# Optionally save a pdf of the final image +if args.savepdf: + pdfout = os.path.splitext(args.outfile)[0] + '.pdf' + im_out.display(cbar_unit=['Tb'],label_type='scale',export_pdf=pdfout) + +# Optionally generate a summary of the final image and associated data consistency metrics +if args.imgsum: + + # Add a large gaussian component to account for the missing flux + # so the final image can be compared with the original data + im_addcmp = im_out.add_zblterm(obs_orig, uv_zblcut, debias=True) + obs_sc_addcmp = eh.selfcal(obs_orig, im_addcmp, method='both', ttype=ttype) + + # Save an image summary sheet + # Note! these chi2 values will not be identical to what is shown in the paper + # because they are combining high and low band + matplotlib.pyplot.close('all') + outimgsum = os.path.splitext(args.outfile)[0] + '_imgsum.pdf' + eh.imgsum(im_addcmp, obs_sc_addcmp, obs_orig, outimgsum ,cp_uv_min=uv_zblcut) From a2474467b8e844b508956eb4fda7d032a6e25eb7 Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Thu, 4 Apr 2024 15:39:00 -0700 Subject: [PATCH 05/11] Set xbit on shell file --- src/smili/run_complete.sh | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 src/smili/run_complete.sh diff --git a/src/smili/run_complete.sh b/src/smili/run_complete.sh old mode 100644 new mode 100755 From 6a20883ff30542f886a1134b1c0a2e46035f6084 Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Thu, 4 Apr 2024 15:41:13 -0700 Subject: [PATCH 06/11] Add missing files --- src/smili/example_driver.py | 132 +++++++++ src/smili/run.sh | 20 ++ src/smili/run_postprocessing.sh | 0 src/smili/smili_imaging_pipeline.py | 439 ++++++++++++++++++++++++++++ 4 files changed, 591 insertions(+) create mode 100755 src/smili/example_driver.py create mode 100755 src/smili/run.sh mode change 100644 => 100755 src/smili/run_postprocessing.sh create mode 100755 src/smili/smili_imaging_pipeline.py diff --git a/src/smili/example_driver.py b/src/smili/example_driver.py new file mode 100755 index 0000000..f5a3152 --- /dev/null +++ b/src/smili/example_driver.py @@ -0,0 +1,132 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +"""Example driver for the SMILI M87 Stokes I Imaging Pipeline for EHT observations in April 2017 + +Authors: The Event Horizon Telescope Collaboration et al. +Date: April 10, 2019 +Primary Reference: The Event Horizon Telescope Collaboration, et al. 2019d, ApJL, 875, L4 (M87 Paper IV) +Data Product Code: 2019-D01-02 + +Brief Description: +The script is an example driver of the SMILI M87 Stokes I Imaging Pipeline +(smili_imaging_pipeline.py) for EHT observations in April 2017 attached in the +same directory. For more detail instructions for smili_imaging_pipeline.py +please read the help document associated in the imaging script +"python smili_imaging_pipeline.py --help". + +This sript reconstructs images on fiducial parameters (see Section 6 of M87 +Paper IV) from calibratd uvfits files released in the Data Product Release +(2019-D01-01; see also M87 Paper III). You can run it by +"python example_driver.py --uvfitsdir --nproc ." +For details, please have a look at the help document by +"python example_driver.py --help". + +References: + - EHT Collaboration Data Portal Website: + https://eventhorizontelescope.org/for-astronomers/data + - The Event Horizon Telescope Collaboration, et al. 2019a, ApJL, 875, L1 (M87 Paper I) + - The Event Horizon Telescope Collaboration, et al. 2019b, ApJL, 875, L2 (M87 Paper II) + - The Event Horizon Telescope Collaboration, et al. 2019c, ApJL, 875, L3 (M87 Paper III) + - The Event Horizon Telescope Collaboration, et al. 2019d, ApJL, 875, L4 (M87 Paper IV) + - The Event Horizon Telescope Collaboration, et al. 2019e, ApJL, 875, L5 (M87 Paper V) + - The Event Horizon Telescope Collaboration, et al. 2019f, ApJL, 875, L6 (M87 Paper VI) + - Akiyama, K., Ikeda, S., Pleau, M., et al. 2017a, ApJ, 838, 1 + - Akiyama, K., Kuramochi, K., Ikeda, S., et al. 2017b, AJ, 153, 159 + - Akiyama, K., Tazaki, F., Moriyama, K., et al. 2019, Zenodo (SMILI version 0.0.0) + - SMILI: https://github.com/astrosmili/smili +""" + + +#------------------------------------------------------------------------------- +# Information of Authors +#------------------------------------------------------------------------------- +__author__ = "The Event Horizon Telescope Collaboration et al." +__copyright__ = "Copyright 2019, the Event Horizon Telescope Collaboration et al." +__license__ = "GPL version 3" +__version__ = "1.0" +__date__ = "April 10 2019" + + +#------------------------------------------------------------------------------- +# Modules +#------------------------------------------------------------------------------- +import os +import argparse +import glob +import tqdm +from multiprocessing import Pool + + +#------------------------------------------------------------------------------- +# Loading command line arguments +#------------------------------------------------------------------------------- +parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) +parser.add_argument('--uvfitsdir',metavar='input uvfits directory',type=str,default='../../data/uvfits/', + help='input uvfits directory of the data product release 2019-XX-XXXX') +parser.add_argument('--outputdir',metavar='output fits file',type=str,default="./smili_reconstructions", + help='output directory. The default is "./smili_reconstructions".') +parser.add_argument('--pipeline',metavar='pipeline script',type=str,default="smili_imaging_pipeline.py", + help='the filename of the pipeline scripts. The default is "smili_imaging_pipeline.py"') +parser.add_argument("--nproc",metavar="Number of parallel processors. Default=4.",type=int,default=4) +args = parser.parse_args() + +nproc = args.nproc +uvfitsdir = args.uvfitsdir +outputdir = args.outputdir +pipeline = args.pipeline + + +#------------------------------------------------------------------------------- +# Preparation +#------------------------------------------------------------------------------- +# create the output directory +if not os.path.isdir(outputdir): + command = "mkdir -p %s"%(outputdir) + print(command) + os.system(command) + +# get uvfits files +uvfitsfiles = sorted(glob.glob(os.path.join(uvfitsdir, "*_hops_netcal_StokesI.uvfits"))) +if len(uvfitsfiles) < 1: + raise ValueError("No input uvfits files in %s"%(uvfitsdir)) + +# copy uvfits files +for uvfitsfile in uvfitsfiles: + uvfitsfilebase = os.path.basename(uvfitsfile) + command = "cp -p %s %s"%(uvfitsfile, os.path.join(outputdir, uvfitsfilebase)) + print(command) + os.system(command) + +# generating commands +commands = [] +for uvfitsfile in uvfitsfiles: + uvfitsfilebase = os.path.basename(uvfitsfile) + + # check observing errors + obsdate = None + for day in [95,96,100,101]: + if "_%03d_"%(day) in uvfitsfilebase: + obsdate = day - 90 + break + if obsdate is None: + raise ValueError("Apparently there is no obsday information in the filename of the uvfits: %s"%(uvfitsfilebase)) + + command = [] + command.append("python") + command.append(pipeline) + command.append("-i %s"%(os.path.join(outputdir, uvfitsfilebase))) + command.append("--day %d"%(obsdate)) + command.append("--nproc 1") + commands.append(" ".join(command)) + +#------------------------------------------------------------------------------- +# Running the pipeline +#------------------------------------------------------------------------------- +# OMP NUM THREADS must be set to 1, so that each command will use only a single processor. +os.environ['OMP_NUM_THREADS'] = '1' + +# running imaging in parallel +pool = Pool(processes=nproc) +for _ in tqdm.tqdm(pool.imap_unordered(os.system,commands),total=len(commands)): + pass +exit() diff --git a/src/smili/run.sh b/src/smili/run.sh new file mode 100755 index 0000000..1474e35 --- /dev/null +++ b/src/smili/run.sh @@ -0,0 +1,20 @@ +#!/usr/bin/env bash +# +# Copyright (C) 2019 The Event Horizon Telescope Collaboration +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +set -e + +./example_driver.py --uvfitsdir ../../data_validation/data/uvfits diff --git a/src/smili/run_postprocessing.sh b/src/smili/run_postprocessing.sh old mode 100644 new mode 100755 diff --git a/src/smili/smili_imaging_pipeline.py b/src/smili/smili_imaging_pipeline.py new file mode 100755 index 0000000..6835a26 --- /dev/null +++ b/src/smili/smili_imaging_pipeline.py @@ -0,0 +1,439 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +"""SMILI M87 Stokes I Imaging Pipeline for EHT observations in April 2017 + +Authors: The Event Horizon Telescope Collaboration et al. +Date: April 10, 2019 +Primary Reference: The Event Horizon Telescope Collaboration, et al. 2019d, ApJL, 875, L4 (M87 Paper IV) +Data Product Code: 2019-D01-02 + +Brief Description: +The pipeline reconstructs an image from a specified uvfits file simultaneously +released in the EHT website (data release ID: 2019-D01-01) using a +python-interfaced imaging package SMILI (Akiyama et al. 2017a,b) version 0.0.0 +(Akiyama et al. 2019).This script assumes to load calibrated data sets at Stokes +I (see M87 Paper III). As described in M87 Paper IV Section 6.2.3, the SMILI +pipeline also uses the eht-imaging library (Chael et al. 2016, 2018) +version 1.1.0 (Chael et al. 2019) solely for pre-imaging calibrations including +time averaging, LMT calibration to use the input visibility data sets consistent +with the eht-imaging imaging pipeline. The script has been tested in Python 2.7 +installed with Anaconda on Ubuntu 18.04 LTS and MacOS 10.13 & 10.14 (with +macport or homebrew). + +Notes: +We note that, as described in 2019-D01-01, released visibility data sets are +slightly different from data sets used in Paper IV for two reasons. + +(1) They have only Stokes I, while Paper IV data sets have dual polarization at +Stokes RR and LL. This slight difference in released data sets will change +self-calibration procedures slightly in this pipeline since R and L gains are +calibrated separately for the latter dual polarization data sets. We find that +this does not produce any significant differences (within < ~5%) in resultant images. + +(2) In Paper IV, this pipeline strictly uses Stoke I data; JCMT which has a +single polarization is flagged when Stokes I visibilities are computed from RR +and LL, while JCMT is included for self-calibration at the corresponding +polarization. On the other hand, in the released data sets, JCMT has pseudo +Stokes RR/LL data, such that Stokes I computed from them is identical to the +original single Stokes. This will make JCMT be used in imaging, which change +images slightly. For reproducibility of Paper IV results using this data set, in +default, we will not use JCMT for imaging, as we did in Paper IV. If you want +to use it to use all of data sets, you can specify --keepsinglepol to keep it. + +The pipeline will output three files. + - image (specified with -o option; assumed to be xxxx.fits here) + - pre-calibrated uvfits file (xxxx.precal.uvfits) + - self-calibrated uvfits file (xxxx.selfcal.uvfits) +For usage and detail parameters of the pipeline, please +read the help document associated in the imaging script +"python smili_imaging_pipeline.py --help". + +References: + - EHT Collaboration Data Portal Website: + https://eventhorizontelescope.org/for-astronomers/data + - The Event Horizon Telescope Collaboration, et al. 2019a, ApJL, 875, L1 (M87 Paper I) + - The Event Horizon Telescope Collaboration, et al. 2019b, ApJL, 875, L2 (M87 Paper II) + - The Event Horizon Telescope Collaboration, et al. 2019c, ApJL, 875, L3 (M87 Paper III) + - The Event Horizon Telescope Collaboration, et al. 2019d, ApJL, 875, L4 (M87 Paper IV) + - The Event Horizon Telescope Collaboration, et al. 2019e, ApJL, 875, L5 (M87 Paper V) + - The Event Horizon Telescope Collaboration, et al. 2019f, ApJL, 875, L6 (M87 Paper VI) + - Akiyama, K., Ikeda, S., Pleau, M., et al. 2017a, ApJ, 838, 1 + - Akiyama, K., Kuramochi, K., Ikeda, S., et al. 2017b, AJ, 153, 159 + - Akiyama, K., Tazaki, F., Moriyama, K., et al. 2019, Zenodo (SMILI version 0.0.0) + - Chael, A., Johnson, M., Narayan, R., et al. 2016, ApJ, 829, 11C + - Chael, A., Johnson, M., Bouman, K., et al. 2018, ApJ, 857, 23C + - Chael, A., Bouman, K., Johnson, M. et al., Zenodo (eht-imaging version 1.1.0) + - Anaconda: https://www.anaconda.com/ + - eht-imaging: https://github.com/achael/eht-imaging + - SMILI: https://github.com/astrosmili/smili +""" + + +#------------------------------------------------------------------------------- +# Information of Authors +#------------------------------------------------------------------------------- +__author__ = "The Event Horizon Telescope Collaboration et al." +__copyright__ = "Copyright 2019, the Event Horizon Telescope Collaboration et al." +__license__ = "GPL version 3" +__version__ = "1.0" +__date__ = "April 10 2019" + + +#------------------------------------------------------------------------------- +# Modules +#------------------------------------------------------------------------------- +from smili import uvdata,imdata,imaging,util +import ehtim as eh +import pandas as pd +import numpy as np +import copy +import os +import argparse + + +#------------------------------------------------------------------------------- +# Loading command line arguments +#------------------------------------------------------------------------------- +parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) +parser.add_argument('-i','--input',metavar='input uvfits file',type=str,required=True, + help='input uvfits file name ') +parser.add_argument('-o','--output',metavar='output fits file',type=str,default=None, + help='output fits file name') +parser.add_argument('--day',metavar='observing day of April 2017',type=int,required=True, + help='Observational Day of April: 5, 6, 10 or 11. This must be specified correctly to apply the correct total flux density used in the network calibration.') +parser.add_argument("--prior",metavar='prior FWHM size',type=float,default=40., + help='Prior FWHM (uas)') +parser.add_argument("--tfd",metavar='total flux density',type=float,default=0.5, + help='Target total flux density (Jy)') +parser.add_argument("--sys",metavar='fractional systematic error',type=float,default=0., + help='Fractional systematic error to be added to visibilities in quadrature (percent)') +parser.add_argument("--keepsinglepol", action='store_true', default=False, + help='This option uses single pol station(s) for imaging; this will cause slight difference in the resultant images from M87 Paper IV.') +parser.add_argument("--lambl1",metavar='lambda l1',type=float,default=1., + help='lambda l1') +parser.add_argument("--lambtv",metavar='lambda tv',type=float,default=1e3, + help='lambda tv') +parser.add_argument("--lambtsv",metavar='lambda tsv',type=float,default=1e4, + help='lambda tsv') +parser.add_argument("--nproc",metavar="Number of parallel processors. Default=1.",type=int,default=1) +args = parser.parse_args() + +# Input uvfits file +inputuvfits = args.input + +# Output directory +outputfits = args.output +if outputfits is None: + outputfits = os.path.splitext(inputuvfits)[0]+".fits" +outputhead = os.path.splitext(outputfits)[0] + +# Observational date +day = args.day +if day not in [5,6,10,11]: + raise ValueError("The observational day (--day integer) must be (April) 5, 6, 10 or 11.") + +# Prior size +prior_fwhm = args.prior + +# Target total flux density +totalflux = args.tfd + +# Systematic error +syserr = args.sys + +# Lambda L1, TV and TSV +lambl1 = args.lambl1 +lambtv = args.lambtv +lambtsv = args.lambtsv + +# Number of processors +nproc = args.nproc +os.environ['OMP_NUM_THREADS'] = '%d'%(nproc) + +#------------------------------------------------------------------------------- +# Other tuning parameters fixed in the paper. +#------------------------------------------------------------------------------- +# The station(s) with single pol data +singlepol = ['JC'] + +# Number of iterative imaging inside a single selfcal cycle +Nweig = 3 + +# Number of closure imaging attempted in this script (must be < Nself) +Nflcl = 2 + +# Number of selfcals +Nself = 4 + +# Imaging pixel size +dx_uas = 2 + +# Number of pixels +nx = 64 + + +#------------------------------------------------------------------------------- +# Set default imaging parameters +#------------------------------------------------------------------------------- +# Definition of the prior image +def gen_prior_gauss(x0=0.,y0=0.,fwhm=40.,dx=dx_uas,nx=nx): + ''' + This function provides a Gaussain image with a specified FWHM in uas. + ''' + # create a blank image + outimage = imdata.IMFITS(dx=dx,nx=nx,angunit="uas") + outimage = outimage.add_gauss(majsize=fwhm, overwrite=True) + return outimage + +# Set a circular Gaussian prior for L1 +l1_prior = gen_prior_gauss(fwhm=prior_fwhm,dx=dx_uas,nx=nx) + +# Default imaging parameters +imprm_init={} + +# Iteration numbers +imprm_init["niter"] = 1000 + +# Regularization parameters +# Flat prior TSV +if lambtsv > 0: + imprm_init["tsv_lambda"] = lambtsv +else: + imprm_init["tsv_lambda"] = -1 + +# Flat prior TV +if lambtv > 0: + imprm_init["tv_lambda"] = lambtv +else: + imprm_init["tv_lambda"] = -1 + +# Weighted L1 +if lambl1 > 0: + imprm_init["l1_lambda"] = lambl1 + imprm_init["l1_prior"] = l1_prior +else: + imprm_init["l1_lambda"] = -1 + +# Maximum Entropy Methods: not used in SMILI pipeline +# GS-MEM +imprm_init["gs_lambda"] = -1 +# KL-MEM +imprm_init["kl_lambda"] = -1 + +# Total flux density constraint +imprm_init["totalflux"] = totalflux # Target total flux density +imprm_init["tfd_lambda"] = 1 # Lambda parameter +imprm_init["tfd_tgterror"] = 1e-2 # Target fractional accuracy + + +#------------------------------------------------------------------------------- +# Step 1: Pre-calibration +# To make sure that we use visibility data sets pre-calibrated consistently with +# two other pipelines, this pipeline also uses eht-imaging library. +#------------------------------------------------------------------------------- +# Load the uvfits file +obs = eh.obsdata.load_uvfits(inputuvfits).switch_polrep('stokes') + +# Rescale short baselines to excite contributions from extended flux. +# Total flux density adopted in the Network Calibration +# see EHT Collaboration et al. 2019c, M87 Paper III +if day < 7: + netcal_tfd = 1.1 # 1.1 Jy was adopted for Apr 5 and 6 +else: + netcal_tfd = 1.2 # 1.2 Jy was adopted for Apr 10 and 11 +# Scaling the total flux density: +# setting totalflux < netcal_tfd assumes there is an extended constant flux +# component of netcal_tfd-totalflux Jy +for j in range(len(obs.data)): + if (obs.data['u'][j]**2 + obs.data['v'][j]**2)**0.5 < 0.1e9: + for k in range(-8,0): + obs.data[j][k] *= totalflux/netcal_tfd + +# Do scan averaging +obs.add_scans() # this seperates the data into scans, if it isn't done so already with an NX table +obs = obs.avg_coherent(0.,scan_avg=True) # average each scan coherantly + +# Order stations. this is to create a minimal set of closure quantities with the highest snr +obs.reorder_tarr_snr() + +# Self-calibrate the LMT to a Gaussian model if LMT is included in the input uvfits file +if "LM" in obs.data["t1"] or "LM" in obs.data["t2"]: + obs_selfcal = obs.copy() + # Add systematic noise for leakage and other non-closing errors + # (reminder: this must be done *after* any averaging) + obs_selfcal = obs_selfcal.add_fractional_noise(syserr) + obs_selfcal = obs_selfcal.switch_polrep('stokes') + # Set a circular Gaussian with a FWHM of 60 uas as the model image + gausspriormodel = eh.image.make_square(obs, int(nx), nx*dx_uas*eh.RADPERUAS) # Create empty image + gausspriormodel = gausspriormodel.add_gauss(0.6, (60*eh.RADPERUAS, 60*eh.RADPERUAS, 0, 0, 0)) # Add gaussian + # Self-calibrate amplitudes + for repeat in range(3): + caltab = eh.self_cal.self_cal( + obs_selfcal.flag_uvdist(uv_max=2e9), + gausspriormodel, + sites=['LM','LM'], + method='vis', + ttype='nfft', + processes=nproc, + caltable=True, + gain_tol=1.0) + obs_selfcal = caltab.applycal(obs_selfcal, interp='nearest', extrapolate=True) + # Save calibrated uvfits data sets + obs_selfcal.save_uvfits(outputhead+".precal.uvfits") + del obs, obs_selfcal, caltab, gausspriormodel, netcal_tfd +else: + obs.save_uvfits(outputhead+".precal.uvfits") + del obs + + +#--------------------------------------------------------------------------- +# Step 2: Generating Data Tables for Imaging +#--------------------------------------------------------------------------- +# Create the full complex visibility table +vtable = uvdata.UVFITS(outputhead+".precal.uvfits").select_stokes("I").make_vistable() +if not args.keepsinglepol: + for stname in singlepol: + vtable = vtable.query("st1name != '"+stname+"' and st2name != '"+stname+"'").reset_index(drop=True) + +# Check the number of data sets +if len(vtable.amp) == 0: + raise ValueError("No data points exist in the input visibility data set.") + +# Create tables for non-trivial/non-redundant bi-spectrum and closure amplitudes +btable = vtable.make_bstable(redundant=[["AA","AP"],["JC","SM"]]) +ctable = vtable.make_catable(redundant=[["AA","AP"],["JC","SM"]]) + +# Check the number of data sets +if len(btable.amp) == 0: + btable = None +if len(ctable.amp) == 0: + ctable = None + +#--------------------------------------------------------------------------- +# Step 3: Imaging +#--------------------------------------------------------------------------- +def edit_image(image,imprm): + ''' + Image editing (shifting, blurring) + ''' + image_edited = copy.deepcopy(image) + orgtotalflux = image_edited.totalflux() + + # Shifting images such that the ring or dominant structure comes to the center + if "vistable" not in imprm.keys(): + image_edited = image_edited.comshift() + + # Blurring with 20 uas Gaussian beam + image_edited = image_edited.convolve_gauss(majsize=20,angunit="uas") + + # Scaling images + if "totalflux" in imprm.keys(): + image_edited.data *= imprm["totalflux"]/image_edited.totalflux() + else: + image_edited.data *= orgtotalflux/image_edited.totalflux() + + image_edited.update_fits() + return image_edited + +# Loop for self-calibration +for iself in range(Nself): + # Set the initial image + # At each iteration of self-calibrations, imaging starts from a circular + # Gaussian. + initimage = imdata.IMFITS( + dx=dx_uas, + nx=nx, + angunit="uas", + uvfits=outputhead+".precal.uvfits", + ) + initimage = initimage.add_gauss(totalflux=totalflux,majsize=prior_fwhm,overwrite=True) + + # Iterative Imaging: + # At each iteration of self-calibration, imaging iterations are attempted + # for Nweig(=5) times. + for iweig in range(Nweig): + # Imaging parameters: copied from the default parameter sets + imprm = copy.deepcopy(imprm_init) + + # Data to be used in default + # Closure phases/amplitudes are used always to avoid over-fitting to + # full complex visibilities + imprm["bstable"]=btable + imprm["catable"]=ctable + + # Use amplitudes or full complex visibilities: + # Parameters depending on self-cal stages + if iself < Nflcl: + # The first half for self-calibration: using closure quantities + + # Get Amplitude Data sets + atable = copy.deepcopy(vtable) + + # Flag the intra-site to avoid conflict with with the total-flux-density regularization. + atable = atable.query("uvdist > 50e6") + + # Adding a fractional error to the data + # For LMT (adding 30% error) + idx = atable["st1name"] == "LM" + idx|= atable["st2name"] == "LM" + atable.loc[idx,"sigma"] = np.sqrt(atable.loc[idx,"sigma"]**2+(atable.loc[idx,"amp"]*0.3)**2) + # For other stations (adding an error specified with amp_err) + idx = idx == False + atable.loc[idx,"sigma"] = np.sqrt(atable.loc[idx,"sigma"]**2+(atable.loc[idx,"amp"]*0.05)**2) + + # Use amplitudes + imprm["amptable"]=atable + else: + # The latter half for self-calibration: using full complex + # visibilities. Note that we still add closure quantities + # to prevent over-fitting to full complex visibilities + vtable_in = vtable.query("uvdist > 50e6") + imprm["vistable"] = vtable_in + + # Imaging + outimage = imaging.lbfgs.imaging(initimage,**imprm) + initimage = edit_image(outimage,imprm) + + # Save image + outimage.to_fits(outputfits) + + # If iself==Nself, self-calibration won't be performed and exit the script. + if iself == Nself-1: + break + + # Self-calibration + # 1. Set the file name + if iself == 0: + uvfitsold = outputhead+".precal.uvfits" + else: + uvfitsold = outputhead+".selfcal.uvfits" + uvfitsnew = outputhead+".selfcal.uvfits" + + # 2. Load uvfits + uvfits = uvdata.UVFITS(uvfitsold) + # + # 3. Edit images + # Before selfcal, we bring the center of the mass to the image center + # and also normalize the image with the target total flux density. + outimage = outimage.comshift() + # Normalize the total flux density + # *** Thanks to the total-flux-density regularization, this will provide only a + # *** slight change in the total flux density by a few pecent or even < 1% + outimage.data[0,0] *= totalflux/outimage.totalflux() + + # 4. Do selfcal + # This is a regularized selfcal function using Gaussian priors + # for gain amplitudes and phases. std_amp is the standard deviation of + # the Gaussian prior for amplitude gains. std_amp=10000 works as almost + # the flat prior, so it is essentially selfcal without any gain constraints. + caltable = uvfits.selfcal(outimage,std_amp=10000) + uvfits = uvfits.apply_cltable(caltable) + + # 5. Save the self-calibrated uvfits file + uvfits.to_uvfits(uvfitsnew) + + # 6. Get a new visibility table + vtable = uvfits.select_stokes("I").make_vistable() + if not args.keepsinglepol: + for stname in singlepol: + vtable = vtable.query("st1name != '"+stname+"' and st2name != '"+stname+"'").reset_index(drop=True) From ac1091f554bd1f1c34d6e8279b63560d709589a2 Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Thu, 4 Apr 2024 15:42:12 -0700 Subject: [PATCH 07/11] Add missing difmap.sh --- src/difmap/difmap.sh | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100755 src/difmap/difmap.sh diff --git a/src/difmap/difmap.sh b/src/difmap/difmap.sh new file mode 100755 index 0000000..33b4601 --- /dev/null +++ b/src/difmap/difmap.sh @@ -0,0 +1,40 @@ +#!/usr/bin/env bash +# +# Copyright (C) 2019 The Event Horizon Telescope Collaboration +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +if [ $# -eq 0 ]; then + echo "usage: $0 [input].uvfits" + exit 0 +fi + +in_name=${1%.uvfits} + +expect <" +send -- "@EHT_Difmap ${in_name},CircMask_r30_x-0.002_y0.022,-10,0.5,0.1,2,-1\r" + +expect "*0>" +send -- "exit\r" + +expect "*quit without saving: " +send -- "\r" + +expect eof +EOF From c0581e3b20a62dd5dfb0b279eea4f2521f209858 Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Thu, 4 Apr 2024 15:42:40 -0700 Subject: [PATCH 08/11] Add missing input files --- src/difmap/CircMask_r30_x-0.002_y0.022.win | 32 +++ src/difmap/afmhot_10us.cmap | 256 +++++++++++++++++++++ 2 files changed, 288 insertions(+) create mode 100644 src/difmap/CircMask_r30_x-0.002_y0.022.win create mode 100755 src/difmap/afmhot_10us.cmap diff --git a/src/difmap/CircMask_r30_x-0.002_y0.022.win b/src/difmap/CircMask_r30_x-0.002_y0.022.win new file mode 100644 index 0000000..d4b2307 --- /dev/null +++ b/src/difmap/CircMask_r30_x-0.002_y0.022.win @@ -0,0 +1,32 @@ +! CLEAN windows written by wwins in difmap. +! Columns specify xmin xmax ymin ymax (mas) of each CLEAN window. + -0.0319833287 0.0279833287 0.0210000000 0.0230000000 + -0.0318496231 0.0278496231 0.0230000000 0.0250000000 + -0.0315803989 0.0275803989 0.0250000000 0.0270000000 + -0.0311719043 0.0271719043 0.0270000000 0.0290000000 + -0.0306181760 0.0266181760 0.0290000000 0.0310000000 + -0.0299105715 0.0259105715 0.0310000000 0.0330000000 + -0.0290370117 0.0250370117 0.0330000000 0.0350000000 + -0.0279807621 0.0239807621 0.0350000000 0.0370000000 + -0.0267184142 0.0227184142 0.0370000000 0.0390000000 + -0.0252163735 0.0212163735 0.0390000000 0.0410000000 + -0.0234242853 0.0194242853 0.0410000000 0.0430000000 + -0.0212613603 0.0172613603 0.0430000000 0.0450000000 + -0.0185831240 0.0145831240 0.0450000000 0.0470000000 + -0.0150766968 0.0110766968 0.0470000000 0.0490000000 + -0.0096811457 0.0056811457 0.0490000000 0.0510000000 + -0.0319833287 0.0279833287 0.0230000000 0.0210000000 + -0.0318496231 0.0278496231 0.0210000000 0.0190000000 + -0.0315803989 0.0275803989 0.0190000000 0.0170000000 + -0.0311719043 0.0271719043 0.0170000000 0.0150000000 + -0.0306181760 0.0266181760 0.0150000000 0.0130000000 + -0.0299105715 0.0259105715 0.0130000000 0.0110000000 + -0.0290370117 0.0250370117 0.0110000000 0.0090000000 + -0.0279807621 0.0239807621 0.0090000000 0.0070000000 + -0.0267184142 0.0227184142 0.0070000000 0.0050000000 + -0.0252163735 0.0212163735 0.0050000000 0.0030000000 + -0.0234242853 0.0194242853 0.0030000000 0.0010000000 + -0.0212613603 0.0172613603 0.0010000000 -0.0010000000 + -0.0185831240 0.0145831240 -0.0010000000 -0.0030000000 + -0.0150766968 0.0110766968 -0.0030000000 -0.0050000000 + -0.0096811457 0.0056811457 -0.0050000000 -0.0070000000 diff --git a/src/difmap/afmhot_10us.cmap b/src/difmap/afmhot_10us.cmap new file mode 100755 index 0000000..e8ef6ee --- /dev/null +++ b/src/difmap/afmhot_10us.cmap @@ -0,0 +1,256 @@ +0.092633 0.067057 0.061824 1 +0.113669 0.062030 0.051540 1 +0.140569 0.050091 0.033232 1 +0.165165 0.033736 0.017454 1 +0.188428 0.014065 0.005529 1 +0.205142 0.000004 0.000001 1 +0.211265 0.000002 0.000001 1 +0.217357 0.000006 0.000002 1 +0.223424 0.000007 0.000002 1 +0.229468 0.000006 0.000002 1 +0.235490 0.000001 0.000000 1 +0.241483 0.000005 0.000002 1 +0.247456 0.000007 0.000002 1 +0.253410 0.000006 0.000002 1 +0.259346 0.000002 0.000001 1 +0.265259 0.000004 0.000001 1 +0.271155 0.000007 0.000002 1 +0.277036 0.000007 0.000002 1 +0.282902 0.000002 0.000001 1 +0.288750 0.000005 0.000001 1 +0.294582 0.000007 0.000002 1 +0.300004 0.000882 0.000282 1 +0.304945 0.002842 0.000922 1 +0.309881 0.004862 0.001598 1 +0.314817 0.006939 0.002307 1 +0.319751 0.009070 0.003050 1 +0.324686 0.011255 0.003824 1 +0.329623 0.013488 0.004629 1 +0.334563 0.015767 0.005462 1 +0.339508 0.018090 0.006322 1 +0.344458 0.020451 0.007207 1 +0.349415 0.022851 0.008117 1 +0.354379 0.025285 0.009048 1 +0.359352 0.027748 0.009999 1 +0.364335 0.030239 0.010968 1 +0.369329 0.032752 0.011953 1 +0.374334 0.035286 0.012952 1 +0.379352 0.037836 0.013964 1 +0.384383 0.040400 0.014985 1 +0.389428 0.042887 0.016014 1 +0.394488 0.045292 0.017048 1 +0.399565 0.047620 0.018086 1 +0.404658 0.049875 0.019124 1 +0.409768 0.052058 0.020162 1 +0.414895 0.054174 0.021197 1 +0.420042 0.056222 0.022226 1 +0.425208 0.058204 0.023248 1 +0.430393 0.060122 0.024259 1 +0.435600 0.061976 0.025257 1 +0.440827 0.063768 0.026241 1 +0.446075 0.065500 0.027209 1 +0.451346 0.067170 0.028158 1 +0.456639 0.068780 0.029086 1 +0.461954 0.070329 0.029991 1 +0.467294 0.071817 0.030870 1 +0.472657 0.073245 0.031721 1 +0.477999 0.074767 0.032190 1 +0.483310 0.076423 0.032180 1 +0.488575 0.078253 0.031579 1 +0.493868 0.080011 0.030948 1 +0.499185 0.081709 0.030255 1 +0.504533 0.083330 0.029546 1 +0.509905 0.084892 0.028776 1 +0.515305 0.086387 0.027966 1 +0.520733 0.087815 0.027114 1 +0.526186 0.089186 0.026196 1 +0.531658 0.090523 0.025161 1 +0.537149 0.091825 0.024017 1 +0.542657 0.093097 0.022742 1 +0.548182 0.094341 0.021334 1 +0.553727 0.095546 0.019813 1 +0.559287 0.096731 0.018140 1 +0.564868 0.097876 0.016358 1 +0.570467 0.098993 0.014437 1 +0.576088 0.100068 0.012387 1 +0.581750 0.101056 0.010171 1 +0.587467 0.101921 0.007745 1 +0.593246 0.102644 0.005128 1 +0.598420 0.104881 0.003606 1 +0.602699 0.109285 0.003830 1 +0.606968 0.113646 0.004025 1 +0.611229 0.117960 0.004211 1 +0.615481 0.122237 0.004377 1 +0.619725 0.126477 0.004528 1 +0.623961 0.130683 0.004664 1 +0.628190 0.134862 0.004781 1 +0.632411 0.139011 0.004889 1 +0.636624 0.143139 0.004972 1 +0.640831 0.147240 0.005052 1 +0.645030 0.151324 0.005105 1 +0.649223 0.155386 0.005154 1 +0.653408 0.159434 0.005180 1 +0.657587 0.163463 0.005201 1 +0.661759 0.167480 0.005201 1 +0.665926 0.171482 0.005195 1 +0.670086 0.175474 0.005171 1 +0.674240 0.179454 0.005139 1 +0.678388 0.183424 0.005093 1 +0.682531 0.187385 0.005037 1 +0.686668 0.191339 0.004970 1 +0.690800 0.195285 0.004892 1 +0.694926 0.199224 0.004805 1 +0.699047 0.203159 0.004707 1 +0.703163 0.207087 0.004603 1 +0.707273 0.211012 0.004487 1 +0.711380 0.214933 0.004366 1 +0.715481 0.218851 0.004235 1 +0.719577 0.222766 0.004100 1 +0.723669 0.226679 0.003956 1 +0.727757 0.230590 0.003809 1 +0.731840 0.234500 0.003653 1 +0.735920 0.238409 0.003496 1 +0.739995 0.242319 0.003332 1 +0.744066 0.246227 0.003167 1 +0.748133 0.250137 0.002997 1 +0.752196 0.254047 0.002828 1 +0.756256 0.257959 0.002654 1 +0.760312 0.261871 0.002481 1 +0.764366 0.265784 0.002301 1 +0.768421 0.269694 0.002110 1 +0.772475 0.273603 0.001905 1 +0.776530 0.277510 0.001690 1 +0.780585 0.281416 0.001460 1 +0.784640 0.285322 0.001221 1 +0.788697 0.289227 0.000968 1 +0.792753 0.293133 0.000706 1 +0.796811 0.297038 0.000432 1 +0.800868 0.300945 0.000148 1 +0.804892 0.304893 0.000003 1 +0.808882 0.308882 0.000001 1 +0.812871 0.312872 0.000003 1 +0.816862 0.316862 0.000001 1 +0.820854 0.320854 0.000002 1 +0.824847 0.324847 0.000002 1 +0.828841 0.328841 0.000002 1 +0.832837 0.332837 0.000002 1 +0.836834 0.336835 0.000002 1 +0.840834 0.340834 0.000002 1 +0.844835 0.344836 0.000002 1 +0.848837 0.348842 0.000012 1 +0.852838 0.352854 0.000042 1 +0.856837 0.356873 0.000091 1 +0.860835 0.360898 0.000162 1 +0.864831 0.364930 0.000255 1 +0.868827 0.368969 0.000372 1 +0.872822 0.373015 0.000514 1 +0.876816 0.377068 0.000681 1 +0.880810 0.381128 0.000876 1 +0.884803 0.385196 0.001098 1 +0.888795 0.389270 0.001352 1 +0.892787 0.393353 0.001635 1 +0.896778 0.397443 0.001952 1 +0.900770 0.401542 0.002300 1 +0.904761 0.405647 0.002686 1 +0.908752 0.409761 0.003105 1 +0.912743 0.413883 0.003565 1 +0.916735 0.418014 0.004061 1 +0.920727 0.422152 0.004601 1 +0.924718 0.426299 0.005179 1 +0.928711 0.430454 0.005806 1 +0.932704 0.434618 0.006473 1 +0.936698 0.438791 0.007191 1 +0.940692 0.442972 0.007953 1 +0.944688 0.447162 0.008769 1 +0.948684 0.451362 0.009633 1 +0.952682 0.455570 0.010553 1 +0.956680 0.459787 0.011526 1 +0.960680 0.464014 0.012556 1 +0.964682 0.468250 0.013645 1 +0.968685 0.472495 0.014793 1 +0.972689 0.476749 0.016004 1 +0.976695 0.481014 0.017278 1 +0.980703 0.485287 0.018619 1 +0.984713 0.489571 0.020025 1 +0.988725 0.493864 0.021504 1 +0.992739 0.498168 0.023050 1 +0.996217 0.503024 0.021093 1 +0.999068 0.508507 0.015182 1 +1.000000 0.515619 0.015585 1 +1.000000 0.523413 0.023380 1 +1.000000 0.531144 0.031112 1 +1.000000 0.538814 0.038828 1 +1.000000 0.546423 0.046496 1 +1.000000 0.553969 0.054033 1 +1.000000 0.561455 0.061509 1 +1.000000 0.568883 0.068927 1 +1.000000 0.576257 0.076291 1 +1.000000 0.583579 0.083603 1 +1.000000 0.590852 0.090864 1 +1.000000 0.598078 0.098079 1 +1.000000 0.605258 0.105269 1 +1.000000 0.612396 0.112414 1 +1.000000 0.619494 0.119517 1 +1.000000 0.626554 0.126579 1 +1.000000 0.633578 0.133603 1 +1.000000 0.640568 0.140590 1 +1.000000 0.647525 0.147543 1 +1.000000 0.654451 0.154464 1 +1.000000 0.661349 0.161354 1 +1.000000 0.668217 0.168220 1 +1.000000 0.675057 0.175067 1 +1.000000 0.681872 0.181885 1 +1.000000 0.688663 0.188677 1 +1.000000 0.695432 0.195445 1 +1.000000 0.702180 0.202190 1 +1.000000 0.708907 0.208913 1 +1.000000 0.715616 0.215616 1 +1.000000 0.722303 0.222309 1 +1.000000 0.728974 0.228982 1 +1.000000 0.735629 0.235638 1 +1.000000 0.742268 0.242277 1 +1.000000 0.748894 0.248900 1 +1.000000 0.755506 0.255509 1 +1.000000 0.762105 0.262107 1 +1.000000 0.768690 0.268695 1 +1.000000 0.775265 0.275271 1 +1.000000 0.781830 0.281836 1 +1.000000 0.788385 0.288390 1 +1.000000 0.794932 0.294934 1 +1.000000 0.801470 0.301471 1 +1.000000 0.807999 0.308002 1 +1.000000 0.814521 0.314526 1 +1.000000 0.821037 0.321042 1 +1.000000 0.827549 0.327552 1 +1.000000 0.834055 0.334056 1 +1.000000 0.840556 0.340557 1 +1.000000 0.847052 0.347055 1 +1.000000 0.853546 0.353549 1 +1.000000 0.860036 0.360040 1 +1.000000 0.866525 0.366528 1 +1.000000 0.873012 0.373013 1 +1.000000 0.879497 0.379499 1 +1.000000 0.885981 0.385984 1 +1.000000 0.892465 0.392468 1 +1.000000 0.898949 0.398952 1 +0.999544 0.905499 0.408416 1 +0.998740 0.912075 0.420422 1 +0.997994 0.918610 0.432492 1 +0.997305 0.925105 0.444650 1 +0.996685 0.931559 0.456866 1 +0.996137 0.937970 0.469141 1 +0.995668 0.944336 0.481489 1 +0.995280 0.950658 0.493931 1 +0.994988 0.956931 0.506436 1 +0.994794 0.963156 0.519004 1 +0.994703 0.969333 0.531642 1 +0.994720 0.975458 0.544388 1 +0.994855 0.981532 0.557202 1 +0.995115 0.987553 0.570083 1 +0.995508 0.993518 0.583031 1 +0.998300 0.998505 0.596729 1 +1.000000 1.000000 0.670757 1 +0.999999 1.000000 0.766032 1 +1.000000 1.000000 0.850975 1 +0.999997 1.000000 0.928392 1 +1.000000 1.000000 1.000000 1 From 5b560feb071fdfa96542554b78539c71409c599a Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Thu, 4 Apr 2024 16:38:40 -0700 Subject: [PATCH 09/11] Use version from inside the container --- src/smili/smili_postprocessing.py | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/src/smili/smili_postprocessing.py b/src/smili/smili_postprocessing.py index 13a9c1c..b1a367a 100644 --- a/src/smili/smili_postprocessing.py +++ b/src/smili/smili_postprocessing.py @@ -3,15 +3,17 @@ import os import argparse import ehtim as eh -import ehtplot.color +import numpy as np -# Load and parse command-line arguments +#------------------------------------------------------------------------------- +# Load command-line arguments +#------------------------------------------------------------------------------- parser = argparse.ArgumentParser(description="Script to convert EHT_Difmap fits files to pdf images") -parser.add_argument('-i', '--infile' , default="" , help="input FITS file") -parser.add_argument('-o', '--outfile' , default="" , help="output PDF file") -parser.add_argument('-s', '--scale' , default=False, help="display scale in output" , action='store_true') -parser.add_argument('-b', '--beam' , default=False, help="display beam size in output" , action='store_true') -parser.add_argument('-l', '--blur' , default=False, help="apply Gaussian blur to image" , action='store_true') +parser.add_argument('-i', '--infile' , default="" , help="input FITS file") +parser.add_argument('-o', '--outfile', default="" , help="output PDF file") +parser.add_argument('-s', '--scale' , default=False, help="display scale in output" , action='store_true') +parser.add_argument('-b', '--beam' , default=False, help="display beam size in output" , action='store_true') +parser.add_argument('-l', '--blur' , default=False, help="apply Gaussian blur to image" , action='store_true') parser.add_argument('-u', '--afmhot10us', default=False, help="use the afmhot_10us colormap" , action='store_true') parser.add_argument('-t', '--notitle' , default=False, help="display with no title or padding" , action='store_true') parser.add_argument('-a', '--all' , default=False, help="perform all post-processing steps", action='store_true') @@ -30,12 +32,16 @@ # Set params for display based on user args scale = 'scale' if (args.scale or args.all) else 'none' -cmap = 'afmhot_10us' if (args.afmhot10us or args.all) else 'afmhot' +#cmap = 'afmhot_10us' if (args.afmhot10us or args.all) else 'afmhot' title = False if (args.notitle or args.all) else True -# Load fits file into image object and set params for beam +# Load fits file into image object im_obj = eh.image.load_fits(args.infile) -params = [9.018e-11, 9.018e-11, 0] # This is for SMILI (18.6 uas) in Paper IV +params = [9.018e-11, 9.018e-11, 0] # This is for SMILI (18.6 uas) + +# Create color map of afmhot_10us, vals copied from ehtplot/color/ctabs/afmhot_10us.ctab file +colors = np.loadtxt("afmhot_10us.cmap", delimiter=" ", unpack=False) +cmap = matplotlib.colors.ListedColormap(colors) if(args.afmhot10us or args.all) else 'afmhot' # Blur the image object using Gaussian blur if specified if(args.blur or args.all): im_obj = im_obj.blur_gauss(params) @@ -46,3 +52,5 @@ export_pdf=args.outfile) else: im_obj.display(cbar_unit=['Tb'], has_title=title, cfun=cmap, cbar_orientation='horizontal', label_type=scale, export_pdf=args.outfile) + +print("================================") From 97008d8e6dc2fa62e47d72dcec0aaeb6335e9f99 Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Fri, 5 Apr 2024 10:39:42 -0700 Subject: [PATCH 10/11] Create output dir beforee running code --- src/smili/run_postprocessing.sh | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/smili/run_postprocessing.sh b/src/smili/run_postprocessing.sh index 6cc9c9b..51ccb16 100755 --- a/src/smili/run_postprocessing.sh +++ b/src/smili/run_postprocessing.sh @@ -10,6 +10,8 @@ # ('-t', '--notitle' , default=False, help="display with no title") # ('-a', '--all' , default=False, help="perform all post-processing steps") +mkdir -p post + # Traverse output dir and apply post-processing to each .fits file for d in 095 096 100 101; do python smili_postprocessing.py \ @@ -28,9 +30,9 @@ done #python smili_postprocessing.py \ # -i ./smili_reconstructions/SR1_M87_2017_101_hi_hops_netcal_StokesI.fits \ # -o ./post/SR1_M87_2017_101_afmhot10us.pdf \ -# --afmhot10us --notitle +# --afmhot10us --notitle #python smili_postprocessing.py \ # -i ./smili_reconstructions/SR1_M87_2017_101_hi_hops_netcal_StokesI.fits \ # -o ./post/SR1_M87_2017_101_afmhot10us_blur.pdf \ -# --blur --notitle --afmhot10us --beam +# --blur --notitle --afmhot10us --beam From cf999de1ee6e2c0d03e88c1dada6c6ecc8aff263 Mon Sep 17 00:00:00 2001 From: Rajiv Mayani Date: Fri, 5 Apr 2024 10:40:29 -0700 Subject: [PATCH 11/11] Add link to afmhot_10us.cmap as it is used by smili as well --- src/smili/afmhot_10us.cmap | 1 + 1 file changed, 1 insertion(+) create mode 120000 src/smili/afmhot_10us.cmap diff --git a/src/smili/afmhot_10us.cmap b/src/smili/afmhot_10us.cmap new file mode 120000 index 0000000..69e2e51 --- /dev/null +++ b/src/smili/afmhot_10us.cmap @@ -0,0 +1 @@ +../difmap/afmhot_10us.cmap \ No newline at end of file