Commit bd635273 authored by alexander.dornheim's avatar alexander.dornheim
Browse files

JWST computation is now double precision, changed type of variables in...

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)
parent 6fc25f2a
......@@ -7,5 +7,5 @@ from phase import Phase
JWST = Phase(JWST_in.new_config)
JWST.solve()
#JWST.compare_to_matlab()
JWST.compare_to_matlab()
#JWST.show()
......@@ -39,6 +39,7 @@ def JWST_data_processor(config):
with open('../InputData/phase_retrieval/pupil.pmod','r') as fid:
# read lower endian float <f
Xi_A = np.fromfile(fid, dtype='<f');
Xi_A = Xi_A.astype(np.float64)
Xi_A = Xi_A.reshape((512,512)).T;
diversity=3;
......@@ -46,11 +47,13 @@ def JWST_data_processor(config):
with open('../InputData/phase_retrieval/phase_p37.pmod','r') as fid:
# read lower endian float <f
temp1 = np.fromfile(fid, dtype='<f');
temp1 = temp1.astype(np.float64)
temp1 = temp1.reshape(512,512).T;
with open('../InputData/phase_retrieval/phase_m37.pmod','r') as fid:
# read lower endian float <f
temp2 = np.fromfile(fid, dtype='<f');
temp2 = temp2.astype(np.float64)
temp2 = temp2.reshape(512,512).T;
defocus=(temp1-temp2)/2;
......@@ -73,11 +76,10 @@ def JWST_data_processor(config):
rt_k= np.zeros((newres,newres,diversity));
for j in range(diversity):
Pj = np.multiply(Xi_A,exp(1j*(aberration[:,:,j]+theta)));
Pj = Xi_A * exp(1j*(aberration[:,:,j]+theta));
k[:,:,j] = np.square(abs(Utilities.FFT(Pj)));
#clear Pj temp1 temp2 defocus Automatic garbage collection?
epsilon=0;
norm_rt_data=np.sqrt(sum(sum(Xi_A)));
data_zeros= [[],[],[]];
......@@ -130,7 +132,7 @@ def JWST_data_processor(config):
# initial guess
# config.u_0 = Xi_A.*exp(1i*2*pi*rand(newres));
config['u_0']= zeros((newres,newres,config['product_space_dimension']),dtype = np.complex64);
config['u_0']= zeros((newres,newres,config['product_space_dimension']),dtype = np.complex128);
for j in range(config['product_space_dimension']):
config['u_0'][:,:,j] = np.multiply(config['abs_illumination'],exp(1j*2*pi*0.5*np.ones(newres)));#rand(newres)));
config['truth'] = true_object;
......
......@@ -72,7 +72,7 @@ def magproj(u, constr):
See LukeBurkeLyon, SIREV 2002
"""
eps = 1e-23
eps = 1e-30
modsq_u = conj(u) * u
denom = modsq_u+eps
......
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