Imaging-related fixes

hyperion/model/image.py

@@ -132,6 +132,17 @@ def unit(self, value):
         else:
             raise ValueError("unit should be a string")
 
+    @property
+    def units(self):
+        """
+        The units of the image values (alias for the ``.unit`` property).
+        """
+        return self.unit
+
+    @units.setter
+    def units(self, value):
+        self.unit = value
+
     @property
     def wav(self):
         """

hyperion/model/sed.py

@@ -123,6 +123,17 @@ def unit(self, value):
         else:
             raise ValueError("unit should be a string")
 
+    @property
+    def units(self):
+        """
+        The units of the SED values (alias for the ``.unit`` property).
+        """
+        return self.unit
+
+    @units.setter
+    def units(self, value):
+        self.unit = value
+
     @property
     def wav(self):
         """

hyperion/model/tests/test_image.py

@@ -728,3 +728,189 @@ def test_get_image_stokes(self, stokes):
         with pytest.raises(ValueError) as exc:
             self.m2.get_image(stokes=stokes)
         assert exc.value.args[0] == "Only the Stokes I value was stored for this image"
+
+
+def _inside_observer_direction(theta, phi):
+    t, p = np.radians(theta), np.radians(phi)
+    return np.array([np.sin(t) * np.cos(p), np.sin(t) * np.sin(p), np.cos(t)])
+
+
+def _inside_observer_expected_lonlat(d, theta_v, phi_v):
+    # Map (longitude, latitude) of a photon arriving from direction d, for an
+    # observer looking towards (theta_v, phi_v). This is d expressed in the
+    # local spherical basis (r_hat, phi_hat, -theta_hat) at the viewing
+    # direction, converted to polar coordinates - i.e. an independent
+    # reimplementation of the transform in images_peeled.f90.
+    t, p = np.radians(theta_v), np.radians(phi_v)
+    st, ct, cp, sp = np.sin(t), np.cos(t), np.cos(p), np.sin(p)
+    vx = np.dot([st * cp, st * sp, ct], d)     # r_hat . d      -> map centre / +x
+    vy = np.dot([-sp, cp, 0.], d)              # phi_hat . d    -> +longitude
+    vz = np.dot([-ct * cp, -ct * sp, st], d)   # -theta_hat . d -> +latitude
+    lon = np.degrees(np.arctan2(vy, vx))
+    lat = np.degrees(np.arctan2(np.hypot(vx, vy), vz)) - 90.
+    return lon, lat
+
+
+@pytest.mark.requires_hyperion_binaries
+def test_inside_observer_sky_coordinates(tmpdir):
+    # Regression test for the inside-observer sky-coordinate transform. A point
+    # source at cartesian position P, seen by an observer at the origin, sends
+    # photons that arrive from direction d = -P/|P|. On the all-sky map for a
+    # viewing direction (theta_v, phi_v) the source must appear at the position
+    # obtained by expressing d in the local spherical frame of the viewing
+    # direction. A previous implementation was discontinuous through
+    # theta_v = 90 deg and ignored phi_v; the viewing angles below include
+    # several cases it placed incorrectly (large phi at theta=90, either side
+    # of theta=90, and general off-axis directions).
+
+    # Place the source so that its photon arrival direction is +x.
+    photon_dir = _inside_observer_direction(90., 0.)
+    P = tuple(-0.5 * photon_dir)
+
+    views = [(90., 0.), (90., 90.), (90., 150.), (60., 0.), (120., 0.),
+             (89., 0.), (91., 0.), (45., 30.), (135., 200.)]
+
+    m = Model()
+    m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.])
+    s = m.add_point_source()
+    s.position = P
+    s.luminosity = 1.
+    s.temperature = 6000.
+    i = m.add_peeled_images(sed=False, image=True)
+    i.set_inside_observer((0., 0., 0.))
+    i.set_image_limits(180., -180., -90., 90.)
+    i.set_image_size(360, 180)
+    i.set_viewing_angles([v[0] for v in views], [v[1] for v in views])
+    i.set_wavelength_range(1, 1., 1000.)
+    m.set_n_initial_iterations(0)
+    m.set_n_photons(imaging=20000)
+    m.set_seed(-1)
+    m.write(tmpdir.join(random_id()).strpath)
+    out = m.run(tmpdir.join(random_id()).strpath)
+
+    val = out.get_image(group=0).val[:, :, :, 0]  # (n_view, n_y, n_x)
+    n_y, n_x = val.shape[1], val.shape[2]
+    xmin, xmax, ymin, ymax = 180., -180., -90., 90.
+
+    for iv, (theta_v, phi_v) in enumerate(views):
+        img = val[iv]
+        ys, xs = np.nonzero(img > 0)
+        assert len(xs) > 0, "source not found for view (%.0f, %.0f)" % (theta_v, phi_v)
+        w = img[ys, xs]
+        lon = ((xmin + (xs + 0.5) * (xmax - xmin) / n_x) * w).sum() / w.sum()
+        lat = ((ymin + (ys + 0.5) * (ymax - ymin) / n_y) * w).sum() / w.sum()
+
+        # (1) Convention-independent invariant: the great-circle distance on the
+        # map from the centre (0, 0) to the source must equal the true angle
+        # between the photon direction and the viewing direction (any rigid
+        # rotation preserves angles). This catches the tearing and the ignored
+        # phi of the old code.
+        map_sep = np.degrees(np.arccos(np.clip(
+            np.cos(np.radians(lat)) * np.cos(np.radians(lon)), -1., 1.)))
+        true_sep = np.degrees(np.arccos(np.clip(
+            np.dot(photon_dir, _inside_observer_direction(theta_v, phi_v)), -1., 1.)))
+        assert abs(map_sep - true_sep) < 2., \
+            "view (%.0f, %.0f): map separation %.1f deg != true %.1f deg" % (
+                theta_v, phi_v, map_sep, true_sep)
+
+        # (2) Full position matches the independent reference (catches roll and
+        # reflection errors that preserve the angular separation).
+        exp_lon, exp_lat = _inside_observer_expected_lonlat(photon_dir, theta_v, phi_v)
+        dlon = (lon - exp_lon + 180.) % 360. - 180.
+        assert abs(dlon) < 2. and abs(lat - exp_lat) < 2., \
+            "view (%.0f, %.0f): (%.1f, %.1f) != expected (%.1f, %.1f)" % (
+                theta_v, phi_v, lon, lat, exp_lon, exp_lat)
+
+
+def test_image_and_sed_units_constructor():
+    # Regression test: the ``units`` constructor argument of Image and SED must
+    # be routed through the validated ``unit`` property, so that it is stored
+    # (readable via ``.unit``) and invalid values are rejected. Previously it
+    # was assigned to an unvalidated attribute and left ``.unit`` unset.
+    from ..sed import SED
+    for cls in (Image, SED):
+        obj = cls(nu=np.array([1.e10, 1.e11]), units='Jy')
+        assert obj.unit == 'Jy'
+        assert obj.units == 'Jy'
+        with pytest.raises(ValueError):
+            cls(nu=np.array([1.e10, 1.e11]), units=42)
+
+
+@pytest.mark.requires_hyperion_binaries
+def test_inside_observer_flux_dilution(tmpdir):
+    # Regression test for the inside-observer 1/d^2 flux dilution. Two point
+    # sources of equal luminosity at distances d1 and d2 from the observer must
+    # have observed fluxes in the ratio (d2/d1)**2. A previous implementation
+    # divided by (d - d_min)**2 instead of d**2, giving the wrong ratio whenever
+    # a nonzero depth minimum was set.
+    d1, d2, d_min = 2., 5., 1.
+    m = Model()
+    m.set_cartesian_grid([-8., 8.], [-8., 8.], [-8., 8.])
+    for pos in [(d1, 0., 0.), (0., d2, 0.)]:
+        s = m.add_point_source()
+        s.position = pos
+        s.luminosity = 1.
+        s.temperature = 6000.
+    i = m.add_peeled_images(sed=False, image=True)
+    i.set_inside_observer((0., 0., 0.))
+    i.set_image_limits(180., -180., -90., 90.)
+    i.set_image_size(360, 180)
+    i.set_viewing_angles([90.], [0.])
+    i.set_wavelength_range(1, 1., 1000.)
+    i.set_depth(d_min, 20.)   # nonzero d_min exposes the bug
+    m.set_n_initial_iterations(0)
+    m.set_n_photons(imaging=200000)
+    m.set_seed(-1)
+    m.write(tmpdir.join(random_id()).strpath)
+    out = m.run(tmpdir.join(random_id()).strpath)
+
+    val = out.get_image(group=0).val[0, :, :, 0]   # single view, single wavelength
+    brightest = np.sort(val[val > 0])[::-1]
+    assert len(brightest) >= 2, "expected two point sources in the image"
+    # the closer source (d1) is the brighter one
+    ratio = brightest[0] / brightest[1]
+    np.testing.assert_allclose(ratio, (d2 / d1) ** 2, rtol=0.08)
+
+
+@pytest.mark.requires_hyperion_binaries
+def test_inside_observer_peeloff_optical_depth(tmpdir):
+    # Regression test for integrating the peeloff optical depth all the way to
+    # the (inside) observer. A point source at distance d sits in uniform,
+    # purely-absorbing dust of flat opacity chi and density rho, so its direct
+    # light is attenuated by exp(-chi*rho*d) over the full path to the observer.
+    # A previous implementation integrated the optical depth only to d - d_min,
+    # under-attenuating the source whenever a nonzero depth minimum was set
+    # (comparing the dust/no-dust flux ratio at fixed distance cancels the
+    # 1/d^2 dilution, isolating the optical-depth path).
+    from ...dust import IsotropicDust
+
+    d, chi, rho, d_min = 1.e16, 1., 1.e-16, 0.5e16
+    tau_full = chi * rho * d   # = 1.0
+
+    def peak_flux(with_dust):
+        m = Model()
+        m.set_cartesian_grid([-2.e16, 2.e16], [-2.e16, 2.e16], [-2.e16, 2.e16])
+        if with_dust:
+            dust = IsotropicDust([3.e9, 3.e16], [0., 0.], [chi, chi])
+            dust.set_lte_emissivities(n_temp=10, temp_min=0.1, temp_max=1.e4)
+            m.add_density_grid(np.ones((1, 1, 1)) * rho, dust)
+        s = m.add_point_source()
+        s.position = (d, 0., 0.)
+        s.luminosity = 1.
+        s.temperature = 6000.
+        i = m.add_peeled_images(sed=False, image=True)
+        i.set_inside_observer((0., 0., 0.))
+        i.set_image_limits(180., -180., -90., 90.)
+        i.set_image_size(360, 180)
+        i.set_viewing_angles([90.], [0.])
+        i.set_wavelength_range(1, 1., 1000.)
+        i.set_depth(d_min, 3.e16)   # nonzero d_min previously truncated the tau path
+        m.set_n_initial_iterations(0)
+        m.set_n_photons(imaging=100000)
+        m.set_seed(-1)
+        m.write(tmpdir.join(random_id()).strpath)
+        out = m.run(tmpdir.join(random_id()).strpath)
+        return out.get_image(group=0).val[0, :, :, 0].max()
+
+    ratio = peak_flux(with_dust=True) / peak_flux(with_dust=False)
+    np.testing.assert_allclose(ratio, np.exp(-tau_full), rtol=0.05)

