Refactoring the C++ extension codes, introducing pointing flags with code optimization

[anand-avinash]

This PR introduces the following breaking changes in the code base:

• Refactoring the C++ extension codes
  • ProcessTimeSamples
  • SparseLO --> PointingLO
  • BlockDiagonalPreconditionerLO
• Refactoring the Python interface
  • brahmap/utilities/process_time_samples.py
  • brahmap/interfaces/linearoperators.py
• Introducing the pointing flags
• Optimization of the for-loops with nested conditionals

[anand-avinash]

In the last commit, I:

• introduced the pointing flags.

• refactored brahmap/utilities/process_time_samples.py, its C++ extensions and relevant tests

• Attempted to optimize the for-loops structures like this:

  for (int idx=0; idx<nsamples; ++idx) {
    int pixel = pointings[idx];
    if (flag[idx]) { // boolean flag
      quantity[pixel] += set1[idx]*set2[idx];
    } // if
  } // for

  with this:

  for (int idx=0; idx<nsamples; ++idx) {
    int pixel = pointings[idx];
    quantity[pixel] += (flag[idx)*set1[idx]*set2[idx];
  } // for

  Benchmarking is pending for this attempt.

[anand-avinash]

In the _rmult_* methods of PointingLO, we had loops like this:

for i in range(Nrows):
    if pointings[i] == -1:
        continue
    vec_out[2 * pointings[i]] += v[i] * cos[i]
    vec_out[2 * pointings[i] + 1] += v[i] * sin[i]

and in _mult_*, we had loops like this:

for i in range(Nrows):
    if pointings[i] == -1:
        continue
    x[i] += v[2 * pointings[i]] * cos[i]
    +v[2 * pointings[i] + 1] * sin[i]

We can translate _mult_* above with pointings_flag as following:

for i in range(Nrows):
    pixel = pointings[i]
    x[i] += pointings_flag[i] * (v[2*pixel] * cos[i]
            + v[2*pixel+1] * sin[i])

However, since PointingLO is now based on new_npix, v[2*pixel] can simply go out of range of v. In the original version of the code, the pointings for the flagged pixels were changed to -1, thus avoiding the out of range error. In our case, this issue can be solved by changing the in the pointings array the pathological pixels to 0 while flagging the sample sample. But changing the pointings array and pointings_flag seems little expensive. So we introduced in ProcessTimeSamples the pixel_flag given by

pixel_flag = np.zeros(npix, dtype=bool)
for pixel in pixel_mask:
    pixel_flag[pixel] = True

With this, we can now write the _mult_* as

for idx in range(Nrows):
    pixel = pointings[idx]
    pixflag = pixel_flag[pixel]
    pixel *= pixflag   # This marks pathological pixel to pixel index 0
    pointflag = pixflag && pointings_flag[idx]  # This flags the sample corresponding to pathological pixel to zero

    x[i] += pointings_flag[i] * (v[2*pixel] * cos[i]
            + v[2*pixel+1] * sin[i])

The same has been implemented with the previous commit.

[anand-avinash]

Picking from the last comment, we had the for loop for _mult_* given by

for idx in range(Nrows):
    pixel = pointings[idx]
    pixflag = pixel_flag[pixel]
    pixel *= pixflag   # This marks pathological pixel to pixel index 0
    pointflag = pixflag && pointings_flag[idx]  # This flags the sample corresponding to pathological pixel to zero

    x[i] += pointings_flag[i] * (v[2*pixel] * cos[i]
            + v[2*pixel+1] * sin[i])

While it successfully deals with pathological pixels stored in the pointings array, v[2*pixel] can still go out of range. Note that the pixels in the pointings array range from 0 to npix while the shape of PointingLO is now nrows = nsamp, ncols = new_npix x solver_type. So the pixel in this loop must range from 0 to new_npix x solver_type. We could not come up with a solution to this problem without relying on if-statements. In the end, we decided to update the pointings and pointings_flag when ProcessTimeSamples is called such that,

1. pointings[pathological_pixels] = 0: This is just to avoid out-of-range errors/warnings when the pointing sample corresponding to a pathological pixel is accessed.
2. pointings[idx] = old2new_pixel[pointings[idx]]: Once we flagged the pathological pixels, we can renew the pixel index such that they vary from 0 to new_npix (where new_npix = npix - n_pathological_pixels). The array old2new_pixel simply maps the old pixel indices to new pixel indices: new_pix_idx = old2new_pixel[old_pix_idx]. In fact, this was already in place in the earlier versions of our code.
3. pointings_flag[pathological_pixels] = False: To ensure that pointing samples corresponding to pathological pixels are excluded from PointingLO operator and eventually the entire map-making pipeline. This makes sure that the pointing samples set to zero (see the first point) will not be accounted in the map-making.

With these changes, the loop becomes

for idx in range(Nrows):
    pixel = pointings[idx]
    pointflag = pointings_flag[idx]

    x[i] += pointings_flag[i] * (v[2*pixel] * cos[i]
            + v[2*pixel+1] * sin[i])

[giuspugl]

Shouldn’t the flag be applied to both operations ,e.g.  x[i] += pointings_flag[i] * (v[2*pixel] * cos[i] + v[2*pixel+1] * sin[i])

[anand-avinash]

@giuspugl Yes, it has already been applied to both operators.

[anand-avinash]

With respect to the refactoring done so far, the previous class/function interfaces are now incompatible with the new one. In light of this, I have now implemented the function brahmap.mapmakers.compute_GLS_maps() to provide a general interface for GLS map-maker. Another function is brahmap.mapmakers.LBSim_compute_GLS_maps() that provides the interface for GLS map-maker from litebird_sim.

[anand-avinash]

For noise-less time-ordered data (TOD) from litebird_sim in the absence of half-wave plate (hwp), we expect the preconditioned conjugate gradient (PCG) solver to converge in a single iteration, as the preconditioner linear operator would be a diagonal operator in this case. However, we found that the PCG was never converging. Digging in the linear operators, we found that few 3x3 diagonal blocks in the preconditioner operator were zero, that was not what we expect. Turned out it was due to the test we were performing during multiplication between BlockDiagonalPreconditionerLO and a vector. The test was intended to identify and ignore the singular $3\times3$ matrix blocks.

To fix this issue, we revisited the function get_pixel_mask_pol(). In this function, we flagged all the pixels where the condition number of the following matrix was above certain threshold:

$$\begin{bmatrix} \sum\cos^2 2\phi_i/w_i & \sum\sin2\phi_i\cos2\phi_i/w_i\\ \sum\sin2\phi_i\cos2\phi_i/w_i & \sum\sin^2 2\phi_i/w_i \end{bmatrix}$$

We decided to replace this flagging scheme with the one where we flag all the pixels that might result in singular matrix block in the BlockDiagonalPreconditionerLO. For QU solver, this block matrix will be

$$\begin{bmatrix}\sum\cos^2 2\phi_i/w_i & \sum\sin2\phi_i\cos2\phi_i/w_i\\ \sum\sin2\phi_i\cos2\phi_i/w_i & \sum\sin^2 2\phi_i/w_i \end{bmatrix}$$

whereas for IQU solver, this matrix block will be

$$\begin{bmatrix} \sum 1/w_i & \sum\cos2\phi_i/w_i & \sum\sin2\phi_i/w_i \\ \sum\cos2\phi_i/w_i & \sum\cos^2 2\phi_i/w_i & \sum\sin2\phi_i\cos2\phi_i/w_i\\ \sum\sin2\phi_i/w_i & \sum\sin2\phi_i\cos2\phi_i/w_i & \sum\sin^2 2\phi_i/w_i \end{bmatrix}$$

Eventually in the new pixel flagging scheme, we are computing the determinant of the matrix blocks in get_pixel_mask_pol() and then flagging it if the determinant is below certain threshold. Since the determinant is now a member of ProcessTimeSamples class, now we won't compute it every time the matrix vector multiplication of BlockDiagonalPreconditionerLO is called.

[giuspugl]

Everything looks good, few observations :

1. i think even in presence of HWP signal only sims should converge with 1 iteration
2. please remember the factor of 2 in the trigonometric functions , e.g. $\sin 2 \phi $

[anand-avinash]

i think even in presence of HWP signal only sims should converge with 1 iteration

Yes, it does converge with one iteration in the presence of HWP

please remember the factor of 2 in the trigonometric functions , e.g. $\sin2\phi$

Thanks, I have now updated the comment!