Dear Gitlab Users, for our upcoming upgrade to Gitlab v14, Gitlab will be unavailable on Thursday, 05.08.2021 from 5:00 pm to approximately 7:00 pm. Note that with v14, certain breaking changes will be introduced (https://about.gitlab.com/blog/2021/06/04/gitlab-moving-to-14-breaking-changes/).

Commit c12ed546 authored by Russell Luke's avatar Russell Luke
Browse files

Merge branch 'ElserII' into 'master'

ElserII DataIO

See merge request !15
parents 604f42cd 23bb2ee0
import pdb
import SetProxPythonPath
from proxtoolbox.experiments.phase.Elser_Experiment import Elser_Experiment
Elser = Elser_Experiment(Atoms=400, category='H', algorithm='RRR', lambda_0=0.95, lambda_max=0.95, anim=True)
Elser.run()
Elser.show()
......@@ -69,7 +69,7 @@ A concrete experiment class needs to override the
GetData.getData('Phase')
print('Loading data file CDI_intensity')
f = loadmat('../InputData/Phase/CDI_intensity.mat')
f = loadmat(datadir/'Phase'/'CDI_intensity.mat')
# diffraction pattern
dp = f['intensity']['img'][0, 0]
orig_res = max(dp.shape[0], dp.shape[1]) # actual data size
......
......@@ -6,6 +6,7 @@ from proxtoolbox.utils.graphics import addColorbar
#for downloading data
import proxtoolbox.utils.GetData as GetData
from proxtoolbox.utils.GetData import datadir
import numpy as np
from scipy.io import loadmat
......@@ -74,7 +75,7 @@ class ART_Experiment(Experiment):
# load data
if not self.silent:
print('Loading data file ART_SheppLogan.mat ')
data_Shepp = loadmat('../InputData/CT/ART_SheppLogan.mat')
data_Shepp = loadmat(datadir/'CT'/'ART_SheppLogan.mat')
N = data_Shepp['N'].item()
p = data_Shepp['p'].item()
theta = data_Shepp['theta']
......
......@@ -5,14 +5,16 @@ from proxtoolbox.experiments.orbitaltomography.planar_molecule import PlanarMole
from proxtoolbox.utils.orbitaltomog import shifted_fft, fourier_interpolate
from proxtoolbox.utils.visualization.complex_field_visualization import complex_to_rgb
#for downloading data
import proxtoolbox.utils.GetData as GetData
from proxtoolbox.utils.GetData import datadir
class DegenerateOrbital(PlanarMolecule):
@staticmethod
def getDefaultParameters():
defaultParams = {
'experiment_name': '2D ARPES',
'data_filename': '..\\Inputdata\\OrbitalTomog'
+ '\\2020_10_27_coronene_Homo1+2_ARPES_2e6counts_corrected_80x80.tif',
'data_filename': datadir/'OrbitalTomog'/'2020_10_27_coronene_Homo1+2_ARPES_2e6counts_corrected_80x80.tif',
'from_intensity_data': True,
'object': 'real',
'degeneracy': 2, # Number of degenerate states to reconstruct
......@@ -59,6 +61,10 @@ class DegenerateOrbital(PlanarMolecule):
- threshold_for_support: float, in range [0,1], fraction of the maximum at which to threshold when
determining support or support for sparsity
"""
# make sure input data can be found, otherwise download it
GetData.getData('OrbitalTomog')
super(DegenerateOrbital, self).loadData()
def createRandomGuess(self):
......
......@@ -7,6 +7,9 @@ from proxtoolbox.experiments.orbitaltomography.orbitExperiment import OrbitalTom
from proxtoolbox.utils.visualization.stack_viewer import XYZStackViewer
from proxtoolbox.utils.orbitaltomog import shifted_fft
#for downloading data
import proxtoolbox.utils.GetData as GetData
from proxtoolbox.utils.GetData import datadir
class Molecule3D(OrbitalTomographyExperiment):
@staticmethod
......@@ -56,6 +59,9 @@ class Molecule3D(OrbitalTomographyExperiment):
"""
Load data and set in the correct format for reconstruction
"""
# make sure input data can be found, otherwise download it
GetData.getData('OrbitalTomog')
# TODO: copy from data processor
raise NotImplementedError
......
......@@ -6,6 +6,9 @@ from glob import glob
from proxtoolbox.utils.h5 import write_dict_to_h5
import numpy as np
#for downloading data
import proxtoolbox.utils.GetData as GetData
from proxtoolbox.utils.GetData import datadir
class OrbitalTomographyExperiment(Experiment):
@staticmethod
......@@ -30,6 +33,10 @@ class OrbitalTomographyExperiment(Experiment):
self.output_ignore_keys = ['u1', 'u2', 'plots']
def loadData(self):
# make sure input data can be found, otherwise download it
GetData.getData('OrbitalTomog')
"""
Load or generate the dataset that will be used for
this experiment. Create the initial iterate.
......
......@@ -7,6 +7,10 @@ from proxtoolbox.experiments.orbitaltomography.planar_molecule import PlanarMole
from proxtoolbox.utils.orbitaltomog import shifted_fft, fourier_interpolate, bin_array, shifted_ifft
from proxtoolbox.utils.visualization.complex_field_visualization import complex_to_rgb
#for downloading data
import proxtoolbox.utils.GetData as GetData
from proxtoolbox.utils.GetData import datadir
class OrbitalMomentumMicroscope(PlanarMolecule):
"""
......@@ -17,8 +21,7 @@ class OrbitalMomentumMicroscope(PlanarMolecule):
def getDefaultParameters():
defaultParams = {
'experiment_name': 'Momentum microscopy',
'data_filename': '..\\Inputdata\\OrbitalTomog'
+ '\\2020_11_05_Coronene_HOMO_near_degenerate_noisefree.tif',
'data_filename': datadir/'OrbitalTomog'/'2020_11_05_Coronene_HOMO_near_degenerate_noisefree.tif',
'from_intensity_data': False,
'object': 'real',
'constraint': 'sparse real',
......@@ -54,6 +57,10 @@ class OrbitalMomentumMicroscope(PlanarMolecule):
self.momentum_axes = (-1, -2)
def loadData(self):
# make sure input data can be found, otherwise download it
GetData.getData('OrbitalTomog')
"""
Load data and set in the correct format for reconstruction
Parameters are taken from experiment class (self) properties, which must include::
......
......@@ -7,13 +7,18 @@ from proxtoolbox.experiments.orbitaltomography.planar_molecule import PlanarMole
from proxtoolbox.utils.orbitaltomog import shifted_fft, fourier_interpolate, bin_array, shifted_ifft
from proxtoolbox.utils.visualization.complex_field_visualization import complex_to_rgb
#for downloading data
import proxtoolbox.utils.GetData as GetData
from proxtoolbox.utils.GetData import datadir
class OrthogonalOrbitals(PlanarMolecule):
@staticmethod
def getDefaultParameters():
defaultParams = {
'experiment_name': '2D ARPES',
'data_filename': '..\\Inputdata\\OrbitalTomog'
'data_filename': datadir
+ '\OrbitalTomog'
+ '\\2020_10_27_coronene_Homo_stack_ARPES_2e6counts_corrected_80x80.tif',
'from_intensity_data': False,
'object': 'real',
......@@ -52,6 +57,10 @@ class OrthogonalOrbitals(PlanarMolecule):
- threshold_for_support: float, in range [0,1], fraction of the maximum at which to threshold when
determining support or support for sparsity
"""
# make sure input data can be found, otherwise download it
GetData.getData('OrbitalTomog')
# load data
if self.data_filename is None:
self.data_filename = input('Please enter the path to the datafile: ')
......
......@@ -7,6 +7,9 @@ from proxtoolbox.experiments.orbitaltomography.orbitExperiment import OrbitalTom
from proxtoolbox.utils.visualization.complex_field_visualization import complex_to_rgb
from proxtoolbox.utils.orbitaltomog import bin_array, shifted_fft, shifted_ifft, fourier_interpolate, roll_to_pos
#for downloading data
import proxtoolbox.utils.GetData as GetData
from proxtoolbox.utils.GetData import datadir
class PlanarMolecule(OrbitalTomographyExperiment):
@staticmethod
......@@ -14,7 +17,7 @@ class PlanarMolecule(OrbitalTomographyExperiment):
# TODO: optimize parameters and proxoperators to get good & consistent phase retrieval using the demo
defaultParams = {
'experiment_name': 'noisy 2D ARPES', # '2D ARPES', #
'data_filename': '../InputData/OrbitalTomog/coronene_homo1.tif',
'data_filename': datadir/'OrbitalTomog'/'coronene_homo1.tif',
'from_intensity_data': False,
'object': 'real',
'constraint': 'sparse real',
......@@ -74,6 +77,9 @@ class PlanarMolecule(OrbitalTomographyExperiment):
- threshold_for_support: float, in range [0,1], fraction of the maximum at which to threshold when
determining support or support for sparsity
"""
# make sure input data can be found, otherwise download it
GetData.getData('OrbitalTomog')
# load data
if self.data_filename is None:
self.data_filename = input('Please enter the path to the datafile: ')
......
......@@ -7,6 +7,7 @@ from proxtoolbox.utils.graphics import addColorbar
#for downloading data
import proxtoolbox.utils.GetData as GetData
from proxtoolbox.utils.GetData import datadir
import numpy as np
from numpy import exp, sqrt, log2, log10, ceil, floor, unravel_index, argmax, zeros
......@@ -77,7 +78,7 @@ class CDI_Experiment(PhaseExperiment):
GetData.getData('Phase')
if not self.silent:
print('Loading data file CDI_intensity')
f = loadmat('../InputData/Phase/CDI_intensity.mat')
f = loadmat(datadir/'Phase'/'CDI_intensity.mat')
# diffraction pattern
dp = f['intensity']['img'][0, 0]
orig_res = max(dp.shape[0], dp.shape[1]) # actual data size
......
......@@ -13,6 +13,11 @@ from numpy.linalg import norm
from numpy.fft import fft2, ifft2
import time
#for downloading data
import proxtoolbox.utils.GetData as GetData
from proxtoolbox.utils.GetData import datadir
print(datadir)
class CDP_Experiment(PhaseExperiment):
"""
......@@ -69,6 +74,9 @@ class CDP_Experiment(PhaseExperiment):
"""
Load CDP dataset. Create the initial iterate.
"""
#make sure input data can be found, otherwise download it
GetData.getData('Phase')
# Implementation of the Wirtinger Flow (WF) algorithm presented in the paper
# "Phase Retrieval via Wirtinger Flow: Theory and algorithms"
# by E. J. Candes, X. Li, and M. Soltanolkotabi
......@@ -91,18 +99,18 @@ class CDP_Experiment(PhaseExperiment):
# make image
if debug:
if n2 == 1:
x_dict = loadMatFile('../InputData/Phase/CDP_1D_x.mat')
x_dict = loadMatFile(datadir/'Phase'/'CDP_1D_x.mat')
debug_image = x_dict['x']
masks_dict = loadMatFile('../InputData/Phase/CDP_1D_Masks.mat')
masks_dict = loadMatFile(datadir/'Phase'/'CDP_1D_Masks.mat')
debug_masks = masks_dict['Masks']
z0_dict = loadMatFile('../InputData/Phase/CDP_1D_z0.mat')
z0_dict = loadMatFile(datadir/'Phase'/'CDP_1D_z0.mat')
debug_z0 = z0_dict['z0']
elif n2 == 256:
x_dict = loadMatFile('../InputData/Phase/CDP_2D_x.mat')
x_dict = loadMatFile(datadir/'Phase'/'CDP_2D_x.mat')
debug_image = x_dict['x']
masks_dict = loadMatFile('../InputData/Phase/CDP_2D_Masks.mat')
masks_dict = loadMatFile(datadir/'Phase'/'CDP_2D_Masks.mat')
debug_masks = masks_dict['Masks']
z0_dict = loadMatFile('../InputData/Phase/CDP_2D_z0.mat')
z0_dict = loadMatFile(datadir/'Phase'/'CDP_2D_z0.mat')
debug_z0 = z0_dict['z0']
x = debug_image
else:
......
import SetProxPythonPath
from proxtoolbox.experiments.phase.phaseExperiment import PhaseExperiment
from proxtoolbox import proxoperators
import numpy as np
from numpy import fromfile, exp, nonzero, zeros, pi, resize, real, angle
from numpy.random import rand
from numpy.linalg import norm
from numpy.fft import fftshift, fft2, ifft2
import proxtoolbox.utils as utils
from proxtoolbox.utils.cell import Cell, isCell
#for downloading data
import proxtoolbox.utils.GetData as GetData
from proxtoolbox.utils.GetData import datadir
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.pyplot import subplots, show, figure
class Elser_Experiment(PhaseExperiment):
'''
Elser's Crystallography experiment class
'''
@staticmethod
def getDefaultParameters():
defaultParams = {
'experiment_name' : 'Elser',
'object': 'complex',
'constraint': 'atom count',
'distance': 'far field',
'Nx': 128,
'Ny': 128,
'Nz': 1,
'Atoms':100,
'category': 'E',
'product_space_dimension': 2,
'MAXIT': 10000,
'TOL': 5e-5,
'lambda_0': 0.85,
'lambda_max': 0.85,
'lambda_switch': 20,
'data_ball': 1e-30,
'diagnostic': True,
'iterate_monitor_name': 'FeasibilityIterateMonitor',
'verbose': 0,
'graphics': 1
}
return defaultParams
def __init__(self, Atoms=100, category='E', **kwargs):
super(Elser_Experiment, self).__init__(**kwargs)
self.fresnel_nr = None
self.use_farfield_formula = None
self.data_zeros = None
self.supp_ampl = None
self.indicator_ampl = None
self.abs_illumination = None
self.illumination_phase = None
self.FT_conv_kernel = Cell(1)
self.Atoms=Atoms
self.category=category
def loadData(self):
# def loadElserData(filename, M):
"""
Load Elser dataset. Create the initial iterate.
"""
#make sure input data can be found, otherwise download it
GetData.getData('Elser')
#defaultParams = getDefaultParameters()
# So far only Nx to set the desired
# resolution
newres = self.Nx
# The following value for snr does not work well in Python
# as this creates numerical issues with Poisson noise
# (data_ball is far too small)
# snr = self.data_ball
# Instead, we choose the following bigger value which
# works well.
snr = 1e-8
print('Loading data.')
fn="".join([str(datadir/"Elser"/"data"), str(self.Atoms), self.category])
# read text data file as a matrix of numbers
data = np.loadtxt(fn, delimiter="\t")
data = np.sqrt(data)/self.Nx
data = np.c_[data, np.zeros(self.Nx)] # add extra column
print(data)
data = data.flatten() # make it a one-dimensional array
M = fftshift(data)
# have to look at what Elser is actually doing with his data...
self.magn = 1 # magnification factor
self.farfield = True
print('Using farfield formula.')
self.rt_data = Cell(1)
self.rt_data[0] = M
# standard for the main program is that
# the data field is the magnitude SQUARED
# in Luke.m this is changed to the magnitude.
self.data = Cell(1)
self.data[0] = M**2
self.norm_rt_data = np.linalg.norm(ifft2(M), 'fro')
self.data_zeros = np.where(M == 0)
self.support_idx = np.nonzero(S)
self.sets = 2
# The support constraint is sparsity determined by the number of atoms indicated in the
# data file name...this needs to be installed
# initial guess: follow Elser's initialization...
def setupProxOperators(self):
"""
Determine the prox operators to be used for this experiment
"""
super(Elser_Experiment, self).setupProxOperators() # call parent's method
self.propagator = 'Propagator_FreFra'
self.inverse_propagator = 'InvPropagator_FreFra'
# remark: self.farfield is always true for these problems
# self.proxOperators already contains a prox operator at slot 0.
# Here, we add the second one.
if self.constraint == 'phaselift':
self.proxOperators.append('P_Rank1')
elif self.constraint == 'phaselift2':
self.proxOperators.append('P_rank1_SR')
else:
self.proxOperators.append('Approx_Pphase_FreFra_Poisson')
self.nProx = self.sets
# Now we adjust data structures and prox operators according to the algorithm
# What follows wa set up for JWST and works in that context. Not
# sure how much of this is trasferrable to Elsers problems (DRL 14.08.2020)
if self.algorithm_name == 'ADMM':
# set-up variables in cells for block-wise algorithms
# ADMM has three blocks of variables, primal, auxilliary (primal
# variables in the image space of some linear mapping) and dual.
# we sort these into a 3D cell.
u0 = self.u0
self.u0 = Cell(3)
self.proxOperators = []
self.productProxOperators = []
K = self.product_space_dimension - 1
self.product_space_dimension = K
self.u0[0] = u0[K]
prox_FreFra_Poisson_class = getattr(proxoperators, 'Approx_Pphase_FreFra_Poisson')
prox_FreFra_Poisson = prox_FreFra_Poisson_class(self)
self.u0[1] = Cell(K)
self.u0[2] = Cell(K)
for j in range(K):
tmp_u0 = prox_FreFra_Poisson.eval(u0[j], j)
self.u0[1][j] = fft2(self.FT_conv_kernel[j]*tmp_u0) / (self.Nx*self.Ny)
self.u0[2][j] = self.u0[1][j] / self.lambda_0
self.proxOperators.append('Prox_primal_ADMM_indicator')
self.proxOperators.append('Approx_Pphase_FreFra_ADMM_Poisson')
elif self.algorithm_name == 'AvP2':
# u_0 should be a cell...we change it into a cell of cells
u0 = self.u0
self.u0 = Cell(3)
prox2_class = getattr(proxoperators, self.proxOperators[1])
prox2 = prox2_class(self)
self.u0[2] = prox2.eval(u0)
P_Diag_class = getattr(proxoperators, 'P_diag')
P_Diag_prox = P_Diag_class(self)
tmp_u = P_Diag_prox.eval(self.u0[2])
self.u0[1] = tmp_u[0]
self.u0[0] = u0[self.product_space_dimension-1]
elif self.algorithm_name == 'PHeBIE':
# set up variables in cells for block-wise, implicit-explicit algorithms
u0 = self.u0
self.u0 = Cell(2)
self.product_space_dimension = self.product_space_dimension - 1
self.u0[0] = u0[self.product_space_dimension].copy() # copy last element of prev u0 into first slot
self.u0[1] = Cell(self.product_space_dimension)
for j in range(self.product_space_dimension):
self.u0[1][j] = u0[j]
self.proxOperators = []
self.nProx = 2
self.proxOperators.append('P_amp_support') # project onto the support/amplitude constraints
self.proxOperators.append('Prox_product_space')
self.productProxOperators = []
self.n_product_Prox = self.product_space_dimension
for _j in range(self.n_product_Prox):
self.productProxOperators.append('Approx_Pphase_FreFra_Poisson')
elif self.algorithm_name == 'Wirtinger':
# Contrary to the Matlab code, we use cells instead of 3D arrays
L = len(self.FT_conv_kernel)
self.product_space_dimension = L
self.masks = Cell(L)
for j in range(L):
self.masks[j] = self.indicator_ampl * self.FT_conv_kernel[j]
self.data_zeros = None
self.FT_conv_kernel = self.masks
self.proxOperators[1] = 'Approx_Pphase_JWST_Wirt'
self.use_farfield_formula = True
self.fresnel_nr = 0
def show(self):
"""
Generate graphical output from the solution
"""
# display plot of far field data
# figure(123)
f, (ax1) = subplots(1, 1,
figsize=(self.figure_width, self.figure_height),
dpi=self.figure_dpi)
im = ax1.imshow(log10(self.dp + 1e-15))
f.colorbar(im, ax=ax1, fraction=0.046, pad=0.04)
ax1.set_title('Far field data')
plt.subplots_adjust(wspace=0.3) # adjust horizontal space (width)
# between subplots (default = 0.2)
f.suptitle('Elser Data')
# call parent to display the other plots
super(Elser_Experiment, self).show()
u_m = self.output['u_monitor']
if isCell(u_m):
u = u_m[0]
if isCell(u):
u = u[0]
u2 = u_m[len(u_m)-1]
if isCell(u2):
u2 = u2[len(u2)-1]
else:
u2 = u_m
if u2.ndim > 2:
u2 = u2[:,:,0]
u = self.output['u']
if isCell(u):
u = u[0]
elif u.ndim > 2:
u = u[:,:,0]
algo_desc = self.algorithm.getDescription()
title = "Algorithm " + algo_desc
# figure 904
titles = ["Best approximation amplitude - physical constraint satisfied",
"Best approximation phase - physical constraint satisfied",
"Best approximation amplitude - Fourier constraint satisfied",
"Best approximation phase - Fourier constraint satisfied"]
f = self.createFourImageFigure(u, u2, titles)
f.suptitle(title)
plt.subplots_adjust(hspace = 0.3) # adjust vertical space (height) between subplots (default = 0.2)
# figure 900
changes = self.output['stats']['changes']
time = self.output['stats']['time']
time_str = "{:.{}f}".format(time, 5) # 5 is precision
label = "Iterations (time = " + time_str + " s)"
f, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2, \
figsize = (self.figure_width, self.figure_height),
dpi = self.figure_dpi)
self.createImageSubFigure(f, ax1, abs(u), titles[0])
self.createImageSubFigure(f, ax2, real(u), titles[1])
ax3.semilogy(changes)
ax3.set_xlabel(label)
ax3.set_ylabel('Log of change in iterates')
if 'gaps' in self.output['stats']:
gaps = self.output['stats']['gaps']
ax4.semilogy(gaps)
ax4.set_xlabel(label)
ax4.set_ylabel('Log of the gap distance')
f.suptitle(title)
plt.subplots_adjust(hspace = 0.3) # adjust vertical space (height) between subplots (default = 0.2)
plt.subplots_adjust(wspace = 0.3) # adjust horizontal space (width) between subplots (default = 0.2)
show()
if __name__ == "__main__":
Elser = Elser_Experiment(Atoms=200, category='E', algorithm='RRR', lambda_0=0.95, lambda_max=0.95, anim=True)
Elser.run()
Elser.show()
load('../../../InputData/Phase/Elser/data100E')
x=data100E;
datapow=0;
M=max(size(x));
fmag=zeros(M, M/2);
Fmag=zeros(M,M);
for(i=1:M/2)
j=1;
datapow=x(i,j)+datapow;
fmag(i,j)=sqrt(x(i,j))/M;
Fmag(i,j) = fmag(i,j);
for(j=2:M/2)
datapow=2*x(i,j)+datapow;
fmag(i,j)=sqrt(x(i,j))/M;
Fmag(i,j) = fmag(i,j);
if(i>1)
Fmag(M-(j-2),M-(i-2)) = fmag(i,j);
end
end
end
for(i=M/2+1:M)
j=1;
datapow=x(i,j)+datapow;
fmag(i,