cbasspy package¶
Subpackages¶
- cbasspy.constantsPlotting package
- cbasspy.groundpy package
- cbasspy.jackknives package
- cbasspy.mapManipulations package
- Submodules
- cbasspy.mapManipulations.beam2bl module
- cbasspy.mapManipulations.beam2bl_I module
- cbasspy.mapManipulations.beam2bl_P module
- cbasspy.mapManipulations.deconvolveMap module
- cbasspy.mapManipulations.integrate module
- cbasspy.mapManipulations.pix_QU_rot module
- cbasspy.mapManipulations.rotate_QU module
- cbasspy.mapManipulations.wl_anafast module
- cbasspy.pixelSimsData package
Submodules¶
cbasspy.angles module¶
Common angle conversions
-
latLon2ThetaPhi
(lat, lon, deg=False)¶ Convert latitude and longitude to colatitude (theta) and longitude (phi) Default angles all in radians
Parameters: lat : array_like
Latitude
lon : array_like
Longitude
deg : optional, {True,False}
Whether angles in deg or rad
-
raDec2GlatGlon
(RA, DEC, deg=False, EPOCH=36525.0)¶ Convert RA and DEC to galactic lat and lon. Default angles in radians Default epoch J2000
Parameters: RA : float
Right ascension
DEC : float
Declination
deg : optional, {True,False}
Whether angles in deg or rad
EPOCH : ephem object, optional {ephem.J2000}
Epoch of conversions
Returns: glat : float
Galactic latitude
glon : float
Galactic longitude
cbasspy.bin_llcl module¶
Function written by PolSPICE
-
bin_llcl
(llcl_in, bin, flatten=False, uniform=False)¶ - turns llcl_in (= continuous l*(l+1)*cl/2Pi) into a binned version with
- a constant or variable binwidth ‘bin’
Parameters: llcl_in : array_like
input l*(l+1)*Cl/2Pi, 1D vector, defined for each l from l=0
bin : array_like,
can be either a scalar = dl or a vector defining the bins boundaries : [low0, low1, low2, …,low(n-2), low(n-1)+1]
flatten : optional, {False,True}
If set, the input C(l) is multiplied internally by l*(l+1)/2Pi before being binned. By default, the input C(l) is binned as is.
uniform : optional, {False, True}
if set, each l is given the same weight in the bin. By default, a weighting propto (2*l+1) (inverse cosmic variance) is applied to each l. In any cases, the output x is the same
Returns: x : array_like
center of bins
y : array_like
binned l*(l+1)*Cl/2Pi
dx : float (?)
width of each bin
dy : float
- returns on output the rms of C(l) for a full sky observation
= C(l) * sqrt( 2/ 2l+1 / dl)
cbasspy.clusteringCode module¶
Python module to cluster pixels in the sky into regions.
-
compute_naive_betas
(map1_name, map2_name, freq1, freq2, offset1=0.0, offset2=0.0, make_P_from_inputs=False)¶ Compute the simple pixel-by-pixel spectral index. Input maps should be in the same coordinate system, have the same Nside and be in the same units.
Parameters: map1_name : str
Name of first map file
map2_name : str
Name of second map file
freq1 : float
Frequency of first map [GHz]
freq2 : float
Frequency of second map [GHz]
offset1 : float, 0.
Zero level offset to add to map1
offset2 : float, 0.
Zero level offset to add to map2
Returns: beta : array_like
Healpix map of spectral indices
-
compute_pixel_vectors
(Nside)¶ Computes the vectors to the pixel centres for every pixel in a map.
Returns: vec : array_like
Array of vectors to the pixel centres
-
find_meanshift_clusters
(beta12, beta23, vec, goodIdx, bandwidth=0.13, min_bin_freq=20)¶ Assign each pixel to a region using the meanshift algorithm.
Parameters: X : array_like
Array of values to cluster on
goodIdx : array_like
The indices to use
bandwidth : float, 0.13
Bandwidth to pass to clustering algorith
min_bin_freq : int, 20
Minimum number of pixels in each region
Returns: labels_map : array_like
Healpix map, with each pixel assigned a region.
-
find_voronoi_clusters
(S_map_name, N_map=None, S_map_offset=0, threshold_score=2.5, min_pixels=20)¶ Assign each pixel to a region using the Voronoi binning algorithm.
Parameters: S_map_name : str
Name of input signal map
N_map : str, None
Name of input noise map, if None then noise is set to unity.
S_map_offset : float, 0
Zero level offset to add to signal map
threshold_score : float, 2.5
Target S/N for each region
min_pixels : int, 20
Target minimum number of pixels for each region
Returns: voronoi_map : array_like
Healpix map, with each pixel assigned a region.
distance_map : array_like
Healpix map, with the distance from each pixel to its region centroid.
cbasspy.const module¶
-
TCMB
= <cbasspy.const.constant object at 0xdaa5450>¶ CMB mean temperature
2.7255 K
-
ap_eff
= <cbasspy.const.constant object at 0xdaa5510>¶ Aperture efficiency
0.55
-
area_eff
= <cbasspy.const.constant object at 0xdaa5590>¶ Effective dish area
ap_eff.value*dish_area.value m**2
-
c
= <cbasspy.const.constant object at 0xdaa5490>¶ Speed of light in vacuum
2.99792458e8 m s**-1
-
class
constant
(value, unit_string, name)¶ Bases:
object
A numerical constant
Methods
-
displayValue
()¶ Print the value of the constant and it’s units
-
-
dish_area
= <cbasspy.const.constant object at 0xdaa5550>¶ Dish area
(6.1/2)**2*np.pi m**2
-
h
= <cbasspy.const.constant object at 0xdaa53d0>¶ Planck constant.
6.62607004e-34 m**2 kg s**-1
-
kb
= <cbasspy.const.constant object at 0xdaa5410>¶ Boltzmann constant
1.38064852e-23 m**2 kg s**-2 K**-1
-
sample_freq
= <cbasspy.const.constant object at 0xdaa54d0>¶ Telescope sampling frequency
- Hz
cbasspy.conv module¶
Module containing functions to convert between brightness, temperature, SI units and Jy
-
brightness2Trj
(B, nu, SI=True)¶ Convert from brightness to brightness temperature in K_RJ.
Brightness temperature is also known as Reyleigh Jeans temperature and sometimes Antenna temperature By default B is assumed to be in SI units, set SI=False for Jy/sr nu in Hz.
Parameters: B : float
Brightness, in SI units by default.
nu : float
Frequency, in Hz
SI : {True, False}
If True, assumes brightness in SI units. If False, assumes brightness in Jy/sr.
Returns: T_RJ : Float
Rayleight-Jeans brightness temperature in Kelvin.
See also
trj2Brightness
- Inverse conversion
-
fluxDensityJy2Si
(F_Jy)¶ Converts a flux density from Jy to SI units W m^{-2} Hz^{-1}
Parameters: F_Jy : Float
Flux density in Jy.
Returns: F_SI : Float
Flux density in W m^{-2} Hz^{-1}
See also
fluxDensitySi2Jy
- Inverse conversion
-
fluxDensitySi2Jy
(F_SI)¶ Converts a flux density from SI units W m^{-2} Hz^{-1} to Jy
Parameters: F_SI : Float
Flux density in W m^{-2} Hz^{-1}
Returns: F_Jy : Float
Flux density in Jy.
See also
fluxDensityJy2Si
- Inverse conversion
-
tcmb2Trj
(T_CMB, nu)¶ Convert from CMB thermodynamic temperature, K_CMB to brightness temperature in K_RJ.
Sometimes K_CMB is called thermodynamic temperature but this could be confusing nu in Hz.
Parameters: T_CMB : Float
Thermodynamic or CMB temperature in Kelvin.
nu : Float
Frequency in Hz.
Returns: T_RJ : Float
Rayleigh-Jeans brightness temperature in Kelvin.
See also
trj2Tcmb
- Inverse conversion
-
trj2Brightness
(T_RJ, nu, SI=True)¶ Convert from brightness temperature in K_RJ to brightness.
By default returns brightness in SI units, set SI=False for Jy/sr Brightness temperature also known as Reyleigh Jeans temperature and sometimes Antenna temperature By default B is returned in SI units, set SI=False for Jy/sr nu in Hz.
Parameters: T_RJ : float
Rayleight-Jeans brightness temperature in Kelvin.
nu : float
Frequency, in Hz
SI : {True, False}
If True, returns brightness in SI units. If False, returns brightness in Jy/sr.
Returns: B : Float
Brightness, in SI units by default.
See also
brightness2Trj
- Inverse conversion
-
trj2Tcmb
(T_RJ, nu)¶ Convert from brightness temperature in K_RJ to K_CMB.
Sometimes K_CMB is called thermodynamic temperature but this could be confusing. nu in Hz.
Parameters: T_RJ : Float
Rayleigh-Jeans brightness temperature in Kelvin.
nu : Float
Frequency in Hz.
Returns: T_CMB : Float
CMB or thermodynamic temperature in Kelvin.
See also
tcmb2Trj
- Inverse conversion.
cbasspy.createClusters module¶
cbasspy.diffuseComponentModels module¶
Python module to simulate the frequency spectra of various emission mechanisms.
-
add_white_noise
(T, rms)¶ Adds white noise to temperatures
Parameters: T : array_like
signal to add noise to
rms : array_like
sigma of gausian from which to draw noise
Returns: Tsim : array_like
T+noise realization
-
cmb_model
(amp, T0, f)¶ CMB model
T_RJ = amp * x^2*exp(x)/(exp(x)-1^2 where x = h*nu/(k*T0)
Parameters: amp : array_like
amplitude of CMB fluctuation [K_CMB].
T_0 : array_like
Temperature of CMB [K]
f : array_like
frequency at which to return amplitude [Hz]
Returns: T_RJ : array_like
[K_RJ] at frequency f.
-
curved_power_law_model
(amp, beta, C, f0, f, factorHalf=False)¶ Curved Power Law model
- if factorHalf==False:
- T_RJ = amp * (nu/nu0)^(beta+C*log(f/f0))
- elif factorHalf==True:
- T_RJ = amp * (nu/nu0)^(beta+0.5*C*log(f/f0))
Parameters: amp : array_like
amplitude at reference frequency
beta : array_like
spectral index
C : array_like
Curvature
f0 : array_like
reference frequency [Hz]
f : array_like
frequency at which to return amplitude [Hz]
factorHalf : {False, True}
Whether to multiply C by 0.5.
Returns: T_RJ : array_like
Amplitude at frequency f in same units as amp.
-
free_free_model
(Te, EM, f)¶ Free-free model
Two parameter model of Draine - put in a reference!
Parameters: Te : array_like
Effective electron temperature in beam [K]
EM : array_like
Effective emission measure in beam [cm^6 pc^-1]
f : array_like
Frequency at which to return amplitude [Hz]
Returns: T_RJ : array_like
[K] at frequency f.
-
interpolater
= <scipy.interpolate.interpolate.interp1d object at 0xdb8d470>¶ scipy interp1d function for the SPDUST2 model
-
power_law_model
(amp, beta, f0, f)¶ Power Law model
T_RJ = amp * (nu/nu0)^beta
Parameters: amp : array_like
amplitude at reference frequency
beta : array_like
spectral index
f0 : array_like
reference frequency [Hz]
f : array_like
frequency at which to return amplitude [Hz]
Returns: T_RJ : array_like
Amplitude at frequency f in same units as amp.
-
spdust2_model
(fPeak, amp, f0, f)¶ Spinning dust model
SPDUST2 in a cold neutral medium shifted in frequency and amplitude space. -Add a reference-
Parameters: fPeak : float
Frequency of peak emission [Hz]
amp : float
Amplitude at reference frequency f0 [K]
f : array_like
Freq at which to return amplitude [Hz]
Returns: T_RJ : array_like
[K] at frequency f.
-
synch_moment_expansion
(amp, beta, moments, f0, f)¶ Created – Richard David Paradorn Grumitt (09/05/2018)
Adapted by Jaz (22/05/2018)
This model implements the moment expansion SED for a power law superposition described in J. Chluba et. al. 2017. Model employed here assumes power law amplitudes and spectral indices are not correlated.
Adapted from the PySM synchrotron class method, which can be used to generate whole sky simulations with instrumental and noise effects.
Parameters: amp : float
Amplitude
beta : float
spectral index
moments : array_like
moments of the synchrotron SED up to desired order, as defined in J. Chluba et. al. 2017. – numpy.ndarray
f0 : float
Reference frequency in Hz
f : array_like
frequency at which to evualate spectrum
Returns: T_RJ : array_like
Power law moment expansion SED for I, Q and U, evaluated at frequency nu
-
thermal_dust_model
(amp, beta, T, f0, f)¶ Thermal dust model
T_RJ = amp*(f/f0)^(beta+1)*(exp(gamma*f0)-1)/(exp(gamma*f)-1.) where gamma = h/(kb*T)
Parameters: amp : array_like
amplitude of dust [X] at reference frequency
beta : array_like
emissivity index
T : array_like
Dust temperature [K]
f0 :array_like
Reference frequency [Hz]
f : array_like
frequency at which to return amplitude [Hz]
Returns: T_RJ : array_like
[X] at frequency f.
cbasspy.fileListBash module¶
A slightly ugly set of wrappers that call commonly useful bash commands for dealing with file lists.
-
createDailyFileLists
(start_date_str, number_of_days, fits_dir, suffix='')¶ Create daily file lists from all the fits files in one directory.
Outputs files as DD-MMM-YYYY_files[suffix].txt
Parameters: start_date_str : str
start date in format DD-MMM-YYYY, e.g. ‘10-Sep-2010’
number_of_days : int
Number of subsequent days to create file lists for, e.g. 20
fits_dir : str
Directory containing TOD .fits files
suffix : str, optional {‘’}
extra suffix to add to output file names
-
createMonthlyFileLists
(start_date_str, number_of_months, fits_dir, suffix='')¶ Create monthly file lists from all the fits files in one directory.
Outputs files as MMM-YYYY_files[suffix].txt
Parameters: start_date_str : str
start date in format MMM-YYYY, e.g. Sep-2010
number_of_months : int
Number of subsequent months to create file lists for, e.g. 20
fits_dir : str
Directory containing TOD .fits files
suffix : str, optional {‘’}
extra suffix to add to output file names
-
excludeDateRange
(file_list_name, output_name, exclusion_min_str, exclusion_max_str)¶ Exclude a date range from a file list.
Parameters: file_list_name : str
Name of file list
output_name : str
Where to save the outputs. Don’t include the file extension .txt, the code adds it.
exclusion_min_str : str
Earliest date to exclude, format e.g. 01-Aug-2014:00:00:01
exclusion_max_str : str
Latest date to exclude, formate e.g. 01-Oct-2014:00:00:01
-
excludeLinesWithMatches
(fileListName, matchesListName, outListName, verbose=False)¶ Function to remove files from fileListName that contain any string in the file matchesListName.
Parameters: file_list_name : str
Name of file list
matchesListName : str
Name of list of matches
outListName: str
What to save the output to
verbose : bool, False
Whether to print messages to screen.
-
makeRaDecList
(fileList, RAmin, RAmax, DECmin, DECmax, SAMPLESmin, verbose=0)¶ Give functiona a list of files and an RA/DEC region (in radians). Fits files with more than SAMPLESmin number of samples within the region are outputted in a goodFileList.
Parameters: fileList : list
list of file names
RAmin : float
Minimum RA to accept [radians]
RAmax : float
Maximum RA to accept [radians]
DECmin : float
Minimum DEC to accept [radians]
DECmax : float
Maximum DEC to accept [radians]
SAMPLESmin : int
Minimum number of samples within the region for the file to be passed to the good list.
verbose : optional {0}
If >=1 then print summary to screen.
Returns: goodFileList : list
List of files that have met the criteria
-
openFileList
(fileListName)¶ Open the .txt file fileListName. Return a list of strings, each element a line from the file
This is only useful if you can’t remember how to open a txt file with python!
Parameters: fileListName : str
Name of file to open.
Returns: fileList : list
List of file names
See also
-
prependString
(string2prepend, inFileName, outFileName, verbose=1)¶ prependString will prepend the string string2prepend to the start of each line of a file in inFileName. The output is saved to outFileName. If verbose>0 then a summary will be printed to the screen.
Parameters: string2prepend : str
String to prepend to each line of input file
inFileName : str
Name (including path) of file to append string to
outFileName : str
Name (including path) of file to save output to
verbose : optional {1,0}
If >1 then a summary will be printed to screen.
-
saveCommonLines
(inFileName1, inFileName2, outFileName, verbose=1.0)¶ saveCommonLines looks for matching lines in inFileName1 and inFileName2. Matching lines will be written to outFileName. If verbose>0 then a summary will be printed to the screen.
Parameters: inFileName1 : str
Name (including path) of first file
inFileName2 : str
Name (including path) of second file
outFileName : str
Name (including path) of file to save output to
verbose : optional {1,0}
If >1 then a summary will be printed to screen.
-
saveFileList
(outListName, fileList)¶ Save a list of files to a txt file called outListName
This is only useful it you can’t remember how to save a txt file with python!
Parameters: outListName : str
Where to save the file list
fileList : list
list of lines to save
See also
-
splitFileListCentralDate
(file_list_name, output_name)¶ Split a file list into two (potentially) unequal halves. The central date, halfway between the first and last dates in the file, is the boundary between files that are outputted into one file and another.
Parameters: file_list_name : str
Name of file list
output_name : str
Where to save the outputs. Don’t include the file extension .txt, the code adds it.
-
splitFileListSortDate
(file_list_name, output_name, number_of_splits=2)¶ Split a file list into equal sections. Each containing 1/number_of_splits of the chronologically oredered files.
Parameters: file_list_name : str
Name of file list
output_name : str
Where to save the outputs. Don’t include the file extension .txt, the code adds it.
number_of_splits : int, 2
Number of splits to make to the data
-
splitOddAndEvenLines
(inFileName, evenFileName, oddFileName, verbose=1)¶ Split a file list into two files, one listing even lines and the other listing odd lines.
Parameters: inFileName : str
File to split
evenFileName : str
where to save the even lines
oddFileName : str
where to save the odd lines
verbose : optional {1}
If >1 print a summary to screen
cbasspy.fit1line module¶
-
TruncatedOneOverX_like
(value, a, b)¶ Depricated. In future this will be replaced by
cbasspy.pymcCustomPriors.BoundedOneOverX_like()
.
-
fit1line
(x_obs, y_obs, x_sigma, y_sigma, init_slope_1=0.1, niter=500000.0, nburn=100000.0, nthin=20, ntune=500.0, seed=1234, rangeLower_x=1e-06, rangeUpper_x=2000.0, rangeLower_y=1e-06, rangeUpper_y=2000.0, locLower_x=-2000.0, locUpper_x=1000.0, locLower_y=-2000.0, locUpper_y=1000.0, save2db=False, dbName='mcmcChains_fit1line.db3')¶ Fit a straight line with outliers to data with errors in x and y. The problem is parametrized according to Gull 1989, this method produces unbiased estimates of the parameters of a straight line. It assumes errors on x_obs and y_obs are Gaussian. The line is reparametrized as hat{Y}=pm hat{X}, hat{X}=(x-x0)/Rx and hat{Y}=(y-y0)/Ry the loc parameters x0 and y0 are given uniform priors, the range parameters Rx and Ry are given 1/R priors. The true locations of the x values are marginalised over. (See Gull 1989 for details)
The outlies are modelled using a mixtude model; y = a*x+b OR 2D Gaussian with no covariance. The 2D Gaussian is paramterized by mean (x0_3, y0_3) and covariance ((Rx_3^2,0),(0,Ry_3^2))
Parameters: x_obs : array_like
x values
y_obs : array_like
y values
x_sigma : array_like
errorbar on x
y_sigma : array_like
errorbar on y
init_slope_1 : float {0.1}
Depricated, kept for back compatablility. Parameter was originally Starting guess for slope of line, now does nothing.
niter : int, optional {5e5}
Number of steps in MCMC chain
nburn : int, optional {1e5}
Number of steps to burn
nthin : int, optional {20}
Factor by which to thin the chain
ntune : int, optional {500}
Number of steps between tuning of step size (only during burn in)
seed : int, optional {1234}
Random seed
rangeLower_x : float, optional {1e-6}
Lower limit of x range parameter [Must be >0]
rangeUpper_x : float, optional {2e3}
Upper limit of x range parameter [Must be >0]
rangeLower_y : float, optional {1e-6}
Lower limit of y range parameter [Must be >0]
rangeUpper_y : float, optional {2e3}
Upper limit of y range parameter [Must be >0]
locLower_x : float, optional {-2e3}
Lower limit of x location parameter [float]
locUpper_x : float, optional {1e3}
Upper limit of x location parameter [float]
locLower_y : float, optional {-2e3}
Lower limit of y location parameter [float]
locUpper_y : float, optional {1e3}
Upper limit of y location parameter [float]
save2db : boolean {False}
Whether to save output as a sqlite database
dbName : str {‘mcmcChains_fit1line.db3’}
Name to save the database as
Returns: mcmc : pymc object
The function returns a PyMC object that contains the traces of
- The free parameters in the model,
- The straight line component x0, y0, Rx, Ry
- model_a_1 and model_b_1 are the calculated slope and y-intercept for each step in the chain
- The outlier component x0_3, y0_3, Rx_3, Ry_3, frac_3
- (frac_3 is the fraction of points that belong to the outlier component)
- True x values x
To access these chains,:
mcmc.trace('x0_3')[:]
Notes
- The lower and upper limits on the range and loc parameters should be set far away from the expected parameter range and central values and decided
on before looking at the data. * If the input parameters result in zero probability, None is returned.
cbasspy.fit_mean_var module¶
-
intrinsic_scatter
(obs, obs_err, mean_lower, mean_upper, sd_lower, sd_upper, nsamples=2000, ntune=1000, accept_prob=0.9, jeffreys_prior=False)¶ Function for fitting the mean and additional intrinsic scattter for a dataset.
Returns: trace: pymc3.MultiTrace
Output trace from your fitting - pymc3.MultiTrace
cbasspy.fitexey module¶
This module may be incorrect, it definitely needs checking.
-
fitexey
(x_obs, y_obs, x_sigma, y_sigma, niter=500000.0, nburn=10000.0, nthin=20, ntune=100.0, dbname='./model_db.pickle', seed=1234)¶ Fit a straight line to data with errors in x and y. y = ax+b Returns traces of the slope, a, and intercept, b. This method is unbiased - Gull 1989. x_obs: x values y_obs: y values x_sigma: errorbar on x y_sigma: errorbar on y niter: Number of steps in MCMC chain nburn: Number of steps to burn nthin: Factor by which to thin the chain ntune: Number of steps between tuning of step size (only during burn in) dbname: Not currently used seed: Random seed
cbasspy.masks module¶
-
produce_dec_mask
(NSIDE, rotTo, cutMin, cutMax)¶ Return a declination mask, masking region between cutMin and cutMax Pixels are set to 0 within the mask and 1 outside.
Parameters: NSIDE : int
Nside of mask, e.g. NSIDE=512
rotTo : str, {‘C’,’G’,’E’}
Coordinate system to rotate mask to.
cutMin : float
Lower declination bound, e.g. cutMin=-90.
cutMax : float
Upper declination bound, e.g. cutMax=0.
Returns: fullR : array_like
HEALPix map.
-
produce_frac_mask
(inMap, frac, account_for_masked=False)¶ Return a mask excluding the brightest fraction of the inMap. The masked pixels are set to 0, unmasked pixels to 1.
Parameters: inMap : array_like
HEALPix sky map
frac : float
Fraction of the sky to mask.
account_for_masked : boolean, False
frac of unmasked pixels if True, else frac of all pixels
Returns: mask : array_like
HEALPix map.
-
produce_plane_mask
(NSIDE, rotTo, cutOff)¶ Return a mask of the galactic plane. Pixels are set to 0 if too close to the plane and 1 if far enough away.
Parameters: NSIDE : int
Nside of mask, e.g. NSIDE=512
rotTo : str, {‘C’,’G’,’E’}
Coordinate system to rotate mask to.
cutOff : float
cutOff in degrees from galactic plane.
Returns: fullR : array_like
HEALPix map.
cbasspy.noise module¶
-
fitNoiseModelToTodFast
(freq, I_N, niter=5000.0, verbose=True)¶ Fits a white + red noise model to a periodogram using a normal approximation to the posterior. Draws posterior from samples drawn from the normal approx.
Parameters: freq : array_like
Frequencies [Hz]
I_N : array_like
Periodogram
niter : int
Number of steps in chain
Returns: Mean[alpha] : float
Mean[sigma_r] : float
Mean[sigma_w] : float
Cov[alpha,alpha] : float
Cov[alpha,sigma_r] : float
Cov[alpha,sigma_w] : float
Cov[sigma_r,sigma_r] : float
Cov[sigma_r,sigma_w] : float
Cov[sigma_w,sigma_w] : float
See also
fitNoiseModelToTodFitsFile
- Does full MCMC
-
fitNoiseModelToTodFitsFile
(freq, I_N, niter=50000.0, nburn=25000.0, thinFactor=1, ntune=1000.0, verbose=False)¶ Fits a white + red noise model to a periodogram using the full Whittle Pseudo-Likelihood. Method described in Luke Jew’s Thesis, Chapter 3.
Parameters: Data parameters
^^^^^^^^^^^^^^^
freq : array_like
Frequency in Hz
I_N : array_like
Periodogram
MCMC parameters
^^^^^^^^^^^^^^^
niter : int
Number of steps in chain
nburn : int
Number of initial steps to burn
thinFactor : int
Thinning factor
ntune : int
How often to tune step size
Returns: Mean[alpha] : float
Max posterior alpha : float
Mean[sigma_r] : float
Max posterior sigma_r : float
Mean[sigma_w] : float
Max posterior sigma_w : float
Mean[f_knee] : float
Max posterior f_knee : float
Var[f_knee] : float
Cov[alpha,alpha] : float
Cov[alpha,sigma_r] : float
Cov[alpha,sigma_w] : float
Cov[sigma_r,sigma_r] : float
Cov[sigma_r,sigma_w] : float
Cov[sigma_w,sigma_w] : float
chi2 : float
dof : int
See also
ln_noise_model
- Noise model that is fit to the data.
cbasspy.const
- Fundamental constants
cbasspy.angles
- Angle conversions to subtract sky signal.
-
leastsq_fitNoiseModelToTodFitsFile
(hdulist, chanName, model_sky_map, sample_freq=100.0, minChunkLength=100, savetxt=False, savetxtStub='example')¶ Least squares fitting method. Needs to average lots of chunks together to estimate P_xx with normal error bars
p[0] = alpha p[1] = f_knee p[2] = sigma_wParameters: hdulist : pyfits object
TOD Fits file
chanName : {‘I1’,’I2’,’Q1’,’Q2’,’U1’,’U2’}
Which channel to fit to
model_sky_map : array_like
HEALPix model sky map
sample_freq : float, optional {const.sample_freq.value}
Sample frequency of TOD
minChunkLength : int, optional {100}
Minimum chunk length to fit to
savetxt = {False,True}
Whether to save outputs
savetxtStub : str, optional {‘example’}
Name to attach to output files.
Returns: p : array_like
Best fitting parameters, alpha, f_knee, sigma_w
mesg : str
Fitting string
flag : int
Notes
Depricated
-
ln_noise_model
(p, f, f_0=0.1)¶ Returns the natural logarithm of a white + 1/f noise model power spectrum.
- ln_P_xx = ln_noise_model(p,f,f_0=0.1)
- = np.log(p[2]**2 + p[1]**2*np.power(f/f_0,-p[0])))
Parameters: p : array_like
Noise parameters shape = [3,]. p[0] = alpha, p[1] = sigma_r, p[2] = sigma_w
f : array_like
Frequency in Hz
Returns: ln_P_xx : array_like
Noise power spectrum.
-
minimizeParams
(f, ln_P_xx, p_guess)¶ Parameters needed to minimise the rms value of the difference between ln_P_xx and the power spectrum computed by ln_noise_model. residsFunc is used to calculate the difference.
Parameters: f : array_like
Frequency in Hz
ln_P_xx : array_like
Natural logarithm of power spectrum.
p_guess : array_like
The starting values of the parameters, shape = [3,]. p_guess[0] = alpha p_guess[1] = sigma_r p_guess[2] = sigma_w
See also
residsFunc
- The residuals function.
-
old_ln_noise_model
(p, f)¶ Returns the natural logarithm of a white + 1/f noise model power spectrum.
Please note, this function parametrises the noise model with f_knee - not sigma_r ln_P_xx = ln_noise_model(p,f)
= np.log(p[2]**2 + p[1]**2*np.power(f,-p[0])))Parameters: p : array_like
Noise parameters shape = [3,]. p[0] = alpha, p[1] = f_knee, p[2] = sigma_w
- f : array_like
Frequeny in Hz.
Returns: ln_P_xx : array_like
Noise power spectrum
-
residsFunc
(p, f, ln_P_xx)¶ Returns the difference between ln_P_xx and the power spectrum computed by ln_noise_model(p,f) with params p.
p[0] = alpha p[1] = sigma_r p[2] = sigma_wParameters: p : array_like
Noise parameters
f : array_like
Frequency in Hz
ln_P_xx : array_like
Natural logarithm of power spectrum to compare to.
Returns: diff : array_like
Difference
See also
ln_noise_model
- Function called
cbasspy.noise_fit_12 module¶
-
calcSpectralContaminantChisq
(f, I_N, ps, freqLower, freqUpper)¶
-
noise_fit_12
(outputFileName, fileListName, model_sky_name, convert2Kelvin=1.0, doPol=1.0, niter=20000, nburn=10000, thinFactor=1, ntune=500, cutName='NIGHT', spikeThreshold=0.035, mpi_on=True)¶ Calculates noise stats for I1/2,Q1/2 and U1/2 Previous versions averaged these to I, Q and U. Outputs saved to a sqlite database with filename outputFileName
Parameters: outputFileName : str
Name of db3 file to create and save results in. If outputFileName already exists then a table called noiseStats is created. If noiseStats already exists then it is overwritten
fileListName : str
Name of txt file that lists the TOD fits files
model_sky_name : str
Name of healpix map .fits file that is used to generate an expected sky signal. MAP MUST BE IN GALACTIC COORDINATES
convert2Kelvin : float, (1.)
Scaling factor by which the TOD and map are multiplied by
doPol : int, {1,0}
if 0 only I, if 1 IQU
niter : int
number of steps to take in MCMC chains
nburn : int
number of steps to burn
thinFactor : int
Factor by which to thin the chains
ntune : int
how frequently to tune the step proposal distribution during tuning
mpi_on : boolean, {(True),False}
Whether this should be run on multiple cores.
cutName : str
Name of extension in fits file, default ‘NIGHT’
spikeThreshold : float, (0.035)
If model sky changes by >spikeThreshold between 2 samples then don’t use that chunk. Current default is 0.035, it used to be 0.025
Returns: nothing :
cbasspy.overwriteFitsNoiseStats module¶
-
ASrSw2Fk
(alpha, sigma_r, sigma_w, f_0=0.1)¶ Take alpha, sigma_r, sigma_w and f_0 (default 0.1Hz) Return f_knee
-
findFitsExtNo
(hdulist, cutName)¶ Find the header extension with name cutName.
Parameters: hdulist : dict
Opened fits file
cutName : str
Name of extension
Returns: extNo : int
hdulist[extNo] is the table of interest
-
overwriteFitsNoiseStats_fromSql
(fileListName, inSqlFileName, cutName='NIGHT')¶ Function to overwrite the noise statistics in the fits files listed in fileList with the noise properties in the SQL database.
Parameters: fileListName : str
Name of the .txt file containing the file list
inSqlFileName : str
Name of the SQL databse
cutName : str, ‘NIGHT’
Which cut extension to overwrite
Returns: filesWithoutNoiseStats : array_like
A list of files where the default noise properties were written.
filesWithNoiseStats : array_like
A list of files where a match was found in the database.
-
wavg_12_calc
(inSqlFileName, outSqlFileName)¶ Calculate the weighted average noise statistics for each fits file. Takes a sqlite database as input. Outputs results in a second sqlite database.
Parameters: inSqlFileName : str
Name of input database, most likely created by
cbasspy.noise_fit_12.noise_fit_12()
.outSqlFileName : str
Name of output database.
-
weightedNoiseStatsMean
(aa, ar, aw, rr, rw, ww, alpha_mu, sigma_r_mu, sigma_w_mu, error=False)¶ Calculate the weighted average noise properties, accounting for parameter covariance
Parameters: aa : array_like
Var(alpha)
ar : array_like
Cov(alpha,sigma_r)
aw : array_like
Cov(alpha,sigma_w)
rr : array_like
Var(sigma_r)
rw : array_like
Cov(sigma_r,sigma_w)
ww : array_like
Var(sigma_w)
alpha_mu : array_like
array of alpha values
sigma_r_mu : array_like
array of sigma_r values
sigma_w_mu : array_like
array of sigma_w values
error : Boolean {False}
If False just return weighted average, if True also return covariance.
Returns: means : array_like
Array of parameter averages (only if error==True)
uncertainty : array_like
Array of uncertainties (only if error==True)
alpha_avg : float
Weighted average alpha (only if error==False)
sigma_r_avg : float
Weighted average sigma_r (only if error==False)
sigma_w_avg : float
Weighted average sigma_w (only if error==False)
cbasspy.pymcCustomPriors module¶
-
BoundedNegOneOverX_like
(value, lower, upper)¶
-
BoundedOneOverX_like
(value, lower, upper)¶ - p(X)~1/X when lower<X<upper
- 0 else
Returns log probability
-
JeffCuPrior_like
(value, sigs, nus, nu0, lower, upper, beta)¶
-
JeffEmPrior_like
(value, sigs, nus, lower, upper, Te)¶
-
JeffNupPrior_like
(value, sigs, nus, f0, Y, lower, upper)¶
-
JeffSiPrior_like
(value, sigs, nus, nu0, lower, upper, G)¶ Jeffreys prior for spectral indices when likelihood is L~1/sigma * exp(-(X-A*G*(nu/nu0)^beta)/2sigma^2) Where A is an amplitude and G is a fixed function of frequency e.g. for power law with no curvature. G = 1.0 e.g. for a power law with curvature, G = exp(0.5*C*ln(nu_i/nu_0)) e.g. for modified black body (i.e. thermal dust) G = (nu_i/nu_0) * [exp(h*nu_0/k*T)-1]/[exp(h*nu_i/k*T)-1]
Parameters: sigs : array_like
array of uncertanties
nus: array_like
array of frequencies
nu0 : float
reference frequency
lower : float
lower bound on spectral index
upper : float
upper bound on spectral index
value :float
proposed value of spectral index
G : array_like
either float or an array of floats.
Returns: ln(probability) : float
-
JeffSpectralPrior_like
(value, sigs, nus, nu0, lower, upper, Y, X, factor)¶
-
JeffTdPrior_like
(value, sigs, nus, nu0, lower, upper, Y)¶
-
TruncatedExpX_like
(value, lower, upper)¶
-
TruncatedX_like
(value, lower, upper)¶
cbasspy.scalePolCovsBy2 module¶
-
scalePolCovsBy2
(cov_name, output_str='_polCovCor')¶ Scale the Q and U covariances by 2. Descart thinks that Q1 and Q2 (U1 and U2) are independent - they are not. If you make a map using both 1 and 2 channels then you will need to apply this scaling.
Parameters: cov_name : str
Name of covariance fits file
output_str : str, ‘_polCovCor’
String to attached to output file name. The output file name takes the form:
'{0}{1}.fits'.format(cov_name[0:-5],output_str)
Returns: cov_map : array_like
Corrected healpix map.
cbasspy.simulate module¶
-
groundFromTemplate
(az, templateFileName)¶ Simulates the ground signal at a given azimuth and elevation by taking a sample ground template.
Parameters: az : array_like
azimuths to sample at [degrees]
templateFileName : str
Template filename (and directory).
Returns: signal : array_like
Simulated ground TOD
Notes
LJ - 05/07/18 AJB - 05/07/18
-
overwrite_tod_fits_file
(fitsFile, whichField, values)¶ Overwrites the field [whichField] in extension 1 of the TOD Fits file [fitsFile] with the values in [values].
Parameters: fitsFile : str
File name
whichField : str
Field name
values : float, array_like
values to overwrite field with
-
polSkyFromMap
(RA, DEC, MJD, LAT, LON, Q_sim_map, U_sim_map)¶ Simulate stokes Q and U signals from input Q and U maps ~~in Celestial coordinates~~.
Future work Support input maps in other coordinate systems (e.g. Galactic) Support angles in degrees
Parameters: RA : array_like
Right ascension in radians
DEC : array_like
Declination in radians
MJD : array_like
MJD
LAT : float
Latitude of telescope, in radians
LON : float
Longitude of telescope, in radians
Q_sim_map : array_like
Simulated Q Healpix sky map, in Celestial coordinates
U_sim_map : array_like
Simulated U Healpix sky map, in Celestial coordinates
Returns: Q_a : array_like
Rotated Q signal in telescope coordinates
U_a : array_like
Rotated U signal in telescope coordinates
-
skyFromMap
(glat, glon, sky_sim_map)¶ Simulate sky-signal TOD from the HEALPix map sky_sim_map (in Galactic coordiantes) at coordinates glat and glon (in radians).
Parameters: glat : array_like
array of Galactic latitudes [radians]
glon : array_like
array of Galactic longitudes [radians]
sky_sim_map : array_like
HEALPix sky map.
Returns: signal : array_like
Simulated sky signal
-
whiteAndRedNoiseSim
(sigma_w, f_knee, alpha, sample_freq, N, return_pdf=False)¶ nosie (,f,P) = whiteAndRedNoiseSim(sigma_w,f_knee,alpha,sample_freq,N,return_pdf=False) Simulate red and white noise in TOD.
Parameters: sigma_w : float
white noise level in K
f_knee : float
knee frequency in Hz
alpha : float
power law
sample_freq : float
for C-BASS this is 100 Hz
N : int
length of noise to simulate
Returns: noise : array_like
Simulated noise signal
cbasspy.smart_ud_grade module¶
-
alm_ud_grade
(map_in, nside_out)¶ Upgrade or downgrade map in alm space.
Parameters: map_in : array_like
input map
nside_out : int
output resolution
Returns: map_out : array_like
up/down graded map
-
smart_ud_grade
(map_in, cov_in, nside_out, pess=False)¶ Upgrade or downgrade the resolution of a healpix map
When degrading the resolution, the value of the output superpixel is set to the inverse variance weighted mean of the input children pixels. When upgrading the resolution, the value of the children pixels are set to that of the superpixel.
When degrading the resolution, the input inverse variances are summed to give the superpixel. When upgrading the resolution, the input inverse variance is divided by the number of pixels.
The default healpy.ud_grade calculates the average of children pixels. Multiplying the average by the number of daughter pixels in a superpixel returns their sum. ud_grade has the option to multiply the outputs by the number of children, setting power=-2.
Parameters: map_in : array_like
input map
cov_in : array_like
input covariance
nside_out : int
desired resolution of output map
Returns: map_out : array_like
up/down graded map
cov_out : array_like
up/down graded cov