spacetelescope / webbpsf

James Webb Space Telescope PSF simulation tool
https://webbpsf.readthedocs.io
BSD 3-Clause "New" or "Revised" License
119 stars 63 forks source link

poppy.test() fail on python 2.7 #166

Closed mperrin closed 6 years ago

mperrin commented 6 years ago

Issue by benjaminpope Tuesday Jun 27, 2017 at 06:59 GMT Originally opened as https://github.com/mperrin/webbpsf/issues/166


Hi Marshall,

I am hoping to use poppy for an aperture mask simulation I'm doing. I just ran poppy.test() on a brand new conda 4.3.21 python 2.7 install on Mac Sierra and it seems to fail on 3 tests unexpectedly. I have copied the full output below. Is this a known issue?

Using astropy 1.3.3, scipy 0.19.1, numpy 1.13.0, matplotlib 2.0.2.

All the best,

Ben

=================================== FAILURES =================================== _____ test_optic_resizing __

@pytest.mark.skipif(
    (scipy is None),
    reason='No SciPy installed'
)
def test_optic_resizing():
    '''
    Tests the rescaling functionality of OpticalElement.getPhasor(),
    by first creating an optic with a small pixel scale and then
    creating an optic with a large pixel scale, and checking the returned
    phasor of each has the dimensions of the input wavefront.
    '''

    # diameter 1 meter, pixel scale 2 mm
    inputwf = poppy_core.Wavefront(diam=1.0, npix=500)

    # Test rescaling from finer scales: diameter 1 meter, pixel scale 1 mm
    test_optic_small=fits.HDUList([fits.PrimaryHDU(np.zeros([1000,1000]))])
    test_optic_small[0].header["PUPLSCAL"]=.001
    test_optic_small_element=poppy_core.FITSOpticalElement(transmission=test_optic_small)
  assert(test_optic_small_element.getPhasor(inputwf).shape ==inputwf.shape )

/anaconda/lib/python2.7/site-packages/poppy/tests/test_core.py:348:


/anaconda/lib/python2.7/site-packages/poppy/poppy_core.py:2315: in getPhasor return self.get_phasor(wave)


self = <poppy.poppy_core.FITSOpticalElement object at 0x111470e50> wave = <poppy.poppy_core.Wavefront object at 0x111470cd0>

def get_phasor(self,wave):
    """ Compute a complex phasor from an OPD, given a wavelength.

        The returned value should be the complex phasor array as appropriate for
        multiplying by the wavefront amplitude.

        Parameters
        ----------
        wave : float or obj
            either a scalar wavelength or a Wavefront object

        """
    #_log.info("Pixelscales for %s: wave %f, optic  %f" % (self.name, wave.pixelscale, self.pixelscale))

    if isinstance(wave, Wavefront):
        wavelength=wave.wavelength
    else:
        wavelength=wave
    scale = 2. * np.pi / wavelength.to(u.meter).value

    # set the self.phasor attribute:
    # first check whether we need to interpolate to do this.
    float_tolerance = 0.001  #how big of a relative scale mismatch before resampling?
    if self.pixelscale is not None and hasattr(wave,'pixelscale') and abs(wave.pixelscale -self.pixelscale)/self.pixelscale >= float_tolerance:
        _log.debug("Pixelscales: wave {}, optic {}" .format(wave.pixelscale, self.pixelscale))
        #raise ValueError("Non-matching pixel scale for wavefront and optic! Need to add interpolation / ing ")
        if hasattr(self,'_resampled_scale') and abs(self._resampled_scale-wave.pixelscale)/self._resampled_scale >= float_tolerance:
            # we already did this same resampling, so just re-use it!
            self.phasor = self._resampled_amplitude * np.exp (1.j * self._resampled_opd * scale)
        else:
            #raise NotImplementedError("Need to implement resampling.")
            zoom=(self.pixelscale/wave.pixelscale).decompose().value
            resampled_opd = scipy.ndimage.interpolation.zoom(self.opd, zoom,
                    output=self.opd.dtype, order=self.interp_order)
            resampled_amplitude = scipy.ndimage.interpolation.zoom(self.amplitude,zoom,output=self.amplitude.dtype,order=self.interp_order)
            _log.debug("resampled optic to match wavefront via spline interpolation by a zoom factor of %.3g"%(zoom))
            _log.debug("resampled optic shape: {}   wavefront shape: {}".format(resampled_amplitude.shape, wave.shape))

            lx,ly=resampled_amplitude.shape
            #crop down to match size of wavefront:
            lx_w,ly_w = wave.amplitude.shape
            border_x = np.abs(np.floor((lx-lx_w)/2))
            border_y = np.abs(np.floor((ly-ly_w)/2))
            if (self.pixelscale*self.amplitude.shape[0] < wave.pixelscale*wave.amplitude.shape[0]) or (self.pixelscale*self.amplitude.shape[1] < wave.pixelscale*wave.amplitude.shape[0]):
                #raise ValueError("Optic is smaller than input wavefront")
                _log.warn("Optic"+str(np.shape(resampled_opd))+" is smaller than input wavefront"+str([lx_w,ly_w])+", will attempt to zero-pad the rescaled array")
                self._resampled_opd = np.zeros([lx_w,ly_w])
                self._resampled_amplitude = np.zeros([lx_w,ly_w])

                self._resampled_opd[border_x:border_x+resampled_opd.shape[0],border_y:border_y+resampled_opd.shape[1]] = resampled_opd
                self._resampled_amplitude[border_x:border_x+resampled_opd.shape[0],border_y:border_y+resampled_opd.shape[1]]=resampled_amplitude
                _log.debug("padded an optic with a %i x %i border to optic to match the wavefront"%(border_x,border_y))

            else:
              self._resampled_opd = resampled_opd[border_x:border_x+lx_w,border_y:border_y+ly_w]

