Commit 6b0228d2 authored by alexander.dornheim's avatar alexander.dornheim
Browse files

Fixed a bug where ProxOperator Approx_P_JWST_Poission changed argument

parent 7af0eca0
......@@ -79,8 +79,10 @@ def JWST_data_processor(config):
for j in range(diversity):
Pj = np.multiply(Xi_A,exp(1j*(aberration[:,:,j]+theta)));
print(Pj[81,81]);
k[:,:,j] = np.square(abs(Utilities.FFT(Pj)));
print(Pj[87,97]);
k[:,:,j] = np.square(abs(Utilities.FFT(Pj)));
print(k[87,97,j]);
print(norm(k[:,:,j]));
#clear Pj temp1 temp2 defocus Automatic garbage collection?
......
......@@ -5,7 +5,7 @@ from proxtoolbox import Algorithms
from proxtoolbox import ProxOperators
from proxtoolbox.ProxOperators.proxoperators import ProxOperator
from numpy.linalg import norm
from numpy import square, sqrt
from numpy import square, sqrt, nonzero
class Phase(Problem):
......@@ -158,6 +158,8 @@ class Phase(Problem):
proj2 = self.config['proxoperators'][1](self.config)
u_2 = proj2.work(u_1);
print(norm(self.config['u_0'][:,:,0]));
# estimate the gap in the relevant metric
if self.config['Nx'] ==1 or self.config['Ny']==1 :
tmp_gap = square(norm(u_1-u_2)/self.config['norm_rt_data']);
......@@ -168,6 +170,11 @@ class Phase(Problem):
tmp_gap = tmp_gap+(norm(u_1[:,:,j]-u_2[:,:,j])/square(self.config['norm_rt_data']));
gap_0=sqrt(tmp_gap);
print(gap_0);
print(norm(u_1[:,:,0]))
print(norm(u_2[:,:,0]))
print(proj1)
print(proj2)
# sets the set fattening to be a percentage of the
# initial gap to the unfattened set with
......@@ -197,6 +204,14 @@ class Phase(Problem):
self.u1,self.u2,self.iters,self.change,self.gap = \
self.algorithm.run(self.config['u_0'],self.config['TOL'],self.config['MAXIT'])
print(self.iters)
#print(nonzero(self.u1))
#print(self.u1[self.u1.nonzero()]);
print(norm(self.u1));
print(norm(self.u2));
#print(self.gap);
#print(self.change);
def _postsolve(self):
......
......@@ -255,11 +255,13 @@ class Approx_P_JWST_Poisson(ProxOperator):
data = self.data;
product_space_dimension = self.product_space_dimension;
abs_illumination = self.abs_illumination
u_new = zeros(u.shape,u.dtype)
for j in range(product_space_dimension-1):
#print(u[:,:,j].nonzero());
#print(u[6,54,j]);
#print(norm(data_sq[:,:,j]))
U = fft2( indicator_ampl *exp(1j*illumination_phase [:,:,j]) * u[:,:,j]);
U_sq = U * conj(U);
tmp = U_sq/data_sq[:,:,j];
......@@ -278,17 +280,20 @@ class Approx_P_JWST_Poisson(ProxOperator):
hU = sum(sum(U_sq *tmp + data_sq[:,:,j] - U_sq));
if hU>=epsilon+TOL2:
U0 = magproj(U,data[:,:,j]); #argument order changed compared to matlab implementation!!!
u[:,:,j] = indicator_ampl *exp(-1j*illumination_phase[:,:,j]) *ifft2(U0);
u_new[:,:,j] = indicator_ampl *exp(-1j*illumination_phase[:,:,j]) *ifft2(U0);
else:
u_new[:,:,j] = u[:,:,j]
#else:
# no change
#print(norm(u[:,:,j]));
# now project onto the pupil constraint.
# this is a qualitative constraint, so no
# noise is taken into account.
j=product_space_dimension-1;
u[:,:,j] = magproj(u[:,:,j],abs_illumination);#argument order changed compared to matlab implementation!!!
u_new[:,:,j] = magproj(u[:,:,j],abs_illumination);#argument order changed compared to matlab implementation!!!
return u;
return u_new;
#original matlab comment
# P_amp.m
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment