Noise in JWST experiment is now working. Added some noise demos

demos/JWST_RAAR.py → demos/JWST_AvP.py

-# old demo
 
 import SetProxPythonPath
 from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
 
-JWST = JWST_Experiment(algorithm = 'RAAR', lambda_0 = 0.85, lambda_max = 0.85, MAXIT = 1000)
+JWST = JWST_Experiment(algorithm='AvP')
 JWST.run()
 JWST.show()

demos/JWST_AvP_noise.py

+
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
+
+JWST = JWST_Experiment(algorithm='AvP', MAXIT=2000, TOL=5e-15, noise=True)
+JWST.run()
+JWST.show()

demos/JWST_CADRl.py

-# old demo
 
 import SetProxPythonPath
 from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
 
-JWST = JWST_Experiment(algorithm = 'CADRl', lambda_0 = 1.0, lambda_max = 1.0, MAXIT = 1000)
+JWST = JWST_Experiment(algorithm='CADRl', formulation='cyclic', 
+                       lambda_0=1.0, lambda_max=1.0, MAXIT=6000, rotate=True)
 JWST.run()
 JWST.show()

demos/JWST_CADRl_noise.py

+
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
+
+JWST = JWST_Experiment(algorithm='CADRl', formulation='cyclic', 
+                       lambda_0=1.0, lambda_max=1.0, MAXIT=6000, 
+                       noise=True, rotate=True)
+JWST.run()
+JWST.show()

demos/JWST_CDRl.py

+
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
+
+JWST = JWST_Experiment(algorithm='CDRl', formulation='cyclic',
+                       lambda_0=0.25, lambda_max=0.25, MAXIT=600, rotate=True)
+JWST.run()
+JWST.show()

demos/JWST_CDRl_noise.py

+
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
+
+JWST = JWST_Experiment(algorithm='CDRl', formulation='cyclic', 
+                       lambda_0=0.25, lambda_max=0.25, 
+                       noise=True, rotate=True)
+JWST.run()
+JWST.show()

demos/JWST_CP.py

-# old demo
 
 import SetProxPythonPath
 from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
 
-JWST = JWST_Experiment(algorithm = 'CP', formulation = 'cyclic')
+JWST = JWST_Experiment(algorithm='CP', formulation='cyclic', MAXIT=1000,
+                       TOL=5e-15, rotate=True)
 JWST.run()
 JWST.show()

demos/JWST_CP_noise.py

+
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
+
+JWST = JWST_Experiment(algorithm='CP', formulation='cyclic', MAXIT=1000,
+                       noise=True, rotate=True)
+JWST.run()
+JWST.show()

demos/JWST_DRAP.py

@@ -2,6 +2,6 @@
 import SetProxPythonPath
 from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
 
-JWST = JWST_Experiment(algorithm = 'DRAP', lambda_0 = 0.65, lambda_max = 0.65)
+JWST = JWST_Experiment(algorithm='DRAP', lambda_0=0.65, lambda_max=0.65)
 JWST.run()
 JWST.show()

demos/JWST_DRAP_noise.py

+
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
+
+JWST = JWST_Experiment(algorithm='DRAP', lambda_0=0.25, lambda_max=0.25,
+                       noise=True)
+JWST.run()
+JWST.show()

demos/JWST_DRl.py

+
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
+
+JWST = JWST_Experiment(algorithm='DRl', lambda_0=0.55, lambda_max=0.55)
+JWST.run()
+JWST.show()

demos/JWST_DRl_noise.py

+
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
+
+JWST = JWST_Experiment(algorithm='DRl', lambda_0=0.55, lambda_max=0.55, noise=True)
+JWST.run()
+JWST.show()

demos/JWST_QNAvP_noise.py

+
+import SetProxPythonPath
+from proxtoolbox.experiments.phase.JWST_Experiment import JWST_Experiment
+
+JWST = JWST_Experiment(algorithm = 'QNAvP', noise=True)
+JWST.run()
+JWST.show()

proxtoolbox/experiments/experiment.py

