diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index d52b4e8..9e54c21 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -25,6 +25,7 @@ jobs: scipy-version: ['<1.6', '<1.7'] speclite-version: ['==0.19'] desiutil-version: ['3.4.2'] + desimodel-version: ['0.19.1'] steps: - name: Checkout code @@ -46,7 +47,7 @@ jobs: python -m pip install 'numpy${{ matrix.numpy-version }}' python -m pip install 'speclite${{ matrix.speclite-version }}' python -m pip install git+https://github.com/desihub/desiutil.git@${{ matrix.desiutil-version }} - python -m pip install git+https://github.com/desihub/desimodel.git@0.09.1 + python -m pip install git+https://github.com/desihub/desimodel.git@${{ matrix.desimodel-version }} - name: Run the test run: pytest @@ -77,7 +78,7 @@ jobs: python -m pip install pytest pytest-astropy coveralls python -m pip install pyyaml numpy\<1.21 scipy\<1.6 astropy\<5.0 speclite==0.19 python -m pip install git+https://github.com/desihub/desiutil.git@3.4.2 - python -m pip install git+https://github.com/desihub/desimodel.git@0.09.1 + python -m pip install git+https://github.com/desihub/desimodel.git@0.19.1 - name: Run the test with coverage run: pytest --cov @@ -141,7 +142,7 @@ jobs: run: | python -m pip install --upgrade pip wheel setuptools flake8 - - name: Test the style; failures are allowed + - name: Test the style; failures are not allowed # This is equivalent to an allowed falure. - continue-on-error: true + # continue-on-error: true run: flake8 specsim --count --max-line-length=100 diff --git a/specsim/conftest.py b/specsim/conftest.py index dbdfd0c..128660e 100644 --- a/specsim/conftest.py +++ b/specsim/conftest.py @@ -7,9 +7,9 @@ # With older versions of Astropy, we actually need to import the pytest # plugins themselves in order to make them discoverable by pytest. from astropy.tests.pytest_plugins import * - + # As of Astropy 5.1 astropy.tests.plugins.display (where PYTEST_HEADER_MODULES -# and TESTED_VERSIONS lived) as been deprecated and removed entirely. +# and TESTED_VERSIONS lived) as been deprecated and removed entirely. ## Uncomment the following lines to treat all DeprecationWarnings as ## exceptions diff --git a/specsim/fiberloss.py b/specsim/fiberloss.py index ab1a318..19b22cd 100644 --- a/specsim/fiberloss.py +++ b/specsim/fiberloss.py @@ -356,7 +356,7 @@ def calculate_fiber_acceptance_fraction( if len(focal_y) != num_fibers: raise ValueError('Arrays focal_x and focal_y must have same length.') - + # Use pre-tabulated fiberloss fractions when requested. if instrument.fiberloss_method == 'table': @@ -371,7 +371,7 @@ def calculate_fiber_acceptance_fraction( floss[i] = instrument.fiber_acceptance_dict[type_name] return floss - # Otherwise, use GalSim or the fastfiberacceptance calibrated on galsim + # Otherwise, use GalSim or the fastfiberacceptance calibrated on galsim # to calculate fiberloss fractions on the fly... # Initialize the grid of wavelengths where the fiberloss will be @@ -380,14 +380,14 @@ def calculate_fiber_acceptance_fraction( wlen_unit = wavelength.unit wlen_grid = np.linspace(wavelength.value[0], wavelength.value[-1], num_wlen) * wlen_unit - + # Calculate the focal-plane optics at the fiber locations. scale, blur, offset = instrument.get_focal_plane_optics( focal_x, focal_y, wlen_grid) # Calculate the atmospheric seeing at each wavelength. seeing_fwhm = atmosphere.get_seeing_fwhm(wlen_grid).to(u.arcsec).value - + # Replicate source parameters from the source config if they are not # provided via args. If they are provided, check for the expected shape. if source_fraction is None: @@ -416,32 +416,32 @@ def calculate_fiber_acceptance_fraction( [num_fibers, 1]) elif source_position_angle.shape != (num_fibers, 2): raise ValueError('Unexpected shape for source_position_angle.') - + fiberloss_grid = None - + # choose here how to compute things if instrument.fiberloss_method == 'fastsim': scale_um_per_arcsec = scale.to(u.um / u.arcsec).value blur_um = blur.to(u.um).value - + sigma = np.zeros((offset.shape[0],offset.shape[1])) disk_half_light_radius = np.zeros(sigma.shape) bulge_half_light_radius = np.zeros(sigma.shape) disk_frac = np.zeros(sigma.shape) bulge_frac = np.zeros(sigma.shape) - + for i in range(offset.shape[0]) : sigma[i] = np.sqrt( (seeing_fwhm/2.35482)**2*scale_um_per_arcsec[i,0]*scale_um_per_arcsec[i,1]+blur_um[i]**2 ) disk_half_light_radius[i] = source_half_light_radius[i,0]*np.ones(offset.shape[1]) bulge_half_light_radius[i] = source_half_light_radius[i,1]*np.ones(offset.shape[1]) disk_frac[i] = source_fraction[i,0]*np.ones(offset.shape[1]) bulge_frac[i] = source_fraction[i,1]*np.ones(offset.shape[1]) - + point_frac = 1 - disk_frac - bulge_frac - + offset_um = offset.to(u.um).value delta = np.sqrt(offset_um[:,:,0]**2+offset_um[:,:,1]**2) - + fiberloss_grid = np.zeros(sigma.shape) if np.sum(point_frac)>0 : fiberloss_grid += point_frac * instrument.fast_fiber_acceptance.value("POINT",sigma,delta) @@ -449,9 +449,9 @@ def calculate_fiber_acceptance_fraction( fiberloss_grid += disk_frac * instrument.fast_fiber_acceptance.value("DISK",sigma,delta,disk_half_light_radius) if np.sum(bulge_frac)>0 : fiberloss_grid += bulge_frac * instrument.fast_fiber_acceptance.value("BULGE",sigma,delta,bulge_half_light_radius) - + else : - + # Initialize a new calculator. calc = GalsimFiberlossCalculator( instrument.fiber_diameter.to(u.um).value, diff --git a/specsim/fitgalsim.py b/specsim/fitgalsim.py index ee77664..e6e9f87 100644 --- a/specsim/fitgalsim.py +++ b/specsim/fitgalsim.py @@ -23,7 +23,7 @@ def generate_fiber_positions(nfiber, seed, desi): """ returns random fiber location on focal surface - + Args: nfibers (int) number of fiber seed (int) random seed @@ -50,17 +50,17 @@ def main(args=None): The method consists in first setting sigma and offset grid, and reverse engineering the atmospheric seeing to retrieve the correct effective sigma on the grid given the telescope blur. - + For each point in the output parameter grid (source type , sigma, offset, - source radius), several calculations are done with random angular + source radius), several calculations are done with random angular orientation of fiber and source to account for the fiber ellipticity (due to anisotropic plate scale) and source ellipticity. - The output file has to be saved in + The output file has to be saved in $DESIMODEL/data/throughput/galsim-fiber-acceptance.fits to be used for fast fiber acceptance computation. This idea is to compute accurate and correlated values of offset, blur, - scale, atmospheric seeing from the fiber location and wavelength, + scale, atmospheric seeing from the fiber location and wavelength, compute the effective sigma and read with a ND interpolation the fiber acceptance value from the file. """ @@ -78,12 +78,12 @@ def main(args=None): axis_ratio = 0.7 # a fixed axis ratio is used for DISK and BULGE (position angles are random) sources = ["POINT","DISK","BULGE"] ######################################################################## - + print("init simulator") - + desi = specsim.simulator.Simulator('desi') wave = np.linspace(6000.,6001.,nsigma) # Angstrom , wavelength are not used - + # optics with random fiber positions to get the range of scale and blur x,y = generate_fiber_positions(nrand, seed, desi) x=x.to(u.um).value @@ -97,24 +97,24 @@ def main(args=None): offset = np.linspace(0,max_offset,noffset)*mscale # this is the final offset array I will save, um sigma = np.linspace(min_sigma,max_sigma,nsigma)*mscale # this is the final sigma array I will save, um - + # random orientations of sources (account for source ellipticity) position_angle_source_deg = 360.*np.random.uniform(size=nrand) - + # random orientations of offsets (account for plate scale asymetry) theta = 2*np.pi*np.random.uniform(size=nrand) rcos_offset = np.cos(theta) rsin_offset = np.sin(theta) - + # init galsim calculator fiber_diameter = desi.instrument.fiber_diameter.to(u.um).value calc = specsim.fiberloss.GalsimFiberlossCalculator(fiber_diameter=fiber_diameter,wlen_grid=wave, num_pixels=16,oversampling=32,moffat_beta=3.5) hdulist=None - + for source in sources : - + nfibers=noffset disk_bulge_fraction = np.zeros((nfibers,2)) # fraction of disk and bulge minor_major_axis_ratio = axis_ratio*np.ones((nfibers,2)) # minor / major axis ratio , for disk and bulge component, respectively @@ -139,21 +139,21 @@ def main(args=None): print("computing fiberloss for",source,"hlr=",half_light_radius_value,"arcsec") sys.stdout.flush() - + sloss=np.zeros((noffset,nsigma)) # sum of loss values sloss2=np.zeros((noffset,nsigma)) # sum of loss2 values for r in range(nrand) : blur2=np.mean(blur[r,:]**2) # scalar, mean blur scale2=scale[r,0]*scale[r,1] # scalar ,sigmax*sigmay - + # we artificially set the seeing array to span the seeing range instead of following # evolution with wavelength # # this is the inverse of (in fiberloss.py) : # sigma[i] = np.sqrt( (seeing_fwhm/2.35482)**2*scale_um_per_arcsec[i,0]*scale_um_per_arcsec[i,1]+blur_um[i]**2 ) - + atmospheric_seeing = np.sqrt((sigma**2 - blur2)/scale2)/fwhm_to_sigma # size nsigma , arcsec, fwhm - + galsim_scale = np.zeros((noffset,2)) galsim_offset = np.zeros((noffset,nsigma,2)) galsim_blur = np.zeros((noffset,nsigma)) @@ -163,18 +163,18 @@ def main(args=None): galsim_blur[i,:] = blur[r,:] # same blur as a function of wavelength galsim_offset[i,:,0] = offset[i]*rcos_offset[r] # apply a random angle to the offset galsim_offset[i,:,1] = offset[i]*rsin_offset[r] # apply a random angle to the offset - position_angle[:,:] = position_angle_source_deg[r] # degrees - + position_angle[:,:] = position_angle_source_deg[r] # degrees + # calculate using galsim (that's long) loss = calc.calculate(seeing_fwhm=atmospheric_seeing,scale=galsim_scale, offset=galsim_offset,blur_rms=galsim_blur, source_fraction=disk_bulge_fraction,source_half_light_radius=half_light_radius, source_minor_major_axis_ratio=minor_major_axis_ratio, source_position_angle=position_angle) - + sloss += loss sloss2 += loss**2 - mloss=sloss/nrand + mloss=sloss/nrand mean_loss.append( mloss.T ) rms_loss.append( np.sqrt(sloss2/nrand-mloss**2).T ) @@ -201,7 +201,7 @@ def main(args=None): moffat_beta=3.5 header.add_comment("galsim Atmospheric PSF: Moffat with beta=%2.1f"%moffat_beta) header.add_comment("galsim Telescope blur PSF: Gaussian") - + if len(mean_loss)>1 : mean_loss=np.array(mean_loss) diff --git a/specsim/observation.py b/specsim/observation.py index ea41851..1fbe734 100644 --- a/specsim/observation.py +++ b/specsim/observation.py @@ -69,7 +69,7 @@ def _update_model(self): # Initialize an observing model at the middle of the exposure and # at the central wavelength of the simulation, i.e., ignore temporal # and chromatic variations (for now). - + # This calculation can raise a non-catastrophic "overflow encountered in # double_scalars" numpy error; catch it here. with np.errstate(all='ignore'): diff --git a/specsim/simulator.py b/specsim/simulator.py index 57f69c2..6ebb713 100644 --- a/specsim/simulator.py +++ b/specsim/simulator.py @@ -587,7 +587,7 @@ def generate_random_noise(self, random_state=None, use_poisson=True): if None is specified. use_poisson : bool If False, use numpy.random.normal instead of numpy.random.poisson. - This is useful for simulations where one wants the same noise + This is useful for simulations where one wants the same noise realization for varying average flux (numpy.random.poisson uses a varying number of random numbers depending on the mean). """