GitLab

Cleaned up some output statements

Fixed 2 bugs in data_processor (S, matrix multiplication)
Fixed one bug in Approx_PM_Gaussion (use iift2 instead of ifft)

Tried naive implementation of magproj, should work now, but commented out
Replaced u*conj(u) by u.real**2+u.imag**2 in magproj. Computation time for 500 JWST iterations is down from 1.3s to 0.5 s (much faster!!!). Computation before was complex!!!
Did the same in Approx_P_JWST_Poisson (also somewhat faster, since using real muliplication, log now).
Implemented P_S

Phase was adjusted for new output
This should allow different algorithms to return different arguments
Still need to adjust other algorithms, Sudoku, pythography...

Added test for tasse with AP, RAAR with beta_0=beta_max=0.5

Added another test

Added a test function to check CDI_data_processors to phase

Updated RAAR algortihm (should now work identically to matlab version), small changes in other files

Moved definition of use_farfield_fromula and fresnel_nr from JWST_in to phase since if statements are not possible in dictonary and user only needs to specify distance

JWST computation is now double precision, changed type of variables in JWST_data_processor to float64, complex128

Compared result to matlab result:

After 1 iteration:

Maximum norm of difference:
2.25868816441e-12
Frobenius norm of difference:
1.9427347826e-11

After 500 iterations:

Maximum norm of difference:
8.35619946335e-05
Frobenius norm of difference:
0.00514996012436

So JWST seems to be working now...
Also epsilon in magproj was changed back to 1e-30 (no overflow with double precision)

Do some precalculations in Approx_P_JWST_Poisson for better performance. This modification results in ~20% less computation time for JWST (500 Iterations).

In phase moved comparision with matlab results to method compare_to_matlab. This is only intended to be used for testing. Note u_0 needs to be fixed in a specific way, noise should be switched off.
Added test for 1 Iteration: After 1 Iteration maximum norm of difference is 5.4e-06, Forbenius norm of difference is 2.5e-05.

ProxPython is now using the built-in random generator numpy.random.poisson. According to the comments in PossionRan the argument is the mean (often called lambda). I also checked Numerical Recipes p. 294. Therefore calling

numpy.random.poisson(x)

should yield the same as calling

PoissonRan(x).

This resolves the overflow issue with PoissonRan, but I still get the following warnings:

../proxtoolbox/ProxOperators/proxoperators.py:265: RuntimeWarning: divide by zero encountered in true_divide
tmp = U_sq/data_sq[:,:,j];
../proxtoolbox/ProxOperators/proxoperators.py:265: RuntimeWarning: invalid value encountered in true_divide
tmp = U_sq/data_sq[:,:,j];

Therefore I had a look at the following line in JWST_data_processor:

k[:,:,i]=np.round(k[:,:,i])

From my understanding it is very likely to get zeros when rounding.

I noted that after commenting out above line, there were no further warnings and the algorithm did converge nicely.

Adjusted Input, import paths in JWST_data_processor