RAAR.py

# -*- coding: utf-8 -*-
"""
Created on Mon Dec 14 12:48:26 2015

@author: rebecca
"""

from math import exp, sqrt
from numpy import zeros

from scipy.linalg import norm

from .algorithms import Algorithm

class RAAR(Algorithm):
    """
    Relaxed Averaged Alternating Reflection algorithm
    """
    
    def __init__(self,config):
        """
        Parameters
        ----------
        config : dict        
                 Dictionary containing the problem configuration.
                 It must contain the following mappings:
            
                proxoperators: 2 ProxOperators
                    Tuple of ProxOperators (the class, no instance)
                beta0: number
                    Starting relaxation parmater
                beta_max: number
                    Maximum relaxation parameter
                beta_switch: int
                    Iteration at which beta moves from beta0 -> beta_max
                normM: number
                    ?
                Nx: int
                    Row-dim of the product space elements
                Ny: int
                    Column-dim of the product space elements
                Nz: int
                    Depth-dim of the product space elements
                dim: int
                    Size of the product space
        """
        self.proj1 = config['proxoperators'][0](config)
        self.proj2 = config['proxoperators'][1](config)
        self.normM = config['normM']
        self.beta0 = config['beta0']
        self.beta_max = config['beta_max']
        self.beta_switch = config['beta_switch']
        self.Nx = config['Nx']; self.Ny = config['Ny']; self.Nz = config['Nz'];
        self.dim = config['dim']
        self.iters = 0

    def run(self, u, tol, maxiter):
        """
        Runs the algorithm for the specified input data
        """
        
        ##### PREPROCESSING
        normM = self.normM
        
        beta = self.beta0
        iters = self.iters
        change = zeros(maxiter+1,dtype=u.dtype)
        change[0] = 999
        gap = change.copy()
        
        tmp1 = 2*self.proj2.work(u) - u
        
        
        ##### LOOP
        while iters < maxiter and change[iters] >= tol:
            tmp = exp((-iters/self.beta_switch)**3);
            beta = (tmp*self.beta0) + ((1-tmp)*self.beta_max);
            iters += 1;
            
            tmp3 = self.proj1.work(tmp1);
            tmp_u = ((beta*(2*tmp3-tmp1)) + ((1-beta)*tmp1) + u)/2;
            tmp2 = self.proj2.work(tmp_u);
            
            tmp3 = self.proj1.work(tmp2);
            
            tmp_change = 0; tmp_gap = 0;
            if self.Ny == 1 or self.Nx == 1:
                tmp_change = (norm(u-tmp_u,'fro')/normM)**2;
                tmp_gap = (norm(tmp3-tmp2,'fro')/normM)**2;
            elif self.Nz == 1:
                for j in range(self.dim):
                    tmp_change += (norm(u[:,:,j]-tmp_u[:,:,j],'fro')/normM)**2;
                    tmp_gap += (norm(tmp3[:,:,j]-tmp2[:,:,j],'fro')/normM)**2;
            else:
                for j in range(self.dim):
                    for k in range(self.Nz):
                        tmp_change += (norm(u[:,:,k,j]-tmp_u[:,:,k,j],'fro')/normM)**2;
                        tmp_gap += (norm(tmp3[:,:,k,j]-tmp2[:,:,k,j],'fro')/normM)**2;
            
            change[iters] = sqrt(tmp_change);
            gap[iters] = sqrt(tmp_gap);
            
            u = tmp_u;
            tmp1 = (2*tmp2) - tmp_u;
            
        
        ##### POSTPROCESSING
        u = tmp2;
        tmp = self.proj1.work(u);
        tmp2 = self.proj2.work(u);
        if self.Ny == 1:
            u1 = tmp[:,1];
            u2 = tmp2[:,1];
        elif self.Nx == 1:
            u1 = tmp[1,:];
            u2 = tmp2[1,:];
        elif self.Nz == 1:
            u1 = tmp[:,:,1];
            u2 = tmp2[:,:,1];
        else:
            u1 = tmp;
            u2 = tmp2;
        change = change[1:iters+1];
        gap = gap[1:iters+1];
        
        return u1, u2, iters, change, gap