Fixed a bug in Approx_P_Poission_JWST, changed structe of data_zeros to list of tuples

proxtoolbox/Problems/Phase/JWST_data_processor.py

@@ -85,7 +85,7 @@ def JWST_data_processor(config):
 
     epsilon=0;
     norm_rt_data=np.sqrt(sum(sum(Xi_A)));
-    data_zeros= np.zeros((newres*newres,3));
+    data_zeros= [[],[],[]];
     for i in range(diversity):
       rt_k[:,:,i]= np.sqrt(k[:,:,i]); # *newres ?
       # normalize the data so that the corresponding 
@@ -111,8 +111,8 @@ def JWST_data_processor(config):
       # fftshift and scale the data for use with fft:
       rt_k[:,:,i] = fftshift(rt_k[:,:,i])*newres;
       k[:,:,i] = fftshift(k[:,:,i])*newres;
-      temp = np.where(rt_k[:,:,i]==0)[1];
-      data_zeros[range(temp.size),i] = temp;
+      data_zeros[i] = rt_k[:,:,i].nonzero() #np.where(rt_k[:,:,i]==0)[1];
+      #data_zeros[range(temp.size),i] = temp;
 
     config['rt_data'] = rt_k;
     config['data'] = k;
@@ -142,7 +142,7 @@ def JWST_data_processor(config):
     config['truth_dim'] = true_object.shape;
     config['norm_truth'] = norm(true_object);
     
-    print(config['norm_rt_data']);
+    #print(config);
 
 
 

proxtoolbox/Problems/Phase/JWST_in.py

@@ -60,7 +60,7 @@ new_config = {
 
 
 # What are the noise characteristics (Poisson or Gaussian or none)?
-    'noise' : 'Poisson', #'Poisson',
+   'noise' : 'Poisson', #'Poisson',
 
 #==========================================
 # Algorithm parameters

proxtoolbox/Problems/Phase/phase.py

@@ -135,7 +135,47 @@ class Phase(Problem):
         # and adjust your input files and projectors accordingly. 
         # you could also do this within the data processor
 
-        self.config['TOL2'] = 1e-15; 
+        self.config['TOL2'] = 1e-15;
+
+        #To estimate the gap in the sequential formulation, we build the
+        # appropriate point in the product space. This allows for code reuse.
+        # Note for sequential diversity diffraction, input.Proj1 is the "RCAAR"
+        # version of the function.
+        if formulation == 'sequential':
+            for j in range(self.config['product_space_dimension']):
+                self.config['proj_iter'] =j;
+                proj1 = self.config['proxoperators'][0](self.config)
+                u_1[:,:,j]= proj1.work(self.config['u_0']);
+                self.config['proj_iter'] = mod(j,config['product_space_dimension'])+1;
+                proj1 = self.config['proxoperators'][0](self.config)
+                u_1[:,:,j]= proj1.work(self.config['u_0']);
+            end;
+        else: #i.e. formulation=='product space'
+            proj1 = self.config['proxoperators'][0](self.config)
+            u_1 = proj1.work(self.config['u_0']);
+            proj2 = self.config['proxoperators'][1](self.config)
+            u_2 = proj2.work(self.config['u_0']);
+
+        # estimate the gap in the relevant metric
+        if self.config['Nx'] ==1 or self.config['Ny']==1 :
+            tmp_gap = (norm(u_1-u_2)/input.norm_rt_data)^2;
+        else:
+            tmp_gap=0;
+            for j in range(self.config['product_space_dimension']):
+                # compute (||P_Sx-P_Mx||/normM)^2:
+                tmp_gap = tmp_gap+(norm(u_1[:,:,j]-u_2[:,:,j])/config['norm_rt_data'])^2;
+
+        gap_0=sqrt(tmp_gap);
+
+        # sets the set fattening to be a percentage of the
+        # initial gap to the unfattened set with 
+        # respect to the relevant metric (KL or L2), 
+        # that percentage given by
+        # input.data_ball input by the user.
+        input.data_ball=input.data_ball*gap_0;
+        # the second tolerance relative to the oder of 
+        # magnitude of the metric
+        input.TOL2 = input.data_ball*1e-15; 
 
         print(self.config);
 
@@ -147,7 +187,7 @@ class Phase(Problem):
         """
         Prepares argument for actual solving routine
         """
-    
+        
     
     def _solve(self):
         """

proxtoolbox/ProxOperators/proxoperators.py

@@ -7,7 +7,9 @@ Created on Mon Sep 22 13:23:26 2014
 The "ProxOperators"-module contains various specific operators that do the actual calculations within the ProxToolbox.
 """
 
-from numpy import conj, dot, empty, ones, sqrt, sum, zeros
+import numpy as np
+from numpy import conj, dot, empty, ones, sqrt, sum, zeros, exp, nonzero, log
+from numpy.fft import fft2
 
 __all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson", "P_amp"]
 
@@ -237,6 +239,7 @@ class Approx_P_JWST_Poisson(ProxOperator):
         self.illumination_phase = config['illumination_phase'];
         self.data_zeros = config['data_zeros'];
         self.data_sq = config['data_sq'];
+        self.product_space_dimension = config['product_space_dimension'];
 
     def work(self,u):
         TOL2 = self.TOL2;
@@ -245,19 +248,23 @@ class Approx_P_JWST_Poisson(ProxOperator):
         illumination_phase = self.illumination_phase;
         data_zeros = self.data_zeros;
         data_sq = self.data_sq;
-        
+        product_space_dimension = self.product_space_dimension;
 
         for j in range(product_space_dimension-1):
             U = fft2( indicator_ampl *exp(1j*illumination_phase [:,:,j]) * u[:,:,j]);
             U_sq = U * conj(U);
             tmp = U_sq/data_sq[:,:,j];
-            tmp2= data_zeros[:,j] != 0;
-            tmpI = data_zeros[tmp2,j];
-            tmp[tmpI]=1;
-            U_sq[tmpI]=0;
+            #tmp2= nonzero(data_zeros[:,j]);
+            #tmpI = data_zeros[data_zeros[:,j] != 0, j];
+            mask = ones(tmp.shape, np.bool);
+            print(mask.shape);
+            mask[data_zeros[j]] = 0;
+            tmp[mask]=1;
+            U_sq[mask]=0;
             IU= tmp==0;
             tmp[IU]=1;
             tmp=log(tmp);
+            print(tmp);
             hU = sum(sum(U_sq *tmp + data_sq[:,:,j] - U_sq));
             if hU>=epsilon+TOL2:
                 U0 = magproj(U,input.data[:,:,j]); #argument order changed compared to matlab implementation!!!
@@ -268,7 +275,7 @@ class Approx_P_JWST_Poisson(ProxOperator):
         # now project onto the pupil constraint.
         # this is a qualitative constraint, so no 
         # noise is taken into account.
-        j=input.product_space_dimension;
+        j=product_space_dimension;
         u[:,:,j] = magproj(u[:,:,j],abs_illumination);#argument order changed compared to matlab implementation!!!
         
         return u;