search

mousset / qubic

Data analysis tools for the QUBIC experiment.
3 stars 2 forks source link

Divide by zero in multiacquisition.py #7

Closed planck2007 closed 4 years ago

planck2007 commented 4 years ago

During the computation of a map via tod2map the following warning is issued: multiacquisition.py:63: RuntimeWarning: divide by zero encountered in true_divide for (mi, ma) in self.bands]) The computation continues but the warning might show a problem. I attach the dict file used to make the simulation TD_spectroimaging.dict.txt

mousset commented 4 years ago

Yes this warning is always here, I think it doesn't depend on the dictionary used. Don't know why...

jchamilton75 commented 4 years ago

I've been looking at the origin of the warning... it's rather harmless. It happens in the function that calculates the preconditionner for the conjugate gradient to converge faster (so it has no - or very marginal - effect on the final result, only on the convergence efficiency). What happens is that we calculate the inverse of the coverage in each sub-bands and then average them int each reconstruction sub-band. So the coverage maps contain zeros in the non observed pixels resulting in a division by zero. Then the next line of code fixes the inf created by setting them to zero. The way it is written is rather compact and avoiding the division by zero would make the code less compact... So we have two choices here:

Now I have a worry here... I am not sure it really makes sense to invert the coverage and then average the inverse coverage... it seems to me that we should instead average the coverage before inversion. So this could be the opportunity to improve a little bit this part...

Let's discuss it a bit and then decide what we change in the code.

mousset commented 4 years ago

We chose the 2nd option: implement the code that avoids divisions by zero in order to clean the output.

As Martin realized, about the worry that JC had, we checked the code and the way it is written is actually fine.

jchamilton75 commented 4 years ago

Here it is. I have modified the code so that the warning does not appear anymore. The results is exactly the same. Issue solved.