Commit 4a9d315f authored by alexander.dornheim's avatar alexander.dornheim
Browse files

Do some precalculations in Approx_P_JWST_Poisson for better performance. This...

Do some precalculations in Approx_P_JWST_Poisson for better performance. This modification results in ~20% less computation time for JWST (500 Iterations).
parent b8068182
......@@ -7,5 +7,5 @@ from phase import Phase
JWST = Phase(JWST_in.new_config)
JWST.solve()
JWST.compare_to_matlab()
JWST.show()
#JWST.compare_to_matlab()
#JWST.show()
......@@ -87,7 +87,7 @@ new_config = {
#maximum number of iterations and tolerances
'MAXIT' : 1,
'MAXIT' : 500,
'TOL' : 1e-8,
# relaxaton parameters in RAAR, HPR and HAAR
......
......@@ -8,7 +8,7 @@ The "ProxOperators"-module contains various specific operators that do the actua
"""
import numpy as np
from numpy import conj, dot, empty, ones, sqrt, sum, zeros, exp, nonzero, log
from numpy import conj, dot, empty, ones, sqrt, sum, zeros, exp, nonzero, log, tile
from numpy.fft import fft2, ifft2
__all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson", "P_amp"]
......@@ -257,6 +257,12 @@ class Approx_P_JWST_Poisson(ProxOperator):
self.product_space_dimension = config['product_space_dimension'];
self.abs_illumination = config['abs_illumination'];
#calculate the following once for better performance (speedup of JWST ~20% for 500 Iterations)
self.exp_illu = exp(1j*self.illumination_phase)*tile(self.indicator_ampl[...,None],(1,1,self.product_space_dimension-1))
self.exp_illu_minus = exp(-1j*self.illumination_phase)*tile(self.indicator_ampl[...,None],(1,1,self.product_space_dimension-1))
#@profile
def work(self,u):
TOL2 = self.TOL2;
epsilon = self.epsilon;
......@@ -267,12 +273,14 @@ class Approx_P_JWST_Poisson(ProxOperator):
data = self.data;
product_space_dimension = self.product_space_dimension;
abs_illumination = self.abs_illumination
exp_illu = self.exp_illu
exp_illu_minus = self.exp_illu_minus
u_new = zeros(u.shape,u.dtype)
for j in range(product_space_dimension-1):
U = fft2( indicator_ampl *exp(1j*illumination_phase [:,:,j]) * u[:,:,j]);
U = fft2( exp_illu[:,:,j] * u[:,:,j]);
U_sq = U * conj(U);
tmp = U_sq/data_sq[:,:,j];
mask = ones(tmp.shape, np.bool);
......@@ -285,7 +293,7 @@ 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_new[:,:,j] = indicator_ampl *exp(-1j*illumination_phase[:,:,j]) *ifft2(U0);
u_new[:,:,j] = exp_illu_minus[:,:,j] *ifft2(U0);
else:
u_new[:,:,j] = u[:,:,j]
#else:
......
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