HERA-Team / hera_sim

Simple simulation code for HERA-like redundant interferometric arrays
Other
16 stars 8 forks source link

VisibilitySimulator incorrectly errors when passed some obsparam files #87

Closed Jackmastr closed 3 years ago

Jackmastr commented 4 years ago

Description

A VisibilitySimulator currently cannot be tested against all of the pyuvsim reference simulations because it is are too restrictive in the format for its obsparam parameter. It should assume that there is zero flux at a given frequency if there are no sources that emit at that frequency. Instead it throws a ValueError.

To Reproduce

Pass a VisibilitySimulator such as VisCPU an obsparam file with Nfreq > 1 and a catalog that only lists sources that emit at 1 frequency. pyuvsim's reference_simulations/obsparam_ref_1.3_gauss.yaml file is one such example with Nfreq = 6250.

Files

simulators.py

lines 91-94
if obsparams:
    (self.uvdata,
     self.beams,
     self.beam_ids) = initialize_uvdata_from_params(obsparams)
lines 101-103
catalog = initialize_catalog_from_params(obsparams)[0]
point_source_pos = np.array([catalog['ra_j2000'], catalog['dec_j2000']]).T * np.pi/180.
point_source_flux = np.atleast_2d(catalog['flux_density_I'])
lines 136-137
if sky_freqs is None:
        self.sky_freqs = np.unique(self.uvdata.freq_array)
lines 164-167
if self.point_source_flux is not None:
        if self.point_source_flux.shape[0] != self.sky_freqs.shape[0]:
            raise ValueError("point_source_flux must have the same number "
                             "of freqs as sky_freqs.")

pyuvsim: reference_simulations/obsparam_ref_1.3_gauss.yaml

filing:
  outdir: "./simulation_results/"
  outfile_name: 'ref_1.3_gauss'
  output_format: 'uvh5'
freq:
        #  Nfreqs:      6144
        #  channel_width: 32552.0
        #  start_freq:  50000000.0
        #  end_freq:   250000000.0
  Nfreqs:      6250.
  channel_width: 32000.0
  start_freq:  50016000.0
  end_freq:   249984000.0
sources:
  catalog: "catalog_files/letter_R_12pt_2458098.38824015.txt"
telescope:
  array_layout: telescope_config/baseline_lite.csv
  telescope_config_name: telescope_config/bl_lite_gauss.yaml
time:
  Ntimes: 2
  integration_time: 11.0
  start_time: 2458098.38824015

pyuvsim: reference_simulations/catalog_files/letter_R_12pt_2458098.38824015.txt

SOURCE_ID   RA_J2000 [deg]  Dec_J2000 [deg] Flux [Jy]   Frequency [Hz]
HERATEST5   59.37045    -28.778843  1   100000000.0
HERATEST6   57.08925    -28.74223   1   100000000.0
HERATEST12  59.38125    -27.778828  1   100000000.0
HERATEST13  58.25100    -27.765359  1   100000000.0
HERATEST21  59.39115    -26.779049  1   100000000.0
HERATEST22  58.27125    -26.765736  1   100000000.0
HERATEST23  57.15150    -26.743624  1   100000000.0
HERATEST30  59.40120    -25.779269  1   100000000.0
HERATEST31  57.18090    -25.744495  1   100000000.0
HERATEST39  59.41035    -24.779242  1   100000000.0
HERATEST40  58.30965    -24.766704  1   100000000.0
HERATEST41  57.20820    -24.744905  1   100000000.0
steven-murray commented 4 years ago

Thanks @Jackmastr for the nice detailed report.

@aelanman: what is the intended behaviour of the input catalog above? Is the idea that each source has flat spectral index?

aelanman commented 4 years ago

@steven-murray Yes. The reference simulation catalogs are all flat-spectrum.

Jackmastr commented 4 years ago

@aelanman should this always be the assumption though? Why does the catalog mention one frequency if we assume a constant flux across all frequencies?

aelanman commented 4 years ago

Recent updates to pyuvsim do support non-flat spectra, and we're working on a new set of reference simulations that take advantage of this. The frequency mentioned in the catalog was a placeholder for future frequency support, which ultimately never got used.

steven-murray commented 4 years ago

I'm trying to think what the best way forward here is. Should hera_sim assume that if the catalog file defines only a single frequency, that it is flat-spectrum, and extend the frequencies according to the obsparam file?

aelanman commented 4 years ago

The freq column in catalogs is deprecated. It was never used in the first place, and with the addition of full-frequency HEALPix maps, we needed multiple reference frequencies for each source anyway. What we're moving toward in pyradiosky is to have a freq_array attribute, with the same length as the fluxes for each source. So far (to my knowledge), only two spectral definitions are supported:

We're adding support for spectral indices, interpolating between subbands, etc.

I think this further indicates the need for new reference simulations...

steven-murray commented 3 years ago

I'm pretty sure this is fixed now. Feel free to re-open if you still find it buggy, @Jackmastr