hyperion/model/tests/test_sed.py

@@ -470,3 +470,34 @@ def test_sed_negative_group(tmpdir):
     m.write(tmpdir.join(random_id()).strpath)
     mo = m.run()
     np.testing.assert_allclose(mo.get_sed(group=-1).val, mo.get_sed(group=0).val)
+
+
+@pytest.mark.requires_hyperion_binaries
+def test_sed_uncertainty_sum_of_squares(tmpdir):
+    # Regression test for the Monte-Carlo uncertainty estimator. For a single
+    # isotropic point source with no dust, every imaging photon contributes the
+    # same weight w to the single SED bin, so the sum-of-squares estimator gives
+    # sigma = sqrt(sum w**2) = sqrt(N) w against a total flux F = N w, i.e. a
+    # relative uncertainty of exactly 1/sqrt(N). A previous (sign-flipped)
+    # sample-variance formula over-estimated this by a factor sqrt(2).
+    n_photons = 10000
+    m = Model()
+    m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.])
+    s = m.add_point_source()
+    s.luminosity = 1.
+    s.temperature = 6000.
+    i = m.add_peeled_images(sed=True, image=False)
+    i.set_viewing_angles([45.], [45.])
+    i.set_aperture_radii(1, 1.e10, 1.e10)
+    i.set_wavelength_range(1, 0.01, 5000.)
+    i.set_uncertainties(True)
+    m.set_n_initial_iterations(0)
+    m.set_n_photons(imaging=n_photons)
+    m.set_seed(-1)
+    m.write(tmpdir.join(random_id()).strpath)
+    out = m.run(tmpdir.join(random_id()).strpath)
+
+    sed = out.get_sed(uncertainties=True)
+    flux = np.nansum(sed.val)
+    sigma = np.sqrt(np.nansum(np.array(sed.unc) ** 2))
+    np.testing.assert_allclose(sigma / flux, 1. / np.sqrt(n_photons), rtol=0.03)

