achael / eht-imaging

Imaging, analysis, and simulation software for radio interferometry

ehtim (eht-imaging)

Python modules for simulating and manipulating VLBI data and producing images with regularized maximum likelihood methods. This version is an early release so please raise an issue, submit a pull request, or email [email protected] if you have trouble or need help for your application.

The package contains several primary classes for loading, simulating, and manipulating VLBI data. The main classes are the Image, Movie, Array, Obsdata, Imager, and Caltable classes, which provide tools for loading images and data, producing simulated data from realistic u-v tracks, calibrating, inspecting, and plotting data, and producing images from data sets in various polarizations using various data terms and regularizing functions.

Installation

The latest stable version (1.2.8) is available on PyPi. Simply install pip and run

pip install ehtim

Incremental updates are developed on the dev branch. To use the very latest (unstable) code, checkout the dev branch, change to the main eht-imaging directory, and run:

pip install .

Installing with pip will update most of the required libraries automatically (numpy, scipy, matplotlib, astropy, ephem, future, h5py, and pandas).

If you want to use fast fourier transforms, you will also need to separately install NFFT and its pynfft wrapper. The simplest way is to use conda to install both:

conda install -c conda-forge pynfft

Alternatively, first install NFFT manually following the instructions on the readme, making sure to use the --enable-openmp flag in compilation. Then install pynfft, with pip, following the readme instructions to link the installation to where you installed NFFT. Finally, reinstall ehtim.

For M1 Macs (OS >= v12.0), install the M1 Mac version of pynfft and follow the instructions on the readme. It has the instructions to install fftw, nfft and then pynfft.

Certain eht-imaging functions require other external packages that are not automatically installed. In addition to pynfft, these include networkx (for image comparison functions), requests (for dynamical imaging), and scikit-image (for a few image analysis functions). However, the vast majority of the code will work without these dependencies.

Documentation and Tutorials

Documentation is here.

A intro to imaging tutorial jupyter notebook can be found in the repo at tutorials/ehtim_tutorial.ipynb

Slides for the included tutorial walk through the basic steps of reconstructing EHT images with the code

Here are some other ways to learn to use the code:

• Start with the script examples/example.py, which contains a series of sample commands to load an image and array, generate data, and produce an image with various imaging algorithms.
• Older Slides from the EHT2016 data generation and imaging workshop contain a tutorial on generating data with the VLBI imaging website, loading into the library, and producing an image.

Citation

If you use ehtim in your publication, please cite Chael+ 2018.

The latest version is also available as a static doi on Zenodo.

Selected publications that use ehtim

Let us know if you use ehtim in your publication and we'll list it here!

• High-Resolution Linear Polarimetric Imaging for the Event Horizon Telescope, Chael et al. 2016
• Computational Imaging for VLBI Image Reconstruction, Bouman et al. 2016
• Stochastic Optics: A Scattering Mitigation Framework for Radio Interferometric Imaging, Johnson 2016
• Reconstructing Video from Interferometric Measurements of Time-Varying Sources, Bouman et al. 2017
• Dynamical Imaging with Interferometry, Johnson et al. 2017
• Interferometric Imaging Directly with Closure Phases and Closure Amplitudes, Chael et al. 2018
• A Model for Anisotropic Interstellar Scattering and its Application to Sgr A*, Psaltis et al. 2018
• The Currrent Ability to Test Theories of Gravity with Black Hole Shadows, Mizuno et al. 2018
• The Scattering and Intrinsic Structure of Sagittarius A* at Radio Wavelengths, Johnson et al. 2018
• Testing GR with the Black Hole Shadow Size and Asymmetry of Sagittarius A*: Limitations from Interstellar Scattering, Zhu et al. 2018
• The Size, Shape, and Scattering of Sagittarius A* at 86 GHz: First VLBI with ALMA, Issaoun et al. 2019a
• First M87 Event Horizon Telescope Results IV: Imaging the Central Supermassive Black Hole, EHTC et al. 2019
• EHT-HOPS Pipeline for Millimeter VLBI Data Reduction, Blackburn et al. 2019
• Discriminating Accretion States via Rotational Symmetry in Simulated Polarimetric Images of M87, Palumbo et al. 2020
• SYMBA: An end-to-end VLBI synthetic data generation pipeline, Roelofs et al. 2020
• Monitoring the Morphology of M87* in 2009-2017 with the Event Horizon Telescope, Wielgus et al. 2020
• EHT imaging of the archetypal blazar 3C 279 at extreme 20 microarcsecond resolution, Kim et al. 2020
• Verification of Radiative Transfer Schemes for the EHT, Gold et al. 2020
• Closure Traces: Novel Calibration-insensitive Quantities for Radio Astronomy, Broderick and Pesce. 2020
• Evaluation of New Submillimeter VLBI Sites for the Event Horizon Telescope, Raymond et al. 2021
• Imaging VGOS Observations and Investigating Source Structure Effects, Xu et al. 2021
• A D-term Modeling Code (DMC) for Simultaneous Calibration and Full-Stokes Imaging of VLBI Data, Pesce et al. 2021
• Using space-VLBI to probe gravity around Sgr A*, Fromm et al. 2021
• Persistent Non-Gaussian Structure in the Image of Sagittarius A* at 86 GHz, Issaoun et al. 2021
• First M87 Event Horizon Telescope Results. VII. Polarization of the Ring, EHTC et al. 2021
• Event Horizon Telescope observations of the jet launching and collimation in Centaurus A, Janssen et al. 2021
• RadioAstron discovers a mini-cocoon around the restarted parsec-scale jet in 3C 84 Savolainen et al. 2021
• Unravelling the Innermost Jet Structure of OJ 287 with the First GMVA+ALMA Observations, Zhao et al. 2022
• First Sagittarius A* Event Horizon Telescope Results. III: Imaging of the Galactic Center Supermassive Black Hole, EHTC et al. 2022
• Resolving the Inner Parsec of the Blazar J1924-2914 with the Event Horizon Telescope, Issaoun et al. 2022
• The Event Horizon Telescope Image of the Quasar NRAO 530, Jorstad et al. 2023
• First M87 Event Horizon Telescope Results. IX. Detection of Near-horizon Circular Polarization, EHTC et al. 2023
• Filamentary structures as the origin of blazar jet radio variability, Fuentes et al. 2023
• The persistent shadow of the supermassive black hole of M 87. I. Observations, calibration, imaging, and analysis, EHTC 2023
• Parsec-scale evolution of the gigahertz-peaked spectrum quasar PKS 0858-279, Kosogorov et al. 2024