@@ -85,7 +85,7 @@ class Experiment(metaclass=ExperimentMetaClass):
                  Nx=None,
                  Ny=None,
                  Nz=None,
-                 noise=None,
+                 noise=False,
                  algorithm=None,
                  numruns=1,
                  MAXIT=1000,

proxtoolbox/experiments/phase/JWST_Experiment.py

@@ -80,7 +80,14 @@ class JWST_Experiment(PhaseExperiment):
         # So far Matlab code used only Nx to set the desired
         # resolution
         newres = self.Nx
-        snr = self.data_ball
+        
+        # The following value for snr does not work well in Python
+        # as this creates numerical issues with Poisson noise 
+        # (data_ball is far too small)
+        # snr = self.data_ball
+        # Instead, we choose the following bigger value which  
+        # works well.
+        snr = 1e-8 
 
         print('Loading data.')
         with open('../InputData/Phase/pupil.pmod','r') as fid:
@@ -137,7 +144,7 @@ class JWST_Experiment(PhaseExperiment):
             Pj = Xi_A*exp(1j*(aberration_j+theta))
             k_j = np.square(abs(utils.FFT(Pj)))
             rt_k_j = np.sqrt(k_j)
-            if self.noise == 'Poisson':
+            if self.noise:
                 # Add Poisson noise according to Ning Lei:
                 # f = alpha*4*4*pi/q/q*abs(cos(qx(iX))*sin(qy(iY))*(sin(q)/q-cos(q)));
                 #        f2 = PoissonRan(f*f*2)/alpha/alpha/2;   # 2 is for I(-q)=I(q)
@@ -147,7 +154,8 @@ class JWST_Experiment(PhaseExperiment):
                 k_j = k_j/snr
                 for ii in range(newres):
                     for jj in range(newres):
-                        k_j[ii,jj]= np.random.poisson(k_j[ii,jj])*snr #use built in numpy possion instead of utils.PoissonRan(k[ii,jj,i])*snr
+                        k_j[ii,jj]= utils.PoissonRan(k_j[ii,jj])*snr # this matches what is currently done in Matlab
+                        #k_j[ii,jj]= np.random.poisson(k_j[ii,jj])*snr # does not work: throw an exception if snr is too small 
                 k_j = np.round(k_j)
                 rt_k_j = np.sqrt(k_j)
 
@@ -156,7 +164,7 @@ class JWST_Experiment(PhaseExperiment):
 
             k[j] = k_j
             rt_k[j] = rt_k_j
-            zeros = np.where(rt_k_j == 0)[0]
+            zeros = np.where(rt_k_j == 0)
             data_zeros[j] = zeros
             self.FT_conv_kernel[j] = self.indicator_ampl*exp(1j*aberration_j)
 
@@ -174,13 +182,13 @@ class JWST_Experiment(PhaseExperiment):
         # self.illumination_phase = aberration -> now input.FT_conv_kernel
 
         # initial guess
-        #self.u0 = zeros((newres,newres,self.product_space_dimension),dtype = np.complex128)
+        perturb = (np.random.rand(newres, newres)).T # to match Matlab 
         if self.formulation == "cyclic":
-            self.u0 = np.multiply(self.abs_illumination, exp(1j*2*pi*np.ones(newres))) #rand(newres)))
+            self.u0 = self.abs_illumination*exp(1j*2*pi*perturb)
         else: # product space
             self.u0 = Cell(self.sets)
             for j in range(self.sets):
-                u0_j = self.abs_illumination*exp(1j*2*pi*np.ones(newres)) #rand(newres)))
+                u0_j = self.abs_illumination*exp(1j*2*pi*perturb) 
                 self.u0[j] = u0_j
 
         self.truth = true_object

proxtoolbox/proxoperators/fourier.py

@@ -60,18 +60,25 @@ class Approx_Pphase_FreFra_Poisson(ProxOperator):
 
         # Propagate to the image plane:
         u_hat = self.propagator.eval(u, prox_idx)
+        
         # Compute the Kullback-Leibler distance of u_hat to the data: 
         if prox_idx is None:
             j = 0
         else:
             j = min(prox_idx, len(self.data)-1)
         u_hat_sq = real(u_hat * conj(u_hat))
