# 2 Usage¶

## 2.1 Basic usage¶

The most common use case for `symlens`

is to get noise curves for lensing, and
to run lensing reconstruction on CMB maps. We will cover this first, and later
come to custom estimators that exploit the full power and flexibility of this
package.

### 2.1.1 Geometry¶

Skip this if you are familiar with the pixell library. Most functions that
return something immediately useful take a `shape`

,``wcs`` pair
as input. In short, this pair specifies the footprint or geometry of the patch of
sky on which calculations are done. If you are playing with simulations or just
want quick forecasts, you can make your own `shape`

,``wcs`` pair as follows:

```
>>> from pixell import enmap
>>> shape,wcs =
enmap.geometry(shape=(2048,2048),res=np.deg2rad(2.0/60.),pos=(0,0))
```

to specify a footprint consisting of 2048 x 2048 pixels of width 2 arcminutes and centered on the origin of a CEA geometry.

All outputs are 2D arrays in Fourier space, so you will need some way to bin it in annuli. A map of the absolute wavenumbers is useful for this

```
>>> modlmap = enmap.modlmap(shape,wcs)
```

To read a map in from disk and get its geometry, you could do

```
>>> imap = enmap.read_map(fits_file_name)
>>> shape,wcs = imap.shape, imap.wcs
```

Please read the documentation for pixell for more information.

### 2.1.2 Noise curves¶

Using the machinery described in `Custom estimators`

below, a number of
pre-defined mode-coupling noise curves have been built. We provide some examples
of using these. All of these require a `shape`

,``wcs`` geometry pair as described
above. Next, you also need to have a Fourier space mask at hand
that enforces what multipoles in the CMB map are included.

A convenience function is provided for generating simple Fourier space masks. e.g.

```
>>> from symlens import utils
>>> tellmin = 100
>>> tellmax = 3000
>>> kmask = utils.mask_kspace(shape,wcs,lmin=tellmin,lmax=tellmax)
```

Finally, you also need a `feed_dict`

, a dictionary which maps names of variables (keys) to
2D arrays containing data, filters, etc. which are fed in at the very final
integration step. With custom estimators described later, you get to choose the
names of your variables. But the convenience of the canned functions described
here comes with the cost of having to learn what variable convention is defined
inside them. We will learn by example.

Lensing noise curves require CMB power spectra. The naming convention for the
`feed_dict`

for these is:

`uC_X_Y`

for CMB XY spectra that go in the lensing response function, e.g.`uC_T_T`

for the TT spectrum (despite the notation, this should be the lensed spectrum, not the unlensed spectrum)`tC_X_Y`

for total CMB XY spectra that go in the lensing filters. e.g.`tC_T_T`

for the total TT spectrum that includes beam-deconvolved noise.

These have to be specified on the 2D Fourier space grid. We can build them like this:

```
>>> feed_dict = {}
>>> feed_dict['uC_T_T'] = utils.interp(ells,ctt)(modlmap)
>>> feed_dict['tC_T_T'] = utils.interp(ells,ctt)(modlmap)+(33.*np.pi/180./60.)**2./utils.gauss_beam(modlmap,7.0)**2.
```

where I’ve used the convenience function `interp`

to interpolate an `ells`

,``cltt``
1D spectrum specification isotropically on to the Fourier space grid, and
created a Planck-like total beam-deconvolved spectrum using the `gauss_beam`

function. That’s it! Now we can get the pre-built Hu Okamoto 2001
(estimator=”hu_ok”) noise for the TT lensing estimator as follows,

```
>>> import symlens as s
>>> nl2d = s.N_l(shape,wcs,feed_dict,"hu_ok","TT",xmask=kmask,ymask=kmask)
```

which can be binned in annuli to obtain a lensing noise curve.

### 2.1.3 Lensing maps¶

To make a lensing map, we need to provide beam deconvolved Fourier maps of the CMB, which for a quadratic estimator <XY> have default variable names of X and Y,

```
>>> feed_dict['X'] = beam_deconvolved_fourier_T_map
>>> feed_dict['Y'] = beam_deconvolved_fourier_T_map
```

An important thing to remember is that by default, the code expects “physical” normalization of FFTs in pixell (not the default normalization in pixell), i.e. you should be passing in Fourier maps that come from something like

```
>>> beam_deconvolved_fourier_T_map = enmap.fft(imap,normalize='phys')
```

or

```
>>> beam_deconvolved_fourier_T_map = enmap.map2harm(imaps,normalize='phys')[0]
```

One can then obtain the unnormalized lensing map simply by doing,

```
>>> ukappa = s.unnormalized_quadratic_estimator(shape,wcs,
feed_dict,"hu_ok","TT",xmask=kmask,ymask=kmask)
```

and also obtain its normalization,

```
>>> norm = s.A_l(shape,wcs,feed_dict,"hu_ok","TT",xmask=kmask,ymask=kmask)
```

and combine into a normalized Fourier space CMB lensing convergence map,

```
>>> fkappa = norm * ukappa
```