oifits Documentation

The oifits_new.py file used for reading/writing .oifits files is a slightly modified version of Paul Boley's package oifits. The oifits read/write functionality in ehtim is still being developed and may not work with all versions of python or astropy.

The documentation is styled after dfm's projects

License

ehtim is licensed under GPLv3. See LICENSE.txt for more details.

eht-imaging's Issues

Dockerfile fails to build

Describe the bug
Building the docker file fails with error

To Reproduce
Steps to reproduce the behavior:

1. Go to project root
2. Run docker build -t eht .
3. Scroll down to step [5/5] RUN pip install -r requirements.txt && conda install -y -c conda-forge pynfft && echo "backend: Agg" >> /opt/conda/lib/python2.7/site-packages/matplotlib/mpl-data/matplotlibrc
4. See error

Expected behavior
The docker build should finish successfully

Screenshots

Desktop (please complete the following information):

• OS: MacOS Big Sur
• Docker Engine Version v19.03.13

Warnings in ehtim_tutorial.ipynb

There are 4 warnings here:
https://github.com/achael/eht-imaging/blob/main/tutorials/ehtim_tutorial.ipynb
colored in red. The first one:

No NX table in uvfits!
/home/achael/anaconda3/lib/python3.8/site-packages/erfa/core.py:154: ErfaWarning: ERFA function "d2dtf" yielded 212584 of "dubious year (Note 5)"
warnings.warn('ERFA function "{}" yielded {}'.format(func_name, wmsg),

Mismatched dimensions

Hi,

I am attempting to make some models for the international baselines of LOFAR using the EHT imager but I have stumbled on an error with the imaging function I do not quite understand.

This is how I call it:

eh.imager_func( obs, gaussprior, gaussprior, zbl, d1='cphase', d2='camp', s1='gs', s2='gs', alpha_d1=50, alpha_d2=50, clipfloor=0.001, maxit=300, stop=0.0001 )

Obs is generated from a measurement set (which was turned into uvfits).

Sorry if this is obscure but it looks like a dimensional issue. Let me know if you need extra information. Hopefully this is a simple issue and a mistake on my behalf.

Thanks!

Joe

Getting bispectra:: type vis, count min, scan 3370/3370 

/home/callingham/.local/lib/python2.7/site-packages/ehtim/imaging/imager_utils.py:2251: RuntimeWarning: invalid value encountered in greater
  snrmask = np.abs(clphasearr['cphase']/clphasearr['sigmacp']) > snrcut
Traceback (most recent call last):
  File "eht_imaging_joe.py", line 371, in <module>
    main( args.vis, closure_tels = args.ctels, npix = args.npix, fov_arcsec = args.fov, maxarcmin = args.maxarcmin, zbl = args.zbl, prior_fwhm_arcsec = args.prior_fwhm, doplots = args.doplots, imfile = args.imfile, remove_tels = args.rtels, niter=args.maxiter, cfloor=args.clipfloor, conv_criteria=args.convg, use_bs=args.use_bs, scratch_model=args.scratch_model, do_diagnostic_plots=args.diagplots, lotss_file=args.lotss_file )
  File "eht_imaging_joe.py", line 320, in main
    out = eh.imager_func( obs, gaussprior, gaussprior, zbl, d1='cphase', d2='camp', s1='gs', s2='gs', alpha_d1=50, alpha_d2=50, clipfloor=cfloor, maxit=niter, stop=conv_criteria )
  File "/home/callingham/.local/lib/python2.7/site-packages/ehtim/imaging/imager_utils.py", line 190, in imager_func
    (data1, sigma1, A1) = chisqdata(Obsdata, Prior, embed_mask, d1, **kwargs)
  File "/home/callingham/.local/lib/python2.7/site-packages/ehtim/imaging/imager_utils.py", line 625, in chisqdata
    (data, sigma, A) = chisqdata_cphase(Obsdata, Prior, mask, **kwargs)
  File "/home/callingham/.local/lib/python2.7/site-packages/ehtim/imaging/imager_utils.py", line 2259, in chisqdata_cphase
    sigma = np.linalg.norm([sigma, systematic_cphase_noise*np.ones(len(sigma))], axis=0)[snrmask]
IndexError: boolean index did not match indexed array along dimension 0; dimension is 70628 but corresponding boolean dimension is 70631
Exception TypeError: TypeError("'NoneType' object is not callable",) in <bound method UmfpackContext.new_del of <scipy.sparse.linalg.dsolve.umfpack.umfpack.UmfpackContext object at 0x4f39910>> ignored

Provide 'extras' for installing optional packages

Looking at https://github.com/achael/eht-imaging/blob/master/setup.py I noticed 4 dependencies commented out due to them only being necessary for optional things. This sort of thing is a perfect use for "extras" installs. It lets users specify the "extras" they want and then only get bonus features if they want them.

No NFFT Installed issue

Describe the bug
Program says that no NFFT is installed after following installation instructions in the README.

To Reproduce
Steps to reproduce the behavior:
1.Install NFFT and pynfft with conda as described in readme
2. Clone repository
3. Build with dockerfile
3. Run dockerfile
4. Run first 35 lines of examples/examples.py in ipython

Expected behavior
Running line 35 causes an error. I ran on both linux and MacOS; in the first case I got an issue that said no NFFT is installed, in the second I got an error that says 'NFFT is not defined.'

Screenshots
If applicable, add screenshots to help explain your problem.

Desktop (please complete the following information):
macOS Catalina
linux centOS Fedora

UVFITS-files of VLA-data exported by CASA cannot be successfully read

Hello Andrew,

I am hoping to use the EHT imaging library to image some data from the VLA, as it might allow me to automate some steps in the process (I have a LOT of data). In this data, Sgr A was observed and both the compact component and Sgr A West are visible. To test this, I exported a single scan at a single frequency (34.256 GHz, 125 MHz BW) from CASA (initially using version 4.7.0) using 'exportuvfits', and I tried loading it into ehtim with 'obsdata.load_uvfits' (using ehtimaging 1.0.0). Unfortunately, I got the following error message:

UnicodeDecodeError: 'ascii' codec can't decode byte 0x90 in position 5: ordinal not in range(128)

I have tried changing the codec to 'utf8', which resulted in a different error message:

UnicodeDecodeError: 'utf8' codec can't decode byte 0x90 in position 5: invalid start byte

I then tried using CASA 5.3.0 to export a UVFITS file, as some things in the export function have changed between 4.7.0 and 5.3.0. With the new UVFITS file, I got a new error message when attempting to import it:

IndexError: index 26 is out of bounds for axis 0 with size 26

This looks like it perhaps has something to do with the antenna table, as I have 26 antennas in the array. Maybe a flagging-related issue? I also tried exporting a UVFITS file from CASA with the setting 'padwithflags=True', but this does not result in any observable difference.

For convenience, I have uploaded a zip file here with the four UVFITS files I have tried so far, using the two CASA versions I have available and various combinations of export settings. All of these files work just fine with Difmap, and the visibilities I see are sensible. Do you know what might be wrong, and how to fix it?

Thank you!

UVFITSfiles.zip

Increase the usage of augmented assignment statements

👀 Some source code analysis tools can help to find opportunities for improving software components.
💭 I propose to increase the usage of augmented assignment statements accordingly.

Would you like to integrate anything from a transformation result which can be generated by a command like the following?
(:point_right: Please check also for questionable change suggestions because of an evolving search pattern.)

lokal$ perl -p -i.orig -0777 -e 's/^(?<indentation>\s+)(?<target>\S+)\s*=\s*\k<target>[ \t]*(?<operator>[+%^@]|-(?!>)|&(?!&)|\|(?!\|)|\*\*?|\/\/?|<<|>>)/$+{indentation}$+{target} $+{operator}=/gm' $(find ~/Projekte/eht-imaging/lokal -name '*.py')

Math behind this project

Hi!
I just discovered this project, it's awesome! Where can I read the math research and all related to this project?

Zenado link not working

Just reporting - this link:
https://zenodo.org/badge/latestdoi/42943499[Zenodo](https://zenodo.org/badge/latestdoi/42943499Zenodo)
is not working.

Please clarify if the license is gpl 3 or later vs gpl 3 only

Please update README.rst at least.

To understand why this is important, see
https://www.gnu.org/licenses/identify-licenses-clearly.en.html

For instructions on how to properly apply the license, read the end of the gpl license.

zero-size array when running imager_func

I'm calling the imager_func like this:
zbl = 0.28
cfloor = 0.0028
niter = 300
conv_criteria = 0.0001
out = eh.imager_func( obs, gaussprior, gaussprior, zbl, d1='cphase', d2='camp', s1='gs', s2='gs', alpha_d1=50, alpha_d2=50, clipfloor=cfloor, maxit=niter, stop=conv_criteria )

The gaussprior is a single gaussian.

I get the following error and traceback:

/home/lmorabit/.local/lib/python2.7/site-packages/ehtim/imaging/imager_utils.pyc in imager_func(Obsdata, InitIm, Prior, flux, d1, d2, d3, alpha_d1, alpha_d2, alpha_d3, s1, s2, s3, alpha_s1, alpha_s2, alpha_s3, alpha_flux, alpha_cm, **kwargs)
294 # Print stats
295 print("Initial S_1: %f S_2: %f S_3: %f" % (reg1(ninit), reg2(ninit), reg3(ninit)))
--> 296 print("Initial Chi^2_1: %f Chi^2_2: %f Chi^2_3: %f" % (chisq1(ninit), chisq2(ninit), chisq3(ninit)))
297 print("Initial Objective Function: %f" % (objfunc(xinit)))
298

/home/lmorabit/.local/lib/python2.7/site-packages/ehtim/imaging/imager_utils.pyc in chisq2(imvec)
199
200 def chisq2(imvec):
--> 201 return chisq(imvec, A2, data2, sigma2, d2, ttype=ttype, mask=embed_mask)
202
203 def chisq2grad(imvec):

/home/lmorabit/.local/lib/python2.7/site-packages/ehtim/imaging/imager_utils.pyc in chisq(imvec, A, data, sigma, dtype, ttype, mask)
373 chisq = chisq_cphase(imvec, A, data, sigma)
374 elif dtype == 'camp':
--> 375 chisq = chisq_camp(imvec, A, data, sigma)
376 elif dtype == 'logcamp':
377 chisq = chisq_logcamp(imvec, A, data, sigma)

/home/lmorabit/.local/lib/python2.7/site-packages/ehtim/imaging/imager_utils.pyc in chisq_camp(imvec, Amatrices, clamp, sigma)
739 def chisq_camp(imvec, Amatrices, clamp, sigma):
740 """Closure Amplitudes (normalized) chi-squared"""
--> 741 print( np.min(imvec) )
742 print( np.min(Amatrices[1] ) )
743 clamp_samples = np.abs(np.dot(Amatrices[0], imvec) * np.dot(Amatrices[1], imvec) / (np.dot(Amatrices[2], imvec) * np.dot(Amatrices[3], imvec)))

/usr/lib64/python2.7/site-packages/numpy/core/fromnumeric.pyc in amin(a, axis, out, keepdims)
2370
2371 return _methods._amin(a, axis=axis,
-> 2372 out=out, **kwargs)
2373
2374

/usr/lib64/python2.7/site-packages/numpy/core/_methods.pyc in _amin(a, axis, out, keepdims)
27
28 def _amin(a, axis=None, out=None, keepdims=False):
---> 29 return umr_minimum(a, axis, None, out, keepdims)
30
31 def _sum(a, axis=None, dtype=None, out=None, keepdims=False):

ValueError: zero-size array to reduction operation minimum which has no identity

ValueError: numpy.dtype size changed

Numpy recently released version 2.0, and I am struggling to get fresh installs working without:

ValueError: numpy.dtype size changed, may indicate binary incompatibility. Expected 96 from C header, got 88 from PyObject

I don't think this is a personal problem, so if anyone else sees this, please chime in.

pyNFFT for M1 Mac (M1 Pro, M1 Max etc.)

conda install -c conda-forge pynfft will not work for newer M1 Mac / M1 Pro / M1 Max

I have forked the pyNFFT git repository here https://github.com/rohandahale/pyNFFT.git and I have written down step-by-step instructions for installing pyNFFT on M1 Macs. I have also made small edit to the setup.py file to use gcc-12 as the C compiler which works the best for installation with openmp.

It could be helpful if this repository link is added as a small note for pynfft installation in the README.md.

Congratulations on the big achievement!

Sorry to open this as an issue, but I really wanted to say congrats to your whole team who made it possible to achieve something truly incredible. I really wish people could admire work like this a billion times than all the celebrities and stuff they usually admire. Brilliant work and kudos to your complete team!

@chanchikwan @lindyblackburn @klbouman @danielpalumbo @jrfarah

:DDDDDDD

Conda fails to install `skyfield` or `paramsurvey[ray]` during project configuration

Describe the bug
See title.

To Reproduce
Open up the project in pycharm and let it try to install all dependencies, then open up requirements.txt

Expected behavior
Both packages should be installed properly.

Desktop (please complete the following information):

• OS: Debian 11 Bullseye

Remediation
In a terminal in the project folder, running pip install skyfield paramsurvey successfully installs the packages.

Additionally, it seems that the packages install by pip are newer than the ones in the conda repos

Lastest numpy 1.24.2 does not work with ehtim

pip install ehtim will install numpy 1.24.2. In a python 3.8 shell when you do import ehtim, you will get an error:

ImportError: cannot import name 'bool' from 'numpy'

The correct way to install ehtim should be:

pip install ehtim numpy==1.21.1

Observations cannot generate tlists(), possibly due to ragged arrays being abandoned in numpy 1.24.

Hi,

After loading an observation, calling obs.tlist() throws the error "ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (72,) + inhomogeneous part." where the 72 is case-specific. Is there a recommended workaround or should I backdate numpy?

Closure phases locked at zero and clip setting

Hi again Andrew,

After that last error, the script is running but there is a couple of features I do not quite understand. For one, the chi2_1 (representing phase) is set to zero and does not change at all in the convergence processes. I do not know why that is the case, the output is below and my call to imag_func.

Also, the script outputs that the Total Pixel #: 16384 and Clipped Pixel #: 16384. That does not make a lot of sense to me, as how can all the pixels be clipped?

Sorry again if I am missing something simple.

Joe

`
Using closure phase and closure amplitude separately.
Warning! Specified flux is > 120% of maximum visibility amplitude!
Getting bispectra:: type vis, count min, scan 3371/3371

Getting closure amps:: type vis camp , count min, scan 3371/3371

Initial S_1: -0.000000 S_2: -0.000000 S_3: 0.000000
Initial Chi^2_1: 0.000000 Chi^2_2: 1244289.391793 Chi^2_3: 1.000000
Initial Objective Function: 326493584.274608
Total Data 1: 70791
Total Data 2: 67420
Total Pixel #: 16384
Clipped Pixel #: 16384

i: 0 chi2_1: 0.00 chi2_2: 1244289.39 chi2_3: 1.00 s_1: 0.00 s_2: 0.00 s_3: 0.00
i: 1 chi2_1: 0.00 chi2_2: 126347.43 chi2_3: 1.00 s_1: 0.03 s_2: 0.03 s_3: 0.00
i: 2 chi2_1: 0.00 chi2_2: 107404.42 chi2_3: 1.00 s_1: 0.03 s_2: 0.03 s_3: 0.00
i: 3 chi2_1: 0.00 chi2_2: 95902.42 chi2_3: 1.00 s_1: 0.03 s_2: 0.03 s_3: 0.00
i: 4 chi2_1: 0.00 chi2_2: 51008.25 chi2_3: 1.00 s_1: 0.04 s_2: 0.04 s_3: 0.00
i: 5 chi2_1: 0.00 chi2_2: 30726.60 chi2_3: 1.00 s_1: 0.05 s_2: 0.05 s_3: 0.00
i: 6 chi2_1: 0.00 chi2_2: 4324.99 chi2_3: 1.00 s_1: 0.31 s_2: 0.31 s_3: 0.00
i: 7 chi2_1: 0.00 chi2_2: 511.05 chi2_3: 1.00 s_1: 5.37 s_2: 5.37 s_3: 0.00
i: 8 chi2_1: 0.00 chi2_2: 774.29 chi2_3: 1.00 s_1: 80.65 s_2: 80.65 s_3: 0.00
i: 9 chi2_1: 0.00 chi2_2: 5.80 chi2_3: 1.00 s_1: 285.23 s_2: 285.23 s_3: 0.00
i: 10 chi2_1: 0.00 chi2_2: 27.71 chi2_3: 1.00 s_1: 292.61 s_2: 292.61 s_3: 0.00
i: 11 chi2_1: 0.00 chi2_2: 49.97 chi2_3: 1.00 s_1: 350.60 s_2: 350.60 s_3: 0.00
i: 12 chi2_1: 0.00 chi2_2: 43.36 chi2_3: 1.00 s_1: 413.17 s_2: 413.17 s_3: 0.00
i: 13 chi2_1: 0.00 chi2_2: 26.81 chi2_3: 1.00 s_1: 495.87 s_2: 495.87 s_3: 0.00
i: 14 chi2_1: 0.00 chi2_2: 28.83 chi2_3: 1.00 s_1: 533.62 s_2: 533.62 s_3: 0.00
i: 15 chi2_1: 0.00 chi2_2: 15.65 chi2_3: 1.00 s_1: 577.21 s_2: 577.21 s_3: 0.00
i: 16 chi2_1: 0.00 chi2_2: 16.63 chi2_3: 1.00 s_1: 608.58 s_2: 608.58 s_3: 0.00
i: 17 chi2_1: 0.00 chi2_2: 9.78 chi2_3: 1.00 s_1: 615.92 s_2: 615.92 s_3: 0.00
i: 18 chi2_1: 0.00 chi2_2: 8.06 chi2_3: 1.00 s_1: 637.82 s_2: 637.82 s_3: 0.00
i: 19 chi2_1: 0.00 chi2_2: 9.04 chi2_3: 1.00 s_1: 641.03 s_2: 641.03 s_3: 0.00
i: 20 chi2_1: 0.00 chi2_2: 11.15 chi2_3: 1.00 s_1: 677.38 s_2: 677.38 s_3: 0.00
i: 21 chi2_1: 0.00 chi2_2: 0.55 chi2_3: 1.00 s_1: 691.50 s_2: 691.50 s_3: 0.00
i: 22 chi2_1: 0.00 chi2_2: 6.49 chi2_3: 1.00 s_1: 694.04 s_2: 694.04 s_3: 0.00
i: 23 chi2_1: 0.00 chi2_2: 8.52 chi2_3: 1.00 s_1: 695.83 s_2: 695.83 s_3: 0.00
i: 24 chi2_1: 0.00 chi2_2: 2.66 chi2_3: 1.00 s_1: 701.40 s_2: 701.40 s_3: 0.00
i: 25 chi2_1: 0.00 chi2_2: 3.82 chi2_3: 1.00 s_1: 706.82 s_2: 706.82 s_3: 0.00
i: 26 chi2_1: 0.00 chi2_2: 1.32 chi2_3: 1.00 s_1: 710.09 s_2: 710.09 s_3: 0.00
i: 27 chi2_1: 0.00 chi2_2: 1.79 chi2_3: 1.00 s_1: 712.76 s_2: 712.76 s_3: 0.00
i: 28 chi2_1: 0.00 chi2_2: 1.21 chi2_3: 1.00 s_1: 714.99 s_2: 714.99 s_3: 0.00
i: 29 chi2_1: 0.00 chi2_2: 1.07 chi2_3: 1.00 s_1: 714.74 s_2: 714.74 s_3: 0.00
i: 30 chi2_1: 0.00 chi2_2: 0.90 chi2_3: 1.00 s_1: 714.46 s_2: 714.46 s_3: 0.00
i: 31 chi2_1: 0.00 chi2_2: 0.65 chi2_3: 1.00 s_1: 714.57 s_2: 714.57 s_3: 0.00
i: 32 chi2_1: 0.00 chi2_2: 0.42 chi2_3: 1.00 s_1: 715.06 s_2: 715.06 s_3: 0.00
i: 33 chi2_1: 0.00 chi2_2: 0.34 chi2_3: 1.00 s_1: 715.62 s_2: 715.62 s_3: 0.00
i: 34 chi2_1: 0.00 chi2_2: 0.26 chi2_3: 1.00 s_1: 716.15 s_2: 716.15 s_3: 0.00
i: 35 chi2_1: 0.00 chi2_2: 0.25 chi2_3: 1.00 s_1: 716.34 s_2: 716.34 s_3: 0.00
i: 36 chi2_1: 0.00 chi2_2: 0.23 chi2_3: 1.00 s_1: 716.70 s_2: 716.70 s_3: 0.00
i: 37 chi2_1: 0.00 chi2_2: 0.21 chi2_3: 1.00 s_1: 717.02 s_2: 717.02 s_3: 0.00
i: 38 chi2_1: 0.00 chi2_2: 0.20 chi2_3: 1.00 s_1: 717.20 s_2: 717.20 s_3: 0.00
time: 464.397219 s
J: 1344.529552
Final Chi^2_1: 0.000000 Chi^2_2: 0.202248 Chi^2_3: 1.000000
CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
Warning! Specified flux is > 120% of maximum visibility amplitude!
Getting bispectra:: type vis, count min, scan 3371/3371

Getting closure amps:: type vis camp , count min, scan 3371/3371

Initial S_1: -0.000000 S_2: -0.000000 S_3: 0.000000
Initial Chi^2_1: 0.000000 Chi^2_2: 0.290918 Chi^2_3: 1.000000
Initial Objective Function: -85.130163
Total Data 1: 70791
Total Data 2: 67420
Total Pixel #: 16384
Clipped Pixel #: 16384

i: 0 chi2_1: 0.00 chi2_2: 0.29 chi2_3: 1.00 s_1: 0.00 s_2: 0.00 s_3: 0.00
i: 1 chi2_1: 0.00 chi2_2: 0.28 chi2_3: 1.00 s_1: 0.00 s_2: 0.00 s_3: 0.00
i: 2 chi2_1: 0.00 chi2_2: 0.28 chi2_3: 1.00 s_1: 0.00 s_2: 0.00 s_3: 0.00
i: 3 chi2_1: 0.00 chi2_2: 0.25 chi2_3: 1.00 s_1: 0.02 s_2: 0.02 s_3: 0.00
i: 4 chi2_1: 0.00 chi2_2: 0.25 chi2_3: 1.00 s_1: 0.02 s_2: 0.02 s_3: 0.00
i: 5 chi2_1: 0.00 chi2_2: 0.24 chi2_3: 1.00 s_1: 0.05 s_2: 0.05 s_3: 0.00
i: 6 chi2_1: 0.00 chi2_2: 0.22 chi2_3: 1.00 s_1: 0.09 s_2: 0.09 s_3: 0.00
i: 7 chi2_1: 0.00 chi2_2: 0.21 chi2_3: 1.00 s_1: 0.19 s_2: 0.19 s_3: 0.00
i: 8 chi2_1: 0.00 chi2_2: 0.20 chi2_3: 1.00 s_1: 0.19 s_2: 0.19 s_3: 0.00
i: 9 chi2_1: 0.00 chi2_2: 0.19 chi2_3: 1.00 s_1: 0.28 s_2: 0.28 s_3: 0.00
i: 10 chi2_1: 0.00 chi2_2: 0.19 chi2_3: 1.00 s_1: 0.35 s_2: 0.35 s_3: 0.00
i: 11 chi2_1: 0.00 chi2_2: 0.18 chi2_3: 1.00 s_1: 0.43 s_2: 0.43 s_3: 0.00
i: 12 chi2_1: 0.00 chi2_2: 0.18 chi2_3: 1.00 s_1: 0.48 s_2: 0.48 s_3: 0.00
i: 13 chi2_1: 0.00 chi2_2: 0.17 chi2_3: 1.00 s_1: 0.56 s_2: 0.56 s_3: 0.00
i: 14 chi2_1: 0.00 chi2_2: 0.16 chi2_3: 1.00 s_1: 0.65 s_2: 0.65 s_3: 0.00
i: 15 chi2_1: 0.00 chi2_2: 0.16 chi2_3: 1.00 s_1: 0.65 s_2: 0.65 s_3: 0.00
i: 16 chi2_1: 0.00 chi2_2: 0.16 chi2_3: 1.00 s_1: 0.71 s_2: 0.71 s_3: 0.00
i: 17 chi2_1: 0.00 chi2_2: 0.16 chi2_3: 1.00 s_1: 0.72 s_2: 0.72 s_3: 0.00
i: 18 chi2_1: 0.00 chi2_2: 0.15 chi2_3: 1.00 s_1: 0.77 s_2: 0.77 s_3: 0.00
i: 19 chi2_1: 0.00 chi2_2: 0.15 chi2_3: 1.00 s_1: 0.77 s_2: 0.77 s_3: 0.00
i: 20 chi2_1: 0.00 chi2_2: 0.15 chi2_3: 1.00 s_1: 0.79 s_2: 0.79 s_3: 0.00
i: 21 chi2_1: 0.00 chi2_2: 0.15 chi2_3: 1.00 s_1: 0.81 s_2: 0.81 s_3: 0.00
i: 22 chi2_1: 0.00 chi2_2: 0.15 chi2_3: 1.00 s_1: 0.80 s_2: 0.80 s_3: 0.00
i: 23 chi2_1: 0.00 chi2_2: 0.15 chi2_3: 1.00 s_1: 0.82 s_2: 0.82 s_3: 0.00
i: 24 chi2_1: 0.00 chi2_2: 0.15 chi2_3: 1.00 s_1: 0.82 s_2: 0.82 s_3: 0.00
i: 25 chi2_1: 0.00 chi2_2: 0.15 chi2_3: 1.00 s_1: 0.82 s_2: 0.82 s_3: 0.00
time: 359.526318 s
J: -90.942385
Final Chi^2_1: 0.000000 Chi^2_2: 0.148237 Chi^2_3: 1.000000
CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
done.
Exception TypeError: TypeError("'NoneType' object is not callable",) in <bound method UmfpackContext.new_del of <scipy.sparse.linalg.dsolve.umfpack.umfpack.UmfpackContext object at 0x44a2910>> ignored`

How to load all IF's and channels without averaging?

Is it possible to load uvfits data without averaging over frequency (IF, channels) with obsdata.load_uvfits()?
My (deleted) example had a different problem: I assumed the data had more than one channels.

Anyway, it seems that obsdata.load_uvfits() cannot load multiple IF's and channels?

Remove requirement for jupyter?

This line seems not to be necessary - the import isn't used and I successfully used the library without it. If removed, eht-imaging will no longer require the installation of the jupyter libraries, which would be nice.

eht-imaging/ehtim/modeling/modeling_utils.py

Line 49 in bbd3aeb

from IPython import display

Missing `*` in eht-imaging/ehtim /image.py

In line 3503 of the file in the title an asterisk is missing. It should be:
factor *= beamarea / (self.psize**2).
Without the asterisk, setting cbar_unit equal to either ['mJy', 'beam'] or ['Jy', 'beam'] yields the same result.

reorder_baselines loses data

The reorder_baselines method of the Obsdata class assumes that any observations obtained on telescopes (t1, t2) and (t2, t1) are conjugates of each other, and throws away one of the two. However, this is not generally they case, in particular if the observations were obtained at different times and/or frequencies.

I discovered this when importing an OIFITS file, only to find that of the 60 or so wavelengths in the file, only a single one remains after reorder_baselines is run.

Merge oifits changes/fixes

Just to let you know, that the (modified?) oifits module used here seems to be an old version. The module is also on github, which is the authoritative source (rather than the old webpage at urfu.ru): https://github.com/pboley/oifits
It might be a good idea to synchronize changes between the two.

avoid numpy deprecation messages

Describe the bug
With the current numpy version, using obsdata.tlist() or obsdata.bllist() returns an ugly deprecation warning.

This deprecation is in numpy 1.19.0 released June 2020
https://numpy.org/doc/stable/release/1.19.0-notes.html#deprecate-automatic-dtype-object-for-ragged-input

To Reproduce

Python 3.8.10 (default, Nov 26 2021, 20:14:08) 
[GCC 9.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> np.array([[1, 2], [3, 4, 5]])
<stdin>:1: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
array([list([1, 2]), list([3, 4, 5])], dtype=object)

Note that the returned array has dtype=object, so adding it does not change the behavior while silencing the warning:

>>> np.array([[1, 2], [3, 4, 5]], dtype=object)
array([list([1, 2]), list([3, 4, 5])], dtype=object)

The patch in the pullreq is two places that are obviously creating ragged arrays, and a comment at a potential 3rd place.

`switch_polrep` gives incorrect results when feeds are missing

I just found a small but silent bug when using eht-imaging to read in polarized data when feeds are missing. Suppose you are reading in data where the R CP feed is missing. If you read in the data using

obs = ehtim.obsdata.load_uvfits(file, polrep="circ")

then everything is fine, and I see the R products are NaN'd. However, if I read in the data using the method

obs = ehtim.obsdata.load_uvfits(file)
obsc = obs.switch_polrep("circ")

I get that the L products are now NaN'd out.

My guess is that the switch_polrep function doesn't know which components are NaN'd when moving from stokes to circular, so it just defaults to the L feeds. I'm unsure if that is the correct behavior because it can lead to some wild d-terms if you get the polarization hands wrong.

`afmhot` is Not A Perceptually Uniform Colormap

Cool picture, guys! But it seems like you've chosen a non-perceptually uniform colormap, which means that humans are likely to slightly misinterpret the image. It'd be great to see this historic image in something more human-comprehensible, like viridis or plasma rather than afmhot!

`interp2d` is now deprecated in `scipy`

The function interp2d is used in ehtim for example in regrid_image, but it is now deprecated in scipy as of version 1.14 and has been removed. See the documentation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html

It is available but throws a deprecation warning as late as version 1.13. Currently ehtim requires scipy but does not specify a version. Alternatively, we could update the usage of interp2d to be compatible with the new scipy. This would involve using RectBivariateSpline or bisplrep/bisplev, per the documentation. This potentially comes with significant drawbacks, including being very slow compared to interp2d. See for example this StackOverflow discussion: https://stackoverflow.com/questions/75427538/regulargridinterpolator-excruciatingly-slow-compared-to-interp2d This is a known issue, as the scipy devs comment and say:

Commenting here for visibility. It is a known performance regression and we are actively working on it. The new interface was needed in term of maintainability moving forward.

Do folks have a preference for enforcing scipy=1.13 vs. updating interp2d to the new scipy functions?

Running example.py throws a NameError

I tried running example.py line by line in a Python terminal but it fails. It throws an error specifically on line 76 of example.py:

out = eh.imager_func(obs, gaussprior, gaussprior, flux, d1 = 'bs', s1 = 'simple', alpha_s1 = 1, alpha_d1 = 100, alpha_flux = 100, alpha_cm = 50, maxit = 100)

Here's the traceback:

Getting bispectra: type: vis count: min scan 107/107 Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/sean/eht/eht-imaging/ehtim/imaging/imager_utils.py", line 178, in imager_func
    (data1, sigma1, A1) = chisqdata(Obsdata, Prior, embed_mask, d1, **kwargs)#ttype=ttype, fft_pad_factor=fft_pad_factor, conv_func=grid_conv, p_rad=grid_prad, order=fft_interp, debias=debias,snrcut=snrcut)
  File "/home/sean/eht/eht-imaging/ehtim/imaging/imager_utils.py", line 583, in chisqdata
    (data, sigma, A) = chisqdata_bs(Obsdata, Prior, mask, **kwargs)
  File "/home/sean/eht/eht-imaging/ehtim/imaging/imager_utils.py", line 1884, in chisqdata_bs
    snrmask = np.abs(biarr['bispec']/data_arr['sigmab']) > snrcut
NameError: name 'data_arr' is not defined

Where should I define the global name data_arr for it to work?

Cannot merge_obs if obsdata have scans.

Describe the bug
The bug is caused by comparing an array to an empty list using ==.

To Reproduce
Here's a simple example using the public EHT data from 2017:

import ehtim as eh
obs = eh.obsdata.load_uvfits('/home/daniel/local_data/m87_data/zblcal+selfcal/hops_hi_3597_M87+zbl-dtcal_selfcal.uvfits')
obs2 = eh.obsdata.load_uvfits('/home/daniel/local_data/m87_data/zblcal+selfcal/hops_hi_3598_M87+zbl-dtcal_selfcal.uvfits')
eh.obsdata.merge_obs([obs,obs2])

Expected behavior
A new merged obsdata.

Error
in merge_obs(obs_List, force_merge)
4731 mjd_offset = obs.mjd - mjd_ref
4732 obs.data['time'] += mjd_offset * 24
-> 4733 if not(obs.scans is None or obs.scans == []):
4734 obs.scans += mjd_offset * 24
4736 # merge the data

ValueError: operands could not be broadcast together with shapes (18,2) (0,)

installation with NFFT on a cluster as a non-sys admin

Step 1: NFFT
git clone https://github.com/NFFT/nfft
Following the installation instructions provided, there are a couple things to watch out for. You need to run configure with openMP enabled, and if you're installing on a cluster then you probably want to do a local installation. This should do the trick:
./configure --enable-all --enable-openmp --prefix=/absolute/path/to/nfft

Step 2: pyNFFT
If you don't have NFFT installed system wide, then you have to run setup.py like this:
ttfamily python setup.py build_ext -I /absolute/path/to/nfft/include -L /absolute/path/to/nfftl/lib -R /absolute/path/to/nfft/lib

However, there is an issue which you have to fix ( https://github.com/ghisvail/pyNFFT/issues/56 ). To do this open pynfft/util.c and either comment out or remove the line #include "nfft3util.h"
You'll probably also have to use the --inplace option while running python setup.py build_ext and then the --user option while running python setup.py install unless you're a system administrator. Note that if the install fails because it can't find the site-packages directory in your installation directory, just add it to your $PYTHONPATH and run the install again.

Step 3: eht-imaging
git clone https://github.com/achael/eht-imaging cd eht-imaging pip install . --user

NB: if you get an error about an unexpected keyword "overwrite" then you need to modify the ehtim/io/save.py script to update hdulist.writeto(fname, overwrite=True) to hdulist.writeto(fname, clobber=True))