From 4122a1fac1076e87194018d8409564d997b568e9 Mon Sep 17 00:00:00 2001 From: Julien Date: Mon, 25 Jan 2021 15:15:13 +0100 Subject: [PATCH 01/15] lint tod.py --- s4cmb/tod.py | 1408 ++++++++++++++++++++++++++++++-------------------- 1 file changed, 856 insertions(+), 552 deletions(-) diff --git a/s4cmb/tod.py b/s4cmb/tod.py index 5a53eeb..74ccbb2 100755 --- a/s4cmb/tod.py +++ b/s4cmb/tod.py @@ -13,6 +13,7 @@ import numpy as np import healpy as hp + # Python3 does not have the cPickle module try: import cPickle as pickle @@ -31,19 +32,47 @@ from s4cmb.xpure import qu_weight_mineig d2r = np.pi / 180.0 -am2rad = np.pi / 180. / 60. +am2rad = np.pi / 180.0 / 60.0 + -class TimeOrderedDataPairDiff(): +class TimeOrderedDataPairDiff: """ Class to handle Time-Ordered Data (TOD) """ - def __init__(self, hardware, scanning_strategy, HealpixFitsMap, - CESnumber, projection='healpix', - nside_out=None, pixel_size=None, width=140., - cut_pixels_outside=True, - array_noise_level=None, array_noise_seed=487587, - array_noise_level2=None, array_noise_seed2=56736, - nclouds=None, corrlength=None, alpha=None, - f0=None, amp_atm=None, mapping_perpair=False, mode='standard', - store_pointing_matrix_input = False, verbose=False, perturb_el = False, perturb_az = False, seed_pointing = 0., mu_pointing = 0., sigma_pointing = 13.,perturb_pol_angs=False,seed_pa=0.,pa_mu=0.,pa_sig=0.,pa_diff=False): + + def __init__( + self, + hardware, + scanning_strategy, + HealpixFitsMap, + CESnumber, + projection="healpix", + nside_out=None, + pixel_size=None, + width=140.0, + cut_pixels_outside=True, + array_noise_level=None, + array_noise_seed=487587, + array_noise_level2=None, + array_noise_seed2=56736, + nclouds=None, + corrlength=None, + alpha=None, + f0=None, + amp_atm=None, + mapping_perpair=False, + mode="standard", + store_pointing_matrix_input=False, + verbose=False, + perturb_el=False, + perturb_az=False, + seed_pointing=0.0, + mu_pointing=0.0, + sigma_pointing=13.0, + perturb_pol_angs=False, + seed_pa=0.0, + pa_mu=0.0, + pa_sig=0.0, + pa_diff=False, + ): """ C'est parti! @@ -121,25 +150,26 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, Store pointing matrix used to scan the input map perturb_el : boolean, optional If True, perturb one of a few pointing instances. - If False, sets error arrays to 0 so that the additional pointing instances - are the same as the original ones. + If False, sets error arrays to 0 so that the additional + pointing instances are the same as the original ones. perturb_az : boolean, optional If True, perturb one of a few pointing instances. - If False, sets error arrays to 0 so that the additional pointing instances - are the same as the original ones. + If False, sets error arrays to 0 so that the additional + pointing instances are the same as the original ones. seed_pointing : integer, optional Seed for boresight pointing perturbation. Doesn't affect anything if perturb_el == False or perturb_az == False. mu_pointing : integer, optional - Mean for boresight error distribution in arcseconds. Doesn't affect anything - perturb_el == False or perturb_az == False. + Mean for boresight error distribution in arcseconds. + Doesn't affect anything perturb_el == False or perturb_az == False. sigma_pointing : integer, optional - Width for boresight error distribution in arcseconds. Doesn't affect anything - perturb_el == False or perturb_az == False. + Width for boresight error distribution in arcseconds. + Doesn't affect anything perturb_el == False or perturb_az == False. perturb_pol_angs : bool, optional If True, perturbing polarization angles. seed_pa : int, optional - seed for perturbing polarization angle (ch is added to it for randomization) + seed for perturbing polarization angle + (ch is added to it for randomization) pa_mu : float, optional mean of pa syst perturbation. pa_sig : float, optional @@ -147,7 +177,7 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, pa_diff : bool, optional if True, doing differential pa perturbation """ - ## Initialise args + # Initialise args self.verbose = verbose self.hardware = hardware self.scanning_strategy = scanning_strategy @@ -164,98 +194,122 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, self.perturb_el = perturb_el self.perturb_az = perturb_az - ## Check if you can run dichroic detectors + # Check if you can run dichroic detectors self.mode = mode - if self.mode == 'dichroic' and ( - not hasattr(self.HealpixFitsMap, 'I2')): + if self.mode == "dichroic" and ( + not hasattr(self.HealpixFitsMap, "I2") + ): raise IOError("You need two sets of maps for dichroic detectors!") - if self.mode == 'dichroic' and ( - not hasattr(self.hardware, 'focal_plane2')): + if self.mode == "dichroic" and ( + not hasattr(self.hardware, "focal_plane2") + ): raise IOError("You need two sets of det for dichroic detectors!") - if self.mode == 'dichroic' and ( - not hasattr(self.hardware, 'beam_model2')): + if self.mode == "dichroic" and ( + not hasattr(self.hardware, "beam_model2") + ): raise IOError("You need two sets of det for dichroic detectors!") self.width = width self.cut_pixels_outside = cut_pixels_outside self.store_pointing_matrix_input = store_pointing_matrix_input - ## Check the projection + # Check the projection self.projection = projection - assert self.projection in ['healpix', 'flat'], \ - ValueError("Projection <{}> for ".format(self.projection) + - "the output map not understood! " + - "Choose among ['healpix', 'flat'].") + assert self.projection in ["healpix", "flat"], ValueError( + """Projection <{}> for + the output map not understood! + Choose among ['healpix', 'flat'].""" + ).format(self.projection) - ## Check if the number of CES is consistent with the scanning strategy + # Check if the number of CES is consistent with the scanning strategy self.CESnumber = CESnumber - assert self.CESnumber < self.scanning_strategy.nces, \ - ValueError("The scan index must be between 0 and {}.".format( + assert self.CESnumber < self.scanning_strategy.nces, ValueError( + "The scan index must be between 0 and {}.".format( self.scanning_strategy.nces - 1 - )) - - ## Initialise internal parameters - self.scan = getattr(self.scanning_strategy, 'scan{}'.format( - self.CESnumber)) - self.nsamples = self.scan['nts'] + ) + ) + + # Initialise internal parameters + self.scan = getattr( + self.scanning_strategy, "scan{}".format(self.CESnumber) + ) + self.nsamples = self.scan["nts"] self.npair = self.hardware.focal_plane.npair self.pair_list = np.reshape( - self.hardware.focal_plane.bolo_index_in_fp, (self.npair, 2)) + self.hardware.focal_plane.bolo_index_in_fp, (self.npair, 2) + ) # boresight pointing perturbation errors if self.perturb_az or self.perturb_el: - mu_pointing = mu_pointing / 3600 * np.pi/180 # arcsec to radians - sigma_pointing = sigma_pointing / (3600*np.sqrt(2)) * np.pi/180 # arcsec to degrees and also dividing by sqrt(2) so that the overal error on position is np.sqrt(\delta_az^2+\delta_el^2) + mu_pointing = mu_pointing / 3600 * np.pi / 180 # arcsec to radians + # arcsec to degrees and also dividing by sqrt(2) so that + # the overal error on position is np.sqrt(\delta_az^2+\delta_el^2) + sigma_pointing = ( + sigma_pointing / (3600 * np.sqrt(2)) * np.pi / 180 + ) if self.perturb_az: # perturb azimuth seed_pointing_az = seed_pointing - state_for_pointing_errors_az = np.random.RandomState(seed_pointing_az) - self.err_azimuth = state_for_pointing_errors_az.normal(mu_pointing, sigma_pointing, self.scan['nts']) + state_for_pointing_errors_az = np.random.RandomState( + seed_pointing_az + ) + self.err_azimuth = state_for_pointing_errors_az.normal( + mu_pointing, sigma_pointing, self.scan["nts"] + ) else: - self.err_azimuth = np.zeros(self.scan['nts']) + self.err_azimuth = np.zeros(self.scan["nts"]) if self.perturb_el: # perturb elevation - seed_pointing_el = seed_pointing + 1234567890 # different seed for each perturbation; tried to avoid having another seed input - state_for_pointing_errors_el = np.random.RandomState(seed_pointing_el) - self.err_elevation = state_for_pointing_errors_el.normal(mu_pointing, sigma_pointing, self.scan['nts']) + seed_pointing_el = ( + seed_pointing + 1234567890 + ) + # different seed for each perturbation; + # tried to avoid having another seed input + state_for_pointing_errors_el = np.random.RandomState( + seed_pointing_el + ) + self.err_elevation = state_for_pointing_errors_el.normal( + mu_pointing, sigma_pointing, self.scan["nts"] + ) else: - self.err_elevation = np.zeros(self.scan['nts']) + self.err_elevation = np.zeros(self.scan["nts"]) - ## Pre-compute boresight pointing objects + # Pre-compute boresight pointing objects self.get_boresightpointing() - ## Polarisation angles: intrinsic and HWP angles + # Polarisation angles: intrinsic and HWP angles self.get_angles() - ## Position of bolometers in the focal plane - ## TODO move that elsewhere... + # Position of bolometers in the focal plane + # TODO move that elsewhere... self.ypos = self.hardware.beam_model.ypos self.xpos = self.hardware.beam_model.xpos self.xpos = self.xpos / np.cos(self.ypos) - ## Initialise pointing matrix, that is the matrix to go from time - ## to map domain, for all pairs of detectors. + # Initialise pointing matrix, that is the matrix to go from time + # to map domain, for all pairs of detectors. if not self.mapping_perpair: self.point_matrix = np.zeros( - (self.npair, self.nsamples), dtype=np.int32) + (self.npair, self.nsamples), dtype=np.int32 + ) else: self.point_matrix = np.zeros((1, self.nsamples), dtype=np.int32) - # If set, stores the pointing matrix used to scan the map if self.store_pointing_matrix_input: - self.point_matrix_input = np.zeros(self.point_matrix.shape, dtype=np.int32) - + self.point_matrix_input = np.zeros( + self.point_matrix.shape, dtype=np.int32 + ) - ## Initialise the mask for timestreams + # Initialise the mask for timestreams self.wafermask_pixel = self.get_timestream_masks() - ## Boundaries for subscans (t_beg, t_end) - self.subscans = self.scan['subscans'] + # Boundaries for subscans (t_beg, t_end) + self.subscans = self.scan["subscans"] - ## Get observed pixels in the input map + # Get observed pixels in the input map if nside_out is None: self.nside_out = self.HealpixFitsMap.nside else: @@ -268,15 +322,16 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, self.obspix, self.npixsky = self.get_obspix( self.width, self.scanning_strategy.ra_mid, - self.scanning_strategy.dec_mid) + self.scanning_strategy.dec_mid, + ) - ## Get timestream weights + # Get timestream weights self.sum_weight, self.diff_weight = self.get_weights() - ## Get detector gains + # Get detector gains self.set_detector_gains() - ## Prepare noise simulator if needed + # Prepare noise simulator if needed self.array_noise_level = array_noise_level self.array_noise_seed = array_noise_seed self.nclouds = nclouds @@ -287,13 +342,14 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, if self.array_noise_level is not None and self.alpha is None: self.noise_generator = WhiteNoiseGenerator( array_noise_level=self.array_noise_level, - ndetectors=2*self.npair, + ndetectors=2 * self.npair, ntimesamples=self.nsamples, - array_noise_seed=self.array_noise_seed) + array_noise_seed=self.array_noise_seed, + ) elif self.array_noise_level is not None and self.alpha is not None: self.noise_generator = CorrNoiseGenerator( array_noise_level=self.array_noise_level, - ndetectors=2*self.npair, + ndetectors=2 * self.npair, ntimesamples=self.nsamples, array_noise_seed=self.array_noise_seed, nclouds=self.nclouds, @@ -301,7 +357,8 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, amp_atm=self.amp_atm, corrlength=self.corrlength, alpha=self.alpha, - sampling_freq=self.scanning_strategy.sampling_freq) + sampling_freq=self.scanning_strategy.sampling_freq, + ) else: self.noise_generator = None @@ -310,13 +367,14 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, if self.array_noise_level2 is not None and self.alpha is None: self.noise_generator2 = WhiteNoiseGenerator( array_noise_level=self.array_noise_level2, - ndetectors=2*self.npair, + ndetectors=2 * self.npair, ntimesamples=self.nsamples, - array_noise_seed=self.array_noise_seed2) + array_noise_seed=self.array_noise_seed2, + ) elif self.array_noise_level2 is not None and self.alpha is not None: self.noise_generator2 = CorrNoiseGenerator( array_noise_level=self.array_noise_level2, - ndetectors=2*self.npair, + ndetectors=2 * self.npair, ntimesamples=self.nsamples, array_noise_seed=self.array_noise_seed2, nclouds=self.nclouds, @@ -324,7 +382,8 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, amp_atm=self.amp_atm, corrlength=self.corrlength, alpha=self.alpha, - sampling_freq=self.scanning_strategy.sampling_freq) + sampling_freq=self.scanning_strategy.sampling_freq, + ) else: self.noise_generator2 = None @@ -334,25 +393,25 @@ def get_angles(self): and initialise total polarisation angle. """ self.hwpangle = self.hardware.half_wave_plate.compute_HWP_angles( - sample_rate=self.scan['sample_rate'], - size=self.nsamples) + sample_rate=self.scan["sample_rate"], size=self.nsamples + ) self.intrinsic_polangle = self.hardware.focal_plane.bolo_polangle self.intrinsic_polangle2 = None - if self.mode == 'dichroic': + if self.mode == "dichroic": self.intrinsic_polangle2 = self.hardware.focal_plane2.bolo_polangle - ## Will contain the total polarisation angles for all bolometers - ## That is PA + intrinsic + 2 * HWP + # Will contain the total polarisation angles for all bolometers + # That is PA + intrinsic + 2 * HWP if not self.mapping_perpair: self.pol_angs = np.zeros((self.npair, self.nsamples)) self.pol_angs2 = None - if self.mode == 'dichroic': + if self.mode == "dichroic": self.pol_angs2 = np.zeros((self.npair, self.nsamples)) else: self.pol_angs = np.zeros((1, self.nsamples)) self.pol_angs2 = None - if self.mode == 'dichroic': + if self.mode == "dichroic": self.pol_angs2 = np.zeros((1, self.nsamples)) def get_timestream_masks(self): @@ -406,48 +465,51 @@ def get_obspix(self, width, ra_src, dec_src): >>> tod = TimeOrderedDataPairDiff(inst, scan, sky_in, CESnumber=0, ... projection='flat') >>> obspix, npix = tod.get_obspix(10., 0., 0.) - >>> print(obspix) ## the same as healpix + >>> print(obspix) # the same as healpix [1376 1439 1440 1504 1567 1568 1632 1695] >>> print(npix, len(obspix)) 16 8 """ - ## Change to radian + # Change to radian ra_src = ra_src * d2r dec_src = dec_src * d2r - ## TODO implement the first line correctly... + # TODO implement the first line correctly... try: xmin, xmax, ymin, ymax = np.array(width) * d2r except TypeError: - xmin = xmax = ymin = ymax = d2r * width / 2. + xmin = xmax = ymin = ymax = d2r * width / 2.0 # If map bound crosses zero make coordinates - ## bounds monotonic in [-pi,pi] - ra_min = (ra_src - xmin) + # bounds monotonic in [-pi,pi] + ra_min = ra_src - xmin if (ra_src + xmax) >= 2 * np.pi: ra_max = (ra_src + xmax) % (2 * np.pi) ra_min = ra_min if ra_min <= np.pi else ra_min - 2 * np.pi else: - ra_min = (ra_src - xmin) - ra_max = (ra_src + xmax) + ra_min = ra_src - xmin + ra_max = ra_src + xmax - dec_min = max([(dec_src - ymin), -np.pi/2]) - dec_max = min([dec_src + ymax, np.pi/2]) + dec_min = max([(dec_src - ymin), -np.pi / 2]) + dec_max = min([dec_src + ymax, np.pi / 2]) self.xmin = ra_min self.xmax = ra_max self.ymin = dec_min self.ymax = dec_max - obspix = input_sky.get_obspix(ra_min, ra_max, - dec_min, dec_max, - self.nside_out) - - if self.projection == 'flat': - npixsky = int( - round( - (self.xmax - self.xmin + self.pixel_size) / - self.pixel_size))**2 - elif self.projection == 'healpix': + obspix = input_sky.get_obspix( + ra_min, ra_max, dec_min, dec_max, self.nside_out + ) + + if self.projection == "flat": + npixsky = ( + int( + round( + (self.xmax - self.xmin + self.pixel_size) / self.pixel_size + ) + ) ** 2 + ) + elif self.projection == "healpix": npixsky = len(obspix) return obspix, npixsky @@ -512,20 +574,23 @@ def set_detector_gains(self, new_gains=None, new_gains2=None): >>> tod.set_detector_gains(new_gains=new_gains, new_gains2=new_gains2) """ if new_gains is not None: - assert len(new_gains) == 2 * self.npair, \ - ValueError("You have to provide {} new gain values!".format( - 2 * self.npair)) + assert len(new_gains) == 2 * self.npair, ValueError( + "You have to provide {} new gain values!".format( + 2 * self.npair + ) + ) self.gain = new_gains else: self.gain = np.ones(2 * self.npair) self.gain2 = None - if self.mode == 'dichroic': + if self.mode == "dichroic": if new_gains2 is not None: - assert len(new_gains2) == 2 * self.npair, \ - ValueError( - "You have to provide {} new gain values!".format( - 2 * self.npair)) + assert len(new_gains2) == 2 * self.npair, ValueError( + "You have to provide {} new gain values!".format( + 2 * self.npair + ) + ) self.gain2 = new_gains2 else: self.gain2 = np.ones(2 * self.npair) @@ -574,18 +639,20 @@ def set_detector_gains_perpair(self, new_gains=None, new_gains2=None): """ if new_gains is not None: msg = "You have to provide ({}, {}) new gain values!" - assert new_gains.shape == (2, self.nsamples), \ - ValueError(msg.format(2, self.nsamples)) + assert new_gains.shape == (2, self.nsamples), ValueError( + msg.format(2, self.nsamples) + ) self.gain = new_gains else: self.gain = np.ones((2, self.nsamples)) self.gain2 = None - if self.mode == 'dichroic': + if self.mode == "dichroic": if new_gains2 is not None: msg = "You have to provide ({}, {}) new gain values!" - assert new_gains2.shape == (2, self.nsamples), \ - ValueError(msg.format(2, self.nsamples)) + assert new_gains2.shape == (2, self.nsamples), ValueError( + msg.format(2, self.nsamples) + ) self.gain2 = new_gains2 else: self.gain2 = np.ones((2, self.nsamples)) @@ -602,48 +669,58 @@ def get_boresightpointing(self): This is to avoid projection artifact by operating a rotation of the coordinates to (0, 0) in flat projection (scan around equator). """ - lat = float( - self.scanning_strategy.telescope_location.lat) * 180. / np.pi + lat = float(self.scanning_strategy.telescope_location.lat)\ + * 180.0 / np.pi - if self.projection == 'healpix': + if self.projection == "healpix": ra_src = 0.0 dec_src = 0.0 - elif self.projection == 'flat': - ra_src = self.scanning_strategy.ra_mid * np.pi / 180. - dec_src = self.scanning_strategy.dec_mid * np.pi / 180. + elif self.projection == "flat": + ra_src = self.scanning_strategy.ra_mid * np.pi / 180.0 + dec_src = self.scanning_strategy.dec_mid * np.pi / 180.0 - ## Perform a rotation of the input to put the point - ## (ra_src, dec_src) at (0, 0). + # Perform a rotation of the input to put the point + # (ra_src, dec_src) at (0, 0). # r = hp.Rotator(rot=[ra_src, self.scanning_strategy.dec_mid]) # theta, phi = hp.pix2ang(self.HealpixFitsMap.nside, # range(12 * self.HealpixFitsMap.nside**2)) # t, p = r(theta, phi, inv=True) # pix = hp.ang2pix(self.HealpixFitsMap.nside, t, p) # - # ## Apply the rotation to our maps + # # Apply the rotation to our maps # self.HealpixFitsMap.I = self.HealpixFitsMap.I[pix] # self.HealpixFitsMap.Q = self.HealpixFitsMap.Q[pix] # self.HealpixFitsMap.U = self.HealpixFitsMap.U[pix] self.pointing = Pointing( - az_enc=self.scan['azimuth'], # degrees, size of nts - el_enc=self.scan['elevation'], # degrees, size of nts - time=self.scan['clock-utc'], + az_enc=self.scan["azimuth"], # degrees, size of nts + el_enc=self.scan["elevation"], # degrees, size of nts + time=self.scan["clock-utc"], value_params=self.hardware.pointing_model.value_params, allowed_params=self.hardware.pointing_model.allowed_params, ut1utc_fn=self.scanning_strategy.ut1utc_fn, - lat=lat, ra_src=ra_src, dec_src=dec_src) + lat=lat, + ra_src=ra_src, + dec_src=dec_src, + ) # boresight pointing systematics if self.perturb_el or self.perturb_az: self.pointing_perturbed = Pointing( - az_enc = (self.scan['azimuth']+self.err_azimuth), # degrees, size of nts - el_enc = (self.scan['elevation']+self.err_elevation), # degrees, size of nts - time=self.scan['clock-utc'], + az_enc=( + self.scan["azimuth"] + self.err_azimuth + ), # degrees, size of nts + el_enc=( + self.scan["elevation"] + self.err_elevation + ), # degrees, size of nts + time=self.scan["clock-utc"], value_params=self.hardware.pointing_model.value_params, allowed_params=self.hardware.pointing_model.allowed_params, ut1utc_fn=self.scanning_strategy.ut1utc_fn, - lat=lat, ra_src=ra_src, dec_src=dec_src) + lat=lat, + ra_src=ra_src, + dec_src=dec_src, + ) else: self.pointing_perturbed = None @@ -693,11 +770,14 @@ def compute_simpolangle(self, ch, parallactic_angle, polangle_err=False): if polangle_err: state_for_polang = np.random.RandomState(self.seed_pa) - - if self.pa_sig==0: - dpolang = np.ones(len(self.intrinsic_polangle))*self.pa_mu # in deg + if self.pa_sig == 0: + dpolang = ( + np.ones(len(self.intrinsic_polangle)) * self.pa_mu + ) # in deg else: - dpolang = state_for_polang.normal(self.pa_mu, self.pa_sig,len(self.intrinsic_polangle)) # in deg + dpolang = state_for_polang.normal( + self.pa_mu, self.pa_sig, len(self.intrinsic_polangle) + ) # in deg if self.pa_diff: # differential pol ang symmetrically @@ -710,29 +790,28 @@ def compute_simpolangle(self, ch, parallactic_angle, polangle_err=False): else: intrinsic_polangle = self.intrinsic_polangle - # calc total pa (with or without perturbations) ang_pix = (90.0 - intrinsic_polangle[ch]) * d2r - ## Demodulation or pair diff use different convention - ## for the definition of the angle. - if not hasattr(self, 'dm'): + # Demodulation or pair diff use different convention + # for the definition of the angle. + if not hasattr(self, "dm"): pol_ang = parallactic_angle + ang_pix + 2.0 * self.hwpangle else: pol_ang = parallactic_angle - ang_pix - 2.0 * self.hwpangle pol_ang2 = None # TODO: add polarization angle perturbations for the second frequency. - if self.mode == 'dichroic': + if self.mode == "dichroic": ang_pix2 = (90.0 - self.intrinsic_polangle2[ch]) * d2r - ## Demodulation or pair diff use different convention - ## for the definition of the angle. - if not hasattr(self, 'dm'): - pol_ang2 = parallactic_angle + ang_pix2 + 2.0*self.hwpangle + # Demodulation or pair diff use different convention + # for the definition of the angle. + if not hasattr(self, "dm"): + pol_ang2 = parallactic_angle + ang_pix2 + 2.0 * self.hwpangle else: - pol_ang2 = parallactic_angle - ang_pix2 - 2.0*self.hwpangle + pol_ang2 = parallactic_angle - ang_pix2 - 2.0 * self.hwpangle return pol_ang, pol_ang2 @@ -778,44 +857,50 @@ def return_parallactic_angle(self, ch, frequency_channel=1): elif frequency_channel == 2: ang_pix = (90.0 - self.intrinsic_polangle2[ch]) * d2r - if not hasattr(self, 'dm'): - return pa - ang_pix - 2*self.hwpangle + if not hasattr(self, "dm"): + return pa - ang_pix - 2 * self.hwpangle else: return pa + ang_pix else: if frequency_channel == 1: - ang_pix = (90.0 - self.intrinsic_polangle[2*ch]) * d2r + ang_pix = (90.0 - self.intrinsic_polangle[2 * ch]) * d2r elif frequency_channel == 2: - ang_pix = (90.0 - self.intrinsic_polangle2[2*ch]) * d2r - if not hasattr(self, 'dm'): - return pa[ch] - ang_pix - 2*self.hwpangle + ang_pix = (90.0 - self.intrinsic_polangle2[2 * ch]) * d2r + if not hasattr(self, "dm"): + return pa[ch] - ang_pix - 2 * self.hwpangle else: return pa[ch] + ang_pix - def get_pixel_indices(self,ra,dec): - ## Retrieve corresponding pixels on the sky, and their index locally. - if self.projection == 'flat': - ## ?? - xmin = - self.width/2.*np.pi/180. - ymin = - self.width/2.*np.pi/180. + def get_pixel_indices(self, ra, dec): + # Retrieve corresponding pixels on the sky, and their index locally. + if self.projection == "flat": + # ?? + xmin = -self.width / 2.0 * np.pi / 180.0 + ymin = -self.width / 2.0 * np.pi / 180.0 index_global, index_local = build_pointing_matrix( - ra, dec, nside_in=self.HealpixFitsMap.nside, + ra, + dec, + nside_in=self.HealpixFitsMap.nside, nside_out=self.nside_out, xmin=xmin, ymin=ymin, pixel_size=self.pixel_size, npix_per_row=int(np.sqrt(self.npixsky)), projection=self.projection, - cut_pixels_outside=self.cut_pixels_outside) - elif self.projection == 'healpix': + cut_pixels_outside=self.cut_pixels_outside, + ) + elif self.projection == "healpix": index_global, index_local = build_pointing_matrix( - ra, dec, nside_in=self.HealpixFitsMap.nside, + ra, + dec, + nside_in=self.HealpixFitsMap.nside, nside_out=self.nside_out, obspix=self.obspix, ext_map_gal=self.HealpixFitsMap.ext_map_gal, projection=self.projection, - cut_pixels_outside=self.cut_pixels_outside) - return index_global,index_local + cut_pixels_outside=self.cut_pixels_outside, + ) + return index_global, index_local def map2tod(self, ch): """ @@ -860,21 +945,28 @@ def map2tod(self, ch): ... mode='dichroic', CESnumber=1) >>> d = tod.map2tod(0) """ - ## Use bolometer beam offsets. + # Use bolometer beam offsets. azd, eld = self.xpos[ch], self.ypos[ch] - ## Compute pointing for detector ch + # Compute pointing for detector ch ra, dec, pa = self.pointing.offset_detector(azd, eld) - ## Compute pointing matrix + # Compute pointing matrix index_global, index_local = self.get_pixel_indices(ra, dec) # Perturbed boresight pointing values if self.pointing_perturbed is not None: - ra_perturbed, dec_perturbed, pa_perturbed = self.pointing_perturbed.offset_detector(azd, eld) + ( + ra_perturbed, + dec_perturbed, + pa_perturbed, + ) = self.pointing_perturbed.offset_detector(azd, eld) # Perturbed boresight pointing values - index_global_perturbed, index_local_perturbed = self.get_pixel_indices(ra_perturbed,dec_perturbed) + ( + index_global_perturbed, + index_local_perturbed, + ) = self.get_pixel_indices(ra_perturbed, dec_perturbed) else: ra_perturbed = ra @@ -883,32 +975,35 @@ def map2tod(self, ch): index_global_perturbed = index_global index_local_perturbed = index_local - if ((self.projection == 'healpix') & (index_local is None)) : + if (self.projection == "healpix") & (index_local is None): # Using a pointer not to increase memory usage index_local = index_global # Perturbed boresight pointing values index_local_perturbed = index_global_perturbed - - ## For flat projection, one needs to flip the sign of U - ## w.r.t to the full-sky basis (IAU angle convention) - if self.projection == 'flat': - sign = -1. - elif self.projection == 'healpix': sign = 1. + # For flat projection, one needs to flip the sign of U + # w.r.t to the full-sky basis (IAU angle convention) + if self.projection == "flat": + sign = -1.0 + elif self.projection == "healpix": + sign = 1.0 self.Usign = sign - ## Store pointing matrix and parallactic angle for tod2map operations + # Store pointing matrix and parallactic angle for tod2map operations if ch % 2 == 0: - if ((self.xpos[ch]!=self.xpos[ch+1]) or (self.ypos[ch]!=self.ypos[ch+1])): - ## Compute pointing matrix for the pair center if beam - ## offsets for top and bottom detector do not coincide. - azd = 0.5*(self.xpos[ch]+self.xpos[ch+1]) - eld = 0.5*(self.ypos[ch]+self.ypos[ch+1]) + if (self.xpos[ch] != self.xpos[ch + 1]) or ( + self.ypos[ch] != self.ypos[ch + 1] + ): + # Compute pointing matrix for the pair center if beam + # offsets for top and bottom detector do not coincide. + azd = 0.5 * (self.xpos[ch] + self.xpos[ch + 1]) + eld = 0.5 * (self.ypos[ch] + self.ypos[ch + 1]) ra, dec, pa_pair = self.pointing.offset_detector(azd, eld) - index_global_pair, index_local = self.get_pixel_indices(ra,dec) - #print("even",ch, self.xpos[ch],self.xpos[ch+1],self.ypos[ch],self.ypos[ch+1],pa,pa_pair) - if index_local is None : + index_global_pair, index_local = self.get_pixel_indices( + ra, dec + ) + if index_local is None: # Using a pointer not to increase memory usage index_local = index_global_pair else: @@ -916,23 +1011,22 @@ def map2tod(self, ch): # introducing differential beam systematics. index_global_pair = index_global - ## Store list of pixels to be mapped only for top bolometers or pair - ## center in presence of differntial pointing + # Store list of pixels to be mapped only for top bolometers or pair + # center in presence of differntial pointing if ch % 2 == 0 and not self.mapping_perpair: - self.point_matrix[int(ch/2)] = index_local + self.point_matrix[int(ch / 2)] = index_local if self.store_pointing_matrix_input: - index_global_pair[index_local==-1] = -1 - self.point_matrix_input[int(ch/2)] = index_global_pair + index_global_pair[index_local == -1] = -1 + self.point_matrix_input[int(ch / 2)] = index_global_pair elif ch % 2 == 0 and self.mapping_perpair: self.point_matrix[0] = index_local if self.store_pointing_matrix_input: - index_global_pair[index_local==-1] = -1 + index_global_pair[index_local == -1] = -1 self.point_matrix_input[0] = index_global_pair - - ## Default gain for a detector is 1., - ## but you can change it using set_detector_gains or - ## set_detector_gains_perpair. + # Default gain for a detector is 1., + # but you can change it using set_detector_gains or + # set_detector_gains_perpair. if len(self.gain) == self.npair * 2: norm = self.gain[ch] elif len(self.gain) == 2: @@ -942,7 +1036,7 @@ def map2tod(self, ch): else: norm = self.gain[1] - ## Noise simulation + # Noise simulation if self.noise_generator is not None: noise = self.noise_generator.simulate_noise_one_detector(ch) else: @@ -954,81 +1048,93 @@ def map2tod(self, ch): noise2 = 0.0 if self.HealpixFitsMap.do_pol: - ## pol_ang2 is None if mode == 'standard' - pol_ang, pol_ang2 = self.compute_simpolangle(ch, pa, polangle_err=False) + # pol_ang2 is None if mode == 'standard' + pol_ang, pol_ang2 = self.compute_simpolangle( + ch, pa, polangle_err=False + ) if self.perturb_pol_angs: do_polangle_err = True else: do_polangle_err = False - pol_ang_perturbed, pol_ang2_perturbed = self.compute_simpolangle(ch, pa, polangle_err=do_polangle_err) - + pol_ang_perturbed, pol_ang2_perturbed = self.compute_simpolangle( + ch, pa, polangle_err=do_polangle_err + ) - ## Use pa of pair center for tod2map if differential pointing - ## is used. Otherwise the polarization angle of a pair is the - ## same as the top and bottom detectors. + # Use pa of pair center for tod2map if differential pointing + # is used. Otherwise the polarization angle of a pair is the + # same as the top and bottom detectors. try: pol_ang_pair, pol_ang2_pair = self.compute_simpolangle( - ch, pa_pair, polangle_err=False) - except: + ch, pa_pair, polangle_err=False + ) + except Exception as e: + print(str(e)) pol_ang_pair = pol_ang pol_ang2_pair = pol_ang2 - ## For demodulation, HWP angles are not included at the level - ## of the pointing matrix (convention). - if hasattr(self, 'dm'): + # For demodulation, HWP angles are not included at the level + # of the pointing matrix (convention). + if hasattr(self, "dm"): pol_ang_out = pol_ang_pair + 2.0 * self.hwpangle else: pol_ang_out = pol_ang_pair - ## Store list of polangle only for top bolometers + # Store list of polangle only for top bolometers if ch % 2 == 0 and not self.mapping_perpair: - self.pol_angs[int(ch/2)] = pol_ang_out + self.pol_angs[int(ch / 2)] = pol_ang_out elif ch % 2 == 0 and self.mapping_perpair: self.pol_angs[0] = pol_ang_out - ts1 = ( - self.HealpixFitsMap.I[index_global] + - self.HealpixFitsMap.Q[index_global] * np.cos(2 * pol_ang) + - sign * self.HealpixFitsMap.U[index_global] * - np.sin(2 * pol_ang) + noise) * norm + ts1 = self.HealpixFitsMap.I[index_global]\ + + self.HealpixFitsMap.Q[index_global] * np.cos(2 * pol_ang)\ + + sign * self.HealpixFitsMap.U[index_global] * np.sin(2 * pol_ang)\ + + noise + ts1 = ts1 * norm - if self.mode == 'standard': + if self.mode == "standard": return ts1 - elif self.mode == 'dichroic': + elif self.mode == "dichroic": # For demodulation, HWP angles are not included at the level - ## of the pointing matrix (convention). - if hasattr(self, 'dm'): + # of the pointing matrix (convention). + if hasattr(self, "dm"): pol_ang_out2 = pol_ang2_pair + 2.0 * self.hwpangle else: pol_ang_out2 = pol_ang2_pair - ## Store list polangle only for top bolometers + # Store list polangle only for top bolometers if ch % 2 == 0 and not self.mapping_perpair: - self.pol_angs2[int(ch/2)] = pol_ang_out2 + self.pol_angs2[int(ch / 2)] = pol_ang_out2 elif ch % 2 == 0 and self.mapping_perpair: self.pol_angs2[0] = pol_ang_out2 - ts2 = ( - self.HealpixFitsMap.I2[index_global_perturbed] + - self.HealpixFitsMap.Q2[index_global_perturbed] * np.cos(2*pol_ang2_perturbed) + - sign * self.HealpixFitsMap.U2[index_global_perturbed] * - np.sin(2 * pol_ang2_perturbed) + noise2) * norm + ts2 = self.HealpixFitsMap.I2[index_global_perturbed]\ + + self.HealpixFitsMap.Q2[index_global_perturbed]\ + * np.cos(2 * pol_ang2_perturbed)\ + + sign\ + * self.HealpixFitsMap.U2[index_global_perturbed]\ + * np.sin(2 * pol_ang2_perturbed)\ + + noise2 + ts2 = ts2 * norm return np.array([ts1, ts2]) else: - ts1 = norm * (self.HealpixFitsMap.I[index_global_perturbed] + noise) - if self.mode == 'standard': + ts1 = norm * ( + self.HealpixFitsMap.I[index_global_perturbed] + noise + ) + if self.mode == "standard": return ts1 - elif self.mode == 'dichroic': - ts2 = norm * (self.HealpixFitsMap.I2[index_global_perturbed] + noise2) + elif self.mode == "dichroic": + ts2 = norm * ( + self.HealpixFitsMap.I2[index_global_perturbed] + noise2 + ) return np.array([ts1, ts2]) - def tod2map(self, waferts, output_maps, - gdeprojection=False, - frequency_channel=1): + def tod2map( + self, waferts, output_maps, gdeprojection=False, frequency_channel=1 + ): """ Project time-ordered data into sky maps for the whole array. Maps are updated on-the-fly. Massive speed-up thanks to the @@ -1104,8 +1210,8 @@ def tod2map(self, waferts, output_maps, ... nside=tod.nside_out, obspix=tod.obspix) >>> for pair in tod.pair_list: ... d = np.array([tod.map2tod(det) for det in pair]) - ... d1 = d[:, 0] ## first frequency channel - ... d2 = d[:, 1] ## second frequency channel + ... d1 = d[:, 0] # first frequency channel + ... d2 = d[:, 1] # second frequency channel ... tod.tod2map(d1, m1) ... tod.tod2map(d2, m2) @@ -1129,17 +1235,19 @@ def tod2map(self, waferts, output_maps, npixfp = nbolofp / 2 nt = int(waferts.shape[-1]) - ## Check sizes - msg = 'Most likely you set mapping_perpair wrongly when ' + \ - 'initialising your TOD.' + \ - 'The way you loop over focal plane pixel to do the mapmaking' + \ - 'depends on the value of mapping_perpair parameter. ' + \ - 'If mapping_perpair=False, one first load all the pixel and ' + \ - 'then you need to perform the mapmaking with all the pixels ' + \ - 'at once. If mapping_perpair=True, one should load ' + \ - 'pair-by-pair and the mapmaking is done pair-by-pair.' + \ - 'See so_MC_crosstalk.py vs so_MC_gain_drift.py to see both ' + \ - 'approaches (s4cmb-resources/Part2), and example in doctest above.' + # Check sizes + msg = """ + Most likely you set mapping_perpair wrongly when + initialising your TOD. + The way you loop over focal plane pixel to do the mapmaking + depends on the value of mapping_perpair parameter. + If mapping_perpair=False, one first load all the pixel and + then you need to perform the mapmaking with all the pixels + at once. If mapping_perpair=True, one should load + pair-by-pair and the mapmaking is done pair-by-pair. + See so_MC_crosstalk.py vs so_MC_gain_drift.py to see both + approaches (s4cmb-resources/Part2), and example in doctest above. + """ assert npixfp == self.point_matrix.shape[0], msg assert nt == self.point_matrix.shape[1], msg @@ -1156,44 +1264,95 @@ def tod2map(self, waferts, output_maps, sum_weight = self.sum_weight.flatten() wafermask_pixel = self.wafermask_pixel.flatten() - if (hasattr(self, 'dm') and (gdeprojection is False)): + if hasattr(self, "dm") and (gdeprojection is False): tod_f.tod2map_hwp_f( - output_maps.d0, output_maps.d4r, output_maps.d4i, - output_maps.w0, output_maps.w4, output_maps.nhit, - point_matrix, pol_angs, waferts, - diff_weight, sum_weight, nt, - wafermask_pixel, npixfp, self.npixsky) - elif ((hasattr(output_maps, 'dm') is False) and gdeprojection): - print('gdep!') + output_maps.d0, + output_maps.d4r, + output_maps.d4i, + output_maps.w0, + output_maps.w4, + output_maps.nhit, + point_matrix, + pol_angs, + waferts, + diff_weight, + sum_weight, + nt, + wafermask_pixel, + npixfp, + self.npixsky, + ) + elif (hasattr(output_maps, "dm") is False) and gdeprojection: + print("gdep!") tod_f.tod2map_pair_gdeprojection_f( - output_maps.d, output_maps.w, - output_maps.dm, output_maps.dc, output_maps.ds, - output_maps.wm, output_maps.cc, output_maps.cs, - output_maps.ss, output_maps.c, output_maps.s, output_maps.nhit, - point_matrix, pol_angs, waferts, - diff_weight, sum_weight, nt, - wafermask_pixel, npixfp, self.npixsky) + output_maps.d, + output_maps.w, + output_maps.dm, + output_maps.dc, + output_maps.ds, + output_maps.wm, + output_maps.cc, + output_maps.cs, + output_maps.ss, + output_maps.c, + output_maps.s, + output_maps.nhit, + point_matrix, + pol_angs, + waferts, + diff_weight, + sum_weight, + nt, + wafermask_pixel, + npixfp, + self.npixsky, + ) else: tod_f.tod2map_pair_f( - output_maps.d, output_maps.w, output_maps.dc, - output_maps.ds, output_maps.cc, output_maps.cs, - output_maps.ss, output_maps.nhit, - point_matrix, pol_angs, waferts, - diff_weight, sum_weight, nt, - wafermask_pixel, npixfp, self.npixsky) + output_maps.d, + output_maps.w, + output_maps.dc, + output_maps.ds, + output_maps.cc, + output_maps.cs, + output_maps.ss, + output_maps.nhit, + point_matrix, + pol_angs, + waferts, + diff_weight, + sum_weight, + nt, + wafermask_pixel, + npixfp, + self.npixsky, + ) # Garbage collector guard wafermask_pixel class TimeOrderedDataDemod(TimeOrderedDataPairDiff): """ Class to """ - def __init__(self, hardware, scanning_strategy, HealpixFitsMap, - CESnumber, projection='healpix', - nside_out=None, pixel_size=None, width=140., - cut_pixels_outside=True, - array_noise_level=None, array_noise_seed=487587, - array_noise_level2=None, array_noise_seed2=56736, - mapping_perpair=False, mode='standard', verbose=False): + + def __init__( + self, + hardware, + scanning_strategy, + HealpixFitsMap, + CESnumber, + projection="healpix", + nside_out=None, + pixel_size=None, + width=140.0, + cut_pixels_outside=True, + array_noise_level=None, + array_noise_seed=487587, + array_noise_level2=None, + array_noise_seed2=56736, + mapping_perpair=False, + mode="standard", + verbose=False, + ): """ C'est parti! @@ -1271,9 +1430,15 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, >>> tod.tod2map(d2, m2, frequency_channel=2) """ TimeOrderedDataPairDiff.__init__( - self, hardware, scanning_strategy, HealpixFitsMap, - CESnumber, projection=projection, - nside_out=nside_out, pixel_size=pixel_size, width=width, + self, + hardware, + scanning_strategy, + HealpixFitsMap, + CESnumber, + projection=projection, + nside_out=nside_out, + pixel_size=pixel_size, + width=width, cut_pixels_outside=cut_pixels_outside, array_noise_level=array_noise_level, array_noise_seed=array_noise_seed, @@ -1281,16 +1446,18 @@ def __init__(self, hardware, scanning_strategy, HealpixFitsMap, array_noise_seed2=array_noise_seed2, mapping_perpair=mapping_perpair, mode=mode, - verbose=verbose) + verbose=verbose, + ) - ## Prepare the demodulation of timestreams + # Prepare the demodulation of timestreams self.dm = Demodulation( self.hardware.half_wave_plate.freq_hwp, self.scanning_strategy.sampling_freq, self.hwpangle, - self.verbose) + self.verbose, + ) - ## Prepare the filters use for the demodulation + # Prepare the filters use for the demodulation self.dm.prepfilter([4], [1.9]) self.dm.prepfftedfilter(nt=self.nsamples) @@ -1320,11 +1487,11 @@ def demodulate_timestreams(self, ts): # dm.b.copy() self.dm.br = ts - ## Do temperature + # Do temperature self.dm.demod(0) newts[:, 0, :] = self.dm.bm.real - ## Do 4f component (effectively polarisation) + # Do 4f component (effectively polarisation) self.dm.demod(4) newts[:, 1, :] = self.dm.bm.real newts[:, 2, :] = self.dm.bm.imag @@ -1332,8 +1499,9 @@ def demodulate_timestreams(self, ts): return newts -class Demodulation(): +class Demodulation: """ Class to handle demodulation of timestreams """ + def __init__(self, hwp_freq, sampling_freq, hwp_angles, verbose=False): """ Demodulation of timestreams means here that we estimate I, Q and U @@ -1363,12 +1531,12 @@ def __init__(self, hwp_freq, sampling_freq, hwp_angles, verbose=False): self.hwp_angles = hwp_angles self.verbose = verbose - ## In Hz - self.speed = np.median( - np.diff(self.hwp_angles)) * self.sampling_freq / (2 * np.pi) + # In Hz + med = np.median(np.diff(self.hwp_angles)) + self.speed = med * self.sampling_freq / (2 * np.pi) - ## Nyquist frequency at half the sampling - self.nyq = self.sampling_freq / 2. + # Nyquist frequency at half the sampling + self.nyq = self.sampling_freq / 2.0 def prepfilter(self, modes, bands=None, numtaps=None, relative=True): """ @@ -1394,7 +1562,7 @@ def prepfilter(self, modes, bands=None, numtaps=None, relative=True): """ if numtaps is None: - if self.nyq < 100./6: + if self.nyq < 100.0 / 6: numtaps = 255 else: numtaps = 1023 @@ -1412,24 +1580,26 @@ def prepfilter(self, modes, bands=None, numtaps=None, relative=True): self.bands *= self.speed self.hwp_freq *= self.speed - bad = (self.modes * self.speed + self.bands > self.nyq) + bad = self.modes * self.speed + self.bands > self.nyq if bad.any(): - print('You need freq = modes * speed + bands < nyquist') - print('MODES', self.modes) - print('SPEED', self.speed) - print('BANDS', self.bands) - print('FREQ', self.modes * self.speed + self.bands) - print('NYQUIST', self.nyq) - print('mode', self.modes[bad], 'are over nyq.') - print('You need either to increase the sampling frequency of', - 'your detectors, or decrease the spin frequency of the HWP.') + print("You need freq = modes * speed + bands < nyquist") + print("MODES", self.modes) + print("SPEED", self.speed) + print("BANDS", self.bands) + print("FREQ", self.modes * self.speed + self.bands) + print("NYQUIST", self.nyq) + print("mode", self.modes[bad], "are over nyq.") + print( + "You need either to increase the sampling frequency of", + "your detectors, or decrease the spin frequency of the HWP.", + ) self.modes = self.modes[~bad] self.bands = self.bands[~bad] self.lpf0 = firwin(numtaps, self.hwp_freq, nyq=self.nyq) - ## Construct the filters + # Construct the filters self.lpfs = [] self.bpfs = [] for mode, band in zip(self.modes, self.bands): @@ -1438,10 +1608,11 @@ def prepfilter(self, modes, bands=None, numtaps=None, relative=True): numtaps, [self.speed * mode - band, self.speed * mode + band], nyq=self.nyq, - pass_zero=False) + pass_zero=False, + ) fbpf = fftpack.fft(bpfreal) - fbpf[: int((numtaps + 1) / 2)] = 0. + fbpf[: int((numtaps + 1) / 2)] = 0.0 self.bpfs.append(fftpack.ifft(fbpf)) self.bpfs = np.array(self.bpfs) @@ -1488,14 +1659,14 @@ def demod(self, mode, bpf=True, lpf=True): lpf : bool, optional If True, perform a Low Pass Filtering. Default is True. """ - ## Temperature or not understood! + # Temperature or not understood! if mode not in self.modes: if mode == 0: # if (self.br == 0.).all(): # self.bm = np.zeros(self.br.shape, dtype=self.br.dtype) # return self.bm - if hasattr(self, 'flpf0'): + if hasattr(self, "flpf0"): flpf0 = self.flpf0 else: flpf0 = None @@ -1503,15 +1674,15 @@ def demod(self, mode, bpf=True, lpf=True): self.bm = convolvefilter(self.br, self.lpf0, flpf0) return self.bm else: - print('Filters for', mode, 'is not prepared!') + print("Filters for", mode, "is not prepared!") return None - ## Polarisation + # Polarisation else: i = np.where(self.modes == mode)[0][0] - if hasattr(self, 'bm'): - del(self.bm) + if hasattr(self, "bm"): + del self.bm # if (self.br == 0.).all(): # self.bm = np.zeros(self.br.shape, dtype=self.br.dtype) + 0.j @@ -1521,7 +1692,7 @@ def demod(self, mode, bpf=True, lpf=True): bpf = self.bpfs[i] try: fbpf = self.fbpfs[i] - except(AttributeError, IndexError): + except (AttributeError, IndexError): fbpf = None else: bpf = None @@ -1529,13 +1700,19 @@ def demod(self, mode, bpf=True, lpf=True): lpf = self.lpfs[i] try: flpf = self.flpfs[i] - except(AttributeError, IndexError): + except (AttributeError, IndexError): flpf = None else: lpf = None - self.bm = demod(self.br, np.exp(1.j * mode * self.hwp_angles), - bpf=bpf, lpf=lpf, fbpf=fbpf, flpf=flpf) + self.bm = demod( + self.br, + np.exp(1.0j * mode * self.hwp_angles), + bpf=bpf, + lpf=lpf, + fbpf=fbpf, + flpf=flpf, + ) return self.bm @@ -1571,7 +1748,7 @@ def demod(x, e, bpf=None, lpf=None, fbpf=None, flpf=None): For examples, see TimeOrderedDataDemod. """ - u = np.ones(list(x.shape[:-1]) + [1], x.dtype) * 2. * e + u = np.ones(list(x.shape[:-1]) + [1], x.dtype) * 2.0 * e if bpf is not None: u *= convolvefilter(x, bpf, fbpf) else: @@ -1580,6 +1757,7 @@ def demod(x, e, bpf=None, lpf=None, fbpf=None, flpf=None): u = convolvefilter(u, lpf, flpf) return u + def convolvefilter(x, f, ff=None, isreal=False): """ Apply filter f to a timestream x (convolution) @@ -1607,18 +1785,18 @@ def convolvefilter(x, f, ff=None, isreal=False): # since f might not be numpy.ndarray assert np.ndim(f) == 1 - ## Get shapes + # Get shapes init_shape = x.shape n = f.size nt = init_shape[-1] - ## How many FFTs to perform + # How many FFTs to perform fftsize = int(2 ** np.ceil(np.log2(nt + 3 * n - 1))) fftslice = (slice(None), slice((3 * n - 1) // 2, (3 * n - 1) // 2 + nt)) x2 = x.reshape(-1, nt) - u = np.zeros((init_shape[0], fftsize), (np.ones(1, x.dtype) + 0.j).dtype) - (u[:, : n].T)[:] = x2[:, 0] + u = np.zeros((init_shape[0], fftsize), (np.ones(1, x.dtype) + 0.0j).dtype) + (u[:, :n].T)[:] = x2[:, 0] u[:, n: nt + n] = x2 (u[:, nt + n: nt + n * 2].T)[:] = x2[:, -1] @@ -1630,7 +1808,7 @@ def convolvefilter(x, f, ff=None, isreal=False): u1 *= ff fftpack.ifft(u1, fftsize, overwrite_x=True) - ## Condition to return the real part + # Condition to return the real part cond1 = x.dtype == float or x.dtype == np.float32 cond2 = x.dtype == int or x.dtype == np.int32 cond3 = f.dtype == float or f.dtype == np.float32 @@ -1641,10 +1819,12 @@ def convolvefilter(x, f, ff=None, isreal=False): return u[fftslice].reshape(init_shape) -class WhiteNoiseGenerator(): +class WhiteNoiseGenerator: """ Class to handle white noise """ - def __init__(self, array_noise_level, ndetectors, ntimesamples, - array_noise_seed): + + def __init__( + self, array_noise_level, ndetectors, ntimesamples, array_noise_seed + ): """ This class is used to simulate time-domain noise. Usually, it is used in combination with map2tod to insert noise @@ -1668,9 +1848,10 @@ def __init__(self, array_noise_level, ndetectors, ntimesamples, self.ndetectors = ndetectors self.ntimesamples = ntimesamples - ## Noise level for one detector - self.detector_noise_level = self.array_noise_level * \ - np.sqrt(self.ndetectors) + # Noise level for one detector + self.detector_noise_level = self.array_noise_level * np.sqrt( + self.ndetectors + ) self.array_noise_seed = array_noise_seed state = np.random.RandomState(self.array_noise_seed) @@ -1703,11 +1884,23 @@ def simulate_noise_one_detector(self, ch): return self.detector_noise_level * vec + class CorrNoiseGenerator(WhiteNoiseGenerator): """ """ - def __init__(self, array_noise_level, ndetectors, ntimesamples, - array_noise_seed, nclouds=10, f0=0.1, alpha=-4, amp_atm=1e2, - corrlength=300, sampling_freq=8): + + def __init__( + self, + array_noise_level, + ndetectors, + ntimesamples, + array_noise_seed, + nclouds=10, + f0=0.1, + alpha=-4, + amp_atm=1e2, + corrlength=300, + sampling_freq=8, + ): """ This class is used to simulate time-domain correlated noise. Usually, it is used in combination with map2tod to insert noise @@ -1774,8 +1967,8 @@ def __init__(self, array_noise_level, ndetectors, ntimesamples, """ WhiteNoiseGenerator.__init__( - self, array_noise_level, ndetectors, - ntimesamples, array_noise_seed) + self, array_noise_level, ndetectors, ntimesamples, array_noise_seed + ) self.nclouds = nclouds self.alpha = alpha self.sampling_freq = sampling_freq @@ -1783,7 +1976,7 @@ def __init__(self, array_noise_level, ndetectors, ntimesamples, self.f0 = f0 self.amp_atm = amp_atm - ## Bolometers in a pair get the same seed for correlated noise + # Bolometers in a pair get the same seed for correlated noise self.pixel_noise_seeds = np.repeat(self.noise_seeds[::2], 2) def simulate_noise_one_detector(self, ch): @@ -1811,44 +2004,41 @@ def simulate_noise_one_detector(self, ch): [ -7536.5882971 -224.58319073 -10795.19644268 ..., -5528.66256308 -3161.93996673 -5174.84161989] """ - ## White noise part + # White noise part state = np.random.RandomState(self.noise_seeds[ch]) vec = state.normal(size=self.ntimesamples) wnoise = self.detector_noise_level * vec - ## Correlated part + # Correlated part state = np.random.RandomState(self.pixel_noise_seeds[ch]) # amps = 2 * (-0.5 + state.uniform(size=1)) amps = state.uniform(size=1) corrdet = int(self.ndetectors / self.nclouds) - state = np.random.RandomState( - self.array_noise_seed + ch // corrdet) + state = np.random.RandomState(self.array_noise_seed + ch // corrdet) phases = 2 * np.pi * state.rand(self.ntimesamples) ts_corr = np.zeros(self.ntimesamples) for i in range(0, self.ntimesamples, self.corrlength): - ## Check that you have enough samples + # Check that you have enough samples if self.ntimesamples - i < self.corrlength: step = self.ntimesamples - i else: step = self.corrlength - ## Get the PSD and the frequency range - fs = fftfreq(step, 1. / self.sampling_freq) + # Get the PSD and the frequency range + fs = fftfreq(step, 1.0 / self.sampling_freq) psd = np.zeros_like(fs) - ## Avoid zero frequency - psd[1:] = self.amp_atm * (1 + (fs[1:]/self.f0)**self.alpha) + # Avoid zero frequency + psd[1:] = self.amp_atm * (1 + (fs[1:] / self.f0) ** self.alpha) - ## Get the TOD from the PSD - ts_corr[i: i+step] = corr_ts( - PSD=psd, - N=step, - amp=amps, - phase=phases[i: i+step]) + # Get the TOD from the PSD + ts_corr[i: i + step] = corr_ts( + PSD=psd, N=step, amp=amps, phase=phases[i: i + step] + ) - ## remove PSD normalisation and add white noise! + # remove PSD normalisation and add white noise! return ts_corr / np.sqrt(self.sampling_freq) + wnoise @@ -1888,18 +2078,19 @@ def corr_ts(PSD, N, amp, phase): -0.49232085 0.07355544 0.68893301 0.46171998 -0.82563144] """ Nf = len(PSD) - PSD[0] = 0. + PSD[0] = 0.0 A = amp * N * np.sqrt(PSD) - FFT = A * np.exp(1j*phase[: Nf]) + FFT = A * np.exp(1j * phase[:Nf]) ts = np.fft.ifft(FFT, n=N) return np.real(ts) + def psdts(ts, sample_rate, NFFT=4096): - ''' + """ Compute the power spectrum (or power spectral density) of a timestream ts. Parameters @@ -1925,15 +2116,16 @@ def psdts(ts, sample_rate, NFFT=4096): >>> fs, psd = psdts(ts, sample_rate=1) >>> print(round(fs[-1], 2), round(psd[-1], 2)) 0.49 0.5 - ''' - ## Remove the mean + """ + # Remove the mean ts -= np.mean(ts) fs, asd = compute_asd(ts, sample_rate=sample_rate, NFFT=NFFT) - psd = asd**2 + psd = asd ** 2 return fs, psd + def compute_asd(x, sample_rate, NFFT=2048, is_complex=False): """ Compute the amplitude spectral density (PSD = ASD**2). @@ -1957,25 +2149,25 @@ def compute_asd(x, sample_rate, NFFT=2048, is_complex=False): asd : 1d array Amplitude spectral density of x. """ - ## Sanity check in the case of no chunking + # Sanity check in the case of no chunking maxnfft = x.shape[0] - (x.shape[0] % 2) - if(NFFT > maxnfft or NFFT < 0): + if NFFT > maxnfft or NFFT < 0: NFFT = maxnfft period = 1.0 / sample_rate window = np.blackman(NFFT) - window_norm = 1.0 / np.average(window**2) + window_norm = 1.0 / np.average(window ** 2) n = len(x) - ## Truncates and rounds down + # Truncates and rounds down nchunks = int((2 * n) / NFFT) if is_complex: PSD = np.zeros(NFFT) else: - PSD = np.zeros(int(NFFT/2)) + PSD = np.zeros(int(NFFT / 2)) - for i in range(0, nchunks-1): + for i in range(0, nchunks - 1): chunk = window * x[i * int(NFFT / 2): (i + 2) * int(NFFT / 2)] fch = fft(chunk) @@ -1991,10 +2183,10 @@ def compute_asd(x, sample_rate, NFFT=2048, is_complex=False): fs = fftfreq(NFFT, period) if not is_complex: - fs = fs[0: int(NFFT/2)] + fs = fs[0: int(NFFT / 2)] if nchunks != 1: - PSD /= (nchunks-1) + PSD /= nchunks - 1 # Numpy normalizes inverse transform PSD /= NFFT @@ -2010,13 +2202,21 @@ def compute_asd(x, sample_rate, NFFT=2048, is_complex=False): PSD = fftshift(PSD) # Convert to amplitude/rtHz - return fs, PSD**0.5 + return fs, PSD ** 0.5 + -class OutputSkyMap(): +class OutputSkyMap: """ Class to handle sky maps generated by tod2map """ - def __init__(self, projection, - obspix=None, npixsky=None, - nside=None, pixel_size=None, demodulation=False): + + def __init__( + self, + projection, + obspix=None, + npixsky=None, + nside=None, + pixel_size=None, + demodulation=False, + ): """ Initialise all maps: weights, projected TOD, and Stokes parameters. @@ -2048,30 +2248,35 @@ def __init__(self, projection, self.pixel_size = pixel_size self.demodulation = demodulation - if self.projection == 'healpix': - assert self.obspix is not None, \ - ValueError("You need to provide the list (obspix) " + - "of observed pixel if projection=healpix!") - assert self.nside is not None, \ - ValueError("You need to provide the resolution (nside) " + - "of the map if projection=healpix!") + if self.projection == "healpix": + assert self.obspix is not None, ValueError( + """You need to provide the list (obspix) + of observed pixel if projection=healpix!""" + ) + assert self.nside is not None, ValueError( + """You need to provide the resolution (nside) + of the map if projection=healpix!""" + ) self.npixsky = len(self.obspix) self.pixel_size = hp.nside2resol(nside, arcmin=True) - elif self.projection == 'flat': - assert self.npixsky is not None, \ - ValueError("You need to provide the number of " + - "observed pixels (npixsky) if projection=flat.") - assert self.pixel_size is not None, \ - ValueError("You need to provide the size of " + - "pixels (pixel_size) in arcmin if projection=flat.") + elif self.projection == "flat": + assert self.npixsky is not None, ValueError( + """You need to provide the number of + observed pixels (npixsky) if projection=flat.""" + ) + assert self.pixel_size is not None, ValueError( + """You need to provide the size of + pixels (pixel_size) in arcmin if projection=flat. + """ + ) self.initialise_sky_maps() if not self.demodulation: - self.to_coadd = 'd dc ds w cc cs ss nhit' + self.to_coadd = "d dc ds w cc cs ss nhit" else: - self.to_coadd = 'd0 d4r d4i w0 w4 nhit' + self.to_coadd = "d0 d4r d4i w0 w4 nhit" def initialise_sky_maps(self): """ @@ -2120,15 +2325,15 @@ def set_idet(self): """ testcc = self.cc * self.ss - self.cs * self.cs idet = np.zeros(testcc.shape) - inonzero = (testcc != 0.) - idet[inonzero] = 1. / testcc[inonzero] + inonzero = testcc != 0.0 + idet[inonzero] = 1.0 / testcc[inonzero] thresh = np.finfo(np.float32).eps try: - izero = (np.abs(testcc) < thresh) + izero = np.abs(testcc) < thresh except FloatingPointError: izero = inan = np.isnan(testcc) - izero[~inan] = (np.abs(testcc[~inan]) < thresh) + izero[~inan] = np.abs(testcc[~inan]) < thresh idet[izero] = 0.0 self.idet = idet @@ -2157,10 +2362,10 @@ def get_I(self): hit = self.w > 0 I = np.zeros_like(self.d) - I[hit] = self.d[hit]/self.w[hit] + I[hit] = self.d[hit] / self.w[hit] return I - def get_QU(self,force=True): + def get_QU(self, force=True): """ Solve for the polarisation maps from projected difference timestream maps and weights: @@ -2190,7 +2395,7 @@ def get_QU(self,force=True): if self.demodulation: return self.get_QU_demod() - if ((not hasattr(self, 'idet')) or force): + if (not hasattr(self, "idet")) or force: self.set_idet() Q = self.idet * (self.ss * self.dc - self.cs * self.ds) @@ -2293,13 +2498,14 @@ def coadd(self, other, to_coadd=None): >>> print(m1.nhit) [ 2. 2. 2. 2.] """ - assert np.all(self.obspix == other.obspix), \ - ValueError("To add maps together, they must have the same obspix!") + assert np.all(self.obspix == other.obspix), ValueError( + "To add maps together, they must have the same obspix!" + ) if to_coadd is None: to_coadd = self.to_coadd - to_coadd_split = to_coadd.split(' ') + to_coadd_split = to_coadd.split(" ") for k in to_coadd_split: a = getattr(self, k) b = getattr(other, k) @@ -2326,19 +2532,23 @@ def coadd_MPI(self, other, MPI, to_coadd=None): >>> from mpi4py import MPI >>> m = OutputSkyMap(projection='healpix', ... nside=16, obspix=np.array([0, 1, 2, 3])) - >>> ## do whatever you want with the maps + >>> # do whatever you want with the maps >>> m.coadd_MPI(m, MPI) """ if to_coadd is None: to_coadd = self.to_coadd - to_coadd_split = to_coadd.split(' ') + to_coadd_split = to_coadd.split(" ") for k in to_coadd_split: - setattr(self, k, MPI.COMM_WORLD.allreduce( - getattr(other, k), op=MPI.SUM)) + setattr( + self, + k, + MPI.COMM_WORLD.allreduce(getattr(other, k), op=MPI.SUM), + ) - def pickle_me(self, fn, shrink_maps=True, crop_maps=False, - epsilon=0., verbose=False): + def pickle_me( + self, fn, shrink_maps=True, crop_maps=False, epsilon=0.0, verbose=False + ): """ Save data into pickle file. Work with pair differenced data only. See pickle_me_demod for @@ -2359,37 +2569,58 @@ def pickle_me(self, fn, shrink_maps=True, crop_maps=False, """ if self.demodulation: self.pickle_me_demod( - fn, shrink_maps=shrink_maps, crop_maps=crop_maps, - epsilon=epsilon, verbose=verbose) + fn, + shrink_maps=shrink_maps, + crop_maps=crop_maps, + epsilon=epsilon, + verbose=verbose, + ) else: try: I, Q, U = self.get_IQU() - wP = qu_weight_mineig(self.cc, self.cs, self.ss, - epsilon=epsilon, verbose=verbose) - - data = {'I': I, 'Q': Q, 'U': U, - 'wI': self.w, 'wP': wP, 'nhit': self.nhit, - 'projection': self.projection, - 'nside': self.nside, 'pixel_size': self.pixel_size, - 'obspix': self.obspix} + wP = qu_weight_mineig( + self.cc, self.cs, self.ss, epsilon=epsilon, verbose=verbose + ) + + data = { + "I": I, + "Q": Q, + "U": U, + "wI": self.w, + "wP": wP, + "nhit": self.nhit, + "projection": self.projection, + "nside": self.nside, + "pixel_size": self.pixel_size, + "obspix": self.obspix, + } except Exception as e: - print('Exception error: ', e) - I, G, Q, U = self.get_IQU() # if using IGQU class - wP = qu_weight_mineig(self.cc, self.cs, self.ss, - epsilon=epsilon, verbose=verbose) - - data = {'I': I, 'G': G, 'Q': Q, 'U': U, - 'wI': self.w, 'wP': wP, 'nhit': self.nhit, - 'projection': self.projection, - 'nside': self.nside, 'pixel_size': self.pixel_size, - 'obspix': self.obspix} - - if shrink_maps and self.projection == 'flat': - data = shrink_me(data, based_on='wP') - elif crop_maps is not False and self.projection == 'flat': - data = crop_me(data, based_on='wP', npix_per_row=crop_maps) - - with open(fn, 'wb') as f: + print("Exception error: ", e) + I, G, Q, U = self.get_IQU() # if using IGQU class + wP = qu_weight_mineig( + self.cc, self.cs, self.ss, epsilon=epsilon, verbose=verbose + ) + + data = { + "I": I, + "G": G, + "Q": Q, + "U": U, + "wI": self.w, + "wP": wP, + "nhit": self.nhit, + "projection": self.projection, + "nside": self.nside, + "pixel_size": self.pixel_size, + "obspix": self.obspix, + } + + if shrink_maps and self.projection == "flat": + data = shrink_me(data, based_on="wP") + elif crop_maps is not False and self.projection == "flat": + data = crop_me(data, based_on="wP", npix_per_row=crop_maps) + + with open(fn, "wb") as f: pickle.dump(data, f, protocol=2) def initialise_sky_maps_demod(self): @@ -2441,7 +2672,7 @@ def get_I_demod(self): """ hit = self.w0 > 0 I = np.zeros_like(self.d0) - I[hit] = self.d0[hit]/self.w0[hit] + I[hit] = self.d0[hit] / self.w0[hit] return I def get_QU_demod(self): @@ -2475,13 +2706,14 @@ def get_QU_demod(self): Q = np.zeros_like(self.d4r) U = np.zeros_like(self.d4i) - Q[hit] = self.d4r[hit]/self.w4[hit] - U[hit] = self.d4i[hit]/self.w4[hit] + Q[hit] = self.d4r[hit] / self.w4[hit] + U[hit] = self.d4i[hit] / self.w4[hit] return Q, U - def pickle_me_demod(self, fn, shrink_maps=True, crop_maps=False, - epsilon=0., verbose=False): + def pickle_me_demod( + self, fn, shrink_maps=True, crop_maps=False, epsilon=0.0, verbose=False + ): """ Save data into pickle file. Work with demodulated data only. @@ -2501,26 +2733,39 @@ def pickle_me_demod(self, fn, shrink_maps=True, crop_maps=False, """ I, Q, U = self.get_IQU() - data = {'I': I, 'Q': Q, 'U': U, - 'wI': self.w0, 'wP': self.w4, 'nhit': self.nhit, - 'projection': self.projection, - 'nside': self.nside, 'pixel_size': self.pixel_size, - 'obspix': self.obspix} - - if shrink_maps and self.projection == 'flat': - data = shrink_me(data, based_on='wP') - elif crop_maps is not False and self.projection == 'flat': - data = crop_me(data, based_on='wP', npix_per_row=crop_maps) - - with open(fn, 'wb') as f: + data = { + "I": I, + "Q": Q, + "U": U, + "wI": self.w0, + "wP": self.w4, + "nhit": self.nhit, + "projection": self.projection, + "nside": self.nside, + "pixel_size": self.pixel_size, + "obspix": self.obspix, + } + + if shrink_maps and self.projection == "flat": + data = shrink_me(data, based_on="wP") + elif crop_maps is not False and self.projection == "flat": + data = crop_me(data, based_on="wP", npix_per_row=crop_maps) + + with open(fn, "wb") as f: pickle.dump(data, f, protocol=2) class OutputSkyMapIGQU(OutputSkyMap): """ Class to handle sky maps generated by tod2map + G deprojection """ - def __init__(self, projection, - obspix=None, npixsky=None, - nside=None, pixel_size=None): + + def __init__( + self, + projection, + obspix=None, + npixsky=None, + nside=None, + pixel_size=None, + ): """ Initialise all maps: weights, projected TOD, and Stokes parameters. Suitable if you want to perform a deprojection of a constant @@ -2548,22 +2793,25 @@ def __init__(self, projection, OutputSkyMap.__init__( self, projection=projection, - obspix=obspix, npixsky=npixsky, - nside=nside, pixel_size=pixel_size, - demodulation=False) + obspix=obspix, + npixsky=npixsky, + nside=nside, + pixel_size=pixel_size, + demodulation=False, + ) - self.to_coadd = self.to_coadd + ' dm wm c s' + self.to_coadd = self.to_coadd + " dm wm c s" - ## Map and weight of deprojected spurious signal + # Map and weight of deprojected spurious signal self.dm = np.zeros(self.npixsky) self.wm = np.zeros(self.npixsky) - ## cosine and sine (weights for cross term GxQ and GxU) + # cosine and sine (weights for cross term GxQ and GxU) self.c = np.zeros(self.npixsky) self.s = np.zeros(self.npixsky) def buildV(self, ipix): - """ Build vector of projected polarisation TOD for a given sky pixel + """Build vector of projected polarisation TOD for a given sky pixel It includes G template for deprojection (dm). @@ -2586,7 +2834,7 @@ def buildV(self, ipix): return np.array([self.dm[ipix], self.dc[ipix], self.ds[ipix]]) def buildP(self, ipix): - """ Build pixel weight matrix for polarisation. + """Build pixel weight matrix for polarisation. The [3x3] matrix contains the weights for G, Q and U components of the mapmaking: @@ -2634,7 +2882,7 @@ def set_goodpix(self): ... nside=16, obspix=np.array(range(12*16**2))) >>> m1.set_goodpix() """ - ## We have 4 components: I, Q, U and G. + # We have 4 components: I, Q, U and G. inonzero = [pix for pix in range(self.npixsky) if self.nhit[pix] > 10] self.goodpix = np.zeros((self.npixsky)) self.goodpix[inonzero] = 1 @@ -2676,10 +2924,11 @@ def get_QU(self): ... nside=16, obspix=np.array([0, 1, 2, 3])) >>> G, Q, U = m1.get_QU() """ - if not hasattr(self, 'goodpix'): + if not hasattr(self, "goodpix"): self.set_goodpix() inonzero = [ - pix for pix in range(self.npixsky) if self.goodpix[pix] != 0] + pix for pix in range(self.npixsky) if self.goodpix[pix] != 0 + ] G = np.zeros((self.npixsky)) Q = np.zeros((self.npixsky)) @@ -2772,10 +3021,11 @@ def shrink_me(dic, based_on): [[ 6 7] [10 11]] """ - assert based_on in dic, \ - KeyError("{} not in input dictionary!".format(based_on)) + assert based_on in dic, KeyError( + "{} not in input dictionary!".format(based_on) + ) - npixr = int(len(dic[based_on])**.5) + npixr = int(len(dic[based_on]) ** 0.5) halfnpixr = int(npixr / 2) mask = dic[based_on].reshape((npixr, npixr)) idx = np.where(mask > 0) @@ -2788,21 +3038,24 @@ def shrink_me(dic, based_on): halfdxy = int(dxy / 2) for k in dic.keys(): - ## Filter out fields which aren't arrays + # Filter out fields which aren't arrays if type(dic[k]) == np.ndarray: - ## Filter out fields which are arrays but not like based_on. - npixr_loc = int(len(dic[k])**.5) + # Filter out fields which are arrays but not like based_on. + npixr_loc = int(len(dic[k]) ** 0.5) if npixr_loc == npixr: dic[k] = np.array( - [i[halfnpixr - halfdxy - 1: halfnpixr + halfdxy+1] for i in - dic[k].reshape( - (npixr, npixr))[ - halfnpixr - halfdxy - 1: - halfnpixr + halfdxy + 1]]).flatten() + [ + i[halfnpixr - halfdxy - 1: halfnpixr + halfdxy + 1] + for i in dic[k].reshape((npixr, npixr))[ + halfnpixr - halfdxy - 1: halfnpixr + halfdxy + 1 + ] + ] + ).flatten() return dic -def crop_me(dic, based_on, npix_per_row=2**12): + +def crop_me(dic, based_on, npix_per_row=2 ** 12): """ Crop maps to a chosen size. Maps have to be squared (so work only for flat sky). @@ -2838,32 +3091,38 @@ def crop_me(dic, based_on, npix_per_row=2**12): [[ 6 7] [10 11]] """ - assert based_on in dic, \ - KeyError("{} not in input dictionary!".format(based_on)) + assert based_on in dic, KeyError( + "{} not in input dictionary!".format(based_on) + ) - npixr = int(len(dic[based_on])**.5) + npixr = int(len(dic[based_on]) ** 0.5) halfnpixr = int(npixr / 2) halfnpix_per_row = int(npix_per_row / 2) for k in dic.keys(): - ## Filter out fields which aren't arrays + # Filter out fields which aren't arrays if type(dic[k]) == np.ndarray: - ## Filter out fields which are arrays but not like based_on. - npixr_loc = int(len(dic[k])**.5) + # Filter out fields which are arrays but not like based_on. + npixr_loc = int(len(dic[k]) ** 0.5) if npixr_loc == npixr: - ## Check that we have enough pixels to start with - assert npixr_loc >= npix_per_row, \ - ValueError("Map too small to " + - "be cropped! ({} vs {})".format( - npixr_loc, npix_per_row)) + # Check that we have enough pixels to start with + assert npixr_loc >= npix_per_row, ValueError( + "Map too small to be cropped! ({} vs {})".format( + npixr_loc, + npix_per_row + ) + ) dic[k] = np.array( - [i[halfnpixr - halfnpix_per_row: - halfnpixr + halfnpix_per_row] for i in - dic[k].reshape( - (npixr, npixr))[halfnpixr - halfnpix_per_row: - halfnpixr + - halfnpix_per_row]]).flatten() + [ + i[ + halfnpixr - halfnpix_per_row: halfnpixr + halfnpix_per_row + ] + for i in dic[k].reshape((npixr, npixr))[ + halfnpixr - halfnpix_per_row: halfnpixr + halfnpix_per_row + ] + ] + ).flatten() return dic @@ -2897,15 +3156,25 @@ def partial2full(partial_obs, obspix, nside, fill_with=0.0): >>> obspix = np.arange(12 * nside**2, dtype=int)[30:40] >>> fullsky = partial2full(data, obspix, nside) """ - fullsky = np.zeros(12 * nside**2) * fill_with + fullsky = np.zeros(12 * nside ** 2) * fill_with fullsky[obspix] = partial_obs return fullsky -def build_pointing_matrix(ra, dec, nside_in, nside_out=None, - projection='healpix', obspix=None, ext_map_gal=False, - xmin=None, ymin=None, - pixel_size=None, npix_per_row=None, - cut_pixels_outside=True): + +def build_pointing_matrix( + ra, + dec, + nside_in, + nside_out=None, + projection="healpix", + obspix=None, + ext_map_gal=False, + xmin=None, + ymin=None, + pixel_size=None, + npix_per_row=None, + cut_pixels_outside=True, +): """ Given pointing coordinates (RA/Dec), retrieve the corresponding healpix pixel index for a full sky map. This acts effectively as an operator @@ -2991,12 +3260,12 @@ def build_pointing_matrix(ra, dec, nside_in, nside_out=None, theta, phi = radec2thetaphi(ra, dec) if ext_map_gal: - r = hp.Rotator(coord=['C', 'G']) + r = hp.Rotator(coord=["C", "G"]) theta, phi = r(theta, phi) index_global = hp.ang2pix(nside_in, theta, phi) - if projection == 'healpix' and obspix is not None: + if projection == "healpix" and obspix is not None: index_global_out = hp.ang2pix(nside_out, theta, phi) index_local = obspix.searchsorted(index_global_out) mask1 = index_local < len(obspix) @@ -3004,25 +3273,29 @@ def build_pointing_matrix(ra, dec, nside_in, nside_out=None, loc[mask1] = obspix[index_local[mask1]] == index_global_out[mask1] outside_pixels = np.invert(loc) - ## Handling annoying cases. - if (np.sum(outside_pixels) and (not cut_pixels_outside)): - msg = "Pixels outside patch boundaries. " + \ - "Patch width insufficient. To avoid this, " + \ - "increase the parameter width while initialising the TOD " + \ - "or set cut_pixels_outside to True to get a cropped map." + # Handling annoying cases. + if np.sum(outside_pixels) and (not cut_pixels_outside): + msg = """ + Pixels outside patch boundaries. + Patch width insufficient. To avoid this, + increase the parameter width while initialising the TOD + or set cut_pixels_outside to True to get a cropped map. + """ raise ValueError(msg) - elif (np.sum(outside_pixels) and cut_pixels_outside): - if (not ('msg_cut' in globals())): + elif np.sum(outside_pixels) and cut_pixels_outside: + if not ("msg_cut" in globals()): global msg_cut - msg_cut = "Pixels outside patch boundaries. " + \ - "Your output map will be cropped. To avoid this, " + \ - "increase the parameter width while initialising the TOD." + msg_cut = """ + Pixels outside patch boundaries. + Your output map will be cropped. To avoid this, + increase the parameter width while initialising the TOD. + """ print(msg_cut) index_local[outside_pixels] = -1 else: pass - elif projection == 'flat': + elif projection == "flat": x, y = input_sky.LamCyl(ra, dec) xminmap = xmin - pixel_size / 2.0 @@ -3033,22 +3306,27 @@ def build_pointing_matrix(ra, dec, nside_in, nside_out=None, index_local = ix * npix_per_row + iy - outside_pixels = (ix < 0) | (ix >= npix_per_row) | \ - (iy < 0) | (iy >= npix_per_row) - - ## Handling annoying cases. - if (np.sum(outside_pixels) and (not cut_pixels_outside)): - msg = "Pixels outside patch boundaries. " + \ - "Patch width insufficient. To avoid this, " + \ - "increase the parameter width while initialising the TOD " + \ - "or set cut_pixels_outside to True to get a cropped map." + outside_pixels = ( + (ix < 0) | (ix >= npix_per_row) | (iy < 0) | (iy >= npix_per_row) + ) + + # Handling annoying cases. + if np.sum(outside_pixels) and (not cut_pixels_outside): + msg = """ + Pixels outside patch boundaries. + Patch width insufficient. To avoid this, + increase the parameter width while initialising the TOD + or set cut_pixels_outside to True to get a cropped map. + """ raise ValueError(msg) - elif (np.sum(outside_pixels) and cut_pixels_outside): - if (not ('msg_cut_flat' in globals())): + elif np.sum(outside_pixels) and cut_pixels_outside: + if not ("msg_cut_flat" in globals()): global msg_cut_flat - msg_cut_flat = "Pixels outside patch boundaries. " + \ - "Your output map will be cropped. To avoid this, " + \ - "increase the parameter width while initialising the TOD." + msg_cut_flat = """ + Pixels outside patch boundaries. + Your output map will be cropped. To avoid this, + increase the parameter width while initialising the TOD. + """ print(msg_cut_flat) index_local[outside_pixels] = -1 @@ -3059,8 +3337,10 @@ def build_pointing_matrix(ra, dec, nside_in, nside_out=None, return index_global, index_local -def load_fake_instrument(nside=16, nsquid_per_mux=1, fwhm_in2=None, - compute_derivatives=False): + +def load_fake_instrument( + nside=16, nsquid_per_mux=1, fwhm_in2=None, compute_derivatives=False +): """ For test purposes. Create instances of HealpixFitsMap, hardware, and @@ -3075,44 +3355,65 @@ def load_fake_instrument(nside=16, nsquid_per_mux=1, fwhm_in2=None, HealpixFitsMap : HealpixFitsMap instance Instance of HealpixFitsMap containing input sky parameters. """ - ## Add paths to load modules - sys.path.insert(0, os.path.realpath(os.path.join(os.getcwd(), '.'))) - sys.path.insert(0, os.path.realpath(os.path.join(os.getcwd(), 's4cmb'))) + # Add paths to load modules + sys.path.insert(0, os.path.realpath(os.path.join(os.getcwd(), "."))) + sys.path.insert(0, os.path.realpath(os.path.join(os.getcwd(), "s4cmb"))) from s4cmb.input_sky import HealpixFitsMap from s4cmb.instrument import Hardware from s4cmb.scanning_strategy import ScanningStrategy - ## Create fake inputs - - ## Sky - sky_in = HealpixFitsMap('s4cmb/data/test_data_set_lensedCls.dat', - do_pol=True, fwhm_in=0.0, fwhm_in2=fwhm_in2, - nside_in=nside, map_seed=48584937, - compute_derivatives=compute_derivatives, - verbose=False, no_ileak=False, no_quleak=False) - - ## Instrument - inst = Hardware(ncrate=1, ndfmux_per_crate=1, - nsquid_per_mux=nsquid_per_mux, npair_per_squid=4, - fp_size=60., fwhm=3.5, - beam_seed=58347, projected_fp_size=3., - pm_name='5params', type_hwp='CRHWP', - freq_hwp=0.2, angle_hwp=0., verbose=False) + # Create fake inputs + + # Sky + sky_in = HealpixFitsMap( + "s4cmb/data/test_data_set_lensedCls.dat", + do_pol=True, + fwhm_in=0.0, + fwhm_in2=fwhm_in2, + nside_in=nside, + map_seed=48584937, + compute_derivatives=compute_derivatives, + verbose=False, + no_ileak=False, + no_quleak=False, + ) + + # Instrument + inst = Hardware( + ncrate=1, + ndfmux_per_crate=1, + nsquid_per_mux=nsquid_per_mux, + npair_per_squid=4, + fp_size=60.0, + fwhm=3.5, + beam_seed=58347, + projected_fp_size=3.0, + pm_name="5params", + type_hwp="CRHWP", + freq_hwp=0.2, + angle_hwp=0.0, + verbose=False, + ) if fwhm_in2 is not None: inst.make_dichroic(fwhm=fwhm_in2) - ## Scanning strategy - scan = ScanningStrategy(nces=2, start_date='2013/1/1 00:00:00', - telescope_longitude='-67:46.816', - telescope_latitude='-22:56.396', - telescope_elevation=5200., - name_strategy='deep_patch', - sampling_freq=8., sky_speed=0.4, - language='fortran') + # Scanning strategy + scan = ScanningStrategy( + nces=2, + start_date="2013/1/1 00:00:00", + telescope_longitude="-67:46.816", + telescope_latitude="-22:56.396", + telescope_elevation=5200.0, + name_strategy="deep_patch", + sampling_freq=8.0, + sky_speed=0.4, + language="fortran", + ) scan.run() return inst, scan, sky_in + def noise_ukam(array_noise_level, fsky, nside, tobs): """ Estimate quickly the noise level in map domain [uK.arcmin] @@ -3146,12 +3447,15 @@ def noise_ukam(array_noise_level, fsky, nside, tobs): >>> print(round(noise, 2), 'uK.arcmin') 3.93 uK.arcmin """ - noise = np.sqrt(array_noise_level**2 * hp.nside2npix(nside) * fsky / tobs) + noise = np.sqrt( + array_noise_level ** 2 * hp.nside2npix(nside) * fsky / tobs + ) return noise * hp.nside2resol(nside, arcmin=True) if __name__ == "__main__": import doctest + if np.__version__ >= "1.14.0": np.set_printoptions(legacy="1.13") doctest.testmod() From 67949e77c26a73d430c47cbf4f816178a9f3b266 Mon Sep 17 00:00:00 2001 From: Julien Date: Mon, 25 Jan 2021 15:18:40 +0100 Subject: [PATCH 02/15] Remove unused code --- s4cmb/tod.py | 28 ---------------------------- 1 file changed, 28 deletions(-) diff --git a/s4cmb/tod.py b/s4cmb/tod.py index 74ccbb2..342b411 100755 --- a/s4cmb/tod.py +++ b/s4cmb/tod.py @@ -954,34 +954,6 @@ def map2tod(self, ch): # Compute pointing matrix index_global, index_local = self.get_pixel_indices(ra, dec) - # Perturbed boresight pointing values - if self.pointing_perturbed is not None: - ( - ra_perturbed, - dec_perturbed, - pa_perturbed, - ) = self.pointing_perturbed.offset_detector(azd, eld) - - # Perturbed boresight pointing values - ( - index_global_perturbed, - index_local_perturbed, - ) = self.get_pixel_indices(ra_perturbed, dec_perturbed) - - else: - ra_perturbed = ra - dec_perturbed = dec - pa_perturbed = pa - index_global_perturbed = index_global - index_local_perturbed = index_local - - if (self.projection == "healpix") & (index_local is None): - # Using a pointer not to increase memory usage - index_local = index_global - - # Perturbed boresight pointing values - index_local_perturbed = index_global_perturbed - # For flat projection, one needs to flip the sign of U # w.r.t to the full-sky basis (IAU angle convention) if self.projection == "flat": From 550a446e639a3078d23ea7afeb4aee4e6e381d59 Mon Sep 17 00:00:00 2001 From: Julien Date: Mon, 25 Jan 2021 15:29:36 +0100 Subject: [PATCH 03/15] TODO: fix this wrong exceptopn --- s4cmb/tod.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/s4cmb/tod.py b/s4cmb/tod.py index 342b411..7bb48a7 100755 --- a/s4cmb/tod.py +++ b/s4cmb/tod.py @@ -1037,12 +1037,13 @@ def map2tod(self, ch): # Use pa of pair center for tod2map if differential pointing # is used. Otherwise the polarization angle of a pair is the # same as the top and bottom detectors. + # TODO: I do not understand this try/except. + # TODO: remove the bare exception!!! try: pol_ang_pair, pol_ang2_pair = self.compute_simpolangle( ch, pa_pair, polangle_err=False ) - except Exception as e: - print(str(e)) + except: pol_ang_pair = pol_ang pol_ang2_pair = pol_ang2 From 9fbac149f3bd2ebcf6e0f734cd30279bac591421 Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 10:18:45 +0100 Subject: [PATCH 04/15] Clean systematics --- s4cmb/systematics.py | 487 ++++++++++++++++++++++++++----------------- 1 file changed, 292 insertions(+), 195 deletions(-) diff --git a/s4cmb/systematics.py b/s4cmb/systematics.py index 6d473a1..9d63acd 100755 --- a/s4cmb/systematics.py +++ b/s4cmb/systematics.py @@ -11,16 +11,25 @@ from scipy import signal from s4cmb.systematics_f import systematics_f -from s4cmb.instrument import construct_beammap, coordinates_on_grid,gauss2d - -arcsecond2rad = np.pi / 180. / 3600. -arcmin2rad = np.pi / 180. / 60. -deg2rad = np.pi / 180. - -def inject_crosstalk_inside_SQUID(bolo_data, squid_ids, bolo_ids, - mu=-3., sigma=1., radius=1, beta=2, - seed=5438765, new_array=None, - language='python'): +from s4cmb.instrument import construct_beammap, coordinates_on_grid, gauss2d + +arcsecond2rad = np.pi / 180.0 / 3600.0 +arcmin2rad = np.pi / 180.0 / 60.0 +deg2rad = np.pi / 180.0 + + +def inject_crosstalk_inside_SQUID( + bolo_data, + squid_ids, + bolo_ids, + mu=-3.0, + sigma=1.0, + radius=1, + beta=2, + seed=5438765, + new_array=None, + language="python", +): """ Introduce leakage between neighboring bolometers within a SQUID. You have to provide the list of bolometers, to which SQUID they @@ -96,9 +105,9 @@ def inject_crosstalk_inside_SQUID(bolo_data, squid_ids, bolo_ids, For large number of bolometers per SQUID, you would prefer fortran to python to perform the loops. Choose python otherwise. """ - ## Make mu and sigma unitless (user provides percentage) - mu = mu / 100. - sigma = sigma / 100. + # Make mu and sigma unitless (user provides percentage) + mu = mu / 100.0 + sigma = sigma / 100.0 combs = {} for bolo in range(len(bolo_data)): @@ -112,35 +121,50 @@ def inject_crosstalk_inside_SQUID(bolo_data, squid_ids, bolo_ids, state = np.random.RandomState(seed) cross_amp = state.normal(mu, sigma, len(bolo_data)) - if language == 'python': + if language == "python": for sq in combs: for ch, i in combs[sq]: for ch2, i2 in combs[sq]: separation_length = abs(ch - ch2) if separation_length > 0 and separation_length <= radius: - tsout[i] += cross_amp[i2] / \ - separation_length**beta * tsout[i2] + tsout[i] += ( + cross_amp[i2] / separation_length ** beta * tsout[i2] + ) - elif language == 'fortran': - ## F2PY convention - tsout = np.array(tsout, order='F') + elif language == "fortran": + # F2PY convention + tsout = np.array(tsout, order="F") for sq in combs: - local_indices = np.array(combs[sq]).flatten()[:: 2] - global_indices = np.array(combs[sq]).flatten()[1:: 2] + local_indices = np.array(combs[sq]).flatten()[::2] + global_indices = np.array(combs[sq]).flatten()[1::2] systematics_f.inject_crosstalk_inside_squid_f( - tsout, local_indices, global_indices, - radius, cross_amp, beta, - len(local_indices), len(bolo_data), len(bolo_data[0])) + tsout, + local_indices, + global_indices, + radius, + cross_amp, + beta, + len(local_indices), + len(bolo_data), + len(bolo_data[0]), + ) if new_array is not None: new_array[:] = tsout else: bolo_data[:] = tsout -def inject_crosstalk_SQUID_to_SQUID(bolo_data, squid_ids, bolo_ids, - mu=-3., sigma=1., - squid_attenuation=100., - seed=5438765, new_array=None): + +def inject_crosstalk_SQUID_to_SQUID( + bolo_data, + squid_ids, + bolo_ids, + mu=-3.0, + sigma=1.0, + squid_attenuation=100.0, + seed=5438765, + new_array=None, +): """ Introduce leakage between bolometers from different SQUIDs. You have to provide the list of bolometers, to which SQUID they @@ -207,9 +231,9 @@ def inject_crosstalk_SQUID_to_SQUID(bolo_data, squid_ids, bolo_ids, 40.948 40.716 """ - ## Make mu and sigma unitless (user provides percentage) - mu = mu / 100. - sigma = sigma / 100. + # Make mu and sigma unitless (user provides percentage) + mu = mu / 100.0 + sigma = sigma / 100.0 combs = {} for bolo in range(len(bolo_data)): @@ -223,15 +247,15 @@ def inject_crosstalk_SQUID_to_SQUID(bolo_data, squid_ids, bolo_ids, state = np.random.RandomState(seed) cross_amp = state.normal(mu, sigma, len(bolo_data)) - ## Take one squid + # Take one squid for sq1 in combs: - ## Take all bolometers in that SQUID + # Take all bolometers in that SQUID for ch, i in combs[sq1]: - ## Take a second SQUID + # Take a second SQUID for sq2 in combs: if sq2 == sq1: continue - ## Take all bolometers in that SQUID + # Take all bolometers in that SQUID for ch2, i2 in combs[sq2]: tsout[i] += cross_amp[i2] / squid_attenuation * tsout[i2] @@ -240,8 +264,14 @@ def inject_crosstalk_SQUID_to_SQUID(bolo_data, squid_ids, bolo_ids, else: bolo_data[:] = tsout -def modify_beam_offsets(bolometer_xpos, bolometer_ypos, - mu_diffpointing=10., sigma_diffpointing=5., seed=5847): + +def modify_beam_offsets( + bolometer_xpos, + bolometer_ypos, + mu_diffpointing=10.0, + sigma_diffpointing=5.0, + seed=5847, +): """ Modify the beam offsets (inject differential pointing between two pixel-pair bolometers). The model is the following: @@ -286,22 +316,23 @@ def modify_beam_offsets(bolometer_xpos, bolometer_ypos, [ 0.92369044 -0.92369044 -1.10310547 1.10310547] """ - assert len(bolometer_xpos) == len(bolometer_ypos), \ - ValueError("x and y bolometer coordinates should have the same length") + assert len(bolometer_xpos) == len(bolometer_ypos), ValueError( + "x and y bolometer coordinates should have the same length" + ) state = np.random.RandomState(seed) npair = int(len(bolometer_xpos) / 2) - ## rho is defined by the user [xpos, ypos in radians - convert it!] - rho = state.normal(mu_diffpointing * arcsecond2rad, - sigma_diffpointing * arcsecond2rad, - npair) + # rho is defined by the user [xpos, ypos in radians - convert it!] + rho = state.normal( + mu_diffpointing * arcsecond2rad, sigma_diffpointing * arcsecond2rad, npair + ) - ## Angle are uniformly drawn from 0 and 360 degree + # Angle are uniformly drawn from 0 and 360 degree theta = state.uniform(0, 360, npair) * deg2rad - rhox = rho / 2. * np.cos(theta) - rhoy = rho / 2. * np.sin(theta) + rhox = rho / 2.0 * np.cos(theta) + rhoy = rho / 2.0 * np.sin(theta) bolometer_xpos[::2] += rhox bolometer_xpos[1::2] += -rhox @@ -310,9 +341,15 @@ def modify_beam_offsets(bolometer_xpos, bolometer_ypos, return bolometer_xpos, bolometer_ypos -def inject_beam_ellipticity(sigma_gaussian, mu_beamellipticity, - sigma_beamellipticity, nbolo, - do_diffbeamellipticity=True, seed=54875): + +def inject_beam_ellipticity( + sigma_gaussian, + mu_beamellipticity, + sigma_beamellipticity, + nbolo, + do_diffbeamellipticity=True, + seed=54875, +): """ Inject differential beam ellipticity (and beam ellipticity) Starting from the definition of the beam ellipticity: @@ -382,29 +419,31 @@ def inject_beam_ellipticity(sigma_gaussian, mu_beamellipticity, state = np.random.RandomState(seed) eps = state.normal( - mu_beamellipticity/100., sigma_beamellipticity/100., nbolo) + mu_beamellipticity / 100.0, sigma_beamellipticity / 100.0, nbolo + ) if do_diffbeamellipticity: - d_plus = 2 * sigma_gaussian / eps * (1. + np.sqrt((1 - eps**2))) - d_minus = 2 * sigma_gaussian / eps * (1. - np.sqrt((1 - eps**2))) + # d_plus = 2 * sigma_gaussian / eps * (1.0 + np.sqrt((1 - eps ** 2))) + d_minus = 2 * sigma_gaussian / eps * (1.0 - np.sqrt((1 - eps ** 2))) d = d_minus else: - ## Bolometers in the same pair will have the same ellipticity. - ## This assumes that bolometers i and i+1 belong - ## to the same pair (default behaviour) + # Bolometers in the same pair will have the same ellipticity. + # This assumes that bolometers i and i+1 belong + # to the same pair (default behaviour) eps = np.repeat(eps[::2], 2) - d_plus = 2 * sigma_gaussian / eps * (1. + np.sqrt((1 - eps**2))) - d_minus = 2 * sigma_gaussian / eps * (1. - np.sqrt((1 - eps**2))) + # d_plus = 2 * sigma_gaussian / eps * (1.0 + np.sqrt((1 - eps ** 2))) + d_minus = 2 * sigma_gaussian / eps * (1.0 - np.sqrt((1 - eps ** 2))) d = d_minus - sig_1 = np.ones(nbolo) * sigma_gaussian + d/2. - sig_2 = np.ones(nbolo) * sigma_gaussian - d/2. + sig_1 = np.ones(nbolo) * sigma_gaussian + d / 2.0 + sig_2 = np.ones(nbolo) * sigma_gaussian - d / 2.0 - ## Angle between ellipses + # Angle between ellipses ellip_ang = state.uniform(-90, 90, nbolo) return sig_1, sig_2, ellip_ang + def modify_pointing_parameters(values, errors): """ This routine is not realistic at all for the moment! @@ -425,8 +464,8 @@ def modify_pointing_parameters(values, errors): values_mod = [p + err for p, err in zip(values, errors)] return values_mod -def step_function(nbolos, nsamples, mean=1, std=0.05, - nbreaks=1, sign='same', seed=0): + +def step_function(nbolos, nsamples, mean=1, std=0.05, nbreaks=1, sign="same", seed=0): """ Generate step functions for each bolometer from 1 to a values drawn from N(mean, std). The full timestream is broken into nbreaks @@ -468,37 +507,37 @@ def step_function(nbolos, nsamples, mean=1, std=0.05, >>> gains = step_function(nbolos, nsamples, nbreaks=1, sign='opposite') >>> assert gains[0][0] == 2 - gains[1][0] """ - ## Fix the seed + # Fix the seed state = np.random.RandomState(seed) - ## Initialise gains to ones + # Initialise gains to ones gains = np.ones((nbolos, nsamples)) - ## Length of each break + # Length of each break length = nsamples // nbreaks - sublength = nsamples // (2*nbreaks) + sublength = nsamples // (2 * nbreaks) for pos in range(nbreaks): - ## 2 bolo in a pair will have the same break to avoid differential gain - end_points = state.normal(mean, std, size=int(nbolos/2)) - end_points = np.tile( - end_points, (2, 1)).T.reshape((2*len(end_points), 1)) + # 2 bolo in a pair will have the same break to avoid differential gain + end_points = state.normal(mean, std, size=int(nbolos / 2)) + end_points = np.tile(end_points, (2, 1)).T.reshape((2 * len(end_points), 1)) - ## Assign values + # Assign values shift = pos * length if pos * length < nsamples: - gains[:, shift + sublength:shift + 2 * sublength] = \ - end_points.reshape((len(end_points), 1)) + gains[:, shift + sublength: shift + 2 * sublength] = end_points.reshape( + (len(end_points), 1) + ) else: continue - if sign == 'opposite': + if sign == "opposite": gains[1] = 2 - gains[1] return gains -def step_function_gen(nsamples, mean=1, std=0.05, - nbreaks=1, sign='same', seed=0): + +def step_function_gen(nsamples, mean=1, std=0.05, nbreaks=1, sign="same", seed=0): """ Generator of step functions for each bolometer from 1 to a values drawn from N(mean, std). The full timestream is broken into nbreaks @@ -540,38 +579,41 @@ def step_function_gen(nsamples, mean=1, std=0.05, >>> g = next(gains_gen) >>> assert g[0][0] == 2 - g[1][0] """ - ## Fix the seed + # Fix the seed state = np.random.RandomState(seed) - ## Initialise gains to ones + # Initialise gains to ones gains = np.ones((2, nsamples)) - ## Length of each break + # Length of each break length = nsamples // nbreaks - sublength = nsamples // (2*nbreaks) + sublength = nsamples // (2 * nbreaks) while 1: for pos in range(nbreaks): - ## 2 bolo in a pair will have the same break - ## to avoid differential gain + # 2 bolo in a pair will have the same break + # to avoid differential gain end_points = state.normal(mean, std, size=1) end_points = np.tile(end_points, (2, 1)).T.reshape((2, 1)) - ## Assign values + # Assign values shift = pos * length if pos * length < nsamples: - gains[:, shift + sublength:shift + 2 * sublength] = \ - end_points.reshape((len(end_points), 1)) + gains[ + :, shift + sublength: shift + 2 * sublength + ] = end_points.reshape((len(end_points), 1)) else: continue - if sign == 'opposite': + if sign == "opposite": gains[1] = 2 - gains[1] yield gains -def linear_function(nbolos, nsamples, mean=1, std=0.05, - nbreaks=1, sign='same', seed=0): + +def linear_function( + nbolos, nsamples, mean=1, std=0.05, nbreaks=1, sign="same", seed=0 +): """ Generate linear functions for each bolometer from 1 to a values drawn from N(mean, std). The full timestream is broken into nbreaks @@ -613,38 +655,43 @@ def linear_function(nbolos, nsamples, mean=1, std=0.05, >>> gains = linear_function(nbolos, nsamples, nbreaks=1, sign='opposite') >>> assert gains[0][0] == 2 - gains[1][0] """ - ## Fix the seed + # Fix the seed state = np.random.RandomState(seed) - ## Initialise gains to ones + # Initialise gains to ones gains = np.ones((nbolos, nsamples)) - ## Length of each break + # Length of each break length = nsamples // nbreaks for pos in range(nbreaks): - ## 2 bolo in a pair will have the same break to avoid differential gain - end_points = state.normal(mean, std, size=int(nbolos/2)) - end_points = np.tile( - end_points, (2, 1)).T.reshape((2*len(end_points), 1)) + # 2 bolo in a pair will have the same break to avoid differential gain + end_points = state.normal(mean, std, size=int(nbolos / 2)) + end_points = np.tile(end_points, (2, 1)).T.reshape((2 * len(end_points), 1)) - ## Assign values + # Assign values shift = pos * length if pos * length < nsamples: - gains[:, shift:shift + length] = np.array([np.interp( - range(shift, shift + length), - [shift, shift + length - 1], - [1, end[0]]) for end in end_points]) + gains[:, shift: shift + length] = np.array( + [ + np.interp( + range(shift, shift + length), + [shift, shift + length - 1], + [1, end[0]], + ) + for end in end_points + ] + ) else: continue - if sign == 'opposite': + if sign == "opposite": gains[1] = 2 - gains[1] return gains -def linear_function_gen(nsamples, mean=1, std=0.05, - nbreaks=1, sign='same', seed=0): + +def linear_function_gen(nsamples, mean=1, std=0.05, nbreaks=1, sign="same", seed=0): """ Generator returning linear functions for each bolometer from 1 to a values drawn from N(mean, std). The full timestream is broken into nbreaks @@ -687,38 +734,47 @@ def linear_function_gen(nsamples, mean=1, std=0.05, >>> assert g[0][0] == 2 - g[1][0] """ - ## Fix the seed + # Fix the seed state = np.random.RandomState(seed) - ## Initialise gains to ones + # Initialise gains to ones gains = np.ones((2, nsamples)) - ## Length of each break + # Length of each break length = nsamples // nbreaks while 1: for pos in range(nbreaks): - ## 2 bolo in a pair will have the same - ## break to avoid differential gain + # 2 bolo in a pair will have the same + # break to avoid differential gain end_points = state.normal(mean, std, size=1) end_points = np.tile(end_points, (2, 1)).T.reshape((2, 1)) - ## Assign values + # Assign values shift = pos * length if pos * length < nsamples: - gains[:, shift:shift + length] = np.array([np.interp( - range(shift, shift + length), - [shift, shift + length - 1], - [1, end[0]]) for end in end_points]) + gains[:, shift: shift + length] = np.array( + [ + np.interp( + range(shift, shift + length), + [shift, shift + length - 1], + [1, end[0]], + ) + for end in end_points + ] + ) else: continue - if sign == 'opposite': + if sign == "opposite": gains[1] = 2 - gains[1] yield gains -def get_kernel_coefficients(beamprm, pairlist,kernel_type='diff', nx=128, pix_size=None,basis='circular'): + +def get_kernel_coefficients( + beamprm, pairlist, kernel_type="diff", nx=128, pix_size=None, basis="circular" +): """ Generate beam expansion coefficients. @@ -801,21 +857,29 @@ def get_kernel_coefficients(beamprm, pairlist,kernel_type='diff', nx=128, pix_si >>> coeffs = get_kernel_coefficients(bm, pairlist, nx=32) """ - assert basis in ['circular','sumbeam'] - assert kernel_type in ['sum','diff','sumdiff','diffsum'] + assert basis in ["circular", "sumbeam"] + assert kernel_type in ["sum", "diff", "sumdiff", "diffsum"] if pix_size is None: - ## Go from sigma to FWHM - mean_FWHM_x = np.mean(beamprm.sig_1) / np.sqrt(8*np.log(2)) - mean_FWHM_y = np.mean(beamprm.sig_2) / np.sqrt(8*np.log(2)) + # Go from sigma to FWHM + mean_FWHM_x = np.mean(beamprm.sig_1) / np.sqrt(8 * np.log(2)) + mean_FWHM_y = np.mean(beamprm.sig_2) / np.sqrt(8 * np.log(2)) - ## 1/7th of the beam size. - pix_size = (mean_FWHM_x + mean_FWHM_y) / 2. / 7. + # 1/7th of the beam size. + pix_size = (mean_FWHM_x + mean_FWHM_y) / 2.0 / 7.0 - if basis=='circular': - ## Creates the a circular beam used to develop the sum/diff beam - fwhm=beamprm.fwhm*np.pi/180/60 + if basis == "circular": + # Creates the a circular beam used to develop the sum/diff beam + fwhm = beamprm.fwhm * np.pi / 180 / 60 xy2f = coordinates_on_grid(pix_size=pix_size, nx=nx) - basis_beam = gauss2d(xy2f,0,0, 1, fwhm/np.sqrt(8*np.log(2)), fwhm/np.sqrt(8*np.log(2)),0).reshape((nx, nx)) + basis_beam = gauss2d( + xy2f, + 0, + 0, + 1, + fwhm / np.sqrt(8 * np.log(2)), + fwhm / np.sqrt(8 * np.log(2)), + 0, + ).reshape((nx, nx)) out_diff = [] out_sum = [] @@ -823,21 +887,26 @@ def get_kernel_coefficients(beamprm, pairlist,kernel_type='diff', nx=128, pix_si for ct, cb in pairlist: summap, diffmap = construct_beammap(beamprm, ct, cb, nx, pix_size) if summap is not None and diffmap is not None: - if basis is not 'circular': basis_beam = summap - if 'diff' in kernel_type: + if basis != "circular": + basis_beam = summap + if "diff" in kernel_type: kernel_diff = split_deriv(basis_beam, diffmap, pix_size) - if 'sum' in kernel_type: + if "sum" in kernel_type: kernel_sum = split_deriv(basis_beam, summap, pix_size) else: - kernel = None kernel_sum = None - if 'diff' in kernel_type: out_diff.append(kernel_diff) - if 'sum' in kernel_type: out_sum.append(kernel_sum) - if (('diff' in kernel_type) and ('sum' in kernel_type)): + if "diff" in kernel_type: + out_diff.append(kernel_diff) + if "sum" in kernel_type: + out_sum.append(kernel_sum) + if ("diff" in kernel_type) and ("sum" in kernel_type): return np.array(out_diff), np.array(out_sum) else: - if 'diff'==kernel_type: return np.array(out_diff) - if 'sum'== kernel_type: return np.array(out_sum) + if "diff" == kernel_type: + return np.array(out_diff) + if "sum" == kernel_type: + return np.array(out_sum) + def split_deriv(sumbeam, diffbeam, pix_size): """ @@ -888,6 +957,7 @@ def split_deriv(sumbeam, diffbeam, pix_size): return x + def xderiv(m, pix_size): """ Perform the derivative of a 2d map with respect to x coordinate. @@ -909,10 +979,11 @@ def xderiv(m, pix_size): >>> m = np.ones((10, 10)) * np.arange(10) >>> dmx = xderiv(m, pix_size=0.5/60.*np.pi/180.) """ - kernel = np.array([[1, 0, -1]]) / (2.0*pix_size) - v = signal.convolve2d(m, kernel, mode='same') + kernel = np.array([[1, 0, -1]]) / (2.0 * pix_size) + v = signal.convolve2d(m, kernel, mode="same") return -v + def yderiv(m, pix_size): """ Perform the derivative of a 2d map with respect to y coordinate. @@ -934,10 +1005,11 @@ def yderiv(m, pix_size): >>> m = np.ones((10, 10)) * np.arange(10) >>> dmy = yderiv(m.T, pix_size=0.5/60.*np.pi/180.) """ - kernel = np.array([[1, 0, -1]]).T / (2.0*pix_size) - v = signal.convolve2d(m, kernel, mode='same') + kernel = np.array([[1, 0, -1]]).T / (2.0 * pix_size) + v = signal.convolve2d(m, kernel, mode="same") return -v + def derivs(m, pix_size): """ Compute full 1st and 2nd derivatives of a 2d map (flat sky). @@ -969,6 +1041,7 @@ def derivs(m, pix_size): return np.array((m00, m10, m01, m11, m20, m02)) + def fixspin(K, spin): """ Selects derivatives of only a given spin (by default we do not compute @@ -1011,21 +1084,21 @@ def fixspin(K, spin): e20 = 0.0 e02 = 0.0 - if '0' in spin: + if "0" in spin: e00 += d00 e10 += 0.0 e01 += 0.0 e11 += 0.0 e20 += x e02 += x - if '1' in spin: + if "1" in spin: e00 += 0.0 e10 += d10 e01 += d01 e11 += 0.0 e20 += 0.0 e02 += 0.0 - if '2' in spin: + if "2" in spin: e00 += 0.0 e10 += 0.0 e01 += 0.0 @@ -1035,6 +1108,7 @@ def fixspin(K, spin): return np.array((e00, e10, e01, e11, e20, e02)) + def rotate_deriv(K, theta): """ Rotate a vector of derivative kernel coefficients K by an angle theta. @@ -1082,24 +1156,35 @@ def rotate_deriv(K, theta): return es -def fix_kernel_type(K,kernel_type='diffnomonopole'): - if 'sum' in kernel_type: - ## subtracts monopole part already included - ## in the input timestream + +def fix_kernel_type(K, kernel_type="diffnomonopole"): + if "sum" in kernel_type: + # subtracts monopole part already included + # in the input timestream K[:, 0] -= 1 - if 'nomonopole' in kernel_type: - ## Force monopole to 0 when using leakage beam + if "nomonopole" in kernel_type: + # Force monopole to 0 when using leakage beam K[:, 0] = 0.0 - ## Flips the sign of dtheta terms. + # Flips the sign of dtheta terms. K[:, 1] *= -1 K[:, 3] *= -1 return K -def waferts_add_diffbeam(waferts, point_matrix, beam_orientation, - intensity_derivatives, diffbeam_kernels, - pairlist, spins='012',kernel_type='diff_nomonopole', - pol_derivatives=None,pol_angle=None,spins_pol='012', - Usign=1): + +def waferts_add_diffbeam( + waferts, + point_matrix, + beam_orientation, + intensity_derivatives, + diffbeam_kernels, + pairlist, + spins="012", + kernel_type="diff_nomonopole", + pol_derivatives=None, + pol_angle=None, + spins_pol="012", + Usign=1, +): """ Modify timestreams by injecting T->P leakage from beam mismatch. Note that timestreams are modified on-the-fly. @@ -1183,15 +1268,15 @@ def waferts_add_diffbeam(waferts, point_matrix, beam_orientation, ... pairlist, spins='012') """ - assert kernel_type in ['sum','diff','sum_nomonopole','diff_nomonopole'] + assert kernel_type in ["sum", "diff", "sum_nomonopole", "diff_nomonopole"] - ## Just to preserve it - K = diffbeam_kernels + 0. + # Just to preserve it + K = diffbeam_kernels + 0.0 - K = fix_kernel_type(K,kernel_type) + K = fix_kernel_type(K, kernel_type) - if 'sum' in kernel_type: - ## adds I->I in sum timestream and P->P in diff timestream + if "sum" in kernel_type: + # adds I->I in sum timestream and P->P in diff timestream bottom_sign = -1 else: bottom_sign = 1 @@ -1199,41 +1284,51 @@ def waferts_add_diffbeam(waferts, point_matrix, beam_orientation, for i in range(len(K)): K[i] = fixspin(K[i], spins) - diffbeamleak = np.zeros((int(waferts.shape[0]/2), waferts.shape[1])) + diffbeamleak = np.zeros((int(waferts.shape[0] / 2), waferts.shape[1])) - ## adds I leakage. if kernel_type is "diff" adds I->P leakage, otherwise - ## adds I->I leakage + # adds I leakage. if kernel_type is "diff" adds I->P leakage, otherwise + # adds I->I leakage diffbeam_map2tod( - diffbeamleak, - intensity_derivatives, - point_matrix, - beam_orientation, - K) + diffbeamleak, intensity_derivatives, point_matrix, beam_orientation, K + ) waferts[::2] += diffbeamleak - waferts[1::2] -= bottom_sign*diffbeamleak - - ## adds P leakage. If kernel_type is "diff" adds P->I leakage, otherwise - ## adds P->P leakage - if ((pol_derivatives is not None) and (pol_angle is not None)): - if ((spins_pol!=spins) and (pol_derivatives is not None)): - ## initialize a new copy of polarization kernels withot modifying - ## the input - K_pol = fix_kernel_type(diffbeam_kernels+0.,kernel_type) + waferts[1::2] -= bottom_sign * diffbeamleak + + # adds P leakage. If kernel_type is "diff" adds P->I leakage, otherwise + # adds P->P leakage + if (pol_derivatives is not None) and (pol_angle is not None): + if (spins_pol != spins) and (pol_derivatives is not None): + # initialize a new copy of polarization kernels withot modifying + # the input + K_pol = fix_kernel_type(diffbeam_kernels + 0.0, kernel_type) K_pol[i] = fixspin(K_pol[i], spins_pol) - else: K_pol = K - diffbeamleak[:,:] = 0.0 + else: + K_pol = K + diffbeamleak[:, :] = 0.0 diffbeam_map2tod( diffbeamleak, pol_derivatives, point_matrix, beam_orientation, - K_pol,pol_angle,Usign) + K_pol, + pol_angle, + Usign, + ) waferts[::2] += diffbeamleak - waferts[1::2] += bottom_sign*diffbeamleak - -def diffbeam_map2tod(out, signal_derivatives,point_matrix, beam_orientation, diffbeam_kernels, pol_angle=None,Usign=1): + waferts[1::2] += bottom_sign * diffbeamleak + + +def diffbeam_map2tod( + out, + signal_derivatives, + point_matrix, + beam_orientation, + diffbeam_kernels, + pol_angle=None, + Usign=1, +): """ Scan the leakage maps to generate polarisation timestream for channel ch, using the differential beam model kernel. @@ -1276,28 +1371,30 @@ def diffbeam_map2tod(out, signal_derivatives,point_matrix, beam_orientation, dif i1d = point_matrix[ipix, :] okpointing = i1d != -1 io = i1d[okpointing] - k = rotate_deriv( - diffbeam_kernels[ipix], - -beam_orientation[ipix, okpointing]) - if (pol_angle is not None): - #compute cosine and sine polarization angles to save time - psi = pol_angle[ipix,:] - c2pol_ang = np.cos(2*psi[okpointing]) - s2pol_ang = np.sin(2*psi[okpointing]) + k = rotate_deriv(diffbeam_kernels[ipix], -beam_orientation[ipix, okpointing]) + if pol_angle is not None: + # compute cosine and sine polarization angles to save time + psi = pol_angle[ipix, :] + c2pol_ang = np.cos(2 * psi[okpointing]) + s2pol_ang = np.sin(2 * psi[okpointing]) if signal_derivatives is not None: for coeff in range(len(k)): # Q derivatives - out[ipix, okpointing] += (signal_derivatives[0][coeff, io]*k[coeff])*c2pol_ang + der0 = signal_derivatives[0][coeff, io] + out[ipix, okpointing] += der0 * k[coeff] * c2pol_ang # U derivatives - out[ipix, okpointing] += (signal_derivatives[1][coeff, io]*k[coeff])*s2pol_ang*Usign - ## Add the leakage on-the-fly + der1 = signal_derivatives[1][coeff, io] + out[ipix, okpointing] += der1 * k[coeff] * s2pol_ang * Usign + # Add the leakage on-the-fly else: for coeff in range(len(k)): # T derivatives - out[ipix, okpointing] += signal_derivatives[coeff, io]*k[coeff] + out[ipix, okpointing] += signal_derivatives[coeff, io] * k[coeff] + if __name__ == "__main__": import doctest + if np.__version__ >= "1.14.0": np.set_printoptions(legacy="1.13") doctest.testmod() From b548ae96c2cd710f1408891269d215960f58e133 Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 10:24:57 +0100 Subject: [PATCH 05/15] Clean detector_pointing --- s4cmb/detector_pointing.py | 285 +++++++++++++++++++++++-------------- 1 file changed, 182 insertions(+), 103 deletions(-) diff --git a/s4cmb/detector_pointing.py b/s4cmb/detector_pointing.py index 0b95532..99ffa14 100755 --- a/s4cmb/detector_pointing.py +++ b/s4cmb/detector_pointing.py @@ -13,24 +13,26 @@ """ from __future__ import division, absolute_import, print_function -import healpy as hp +import os import numpy as np from numpy import cos from numpy import sin from numpy import tan from pyslalib import slalib from astropy.utils import iers +from astropy.time import Time import warnings from s4cmb.detector_pointing_f import detector_pointing_f -sec2deg = 360.0/86400.0 +sec2deg = 360.0 / 86400.0 d2r = np.pi / 180.0 ASTROMETRIC_GEOCENTRIC = 0 APPARENT_GEOCENTRIC = 1 APPARENT_TOPOCENTRIC = 2 + def get_ut1utc(ut1utc_fn, mjd): - """ Return the time correction to UTC. + """Return the time correction to UTC. Those are computed from 01 01 2012 (MJD=55927) to the latest MJD available in the s4cmb/data/ut1utc.ephem file. @@ -70,17 +72,32 @@ def get_ut1utc(ut1utc_fn, mjd): mjd_tmp = mjd warn_me = False if warn_me: - warnings.warn('MJD outside ut1utc file date range %.6f -- %.6f. Default corrections applied. Consider updating the ut1utc file.'%(umjds[0],umjds[-1])) + warnings.warn( + """ + MJD outside ut1utc file date range {:.6f} -- {:.6f}. + Default corrections applied. Consider updating the ut1utc file. + """.format(umjds[0], umjds[-1]) + ) uindex = np.searchsorted(umjds, mjd_tmp) ut1utc = ut1utcs[uindex] return ut1utc -class Pointing(): + +class Pointing: """ Class to handle detector pointing """ - def __init__(self, az_enc, el_enc, time, value_params, - allowed_params='ia ie ca an aw', - ra_src=0.0, dec_src=0.0, lat=-22.958, - ut1utc_fn='s4cmb/data/ut1utc.ephem'): + + def __init__( + self, + az_enc, + el_enc, + time, + value_params, + allowed_params="ia ie ca an aw", + ra_src=0.0, + dec_src=0.0, + lat=-22.958, + ut1utc_fn="s4cmb/data/ut1utc.ephem", + ): """ Apply pointing model with parameters `value_params` and names `allowed_params` to encoder az,el. Order of terms is @@ -158,7 +175,7 @@ def __init__(self, az_enc, el_enc, time, value_params, self.ut1utc = get_ut1utc(self.ut1utc_fn, self.time[0]) - ## Initialise the object + # Initialise the object self.az, self.el = self.apply_pointing_model() self.azel2radec() @@ -173,66 +190,94 @@ def apply_pointing_model(self): el : 1d array The corrected elevation in arcminutes. """ - assert len(self.value_params) == len(self.allowed_params.split()), \ - AssertionError("Vector containing parameters " + - "(value_params) has to have the same " + - "length than the vector containing names " + - "(allowed_params).") - - ## Here are many parameters defining a pointing model. - ## Of course, we do not use all of them. They are zero by default, - ## and only those specified by the user will be used. - params = {p: 0.0 for p in ['an', 'aw', 'an2', 'aw2', 'an4', - 'aw4', 'npae', 'ca', 'ia', 'ie', 'tf', - 'tfs', 'ref', 'dt', 'elt', 'ta1', - 'te1', 'sa', 'se', 'sa2', - 'se2', 'sta', 'ste', 'sta2', 'ste2']} + assert len(self.value_params) == len( + self.allowed_params.split() + ), AssertionError( + """Vector containing parameters + (value_params) has to have the same + length than the vector containing names + (allowed_params).""" + ) + + # Here are many parameters defining a pointing model. + # Of course, we do not use all of them. They are zero by default, + # and only those specified by the user will be used. + params = { + p: 0.0 + for p in [ + "an", + "aw", + "an2", + "aw2", + "an4", + "aw4", + "npae", + "ca", + "ia", + "ie", + "tf", + "tfs", + "ref", + "dt", + "elt", + "ta1", + "te1", + "sa", + "se", + "sa2", + "se2", + "sta", + "ste", + "sta2", + "ste2", + ] + } for param in params: if param in self.allowed_params.split(): index = self.allowed_params.split().index(param) params[param] = self.value_params[index] - params['dt'] *= sec2deg + params["dt"] *= sec2deg - ## Azimuth - azd = -params['an'] * sin(self.az_enc) * sin(self.el_enc) - azd -= params['aw'] * cos(self.az_enc) * sin(self.el_enc) + # Azimuth + azd = -params["an"] * sin(self.az_enc) * sin(self.el_enc) + azd -= params["aw"] * cos(self.az_enc) * sin(self.el_enc) - azd -= -params['an2'] * sin(2 * self.az_enc) * sin(self.el_enc) - azd -= params['aw2'] * cos(2 * self.az_enc) * sin(self.el_enc) + azd -= -params["an2"] * sin(2 * self.az_enc) * sin(self.el_enc) + azd -= params["aw2"] * cos(2 * self.az_enc) * sin(self.el_enc) - azd -= -params['an4'] * sin(4 * self.az_enc) * sin(self.el_enc) - azd -= params['aw4'] * cos(4 * self.az_enc) * sin(self.el_enc) + azd -= -params["an4"] * sin(4 * self.az_enc) * sin(self.el_enc) + azd -= params["aw4"] * cos(4 * self.az_enc) * sin(self.el_enc) - azd += params['npae'] * sin(self.el_enc) - azd -= params['ca'] - azd += params['ia'] * cos(self.el_enc) + azd += params["npae"] * sin(self.el_enc) + azd -= params["ca"] + azd += params["ia"] * cos(self.el_enc) - azd += params['dt'] * ( - -sin(self.lat) + cos(self.az_enc) * - cos(self.lat) * tan(self.el_enc)) + azd += params["dt"] * ( + -sin(self.lat) + cos(self.az_enc) * cos(self.lat) * tan(self.el_enc) + ) - ## Elevation - eld = params['an'] * cos(self.az_enc) - eld -= params['aw'] * sin(self.az_enc) - eld -= params['an2'] * cos(2 * self.az_enc) - eld -= params['aw2'] * sin(2 * self.az_enc) - eld -= params['an4'] * cos(4 * self.az_enc) - eld -= params['aw4'] * sin(4 * self.az_enc) + # Elevation + eld = params["an"] * cos(self.az_enc) + eld -= params["aw"] * sin(self.az_enc) + eld -= params["an2"] * cos(2 * self.az_enc) + eld -= params["aw2"] * sin(2 * self.az_enc) + eld -= params["an4"] * cos(4 * self.az_enc) + eld -= params["aw4"] * sin(4 * self.az_enc) - eld -= params['ie'] - eld += params['tf'] * cos(self.el_enc) - eld += params['tfs'] * sin(self.el_enc) - eld -= params['ref'] / tan(self.el_enc) + eld -= params["ie"] + eld += params["tf"] * cos(self.el_enc) + eld += params["tfs"] * sin(self.el_enc) + eld -= params["ref"] / tan(self.el_enc) - eld += -params['dt'] * cos(self.lat) * sin(self.az_enc) + eld += -params["dt"] * cos(self.lat) * sin(self.az_enc) - eld += params['elt'] * (self.time - np.min(self.time)) + eld += params["elt"] * (self.time - np.min(self.time)) - ## Convert back in radian and apply to the encoder values. - azd *= np.pi / (180.0 * 60.) - eld *= np.pi / (180.0 * 60.) + # Convert back in radian and apply to the encoder values. + azd *= np.pi / (180.0 * 60.0) + eld *= np.pi / (180.0 * 60.0) azd /= np.cos(self.el_enc) @@ -265,13 +310,13 @@ def azel2radec(self): self.meanpa = np.median(v_pa) - self.quaternion = Quaternion(v_ra, v_dec, v_pa, - v_ra_src, v_dec_src) + self.quaternion = Quaternion(v_ra, v_dec, v_pa, v_ra_src, v_dec_src) q = self.quaternion.offset_radecpa_makequat() - assert q.shape == (self.az.size, 4), \ - AssertionError("Wrong size for the quaternions!") + assert q.shape == (self.az.size, 4), AssertionError( + "Wrong size for the quaternions!" + ) self.q = q @@ -291,7 +336,7 @@ def azel2radecpa(self): >>> print(round(ra[0], 2), round(dec[0], 2), round(pa[0], 2)) 0.56 0.67 3.13 """ - ## TODO pass lon, lat, etc from the ScanningStrategy module! + # TODO pass lon, lat, etc from the ScanningStrategy module! converter = Azel2Radec(self.time[0], self.ut1utc) vconv = np.vectorize(converter.azel2radecpa) ra, dec, pa = vconv(self.time, self.az, self.el) @@ -312,7 +357,7 @@ def radec2azel(self): >>> assert np.all(np.round(az[2:4],2) == np.round(pointing.az[2:4],2)) >>> assert np.all(np.round(el[2:4],2) == np.round(pointing.el[2:4],2)) """ - ## TODO pass lon, lat, etc from the ScanningStrategy module! + # TODO pass lon, lat, etc from the ScanningStrategy module! converter = Azel2Radec(self.time[0], self.ut1utc) vconv = np.vectorize(converter.radec2azel) az, el = vconv(self.time, self.ra, self.dec) @@ -339,15 +384,25 @@ def offset_detector(self, azd, eld): pa : 1d array Parallactic angle in radian. """ - ra, dec, pa = self.quaternion.offset_radecpa_applyquat( - self.q, -azd, -eld) + ra, dec, pa = self.quaternion.offset_radecpa_applyquat(self.q, -azd, -eld) return ra, dec, pa + class Azel2Radec(object): """ Class to handle az/el <-> ra/dec conversion """ - def __init__(self, mjd, ut1utc, - lon=-67.786, lat=-22.958, height=5200., - pressure=533.29, temp=273.15, humidity=0.1, epequi=2000.0): + + def __init__( + self, + mjd, + ut1utc, + lon=-67.786, + lat=-22.958, + height=5200.0, + pressure=533.29, + temp=273.15, + humidity=0.1, + epequi=2000.0, + ): """ This class is mostly a wrapper around slalib. The default parameters correspond to the Polarbear observation site. @@ -404,10 +459,20 @@ def updateaoprms(self, mjd, lapserate=0.0065, xpm=0.0, ypm=0.0): """ wavelength = (299792458.0 / 150.0e9) * 1e6 - self.aoprms = slalib.sla_aoppa(mjd, self.ut1utc, self.lon, self.lat, - self.height, xpm, ypm, self.temp, - self.pressure, self.humidity, - wavelength, lapserate) + self.aoprms = slalib.sla_aoppa( + mjd, + self.ut1utc, + self.lon, + self.lat, + self.height, + xpm, + ypm, + self.temp, + self.pressure, + self.humidity, + wavelength, + lapserate, + ) def azel2radecpa(self, mjd, az, el): """ @@ -436,9 +501,9 @@ def azel2radecpa(self, mjd, az, el): amprms = slalib.sla_mappa(self.epequi, mjd) self.aoprms = slalib.sla_aoppat(mjd, self.aoprms) - ra_app1, dec_app1 = slalib.sla_oapqk('a', az, zd + 1e-8, self.aoprms) + ra_app1, dec_app1 = slalib.sla_oapqk("a", az, zd + 1e-8, self.aoprms) ra1, dec1 = slalib.sla_ampqk(ra_app1, dec_app1, amprms) - ra_app2, dec_app2 = slalib.sla_oapqk('a', az, zd - 1e-8, self.aoprms) + ra_app2, dec_app2 = slalib.sla_oapqk("a", az, zd - 1e-8, self.aoprms) ra2, dec2 = slalib.sla_ampqk(ra_app2, dec_app2, amprms) pa = slalib.sla_dbear(ra1, dec1, ra2, dec2) ra = 0.5 * (ra1 + ra2) @@ -473,8 +538,10 @@ def radec2azel(self, mjd, ra, dec): el = np.pi / 2 - zd return az, el -class Quaternion(): + +class Quaternion: """ Class to handle quaternions """ + def __init__(self, ra, dec, pa, v_ra_src, v_dec_src): """ Once you have RA/Dec coordinates of the reference detector, one @@ -563,11 +630,12 @@ def offset_radecpa_applyquat(self, q, azd, eld): assert seq.shape[1] == 4, AssertionError("Wrong size!") - n = seq.shape[0] + # n = seq.shape[0] phi, theta, psi = quat_to_radecpa_fortran(seq) return psi, -theta, -phi + def radec2thetaphi(ra, dec): """ Correspondance between RA/Dec and theta/phi coordinate systems. @@ -596,6 +664,7 @@ def radec2thetaphi(ra, dec): phi = ra return theta, phi + def mult(p, q): """ Multiply arrays of quaternions, @@ -639,11 +708,11 @@ def mult(p, q): pq[:, 3] = ps * qs pq[:, 3] -= arraylist_dot(pv, qv).flatten() - pq[:, :3] = ps[:, np.newaxis] * qv + \ - pv * qs[:, np.newaxis] + np.cross(pv, qv) + pq[:, :3] = ps[:, np.newaxis] * qv + pv * qs[:, np.newaxis] + np.cross(pv, qv) return pq + def mult_fortran(p, q): """ Inline version for when p is an array of quaternions @@ -687,6 +756,7 @@ def mult_fortran(p, q): pq = pq.reshape(shape) return pq + def arraylist_dot(a, b): """ Dot product of ndarrays. @@ -713,6 +783,7 @@ def arraylist_dot(a, b): else: return np.sum(a * b, axis=1)[:, np.newaxis] + def euler_quatx(alpha): """ Generate quaternion units along x axis @@ -734,6 +805,7 @@ def euler_quatx(alpha): s = np.sin(alpha * 0.5) return np.array([s, z, z, c]).T + def euler_quaty(alpha): """ Generate quaternion units along y axis @@ -755,6 +827,7 @@ def euler_quaty(alpha): s = np.sin(alpha * 0.5) return np.array([z, s, z, c]).T + def euler_quatz(alpha): """ Generate quaternion units along z axis @@ -776,6 +849,7 @@ def euler_quatz(alpha): s = np.sin(alpha * 0.5) return np.array([z, z, s, c]).T + def quat_to_radecpa_fortran(seq): """ Routine to compute phi/theta/psi from a sequence @@ -800,10 +874,10 @@ def quat_to_radecpa_fortran(seq): theta = np.zeros_like(q0) psi = np.zeros_like(q0) - detector_pointing_f.quat_to_radecpa_fortran_f( - q0, q1, q2, q3, phi, theta, psi, n) + detector_pointing_f.quat_to_radecpa_fortran_f(q0, q1, q2, q3, phi, theta, psi, n) return phi, theta, psi + def quat_to_radecpa_python(seq): """ Routine to compute phi/theta/psi from a sequence @@ -824,13 +898,12 @@ def quat_to_radecpa_python(seq): """ q1, q2, q3, q0 = seq.T - phi = np.arctan2(2 * (q0 * q1 + q2 * q3), - 1. - 2. * (q1 * q1 + q2 * q2)) + phi = np.arctan2(2 * (q0 * q1 + q2 * q3), 1.0 - 2.0 * (q1 * q1 + q2 * q2)) theta = np.arcsin(2 * (q0 * q2 - q3 * q1)) - psi = np.arctan2(2 * (q0 * q3 + q1 * q2), - 1. - 2. * (q2 * q2 + q3 * q3)) + psi = np.arctan2(2 * (q0 * q3 + q1 * q2), 1.0 - 2.0 * (q2 * q2 + q3 * q3)) return phi, theta, psi + def load_fake_pointing(): """ Load fake pointing parameters for testing purposes. @@ -850,17 +923,16 @@ def load_fake_pointing(): Encoder time (UTC) in mjd """ - allowed_params = 'ia ie ca an aw' - value_params = [10.28473073, 8.73953334, -15.59771781, - -0.50977716, 0.10858016] - az_enc = np.array([np.sin(2 * np.pi * i / 100) - for i in range(100)]) + allowed_params = "ia ie ca an aw" + value_params = [10.28473073, 8.73953334, -15.59771781, -0.50977716, 0.10858016] + az_enc = np.array([np.sin(2 * np.pi * i / 100) for i in range(100)]) el_enc = np.ones(100) * 0.5 - time = np.array([56293 + t/84000. for t in range(100)]) + time = np.array([56293 + t / 84000.0 for t in range(100)]) return allowed_params, value_params, az_enc, el_enc, time -def update_ut1utc(fname='ut1utc_user.ephem',begin_date='2012-01-01',end_date=None): + +def update_ut1utc(fname="ut1utc_user.ephem", begin_date="2012-01-01", end_date=None): """ Creates an updated database file with time correction to UTC. @@ -885,31 +957,38 @@ def update_ut1utc(fname='ut1utc_user.ephem',begin_date='2012-01-01',end_date=Non """ iers_table = iers.IERS_Auto.open() begin_time = Time(begin_date) - mask = (iers_table['MJD'].value>=begin_time.mjd) + mask = iers_table["MJD"].value >= begin_time.mjd if end_date is not None: - end_time = Time(end_time) - mask &= (iers_table['MJD'].value<=end_time.mjd) - y=iers_table['year'][mask] - m=iers_table['month'][mask] - d=iers_table['day'][mask] - mjd=iers_table['MJD'].value[mask] - ut1utc=iers_table['UT1_UTC'][mask] + end_time = Time(end_date) + mask &= iers_table["MJD"].value <= end_time.mjd + y = iers_table["year"][mask] + m = iers_table["month"][mask] + d = iers_table["day"][mask] + mjd = iers_table["MJD"].value[mask] + ut1utc = iers_table["UT1_UTC"][mask] # Corrects MJD according to astropy doc - y[mjd<=51543]+=1900 - y[mjd>=51544]+=2000 + y[mjd <= 51543] += 1900 + y[mjd >= 51544] += 2000 try: - fname_path = os.path.join(os.environ['s4cmbPATH'],'s4cmb/data') - except: - fname_path = './' - - np.savetxt(os.path.join(fname_path,fname),np.column_stack((y,m,d,mjd,ut1utc)),fmt="%d-%d-%d\t%.6f\t%.6f") + fname_path = os.path.join(os.environ["s4cmbPATH"], "s4cmb/data") + except KeyError as e: + print(str(e), "Saving the current directory instead") + fname_path = "./" + + np.savetxt( + os.path.join(fname_path, fname), + np.column_stack((y, m, d, mjd, ut1utc)), + fmt="%d-%d-%d\t%.6f\t%.6f", + ) iers_table.close() return + if __name__ == "__main__": import doctest + if np.__version__ >= "1.14.0": np.set_printoptions(legacy="1.13") doctest.testmod() From 4b8e4c79c4eafa8c238d43aea9db72e39a148d23 Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 10:40:42 +0100 Subject: [PATCH 06/15] Clean scanning strategy --- s4cmb/scanning_strategy.py | 807 ++++++++++++++++++++++--------------- 1 file changed, 493 insertions(+), 314 deletions(-) diff --git a/s4cmb/scanning_strategy.py b/s4cmb/scanning_strategy.py index 340c17e..171e1a1 100755 --- a/s4cmb/scanning_strategy.py +++ b/s4cmb/scanning_strategy.py @@ -15,18 +15,28 @@ from pyslalib import slalib -## numerical constants -radToDeg = 180. / np.pi +# numerical constants +radToDeg = 180.0 / np.pi sidDayToSec = 86164.0905 -class ScanningStrategy(): + +class ScanningStrategy: """ Class to handle the scanning strategy of the telescope """ - def __init__(self, nces=12, start_date='2013/1/1 00:00:00', - telescope_longitude='-67:46.816', - telescope_latitude='-22:56.396', telescope_elevation=5200., - name_strategy='deep_patch', sampling_freq=30., sky_speed=0.4, - ut1utc_fn='s4cmb/data/ut1utc.ephem', language='python', - verbose=False): + + def __init__( + self, + nces=12, + start_date="2013/1/1 00:00:00", + telescope_longitude="-67:46.816", + telescope_latitude="-22:56.396", + telescope_elevation=5200.0, + name_strategy="deep_patch", + sampling_freq=30.0, + sky_speed=0.4, + ut1utc_fn="s4cmb/data/ut1utc.ephem", + language="python", + verbose=False, + ): """ A scanning strategy consists in defining the site of observation on earth for which we will make the observation, the region @@ -78,22 +88,28 @@ def __init__(self, nces=12, start_date='2013/1/1 00:00:00', self.ut1utc_fn = ut1utc_fn if not os.path.isfile(self.ut1utc_fn): - url = 'http://tycho.usno.navy.mil/leapsec.html' - msg = 'The path {} does not point to a valid file! see ' + \ - 's4cmb/data/ut1utc.ephem provided with the package. ' + \ - 'For more information, see {}.' + url = "http://tycho.usno.navy.mil/leapsec.html" + msg = """ + The path {} does not point to a valid file! see + s4cmb/data/ut1utc.ephem provided with the package. + For more information, see {}. + """ raise Exception(msg.format(self.ut1utc_fn, url)) self.verbose = verbose self.telescope_location = self.define_telescope_location( - telescope_longitude, telescope_latitude, telescope_elevation) + telescope_longitude, telescope_latitude, telescope_elevation + ) self.define_boundary_of_scan() - def define_telescope_location(self, telescope_longitude='-67:46.816', - telescope_latitude='-22:56.396', - telescope_elevation=5200.): + def define_telescope_location( + self, + telescope_longitude="-67:46.816", + telescope_latitude="-22:56.396", + telescope_elevation=5200.0, + ): """ Routine to define the site of observation on earth for which positions are to be computed. The location of the Polarbear telescope @@ -160,33 +176,87 @@ def define_boundary_of_scan(self): ['deep_patch', 'shallow_patch', 'small_aperture', 'custom'] """ self.allowed_scanning_strategies = [ - 'deep_patch', - 'shallow_patch', - 'small_aperture', - 'custom'] - - if self.name_strategy == 'deep_patch': - self.elevation = [30.0, 45.5226, 47.7448, 49.967, - 52.1892, 54.4114, 56.6336, 58.8558, - 61.078, 63.3002, 65.5226, 35.2126] - - self.az_min = [134.2263, 162.3532, 162.3532, 162.3532, - 162.3532, 162.3532, 162.3532, 162.3532, - 162.3532, 162.3532, 162.3532, 204.7929] - - self.az_max = [154.2263, 197.3532, 197.3532, 197.3532, - 197.3532, 197.3532, 197.3532, 197.3532, - 197.3532, 197.3532, 197.3532, 224.7929] - - self.begin_LST = ['17:07:54.84', '22:00:21.76', '22:00:21.76', - '22:00:21.76', '22:00:21.76', '22:00:21.76', - '22:00:21.76', '22:00:21.76', '22:00:21.76', - '22:00:21.76', '22:00:21.76', '2:01:01.19'] - - self.end_LST = ['22:00:21.76', '02:01:01.19', '02:01:01.19', - '02:01:01.19', '02:01:01.19', '02:01:01.19', - '02:01:01.19', '02:01:01.19', '02:01:01.19', - '02:01:01.19', '02:01:01.19', '6:53:29.11'] + "deep_patch", + "shallow_patch", + "small_aperture", + "custom", + ] + + if self.name_strategy == "deep_patch": + self.elevation = [ + 30.0, + 45.5226, + 47.7448, + 49.967, + 52.1892, + 54.4114, + 56.6336, + 58.8558, + 61.078, + 63.3002, + 65.5226, + 35.2126, + ] + + self.az_min = [ + 134.2263, + 162.3532, + 162.3532, + 162.3532, + 162.3532, + 162.3532, + 162.3532, + 162.3532, + 162.3532, + 162.3532, + 162.3532, + 204.7929, + ] + + self.az_max = [ + 154.2263, + 197.3532, + 197.3532, + 197.3532, + 197.3532, + 197.3532, + 197.3532, + 197.3532, + 197.3532, + 197.3532, + 197.3532, + 224.7929, + ] + + self.begin_LST = [ + "17:07:54.84", + "22:00:21.76", + "22:00:21.76", + "22:00:21.76", + "22:00:21.76", + "22:00:21.76", + "22:00:21.76", + "22:00:21.76", + "22:00:21.76", + "22:00:21.76", + "22:00:21.76", + "2:01:01.19", + ] + + self.end_LST = [ + "22:00:21.76", + "02:01:01.19", + "02:01:01.19", + "02:01:01.19", + "02:01:01.19", + "02:01:01.19", + "02:01:01.19", + "02:01:01.19", + "02:01:01.19", + "02:01:01.19", + "02:01:01.19", + "6:53:29.11", + ] self.dec_min = None self.dec_max = None @@ -194,23 +264,33 @@ def define_boundary_of_scan(self): self.end_RA = None self.orientation = None - ## Center of the patch in RA/Dec - self.ra_mid = 0. + # Center of the patch in RA/Dec + self.ra_mid = 0.0 self.dec_mid = -57.5 - elif self.name_strategy == 'small_aperture': - self.elevation = [50.] * 6 - - self.az_min = [37.1, 282.8, 110.3, - 197.0, 110.3, 197.0] - - self.az_max = [77.7, 322.9, 163.0, - 250.2, 163.0, 250.2] - - self.begin_LST = ['08:30:00.00', '12:30:00.00', '16:30:00.00', - '02:30:00.00', '16:30:00.00', '22:30:00.00'] - - self.end_LST = ['12:29:59.32', '16:29:59.32', '02:29:59.30', - '08:29:59.98', '22:29:59.99', '08:29:59.30'] + elif self.name_strategy == "small_aperture": + self.elevation = [50.0] * 6 + + self.az_min = [37.1, 282.8, 110.3, 197.0, 110.3, 197.0] + + self.az_max = [77.7, 322.9, 163.0, 250.2, 163.0, 250.2] + + self.begin_LST = [ + "08:30:00.00", + "12:30:00.00", + "16:30:00.00", + "02:30:00.00", + "16:30:00.00", + "22:30:00.00", + ] + + self.end_LST = [ + "12:29:59.32", + "16:29:59.32", + "02:29:59.30", + "08:29:59.98", + "22:29:59.99", + "08:29:59.30", + ] self.dec_min = None self.dec_max = None @@ -218,35 +298,59 @@ def define_boundary_of_scan(self): self.end_RA = None self.orientation = None - ## Center of the patch in RA/Dec - self.ra_mid = 0. - self.dec_mid = 0. - elif self.name_strategy == 'shallow_patch': - self.elevation = [30., 30., 30., 30., 30., 30.] - - self.dec_min = ["-5:00:00", "-5:00:00", "-15:00:00", - "-15:00:00", "-60:00:00", "-60:00:00"] - - self.dec_max = ["20:00:00", "20:00:00", "20:00:00", - "20:00:00", "-15:00:00", "-15:00:00"] - - self.begin_RA = ["7:40:00", "7:40:00", "20:00:00", - "20:00:00", "20:00:00", "20:00:00"] - - self.end_RA = ["15:20:00", "15:20:00", "5:40:00", - "5:40:00", "5:40:00", "5:40:00"] - - self.orientation = ['west', 'east', 'west', 'east', 'west', 'east'] + # Center of the patch in RA/Dec + self.ra_mid = 0.0 + self.dec_mid = 0.0 + elif self.name_strategy == "shallow_patch": + self.elevation = [30.0, 30.0, 30.0, 30.0, 30.0, 30.0] + + self.dec_min = [ + "-5:00:00", + "-5:00:00", + "-15:00:00", + "-15:00:00", + "-60:00:00", + "-60:00:00", + ] + + self.dec_max = [ + "20:00:00", + "20:00:00", + "20:00:00", + "20:00:00", + "-15:00:00", + "-15:00:00", + ] + + self.begin_RA = [ + "7:40:00", + "7:40:00", + "20:00:00", + "20:00:00", + "20:00:00", + "20:00:00", + ] + + self.end_RA = [ + "15:20:00", + "15:20:00", + "5:40:00", + "5:40:00", + "5:40:00", + "5:40:00", + ] + + self.orientation = ["west", "east", "west", "east", "west", "east"] self.az_min = None self.az_max = None self.begin_LST = None self.end_LST = None - ## Center of the patch in RA/Dec - self.ra_mid = 0. - self.dec_mid = 0. - elif self.name_strategy == 'custom': + # Center of the patch in RA/Dec + self.ra_mid = 0.0 + self.dec_mid = 0.0 + elif self.name_strategy == "custom": self.elevation = None self.dec_min = None self.dec_max = None @@ -258,14 +362,16 @@ def define_boundary_of_scan(self): self.begin_LST = None self.end_LST = None - ## Center of the patch in RA/Dec - self.ra_mid = 0. - self.dec_mid = 0. + # Center of the patch in RA/Dec + self.ra_mid = 0.0 + self.dec_mid = 0.0 else: - raise ValueError("Only name_strategy = deep_patch or " + - "shallow_patch or custom are " + - "currently available. For another usage " + - "(advanced users), modify this routine.") + raise ValueError( + """Only name_strategy = deep_patch or + shallow_patch or custom are + currently available. For another usage + (advanced users), modify this routine.""" + ) def run_one_scan(self, scan_file, scan_number): """ @@ -284,104 +390,110 @@ def run_one_scan(self, scan_file, scan_number): Returns True if the scan has been generated, and False if the scan already exists on the disk. """ - ## Check if we have too much/enough information to make a scan + # Check if we have too much/enough information to make a scan msg = "You cannot specify azimuth and declination!" - assert (getattr(self, 'az_min') and not getattr(self, 'dec_min')) or \ - (not getattr(self, 'az_min') and getattr(self, 'dec_min')), msg - assert (getattr(self, 'az_max') and not getattr(self, 'dec_max')) or \ - (not getattr(self, 'az_max') and getattr(self, 'dec_max')), msg + assert (getattr(self, "az_min") and not getattr(self, "dec_min")) or ( + not getattr(self, "az_min") and getattr(self, "dec_min") + ), msg + assert (getattr(self, "az_max") and not getattr(self, "dec_max")) or ( + not getattr(self, "az_max") and getattr(self, "dec_max") + ), msg msg = "You cannot specify LST and RA!" - assert (getattr(self, 'begin_LST') and not getattr(self, 'begin_RA')) or \ - (not getattr(self, 'begin_LST') and getattr(self, 'begin_RA')), msg - assert (getattr(self, 'end_LST') and not getattr(self, 'end_RA')) or \ - (not getattr(self, 'end_LST') and getattr(self, 'end_RA')), msg - - if getattr(self, 'az_min') and getattr(self, 'az_max'): - msg = 'You need to define timing bounds!' - assert getattr(self, 'begin_LST') is not None, msg - assert getattr(self, 'end_LST') is not None, msg + assert (getattr(self, "begin_LST") and not getattr(self, "begin_RA")) or ( + not getattr(self, "begin_LST") and getattr(self, "begin_RA") + ), msg + assert (getattr(self, "end_LST") and not getattr(self, "end_RA")) or ( + not getattr(self, "end_LST") and getattr(self, "end_RA") + ), msg + + if getattr(self, "az_min") and getattr(self, "az_max"): + msg = "You need to define timing bounds!" + assert getattr(self, "begin_LST") is not None, msg + assert getattr(self, "end_LST") is not None, msg azLST = True RADEC = False if self.verbose: print("Using azimuth and LST bounds") - elif getattr(self, 'dec_min') and getattr(self, 'dec_max'): - msg = 'You need to define RA bounds!' - assert getattr(self, 'begin_RA') is not None, msg - assert getattr(self, 'end_RA') is not None, msg - msg = 'You need to define orientation of scan (east/west)!' - assert getattr(self, 'orientation') is not None, msg + elif getattr(self, "dec_min") and getattr(self, "dec_max"): + msg = "You need to define RA bounds!" + assert getattr(self, "begin_RA") is not None, msg + assert getattr(self, "end_RA") is not None, msg + msg = "You need to define orientation of scan (east/west)!" + assert getattr(self, "orientation") is not None, msg azLST = False RADEC = True if self.verbose: print("Using RA and Dec bounds") - ## Figure out the elevation to run the scan + # Figure out the elevation to run the scan el = self.elevation[scan_number] - ## Define the sampling rate in Hz + # Define the sampling rate in Hz sampling_freq = self.sampling_freq - ######################################################### - ## Define geometry of the scan - ######################################################### + ############################# + # Define geometry of the scan + ############################# if azLST: - ## Define geometry of the scan by figuring out the azimuth bounds - az_mean = ( - self.az_min[scan_number] + self.az_max[scan_number]) * 0.5 - az_throw = (self.az_max[scan_number] - - self.az_min[scan_number]) / np.cos(el / radToDeg) + # Define geometry of the scan by figuring out the azimuth bounds + az_mean = (self.az_min[scan_number] + self.az_max[scan_number]) * 0.5 + az_throw = (self.az_max[scan_number] - self.az_min[scan_number]) / np.cos( + el / radToDeg + ) elif RADEC: - ## If given bounds in declination, make bounds in azimuth - ## note there is no sanity checking here! - az_array = np.linspace(0., 180., endpoint=True, num=360) + # If given bounds in declination, make bounds in azimuth + # note there is no sanity checking here! + az_array = np.linspace(0.0, 180.0, endpoint=True, num=360) ra_array = np.zeros(az_array.shape) dec_array = np.zeros(az_array.shape) dec_min = ephem.degrees(self.dec_min[scan_number]) dec_max = ephem.degrees(self.dec_max[scan_number]) - if(self.orientation[scan_number] == 'west'): - az_array += 180. + if self.orientation[scan_number] == "west": + az_array += 180.0 for i in range(0, az_array.shape[0]): - ra_array[i], dec_array[i] = \ - self.telescope_location.radec_of( - az_array[i] / radToDeg, el / radToDeg) + ra_array[i], dec_array[i] = self.telescope_location.radec_of( + az_array[i] / radToDeg, el / radToDeg + ) az_allowed = np.asarray( - [az_array[i] for i in range(0, az_array.shape[0]) - if (dec_array[i] > dec_min and dec_array[i] < dec_max)]) - - if (az_allowed.shape[0] < 2): - ms = 'Invalid combination of declination bounds and elevation.' + [ + az_array[i] + for i in range(0, az_array.shape[0]) + if (dec_array[i] > dec_min and dec_array[i] < dec_max) + ] + ) + + if az_allowed.shape[0] < 2: + ms = "Invalid combination of declination bounds and elevation." print(ms) exit(1) az_max = np.max(az_allowed) az_min = np.min(az_allowed) az_mean = (az_min + az_max) * 0.5 - az_throw = (az_max - az_min) + az_throw = az_max - az_min - ######################################################### - ## Define the timing bounds! - ######################################################### + ############################# + # Define the timing bounds! + ############################# if azLST: - LST_now = float( - self.telescope_location.sidereal_time()) / (2 * np.pi) - begin_LST = float( - ephem.hours(self.begin_LST[scan_number])) / (2 * np.pi) - end_LST = float( - ephem.hours(self.end_LST[scan_number])) / (2 * np.pi) + LST_now = float(self.telescope_location.sidereal_time()) / (2 * np.pi) + begin_LST = float(ephem.hours(self.begin_LST[scan_number])) / (2 * np.pi) + end_LST = float(ephem.hours(self.end_LST[scan_number])) / (2 * np.pi) - if (begin_LST > end_LST): - begin_LST -= 1. + if begin_LST > end_LST: + begin_LST -= 1.0 - ## Reset the date to correspond to the sidereal time to start + # Reset the date to correspond to the sidereal time to start self.telescope_location.date -= ( - (LST_now - begin_LST) * sidDayToSec) * ephem.second + (LST_now - begin_LST) * sidDayToSec + ) * ephem.second - ## Figure out how long to run the scan for + # Figure out how long to run the scan for num_pts = int((end_LST - begin_LST) * sidDayToSec * sampling_freq) if RADEC: @@ -393,7 +505,8 @@ def run_one_scan(self, scan_file, scan_number): # Figure out where we are looking now ra_target, dec_target = self.telescope_location.radec_of( - az_mean / radToDeg, el / radToDeg) + az_mean / radToDeg, el / radToDeg + ) # Instantiate the targets target_min_ra._dec = dec_target @@ -404,138 +517,151 @@ def run_one_scan(self, scan_file, scan_number): # Compute initial RA target_min_ra.compute(self.telescope_location) target_max_ra.compute(self.telescope_location) - if(self.orientation[scan_number] == 'east'): - self.telescope_location.date = \ - self.telescope_location.next_rising(target_min_ra) + if self.orientation[scan_number] == "east": + self.telescope_location.date = self.telescope_location.next_rising( + target_min_ra + ) # Recompute coodinates in the light of change of date target_min_ra.compute(self.telescope_location) target_max_ra.compute(self.telescope_location) - ## Update number of time samples for the scan - num_pts = int( - (self.telescope_location.next_rising( - target_max_ra) - self.telescope_location.date) / - ephem.second * sampling_freq) + # Update number of time samples for the scan + ra_max_rising = self.telescope_location.next_rising(target_max_ra) + diff_ra = ra_max_rising - self.telescope_location.date + num_pts = int(diff_ra / ephem.second * sampling_freq) - if(self.orientation[scan_number] == 'west'): - self.telescope_location.date = \ - self.telescope_location.next_setting(target_min_ra) + if self.orientation[scan_number] == "west": + self.telescope_location.date = self.telescope_location.next_setting( + target_min_ra + ) # Recompute coodinates in the light of change of date target_min_ra.compute(self.telescope_location) target_max_ra.compute(self.telescope_location) - ## Update number of time samples for the scan - num_pts = int( - (self.telescope_location.next_setting( - target_max_ra) - self.telescope_location.date) / - ephem.second * sampling_freq) + # Update number of time samples for the scan + ra_max_setting = self.telescope_location.next_setting(target_max_ra) + diff_ra = ra_max_setting - self.telescope_location.date + num_pts = int(diff_ra / ephem.second * sampling_freq) - ## Run the scan! - pb_az_dir = 1. - upper_az = az_mean + az_throw / 2. - lower_az = az_mean - az_throw / 2. + # Run the scan! + pb_az_dir = 1.0 + upper_az = az_mean + az_throw / 2.0 + lower_az = az_mean - az_throw / 2.0 az_speed = self.sky_speed / np.cos(el / radToDeg) running_az = az_mean - ## Initialize arrays + # Initialize arrays pb_az_array = np.zeros(num_pts) pb_mjd_array = np.zeros(num_pts) pb_ra_array = np.zeros(num_pts) pb_dec_array = np.zeros(num_pts) pb_el_array = np.ones(num_pts) * el - ## Subscan boundaries + # Subscan boundaries pb_subscans = [] pb_subscans.append(0) - ## Loop over time samples + # Loop over time samples # begin_lst = str(self.telescope_location.sidereal_time()) # Pad scans 10 seconds on either side time_padding = 10.0 / 86400.0 - ## Start of the scan + # Start of the scan pb_az_array[0] = running_az pb_mjd_array[0] = date_to_mjd(self.telescope_location.date) - ## Initialize the time - scan_file['firstmjd'] = pb_mjd_array[0] - time_padding + # Initialize the time + scan_file["firstmjd"] = pb_mjd_array[0] - time_padding - ## Update before starting the loop + # Update before starting the loop running_az += az_speed * pb_az_dir / sampling_freq self.telescope_location.date += ephem.second / sampling_freq - if self.language == 'python': + if self.language == "python": for t in range(0, num_pts): - ## Set the Azimuth and time + # Set the Azimuth and time pb_az_array[t] = running_az - pb_ra_array[t], pb_dec_array[t] = \ - self.telescope_location.radec_of( - pb_az_array[t] * np.pi / 180., - pb_el_array[t] * np.pi / 180.) + pb_ra_array[t], pb_dec_array[t] = self.telescope_location.radec_of( + pb_az_array[t] * np.pi / 180.0, pb_el_array[t] * np.pi / 180.0 + ) - ## Case to change the direction of the scan - if(running_az > upper_az): - pb_az_dir = -1. + # Case to change the direction of the scan + if running_az > upper_az: + pb_az_dir = -1.0 pb_subscans.append(t) - elif(running_az < lower_az): - pb_az_dir = 1. + elif running_az < lower_az: + pb_az_dir = 1.0 pb_subscans.append(t) running_az += az_speed * pb_az_dir / sampling_freq - ## Increment the time by one second / sampling rate + # Increment the time by one second / sampling rate if t > 0: - pb_mjd_array[t] = pb_mjd_array[t-1] + \ - ephem.second / sampling_freq + pb_mjd_array[t] = ( + pb_mjd_array[t - 1] + ephem.second / sampling_freq + ) - ## Increment the time by one second / sampling rate + # Increment the time by one second / sampling rate self.telescope_location.date += ephem.second / sampling_freq - elif self.language == 'fortran': - second = 1./24./3600. + elif self.language == "fortran": + second = 1.0 / 24.0 / 3600.0 scanning_strategy_f.run_one_scan_f( - pb_az_array, pb_mjd_array, - running_az, upper_az, lower_az, az_speed, pb_az_dir, - second, sampling_freq, num_pts) - - ## Do not use that for precision - it truncates values + pb_az_array, + pb_mjd_array, + running_az, + upper_az, + lower_az, + az_speed, + pb_az_dir, + second, + sampling_freq, + num_pts, + ) + + # Do not use that for precision - it truncates values self.telescope_location.date += num_pts * ephem.second / sampling_freq - ## Save in file - scan_file['nces'] = self.nces - scan_file['CES'] = scan_number - scan_file['sample_rate'] = sampling_freq - scan_file['sky_speed'] = self.sky_speed - scan_file['lastmjd'] = pb_mjd_array[-1] + time_padding + # Save in file + scan_file["nces"] = self.nces + scan_file["CES"] = scan_number + scan_file["sample_rate"] = sampling_freq + scan_file["sky_speed"] = self.sky_speed + scan_file["lastmjd"] = pb_mjd_array[-1] + time_padding - scan_file['azimuth'] = pb_az_array * np.pi / 180 - scan_file['elevation'] = pb_el_array * np.pi / 180 - scan_file['clock-utc'] = pb_mjd_array + scan_file["azimuth"] = pb_az_array * np.pi / 180 + scan_file["elevation"] = pb_el_array * np.pi / 180 + scan_file["clock-utc"] = pb_mjd_array - scan_file['RA'] = pb_ra_array - scan_file['Dec'] = pb_dec_array + scan_file["RA"] = pb_ra_array + scan_file["Dec"] = pb_dec_array - scan_file['nts'] = len(pb_mjd_array) + scan_file["nts"] = len(pb_mjd_array) - scan_file['subscans'] = zip( - pb_subscans[:-1], np.array(pb_subscans[1:]) - 1) + scan_file["subscans"] = zip(pb_subscans[:-1], np.array(pb_subscans[1:]) - 1) if self.verbose: - print('+-----------------------------------+') - print(' CES starts at %s and finishes at %s' % ( - mjd_to_greg(scan_file['firstmjd']), - mjd_to_greg(scan_file['lastmjd']))) - print(' It lasts %.3f hours' % ( - (scan_file['lastmjd'] - scan_file['firstmjd']) * 24)) - print('+-----------------------------------+') - - ## Add one day before the next CES (to avoid conflict of time) + print("+-----------------------------------+") + print( + " CES starts at %s and finishes at %s" + % ( + mjd_to_greg(scan_file["firstmjd"]), + mjd_to_greg(scan_file["lastmjd"]), + ) + ) + print( + " It lasts %.3f hours" + % ((scan_file["lastmjd"] - scan_file["firstmjd"]) * 24) + ) + print("+-----------------------------------+") + + # Add one day before the next CES (to avoid conflict of time) self.telescope_location.date += 24 * ephem.second * 3600 - ## Add the scan into the instance + # Add the scan into the instance # self._update('scan{}'.format(scan_number), scan_file) return True @@ -586,20 +712,31 @@ def run(self): orientations (east/west) """ - ## Initialise the date and loop over CESes + # Initialise the date and loop over CESes self.telescope_location.date = self.start_date for CES_position in range(self.nces): - ## Initialise the starting date of observation - ## It will be updated then automatically - setattr(self, 'scan{}'.format(CES_position), {}) + # Initialise the starting date of observation + # It will be updated then automatically + setattr(self, "scan{}".format(CES_position), {}) # Create the scan strategy self.run_one_scan( - getattr(self, 'scan{}'.format(CES_position)), CES_position) - - def visualize_my_scan(self, nside, reso=6.9, xsize=900, rot=[0, -57.5], - nfid_bolometer=6000, fp_size=180., boost=1., - fullsky=False,flatsky=False,nest=False): + getattr(self, "scan{}".format(CES_position)), CES_position + ) + + def visualize_my_scan( + self, + nside, + reso=6.9, + xsize=900, + rot=[0, -57.5], + nfid_bolometer=6000, + fp_size=180.0, + boost=1.0, + fullsky=False, + flatsky=False, + nest=False, + ): """ Simple map-making: project time ordered data into sky maps for visualisation. It works only in pure python (i.e. if you set @@ -625,102 +762,140 @@ def visualize_my_scan(self, nside, reso=6.9, xsize=900, rot=[0, -57.5], Boost factor to artificially increase the number of hits. It doesn't change the shape of the survey (just the amplitude). fullsky : boolean - If True, visualising full sky. If False, visualising area defined by reso and my xsize. + If True, visualising full sky. If False, visualising area defined by + reso and my xsize. flatsky : boolean - If True, visualising projection in a square array with the specified resolution. If False, visualising using healpy functions. + If True, visualising projection in a square array with the specified + resolution. If False, visualising using healpy functions. nest : boolean If flatsky projection, using the defined nesting in projection. Outputs ---------- - * nhit: Sky map with cumulative hit counts. If flatsky is True, a square numpy array of size xsize x xsize is returned. Else, a 1D healpy array is returned. + nhit: Sky map with cumulative hit counts. + If flatsky is True, a square numpy array of size xsize x xsize is + returned. Else, a 1D healpy array is returned. Examples ---------- >>> import matplotlib - >>> matplotlib.use("Agg") ## remove this if you want to display + >>> matplotlib.use("Agg") # remove this if you want to display >>> scan = ScanningStrategy(sampling_freq=1., nces=1) >>> scan.run() >>> nhit = scan.visualize_my_scan(512) """ import pylab as pl - if self.language != 'python': - raise ValueError("Visualisation is available only in pure " + - "python because we do not provide (yet) " + - "RA and Dec in fortran. Relaunch " + - "using language='python' in the class " + - "ScanningStrategy.") + + if self.language != "python": + raise ValueError( + """Visualisation is available only in pure + python because we do not provide (yet) + RA and Dec in fortran. Relaunch + using language='python' in the class + ScanningStrategy.""" + ) npix = hp.pixelfunc.nside2npix(nside) nhit = np.zeros(npix) for scan_number in range(self.nces): - scan = getattr(self, 'scan{}'.format(scan_number)) + scan = getattr(self, "scan{}".format(scan_number)) - num_pts = len(scan['clock-utc']) + num_pts = len(scan["clock-utc"]) pix_global = hp.pixelfunc.ang2pix( - nside, (np.pi/2.) - scan['Dec'], scan['RA']) + nside, (np.pi / 2.0) - scan["Dec"], scan["RA"] + ) - ## Boresight pointing healpix maps + # Boresight pointing healpix maps nhit_loc = np.zeros(npix, dtype=np.int32) - scanning_strategy_f.mapmaking( - pix_global, nhit_loc, npix, num_pts) + scanning_strategy_f.mapmaking(pix_global, nhit_loc, npix, num_pts) - ## Fake large focal plane with many bolometers for visualisation. - nhit_loc = convolve_focalplane(nhit_loc, nfid_bolometer, - fp_size, boost) + # Fake large focal plane with many bolometers for visualisation. + nhit_loc = convolve_focalplane(nhit_loc, nfid_bolometer, fp_size, boost) nhit += nhit_loc if self.verbose: - print('Stats: nhits = {}/{} (fsky={}%), max hit = {}'.format( - len(nhit[nhit > 0]), - len(nhit), - round(len(nhit[nhit > 0])/len(nhit) * 100, 2), - int(np.max(nhit)))) - + print( + "Stats: nhits = {}/{} (fsky={}%), max hit = {}".format( + len(nhit[nhit > 0]), + len(nhit), + round(len(nhit[nhit > 0]) / len(nhit) * 100, 2), + int(np.max(nhit)), + ) + ) + + title = """ + nbolos = {}, + fp size = {} arcmin, + nhit boost = {} + """.format(nfid_bolometer, fp_size, boost) if flatsky: if fullsky: # projecting full sky with the given resolution onto a squared array. - A_sky_deg2 = 360.**2/np.pi - N_pix = 60./reso * np.sqrt(A_sky_deg2) # number of pixels on a side - N_pix = np.ceil(N_pix) # rounding up, integer number required - ################## + A_sky_deg2 = 360.0 ** 2 / np.pi + N_pix = ( + 60.0 / reso * np.sqrt(A_sky_deg2) + ) # number of pixels on a side + N_pix = np.ceil(N_pix) # rounding up, integer number required + ######### # for full sky: - lonra = [-180,180] - latra = [-90,90] - ################## - hpcp = hp.projector.CartesianProj(xsize=N_pix,ysize=N_pix, lonra=lonra, latra=latra) - f = lambda x,y,z: hp.pixelfunc.vec2pix(nside,x,y,z,nest=nest) - flat_hits = np.flipud(hpcp.projmap(nhit, f)) # flipping array just for aesthetics + lonra = [-180, 180] + latra = [-90, 90] + ######### + hpcp = hp.projector.CartesianProj( + xsize=N_pix, ysize=N_pix, lonra=lonra, latra=latra + ) + flat_hits = np.flipud( + hpcp.projmap( + nhit, + lambda x, y, z: + hp.pixelfunc.vec2pix(nside, x, y, z, nest=nest) + ) + ) # flipping array just for aesthetics # plot - pl.imshow(flat_hits,cmap=pl.cm.viridis) + pl.imshow(flat_hits, cmap=pl.cm.viridis) pl.colorbar() - pl.title('npix = {}, '.format(int(N_pix)) + 'nbolos = {}, '.format(nfid_bolometer) + 'fp size = {} arcmin, '.format(fp_size) + 'nhit boost = {}'.format(boost)) - nhit = flat_hits # will return a flat hits map + title_ = """ + npix = {}, + nbolos = {}, + fp size = {} arcmin, + nhit boost = {} + """.format(int(N_pix), nfid_bolometer, fp_size, boost) + pl.title(title_) + nhit = flat_hits # will return a flat hits map else: nhit[nhit == 0] = hp.UNSEEN - nhit = hp.gnomview(nhit, rot=rot, reso=reso, xsize=xsize, - cmap=pl.cm.viridis, - title='nbolos = {}, '.format(nfid_bolometer) + - 'fp size = {} arcmin, '.format(fp_size) + - 'nhit boost = {}'.format(boost), return_projected_map = True) # will return a flat hits map + nhit = hp.gnomview( + nhit, + rot=rot, + reso=reso, + xsize=xsize, + cmap=pl.cm.viridis, + title=title, + return_projected_map=True, + ) # will return a flat hits map hp.graticule(verbose=self.verbose) else: if fullsky: nhit[nhit == 0] = hp.UNSEEN - hp.mollview(nhit, rot=rot, cmap=pl.cm.viridis, - title='nbolos = {}, '.format(nfid_bolometer) + - 'fp size = {} arcmin, '.format(fp_size) + - 'nhit boost = {}'.format(boost)) + hp.mollview( + nhit, + rot=rot, + cmap=pl.cm.viridis, + title=title + ) hp.graticule(verbose=self.verbose) else: nhit[nhit == 0] = hp.UNSEEN - hp.gnomview(nhit, rot=rot, reso=reso, xsize=xsize, - cmap=pl.cm.viridis, - title='nbolos = {}, '.format(nfid_bolometer) + - 'fp size = {} arcmin, '.format(fp_size) + - 'nhit boost = {}'.format(boost)) + hp.gnomview( + nhit, + rot=rot, + reso=reso, + xsize=xsize, + cmap=pl.cm.viridis, + title=title, + ) hp.graticule(verbose=self.verbose) pl.show() @@ -741,6 +916,7 @@ def _update(self, name, value): """ setattr(self, name, value) + def convolve_focalplane(bore_nhits, nbolos, fp_radius_amin, boost): """ Given a hit count map, @@ -770,11 +946,11 @@ def convolve_focalplane(bore_nhits, nbolos, fp_radius_amin, boost): ---------- >>> bore_nhits = np.zeros(hp.nside2npix(128)) - ## Put a patch in the center + # Put a patch in the center >>> bore_nhits[hp.query_disc(128, hp.ang2vec(np.pi/2, 0.), ... radius=10*np.pi/180.)] = 1. - ## Increase the number of detector (x100) and make a boost (x10) + # Increase the number of detector (x100) and make a boost (x10) >>> conv_bore_nhits = convolve_focalplane(bore_nhits, nbolos=100, ... fp_radius_amin=180, boost=10) >>> print(round(np.max(bore_nhits), 2), round(np.max(conv_bore_nhits), 2)) @@ -786,25 +962,24 @@ def convolve_focalplane(bore_nhits, nbolos, fp_radius_amin, boost): # Resolution of our healpix map nside = hp.npix2nside(focalplane_nhits.shape[0]) resol_amin = hp.nside2resol(nside, arcmin=True) - fp_rad_bins = int(fp_radius_amin * 2. / resol_amin) + fp_rad_bins = int(fp_radius_amin * 2.0 / resol_amin) fp_diam_bins = (fp_rad_bins * 2) + 1 # Build the focal plane model and a list of offsets (x_fp, y_fp) = np.array( - np.unravel_index( - range(0, fp_diam_bins**2), - (fp_diam_bins, fp_diam_bins))).reshape( - 2, fp_diam_bins, fp_diam_bins) - (fp_rad_bins) - fp_map = ((x_fp**2 + y_fp**2) < (fp_rad_bins)**2) + np.unravel_index(range(0, fp_diam_bins ** 2), (fp_diam_bins, fp_diam_bins)) + ).reshape(2, fp_diam_bins, fp_diam_bins) - (fp_rad_bins) + fp_map = (x_fp ** 2 + y_fp ** 2) < (fp_rad_bins) ** 2 bolo_per_pix = nbolos / float(np.sum(fp_map)) - dRA = np.ndarray.flatten( - (x_fp[fp_map].astype(float) * fp_radius_amin) / ( - fp_rad_bins * 60. * (180. / (np.pi)))) - dDec = np.ndarray.flatten( - (y_fp[fp_map].astype(float) * fp_radius_amin) / ( - fp_rad_bins * 60. * (180. / (np.pi)))) + fp_amin_bins = fp_rad_bins * 60.0 * (180.0 / (np.pi)) + dra = np.ndarray.flatten( + (x_fp[fp_map].astype(float) * fp_radius_amin) / fp_amin_bins + ) + ddec = np.ndarray.flatten( + (y_fp[fp_map].astype(float) * fp_radius_amin) / fp_amin_bins + ) pixels_global = np.array(np.where(bore_nhits != 0)[0], dtype=int) for n in pixels_global: @@ -812,21 +987,23 @@ def convolve_focalplane(bore_nhits, nbolos, fp_radius_amin, boost): # Compute pointing offsets (theta_bore, phi_bore) = hp.pix2ang(nside, n) - phi = phi_bore + dRA * np.sin(theta_bore) - theta = theta_bore + dDec + phi = phi_bore + dra * np.sin(theta_bore) + theta = theta_bore + ddec pixels = hp.ang2pix(nside, theta, phi) npix_loc = len(pixels) - ## Necessary because the values in pixels aren't necessarily unique - ## This is a poor design choice and should probably be fixed + # Necessary because the values in pixels aren't necessarily unique + # This is a poor design choice and should probably be fixed scanning_strategy_f.convolve_focalplane_f( - bore_nhits[n], focalplane_nhits, pixels, - bolo_per_pix, boost, npix_loc) + bore_nhits[n], focalplane_nhits, pixels, bolo_per_pix, boost, npix_loc + ) return focalplane_nhits -## Here are a bunch of routines to handle dates... + +# Here are a bunch of routines to handle dates... + def date_to_mjd(date): """ @@ -851,7 +1028,7 @@ def date_to_mjd(date): Examples ---------- >>> e = ephem.Observer() - >>> e.date = 0.0 ## 1899 December 31 12:00 UT + >>> e.date = 0.0 # 1899 December 31 12:00 UT >>> mjd = date_to_mjd(e.date) >>> print('DATE={} ->'.format(round(e.date, 2)), ... 'MJD={}'.format(round(mjd, 2))) @@ -859,6 +1036,7 @@ def date_to_mjd(date): """ return greg_to_mjd(date_to_greg(date)) + def date_to_greg(date): """ Convert date in ephem.date format to gregorian date. @@ -882,7 +1060,7 @@ def date_to_greg(date): Examples ---------- >>> e = ephem.Observer() - >>> e.date = 0.0 ## 1899 December 31 12:00 UT + >>> e.date = 0.0 # 1899 December 31 12:00 UT >>> greg = date_to_greg(e.date) >>> print('DATE={} ->'.format(round(e.date, 2)), ... 'GREG={}'.format(greg)) @@ -890,13 +1068,11 @@ def date_to_greg(date): """ date_ = str(date) date_ = str(date.datetime()) - greg = date_.split('.')[0].replace('-', - '').replace(':', - '').replace(' ', - '_') + greg = date_.split(".")[0].replace("-", "").replace(":", "").replace(" ", "_") return greg + def greg_to_mjd(greg): """ Convert gregorian date into MJD. @@ -931,6 +1107,7 @@ def greg_to_mjd(greg): return mjd + def mjd_to_greg(mjd): """ Convert MJD into gregorian date. @@ -955,19 +1132,21 @@ def mjd_to_greg(mjd): year, month, day, fracday, baddate = slalib.sla_djcl(mjd) if baddate: - raise ValueError(BadMJD) + raise ValueError("Bad MJD") sign, (hour, minute, second, frac) = slalib.sla_dd2tf(2, fracday) - s = '{:4d}{:2d}{:2d}_{:2d}{:2d}{:2d}'.format( - year, month, day, hour, minute, second) - s = s.replace(' ', '0') + s = "{:4d}{:2d}{:2d}_{:2d}{:2d}{:2d}".format( + year, month, day, hour, minute, second + ) + s = s.replace(" ", "0") return s if __name__ == "__main__": import doctest + if np.__version__ >= "1.14.0": np.set_printoptions(legacy="1.13") doctest.testmod() From 7bb98b51862b92c53d93b46c4e7c32737f20122d Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 10:45:58 +0100 Subject: [PATCH 07/15] Clean input sky --- s4cmb/input_sky.py | 273 +++++++++++++++++++++++++++------------------ 1 file changed, 165 insertions(+), 108 deletions(-) diff --git a/s4cmb/input_sky.py b/s4cmb/input_sky.py index 636ebc9..6890d45 100755 --- a/s4cmb/input_sky.py +++ b/s4cmb/input_sky.py @@ -16,18 +16,29 @@ class HealpixFitsMap to handle it easily. import healpy as hp import numpy as np -from astropy.io import fits as pyfits -from s4cmb.config_s4cmb import compare_version_number from s4cmb.tools import alm2map_spin_der1 -class HealpixFitsMap(): + +class HealpixFitsMap: """ Class to handle fits file containing healpix maps """ - def __init__(self, input_filename, - do_pol=True, verbose=False, fwhm_in=0.0, fwhm_in2=None, - nside_in=16, lmax=None, map_seed=53543, no_ileak=False, - no_quleak=False, compute_derivatives=None,derivatives_type='T', - ext_map_gal=False): + + def __init__( + self, + input_filename, + do_pol=True, + verbose=False, + fwhm_in=0.0, + fwhm_in2=None, + nside_in=16, + lmax=None, + map_seed=53543, + no_ileak=False, + no_quleak=False, + compute_derivatives=None, + derivatives_type="T", + ext_map_gal=False, + ): """ Parameters @@ -113,28 +124,30 @@ def __init__(self, input_filename, print("Reading sky maps from alms file...") self.load_healpix_fits_map_from_alms() fromalms = True - elif self.input_filename[-4:] == '.dat': + elif self.input_filename[-4:] == ".dat": if self.verbose: print("Creating sky maps from cl file...") self.create_healpix_fits_map() - elif self.input_filename[-5:] == '.fits': + elif self.input_filename[-5:] == ".fits": if self.verbose: print("Reading sky maps from fits file...") self.load_healpix_fits_map() else: - raise IOError("Input file not understood! Should be either a " + - "fits file containing the sky maps " + - "(data will be loaded), or a CAMB lensed cl file " + - "(.dat) containing lensed power spectra with " + - "order ell, TT, EE, BB, TE " + - "(maps will be created on-the-fly).") + raise IOError( + """Input file not understood! Should be either a + fits file containing the sky maps + (data will be loaded), or a CAMB lensed cl file + (.dat) containing lensed power spectra with + order ell, TT, EE, BB, TE + (maps will be created on-the-fly).""" + ) self.set_leakage_to_zero() if self.compute_derivatives: - if 'T' in self.derivatives_type: + if "T" in self.derivatives_type: self.compute_intensity_derivatives(fromalm=fromalms) - if 'P' in self.derivatives_type: + if "P" in self.derivatives_type: self.compute_pol_derivatives(fromalm=fromalms) def load_healpix_fits_map(self, force=False): @@ -170,10 +183,12 @@ def load_healpix_fits_map(self, force=False): if self.I is None or force: if self.do_pol: self.I, self.Q, self.U = hp.read_map( - self.input_filename, (0, 1, 2), verbose=self.verbose) + self.input_filename, (0, 1, 2), verbose=self.verbose + ) else: self.I = hp.read_map( - self.input_filename, field=0, verbose=self.verbose) + self.input_filename, field=0, verbose=self.verbose + ) self.nside = hp.npix2nside(len(self.I)) else: print("External data already present in memory") @@ -228,17 +243,23 @@ def load_healpix_fits_map_from_alms(self, force=False): [tlm, elm, blm], nside=self.nside_in, pixwin=False, - fwhm=self.fwhm_in / 60. * np.pi / 180., - sigma=None, pol=True, inplace=False, - verbose=self.verbose) + fwhm=self.fwhm_in / 60.0 * np.pi / 180.0, + sigma=None, + pol=True, + inplace=False, + verbose=self.verbose, + ) if self.fwhm_in2 is not None: self.I2, self.Q2, self.U2 = hp.alm2map( [tlm, elm, blm], nside=self.nside_in, pixwin=False, - fwhm=self.fwhm_in2 / 60. * np.pi / 180., - sigma=None, pol=True, inplace=False, - verbose=self.verbose) + fwhm=self.fwhm_in2 / 60.0 * np.pi / 180.0, + sigma=None, + pol=True, + inplace=False, + verbose=self.verbose, + ) else: tlm = hp.read_alm(self.input_filename[0]) @@ -246,17 +267,23 @@ def load_healpix_fits_map_from_alms(self, force=False): tlm, nside=self.nside_in, pixwin=False, - fwhm=self.fwhm_in / 60. * np.pi / 180., - sigma=None, pol=False, inplace=False, - verbose=self.verbose) + fwhm=self.fwhm_in / 60.0 * np.pi / 180.0, + sigma=None, + pol=False, + inplace=False, + verbose=self.verbose, + ) if self.fwhm_in2 is not None: self.I2 = hp.alm2map( tlm, nside=self.nside_in, pixwin=False, - fwhm=self.fwhm_in2 / 60. * np.pi / 180., - sigma=None, pol=False, inplace=False, - verbose=self.verbose) + fwhm=self.fwhm_in2 / 60.0 * np.pi / 180.0, + sigma=None, + pol=False, + inplace=False, + verbose=self.verbose, + ) self.nside = hp.npix2nside(len(self.I)) else: @@ -304,28 +331,32 @@ def create_healpix_fits_map(self, force=False): nside=self.nside_in, FWHM=self.fwhm_in, seed=self.map_seed, - lmax=self.lmax) + lmax=self.lmax, + ) if self.fwhm_in2 is not None: self.I2, self.Q2, self.U2 = create_sky_map( self.input_filename, nside=self.nside_in, FWHM=self.fwhm_in2, seed=self.map_seed, - lmax=self.lmax) + lmax=self.lmax, + ) else: self.I = create_sky_map( self.input_filename, nside=self.nside_in, FWHM=self.fwhm_in, seed=self.map_seed, - lmax=self.lmax) + lmax=self.lmax, + ) if self.fwhm_in2 is not None: self.I2 = create_sky_map( self.input_filename, nside=self.nside_in, FWHM=self.fwhm_in2, seed=self.map_seed, - lmax=self.lmax) + lmax=self.lmax, + ) self.nside = hp.npix2nside(len(self.I)) else: print("External data already present in memory") @@ -356,13 +387,13 @@ def set_leakage_to_zero(self): >>> print(hpmap.I, hpmap.I2) [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] """ - ## Set temperature to zero to avoid I->QU leakage + # Set temperature to zero to avoid I->QU leakage if self.no_ileak: self.I[:] = 0.0 if self.I2 is not None: self.I2[:] = 0.0 - ## Set polarisation to zero to avoid QU leakage + # Set polarisation to zero to avoid QU leakage if self.no_quleak: if self.Q is not None: self.Q[:] = 0.0 @@ -400,17 +431,18 @@ def compute_intensity_derivatives(self, fromalm=False): else: alm = hp.map2alm(self.I, self.lmax) - lmax=hp.Alm.getlmax(alm.size) - if 'T1' in self.derivatives_type: + # lmax = hp.Alm.getlmax(alm.size) + if "T1" in self.derivatives_type: junk, self.dIdt, self.dIdp = hp.alm2map_der1( - alm, self.nside_in, self.lmax) + alm, self.nside_in, self.lmax + ) else: # computes first and second derivative as derivatives of spin-1 # transform of a scalar field with _1Elm=sqrt(l(l+1))Ilm _1Blm=0 - l=np.arange(self.lmax+1) - grad=np.sqrt(l*(l+1)) + l = np.arange(self.lmax + 1) + grad = np.sqrt(l * (l + 1)) curl = np.zeros_like(alm) - dervs = alm2map_spin_der1([hp.almxfl(alm,grad),curl],self.nside_in,1) + dervs = alm2map_spin_der1([hp.almxfl(alm, grad), curl], self.nside_in, 1) self.dIdt = dervs[0][0] self.dIdp = dervs[0][1] self.d2Id2t = dervs[1][0] @@ -442,23 +474,27 @@ def compute_pol_derivatives(self, fromalm=False): Elm = hp.read_alm(self.input_filename[1]) Blm = hp.read_alm(self.input_filename[2]) else: - alm = hp.map2alm([self.I,self.Q,self.U], self.lmax) - Elm=alm[1] - Blm=alm[2] - lmax=hp.Alm.getlmax(Elm.size) - if 'P1' in self.derivatives_type: - out = alm2map_spin_der1([Elm,Blm], self.nside_in, 2) - self.dQdt =out[1][0] - self.dUdt =out[1][1] - self.dQdp =out[2][0] - self.dUdp =out[2][1] + alm = hp.map2alm([self.I, self.Q, self.U], self.lmax) + Elm = alm[1] + Blm = alm[2] + # lmax = hp.Alm.getlmax(Elm.size) + if "P1" in self.derivatives_type: + out = alm2map_spin_der1([Elm, Blm], self.nside_in, 2) + self.dQdt = out[1][0] + self.dUdt = out[1][1] + self.dQdp = out[2][0] + self.dUdp = out[2][1] else: - warnings.warn('Computation of second order polarization derivatives not implemented yet. Set to 0.') - out = alm2map_spin_der1([Elm,Blm], self.nside_in, 2) - self.dQdt =out[1][0] - self.dUdt =out[1][1] - self.dQdp =out[2][0] - self.dUdp =out[2][1] + warnings.warn(""" + Computation of second order polarization derivatives not + implemented yet. Set to 0. + """) + + out = alm2map_spin_der1([Elm, Blm], self.nside_in, 2) + self.dQdt = out[1][0] + self.dUdt = out[1][1] + self.dQdp = out[2][0] + self.dUdp = out[2][1] self.d2Qd2t = np.zeros_like(self.dQdt) self.d2Qd2p = np.zeros_like(self.dQdt) self.d2Qdpdt = np.zeros_like(self.dQdt) @@ -466,6 +502,7 @@ def compute_pol_derivatives(self, fromalm=False): self.d2Ud2p = np.zeros_like(self.dQdt) self.d2Udpdt = np.zeros_like(self.dQdt) + def add_hierarch(lis): """ Convert in correct format for fits header. @@ -488,11 +525,12 @@ def add_hierarch(lis): """ for i, item in enumerate(lis): if len(item) == 3: - lis[i] = ('HIERARCH ' + item[0], item[1], item[2]) + lis[i] = ("HIERARCH " + item[0], item[1], item[2]) else: - lis[i] = ('HIERARCH ' + item[0], item[1]) + lis[i] = ("HIERARCH " + item[0], item[1]) return lis + def get_obspix(xmin, xmax, ymin, ymax, nside): """ Given RA/Dec boundaries, return the observed pixels in the healpix scheme. @@ -523,21 +561,21 @@ def get_obspix(xmin, xmax, ymin, ymax, nside): 18, 19, 20, 21, 26, 27, 28, 29, 30, 34, 35, 36, 37, 42, 43, 44]) """ - theta_min = np.pi / 2. - ymax - theta_max = np.pi / 2. - ymin - fpix, lpix = hp.ang2pix(nside, [theta_min, theta_max], [0., 2.*np.pi]) + theta_min = np.pi / 2.0 - ymax + theta_max = np.pi / 2.0 - ymin + fpix, lpix = hp.ang2pix(nside, [theta_min, theta_max], [0.0, 2.0 * np.pi]) pixs = np.arange(fpix, lpix + 1, dtype=np.int) theta, phi = hp.pix2ang(nside, pixs) if xmin < 0: - phi[phi > np.pi] = (phi[phi > np.pi] - 2 * np.pi) - good = (theta >= theta_min) * (theta <= theta_max) * \ - (phi <= xmax) * (phi >= xmin) + phi[phi > np.pi] = phi[phi > np.pi] - 2 * np.pi + good = (theta >= theta_min) * (theta <= theta_max) * (phi <= xmax) * (phi >= xmin) obspix = pixs[good] obspix.sort() return obspix + def LamCyl(ra, dec): """ Referred to as cylindrical equal-area in the USGS report, assuming @@ -545,18 +583,22 @@ def LamCyl(ra, dec): """ return ra, np.sin(dec) + def SFL(ra, dec): - '''SFL stands for Sanson-Flamsteed. In the USGS report this is + """SFL stands for Sanson-Flamsteed. In the USGS report this is referred to as the Sinusoidal Projection. It is equal-area. Parallels are equally spaced straight lines. Scale is true along central meridian - and all paralles.''' + and all paralles.""" return ra * np.cos(dec), dec -def deSFL(x,y): - return x/np.cos(y),y -def deLamCyl(x,y): - return x,np.arcsin(y) +def deSFL(x, y): + return x / np.cos(y), y + + +def deLamCyl(x, y): + return x, np.arcsin(y) + def create_sky_map(cl_fn, nside=16, FWHM=0.0, seed=548397, lmax=None): """ @@ -587,27 +629,41 @@ def create_sky_map(cl_fn, nside=16, FWHM=0.0, seed=548397, lmax=None): 80.8613084 ] """ if lmax is None: - lmax = 2*nside + lmax = 2 * nside ell, TT, EE, BB, TE = np.loadtxt(cl_fn).T - ## Take out the normalisation... - llp = ell * (ell + 1.) / (2 * np.pi) + # Take out the normalisation... + llp = ell * (ell + 1.0) / (2 * np.pi) - ## Arcmin to rad - FWHM_rad = FWHM / 60. * np.pi / 180. + # Arcmin to rad + FWHM_rad = FWHM / 60.0 * np.pi / 180.0 np.random.seed(seed) I, Q, U = hp.synfast( - [TT / llp, EE / llp, BB / llp, TE / llp], nside, - lmax=lmax, mmax=None, alm=False, - pol=True, pixwin=False, - fwhm=FWHM_rad, sigma=None, new=True, - verbose=False) + [TT / llp, EE / llp, BB / llp, TE / llp], + nside, + lmax=lmax, + mmax=None, + alm=False, + pol=True, + pixwin=False, + fwhm=FWHM_rad, + sigma=None, + new=True, + verbose=False, + ) return I, Q, U -def write_healpix_cmbmap(output_filename, data, fits_IDL=False, - coord=None, colnames=['I', 'Q', 'U'], partial=True, - nest=False): + +def write_healpix_cmbmap( + output_filename, + data, + fits_IDL=False, + coord=None, + colnames=["I", "Q", "U"], + partial=True, + nest=False, +): """ Write healpix fits map in full sky mode or partial sky, Input data have to be a list with n fields to be written. @@ -640,26 +696,25 @@ def write_healpix_cmbmap(output_filename, data, fits_IDL=False, >>> write_healpix_cmbmap('myfits_to_test_.fits', ... data=[I, Q, U], colnames=colnames) """ - ## Write the header + # Write the header extra_header = [] for c in colnames: - extra_header.append(('column_names', c)) + extra_header.append(("column_names", c)) extra_header = add_hierarch(extra_header) - ## Need to introduce this workaround because the last - ## version of healpy introduced non-backward compatibility... - try: - hp.write_map( - output_filename, data, fits_IDL=fits_IDL, - coord=coord, column_names=None, partial=partial, - extra_header=extra_header, overwrite=True) - except: - hp.write_map( - output_filename, data, fits_IDL=fits_IDL, - coord=coord, column_names=None, partial=partial, - extra_header=extra_header) - -def write_dummy_map(filename='myfits_to_test_.fits', nside=16): + hp.write_map( + output_filename, + data, + fits_IDL=fits_IDL, + coord=coord, + column_names=None, + partial=partial, + extra_header=extra_header, + overwrite=True, + ) + + +def write_dummy_map(filename="myfits_to_test_.fits", nside=16): """ Write dummy file on disk for test purposes. @@ -676,10 +731,11 @@ def write_dummy_map(filename='myfits_to_test_.fits', nside=16): """ nside = 16 I, Q, U = np.random.rand(3, hp.nside2npix(nside)) - colnames = ['I', 'Q', 'U'] + colnames = ["I", "Q", "U"] write_healpix_cmbmap(filename, data=[I, Q, U], colnames=colnames) -def remove_test_data(has_id='_to_test_', silent=True): + +def remove_test_data(has_id="_to_test_", silent=True): """ Remove data with name containing the `has_id`. @@ -695,16 +751,17 @@ def remove_test_data(has_id='_to_test_', silent=True): >>> remove_test_data(has_id='_to_erase_', silent=False) Removing files: ['file_to_erase_.txt'] """ - fns = glob.glob('*' + has_id + '*') + fns = glob.glob("*" + has_id + "*") if not silent: - print('Removing files: ', fns) + print("Removing files: ", fns) for fn in fns: os.remove(fn) if __name__ == "__main__": import doctest + if np.__version__ >= "1.14.0": np.set_printoptions(legacy="1.13") doctest.testmod() - remove_test_data(has_id='_to_test_', silent=True) + remove_test_data(has_id="_to_test_", silent=True) From 1242821d2fc5205634130f68b90cab747793e25d Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 10:49:58 +0100 Subject: [PATCH 08/15] CLean instrument --- s4cmb/instrument.py | 505 +++++++++++++++++++++++++++----------------- 1 file changed, 308 insertions(+), 197 deletions(-) diff --git a/s4cmb/instrument.py b/s4cmb/instrument.py index 11a52a6..d14bc69 100755 --- a/s4cmb/instrument.py +++ b/s4cmb/instrument.py @@ -10,15 +10,13 @@ """ from __future__ import division, absolute_import, print_function -import os import copy -import glob -import datetime import numpy as np -def coordinates_on_grid(pix_size=None, row_size=None, - nx=None, nx2=None, - max_points=None): + +def coordinates_on_grid( + pix_size=None, row_size=None, nx=None, nx2=None, max_points=None +): """ Return the x and y coordinates of points on a grid. The grid is centered on (0, 0). @@ -92,18 +90,23 @@ def coordinates_on_grid(pix_size=None, row_size=None, a pixel (pix_size), or the size of a row (row_size). """ if (nx is None and nx2 is None) or (nx is not None and nx2 is not None): - raise AssertionError('You should specify either the ' + - 'number of pixel per row column (nx), or ' + - 'the total number of point in the grid (nx2).\n') - - if (pix_size is None and row_size is None) or \ - (pix_size is not None and row_size is not None): - raise AssertionError('You should specify either the ' + - 'size of a pixel (pix_size), ' + - 'or the size of a row (row_size).\n') + raise AssertionError( + """You should specify either the + number of pixel per row column (nx), or + the total number of point in the grid (nx2).\n""" + ) + + if (pix_size is None and row_size is None) or ( + pix_size is not None and row_size is not None + ): + raise AssertionError( + """You should specify either the + size of a pixel (pix_size), + or the size of a row (row_size).\n""" + ) if nx2 is not None: - ## Look for the closest number with square root being an integer + # Look for the closest number with square root being an integer nx2_tmp = copy.copy(nx2) while True: nx = np.sqrt(nx2_tmp) @@ -113,7 +116,7 @@ def coordinates_on_grid(pix_size=None, row_size=None, else: nx2_tmp += 1 else: - nx2 = nx**2 + nx2 = nx ** 2 if max_points is None: max_points = nx2 @@ -124,15 +127,14 @@ def coordinates_on_grid(pix_size=None, row_size=None, row_size = pix_size * nx ix1 = np.arange(nx) - xs1 = (ix1 - (nx - 1.) / 2.) * pix_size + xs1 = (ix1 - (nx - 1.0) / 2.0) * pix_size x2, y2 = np.meshgrid(xs1, xs1) - coordinates = np.array( - (x2.flatten()[:max_points], - y2.flatten()[:max_points])) + coordinates = np.array((x2.flatten()[:max_points], y2.flatten()[:max_points])) return coordinates + def convert_pair_to_bolometer_position(xcoord_pairs, ycoord_pairs): """ Return the position of bolometers given the position of pairs. @@ -167,16 +169,24 @@ def convert_pair_to_bolometer_position(xcoord_pairs, ycoord_pairs): [-15. -15. -15. -15. 15. 15. 15. 15.] """ nbolometer = 2 * len(xcoord_pairs) - xcoord_bolometers = np.dstack( - (xcoord_pairs, xcoord_pairs)).reshape((1, nbolometer))[0] - ycoord_bolometers = np.dstack( - (ycoord_pairs, ycoord_pairs)).reshape((1, nbolometer))[0] + xcoord_bolometers = np.dstack((xcoord_pairs, xcoord_pairs)).reshape( + (1, nbolometer) + )[0] + ycoord_bolometers = np.dstack((ycoord_pairs, ycoord_pairs)).reshape( + (1, nbolometer) + )[0] return xcoord_bolometers, ycoord_bolometers -def show_focal_plane(bolo_xcoord, bolo_ycoord, bolo_polangle=None, - fn_out='plot_hardware_map_test.png', - save_on_disk=True, display=False): + +def show_focal_plane( + bolo_xcoord, + bolo_ycoord, + bolo_polangle=None, + fn_out="plot_hardware_map_test.png", + save_on_disk=True, + display=False, +): """ Show the focal plane of the instrument, split in two panels: top and bottom bolometers @@ -206,8 +216,10 @@ def show_focal_plane(bolo_xcoord, bolo_ycoord, bolo_polangle=None, """ if not display: import matplotlib as mpl - mpl.use('Agg') + + mpl.use("Agg") import matplotlib.pyplot as pl + pl.ioff() else: import matplotlib.pyplot as pl @@ -216,25 +228,49 @@ def show_focal_plane(bolo_xcoord, bolo_ycoord, bolo_polangle=None, bolo_polangle = np.ones_like(bolo_xcoord) fig, ax = pl.subplots(1, 2, figsize=(10, 5)) - ## Top pixel - ax[0].scatter(bolo_xcoord[::2], bolo_ycoord[::2], - c=bolo_polangle[::2], alpha=1, s=30, cmap=pl.cm.jet) - ax[0].scatter(bolo_xcoord[::2], bolo_ycoord[::2], - c='black', s=30, marker='|', - label='Top pixel', alpha=0.6) - ax[0].set_ylabel('y position (cm)') - ax[0].set_xlabel('x position (cm)') - ax[0].set_title('Top pixels') - - ## Bottom pixel - ax[1].scatter(bolo_xcoord[1::2], bolo_ycoord[1::2], - c=bolo_polangle[1::2], alpha=1, s=30, cmap=pl.cm.jet) - ax[1].scatter(bolo_xcoord[1::2], bolo_ycoord[1::2], - c='black', s=30, marker='_', - label='Bottom pixel', alpha=0.6) - ax[1].set_ylabel('y position (cm)') - ax[1].set_xlabel('x position (cm)') - ax[1].set_title('Bottom pixels') + # Top pixel + ax[0].scatter( + bolo_xcoord[::2], + bolo_ycoord[::2], + c=bolo_polangle[::2], + alpha=1, + s=30, + cmap=pl.cm.jet, + ) + ax[0].scatter( + bolo_xcoord[::2], + bolo_ycoord[::2], + c="black", + s=30, + marker="|", + label="Top pixel", + alpha=0.6, + ) + ax[0].set_ylabel("y position (cm)") + ax[0].set_xlabel("x position (cm)") + ax[0].set_title("Top pixels") + + # Bottom pixel + ax[1].scatter( + bolo_xcoord[1::2], + bolo_ycoord[1::2], + c=bolo_polangle[1::2], + alpha=1, + s=30, + cmap=pl.cm.jet, + ) + ax[1].scatter( + bolo_xcoord[1::2], + bolo_ycoord[1::2], + c="black", + s=30, + marker="_", + label="Bottom pixel", + alpha=0.6, + ) + ax[1].set_ylabel("y position (cm)") + ax[1].set_xlabel("x position (cm)") + ax[1].set_title("Bottom pixels") if save_on_disk: pl.savefig(fn_out) @@ -242,6 +278,7 @@ def show_focal_plane(bolo_xcoord, bolo_ycoord, bolo_polangle=None, if display: pl.show() + def convert_cm_to_rad(xcm, ycm, conversion): """ Convert positions in cm of the pairs of bolometers or bolometers @@ -270,7 +307,7 @@ def convert_cm_to_rad(xcm, ycm, conversion): Focal plane of 60 cm diameter and mirror giving a 3 deg projection on the sky by default >>> fp = FocalPlane(verbose=False) - >>> projected_fp_size = 3. ## degrees + >>> projected_fp_size = 3. # degrees >>> xcm, ycm = coordinates_on_grid( ... row_size=fp.fp_size, nx2=fp.npair) >>> print(convert_cm_to_rad(xcm, ycm, @@ -282,6 +319,7 @@ def convert_cm_to_rad(xcm, ycm, conversion): """ return np.array(xcm) * conversion, np.array(ycm) * conversion + def construct_beammap(beamprm, ct, cb, nx, pix_size): """ Construct the pixel beam maps @@ -333,21 +371,32 @@ def construct_beammap(beamprm, ct, cb, nx, pix_size): xy2f = coordinates_on_grid(pix_size=pix_size, nx=nx) - tmap = gauss2d(xy2f, tx, ty, - beamprm.Amp[ct], beamprm.sig_1[ct], - beamprm.sig_2[ct], - beamprm.ellip_ang[ct]).reshape((nx, nx)) - - bmap = gauss2d(xy2f, bx, by, - beamprm.Amp[cb], beamprm.sig_1[cb], - beamprm.sig_2[cb], - beamprm.ellip_ang[cb]).reshape((nx, nx)) + tmap = gauss2d( + xy2f, + tx, + ty, + beamprm.Amp[ct], + beamprm.sig_1[ct], + beamprm.sig_2[ct], + beamprm.ellip_ang[ct], + ).reshape((nx, nx)) + + bmap = gauss2d( + xy2f, + bx, + by, + beamprm.Amp[cb], + beamprm.sig_1[cb], + beamprm.sig_2[cb], + beamprm.ellip_ang[cb], + ).reshape((nx, nx)) summap = 0.5 * (tmap + bmap) diffmap = 0.5 * (tmap - bmap) return summap, diffmap + def gauss2d(xy, x_0, y_0, Amp, sig_xp, sig_yp, psi): """ 2D Gaussian model for beams. @@ -401,11 +450,10 @@ def gauss2d(xy, x_0, y_0, Amp, sig_xp, sig_yp, psi): psi2 = -psi * np.pi / 180.0 # x/y coordinates make an angle psi with the xp/yp input coordinates - R = np.array( - [[np.cos(psi2), -np.sin(psi2)], [np.sin(psi2), np.cos(psi2)]]) + R = np.array([[np.cos(psi2), -np.sin(psi2)], [np.sin(psi2), np.cos(psi2)]]) p = np.dot(R, xy_1) - u = p[0, :]**2 / (2 * sig_xp**2) + p[1, :]**2 / (2 * sig_yp**2) + u = p[0, :] ** 2 / (2 * sig_xp ** 2) + p[1, :] ** 2 / (2 * sig_yp ** 2) # Hide underflow by clipping beam function at -430dB level mask = u < 100 @@ -413,15 +461,26 @@ def gauss2d(xy, x_0, y_0, Amp, sig_xp, sig_yp, psi): return z -class Hardware(): + +class Hardware: """ Class to load all the hardware and models of the instrument in once """ - def __init__(self, - ncrate=1, ndfmux_per_crate=1, nsquid_per_mux=1, - npair_per_squid=4, fp_size=60., - fwhm=3.5, beam_seed=58347, - projected_fp_size=3., - pm_name='5params', - type_hwp='CRHWP', freq_hwp=2., angle_hwp=0., verbose=False): + + def __init__( + self, + ncrate=1, + ndfmux_per_crate=1, + nsquid_per_mux=1, + npair_per_squid=4, + fp_size=60.0, + fwhm=3.5, + beam_seed=58347, + projected_fp_size=3.0, + pm_name="5params", + type_hwp="CRHWP", + freq_hwp=2.0, + angle_hwp=0.0, + verbose=False, + ): """ This class creates the data used to model the instrument: * focal plane @@ -470,19 +529,26 @@ def __init__(self, ---------- >>> instrument = Hardware() """ - self.focal_plane = FocalPlane(ncrate, ndfmux_per_crate, - nsquid_per_mux, npair_per_squid, - fp_size, verbose) - - self.beam_model = BeamModel(self.focal_plane, fwhm, beam_seed, - projected_fp_size, verbose) + self.focal_plane = FocalPlane( + ncrate, + ndfmux_per_crate, + nsquid_per_mux, + npair_per_squid, + fp_size, + verbose, + ) + + self.beam_model = BeamModel( + self.focal_plane, fwhm, beam_seed, projected_fp_size, verbose + ) self.pointing_model = PointingModel(pm_name) self.half_wave_plate = HalfWavePlate(type_hwp, freq_hwp, angle_hwp) - def make_dichroic(self, fwhm=1.8, beam_seed=58347, projected_fp_size=3., - shift_angle=45): + def make_dichroic( + self, fwhm=1.8, beam_seed=58347, projected_fp_size=3.0, shift_angle=45 + ): """ Add a layer of detectors on top of the existing ones to have dichroic detectors. The new detectors will sit on top of the @@ -515,19 +581,28 @@ def make_dichroic(self, fwhm=1.8, beam_seed=58347, projected_fp_size=3., """ self.focal_plane2 = copy.copy(self.focal_plane) - ## Shift the polarisation angles (and clip it btw 0 and 360 deg) + # Shift the polarisation angles (and clip it btw 0 and 360 deg) self.focal_plane2.bolo_polangle = ( - np.array(self.focal_plane2.bolo_polangle) + shift_angle) % 360 + np.array(self.focal_plane2.bolo_polangle) + shift_angle + ) % 360 self.beam_model2 = BeamModel( - self.focal_plane2, fwhm, beam_seed, - projected_fp_size) + self.focal_plane2, fwhm, beam_seed, projected_fp_size + ) + -class FocalPlane(): +class FocalPlane: """ Class to handle the focal plane of the instrument. """ - def __init__(self, - ncrate=1, ndfmux_per_crate=1, nsquid_per_mux=1, - npair_per_squid=4, fp_size=60., verbose=False): + + def __init__( + self, + ncrate=1, + ndfmux_per_crate=1, + nsquid_per_mux=1, + npair_per_squid=4, + fp_size=60.0, + verbose=False, + ): """ Initialise our focal plane. @@ -562,9 +637,12 @@ def __init__(self, self.nsquid_per_mux = nsquid_per_mux self.npair_per_squid = npair_per_squid - ## Total number of pairs and bolometers in the focal plane - self.npair = self.ncrate * self.ndfmux_per_crate * \ - self.nsquid_per_mux * self.npair_per_squid + # Total number of pairs and bolometers in the focal plane + self.npair = self.ncrate\ + * self.ndfmux_per_crate\ + * self.nsquid_per_mux\ + * self.npair_per_squid + self.nbolometer = self.npair * 2 self.fp_size = fp_size @@ -595,17 +673,16 @@ def make_focal_plane(self): >>> fp = FocalPlane(verbose=True) Hardware map generated... """ - ## Retrieve coordinate of the pairs inside the focal plane + # Retrieve coordinate of the pairs inside the focal plane # xcoord, ycoord = self.compute_pairs_coordinates(self.npair) - xcoord, ycoord = coordinates_on_grid(row_size=self.fp_size, - nx2=self.npair) + xcoord, ycoord = coordinates_on_grid(row_size=self.fp_size, nx2=self.npair) - ## Initialise + # Initialise self.crate_id, self.dfmux_id = [], [] self.squid_id, self.bolo_id = [], [] self.bolo_index_in_squid, self.bolo_index_in_fp = [], [] self.bolo_xcoord, self.bolo_ycoord, self.bolo_polangle = [], [], [] - ## Construct the hardware map + # Construct the hardware map max_hit = False while max_hit is False: bolo_index = 0 @@ -613,97 +690,101 @@ def make_focal_plane(self): squid_index = 0 dfmux_index = 0 for crate in range(self.ncrate): - ## CRATE - self.crate_id.append('Cr{:03}'.format(crate)) + # CRATE + self.crate_id.append("Cr{:03}".format(crate)) for dfmux in range(self.ndfmux_per_crate): - ## DFMUX - self.dfmux_id.append('Cr{:03}Df{:03}'.format( - crate, dfmux_index)) + # DFMUX + self.dfmux_id.append("Cr{:03}Df{:03}".format(crate, dfmux_index)) for squid in range(self.nsquid_per_mux): - ## SQUID - self.squid_id.append('Cr{:03}Df{:03}Sq{:03}'.format( - crate, dfmux_index, squid_index)) + # SQUID + self.squid_id.append( + "Cr{:03}Df{:03}Sq{:03}".format( + crate, dfmux_index, squid_index + ) + ) for pair in range(self.npair_per_squid): - ## BOLOMETER + # BOLOMETER - ## Split Q/U - boloQ = 2*pair - boloU = 2*pair + 1 + # Split Q/U + boloQ = 2 * pair + boloU = 2 * pair + 1 - if int(pair_index/np.sqrt(self.npair)) % 2 == 0: - shift = 0. - shiftinv = 45. + if int(pair_index / np.sqrt(self.npair)) % 2 == 0: + shift = 0.0 + shiftinv = 45.0 else: - shift = 45. - shiftinv = 0. + shift = 45.0 + shiftinv = 0.0 - ## Top pixel - ## Position of the bolometer within the SQUID + # Top pixel + # Position of the bolometer within the SQUID self.bolo_index_in_squid.append(boloQ) self.bolo_index_in_fp.append(bolo_index) self.bolo_id.append( - 'Cr{:03}Df{:03}Sq{:03}Bo{:03}t'.format( - crate, dfmux_index, squid_index, boloQ)) + "Cr{:03}Df{:03}Sq{:03}Bo{:03}t".format( + crate, dfmux_index, squid_index, boloQ + ) + ) self.bolo_xcoord.append(xcoord[pair_index]) self.bolo_ycoord.append(ycoord[pair_index]) - ## Q/U pixels + # Q/U pixels # bolo in pairs are separated by 90 deg, # and each quadran is separated by 90 deg. if bolo_index % 4 == 0: - shift_ = 45. - shiftinv_ = 0. + shift_ = 45.0 + shiftinv_ = 0.0 else: - shift_ = 0. - shiftinv_ = 45. - if xcoord[pair_index] < 0 and \ - ycoord[pair_index] >= 0: + shift_ = 0.0 + shiftinv_ = 45.0 + if xcoord[pair_index] < 0 and ycoord[pair_index] >= 0: angle = shift - elif xcoord[pair_index] >= 0 and \ - ycoord[pair_index] >= 0: - angle = 90. + shift_ - elif xcoord[pair_index] >= 0 and \ - ycoord[pair_index] <= 0: - angle = 180. + shiftinv - elif xcoord[pair_index] <= 0 and \ - ycoord[pair_index] <= 0: - angle = 270. + shiftinv_ + elif xcoord[pair_index] >= 0 and ycoord[pair_index] >= 0: + angle = 90.0 + shift_ + elif xcoord[pair_index] >= 0 and ycoord[pair_index] <= 0: + angle = 180.0 + shiftinv + elif xcoord[pair_index] <= 0 and ycoord[pair_index] <= 0: + angle = 270.0 + shiftinv_ self.bolo_polangle.append(angle) - ## Move in bolo space (2 bolos/pair) + # Move in bolo space (2 bolos/pair) bolo_index += 1 - ## Bottom pixel - ## Top pixel - ## Position of the bolometer within the SQUID + # Bottom pixel + # Top pixel + # Position of the bolometer within the SQUID self.bolo_index_in_squid.append(boloU) self.bolo_index_in_fp.append(bolo_index) self.bolo_id.append( - 'Cr{:03}Df{:03}Sq{:03}Bo{:03}b'.format( - crate, dfmux_index, squid_index, boloU)) + "Cr{:03}Df{:03}Sq{:03}Bo{:03}b".format( + crate, dfmux_index, squid_index, boloU + ) + ) self.bolo_xcoord.append(xcoord[pair_index]) self.bolo_ycoord.append(ycoord[pair_index]) - ## 90 degree difference wrt top bolometer + # 90 degree difference wrt top bolometer self.bolo_polangle.append((angle + 90) % 360) - ## Move again in bolo space (2 bolos/pair) + # Move again in bolo space (2 bolos/pair) bolo_index += 1 - ## Move in pair space + # Move in pair space pair_index += 1 - ## Close the job if you hit the maximum number of - ## bolometers or pairs. + # Close the job if you hit the maximum number of + # bolometers or pairs. try: - assert bolo_index < self.nbolometer, \ - 'Hardware map generated...' - assert pair_index < self.npair, \ - 'Hardware map generated...' + assert ( + bolo_index < self.nbolometer + ), "Hardware map generated..." + assert ( + pair_index < self.npair + ), "Hardware map generated..." except AssertionError as e: if self.verbose: print(str(e)) @@ -712,7 +793,7 @@ def make_focal_plane(self): squid_index += 1 dfmux_index += 1 - def get_indices(self, name='Cr'): + def get_indices(self, name="Cr"): """ Returns Cr(ate), Df(mux) or Sq(uid) indices. @@ -738,19 +819,26 @@ def get_indices(self, name='Cr'): [0, 0, 0, 0, 0, 0, 0, 0] """ - assert name in ['Cr', 'Sq', 'Df'], \ - ValueError("name must be in ['Cr', 'Sq', 'Df'].") + assert name in ["Cr", "Sq", "Df"], ValueError( + "name must be in ['Cr', 'Sq', 'Df']." + ) - indices = [ - int(comp.split(name)[-1][:3]) for comp in self.bolo_id] + indices = [int(comp.split(name)[-1][:3]) for comp in self.bolo_id] return indices -class BeamModel(): + +class BeamModel: """ Class to handle the beams of the detectors """ - def __init__(self, - focal_plane, fwhm=3.5, beam_seed=58347, - projected_fp_size=3., verbose=False): + + def __init__( + self, + focal_plane, + fwhm=3.5, + beam_seed=58347, + projected_fp_size=3.0, + verbose=False, + ): """ Parameters ---------- @@ -768,15 +856,15 @@ def __init__(self, verbose : boolean, optional If True, print out a number of useful comments for verboseging. """ - ## Focal plane parameters + # Focal plane parameters self.focal_plane = focal_plane - ## Beam model and mirror parameters + # Beam model and mirror parameters self.fwhm = fwhm self.beam_seed = beam_seed self.projected_fp_size = projected_fp_size - ## Paths and names + # Paths and names self.verbose = verbose self.beamprm = self.generate_beam_parameters() @@ -801,43 +889,56 @@ def generate_beam_parameters(self): -0.01308997 -0.01308997 0.01308997 0.01308997] """ - beamprm_header = ['Amp', 'Amp_err', - 'ellip_ang', 'ellip_ang_err', - 'sig_1', 'sig_1_err', - 'sig_2', 'sig_2_err', - 'xpos', 'xpos_err', - 'ypos', 'ypos_err'] + beamprm_header = [ + "Amp", + "Amp_err", + "ellip_ang", + "ellip_ang_err", + "sig_1", + "sig_1_err", + "sig_2", + "sig_2_err", + "xpos", + "xpos_err", + "ypos", + "ypos_err", + ] for b in beamprm_header: setattr(self, b, np.zeros(self.focal_plane.nbolometer)) - ## Position of the bolometers - ## (bolometers cm -> bolometers radians) + # Position of the bolometers + # (bolometers cm -> bolometers radians) self.xpos, self.ypos = convert_cm_to_rad( - self.focal_plane.bolo_xcoord, self.focal_plane.bolo_ycoord, - conversion=self.projected_fp_size / - self.focal_plane.fp_size * np.pi / 180.) - - ## Generate Gaussian beams. - ## FWHM arcmin -> FWHM rad -> sigma rad - FWHM_rad = self.fwhm / 60. * np.pi / 180. + self.focal_plane.bolo_xcoord, + self.focal_plane.bolo_ycoord, + conversion=( + self.projected_fp_size / self.focal_plane.fp_size * np.pi / 180.0 + ), + ) + + # Generate Gaussian beams. + # FWHM arcmin -> FWHM rad -> sigma rad + FWHM_rad = self.fwhm / 60.0 * np.pi / 180.0 sigma_rad = FWHM_rad / np.sqrt(8 * np.log(2)) self.sig_1 = np.ones(self.focal_plane.nbolometer) * sigma_rad self.sig_2 = np.ones(self.focal_plane.nbolometer) * sigma_rad - ## Angle of rotation for the ellipses. - ## Between -90 and 90 degrees. + # Angle of rotation for the ellipses. + # Between -90 and 90 degrees. state = np.random.RandomState(self.beam_seed) self.ellip_ang = state.uniform(-90, 90, self.focal_plane.nbolometer) - ## Amplitude of the beams - ## Default is one. + # Amplitude of the beams + # Default is one. self.Amp = np.ones(self.focal_plane.nbolometer) -class PointingModel(): + +class PointingModel: """ Class to handle the pointing model of the telescope """ - def __init__(self, pm_name='5params'): + + def __init__(self, pm_name="5params"): """ We focus on a five-parameter pointing model (Mangum 2001) to characterize the relationship between the telescope's encoder @@ -875,30 +976,35 @@ def __init__(self, pm_name='5params'): """ self.pm_name = pm_name - if self.pm_name == '5params': + if self.pm_name == "5params": self.five_parameter_pointing_model() else: - raise ValueError('Only the five-parameter ' + - 'pointing model (Mangum 2001) is implemented ' + - 'for the moment (pm_name = 5params)') + raise ValueError( + """Only the five-parameter + pointing model (Mangum 2001) is implemented + for the moment (pm_name = 5params)""" + ) def five_parameter_pointing_model(self): """ Parameters based on Polarbear configuration. """ - self.allowed_params = 'ia ie ca an aw' + self.allowed_params = "ia ie ca an aw" - self.value_params = np.array([-10.28473073, 8.73953334, - -15.59771781, -0.50977716, 0.10858016]) + self.value_params = np.array( + [-10.28473073, 8.73953334, -15.59771781, -0.50977716, 0.10858016] + ) - ## Set this to zero for the moment + # Set this to zero for the moment self.RMS_AZ = 0.0 self.RMS_EL = 0.0 self.RMS_RESID = 0.0 -class HalfWavePlate(): + +class HalfWavePlate: """ Class to handle the Half-Wave Plate (HWP) """ - def __init__(self, type_hwp='CRHWP', freq_hwp=2., angle_hwp=0.): + + def __init__(self, type_hwp="CRHWP", freq_hwp=2.0, angle_hwp=0.0): """ This class provides routines to compute the HWP angles. This can be use later to generate timestreams. @@ -919,15 +1025,16 @@ def __init__(self, type_hwp='CRHWP', freq_hwp=2., angle_hwp=0.): self.freq_hwp = freq_hwp self.angle_hwp = angle_hwp - if self.type_hwp not in ['CRHWP', 'stepped']: + if self.type_hwp not in ["CRHWP", "stepped"]: raise ValueError("`type_hwp` has to be 'CRHWP' or 'stepped'.") - if self.type_hwp is 'stepped' and freq_hwp != 0.0: - raise AssertionError("You cannot have a stepped HWP and non-" + - "zero frequency! set freq_hwp=0.0 " + - "if you want a stepped HWP.") + if self.type_hwp == "stepped" and freq_hwp != 0.0: + raise AssertionError( + """You cannot have a stepped HWP and non-zero frequency! + set freq_hwp=0.0 if you want a stepped HWP.""" + ) - def compute_HWP_angles(self, sample_rate=1., size=1): + def compute_HWP_angles(self, sample_rate=1.0, size=1): """ Generate HWP angles which can be use later to generate timestreams. @@ -968,11 +1075,14 @@ def compute_HWP_angles(self, sample_rate=1., size=1): AssertionError: You cannot have a stepped HWP and non-zero frequency! set freq_hwp=0.0 if you want a stepped HWP. """ - angle = self.angle_hwp * np.pi / 180. + angle = self.angle_hwp * np.pi / 180.0 HWP_angles = np.array( - [angle + t * (self.freq_hwp / sample_rate) * - 2. * np.pi for t in range(size)]) + [ + angle + t * (self.freq_hwp / sample_rate) * 2.0 * np.pi + for t in range(size) + ] + ) return HWP_angles @@ -1015,6 +1125,7 @@ def update_hardware(self, new_type_hwp, new_freq_hwp, new_angle_hwp): if __name__ == "__main__": import doctest + if np.__version__ >= "1.14.0": np.set_printoptions(legacy="1.13") doctest.testmod() From 10cc7f3da37939f829f25a8163a4524750c45dd1 Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 10:51:02 +0100 Subject: [PATCH 09/15] Clean config --- s4cmb/config_s4cmb.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/s4cmb/config_s4cmb.py b/s4cmb/config_s4cmb.py index 7376018..61e6774 100755 --- a/s4cmb/config_s4cmb.py +++ b/s4cmb/config_s4cmb.py @@ -11,6 +11,7 @@ import importlib import numpy as np + def compare_version_number(version, threshold): """ Compare two version numbers. @@ -36,12 +37,12 @@ def compare_version_number(version, threshold): True """ - ## If the two versions are equal + # If the two versions are equal if version == threshold: return True - version_numbers = version.split('.') - threshold_numbers = threshold.split('.') + version_numbers = version.split(".") + threshold_numbers = threshold.split(".") for v, t in zip(version_numbers, threshold_numbers): v = int(v) t = int(t) @@ -53,6 +54,7 @@ def compare_version_number(version, threshold): return False return True + def import_string_as_module(fn_full): """ Import module from its name given as a string. @@ -76,16 +78,18 @@ def import_string_as_module(fn_full): True """ - ## Import parameters from the user parameter file - fn_short = os.path.basename(fn_full).split('.py')[0] - sys.path.insert(0, os.path.realpath( - os.path.join(os.getcwd(), os.path.dirname(fn_full)))) + # Import parameters from the user parameter file + fn_short = os.path.basename(fn_full).split(".py")[0] + sys.path.insert( + 0, os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(fn_full))) + ) params = importlib.import_module(fn_short) return params if __name__ == "__main__": import doctest + if np.__version__ >= "1.14.0": np.set_printoptions(legacy="1.13") doctest.testmod() From 7f84212815488b1907321e4ab9aecacf7ae35481 Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 11:00:47 +0100 Subject: [PATCH 10/15] Clean xpure --- s4cmb/xpure.py | 999 +++++++++++++++++++++++++++---------------------- 1 file changed, 552 insertions(+), 447 deletions(-) diff --git a/s4cmb/xpure.py b/s4cmb/xpure.py index 42c457d..0fabc82 100755 --- a/s4cmb/xpure.py +++ b/s4cmb/xpure.py @@ -11,9 +11,11 @@ import os import numpy as np +import errno from s4cmb.input_sky import write_healpix_cmbmap + def safe_mkdir(path, verbose=False): """ Create a path and catch the race condition between path exists and mkdir. @@ -47,7 +49,8 @@ def safe_mkdir(path, verbose=False): if verbose: print("Folders {} already exist. Not created.".format(path)) -def qu_weight_mineig(cc, cs, ss, epsilon=0., verbose=False): + +def qu_weight_mineig(cc, cs, ss, epsilon=0.0, verbose=False): """ Create a weight map using the smallest eigenvalue of the polarization matrix. @@ -72,10 +75,10 @@ def qu_weight_mineig(cc, cs, ss, epsilon=0., verbose=False): """ trace = cc + ss trace2 = trace * trace - det = (cc * ss - cs * cs) + det = cc * ss - cs * cs val2 = trace2 - 4 * det - valid = (val2 > 0.0) + valid = val2 > 0.0 val = np.zeros_like(val2) val[valid] = np.sqrt(val2[valid]) @@ -83,19 +86,23 @@ def qu_weight_mineig(cc, cs, ss, epsilon=0., verbose=False): lambda_minus = (trace - val) / 2 if verbose: - print('criterion is', epsilon, '< det < 1/4 (epsilon= 0. by default)') + print("criterion is", epsilon, "< det < 1/4 (epsilon= 0. by default)") - valid = (lambda_minus > (trace - np.sqrt(trace2 - 4*epsilon * trace2)) / 2) + valid = lambda_minus > (trace - np.sqrt(trace2 - 4 * epsilon * trace2)) / 2 if verbose: valid3 = [x for x in valid if x is True] - print('number of pixels kept:', len(valid3), '/', np.sum(trace > 0)) - print('Percentage cut: {:3f} %%'.format( - (1. - float(len(valid3)) / np.sum(trace > 0)) * 100.)) + print("number of pixels kept:", len(valid3), "/", np.sum(trace > 0)) + print( + "Percentage cut: {:3f} %%".format( + (1.0 - float(len(valid3)) / np.sum(trace > 0)) * 100.0 + ) + ) weight[valid] = lambda_minus[valid] return weight + def write_maps_a_la_xpure(OutputSkyMap, name_out, output_path): """ Grab sky data from OutputSkyMap and write them into files readable by @@ -122,21 +129,22 @@ def write_maps_a_la_xpure(OutputSkyMap, name_out, output_path): fits_Q[OutputSkyMap.obspix] = Q fits_U[OutputSkyMap.obspix] = U - full_path = os.path.join( - output_path, name_out, 'IQU_{}.fits'.format(name_out)) + full_path = os.path.join(output_path, name_out, "IQU_{}.fits".format(name_out)) - write_healpix_cmbmap(full_path, - data=[fits_I, fits_Q, fits_U], - fits_IDL=False, - coord=None, - colnames=['I', 'Q', 'U'], - partial=False, - nest=False) + write_healpix_cmbmap( + full_path, + data=[fits_I, fits_Q, fits_U], + fits_IDL=False, + coord=None, + colnames=["I", "Q", "U"], + partial=False, + nest=False, + ) del fits_I, fits_Q, fits_U -def write_weights_a_la_xpure(OutputSkyMap, name_out, output_path, epsilon, - HWP=False): + +def write_weights_a_la_xpure(OutputSkyMap, name_out, output_path, epsilon, HWP=False): """ Grab weight and mask from OutputSkyMap and write them into files readable by the software xpure. @@ -159,7 +167,7 @@ def write_weights_a_la_xpure(OutputSkyMap, name_out, output_path, epsilon, """ safe_mkdir(os.path.join(output_path, name_out)) - ## Intensity + # Intensity weight = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside)) if not HWP: @@ -167,59 +175,69 @@ def write_weights_a_la_xpure(OutputSkyMap, name_out, output_path, epsilon, else: weight[OutputSkyMap.obspix] = OutputSkyMap.w0 - full_path = os.path.join( - output_path, name_out, 'Iw_{}.fits'.format(name_out)) - write_healpix_cmbmap(full_path, - data=weight, - fits_IDL=False, - coord=None, - colnames=['I'], - partial=False, - nest=False) + full_path = os.path.join(output_path, name_out, "Iw_{}.fits".format(name_out)) + write_healpix_cmbmap( + full_path, + data=weight, + fits_IDL=False, + coord=None, + colnames=["I"], + partial=False, + nest=False, + ) binary = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside)) mask = np.where((weight > 0))[0] binary[mask] = 1.0 full_path = os.path.join( - output_path, name_out, 'Iw_{}_norm.fits'.format(name_out)) - write_healpix_cmbmap(full_path, - data=binary, - fits_IDL=False, - coord=None, - colnames=['I'], - partial=False, - nest=False) - - ## Polarisation + output_path, name_out, "Iw_{}_norm.fits".format(name_out) + ) + write_healpix_cmbmap( + full_path, + data=binary, + fits_IDL=False, + coord=None, + colnames=["I"], + partial=False, + nest=False, + ) + + # Polarisation weight = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside)) if not HWP: weight[OutputSkyMap.obspix] = qu_weight_mineig( - OutputSkyMap.cc, OutputSkyMap.cs, OutputSkyMap.ss, epsilon) + OutputSkyMap.cc, OutputSkyMap.cs, OutputSkyMap.ss, epsilon + ) else: weight[OutputSkyMap.obspix] = OutputSkyMap.w4 - full_path = os.path.join( - output_path, name_out, 'Pw_{}.fits'.format(name_out)) - write_healpix_cmbmap(full_path, - data=weight, - fits_IDL=False, - coord=None, - colnames=['P'], - partial=False, - nest=False) + full_path = os.path.join(output_path, name_out, "Pw_{}.fits".format(name_out)) + write_healpix_cmbmap( + full_path, + data=weight, + fits_IDL=False, + coord=None, + colnames=["P"], + partial=False, + nest=False, + ) binary = np.zeros((12 * OutputSkyMap.nside * OutputSkyMap.nside)) mask = np.where((weight > 0))[0] binary[mask] = 1.0 full_path = os.path.join( - output_path, name_out, 'Pw_{}_norm.fits'.format(name_out)) - write_healpix_cmbmap(full_path, - data=binary, - fits_IDL=False, - coord=None, - colnames=['P'], - partial=False, - nest=False) + output_path, name_out, "Pw_{}_norm.fits".format(name_out) + ) + write_healpix_cmbmap( + full_path, + data=binary, + fits_IDL=False, + coord=None, + colnames=["P"], + partial=False, + nest=False, + ) + def create_batch(batch_file, name_out, params_s4cmb, params_xpure): """ @@ -234,416 +252,503 @@ def create_batch(batch_file, name_out, params_s4cmb, params_xpure): params_xpure : NormaliseXpureParser instance Object with xpure parameter values. """ - host = os.environ['NERSC_HOST'] - if host == 'edison': + host = os.environ["NERSC_HOST"] + if host == "edison": params_xpure.nproc_per_node = 24 - elif host == 'cori': + elif host == "cori": params_xpure.nproc_per_node = 32 - with open(batch_file, 'w') as f: - print('#!/bin/bash -l', file=f) - print('#SBATCH -p {}'.format(params_xpure.queue), file=f) - print('#SBATCH -N {}'.format(params_xpure.node), file=f) - print('#SBATCH -t {}'.format(params_xpure.time), file=f) - print('#SBATCH -J {}'.format(name_out), file=f) - if host == 'cori': - print('#SBATCH -C haswell', file=f) - - print(' ', file=f) - print('######## TO BE CHANGED BY USER ########', file=f) - print('## Name for files and folders', file=f) - print('name={}'.format(name_out), file=f) - - print(' ', file=f) - print('## Radius for apodization', file=f) - print('radius={}'.format(params_xpure.radius_apodization), file=f) - - print(' ', file=f) - print('## Parameters for XPURE', file=f) - print('LMAX_USER={}'.format(params_xpure.lmax_user), file=f) - print('NSIDE_USER={}'.format(params_s4cmb.nside_out), file=f) - - print(' ', file=f) - print('## Number of simulation (put 1 if not simulation)', file=f) - print('NSIMU_USER=1', file=f) - - print(' ', file=f) - print('## Number of masks (put one if several maps but one common masks)', file=f) - print('NMASK_USER=1', file=f) - - print(' ', file=f) - print('## Number of maps (if > 1, will do the cross correlations)', file=f) - print('NMAPS_USER=1', file=f) - - print(' ', file=f) - print('## Beam and bins', file=f) - print('BEAMFILE={}'.format(params_xpure.beam_file), file=f) - print('BINTAB={}'.format(params_xpure.bin_file), file=f) - - print(' ', file=f) - print('## XPURE running mode (0=xpol, 1=xpure, 2=hybrid)', file=f) - print('MODE_XPURE={}'.format(params_xpure.xpure_mode), file=f) - - print(' ', file=f) - print('## FULL (0=myapodizemask, create_mll and XPURE) or FAST (1=only XPURE) or SEMI-FAST (=2 create_mll and XPURE)', file=f) - print('FAST={}'.format(params_xpure.fast), file=f) - - print(' ', file=f) - print('######## END TO BE CHANGED BY USER ########', file=f) - - print(' ', file=f) - print('#ENVIRONEMENT', file=f) - print('SCRATCH2=/global/cscratch1/sd/peloton', file=f) - if host == 'edison': - print('BINDIR=${HOME}/mapmaker/xpure/xpure/trunk/build/edison', file=f) - elif host == 'cori': - print('BINDIR=${HOME}/mapmaker/xpure/xpure/trunk/build/cori', file=f) - print('THEORYDIR=${SCRATCH2}/s4cmb/additional_files', file=f) - print('OUTPUTMLL=${SCRATCH2}/s4cmb/xpure/cls/${name}', file=f) - print('OUTPUTCL=${SCRATCH2}/s4cmb/xpure/cls/${name}', file=f) - print('mkdir -p ${OUTPUTMLL}', file=f) - print('mkdir -p ${OUTPUTCL}', file=f) - print('MASKDIR=${SCRATCH2}/s4cmb/xpure/masks/${name}', file=f) - print('MAPDIR=${SCRATCH2}/s4cmb/xpure/maps/${name}', file=f) - - print(' ', file=f) - print('####################################################', file=f) - print('# MASK (COMMON MASK)', file=f) - print('####################################################', file=f) - print('cd $MASKDIR', file=f) - - print(' ', file=f) + with open(batch_file, "w") as f: + print("#!/bin/bash -l", file=f) + print("#SBATCH -p {}".format(params_xpure.queue), file=f) + print("#SBATCH -N {}".format(params_xpure.node), file=f) + print("#SBATCH -t {}".format(params_xpure.time), file=f) + print("#SBATCH -J {}".format(name_out), file=f) + if host == "cori": + print("#SBATCH -C haswell", file=f) + + print(" ", file=f) + print("#### TO BE CHANGED BY USER ####", file=f) + print("# Name for files and folders", file=f) + print("name={}".format(name_out), file=f) + + print(" ", file=f) + print("# Radius for apodization", file=f) + print("radius={}".format(params_xpure.radius_apodization), file=f) + + print(" ", file=f) + print("# Parameters for XPURE", file=f) + print("LMAX_USER={}".format(params_xpure.lmax_user), file=f) + print("NSIDE_USER={}".format(params_s4cmb.nside_out), file=f) + + print(" ", file=f) + print("# Number of simulation (put 1 if not simulation)", file=f) + print("NSIMU_USER=1", file=f) + + print(" ", file=f) + print( + "# Number of masks (put one if several maps but one common masks)", + file=f, + ) + print("NMASK_USER=1", file=f) + + print(" ", file=f) + print("# Number of maps (if > 1, will do the cross correlations)", file=f) + print("NMAPS_USER=1", file=f) + + print(" ", file=f) + print("# Beam and bins", file=f) + print("BEAMFILE={}".format(params_xpure.beam_file), file=f) + print("BINTAB={}".format(params_xpure.bin_file), file=f) + + print(" ", file=f) + print("# XPURE running mode (0=xpol, 1=xpure, 2=hybrid)", file=f) + print("MODE_XPURE={}".format(params_xpure.xpure_mode), file=f) + + print(" ", file=f) + txt = "# FULL (0=myapodizemask, create_mll and XPURE) or FAST "\ + + "(1=only XPURE) or SEMI-FAST (=2 create_mll and XPURE)" + print( + txt, + file=f, + ) + print("FAST={}".format(params_xpure.fast), file=f) + + print(" ", file=f) + print("#### END TO BE CHANGED BY USER ####", file=f) + + print(" ", file=f) + print("#ENVIRONEMENT", file=f) + print("SCRATCH2=/global/cscratch1/sd/peloton", file=f) + if host == "edison": + print("BINDIR=${HOME}/mapmaker/xpure/xpure/trunk/build/edison", file=f) + elif host == "cori": + print("BINDIR=${HOME}/mapmaker/xpure/xpure/trunk/build/cori", file=f) + print("THEORYDIR=${SCRATCH2}/s4cmb/additional_files", file=f) + print("OUTPUTMLL=${SCRATCH2}/s4cmb/xpure/cls/${name}", file=f) + print("OUTPUTCL=${SCRATCH2}/s4cmb/xpure/cls/${name}", file=f) + print("mkdir -p ${OUTPUTMLL}", file=f) + print("mkdir -p ${OUTPUTCL}", file=f) + print("MASKDIR=${SCRATCH2}/s4cmb/xpure/masks/${name}", file=f) + print("MAPDIR=${SCRATCH2}/s4cmb/xpure/maps/${name}", file=f) + + print(" ", file=f) + print("##########################", file=f) + print("# MASK (COMMON MASK)", file=f) + print("##########################", file=f) + print("cd $MASKDIR", file=f) + + print(" ", file=f) print('binary_name=$(basename $(find -name "Iw*${name}_norm.fits"))', file=f) print('weight_name=$(basename $(find -name "Iw*${name}.fits"))', file=f) - print('BINARY_MASK_I1=${MASKDIR}/${binary_name}', file=f) - print('WEIGHT_I1=${MASKDIR}/${weight_name}', file=f) - print('APODIZED_MASK_I1=${MASKDIR}/ap30_${weight_name}', file=f) + print("BINARY_MASK_I1=${MASKDIR}/${binary_name}", file=f) + print("WEIGHT_I1=${MASKDIR}/${weight_name}", file=f) + print("APODIZED_MASK_I1=${MASKDIR}/ap30_${weight_name}", file=f) - print(' ', file=f) + print(" ", file=f) print('binary_name=$(basename $(find -name "Pw*${name}_norm.fits"))', file=f) print('weight_name=$(basename $(find -name "Pw*${name}.fits"))', file=f) - print('BINARY_MASK_P1=${MASKDIR}/${binary_name}', file=f) - print('WEIGHT_P1=${MASKDIR}/${weight_name}', file=f) - print('APODIZED_MASK_P1=${MASKDIR}/ap30_${weight_name}', file=f) - - print(' ', file=f) - print('### Intensity', file=f) - print('output_spin0_I1=${MASKDIR}/spin0_I_${name}.fits', file=f) - print('output_spin1_I1=${MASKDIR}/spin1_I_${name}.fits', file=f) - print('output_spin2_I1=${MASKDIR}/spin2_I_${name}.fits', file=f) - - print(' ', file=f) - print('### Polarization', file=f) - print('output_spin0_P1=${MASKDIR}/spin0_P_${name}.fits', file=f) - print('output_spin1_P1=${MASKDIR}/spin1_P_${name}.fits', file=f) - print('output_spin2_P1=${MASKDIR}/spin2_P_${name}.fits', file=f) - - print(' ', file=f) - print('####################################################', file=f) - print('# MAPS', file=f) - print('####################################################', file=f) - print('cd ${MAPDIR}', file=f) - - print(' ', file=f) - print('for (( i=1; i<=${NMAPS_USER}; i++)); do', file=f) + print("BINARY_MASK_P1=${MASKDIR}/${binary_name}", file=f) + print("WEIGHT_P1=${MASKDIR}/${weight_name}", file=f) + print("APODIZED_MASK_P1=${MASKDIR}/ap30_${weight_name}", file=f) + + print(" ", file=f) + print("## Intensity", file=f) + print("output_spin0_I1=${MASKDIR}/spin0_I_${name}.fits", file=f) + print("output_spin1_I1=${MASKDIR}/spin1_I_${name}.fits", file=f) + print("output_spin2_I1=${MASKDIR}/spin2_I_${name}.fits", file=f) + + print(" ", file=f) + print("## Polarization", file=f) + print("output_spin0_P1=${MASKDIR}/spin0_P_${name}.fits", file=f) + print("output_spin1_P1=${MASKDIR}/spin1_P_${name}.fits", file=f) + print("output_spin2_P1=${MASKDIR}/spin2_P_${name}.fits", file=f) + + print(" ", file=f) + print("##########################", file=f) + print("# MAPS", file=f) + print("##########################", file=f) + print("cd ${MAPDIR}", file=f) + + print(" ", file=f) + print("for (( i=1; i<=${NMAPS_USER}; i++)); do", file=f) print(' MAPS[${i}]=$(basename $(find -name "IQU*${name}.fits"))', file=f) - print('done', file=f) - - print(' ', file=f) - print('## CHECK', file=f) - print('for (( i=1; i<=${NMAPS_USER}; i++)); do', file=f) - print(' echo ${MAPS[${i}]}', file=f) - print('done', file=f) - - print(' ', file=f) - print('cd $PBS_O_WORKDIR', file=f) - - print(' ', file=f) - print('#########################################################################################', file=f) - print('# From binary mask to binary apodized mask', file=f) - print('#########################################################################################', file=f) + print("done", file=f) + + print(" ", file=f) + print("# CHECK", file=f) + print("for (( i=1; i<=${NMAPS_USER}; i++)); do", file=f) + print(" echo ${MAPS[${i}]}", file=f) + print("done", file=f) + + print(" ", file=f) + print("cd $PBS_O_WORKDIR", file=f) + + print(" ", file=f) + print( + "################################", + file=f, + ) + print("# From binary mask to binary apodized mask", file=f) + print( + "################################", + file=f, + ) print('if [ "${FAST}" -eq "0" ]', file=f) - print(' then', file=f) - print(' time srun -N 1 -n {} ${{BINDIR}}/myapodizemask ${{BINARY_MASK_I1}} ${{APODIZED_MASK_I1}} -minpix 1 -inside 1 -radius ${{radius}} & time srun -N 1 -n {} ${{BINDIR}}/myapodizemask ${{BINARY_MASK_P1}} ${{APODIZED_MASK_P1}} -minpix 1 -inside 1 -radius ${{radius}}'.format(params_xpure.nproc_apo, params_xpure.nproc_apo), file=f) - print(' wait', file=f) - print('else', file=f) + print(" then", file=f) + txt = " time srun -N 1 -n {} ${{BINDIR}}/myapodizemask "\ + + "${{BINARY_MASK_I1}} ${{APODIZED_MASK_I1}} -minpix 1 -inside 1 "\ + + "-radius ${{radius}} & time srun -N 1 -n {} "\ + + "${{BINDIR}}/myapodizemask ${{BINARY_MASK_P1}} "\ + + "${{APODIZED_MASK_P1}} -minpix 1 -inside 1 -radius ${{radius}}" + print( + txt.format( + params_xpure.nproc_apo, params_xpure.nproc_apo + ), + file=f, + ) + print(" wait", file=f) + print("else", file=f) print(' echo "Go fast - skip myapodizemask"', file=f) - print('fi', file=f) - - print(' ', file=f) - print('#########################################################################################', file=f) - print('# Computation of spin window function', file=f) - print('#########################################################################################', file=f) - - print(' ', file=f) - print('cat > param_all_I11${name}.par << EOF', file=f) - print('##healpix parameters', file=f) - print('###################', file=f) - print('nside = ${NSIDE_USER}', file=f) - print('lmax = ${LMAX_USER} #should not exceed 2*Nside', file=f) - - print(' ', file=f) - print('##input file parameters', file=f) - print('#######################', file=f) - print('maskBinary = ${BINARY_MASK_I1}', file=f) - print('window_spin0 = ${APODIZED_MASK_I1}', file=f) - print('inverseNoise = ${WEIGHT_I1}', file=f) - - print(' ', file=f) - print('##output file parameters', file=f) - print('########################', file=f) - print('output_spin0 = ${output_spin0_I1}', file=f) - print('output_spin1 = ${output_spin1_I1}', file=f) - print('output_spin2 = ${output_spin2_I1}', file=f) - print('EOF', file=f) - - print(' ', file=f) - print('cat > param_all_P11${name}.par << EOF', file=f) - print('##healpix parameters', file=f) - print('###################', file=f) - print('nside = ${NSIDE_USER}', file=f) - print('lmax = ${LMAX_USER} #should not exceed 2*Nside', file=f) - - print(' ', file=f) - print('##input file parameters', file=f) - print('#######################', file=f) - print('maskBinary = ${BINARY_MASK_P1}', file=f) - print('window_spin0 = ${APODIZED_MASK_P1}', file=f) - print('inverseNoise = ${WEIGHT_P1}', file=f) - - print(' ', file=f) - print('##output file parameters', file=f) - print('########################', file=f) - print('output_spin0 = ${output_spin0_P1}', file=f) - print('output_spin1 = ${output_spin1_P1}', file=f) - print('output_spin2 = ${output_spin2_P1}', file=f) - print('EOF', file=f) - - print(' ', file=f) + print("fi", file=f) + + print(" ", file=f) + print( + "################################", + file=f, + ) + print("# Computation of spin window function", file=f) + print( + "################################", + file=f, + ) + + print(" ", file=f) + print("cat > param_all_I11${name}.par << EOF", file=f) + print("#healpix parameters", file=f) + print("##########", file=f) + print("nside = ${NSIDE_USER}", file=f) + print("lmax = ${LMAX_USER} #should not exceed 2*Nside", file=f) + + print(" ", file=f) + print("#input file parameters", file=f) + print("############", file=f) + print("maskBinary = ${BINARY_MASK_I1}", file=f) + print("window_spin0 = ${APODIZED_MASK_I1}", file=f) + print("inverseNoise = ${WEIGHT_I1}", file=f) + + print(" ", file=f) + print("#output file parameters", file=f) + print("############", file=f) + print("output_spin0 = ${output_spin0_I1}", file=f) + print("output_spin1 = ${output_spin1_I1}", file=f) + print("output_spin2 = ${output_spin2_I1}", file=f) + print("EOF", file=f) + + print(" ", file=f) + print("cat > param_all_P11${name}.par << EOF", file=f) + print("#healpix parameters", file=f) + print("##########", file=f) + print("nside = ${NSIDE_USER}", file=f) + print("lmax = ${LMAX_USER} #should not exceed 2*Nside", file=f) + + print(" ", file=f) + print("#input file parameters", file=f) + print("############", file=f) + print("maskBinary = ${BINARY_MASK_P1}", file=f) + print("window_spin0 = ${APODIZED_MASK_P1}", file=f) + print("inverseNoise = ${WEIGHT_P1}", file=f) + + print(" ", file=f) + print("#output file parameters", file=f) + print("############", file=f) + print("output_spin0 = ${output_spin0_P1}", file=f) + print("output_spin1 = ${output_spin1_P1}", file=f) + print("output_spin2 = ${output_spin2_P1}", file=f) + print("EOF", file=f) + + print(" ", file=f) print('if [ "${FAST}" -eq "0" ]', file=f) - print(' then', file=f) - print(' time srun -N {} -n {} ${{BINDIR}}/scalar2spin param_all_I11${{name}}.par >& output_scalar2spinI11${{name}} & time srun -N {} -n {} ${{BINDIR}}/scalar2spin param_all_P11${{name}}.par >& output_scalar2spinP11${{name}}'.format(int(params_xpure.nproc_scalar_to_spin // params_xpure.nproc_per_node), params_xpure.nproc_scalar_to_spin, int(params_xpure.nproc_scalar_to_spin // params_xpure.nproc_per_node), params_xpure.nproc_scalar_to_spin), file=f) - print(' wait', file=f) - print('else', file=f) + print(" then", file=f) + txt = " time srun -N {} -n {} ${{BINDIR}}/scalar2spin "\ + + "param_all_I11${{name}}.par >& output_scalar2spinI11${{name}} "\ + + "& time srun -N {} -n {} "\ + + "${{BINDIR}}/scalar2spin param_all_P11${{name}}.par "\ + + ">& output_scalar2spinP11${{name}}" + print( + txt.format( + int(params_xpure.nproc_scalar_to_spin // params_xpure.nproc_per_node), + params_xpure.nproc_scalar_to_spin, + int(params_xpure.nproc_scalar_to_spin // params_xpure.nproc_per_node), + params_xpure.nproc_scalar_to_spin, + ), + file=f, + ) + print(" wait", file=f) + print("else", file=f) print(' echo "Go fast - skip scalar2spin"', file=f) - print('fi', file=f) - - print(' ', file=f) - print('#########################################################################################', file=f) - print('# X2pure', file=f) - print('#########################################################################################', file=f) - print('cd ${OUTPUTMLL}', file=f) - - print(' ', file=f) - print('#########################################################################################', file=f) - print('# Mode-mode mixing matrix', file=f) - print('#########################################################################################', file=f) - print('cat > createMll.par << EOF', file=f) - - print(' ', file=f) - print('######### MODE #############', file=f) - print('# 0 : Standard formalism', file=f) - print('# 1 : Pure formalism', file=f) - print('# 2 : Hybrid formalism', file=f) - print('############################', file=f) - print('mode = ${MODE_XPURE}', file=f) - - print(' ', file=f) - print('############ SETUP #########', file=f) - print('nside = ${NSIDE_USER}', file=f) - print('lmax = ${LMAX_USER}', file=f) - print('nmask = ${NMASK_USER}', file=f) - - print(' ', file=f) - print('EOF', file=f) - - print(' ', file=f) - print('for (( i=1; i<=${NMASK_USER}; i++)); do', file=f) - print(' ind=$(($i - 1))', file=f) - print(' cat >> createMll.par << EOF', file=f) - - print(' ', file=f) - print('maskfile${i}_T = ${output_spin0_I1}', file=f) - - print(' ', file=f) - print('maskfile${i}_E_spin0 = ${output_spin0_P1}', file=f) - print('maskfile${i}_E_spin1 = ${output_spin1_P1}', file=f) - print('maskfile${i}_E_spin2 = ${output_spin2_P1}', file=f) - - print(' ', file=f) - print('maskfile${i}_B_spin0 = ${output_spin0_P1}', file=f) - print('maskfile${i}_B_spin1 = ${output_spin1_P1}', file=f) - print('maskfile${i}_B_spin2 = ${output_spin2_P1}', file=f) - - print(' ', file=f) - print('mllfile_TT_TT_${i} = ${OUTPUTMLL}/mll_TT_TT_BinMask${i}.fits', file=f) - - print(' ', file=f) - print('mllfile_EE_EE_${i} = ${OUTPUTMLL}/mll_spinEE_EE_pcg${i}.fits', file=f) - print('mllfile_EE_BB_${i} = ${OUTPUTMLL}/mll_spinEE_BB_pcg${i}.fits', file=f) - print('mllfile_EE_EB_${i} = ${OUTPUTMLL}/mll_spinEE_EB_pcg${i}.fits', file=f) - print('mllfile_BB_BB_${i} = ${OUTPUTMLL}/mll_spinBB_BB_pcg${i}.fits', file=f) - print('mllfile_BB_EE_${i} = ${OUTPUTMLL}/mll_spinBB_EE_pcg${i}.fits', file=f) - print('mllfile_BB_EB_${i} = ${OUTPUTMLL}/mll_spinBB_EB_pcg${i}.fits', file=f) - - print(' ', file=f) - print('mllfile_TE_TE_${i} = ${OUTPUTMLL}/mll_spinTE_TE_pcg${i}.fits', file=f) - print('mllfile_TE_TB_${i} = ${OUTPUTMLL}/mll_spinTE_TB_pcg${i}.fits', file=f) - print('mllfile_TB_TE_${i} = ${OUTPUTMLL}/mll_spinTB_TE_pcg${i}.fits', file=f) - print('mllfile_TB_TB_${i} = ${OUTPUTMLL}/mll_spinTB_TB_pcg${i}.fits', file=f) - - print(' ', file=f) - print('mllfile_EB_EB_${i} = ${OUTPUTMLL}/mll_spinEB_EB_pcg${i}.fits', file=f) - print('mllfile_EB_EE_${i} = ${OUTPUTMLL}/mll_spinEB_EE_pcg${i}.fits', file=f) - print('mllfile_EB_BB_${i} = ${OUTPUTMLL}/mll_spinEB_BB_pcg${i}.fits', file=f) - - print(' ', file=f) - print('EOF', file=f) - - print(' ', file=f) - print('done', file=f) - - print(' ', file=f) + print("fi", file=f) + + print(" ", file=f) + print( + "################################", + file=f, + ) + print("# X2pure", file=f) + print( + "################################", + file=f, + ) + print("cd ${OUTPUTMLL}", file=f) + + print(" ", file=f) + print( + "################################", + file=f, + ) + print("# Mode-mode mixing matrix", file=f) + print( + "################################", + file=f, + ) + print("cat > createMll.par << EOF", file=f) + + print(" ", file=f) + print("##### MODE #######", file=f) + print("# 0 : Standard formalism", file=f) + print("# 1 : Pure formalism", file=f) + print("# 2 : Hybrid formalism", file=f) + print("##############", file=f) + print("mode = ${MODE_XPURE}", file=f) + + print(" ", file=f) + print("###### SETUP #####", file=f) + print("nside = ${NSIDE_USER}", file=f) + print("lmax = ${LMAX_USER}", file=f) + print("nmask = ${NMASK_USER}", file=f) + + print(" ", file=f) + print("EOF", file=f) + + print(" ", file=f) + print("for (( i=1; i<=${NMASK_USER}; i++)); do", file=f) + print(" ind=$(($i - 1))", file=f) + print(" cat >> createMll.par << EOF", file=f) + + print(" ", file=f) + print("maskfile${i}_T = ${output_spin0_I1}", file=f) + + print(" ", file=f) + print("maskfile${i}_E_spin0 = ${output_spin0_P1}", file=f) + print("maskfile${i}_E_spin1 = ${output_spin1_P1}", file=f) + print("maskfile${i}_E_spin2 = ${output_spin2_P1}", file=f) + + print(" ", file=f) + print("maskfile${i}_B_spin0 = ${output_spin0_P1}", file=f) + print("maskfile${i}_B_spin1 = ${output_spin1_P1}", file=f) + print("maskfile${i}_B_spin2 = ${output_spin2_P1}", file=f) + + print(" ", file=f) + print("mllfile_TT_TT_${i} = ${OUTPUTMLL}/mll_TT_TT_BinMask${i}.fits", file=f) + + print(" ", file=f) + print("mllfile_EE_EE_${i} = ${OUTPUTMLL}/mll_spinEE_EE_pcg${i}.fits", file=f) + print("mllfile_EE_BB_${i} = ${OUTPUTMLL}/mll_spinEE_BB_pcg${i}.fits", file=f) + print("mllfile_EE_EB_${i} = ${OUTPUTMLL}/mll_spinEE_EB_pcg${i}.fits", file=f) + print("mllfile_BB_BB_${i} = ${OUTPUTMLL}/mll_spinBB_BB_pcg${i}.fits", file=f) + print("mllfile_BB_EE_${i} = ${OUTPUTMLL}/mll_spinBB_EE_pcg${i}.fits", file=f) + print("mllfile_BB_EB_${i} = ${OUTPUTMLL}/mll_spinBB_EB_pcg${i}.fits", file=f) + + print(" ", file=f) + print("mllfile_TE_TE_${i} = ${OUTPUTMLL}/mll_spinTE_TE_pcg${i}.fits", file=f) + print("mllfile_TE_TB_${i} = ${OUTPUTMLL}/mll_spinTE_TB_pcg${i}.fits", file=f) + print("mllfile_TB_TE_${i} = ${OUTPUTMLL}/mll_spinTB_TE_pcg${i}.fits", file=f) + print("mllfile_TB_TB_${i} = ${OUTPUTMLL}/mll_spinTB_TB_pcg${i}.fits", file=f) + + print(" ", file=f) + print("mllfile_EB_EB_${i} = ${OUTPUTMLL}/mll_spinEB_EB_pcg${i}.fits", file=f) + print("mllfile_EB_EE_${i} = ${OUTPUTMLL}/mll_spinEB_EE_pcg${i}.fits", file=f) + print("mllfile_EB_BB_${i} = ${OUTPUTMLL}/mll_spinEB_BB_pcg${i}.fits", file=f) + + print(" ", file=f) + print("EOF", file=f) + + print(" ", file=f) + print("done", file=f) + + print(" ", file=f) print('if [ "${FAST}" -eq "0" ]', file=f) - print(' then', file=f) - print(' time srun -N {} -n {} ${{BINDIR}}/x2pure_create_mll createMll.par'.format(int(params_xpure.nproc_mll//params_xpure.nproc_per_node), params_xpure.nproc_mll), file=f) - - print(' rm -f createMll.par', file=f) + print(" then", file=f) + txt = " time srun -N {} -n {} "\ + + "${{BINDIR}}/x2pure_create_mll createMll.par" + print( + txt.format( + int(params_xpure.nproc_mll // params_xpure.nproc_per_node), + params_xpure.nproc_mll, + ), + file=f, + ) + + print(" rm -f createMll.par", file=f) print('elif [ "${FAST}" -eq "2" ]', file=f) - print(' then', file=f) - print(' time srun -N {} -n {} ${{BINDIR}}/x2pure_create_mll createMll.par'.format(int(params_xpure.nproc_mll//params_xpure.nproc_per_node), params_xpure.nproc_mll), file=f) - print(' rm -f createMll.par', file=f) - print('else', file=f) + print(" then", file=f) + txt = " time srun -N {} -n {} "\ + + "${{BINDIR}}/x2pure_create_mll createMll.par" + print( + txt.format( + int(params_xpure.nproc_mll // params_xpure.nproc_per_node), + params_xpure.nproc_mll, + ), + file=f, + ) + print(" rm -f createMll.par", file=f) + print("else", file=f) print(' echo "Go fast - skip x2pure_create_mll"', file=f) - print('fi', file=f) - - print(' ', file=f) - print('#########################################################################################', file=f) - print('# X2pure', file=f) - print('#########################################################################################', file=f) - print('for (( n=0; n<${NSIMU_USER}; n++ )); do', file=f) - - print(' ', file=f) - print(' echo "************************ simu $n ************************"', file=f) - - print(' ', file=f) - print(' num=$n', file=f) - - print(' ', file=f) - print(' #CREATE PARAMETER FILE', file=f) - print(' cat > xpure.par << _EOF_', file=f) - - print(' ', file=f) - print('mode = ${MODE_XPURE}', file=f) - - print(' ', file=f) - print('nside = ${NSIDE_USER}', file=f) - print('nmaps = ${NMAPS_USER}', file=f) - print('nmasks = ${NMASK_USER}', file=f) - - print(' ', file=f) - print('_EOF_', file=f) - - print(' ', file=f) - print(' for ((j=1; j<=${NMAPS_USER}; j++)); do', file=f) - - print(' ', file=f) - print(' cat >> xpure.par << _EOF_', file=f) - - print(' ', file=f) - print('bellfile${j} = ${THEORYDIR}/${BEAMFILE}', file=f) - print('mapfile${j} = ${MAPDIR}/${MAPS[${j}]}', file=f) - - print(' ', file=f) - print('_EOF_', file=f) - - print(' ', file=f) - print(' done', file=f) - - print(' ', file=f) - print(' cat >> xpure.par << _EOF_', file=f) - - print(' ', file=f) - print('lmaxSim = ${LMAX_USER}', file=f) - - print(' ', file=f) - print('_EOF_', file=f) - - print(' ', file=f) - print(' for(( i=1; i<=${NMASK_USER}; i++)); do', file=f) - print(' ind=$(($i - 1))', file=f) - print(' cat >> xpure.par << EOF', file=f) - - print(' ', file=f) - print('maskfile${i}_T = ${output_spin0_I1}', file=f) - - print(' ', file=f) - print('maskfile${i}_E_spin0 = ${output_spin0_P1}', file=f) - print('maskfile${i}_E_spin1 = ${output_spin1_P1}', file=f) - print('maskfile${i}_E_spin2 = ${output_spin2_P1}', file=f) - - print(' ', file=f) - print('maskfile${i}_B_spin0 = ${output_spin0_P1}', file=f) - print('maskfile${i}_B_spin1 = ${output_spin1_P1}', file=f) - print('maskfile${i}_B_spin2 = ${output_spin2_P1}', file=f) - - print(' ', file=f) - print('mllfile_TT_TT_${i} = ${OUTPUTMLL}/mll_TT_TT_BinMask${i}.fits', file=f) - - print(' ', file=f) - print('mllfile_EE_EE_${i} = ${OUTPUTMLL}/mll_spinEE_EE_pcg${i}.fits', file=f) - print('mllfile_EE_BB_${i} = ${OUTPUTMLL}/mll_spinEE_BB_pcg${i}.fits', file=f) - print('mllfile_EE_EB_${i} = ${OUTPUTMLL}/mll_spinEE_EB_pcg${i}.fits', file=f) - print('mllfile_BB_BB_${i} = ${OUTPUTMLL}/mll_spinBB_BB_pcg${i}.fits', file=f) - print('mllfile_BB_EE_${i} = ${OUTPUTMLL}/mll_spinBB_EE_pcg${i}.fits', file=f) - print('mllfile_BB_EB_${i} = ${OUTPUTMLL}/mll_spinBB_EB_pcg${i}.fits', file=f) - - print(' ', file=f) - print('mllfile_TE_TE_${i} = ${OUTPUTMLL}/mll_spinTE_TE_pcg${i}.fits', file=f) - print('mllfile_TE_TB_${i} = ${OUTPUTMLL}/mll_spinTE_TB_pcg${i}.fits', file=f) - print('mllfile_TB_TE_${i} = ${OUTPUTMLL}/mll_spinTB_TE_pcg${i}.fits', file=f) - print('mllfile_TB_TB_${i} = ${OUTPUTMLL}/mll_spinTB_TB_pcg${i}.fits', file=f) - - print(' ', file=f) - print('mllfile_EB_EB_${i} = ${OUTPUTMLL}/mll_spinEB_EB_pcg${i}.fits', file=f) - print('mllfile_EB_EE_${i} = ${OUTPUTMLL}/mll_spinEB_EE_pcg${i}.fits', file=f) - print('mllfile_EB_BB_${i} = ${OUTPUTMLL}/mll_spinEB_BB_pcg${i}.fits', file=f) - - print(' ', file=f) - print('EOF', file=f) - print(' done', file=f) - - print(' ', file=f) - print(' cat >> xpure.par << _EOF_', file=f) - print('noise_biasT_1 = 0.', file=f) - print('noise_biasT_2_2 = 0.', file=f) - print('noise_biasT_1_2 = 0.', file=f) - - print(' ', file=f) - print('noise_biasP_1 = 0.', file=f) - print('noise_biasP_2_2 = 0.', file=f) - print('noise_biasP_1_2 = 0.', file=f) - - print(' ', file=f) - print('bintab = ${THEORYDIR}/${BINTAB}', file=f) - print('#mask_list = ${INPUT}/bin2mask_43bins.fits', file=f) - print('_EOF_', file=f) - - print(' ', file=f) - print('cat >> xpure.par << _EOF_', file=f) - print('pseudofile = ${OUTPUTCL}/pseudopure_pcg_${name}_$num', file=f) - print('cellfile = ${OUTPUTCL}/cellpure_pbear_${name}_$num', file=f) - print('lmax = ${LMAX_USER}', file=f) - print('_EOF_', file=f) - print(' #RUN', file=f) - print(' time srun -N {} -n {} ${{BINDIR}}/x2pure xpure.par'.format(int(params_xpure.nproc_xpure//params_xpure.nproc_per_node), params_xpure.nproc_xpure), file=f) - - print(' ', file=f) - print('done', file=f) + print("fi", file=f) + + print(" ", file=f) + print( + "#################################", + file=f, + ) + print("# X2pure", file=f) + print( + "################################", + file=f, + ) + print("for (( n=0; n<${NSIMU_USER}; n++ )); do", file=f) + + print(" ", file=f) + print( + ' echo "************************ simu $n ************************"', + file=f, + ) + + print(" ", file=f) + print(" num=$n", file=f) + + print(" ", file=f) + print(" #CREATE PARAMETER FILE", file=f) + print(" cat > xpure.par << _EOF_", file=f) + + print(" ", file=f) + print("mode = ${MODE_XPURE}", file=f) + + print(" ", file=f) + print("nside = ${NSIDE_USER}", file=f) + print("nmaps = ${NMAPS_USER}", file=f) + print("nmasks = ${NMASK_USER}", file=f) + + print(" ", file=f) + print("_EOF_", file=f) + + print(" ", file=f) + print(" for ((j=1; j<=${NMAPS_USER}; j++)); do", file=f) + + print(" ", file=f) + print(" cat >> xpure.par << _EOF_", file=f) + + print(" ", file=f) + print("bellfile${j} = ${THEORYDIR}/${BEAMFILE}", file=f) + print("mapfile${j} = ${MAPDIR}/${MAPS[${j}]}", file=f) + + print(" ", file=f) + print("_EOF_", file=f) + + print(" ", file=f) + print(" done", file=f) + + print(" ", file=f) + print(" cat >> xpure.par << _EOF_", file=f) + + print(" ", file=f) + print("lmaxSim = ${LMAX_USER}", file=f) + + print(" ", file=f) + print("_EOF_", file=f) + + print(" ", file=f) + print(" for(( i=1; i<=${NMASK_USER}; i++)); do", file=f) + print(" ind=$(($i - 1))", file=f) + print(" cat >> xpure.par << EOF", file=f) + + print(" ", file=f) + print("maskfile${i}_T = ${output_spin0_I1}", file=f) + + print(" ", file=f) + print("maskfile${i}_E_spin0 = ${output_spin0_P1}", file=f) + print("maskfile${i}_E_spin1 = ${output_spin1_P1}", file=f) + print("maskfile${i}_E_spin2 = ${output_spin2_P1}", file=f) + + print(" ", file=f) + print("maskfile${i}_B_spin0 = ${output_spin0_P1}", file=f) + print("maskfile${i}_B_spin1 = ${output_spin1_P1}", file=f) + print("maskfile${i}_B_spin2 = ${output_spin2_P1}", file=f) + + print(" ", file=f) + print("mllfile_TT_TT_${i} = ${OUTPUTMLL}/mll_TT_TT_BinMask${i}.fits", file=f) + + print(" ", file=f) + print("mllfile_EE_EE_${i} = ${OUTPUTMLL}/mll_spinEE_EE_pcg${i}.fits", file=f) + print("mllfile_EE_BB_${i} = ${OUTPUTMLL}/mll_spinEE_BB_pcg${i}.fits", file=f) + print("mllfile_EE_EB_${i} = ${OUTPUTMLL}/mll_spinEE_EB_pcg${i}.fits", file=f) + print("mllfile_BB_BB_${i} = ${OUTPUTMLL}/mll_spinBB_BB_pcg${i}.fits", file=f) + print("mllfile_BB_EE_${i} = ${OUTPUTMLL}/mll_spinBB_EE_pcg${i}.fits", file=f) + print("mllfile_BB_EB_${i} = ${OUTPUTMLL}/mll_spinBB_EB_pcg${i}.fits", file=f) + + print(" ", file=f) + print("mllfile_TE_TE_${i} = ${OUTPUTMLL}/mll_spinTE_TE_pcg${i}.fits", file=f) + print("mllfile_TE_TB_${i} = ${OUTPUTMLL}/mll_spinTE_TB_pcg${i}.fits", file=f) + print("mllfile_TB_TE_${i} = ${OUTPUTMLL}/mll_spinTB_TE_pcg${i}.fits", file=f) + print("mllfile_TB_TB_${i} = ${OUTPUTMLL}/mll_spinTB_TB_pcg${i}.fits", file=f) + + print(" ", file=f) + print("mllfile_EB_EB_${i} = ${OUTPUTMLL}/mll_spinEB_EB_pcg${i}.fits", file=f) + print("mllfile_EB_EE_${i} = ${OUTPUTMLL}/mll_spinEB_EE_pcg${i}.fits", file=f) + print("mllfile_EB_BB_${i} = ${OUTPUTMLL}/mll_spinEB_BB_pcg${i}.fits", file=f) + + print(" ", file=f) + print("EOF", file=f) + print(" done", file=f) + + print(" ", file=f) + print(" cat >> xpure.par << _EOF_", file=f) + print("noise_biasT_1 = 0.", file=f) + print("noise_biasT_2_2 = 0.", file=f) + print("noise_biasT_1_2 = 0.", file=f) + + print(" ", file=f) + print("noise_biasP_1 = 0.", file=f) + print("noise_biasP_2_2 = 0.", file=f) + print("noise_biasP_1_2 = 0.", file=f) + + print(" ", file=f) + print("bintab = ${THEORYDIR}/${BINTAB}", file=f) + print("#mask_list = ${INPUT}/bin2mask_43bins.fits", file=f) + print("_EOF_", file=f) + + print(" ", file=f) + print("cat >> xpure.par << _EOF_", file=f) + print("pseudofile = ${OUTPUTCL}/pseudopure_pcg_${name}_$num", file=f) + print("cellfile = ${OUTPUTCL}/cellpure_pbear_${name}_$num", file=f) + print("lmax = ${LMAX_USER}", file=f) + print("_EOF_", file=f) + print(" #RUN", file=f) + print( + " time srun -N {} -n {} ${{BINDIR}}/x2pure xpure.par".format( + int(params_xpure.nproc_xpure // params_xpure.nproc_per_node), + params_xpure.nproc_xpure, + ), + file=f, + ) + + print(" ", file=f) + print("done", file=f) if __name__ == "__main__": import doctest + if np.__version__ >= "1.14.0": np.set_printoptions(legacy="1.13") doctest.testmod() From cebfc1f31a7c06df555007f68aee510cf7c94ac7 Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 11:01:37 +0100 Subject: [PATCH 11/15] Clean celestial_body_trajectories --- s4cmb/celestial_body_trajectories.py | 45 ++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 12 deletions(-) diff --git a/s4cmb/celestial_body_trajectories.py b/s4cmb/celestial_body_trajectories.py index a2c7480..e642097 100755 --- a/s4cmb/celestial_body_trajectories.py +++ b/s4cmb/celestial_body_trajectories.py @@ -11,10 +11,19 @@ from datetime import datetime, date, time, timedelta import numpy as np -class celestial_trajectory(): + +class celestial_trajectory: """Predict the trajectory of a body (Sun, Moon, ...) """ - def __init__(self, body, lon_observer, lat_observer, elevation_observer, - year, tz_offset=0.): + + def __init__( + self, + body, + lon_observer, + lat_observer, + elevation_observer, + year, + tz_offset=0.0, + ): """ Parameters ---------- @@ -38,8 +47,20 @@ def __init__(self, body, lon_observer, lat_observer, elevation_observer, self.year = year self.tz_offset = tz_offset - self.months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', - 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + self.months = [ + "Jan", + "Feb", + "Mar", + "Apr", + "May", + "Jun", + "Jul", + "Aug", + "Sep", + "Oct", + "Nov", + "Dec", + ] self.altaz = [] self.radec = [] @@ -52,7 +73,7 @@ def initialise_observer(self): """ Set the Longitute and Latitude of the Observer. """ - ## Define the observer + # Define the observer self.ob = ephem.Observer() self.ob.lat, self.ob.lon = self.lat_observer, self.lon_observer self.ob.elevation = self.elevation_observer @@ -141,20 +162,19 @@ def traj_of_body_oneday(self, date): WTF? """ - ## Update the date - date = datetime.combine(date, time(12)) - timedelta( - hours=self.tz_offset) + # Update the date + date = datetime.combine(date, time(12)) - timedelta(hours=self.tz_offset) self.ob.date = date - ## Alt/Az + # Alt/Az ALTAZ = self.alt_az(date) self.altaz.append(ALTAZ) - ## RA/Dec + # RA/Dec RADEC = self.radec_of(ALTAZ) self.radec.append(RADEC) - ## Theta/Phi + # Theta/Phi THETAPHI = self.radec2thetaphi(RADEC) self.thetaphi.append(THETAPHI) @@ -188,4 +208,5 @@ def traj_of_body_oneyear(self): if __name__ == "__main__": import doctest + doctest.testmod() From c159d536302a478a4399d6314c557d73392f5da2 Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 11:03:45 +0100 Subject: [PATCH 12/15] Clean tools --- s4cmb/tools.py | 76 +++++++++++++++++++++++++++++++++++--------------- 1 file changed, 53 insertions(+), 23 deletions(-) diff --git a/s4cmb/tools.py b/s4cmb/tools.py index 5a56af4..02829b0 100755 --- a/s4cmb/tools.py +++ b/s4cmb/tools.py @@ -2,9 +2,10 @@ import healpy as hp """ -Module including diverse tools to manipulate alms and maps . +Module including diverse tools to manipulate alms and maps . """ + def get_healpix_ring_pixel_layout(nside, th_idx): """Healpix ring layout. @@ -14,47 +15,57 @@ def get_healpix_ring_pixel_layout(nside, th_idx): """ ith = th_idx + 1 nrings = 2 * nside - assert 1 <= ith <= 4 * nside -1, (ith, nrings) + assert 1 <= ith <= 4 * nside - 1, (ith, nrings) if ith > nrings: - startpix, nphi, kphi0, cth, sth = get_healpix_ring_pixel_layout(nside,ith - 2 * (ith - nrings) -1) - return 12 * nside ** 2 - startpix -nphi, nphi, kphi0, -cth, sth - dth1 = 1. / 3. / nside ** 2 - dth2 = 2. / 3. / nside - dst1 = 1. / (np.sqrt(6.) * nside) - if ith < nside: # polar cap (north) - cth = 1. - ith ** 2 * dth1 + startpix, nphi, kphi0, cth, sth = get_healpix_ring_pixel_layout( + nside, ith - 2 * (ith - nrings) - 1 + ) + return 12 * nside ** 2 - startpix - nphi, nphi, kphi0, -cth, sth + dth1 = 1.0 / 3.0 / nside ** 2 + dth2 = 2.0 / 3.0 / nside + dst1 = 1.0 / (np.sqrt(6.0) * nside) + if ith < nside: # polar cap (north) + cth = 1.0 - ith ** 2 * dth1 nphi = 4 * ith kphi0 = 1 - sth = np.sin( 2.0 * np.arcsin(ith * dst1) ) # sin(theta) + sth = np.sin(2.0 * np.arcsin(ith * dst1)) startpix = 2 * ith * (ith - 1) else: cth = (2 * nside - ith) * dth2 nphi = 4 * nside kphi0 = (ith + 1 - nside) % 2 - sth = np.sqrt((1.0 -cth) * (1.0+cth)) #! sin(theta) + sth = np.sqrt((1.0 - cth) * (1.0 + cth)) startpix = 2 * nside * (nside - 1) + (ith - nside) * int(nphi) return startpix, nphi, kphi0, cth, sth + def get_alpha_raise(s, lmax): """Response coefficient of spin-s spherical harmonic to spin raising operator. Author: Julien Carron (j.carron@sussex.ac.uk) """ ret = np.zeros(lmax + 1, dtype=float) - ret[abs(s):] = np.sqrt(np.arange(abs(s) -s, lmax - s + 1) * np.arange(abs(s) + s + 1, lmax + s + 2)) + ret[abs(s):] = np.sqrt( + np.arange(abs(s) - s, lmax - s + 1) * np.arange(abs(s) + s + 1, lmax + s + 2) + ) return ret + def get_alpha_lower(s, lmax): """Response coefficient of spin-s spherical harmonic to spin lowering operator. Author: Julien Carron (j.carron@sussex.ac.uk) """ ret = np.zeros(lmax + 1, dtype=float) - ret[abs(s):] = -np.sqrt(np.arange(s + abs(s), lmax + s + 1) * np.arange(abs(s) - s + 1, lmax - s + 2)) + ret[abs(s):] = -np.sqrt( + np.arange(s + abs(s), lmax + s + 1) * np.arange(abs(s) - s + 1, lmax - s + 2) + ) return ret -def alm2map_spin_der1(gclm, nside, spin, zbounds=(-1., 1.), ret_slice=None): - """Returns spin-s transform '_{s}d' of alm, together with d/dtheta _{s}d and 1/sin tht d/dphi _{s}d. + +def alm2map_spin_der1(gclm, nside, spin, zbounds=(-1.0, 1.0), ret_slice=None): + """Returns spin-s transform '_{s}d' of alm, + together with d/dtheta _{s}d and 1/sin tht d/dphi _{s}d. This crude version has three calls to spin-weight harmonics alm2map_spin. @@ -64,28 +75,47 @@ def alm2map_spin_der1(gclm, nside, spin, zbounds=(-1., 1.), ret_slice=None): assert hp.Alm.getlmax(gclm[0].size) == hp.Alm.getlmax(gclm[1].size) lmax = hp.Alm.getlmax(gclm[0].size) zbounds = np.sort(np.array(zbounds)) - # shape (2, 12 * nside ** 2), first entry = real part, second entry imaginary part. + # shape (2, 12 * nside ** 2), + # first entry = real part, second entry imaginary part. _sd = np.array(hp.alm2map_spin(gclm, nside, spin, lmax)) if spin > 1: - _gclm = [hp.almxfl(gclm[0], get_alpha_lower(spin, lmax)), hp.almxfl(gclm[1], get_alpha_lower(spin, lmax))] + _gclm = [ + hp.almxfl(gclm[0], get_alpha_lower(spin, lmax)), + hp.almxfl(gclm[1], get_alpha_lower(spin, lmax)), + ] _sm1d = np.array(hp.alm2map_spin(_gclm, nside, spin - 1, lmax)) else: - _sm1d = -np.array([hp.alm2map(hp.almxfl(gclm[0], get_alpha_lower(spin, lmax)), nside,verbose=False), - hp.alm2map(hp.almxfl(gclm[1], get_alpha_lower(spin, lmax)), nside,verbose=False)]) + _sm1d = -np.array( + [ + hp.alm2map( + hp.almxfl(gclm[0], get_alpha_lower(spin, lmax)), + nside, + verbose=False, + ), + hp.alm2map( + hp.almxfl(gclm[1], get_alpha_lower(spin, lmax)), + nside, + verbose=False, + ), + ] + ) - _gclm = [hp.almxfl(gclm[0], get_alpha_raise(spin, lmax)), hp.almxfl(gclm[1], get_alpha_raise(spin, lmax))] - _sp1d = np.array(hp.alm2map_spin(_gclm, nside, spin + 1 , lmax)) + _gclm = [ + hp.almxfl(gclm[0], get_alpha_raise(spin, lmax)), + hp.almxfl(gclm[1], get_alpha_raise(spin, lmax)), + ] + _sp1d = np.array(hp.alm2map_spin(_gclm, nside, spin + 1, lmax)) d_dth = -0.5 * (_sp1d + _sm1d) d_dphi_sin0 = 0.5 * np.array([-_sp1d[1] + _sm1d[1], _sp1d[0] - _sm1d[0]]) - for iring in range(4 * nside -1): + for iring in range(4 * nside - 1): startpix, nphi, kphi0, cth, sth = get_healpix_ring_pixel_layout(nside, iring) if zbounds[0] <= cth <= zbounds[1]: slic = slice(startpix, startpix + nphi) d_dphi_sin0[1, slic] -= spin * (cth / sth) * _sd[0, slic] d_dphi_sin0[0, slic] += spin * (cth / sth) * _sd[1, slic] if ret_slice is not None: - return _sd[:, ret_slice], d_dth[:, ret_slice], d_dphi_sin0[:,ret_slice] + return _sd[:, ret_slice], d_dth[:, ret_slice], d_dphi_sin0[:, ret_slice] return _sd, d_dth, d_dphi_sin0 From 720a3b27cfcca0c01a88a6419f10c2c8f8abee16 Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 11:19:05 +0100 Subject: [PATCH 13/15] Fix bug in tod --- s4cmb/tod.py | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/s4cmb/tod.py b/s4cmb/tod.py index 7bb48a7..e84d586 100755 --- a/s4cmb/tod.py +++ b/s4cmb/tod.py @@ -954,6 +954,30 @@ def map2tod(self, ch): # Compute pointing matrix index_global, index_local = self.get_pixel_indices(ra, dec) + # Perturbed boresight pointing values + if self.pointing_perturbed is not None: + ra_perturbed, dec_perturbed, pa_perturbed = \ + self.pointing_perturbed.offset_detector(azd, eld) + + # Perturbed boresight pointing values + index_global_perturbed, index_local_perturbed = self.get_pixel_indices( + ra_perturbed, dec_perturbed + ) + + else: + ra_perturbed = ra + dec_perturbed = dec + # pa_perturbed = pa + index_global_perturbed = index_global + # index_local_perturbed = index_local + + if ((self.projection == 'healpix') & (index_local is None)): + # Using a pointer not to increase memory usage + index_local = index_global + + # Perturbed boresight pointing values + # index_local_perturbed = index_global_perturbed + # For flat projection, one needs to flip the sign of U # w.r.t to the full-sky basis (IAU angle convention) if self.projection == "flat": @@ -1043,7 +1067,7 @@ def map2tod(self, ch): pol_ang_pair, pol_ang2_pair = self.compute_simpolangle( ch, pa_pair, polangle_err=False ) - except: + except ValueError: pol_ang_pair = pol_ang pol_ang2_pair = pol_ang2 From cb6b0aade1435707a20c67301f6f06b2fd8b957f Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 11:33:40 +0100 Subject: [PATCH 14/15] Fix bug in tod --- s4cmb/tod.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/s4cmb/tod.py b/s4cmb/tod.py index e84d586..979d7a1 100755 --- a/s4cmb/tod.py +++ b/s4cmb/tod.py @@ -1067,7 +1067,7 @@ def map2tod(self, ch): pol_ang_pair, pol_ang2_pair = self.compute_simpolangle( ch, pa_pair, polangle_err=False ) - except ValueError: + except UnboundLocalError: pol_ang_pair = pol_ang pol_ang2_pair = pol_ang2 From 35a7263347901d59b04807b6409fbec1fae591aa Mon Sep 17 00:00:00 2001 From: Julien Date: Tue, 26 Jan 2021 11:36:33 +0100 Subject: [PATCH 15/15] Add a new workflow to lint the code --- .github/workflows/lint-s4cmb.yml | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 .github/workflows/lint-s4cmb.yml diff --git a/.github/workflows/lint-s4cmb.yml b/.github/workflows/lint-s4cmb.yml new file mode 100644 index 0000000..7b8a587 --- /dev/null +++ b/.github/workflows/lint-s4cmb.yml @@ -0,0 +1,25 @@ +name: PEP8 + +on: [push, pull_request] + +jobs: + build: + + runs-on: ubuntu-latest + strategy: + matrix: + python-version: [3.6] + + steps: + - uses: actions/checkout@v2 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v2 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install flake8 + - name: Analysing the code with flake8 + run: | + flake8 s4cmb/*.py --count --show-source --statistics --max-line-length 86 --doctests --ignore=E741,W605