Commit adffb0d3 authored by dinh's avatar dinh
Browse files

to check

parent ffd00a66
import sys, os
import SetProxPythonPath
from proxtoolbox.experiments.orbitaltomography import molecule_3d as m3
import random
exp_m3 = m3.Molecule3D(experiment='3D ARPES', # 'noisy 2D ARPES' OR '2D ARPES'
constraint='sparse real',
sparsity_parameter=790,
algorithm = 'AvP',
formulation = 'product space',
crop_data=(0, 4, 4),
use_sparsity_with_support=False,
TOL=1e-14, MAXIT=500000,
rnd_seed=42, verbose=0,
progressbar='tqdm')
#exp_m3.plotInputData()
exp_m3.run()
exp_m3.saveOutput(savedir='./output/', name='AvP')
exp_m3.show(zoom='auto')
rand = random.sample(range(2**0,2**28),150)
for randm in rand:
print('seed = ',randm)
exp_m3 = m3.Molecule3D( data_dir = "/home/dinh/Desktop/PhD_topic/ot3d/Notebook_3d_tomography/created_data/0.10/",
data_filename = 'measured_kspace.tif',
sensitivity_filename = 'sensitivity_lowpass.tif',
exact_filename = 'exact.tif',
experiment='3D ARPES',
algorithm = 'AvP',
constraint= 'sparse real',
sparsity_parameter=1500,
crop_data=(0, 0, 0),
formulation='product space',
use_sparsity_with_support=False,
TOL=1e-14, MAXIT=200,
rnd_seed=randm, verbose=0,
progressbar='tqdm')
# exp_m3.plotInputData()
#exp_m3.show_3d(inpt=True, fft=True) #3d simulation of input in Fourier space
exp_m3.run()
# exp_m3.saveOutput(savedir='./output/', name='AvP')
exp_m3.show() #convergence plot
exp_m3.show_sclices(zoom='auto') #slices
exp_m3.show_3d(inpt= False,fft=False) #3d simulation in real space
exp_m3.show_3d(inpt=False, fft=True, Real=True) #3d simulation in Fourier space, real part
exp_m3.show_3d(inpt=False, fft=True, Real=False) #3d simulation in Fourier space, imaginary part
break
import sys, os
import SetProxPythonPath
from proxtoolbox.experiments.orbitaltomography import molecule_3d as m3
import random
exp_m3 = m3.Molecule3D(experiment='3D ARPES', # 'noisy 2D ARPES' OR '2D ARPES'
constraint='sparse real',
sparsity_parameter=203161,
algorithm = 'AvP2',
formulation='product space',
crop_data=(0, 4, 4),
use_sparsity_with_support=True,
threshold_for_support = 0.17,
TOL=1e-14, MAXIT=100,
lambda_0=0.05, lambda_max=0.1,
rnd_seed=42, verbose=0,
progressbar='tqdm')
rand = random.sample(range(2**0,2**28),150)
for randm in rand:
print('seed = ',randm)
exp_m3 = m3.Molecule3D( data_dir = "/home/dinh/Desktop/PhD_topic/ot3d/Notebook_3d_tomography/created_data/0.10/",
data_filename = 'measured_kspace.tif',
sensitivity_filename = 'sensitivity_lowpass.tif',
exact_filename = 'exact.tif',
experiment='3D ARPES',
algorithm = 'AvP2',
constraint= 'sparse real',
sparsity_parameter=1500,
crop_data=(0, 0, 0),
formulation='product space',
use_sparsity_with_support=False,
TOL=1e-14, MAXIT=200,
lambda_0=0.05, lambda_max=0.1,
rnd_seed=42, verbose=0,
progressbar='tqdm')
#exp_m3.plotInputData()
exp_m3.run()
exp_m3.saveOutput(savedir='./output/', name='AvP2')
......
......@@ -2,28 +2,35 @@
from numpy.core.fromnumeric import diagonal
import SetProxPythonPath
from proxtoolbox.experiments.orbitaltomography import molecule_3d as m3
import random
exp_m3 = m3.Molecule3D( data_dir = "./InputData/OrbitalTomog/3d_arpes_reconstruction/",
experiment='3D ARPES',
constraint= 'sparse real',
sparsity_parameter=293,
algorithm = 'CP',
formulation ='cyclic',
crop_data=(4, 4, 4),
use_sparsity_with_support=False,
TOL=1e-14, MAXIT=500000,
rnd_seed=42, verbose=0,
progressbar='tqdm')
#exp_m3.plotInputData()
#exp_m3.show_3d(inpt=True, fft=True) #3d simulation of input in Fourier space
exp_m3.run()
exp_m3.saveOutput(savedir='./output/', name='CP')
exp_m3.show() #convergence plot
exp_m3.show_sclices(zoom='auto') #slices
exp_m3.show_3d(inpt= False,fft=False) #3d simulation in real space
exp_m3.show_3d(inpt=False, fft=True) #3d simulation in Fourier space
for s in range(500,3000,50):
exp_m3 = m3.Molecule3D( data_dir = "/home/dinh/Desktop/PhD_topic/ot3d/Notebook_3d_tomography/created_data/0.10/",
data_filename = 'measured_kspace.tif',
sensitivity_filename = 'sensitivity_lowpass.tif',
exact_filename = 'exact.tif',
experiment='3D ARPES',
# constraint= 'sparse real',
constraint= 'symmetric sparse real',
sparsity_parameter=s,
algorithm = 'CP',
formulation ='cyclic',
crop_data=(0,0,0),
truth = True,
use_sparsity_with_support=True,
threshold_for_support = 0.0006,
TOL=1e-14, MAXIT=500000,
rnd_seed=42,
verbose=0,
progressbar='tqdm')
# exp_m3.plotInputData()
# exp_m3.show_3d(inpt=True, fft=True) #3d simulation of input in Fourier space
exp_m3.run()
exp_m3.saveOutput(savedir='./output/', name='CP_0.10_S_3hemisph20-24-29_1e5*7')
# exp_m3.show() #convergence plot
# exp_m3.show_sclices(zoom='auto') #slices
# exp_m3.show_3d(inpt= False,fft=False) #3d simulation in real space
# exp_m3.show_3d(inpt=False, fft=True, Real=True) #3d simulation in Fourier space, real part
# break
\ No newline at end of file
import sys, os
import SetProxPythonPath
from proxtoolbox.experiments.orbitaltomography import molecule_3d as m3
import random
exp_m3 = m3.Molecule3D( data_dir = "./InputData/OrbitalTomog/3d_arpes_reconstruction/",
experiment='3D ARPES',
constraint= 'sparse real',
sparsity_parameter=270,
algorithm = 'DRl',
crop_data=(4, 4, 4),
formulation = 'cyclic',
use_sparsity_with_support=False,
TOL=1e-14, MAXIT=6000,
lambda_0=0.1, lambda_max=0.33, lambda_switch=30,
rnd_seed=42, verbose=0,
progressbar='tqdm')
#exp_m3.plotInputData()
#exp_m3.show_3d(inpt=True, fft=True) #3d simulation of input in Fourier space
exp_m3.run()
exp_m3.saveOutput(savedir='./output/', name='DRl')
exp_m3.show() #convergence plot
exp_m3.show_sclices(zoom='auto') #slices
exp_m3.show_3d(inpt= False,fft=False) #3d simulation in real space
exp_m3.show_3d(inpt=False, fft=True) #3d simulation in Fourier space
rand = random.sample(range(2**0,2**28),150)
for randm in rand:
print('seed = ',randm)
exp_m3 = m3.Molecule3D( data_dir = "/home/dinh/Desktop/PhD_topic/ot3d/Notebook_3d_tomography/created_data/0.50/",
data_filename = 'measured_kspace.tif',
sensitivity_filename = 'sensitivity_lowpass.tif',
exact_filename = 'c22h14_HOMO.cub',
experiment='3D ARPES',
#constraint= 'sparse real',
constraint= 'symmetric sparse real',
sparsity_parameter=1500,
crop_data=(0, 0, 0),
formulation = 'cyclic',
use_sparsity_with_support=True,
threshold_for_support = 0.025,
truth = True,
TOL=2.6e-5, MAXIT=200,
lambda_0=0.1, lambda_max=0.33, lambda_switch=30,
rnd_seed=randm, verbose=0,
progressbar='tqdm')
#exp_m3.plotInputData()
# #exp_m3.show_3d(inpt=True, fft=True) #3d simulation of input in Fourier space
exp_m3.run()
#exp_m3.saveOutput(savedir='./output/', name='DRl')
exp_m3.show() #convergence plot
exp_m3.show_sclices(zoom='auto') #slices
exp_m3.show_3d(inpt= False,fft=False) #3d simulation in real space
exp_m3.show_3d(inpt=False, fft=True, Real=True) #3d simulation in Fourier space, real part
exp_m3.show_3d(inpt=False, fft=True, Real=False) #3d simulation in Fourier space, imaginary part
break
......@@ -228,7 +228,7 @@ class Algorithm:
iteration = self.iter + 1 # add 1 to conform with matlab version
# gradually changes lambda from lambda_0 into lambda_max
a = -iteration / self.lambda_switch
a *= a * a # compute a^3
a *= a * a
s = exp(a)
lmbda = s * self.lambda_0 + (1 - s) * self.lambda_max
return lmbda
from numpy import angle, trace, exp, sqrt, matmul, reshape,sum
from numpy import angle, trace, exp, sqrt, matmul, reshape,sum, unique, max, min
from numpy.linalg import norm
import numpy as np
from proxtoolbox.algorithms.feasibilityIterateMonitor import FeasibilityIterateMonitor
from proxtoolbox.utils.cell import Cell, isCell
from proxtoolbox.utils.orbitaltomog.gsmj_tools.arrays import error_measurment
from proxtoolbox.utils.size import size_matlab
from proxtoolbox.utils.orbitaltomog.gsmj_tools import arrays
class Orbit3DIterateMonitor(FeasibilityIterateMonitor):
......@@ -18,9 +20,12 @@ class Orbit3DIterateMonitor(FeasibilityIterateMonitor):
super(Orbit3DIterateMonitor, self).__init__(experiment)
self.algorithm_name = experiment.algorithm_name
self.productProxOperators = experiment.productProxOperators
if self.truth is not None:
self.exact = experiment.exact
self.norm_exact = experiment.norm_exact
self.gaps = None
self.shadow_changes = None
self.rel_errors = None
self.error = None
self.n_prox = len(self.proxOperators)
......@@ -64,8 +69,8 @@ class Orbit3DIterateMonitor(FeasibilityIterateMonitor):
self.shadow_changes = self.changes.copy()
self.shadow_changes[0] = sqrt(999)
if self.truth is not None:
self.rel_errors = self.changes.copy()
self.rel_errors[0] = sqrt(999)
self.error = self.changes.copy()
self.error[0] = sqrt(999)
def updateStatistics(self, alg):
"""
......@@ -115,6 +120,7 @@ class Orbit3DIterateMonitor(FeasibilityIterateMonitor):
tmp_change = 0
tmp_gap = 0
tmp_shadow = 0
tmp_error = 0
if isCell(self.u_monitor[0]):
for j in range(len(self.u_monitor[0])):
......@@ -126,18 +132,34 @@ class Orbit3DIterateMonitor(FeasibilityIterateMonitor):
for j in range(len(self.u_monitor[0])):
tmp_gap = tmp_gap + sum(abs(self.u_monitor[k][j]-u1[j])**2)/self.norm_data**2
tmp_shadow = tmp_shadow + sum(abs(self.u_monitor[k][j]-prev_u_mon[k][j])**2)/self.norm_data**2
if self.truth is not None:
pred = self.u_monitor[1]
tmp_error += error_measurment.error(self.exact,pred)
# phase_offset = angle(sum(self.exact * pred.conj()))
# tmp_error = tmp_error + sum(abs((self.exact - pred * exp(1j * phase_offset)) ** 2)) / self.norm_data ** 2
# #tmp_error = tmp_error + np.sum(abs(self.exact-pred)**2)/self.norm_data**2
else:
tmp_change = tmp_change + sum(abs(self.u_monitor[0]-prev_u_mon[0])**2)/self.norm_data**2
if self.diagnostic:
tmp_shadow = tmp_change
tmp_error = tmp_change
for k in range(1, nProx):
# compute (||P_SP_Mx-P_Mx||/normM)^2:
tmp_gap = tmp_gap + sum(abs(self.u_monitor[k]-u1)**2)/self.norm_data**2
tmp_shadow = tmp_shadow + sum(abs(self.u_monitor[k]-prev_u_mon[k])**2)/self.norm_data**2
if self.truth is not None:
pred = self.u_monitor[1]
tmp_error += error_measurment.error(self.exact,pred)
# phase_offset = angle(sum(self.exact * pred.conj()))
# tmp_error = tmp_error + sum(abs((self.exact - pred * exp(1j * phase_offset)) ** 2)) / self.norm_data ** 2
# # tmp_error = tmp_error + np.sum(abs(self.exact-pred)**2)/self.norm_data**2
self.changes[alg.iter] = sqrt(tmp_change)
#print(self.changes[alg.iter])
if self.diagnostic:
if self.truth is not None:
self.error[alg.iter] = sqrt(tmp_error)
self.gaps[alg.iter] = sqrt(tmp_gap)
self.shadow_changes[alg.iter] = sqrt(tmp_shadow) # this is the Euclidean norm of the gap to
# the unregularized set. To monitor the Euclidean norm of the gap to the
......@@ -184,8 +206,7 @@ class Orbit3DIterateMonitor(FeasibilityIterateMonitor):
if self.diagnostic:
out = self.proxOperators[1].eval(output['u_monitor'])
output['u1'] = out[0] # projection on S_real
output['u2'] = out[1] # projection on M_masked
output['u2'] = out[1] # projection on M_masked
output['stats']['changes'] = self.changes[1:alg.iter + 1]
if self.diagnostic:
......@@ -193,5 +214,5 @@ class Orbit3DIterateMonitor(FeasibilityIterateMonitor):
stats['gaps'] = self.gaps[1:alg.iter + 1]
stats['shadow_changes'] = self.shadow_changes[1:alg.iter + 1]
if self.truth is not None:
stats['rel_errors'] = self.rel_errors[1:alg.iter + 1]
stats['error'] = self.error[1:alg.iter + 1]
return output
\ No newline at end of file
......@@ -8,18 +8,16 @@ import h5py
from glob import glob
from scipy import sparse
from matplotlib.pyplot import show,subplots
from skimage.transform import rescale
from mayavi import mlab
from proxtoolbox.utils.cell import Cell
from proxtoolbox.experiments.orbitaltomography.orbitExperiment import OrbitalTomographyExperiment
from proxtoolbox.utils.visualization.stack_viewer import XYZStackViewer
from proxtoolbox.utils.orbitaltomog import shifted_fft
from proxtoolbox.utils.orbitaltomog import mirror_kspace
from proxtoolbox.utils.orbitaltomog import shifted_fft,mirror_kspace, binning
from proxtoolbox import proxoperators
from proxtoolbox.experiments.orbitaltomography.planar_molecule import support_from_autocorrelation
from proxtoolbox.utils.orbitaltomog.gsmj_tools import fourier_interpolate,roll_to_pos
from proxtoolbox.utils.orbitaltomog.gsmj_tools import fourier_interpolate, roll_to_pos, bin_array, load_cube_file,arrays
from proxtoolbox.utils import GetData
from proxtoolbox.utils.orbitaltomog.gsmj_tools.io import load_cube_file
from proxtoolbox.utils.orbitaltomog.gsmj_tools import roll_to_pos,bin_array
from proxtoolbox.utils.h5 import write_dict_to_h5
from proxtoolbox.utils.cell import Cell, isCell
......@@ -32,7 +30,7 @@ class Molecule3D(OrbitalTomographyExperiment):
'object': 'real',
'constraint': 'sparse real',
'sparsity_parameter': 180,
'use_sparsity_with_support': True, # Apply a loose support
'use_sparsity_with_support': False, # Apply a loose support
'symmetry_type': None,
'symmetry_axis': 0,
'algorithm': 'DRl',
......@@ -49,6 +47,7 @@ class Molecule3D(OrbitalTomographyExperiment):
'from_intensity_data': False, # File gives field amplitudes
'sensitivity_filename': 'sensitivity_lowpass.tif', # Measurement sensitivity file
'support_filename': None,
'exact_filename':None,
'noise':False,
'MAXIT': 500,
'TOL': 1e-10,
......@@ -63,6 +62,7 @@ class Molecule3D(OrbitalTomographyExperiment):
'zoomin_on_result': True, # Default zoom in to 50% of the field of view
'rotate': False,
'truth' : None,
'rnd_seed': 42,
'verbose': 1,
'graphics': 1,
'anim': False,
......@@ -82,10 +82,7 @@ class Molecule3D(OrbitalTomographyExperiment):
self.from_intensity_data = kwargs['from_intensity_data']
self.sensitivity_filename = kwargs['sensitivity_filename']
self.support_filename = kwargs['support_filename'] # optional class argument
self.threshold_for_support = kwargs['threshold_for_support']
self.sparsity_parameter = kwargs['sparsity_parameter']
self.symmetry_type = kwargs['symmetry_type']
self.symmetry_axis = kwargs['symmetry_axis']
self.noise = kwargs['noise']
self.use_sparsity_with_support = kwargs['use_sparsity_with_support']
self.threshold_for_support = kwargs['threshold_for_support'] # to determine support from the autocorrelation.
......@@ -100,12 +97,16 @@ class Molecule3D(OrbitalTomographyExperiment):
self.progressbar = kwargs['progressbar']
self.diagnostic = kwargs['diagnostic']
self.iterate_monitor_name = kwargs['iterate_monitor_name']
self.truth = kwargs['truth']
self.exact_filename = kwargs['exact_filename']
self.rnd_seed = kwargs['rnd_seed']
# the following data members are set by loadData(), in addition to those specified in parent classes
self.sensitivity = None # measurement sensitivity distribution
self.support = None # support
self.sparsity_support = None
self.fmask = None # Fourier mask - masking points in Fourier space
self.data_zeros = None # Necessary for Approx operators
self.exact = None # Exact of uper part molecule
# For self.plotInputData
self.input_plots = []
self.input_ignore_attrs += ['input_plots'] # TODO: others?
......@@ -116,21 +117,27 @@ class Molecule3D(OrbitalTomographyExperiment):
"""
print('Loading data.')
inp = imread(f"{self.data_dir}{self.data_filename}")
inp_shape = inp.shape
# Load measurement sensitivity, determine data_zeros
sensitivity = imread(f"{self.data_dir}{self.sensitivity_filename}")
if self.truth is not None:
homo = imread(f"{self.data_dir}{self.exact_filename}")
# homo,x1,y1,z1 = load_cube_file(f"{self.data_dir}{self.exact_filename}",pad_result=True)
# homo = rescale(homo,(inp_shape[2]/200,inp_shape[1]/200,inp_shape[0]/200*2)) # order (z,x,y) since the measured data has been transposed in measurement process.
if self.crop_data != (0, 0, 0):
cy, cx, cz = self.crop_data
ny, nx, nz = inp.shape
inp = inp[cy:ny - cy, cx:nx - cx, cz:nz - cz]
sensitivity = sensitivity[cy:ny - cy, cx:nx - cx, cz:nz - cz]
cz, cy, cx = self.crop_data
nz, ny, nx = inp.shape
inp = inp[cz:nz - cz, cy:ny - cy, cx:nx - cx]
sensitivity = sensitivity[cz:nz - cz, cy:ny - cy, cx:nx - cx]
if self.truth is not None:
homo = homo[cz:nz - cz, cy:ny - cy, cx:nx - cx]
# Add the negative K_z part to determine a real-valued orbital?
if self.add_negative_kz:
inp = np.moveaxis(inp, self.kz_dimension, 2)
sensitivity = np.moveaxis(sensitivity, self.kz_dimension, 2)
inp, sensitivity = mirror_kspace(inp, sensitivity, square_array=True)
inp = np.moveaxis(inp, 2, self.kz_dimension)
sensitivity = np.moveaxis(sensitivity, 2, self.kz_dimension)
ny, nx, nz = inp.shape
inp_upper = inp[::-1, :, :]
sensitivity_upper = sensitivity[::-1, :, :]
inp = np.vstack((inp_upper,inp))
sensitivity = np.vstack((sensitivity_upper,sensitivity))
nz, ny, nx = inp.shape
snr = 1e-2
if self.noise:
noise_mask = np.random.poisson(inp)*snr
......@@ -138,6 +145,8 @@ class Molecule3D(OrbitalTomographyExperiment):
self.data = inp
self.sensitivity = sensitivity
self.fmask = (self.sensitivity == 0)
if self.truth is not None:
self.exact = homo
self.sets = 2
# Keep the same resolution?
if self.Nx is None:
......@@ -165,6 +174,8 @@ class Molecule3D(OrbitalTomographyExperiment):
# Calculate the norm of the data, estimating a correction for the lack of data given by data_zeros
self.norm_data = np.sqrt(np.sum(self.data ** 2))
self.norm_data /= np.mean((1 - self.fmask))
if self.truth is not None:
self.norm_exact = np.sqrt(np.sum(self.exact ** 2))
# Object support determination
if self.constraint in ['support only', 'real and support',
'nonnegative and support'] or self.use_sparsity_with_support:
......@@ -188,7 +199,7 @@ class Molecule3D(OrbitalTomographyExperiment):
if self.formulation=='cyclic':
self.ph_init = 2 * np.pi * np.random.random_sample(self.data.shape)
self.u0 = abs(self.data) * np.exp(1j * self.ph_init)
self.u0 = fftn(self.u0)
self.u0 = shifted_fft(self.u0)
assert self.u0.shape == self.data.shape, 'Array size of the guess does not match the data'
else: #product space
......@@ -198,7 +209,7 @@ class Molecule3D(OrbitalTomographyExperiment):
ph_init_j = 2 * np.pi * np.random.random_sample(self.data.shape)
u0_j = abs(self.data) * np.exp(1j * ph_init_j)
self.u0[j] = u0_j
self.u0[j] = fftn(self.u0[j])
self.u0[j] = shifted_fft(self.u0[j])
else:
pass
def setupProxOperators(self):
......@@ -312,13 +323,14 @@ class Molecule3D(OrbitalTomographyExperiment):
:param name: sub-directory name
"""
# Create directory (if needed); does not raise an exception if the directory already exists
date_time = datetime.now() # current date and time
date_str = date_time.strftime("%Y_%m_%d")
# date_time = datetime.now() # current date and time
# date_str = date_time.strftime("%Y_%m_%d")
if savedir == '':
savedir = os.curdir
if name != '': # connect date and name using underscore
date_str += '_'
directory = savedir + date_str + name
# if name != '': # connect date and name using underscore
# date_str += '_'
# directory = savedir + date_str + name
directory = savedir + name
os.makedirs(directory, exist_ok=True)
file_id = 1
......@@ -377,7 +389,7 @@ class Molecule3D(OrbitalTomographyExperiment):
zz, zy, zx = tuple(s // 4 for s in self.output['u1'].shape)
else:
zz, zy, zx = zoom
u1 = roll_to_pos(self.output['u1'],pos=(nz//2,ny//2,nx//2)) #shift the mass of image to center
u1 = roll_to_pos(self.output['u1'],pos=(nz//2,ny//2,nx//2),move_maximum=True) #shift the mass of image to center
u2 = roll_to_pos(self.output['u2'],pos=(nz//2,ny//2,nx//2),move_maximum=True)
phys_constraint_plot = XYZStackViewer(u1[zz:nz - zz, zy:ny - zy, zx:nx - zx],
name='Physical constraint satisfied (%s)' % self.constraint)
......@@ -396,7 +408,7 @@ class Molecule3D(OrbitalTomographyExperiment):
shell = np.sqrt(k0**2-(yy-ny/2)**2-(xx-nx/2)**2)
return shell, yy, xx
def show_3d(self,inpt=True, fft=True):
def show_3d(self,inpt=True, fft=True, Real=True):
"""
Create results for 3d reconstruction
"""
......@@ -404,29 +416,36 @@ class Molecule3D(OrbitalTomographyExperiment):
u1 = self.data
else:
u1 = self.output['u1']
if self.constraint == 'sparse complex':
u1 = u1.real
nx0,ny0,nz0 = u1.shape
u1 = roll_to_pos(u1,pos=(nx0//2,ny0//2,nz0//2))
u1_rot = bin_array(u1,1)
nz,ny, nx = u1_rot.shape
zz, zy, zx = tuple(s // 4 for s in u1_rot.shape)
u1_rot=u1_rot[zz:nz - zz, zy:ny - zy, zx:nx - zx]
u1 = roll_to_pos(u1,pos=(nx0//2,ny0//2,nz0//2),move_maximum=True)#shift the mass of image to center
u1 = arrays.u1_process(u1,True)
if fft==True:
fft_u1 = shifted_fft(u1_rot)
fft_u1 = fourier_interpolate(fft_u1,4)
inp = np.transpose(fft_u1.imag,(2,1,0))
inp = u1
nx0,ny0,nz0 = inp.shape
inp = roll_to_pos(inp,pos=(nx0//2,ny0//2,nz0//2))
inp = bin_array(inp,2)
nz,ny, nx = inp.shape
zz, zy, zx = tuple(s // 4 for s in inp.shape)
inp=inp[zz:nz - zz, zy:ny - zy, zx:nx - zx]
fft_u1 = shifted_fft(inp)
fft_u1 = fourier_interpolate(fft_u1,2)
inp = fft_u1.real
edge = 113
src = mlab.pipeline.scalar_field(inp[:,:,5:])
src = mlab.pipeline.scalar_field(inp)
op = 1
mlab.close()
mlab.figure(bgcolor=(1,1,1),size=(1800,1200), fgcolor=(0,0,0))
# show iso surface molecule in Fourier transform
if op != 0:
mlab.pipeline.iso_surface(src, contours=[-0.07*inp.ptp(), ],
mlab.pipeline.iso_surface(src, contours=[-0.04*inp.ptp(), ],
opacity=op, color=(1,0,0))
mlab.pipeline.iso_surface(src, contours=[0.07*inp.ptp(), ],
mlab.pipeline.iso_surface(src, contours=[0.04*inp.ptp(), ],
opacity=op, color=(0,0,1))
# plot box and hemisphericals:
for k0 in [20,30,50]:
for k0 in [20,30,40]:
step = 0.1
s,y,x = self.shell(inp.shape,k0,step)
mlab.surf(y.T[:,:int(edge/step)],
......@@ -439,13 +458,13 @@ class Molecule3D(OrbitalTomographyExperiment):
mlab.show()
# show isosurface of molecule in real space
else:
x1,y1,z1 = u1_rot.shape
data = u1_rot
x1,y1,z1 = u1.shape
data = u1
src = mlab.pipeline.scalar_field(data)
mlab.close()
realspace = mlab.figure(bgcolor=(1,1,1),size=(1800,1200), fgcolor=(0,0,0))
mlab.pipeline.iso_surface(src, contours=[-0.07*data.ptp(), ],opacity=1,color=(1,0,0))
mlab.pipeline.iso_surface(src, contours=[0.07*data.ptp(), ],opacity=0.6,color=(0,0,1))
mlab.pipeline.iso_surface(src, contours=[0.07*data.ptp(), ],opacity=1,color=(0,0,1))
mlab.outline()
#mlab.axes(ranges=rngs)
mlab.show()
......
......@@ -91,7 +91,11 @@ class OrbitalTomographyExperiment(Experiment):
:return:
"""
self.output['plots'] = []
convergence_plots, ax = subplots(1, 2, figsize=(7, 3), num="Convergence plots")
convergence_plots, ax = subplots(1, 3, figsize=(10
, 3), num="Convergence plots")
if 'change' in self.output['stats']:
change_key = 'change'
else:
......@@ -99,13 +103,19 @@ class OrbitalTomographyExperiment(Experiment):
ax[0].semilogy(self.output['stats'][change_key])