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

Made some progress in phase.py on initialization

Added ProxOperator Approx_P_JWST_Poisson: still needs to be tested
parent 4cf66568
......@@ -4,7 +4,5 @@ import JWST_data_processor
import JWST_in
from phase import Phase
print(JWST_in.new_config);
JWST = Phase(JWST_in.new_config);
JWST.solve();
......@@ -23,8 +23,123 @@ class Phase(Problem):
new_config : dict, optional - Parameters to initialize the problem. If unspecified, the default config is used.
"""
self.config.update(new_config)
print(new_config['algorithm']);
self.algorithm = getattr(Algorithms, self.config['algorithm'])(new_config)
#call data processor, read data
module = __import__(self.config['data_filename'])
data_processor = getattr(module, self.config['data_filename'])
data_processor(self.config)
#reshape and rename the data
self.config['data_sq'] = self.config['data'];
self.config['data'] = self.config['rt_data'];
tmp = self.config['data'].shape;
if(tmp[0]==1 or tmp[1]==1):
self.config['data_sq'] = self.config['data_sq'].reshape((self.config['Nx'],self.config['Ny']));
#the projection algorithms work with the square root of the measurement:
self.config['data'] = self.config['data'].reshape((self.config['Nx'],self.config['Ny']));
self.config['norm_data']=self.config['norm_rt_data'];
if 'Nz' not in self.config:
self.config[Nz] = 1;
#If method_config[formulation is does not exist, i.e. not specified in
#the *_in.m file, use the product space as the default.
if 'formulation' in self.config:
formulation = self.config['formulation'];
else:
formulation = 'product space';
# Set the projectors and inputs based on the types of constraints and
# experiments
proxoperators = ['','',''];
if self.config['constraint'] == 'hybrid':
proxoperators[0] = 'P_cP'; # This will be problem specific
elif self.config['constraint'] == 'support only':
proxoperators[0] = 'P_S';
elif self.config['constraint'] == 'real and support':
proxoperators[0] ='P_S_real';
elif self.config['constraint'] =='nonnegative and support':
proxoperators[0] ='P_SP';
elif self.config['constraint'] =='amplitude only':
proxoperators[0] ='P_amp';
elif self.config['constraint'] =='minimum amplitude':
proxoperators[0] = 'P_min_amp';
elif self.config['constraint'] =='sparse':
proxoperators[0] = 'not in yet';
elif self.config['constraint'] =='phaselift':
proxoperators[0] = 'P_mean_SP';
elif self.config['constraint'] =='phaselift2':
proxoperators[0] ='P_liftM';
proxoperators[2] ='Approx_PM_Poisson'; # Patrick: This is just to monitor the change of phases!
if self.config['experiment'] == 'single diffraction':
if self.config[distance] == 'far field':
if self.config['constraint'] == 'phaselift':
proxoperators[1] = 'P_Rank1';
elif self.config['constraint'] == 'phaselift2':
proxoperators[1] = 'P_rank1_SR';
else:
if self.config['noise'] == 'Poisson':
proxoperators[1] ='Approx_PM_Poisson';
else:
proxoperators[1] ='Approx_PM_Gaussian';
else:
proxoperators[1]='P_Fresnel';
# The following selects the projectors for diversity diffraction not
# performed in the product space. So far only used for RCAAR.
elif self.config['experiment'] == 'diversity diffraction' and formulation == 'sequential':
proxoperators[1] = 'Approx_P_RCAAR_JWST_Poisson';
proxoperators[0] = proxoperators[1];
elif self.config['experiment'] == 'JWST':
proxoperators[1] = 'Approx_P_JWST_Poisson';
proxoperators[2] = proxoperators[0];
proxoperators[0] = 'P_diag';
elif self.config['experiment'] == 'CDP':
proxoperators[1] = 'P_CDP';
proxoperators[2] = proxoperators[0];
proxoperators[0] = 'P_diag';
elif self.config['experiment'] == 'ptychography':
proxoperators[1] = 'not in yet';
elif self.config['experiment'] == 'complex':
proxoperators[1] = 'not in yet';
elif self.config['constraint'] == 'phaselift':
proxoperators[1] ='P_PL_lowrank';
self.config['proxoperators'] = [];
for prox in proxoperators:
self.config['proxoperators'].append(getattr(ProxOperators, prox))
# input.Proj1_input.F=F; % is it any more expensive to pass everything
# into the projectors rather than just a selection of data and
# parameters? If not, and we pass everything anyway, there is no need
# to create a new structure element.
if 'product_space_dimension' not in self.config:
self.config['product_space_dimension'] = 1;
# set the animation program:
self.config['animation']='Phase_animation';
#
# if you are only working with two sets but
# want to do averaged projections
# (= alternating projections on the product space)
# or RAAR on the product space (=swarming), then
# you will want to change product_space_dimension=2
# and adjust your input files and projectors accordingly.
# you could also do this within the data processor
self.config['TOL2'] = 1e-15;
print(self.config);
self.algorithm = getattr(Algorithms, self.config['algorithm'])(self.config);
......
......@@ -9,4 +9,4 @@ The "ProxOperators"-module contains all proximal operators that are already impl
from .proxoperators import *
__all__ = ["P_diag","P_parallel","magproj"]
__all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson"]
......@@ -9,7 +9,7 @@ The "ProxOperators"-module contains various specific operators that do the actua
from numpy import conj, dot, empty, ones, sqrt, sum, zeros
__all__ = ["P_diag","P_parallel","magproj"]
__all__ = ["P_diag","P_parallel","magproj", "Approx_P_JWST_Poisson"]
......@@ -207,3 +207,59 @@ class P_parallel(ProxOperator):
return u_parallel
#original matlab comment
# Approx_PM_Poisson.m
# written on Feb. 18, 2011 by
# Russell Luke
# Inst. Fuer Numerische und Angewandte Mathematik
# Universitaet Gottingen
#
# DESCRIPTION: Projection subroutine for projecting onto Fourier
# magnitude constraints
#
# INPUT: input = data structure
# .data_ball is the regularization parameter described in
# D. R. Luke, Nonlinear Analysis 75 (2012) 1531–1546.
# u = function in the physical domain to be projected
#
# OUTPUT:
#
# USAGE: u_epsilon = P_M(input,u)
#
# Nonstandard Matlab function calls: MagProj
class Approx_P_JWST_Poisson(ProxOperator):
def __init__(self,config):
self.TOL2 = config['TOL2'];
self.epsilon = config['data_ball'];
def work(self,u):
for j in range(product_space_dimension-1):
U = fft2( config['indicator_ampl'] *exp(1j*confg['illumination_phase'][:,:,j]) * u[:,:,j]);
U_sq = U * conj(U);
tmp = U_sq/config['data_sq'][:,:,j];
tmp2= config['data_zeros'][:,j] != 0;
tmpI = config['data_zeros'][tmp2,j];
tmp[tmpI]=1;
U_sq[tmpI]=0;
IU= tmp==0;
tmp[IU]=1;
tmp=log(tmp);
hU = sum(sum(U_sq *tmp + config['data_sq'][:,:,j] - U_sq));
if hU>=epsilon+TOL2:
U0 = MagProj(input.data[:,:,j],U);
u[:,:,j] = config['indicator_ampl'] *exp(-1j*config['illumination_phase'][:,:,j]) *ifft2(U0);
#else:
# no change
# now project onto the pupil constraint.
# this is a qualitative constraint, so no
# noise is taken into account.
j=input.product_space_dimension;
u[:,:,j] = magproj(config['abs_illumination'],u[:,:,j]);
return u;
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