E TypeError: slice indices must be integers or None or have an index method

/anaconda/lib/python2.7/site-packages/poppy/poppy_core.py:2289: TypeError __ test_fresnel_optical_system_Hubble __

display = False

def test_fresnel_optical_system_Hubble(display=False):
    """ Test the FresnelOpticalSystem infrastructure
    This is a fairly comprehensive test case using as its
    example a simplified version of the Hubble Space Telescope.

    The specific numeric values used for Hubble are taken from
    the HST example case in the PROPER manual by John Krist,
    version 2.0a available from http://proper-library.sourceforge.net
    This is not intended as a high-fidelity model of Hubble, and this
    test case neglects the primary aperture obscuration as well as the
    specific conic constants of the optics including the infamous
    spherical aberration.

    This function tests the FresnelOpticalSystem functionality including
    assembly of the optical system and propagation of wavefronts,
    intermediate beam sizes through the optical system,
    intermediate and final system focal lengths,
    toggling between different types of optical planes,
    and the properties of the output PSF including FWHM and
    comparison to the Airy function.

    """

    # HST example - Following example in PROPER Manual V2.0 page 49.
    # This is an idealized case and does not correspond precisely to the real telescope
    # Define system with units
    diam = 2.4 * u.m
    fl_pri = 5.52085 * u.m
    d_pri_sec = 4.907028205 * u.m  # This is what's used in the PROPER example
    #d_pri_sec = 4.9069 * u.m      # however Lallo 2012 gives this value, which differs slightly
                                   # from what is used in the PROPER example case.
    fl_sec = -0.6790325 * u.m
    d_sec_to_focus = 6.3916645 * u.m # place focal plane right at the beam waist after the SM

    osamp = 2 #oversampling factor

    hst = fresnel.FresnelOpticalSystem(pupil_diameter=2.4*u.m, beam_ratio=0.25)
    g1 = fresnel.QuadraticLens(fl_pri, name='Primary', planetype=poppy_core.PlaneType.pupil)
    g2 = fresnel.QuadraticLens(fl_sec, name='Secondary')

    hst.add_optic(optics.CircularAperture(radius=diam.value/2))
    hst.add_optic(g1)
    hst.add_optic(g2, distance=d_pri_sec)
    hst.add_optic(optics.ScalarTransmission(planetype=poppy_core.PlaneType.image), distance=d_sec_to_focus)

    # Create a PSF
    psf, waves = hst.calc_psf(wavelength=0.5e-6, display_intermediates=display, return_intermediates=True)

    ### check the beam size is as expected at primary and secondary mirror
    assert(np.allclose(waves[1].spot_radius().value, 1.2))
    # can't find a definitive reference for the beam diam at the SM, but
    # the secondary mirror's radius is 14.05 cm
    # We find that the beam is indeed slightly smaller than that.
    assert(waves[2].spot_radius() > 13*u.cm )
    assert(waves[2].spot_radius() < 14*u.cm )

    ### check the focal length of the overall system is as expected
    expected_system_focal_length = 1./(1./fl_pri + 1./fl_sec - (d_pri_sec)/(fl_pri*fl_sec))
    # n.b. the value calculated here, 57.48 m, is a bit less than the
    # generally stated focal length of Hubble, 57.6 meters. Adjusting the
    # primary-to-secondary spacing by about 100 microns can resolve this
    # discrepancy. We here opt to stick with the values used in the PROPER
    # example, to facilitate cross-checking the two codes.

    assert(not np.isfinite(waves[0].focal_length))  # plane wave after circular aperture
    assert(waves[1].focal_length==fl_pri)           # focal len after primary
    # NB. using astropy.Quantities with np.allclose() doesn't work that well
    # so pull out the values here:
    assert(np.allclose(waves[2].focal_length.to(u.m).value,
           expected_system_focal_length.to(u.m).value)) # focal len after secondary

    ### check the FWHM of the PSF is as expected
    measured_fwhm = utils.measure_fwhm(psf)
    expected_fwhm = 1.028*0.5e-6/2.4*206265
    # we only require this to have < 5% accuracy with respect to the theoretical value
    # given discrete pixelization etc.
    assert(np.abs((measured_fwhm-expected_fwhm)/expected_fwhm) < 0.05)

    ### check the various plane types are as expected, including toggling into angular coordinates
    assert_message = ("Expected FresnelWavefront at plane #{} to have {} == {}, but got {}")
    system_planetypes = [PlaneType.pupil, PlaneType.pupil, PlaneType.intermediate, PlaneType.image]
    for idx, (wavefront, planetype) in enumerate(zip(waves, system_planetypes)):
        assert wavefront.planetype == planetype, assert_message.format(
            idx, "planetype", plane_type, wavefront.planetype
        )

    angular_coordinates_flags = [False, False, False, True]
    for idx, (wavefront, angular_coordinates) in enumerate(zip(waves, angular_coordinates_flags)):
        assert wavefront.angular_coordinates == angular_coordinates, assert_message.format(
            idx, "angular_coordinates", angular_coordinates, wavefront.angular_coordinates
        )

    spherical_flags = [False, True, True, False]
    for idx, (wavefront, spherical) in enumerate(zip(waves, spherical_flags)):
        assert wavefront.spherical == spherical, assert_message.format(
            idx, "spherical", spherical, wavefront.spherical
        )

    ### and check that the resulting function is a 2D Airy function
    #create an airy function matching the center part of this array
    airy = misc.airy_2d(diameter=diam.value, wavelength=0.5e-6,
                              shape=(128,128), pixelscale=psf[0].header['PIXELSCL'],
                             center=(64,64))

    centerpix = hst.npix / hst.beam_ratio / 2
  cutout = psf[0].data[centerpix-64:centerpix+64, centerpix-64:centerpix+64] / psf[0].data[centerpix,centerpix]

E TypeError: slice indices must be integers or None or have an index method

/anaconda/lib/python2.7/site-packages/poppy/tests/test_fresnel.py:354: TypeError _ test_estimatenprocesses

def test_estimate_nprocesses():
    """ Apply some basic functionality tests to the
    estimate nprocesses function.
    """
    osys = poppy_core.OpticalSystem("test")
    osys.addPupil(optics.CircularAperture(radius=1))
    osys.addPupil(optics.CircularAperture(radius=0.5))
    osys.addDetector(pixelscale=0.1, fov_arcsec=2.0)

    answer = utils.estimate_optimal_nprocesses(osys)

    #see if it's an int with a reasonable value
    assert type(answer) is int
  assert answer > 0, "Estimated optimal nprocesses must be positive integer"

E AssertionError: Estimated optimal nprocesses must be positive integer E assert 0 > 0

/anaconda/lib/python2.7/site-packages/poppy/tests/test_multiprocessing.py:105: AssertionError ========= 3 failed, 92 passed, 6 skipped, 1 xfailed in 170.51 seconds ==========

mperrin commented 6 years ago

Comment by mperrin Tuesday Jun 27, 2017 at 15:12 GMT


Hmm, thanks for bringing this to our attention. The first two errors are clearly due to the recent change in numpy handling of indices, but I thought we'd fixed all those cases already. The third test failure's a new one I don't recognize but I suspect it may be some related float/integer issue.

Thanks for bringing this to our attention, and we will try to look into it more sometime soon.

mperrin commented 6 years ago

Comment by josePhoenix Tuesday Jun 27, 2017 at 18:49 GMT


I thought we'd taken care of these, too. @benjaminpope could you say which version of POPPY you're using?

Also, if it's the current version, would you mind opening an issue on mperrin/poppy instead? Thanks :)

mperrin commented 6 years ago

Comment by josePhoenix Tuesday Jun 27, 2017 at 19:07 GMT


Ah, they were fixed on master since the 0.5.1 release. I would recommend installing POPPY from the development sources (see http://pythonhosted.org/poppy/installation.html) until we get a chance to release 0.5.2 or 0.6.0.

mperrin commented 6 years ago

Comment by benjaminpope Wednesday Jun 28, 2017 at 03:44 GMT


Ah, thanks! Dev source worked. Was using standard pip distribution. Can close this issue then.


From: Joseph Long [notifications@github.com] Sent: Tuesday, June 27, 2017 8:07 PM To: mperrin/webbpsf Cc: Benjamin Pope; Mention Subject: Re: [mperrin/webbpsf] poppy.test() fail on python 2.7 (#166)

Ah, they were fixed on master since the 0.5.1 release. I would recommend installing POPPY from the development sources (see http://pythonhosted.org/poppy/installation.html) until we get a chance to release 0.5.2 or 0.6.0.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/mperrin/webbpsf/issues/166#issuecomment-311455259, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AFvsFU0G-RQ2ja6707GE4uXBmJPV0D7Hks5sIVL7gaJpZM4OGN8k.