phase.py

# -*- coding: utf-8 -*-

from proxtoolbox.Problems.problems import Problem
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
from numpy import square, sqrt


class Phase(Problem):
    """
    Phase Problem
    """
    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 with non-empty keys:
            object, constraint, experiment, distance, algorithm
            Ny, Nx
            noise
            MAXIT, TOL
            diagnostic
            data, data_norm, support,
        """
        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)

        if 'Nz' not in self.config:
            self.config['Nz'] = 1

        # 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
        used_proxoperators = ['', '', '']

        # Projector 1 (real / object space)
        if self.config['constraint'] == 'support only':
            used_proxoperators[0] = 'P_S'
        elif self.config['constraint'] == 'real and support':
            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'

        # 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 used_proxoperators:
            if prox != '':
                self.config['proxoperators'].append(getattr(ProxOperators, prox))

        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

        # 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':
            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:  # 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_data']) ** 2  # norm_rt_data
        else:
            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_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)

    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 problem

        Output can be accessed from self.output.
        """
        # 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)