src/images/image_type.f90

@@ -665,13 +665,12 @@ subroutine image_write(img,group)
 
        cube5d = img%sed
 
-       if(img%uncertainties) then
-          where(img%sedn > 1)
-             cube5de = sqrt((img%sed2 + (img%sed)**2 / img%sedn) / (img%sedn - 1)) * sqrt(img%sedn)
-          elsewhere
-             cube5de = 0._dp
-          end where
-       end if
+       ! The uncertainty on the total flux in each bin is estimated with the
+       ! standard Monte-Carlo sum-of-squares estimator sqrt(sum(x_i**2)),
+       ! which accounts for both the scatter in the photon weights and the
+       ! Poisson counting uncertainty.
+
+       if(img%uncertainties) cube5de = sqrt(img%sed2)
 
        if(.not.img%use_exact_nu) then
           cube5d = cube5d / dnunorm
@@ -727,13 +726,9 @@ subroutine image_write(img,group)
 
        cube6d = img%img
 
-       if(img%uncertainties) then
-          where(img%imgn > 1)
-             cube6de = sqrt((img%img2 + (img%img)**2 / img%imgn) / (img%imgn - 1)) * sqrt(img%imgn)
-          elsewhere
-             cube6de = 0._dp
-          end where
-       end if
+       ! See the SED branch above for the derivation of the uncertainties
+
+       if(img%uncertainties) cube6de = sqrt(img%img2)
 
        if(.not.img%use_exact_nu) then
           cube6d = cube6d / dnunorm

