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
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
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
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
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)
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
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
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_complete.sh b/src/smili/run_complete.sh
old mode 100644
new mode 100755
diff --git a/src/smili/run_postprocessing.sh b/src/smili/run_postprocessing.sh
old mode 100644
new mode 100755
index 6cc9c9b..51ccb16
--- 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
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)
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("================================")