Commit f91e1f67 authored by jansen31's avatar jansen31
Browse files

Merge remote-tracking branch 'origin/dornheim' into dornheim

parents 6bcc022b ee7059c3
%% Cell type:code id: tags:
``` python
import numpy as np
from scipy.io import loadmat
%pylab inline
```
%% Output
Populating the interactive namespace from numpy and matplotlib
%% Cell type:code id: tags:
``` python
inp=loadmat('../../../InputData/OrbitalTomog/coronen_homo1_fourier_noise15.mat')
imshow(inp['I2'])
# Berechnung und Speicherung als png-Datei des Supports
# Autokorrelation
threshold_autocorr=0.1;
support0=np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(inp['I2'])));
cond2=(support0 < threshold_autocorr*np.amax(support0)).astype(uint);
subplot(121)
imshow(support0.real)
colorbar()
subplot(122)
imshow(cond2)
colorbar()
# Anfangsbedingungen
support=support0;
F0=inp['I2']**0.5;
ph_init = 2*np.pi*np.random.rand(F0.shape);
# ph_init = np.angle(np.fft.fft2(ph_init));
u0 = F0 * np.exp(1j.*ph_init);
previous = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(u0)));
```
%% Output
<matplotlib.colorbar.Colorbar at 0x7f0dcab00198>
%% Cell type:code id: tags:
``` python
imshow(inp["I2"])
```
%% Output
<matplotlib.image.AxesImage at 0x7f0dca741358>
%% Cell type:code id: tags:
``` python
inp=loadmat('../../../InputData/OrbitalTomog/coronen_homo1_fourier_noise15.mat')
%Parametereinstellung
N_ER=12;
N_HIO=0;
N_RAAR=0;
threshold_autocorr=0.1;
beta=0.9;
N_Ges=N_ER+N_HIO+N_RAAR;
n_Ges=0;
%Berechnung und Speicherung als png-Datei des Supports
%Autokorrelation
support0=fftshift(ifft2(ifftshift(I2)));
cond2=support0 < threshold_autocorr*max(support0(:));
support0=ones(size(support0));
support0(cond2)=0;
figure(1); imagesc(support0);
daspect([1 1 1]);
zoom(1);
print(gcf, '-dpng', 'support0', '-r150');
%%
%Anfangsbedingungen
support=support0;
F0=I2.^0.5;
ph_init = rand(size(F0));
ph_init = angle(fft2(ph_init));
G2 = F0 .* exp(1j.*ph_init);
previous = fftshift(ifft2(ifftshift(G2)));
%%
%Erstellung eines neuen Ordners mit Angaben ueber die gewaehlten Parameter
name_folder= sprintf('ER_%.f_HIO_%.f_RAAR_%.f_beta_%.2f',N_ER, N_HIO, N_RAAR, beta);
mkdir(name_folder)
oldFolder = cd(name_folder);
%Speichern der gewaehlten Parameter
fid = fopen('Parameter.txt', 'w');
fprintf(fid, 'N_ER\t%.f\r\nN_HIO\t%.f\r\nN_RAAR\t%.f\r\nthreshold_autocorr\t%.1f\r\nbeta\t%.2f', N_ER, N_HIO,N_RAAR,threshold_autocorr, beta);
fclose(fid);
%Der Fehler wird nach jeder Iteration in der Datei 'error.txt'
%abgespeichert
fid = fopen('error.txt', 'w');
E = zeros(N_Ges,1);
%%
%HIO-Algorithmus
```
%% Cell type:code id: tags:
``` python
imshow(np.ones((16,16)))
```
%% Output
<matplotlib.image.AxesImage at 0x7f46afb32588>
%% Cell type:code id: tags:
``` python
ff = {'a':1, 'b':3}
```
%% Cell type:code id: tags:
``` python
if 'c' not in ff:
ff['c'] = 2
ff
```
%% Output
{'a': 1, 'b': 3, 'c': 2}
%% Cell type:code id: tags:
``` python
```
......
import sys
sys.path.append('..')
from . import coronene_test
from .phase import Phase
sys.path.append('../../../')
# from .
import coronene_test # base config
import orbitaltomog_data_processor # extends config
from phase import Phase
# sys.path.append('../proxtoolbox/Problems/Phase')
orbit = Phase(coronene_test.new_config)
orbit_config = orbitaltomog_data_processor.data_processor(coronene_test.new_config)
orbit = Phase(orbit_config)
orbit.solve()
orbit.show()
\ No newline at end of file
%% Cell type:code id: tags:
``` python
import numpy as np
from scipy.io import loadmat
%pylab inline
```
%% Output
Populating the interactive namespace from numpy and matplotlib
%% Cell type:code id: tags:
``` python
inp=loadmat('../../../InputData/OrbitalTomog/coronen_homo1_fourier_noise15.mat')
imshow(inp['I2'])
# Berechnung und Speicherung als png-Datei des Supports
# Autokorrelation
threshold_autocorr=0.1;
support0=np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(inp['I2'])));
cond2=(support0 < threshold_autocorr*np.amax(support0)).astype(uint);
subplot(121)
imshow(support0.real)
colorbar()
subplot(122)
imshow(cond2)
colorbar()
# Anfangsbedingungen
support=support0;
F0=inp['I2']**0.5;
ph_init = 2*np.pi*np.random.rand(F0.shape);
# ph_init = np.angle(np.fft.fft2(ph_init));
u0 = F0 * np.exp(1j.*ph_init);
previous = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(u0)));
```
%% Output
<matplotlib.colorbar.Colorbar at 0x7f0dcab00198>
%% Cell type:code id: tags:
``` python
imshow(inp["I2"])
```
%% Output
<matplotlib.image.AxesImage at 0x7f0dca741358>
%% Cell type:code id: tags:
``` python
inp=loadmat('../../../InputData/OrbitalTomog/coronen_homo1_fourier_noise15.mat')
%Parametereinstellung
N_ER=12;
N_HIO=0;
N_RAAR=0;
threshold_autocorr=0.1;
beta=0.9;
N_Ges=N_ER+N_HIO+N_RAAR;
n_Ges=0;
%Berechnung und Speicherung als png-Datei des Supports
%Autokorrelation
support0=fftshift(ifft2(ifftshift(I2)));
cond2=support0 < threshold_autocorr*max(support0(:));
support0=ones(size(support0));
support0(cond2)=0;
figure(1); imagesc(support0);
daspect([1 1 1]);
zoom(1);
print(gcf, '-dpng', 'support0', '-r150');
%%
%Anfangsbedingungen
support=support0;
F0=I2.^0.5;
ph_init = rand(size(F0));
ph_init = angle(fft2(ph_init));
G2 = F0 .* exp(1j.*ph_init);
previous = fftshift(ifft2(ifftshift(G2)));
%%
%Erstellung eines neuen Ordners mit Angaben ueber die gewaehlten Parameter
name_folder= sprintf('ER_%.f_HIO_%.f_RAAR_%.f_beta_%.2f',N_ER, N_HIO, N_RAAR, beta);
mkdir(name_folder)
oldFolder = cd(name_folder);
%Speichern der gewaehlten Parameter
fid = fopen('Parameter.txt', 'w');
fprintf(fid, 'N_ER\t%.f\r\nN_HIO\t%.f\r\nN_RAAR\t%.f\r\nthreshold_autocorr\t%.1f\r\nbeta\t%.2f', N_ER, N_HIO,N_RAAR,threshold_autocorr, beta);
fclose(fid);
%Der Fehler wird nach jeder Iteration in der Datei 'error.txt'
%abgespeichert
fid = fopen('error.txt', 'w');
E = zeros(N_Ges,1);
%%
%HIO-Algorithmus
```
%% Cell type:code id: tags:
``` python
imshow(np.ones((16,16)))
```
%% Output
<matplotlib.image.AxesImage at 0x7f46afb32588>
%% Cell type:code id: tags:
``` python
ff = {'a':1, 'b':3}
```
%% Cell type:code id: tags:
``` python
if 'c' not in ff:
ff['c'] = 2
ff
```
%% Output
{'a': 1, 'b': 3, 'c': 2}
%% Cell type:code id: tags:
``` python
```
......
......@@ -22,8 +22,8 @@ def data_processor(config):
# Anfangsbedingungen
ph_init = 2 * np.pi * np.random.rand(ny, nx)
# ph_init = np.angle(np.fft.fft2(ph_init));
config['u0'] = inp * np.exp(1j * ph_init)
previous = fft.fftshift(fft.ifft2(fft.ifftshift(config['u0'])))
config['u_0'] = inp * np.exp(1j * ph_init)
previous = fft.fftshift(fft.ifft2(fft.ifftshift(config['u_0'])))
plt.figure(figsize=(15, 4))
plt.subplot(131)
......@@ -39,3 +39,10 @@ def data_processor(config):
plt.colorbar()
plt.title("Initial support")
plt.show()
# Other settings
config['fresnel_nr'] = 0
config['FT_conv_kernel'] = 1
config['use_farfield_formula'] = True
return config
......@@ -5,299 +5,155 @@ from proxtoolbox import Algorithms
from proxtoolbox import ProxOperators
from proxtoolbox.ProxOperators.proxoperators import ProxOperator
from proxtoolbox.Problems.Phase import Graphics
# from .proxoperators import ProxOperator
#from phase_retrieval.back_end_utilities import norm
from numpy.linalg import norm
import numpy as np
import h5py
from numpy import square, sqrt, nonzero, size
from numpy import square, sqrt
class Phase(Problem):
"""
Phase Problem
"""
config = { } # This line does not do anything.
def __init__(self, new_config={}):
config = {}
def __init__(self, new_config):
"""
The initialization of a Phase instance takes the default configuration
and updates the parameters with the arguments in new_config.
Parameters
----------
new_config : dict, optional - Parameters to initialize the problem. If unspecified, the default config is used.
new_config : dict with non-empty keys:
object, constraint, experiment, distance, algorithm
Ny, Nx
noise
MAXIT, TOL
diagnostic
data, data_norm, support,
"""
self.config = new_config
self.config = new_config
#self.back_end = back_end
#call data processor, read data
module = __import__(self.config['data_filename'])
data_processor = getattr(module, self.config['data_filename'])
data_processor(self.config)
# module = __import__(self.config['data_filename'])
# data_processor = getattr(module, self.config['data_filename'])
# data_processor(self.config)
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 method_config['formulation'] does not exist, 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 = ['','','']
# Set the projectors and inputs based on the types of constraints and experiments
used_proxoperators = ['', '', '']
# Projector 1 (real / object space)
if self.config['constraint'] == 'support only':
proxoperators[0] = 'P_S'
used_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'] == 'phase on support':
proxoperators[0] ='P_Amod'
elif self.config['constraint'] =='minimum amplitude':
proxoperators[0] = 'P_min_amp'
used_proxoperators[0] = 'P_S_real'
elif self.config['constraint'] == 'nonnegative and support':
used_proxoperators[0] = 'P_SP'
elif self.config['constraint'] == 'amplitude only':
used_proxoperators[0] = 'P_amp'
proxoperators[1]='Approx_P_FreFra_Poisson'
# Projector 2 (k / Fourier space)
# used_proxoperators[1] = 'P_M' # 'Approx_P_FreFra_Poisson'
used_proxoperators[1] = 'Approx_P_FreFra_Poisson'
self.config['proxoperators'] = []
for prox in proxoperators:
for prox in used_proxoperators:
if prox != '':
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['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
#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.
# 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(u_1)
raise NotImplementedError("This needs careful matlab -> python translation")
# 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'
proj_1 = self.config['proxoperators'][0](self.config)#, self.back_end)
u_1 = proj_1.work(self.config['u_0'])
proj_2 = self.config['proxoperators'][1](self.config)#, self.back_end)
u_2 = proj_2.work(u_1)
# 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'])
if self.config['Nx'] == 1 or self.config['Ny'] == 1: # 1D problem
tmp_gap = square(norm(u_1 - u_2) / self.config['norm_data']) # norm_rt_data
elif self.config['product_space_dimension'] == 1:
tmp_gap = (norm(u_1-u_2)/self.config['norm_rt_data'])**2
tmp_gap = (norm(u_1 - u_2) / self.config['norm_data']) ** 2 # norm_rt_data
else:
tmp_gap=0
tmp_gap = 0
for j in range(self.config['product_space_dimension']):
# compute (||P_Sx-P_Mx||/norm_data)^2:
tmp_gap = tmp_gap+(norm(u_1[:,:,j]-u_2[:,:,j])/self.config['norm_rt_data'])**2
gap_0=sqrt(tmp_gap)
tmp_gap = tmp_gap + (norm(u_1[:, :, j] - u_2[:, :, j]) / self.config['norm_data']) ** 2 # norm_rt_data
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.
self.config['data_ball'] = self.config['data_ball'] * gap_0
# the second tolerance relative to the order of magnitude of the metric
self.config['TOL2'] = self.config['data_ball'] * 1e-15
# self.config['proxoperators']
self.algorithm = getattr(algorithms, self.config['algorithm'])(self.config)#, self.back_end)
# 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.
self.config['data_ball']=self.config['data_ball']*gap_0
# the second tolerance relative to the oder of
# magnitude of the metric
self.config['TOL2'] = self.config['data_ball']*1e-15
self.config['proxoperators']
self.algorithm = getattr(Algorithms, self.config['algorithm'])(self.config)
def _presolve(self):
"""
Prepares argument for actual solving routine
"""
pass # Actually nothing for the phase problem
def _solve(self):
"""
Runs the algorithm to solve the given sudoku problem
Runs the algorithm to solve the given problem
Output can be accessed from self.output.
"""
# algorithm = self.config['algorithm'](self.config)
self.output = self.algorithm.run(self.config['u_0'],self.config['TOL'],self.config['MAXIT'])
# Call to the algorithm, specifically run in SimpleAlgorithm
self.output = self.algorithm.run(self.config['u_0'], self.config['TOL'], self.config['MAXIT'])
print('Iterations:' + str(self.output['iter']))
def _postsolve(self):
"""
Processes the solution and generates the output
"""
pass
def show(self):
"""
Generates graphical output from the solution
"""
print("Calculation time:")
print(self.elapsed_time)
graphics = getattr(Graphics,self.config['graphics_display'])
graphics(self.config,self.output)
def compare_to_matlab(self):
"""
Routine to test and verify results by comparing to matlab
Note that this is only for development and should not be used by a normal user
For result to match u_0 should be chosen as np.multiply(config['abs_illumination'],exp(1j*2*pi*0.5*np.ones(newres)))']
"""
print(self.config['proxoperators'])
if self.config['experiment'] == 'JWST':
if self.config['algorithm'] == 'RAAR':
if self.config['MAXIT'] == 1:
f = h5py.File('Phase_test_data/u1_1.mat')
elif self.config['MAXIT'] == 500 :
f = h5py.File('Phase_test_data/u1_500.mat')
else:
print("No file available to compare to.")
return
elif self.config['algorithm'] == 'AP':
f = h5py.File('Phase_test_data/JWST_u1_ap_' + str(self.config['MAXIT']) + '.mat')
u1 = f['u1'].value.view(np.complex)
elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'amplitude':
f = h5py.File('Phase_test_data/siemens_amplitude_u1_' + str(self.config['MAXIT']) + '.mat')
u1 = f['u1'].value.view(np.complex)
elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'nonnegative and support':
f = h5py.File('Phase_test_data/siemens_nonneg_u1_' + str(self.config['MAXIT']) + '.mat')
u1 = f['u1'].value.view(np.float64)
elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'real and support':
f = h5py.File('Phase_test_data/siemens_real_u1_' + str(self.config['MAXIT']) + '.mat')
u1 = f['u1'].value.view(np.float64)
else:
if self.config['algorithm'] == 'RAAR':
if self.config['beta_0'] == 0.95:
if self.config['MAXIT'] == 1000 :
f = h5py.File('Phase_test_data/tasse_u1_1000.mat')
elif self.config['MAXIT'] == 20:
f = h5py.File('Phase_test_data/tasse_u1_20.mat')
elif self.config['MAXIT'] == 1:
f = h5py.File('Phase_test_data/tasse_u1_1.mat')
else:
print("No file available to compare to.")
return
elif self.config['beta_0'] == 0.50:
f = h5py.File('Phase_test_data/tasse_u1_'+ str(self.config['MAXIT']) + '_beta_0_5.mat')
else:
print("No file available to compare to.")
return
elif self.config['algorithm'] == 'AP' and self.config['constraint'] == 'support only':
f = h5py.File('Phase_test_data/tasse_supp_u1_ap_' + str(self.config['MAXIT']) + '.mat')
elif ( self.config['algorithm'] == 'AP' or self.config['algorithm'] == 'AP_expert') and self.config['constraint'] == 'nonnegative and support':
f = h5py.File('Phase_test_data/tasse_u1_ap_' + str(self.config['MAXIT']) + '.mat')
u1 = f['u1'].value.view(np.float64)
u1 =np.array(u1)
u1 = u1.T
print("Compare u1:")
#print("Nonzero indices matlab:")
#print(nonzero(u1))
#print("Nonzero indices python:")
#print(nonzero(self.output['u1']))
print("Nonzero indices equal:")
print(np.array_equal(nonzero(u1),nonzero(self.output['u1'])))
#print("Nonzero values matlab:")
#print(u1[nonzero(u1)])
#print("Nonzero values python:")
#print(self.output['u1'][nonzero(self.output['u1'])])
#print("Difference at nonzero values:")
#nonz = nonzero(u1)
diff = u1 - self.output['u1']
#print(diff[nonz])
print("Maximum norm of difference:")
print(np.amax(abs(diff)))
print("Frobenius norm of difference:")
print(norm(diff))
print("Frobenius norm of matlab u1:")
print(norm(u1))
print("Frobenius norm of python u1:")
print(norm(self.output['u1']))
def compare_data_to_matlab(self):
"""
Routine to test and verify results by comparing to matlab
Note that this is only for development and should not be used by a normal user
For result to match u_0 should be chosen as np.multiply(config['abs_illumination'],exp(1j*2*pi*0.5*np.ones(newres)))'] =
"""
if self.config['data_filename'] == 'CDI_data_processor':
f = h5py.File('Phase_test_data/CDI_data_processor_rt_data.mat')
data_mat = f['rt_data'].value.view(np.float)
data_py = self.config['rt_data']
elif self.config['data_filename'] == 'Siemens_processor':
f = h5py.File('Phase_test_data/Siemens_processor_truth.mat')
data_mat = f['truth'].value.view(np.float)
data_py = self.config['truth']