Commit 09f1a9ba authored by alexander.dornheim's avatar alexander.dornheim
Browse files

Added first files for JWST Phase, not yet working.

parent 17239173
# -*- coding: utf-8 -*-
new_config = {
'''We start very general.
What type of problem is being solved? Classification
is according to the geometry: Affine, Cone, Convex,
Phase, Affine-sparsity, Nonlinear-sparsity, Sudoku'''
'problem_family' : 'Phase',
'''==========================================
Problem parameters
==========================================
What is the name of the data file?'''
'data_filename' : 'JWST_data_processor',
''' What type of object are we working with?
Options are: 'phase', 'real', 'nonnegative', 'complex' '''
'object' : 'complex',
''' What type of constraints do we have?
Options are: 'support only', 'real and support', 'nonnegative and support',
'amplitude only', 'sparse real', 'sparse complex', and 'hybrid' '''
'constraint' : 'amplitude only',
''' What type of measurements are we working with?
Options are: 'single diffraction', 'diversity diffraction',
'ptychography', and 'complex' '''
'experiment' : 'JWST',
''' Next we move to things that most of our users will know
better than we will. Some of these may be overwritten in the
data processor file which the user will most likely write.
Are the measurements in the far field or near field?
Options are: 'far field' or 'near field' '''
'distance' : 'far field',
# What are the dimensions of the measurements?
'Nx' : 128,
'Ny' : 128,
'Nz' : 1, # do not formulate in the product space
#if 'distance' =='near field':
# 'fresnel_nr' : 1*2*pi*config['Nx'],
# 'use_farfield_formula' : 0,
#else:
'fresnel_nr' : 0, #1*2*pi*prbl.Nx;
'use_farfield_formula' : 1,
# What are the noise characteristics (Poisson or Gaussian or none)?
'noise' : 'Poisson', #'Poisson',
'''==========================================
Algorithm parameters
==========================================
Now set some algorithm parameters that the user should be
able to control (without too much damage)'''
# Algorithm:
'method' : 'RAAR', #'Accelerated_AP_product_space';
'numruns' : 1, # the only time this parameter will
# be different than 1 is when we are
# benchmarking...not something a normal user
# would be doing.
# The following are parameters specific to RAAR, HPR, and HAAR that the
# user should be able to set/modify. Surely
# there will be other algorithm specific parameters that a user might
# want to play with. Don't know how best
# to do this. Thinking of a GUI interface, we could hard code all the
# parameters the user might encounter and have the menu options change
# depending on the value of the prbl.method field.
# do different things depending on the chosen algorithm:
#maximum number of iterations and tolerances
'MAXIT' : 500,
'TOL' : 1e-8,
# relaxaton parameters in RAAR, HPR and HAAR
'beta_0' : 1.0, # starting relaxation prameter (only used with
# HAAR, HPR and RAAR)
'beta_max' : 1.0, # maximum relaxation prameter (only used with
# HAAR, RAAR, and HPR)
'beta_switch' : 20, # iteration at which beta moves from beta_0 -> beta_max
# parameter for the data regularization
# need to discuss how/whether the user should
# put in information about the noise
'data_ball' : 1e-15,
# the above is the percentage of the gap
# between the measured data and the
# initial guess satisfying the
# qualitative constraints. For a number
# very close to one, the gap is not expected
# to improve much. For a number closer to 0
# the gap is expected to improve a lot.
# Ultimately the size of the gap depends
# on the inconsistency of the measurement model
# with the qualitative constraints.
#==========================================
# parameters for plotting and diagnostics
#==========================================
'verbose' : 1, # options are 0 or 1
'graphics' : 1, # whether or not to display figures, options are 0 or 1.
# default is 1.
'anim' : 2, # whether or not to disaply ``real time" reconstructions
# options are 0=no, 1=yes, 2=make a movie
# default is 1.
'graphics_display' : 'JWST_graphics' # unless specified, a default
# plotting subroutine will generate
# the graphics. Otherwise, the user
# can write their own plotting subroutine
#======================================================================
# Technical/software specific parameters
#======================================================================
# Given the parameter values above, the following technical/algorithmic
# parameters are automatically set. The user does not need to know
# about these details, and so probably these parameters should be set in
# a module one level below this one.
}
'''
JWST_data_processor.m
written on May 27, 2012 by
Russell Luke
Inst. Fuer Numerische und Angewandte Mathematik
Universitaet Goettingen
DESCRIPTION:
INPUT: input = a data structure
OUTPUT: input = a data structure
USAGE: input = JWST_data_processor(input)
Data loaded:
Nonstandard function calls: Rbin, ZeroPad, Mgen, Resize, OnesPad
data reader/processor
'''
import numpy as np
from numpy import fromfile, exp
def JWST_data_processor(config):
data_ball = config['data_ball'];
newres = config['Nx'];
noise = config['noise'];
snr = config['data_ball'];
with open('pupil.pmod','r') as fid:
# read lower endian float <f
Xi_A = np.fromfile(fid, dtype='<f');
Xi_A.reshape(512,512);
print(Xi_A);
diversity=3;
with open('phase_p37.pmod','r') as fid:
# read lower endian float <f
temp1 = np.fromfile(fid, dtype='<f');
temp1.reshape(512,512);
with open('phase_m37.pmod','r') as fid:
# read lower endian float <f
temp2 = np.fromfile(fid, dtype='<f');
temp2.reshape(512,512);
defocus=(temp1-temp2)/2;
theta=(temp1+temp2)/2;
if newres != 512:
defocus.resize((newres,newres),refcheck =False);
theta.resize((newres,newres),refcheck =False);
Xi_A.resize((newres,newres),refcheck =False);
aberration = np.array( [np.zeros((newres,newres)), defocus, -defocus]);
print(aberration.shape);
aberration = aberration.transpose((1,2,0));
print(aberration.shape);
# aberration[:,:,1]=np.zeros((newres,newres)); # Note this order!!!!!
# aberration[:,:,2]=defocus; # Note this order!!!!!
# aberration[:,:,3]=-defocus; # Note this order!!!!!
true_object = np.multiply(Xi_A,exp(1j*(theta)));
k= np.zeros((newres,newres,diversity));
rt_k= np.zeros((newres,newres,diversity));
for j in range(diversity):
Pj = np.multiply(Xi_A,exp(1j*(aberration[:,:,j]+theta)));
print(Pj.shape);
k[:,:,j] = np.square(abs(np.fft.fft(Pj)));
#clear Pj temp1 temp2 defocus Automatic garbage collection?
epsilon=0;
norm_rt_data=np.sqrt(sum(sum(Xi_A)));
data_zeros= np.zeros(newres^2,3);
for i in range(diversity):
rt_k[:,:,i]= np.sqrt(k[:,:,i]); # *newres ?
# normalize the data so that the corresponding
# pupil function has amplitude 1 everywhere
rt_k[:,:,] = rt_k[:,:,i]/norm(IFFT(rt_k[:,:,i]),'fro')*norm_rt_data;
k[:,:,i] = np.square(rt_k[:,:,i]);
if config['noise'] == 'Poisson':
# Add Poisson noise according to Ning Lei:
# f = alpha*4*4*pi/q/q*abs(cos(qx(iX))*sin(qy(iY))*(sin(q)/q-cos(q)));
# f2 = PoissonRan(f*f*2)/alpha/alpha/2; # 2 is for I(-q)=I(q)
# sigma(iX, iY, iZ) = f/np.sqrt(2)/alpha/alpha/abs(ff)/abs(ff);
# July 15, 2010: Since this is meant to model photon counts (which
# are integers) we add noise and then round down
k[:,:,i] = k[:,:,i]/snr;
for ii in range(newres):
for jj in range(newres):
k[ii,jj,i]= PoissonRan(k[ii,jj,i])*snr;
k[:,:,i]=round(k[:,:,i]);
rt_k[:,:,i]= np.sqrt(k[:,:,i]);
# fftshift and scale the data for use with fft:
rt_k[:,:,i] = fftshift(rt_k[:,:,i])*newres;
k[:,:,i] = fftshift(k[:,:,i])*newres;
temp = find(rt_k[:,:,i]==0);
data_zeros[range(length(temp)),i] = temp;
config['rt_data'] = rt_k;
config['data'] = k;
config['norm_rt_data'] =norm_rt_data; # this is correct since
# is is calculated in the
# object domain
config['data_zeros'] = data_zeros;
config['supp_ampl'] = find(Xi_A != 0);
config['indicator_ampl'] = zeros(size(Xi_A));
config['indicator_ampl'][config['supp_ampl']] = 1;
config['product_space_dimension'] = diversity+1;
# config.Nz = 1; # already set in the config file, but just in case
config['abs_illumination'] = Xi_A;
# Normalize the illumination so that it has the
# same energy as the data.
config['abs_illumination'] = config['abs_illumination']/norm(config['abs_illumination'],'fro')*config['norm_rt_data'][0];
config['illumination_phase'] = aberration;
# initial guess
# config.u_0 = Xi_A.*exp(1i*2*pi*rand(newres));
for j in range(config['product_space_dimension']):
config['u_0'][:,:,j] = np.multiply(config['abs_illumination'],exp(1j*2*pi*rand(newres)));
config['truth'] = true_object;
config['truth_dim'] = size(true_object);
config['norm_truth'] = norm(true_object,'fro');
#clear I Y k rt_k Xi_A temp aberration theta
# -*- coding: utf-8 -*-
new_config = {
'''We start very general.
What type of problem is being solved? Classification
is according to the geometry: Affine, Cone, Convex,
Phase, Affine-sparsity, Nonlinear-sparsity, Sudoku'''
'problem_family' : 'Phase',
'''==========================================
Problem parameters
==========================================
What is the name of the data file?'''
'data_filename' : 'JWST_data_processor',
''' What type of object are we working with?
Options are: 'phase', 'real', 'nonnegative', 'complex' '''
'object' : 'complex',
''' What type of constraints do we have?
Options are: 'support only', 'real and support', 'nonnegative and support',
'amplitude only', 'sparse real', 'sparse complex', and 'hybrid' '''
'constraint' : 'amplitude only',
''' What type of measurements are we working with?
Options are: 'single diffraction', 'diversity diffraction',
'ptychography', and 'complex' '''
'experiment' : 'JWST',
''' Next we move to things that most of our users will know
better than we will. Some of these may be overwritten in the
data processor file which the user will most likely write.
Are the measurements in the far field or near field?
Options are: 'far field' or 'near field' '''
'distance' : 'far field',
# What are the dimensions of the measurements?
'Nx' : 128,
'Ny' : 128,
'Nz' : 1, # do not formulate in the product space
#if 'distance' =='near field':
# 'fresnel_nr' : 1*2*pi*config['Nx'],
# 'use_farfield_formula' : 0,
#else:
'fresnel_nr' : 0, #1*2*pi*prbl.Nx;
'use_farfield_formula' : 1,
# What are the noise characteristics (Poisson or Gaussian or none)?
'noise' : 'Poisson', #'Poisson',
'''==========================================
Algorithm parameters
==========================================
Now set some algorithm parameters that the user should be
able to control (without too much damage)'''
# Algorithm:
'method' : 'RAAR', #'Accelerated_AP_product_space';
'numruns' : 1, # the only time this parameter will
# be different than 1 is when we are
# benchmarking...not something a normal user
# would be doing.
# The following are parameters specific to RAAR, HPR, and HAAR that the
# user should be able to set/modify. Surely
# there will be other algorithm specific parameters that a user might
# want to play with. Don't know how best
# to do this. Thinking of a GUI interface, we could hard code all the
# parameters the user might encounter and have the menu options change
# depending on the value of the prbl.method field.
# do different things depending on the chosen algorithm:
#maximum number of iterations and tolerances
'MAXIT' : 500,
'TOL' : 1e-8,
# relaxaton parameters in RAAR, HPR and HAAR
'beta_0' : 1.0, # starting relaxation prameter (only used with
# HAAR, HPR and RAAR)
'beta_max' : 1.0, # maximum relaxation prameter (only used with
# HAAR, RAAR, and HPR)
'beta_switch' : 20, # iteration at which beta moves from beta_0 -> beta_max
# parameter for the data regularization
# need to discuss how/whether the user should
# put in information about the noise
'data_ball' : 1e-15,
# the above is the percentage of the gap
# between the measured data and the
# initial guess satisfying the
# qualitative constraints. For a number
# very close to one, the gap is not expected
# to improve much. For a number closer to 0
# the gap is expected to improve a lot.
# Ultimately the size of the gap depends
# on the inconsistency of the measurement model
# with the qualitative constraints.
#==========================================
# parameters for plotting and diagnostics
#==========================================
'verbose' : 1, # options are 0 or 1
'graphics' : 1, # whether or not to display figures, options are 0 or 1.
# default is 1.
'anim' : 2, # whether or not to disaply ``real time" reconstructions
# options are 0=no, 1=yes, 2=make a movie
# default is 1.
'graphics_display' : 'JWST_graphics' # unless specified, a default
# plotting subroutine will generate
# the graphics. Otherwise, the user
# can write their own plotting subroutine
#======================================================================
# Technical/software specific parameters
#======================================================================
# Given the parameter values above, the following technical/algorithmic
# parameters are automatically set. The user does not need to know
# about these details, and so probably these parameters should be set in
# a module one level below this one.
}
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 16 10:15:22 2014
@author: stefan
"""
from numpy import amax, argmax, array, float32, int32, zeros, zeros_like
from matplotlib import pyplot
from .problems import Problem
from proxtoolbox import Algorithms
from proxtoolbox import ProxOperators
from proxtoolbox.ProxOperators.proxoperators import ProxOperator
#__all__ = ["Sudoku"]
class ProjGiven(ProxOperator):
"""
Projection onto the given entries in a Sudoku problem
"""
def __init__(self,config):
"""
Initialization method
Parameters
----------
config : dict - Dictionary containing the problem configuration. It must have the key 'sudoku' mapping to the input Sudoku.
"""
self.given = config['sudoku'].copy()
def work(self,u):
"""
Applies the prox operator to the input data
Parameters
----------
u : array_like - Input data for the operator
Returns
-------
array_like - Result of the operation
"""
A = zeros((9,9,9),dtype=u.dtype)
for x in range(9):
for y in range(9):
z = self.given[x,y]
if z > 0:
A[x,y,z-1] = 1
else:
A[x,y,argmax(u[x,y,:])] = 1
return A
class ProjSquare(ProxOperator):
"""
Projection onto the box constraints in a Sudoku problem
"""
def __init__(self,config=None):
"""
Initialization
Parameters
-----------
config : dict, optional - Not used here
"""
return
def work(self,u):
"""
Applies the prox operator to the input data
Parameters
----------
u : array_like - Input data for the operator
Returns
-------
array_like - Result of the operation
"""
Q = zeros((9,9,9),dtype=u.dtype)
for z in range(9):
for x in range(0,9,3):
for y in range(0,9,3):
v = argmax(u[x:(x+3),y:(y+3),z],axis=0)
w = argmax(amax(u[x:(x+3),y:(y+3),z],axis=0))
Q[x+v[w],y+w,z] = 1
return Q
class ProjColumn(ProxOperator):
"""
Projection onto the column constraints in a Sudoku problem
"""
def __init__(self,config):
"""
Initialization
Parameters
----------
config : dict, optional - Not used here
"""
return
def work(self,u):
"""
Applies the prox operator to the input data
Parameters
----------
u : array_like - Input data for the operator
Returns
-------
array_like - Result of the operation
"""
C = zeros((9,9,9),dtype=u.dtype)
for x in range(9):
for z in range(9):
y = argmax(u[x,:,z])
C[x,y,z] = 1
return C
class ProjRow(ProxOperator):
"""
Projection onto the row constraints in a Sudoku problem
"""
def __init__(self,config):
"""
Initialization
Parameters
----------
config : dict, optional - Not used here
"""
return
def work(self,u):
"""
Applies the prox operator to the input data
Parameters
----------
u : array_like - Input data for the operator
Returns
-------
array_like - Result of the operation
"""
R = zeros((9,9,9),dtype=u.dtype)
for y in range(9):
for z in range(9):
x = argmax(u[:,y,z])
R[x,y,z] = 1
return R
class Sudoku(Problem):
"""
Sudoku Problem
The goal of a standard Sudoku puzzle is to fill a 9x9 array with the numbers
from 1 to 9. Every row, every column and each of the 9 3x3 subblocks should
contain each number only once.
Starting point is a grid that already contains several numbers.
Usually, there exists a unique solution.
"""
config = {
# This is the algorithm we use. RAAR and HPR will work.
'algorithm':'RAAR',
# RAAR requires 2 ProxOperators
'proxoperators':('P_diag','P_parallel'),
# P_parallel requires a sequence of projectors
'projectors':('ProjRow','ProjColumn','ProjSquare','ProjGiven'),
# Relaxation parameters for RAAR/HPR
'beta0':1,
'beta_max':1,
'beta_switch':1,
# Any algorithm requires these
'maxiter':2000,
'tol':1e-9,