### 2.1.4 General noise curves¶

To perform more complicated calculations like cross-covariances, noise for
non-optimal estimators, mixed experiment estimators (for gradient cleaning),
split-based lensing N0 curves, etc., we need to learn how to attach field names,
which make the `feed_dict`

expect more variables than what was described
earlier.

Let’s first show how we can obtain a general noise cross-covariance. We can for example obtain the same TT lensing noise curve as above but in a more round-about way by asking what the cross-covariance of the TT estimator is with the TT estimator itself,

```
>>> Nl =
N_l_cross(shape,wcs,feed_dict,alpha_estimator="hu_ok",alpha_XY="TT",
beta_estimator="hu_ok",beta_XY="TT",
xmask=kmask,ymask=kmask)
```

This works just like before. However, what if the instrument noise in the first leg of the estimator is uncorrelated with the noise in the second leg? Then, we need to differentiate between the four fields that appear above. We can do that by providing names for these fields.

```
>>> Nl = N_l_cross(shape,wcs,feed_dict,
alpha_estimator="hu_ok",alpha_XY="TT",
beta_estimator="hu_ok",beta_XY="TT",
xmask=kmask,ymask=kmask,
field_names_alpha=['E1','E2'],
field_names_beta=['E1','E2'])
```

This modifies the total power spectra variable names that feed_dict expects. The
above command will not work unless `tC_E1_T_E1_T`

, `tC_E2_T_E2_T`

,
`tC_E1_T_E2_T`

, `tC_E2_T_E1_T`

are also provided, instead of just the usual
`tC_T_T`

. Specifying these in feed_dict allows one to generalize to a wider
variety of estimators.

### 2.1.5 Other built-in estimators¶

The following are currently available:

- Hu Okamoto 2001 TT, TE, EE, EB, TB
- Hu DeDeo Vale 2007 TT, TE, ET, EE, EB, TB
- Schaan, Ferraro 2018 shear TT

For the shear estimator, the built-in variable scheme also expects duC_T_T , the logarithmic derivative of the lensed CMB temperature,

```
>>> feed_dict['duC_T_T'] =
utils.interp(ells,np.gradient(np.log(lcltt),np.log(ells)))(modlmap)
```

Once this is added to feed_dict, noise curves and shear maps can be obtained as before,

```
>>> Nl = s.N_l(shape,wcs,feed_dict,"shear","TT",
xmask=tmask,ymask=tmask)
>>> Al = s.A_l(shape,wcs,feed_dict,"shear","TT",xmask=tmask,ymask=tmask)
>>> ushear = s.unnormalized_quadratic_estimator(shape,wcs,feed_dict,"shear","TT",xmask=tmask,ymask=tmask)
>>> shear = Al * ushear
```

## 2.2 Custom estimators¶

We can build general factorizable quadratic estimators as follows.

We need to specify the mode coupling form (little f):

and specify the filter form (big F):

For reference, these are related to the quadratic estimator,

and normalization,

where .

The expressions and must be specified in terms of the following special symbols:

- Ldl1 for
- Ldl2 for
- cos2t12 for
- sin2t12 for
- L for

and any other arbitrary symbols which will be replaced with numerical data later on.

The special symbols can be accessed directly from the module, e.g.:

```
>>> import symlens as s
>>> s.Ldl1
>>> s.L
>>> s.cost2t12
```

and arbitrary symbols can be defined either as functions of l1 or of l2, using a wrapper in the module:

```
>>> s.e('X_l1')
>>> s.e('Y_l2')
```

The ‘_l1’ or ‘_l2’ suffix for arbitrary symbols is critical for the factorizer to know. With these, a large variety of estimators and noise functions can be built, including lensing, magnification, shear, birefringence, patchy tau, mixed estimators (for gradient cleaning), split lensing estimators, etc.

e.g., we can build an integrand for the Hu, Okamoto 2001 TT lensing estimator normalization as follows,

```
# Build HuOk TT estimator integrand
>>> f = s.Ldl1 * s.e('uC_T_T_l1') + s.Ldl2 * s.e('uC_T_T_l2')
>>> F = f / 2 / s.e('tC_T_T_l1') / s.e('tC_T_T_l2')
>>> expr1 = f * F # this is the integrand
```

We then provide data arrays for use after factorization in `feed_dict`

. These are lensed TT spectra interpolated on to 2D Fourier space.

```
>>> feed_dict = {}
>>> feed_dict['uC_T_T'] = utils.interp(ells,ctt)(modlmap)
>>> feed_dict['tC_T_T'] = utils.interp(ells,ctt)(modlmap)
```

For the integral to be sensible, we must also mask regions in Fourier space we don’t want to include.

```
>>> tellmin = 10 ; tellmax = 3000
>>> xmask = utils.mask_kspace(shape,wcs,lmin=tellmin,lmax=tellmax)
```

With these in hand, we can call the core function in symlens for the factorized integral.

```
>>> integral = s.integrate(shape,wcs,feed_dict,expr1,xmask=xmask,ymask=xmask).real
```