-        # update data_sq to prevent division by zero
+        # Update data_sq to prevent division by zero.
+        # Note that data_sq_zeros, as defined below,
+        # can be different from self.data_zeros[j]] 
+        # (e.g., JWST experiment).
         data_sq = self.data_sq[j].copy()
-        data_sq[self.data_zeros[j]] = 1
+        data_sq_zeros = np.where(data_sq == 0)
+        data_sq[data_sq_zeros] = 1
+
         tmp = u_hat_sq / data_sq
-        tmp[self.data_zeros[j]] = 1
-        u_hat_sq[self.data_zeros[j]] = 0
+        data_zeros = self.data_zeros[j]
+        tmp[data_zeros] = 1
+        u_hat_sq[data_zeros] = 0
         I_u_hat = tmp == 0
         tmp[I_u_hat] = 1
         tmp = log(tmp)

proxtoolbox/utils/PoissonRan.py

@@ -9,47 +9,74 @@ from scipy.special import gammaln
 import numpy as np
 
 def PoissonRan(xm):
-  oldm = -1;
-  if xm<12:
-    if xm != oldm:
-      oldm = xm;
-      g = exp(-xm);
-
-    em = -1;
-    t = 1;
-
-    em = em+1;
-    t = t*rand(1);
-    while t>g :
-      em = em+1;
-      t = t*rand(1);
-
-  else:
-    if xm != oldm:
-      oldm = xm;
-      sq = sqrt(2.0*xm);
-      alxm = log(xm);
-      g = xm*alxm-gammaln(xm+1);
-
-    y = tan(pi*rand(1));
-    em = sq*y+xm;
-
-    while em < 0:
-      y = tan(pi*rand(1));
-      em = sq*y+xm;
+    """
+    Poisson random number generation with mean of xm
+
+    See numerical recipe C++ version
+    ref: y(i+1) = xm^i*exp(-xm)/factorial(i); see Leo Breiman Probability
+        
+    Parameters
+    ----------
+    xm : real number
+        Mean
+        
+    Returns
+    -------
+    em : real number
+        Sample from the parameterized Poisson distribution
+
+    Notes    
+    -----
+    This code suffers from numerical precision issues. With big 
+    values for xm (for example, xm = 6.910459246361838e+27, as was 
+    initially the case with the JWST experiment) the difference
+    with Matlab is substential. In Python, one exponential 
+    evalutation gives +inf while it gives zero in Maltab.
+    """
+
+    oldm = -1
+
+    if xm<12:
+        if xm != oldm:
+            oldm = xm
+            g = exp(-xm)
+
+        em = -1
+        t = 1
+
+        em = em+1
+        t = t*rand()
+        while t>g:
+            em = em+1
+            t = t*rand()
+    else:
+        if xm != oldm:
+            oldm = xm
+            sq = sqrt(2.0*xm)
+            alxm = log(xm)
+            g = xm*alxm-gammaln(xm+1)
+
+        y = tan(pi*rand())
+        em = sq*y+xm
+
+        while em < 0:
+            y = tan(pi*rand())
+            em = sq*y+xm
   
-    em = floor(em);
-    t = 0.9*(1+y*y)*exp(em*alxm-gammaln(em+1)-g);
+        em = floor(em)
 
-    while rand(1) > t:
-      y = tan(pi*rand(1));
-      em = sq*y+xm;
+        t = 0.9*(1+y*y)*exp(em*alxm-gammaln(em+1)-g)
 
-      while em < 0:
-        y = tan(pi*rand(1));
-        em = sq*y+xm;
+        while rand() > t:
+            y = tan(pi*rand())
+            em = sq*y+xm
 
-      em = floor(em);
-      t = 0.9*(1+y*y)*exp(em*alxm-gammaln(em+1)-g);
+            while em < 0:
+                y = tan(pi*rand())
+                em = sq*y+xm
+
+            em = floor(em)
+            t = 0.9*(1+y*y)*exp(em*alxm-gammaln(em+1)-g)
    
-  return em;
+    return em
+