src/images/images_peeled.f90

@@ -100,8 +100,8 @@ subroutine peeloff_photon(p_orig, polychromatic)
     type(photon) :: p
     real(dp) :: tau
     integer :: ip,ig,iv,id
-    type(angle3d_dp) :: a_req, a_sky, a_diff
-    type(vector3d_dp) :: v_req
+    type(angle3d_dp) :: a_req, a_view
+    type(vector3d_dp) :: v_req, v_a, v_sky
     logical,intent(in) :: polychromatic
     real(dp) :: x_image, y_image
     real(dp) :: tmax
@@ -151,29 +151,53 @@ subroutine peeloff_photon(p_orig, polychromatic)
        ! the xmax wall.
        call place_in_cell(p)
 
+       ! The depth range given by d_min and d_max is only used to select
+       ! which events to include in the image - the optical depth should be
+       ! integrated all the way to the observer (for inside observers) or to
+       ! the edge of the grid (for external observers).
        if(inside_observer(ig)) then
           dr = p%r-r_peeloff(ig)
           d = sqrt(dr.dot.dr)
-          tmax = d - d_min(ig)
+          tmax = d
        else
           d = -(v_req.dot.p%r)
-          tmax = - d_min(ig) + d
+          tmax = huge(1._dp)
        end if
 
        if(d < d_min(ig) .or. d > d_max(ig)) cycle
 
        if(inside_observer(ig)) then
 
-          if(abs(viewing_angles(ip)%cost) .gt. 1.e-10) then
-             call rotate_angle3d(viewing_angles(ip), angle3d_deg(90._dp, 0._dp), a_diff)
-             call difference_angle3d(a_diff, -p%a, a_sky)
-          else
-             a_sky = p%a
-          end if
-
-          ! Convert to angles in degrees
-          x_image = atan2(a_sky%sinp, a_sky%cosp) * rad2deg
-          y_image = atan2(a_sky%sint, a_sky%cost) * rad2deg - 90._dp
+          ! Sky position of the photon for an observer looking towards the
+          ! direction a_view = (theta_v, phi_v). The map is centred on the
+          ! viewing direction, so we express the photon direction p%a in the
+          ! local orthonormal frame at a_view given by the spherical basis
+          ! vectors (r_hat, phi_hat, -theta_hat), writing t, p for theta_v,
+          ! phi_v:
+          !
+          !   r_hat      = ( sin t cos p,  sin t sin p,  cos t )  (viewing dir)
+          !   phi_hat    = (      -sin p,        cos p,      0  )  (local east)
+          !   -theta_hat = (-cos t cos p, -cos t sin p,  sin t )  (local north)
+          !
+          ! The three components of v_sky below are the projections of p%a onto
+          ! these axes, i.e. v_sky = R p%a where R is the rigid rotation whose
+          ! rows are the vectors above and which therefore sends a_view to the
+          ! map centre +x = (theta, phi) = (90, 0). R is orthonormal, so this is
+          ! exactly continuous in (theta_v, phi_v) with no special cases at the
+          ! poles (theta_v = 0, 90, 180) - unlike a frame built from the
+          ! pole-referenced rotate_angle3d / difference_angle3d routines.
+          a_view = viewing_angles(ip)
+          call angle3d_to_vector3d(p%a, v_a)
+          v_sky%x = (v_a%x * a_view%cosp + v_a%y * a_view%sinp) * a_view%sint + v_a%z * a_view%cost
+          v_sky%y = - v_a%x * a_view%sinp + v_a%y * a_view%cosp
+          v_sky%z = - (v_a%x * a_view%cosp + v_a%y * a_view%sinp) * a_view%cost + v_a%z * a_view%sint
+
+          ! Convert to map coordinates: longitude is the azimuth of v_sky in the
+          ! x-y plane, and latitude is (colatitude measured from +z) - 90 deg.
+          ! Both use atan2 on the Cartesian components, so they stay well-defined
+          ! when v_sky is near the poles (where sin(colatitude) -> 0).
+          x_image = atan2(v_sky%y, v_sky%x) * rad2deg
+          y_image = atan2(sqrt(v_sky%x**2 + v_sky%y**2), v_sky%z) * rad2deg - 90._dp
 
           ! Make sure the photon falls inside the image (wrap angles around)
           x_image = peeled_image(ig)%x_max + modulo(x_image - peeled_image(ig)%x_max, 360._dp)
@@ -205,12 +229,11 @@ subroutine peeloff_photon(p_orig, polychromatic)
              end if
           end if
 
-          ! For inside observer, don't want optical depth to escape grid, just to go to observer!
-          ! Need to include 1/d^2!
-
           if(.not.killed) then
 
-             if(inside_observer(ig)) p%s = p%s / (4._dp * pi * tmax**2._dp)
+             ! For inside observers, include the 1/d^2 flux dilution
+
+             if(inside_observer(ig)) p%s = p%s / (4._dp * pi * d**2._dp)
 
              if(polychromatic) then