nikcleju@21: # -*- coding: utf-8 -*-
nikcleju@21: """
nikcleju@21: Created on Thu Oct 13 14:05:22 2011
nikcleju@21: 
nikcleju@21: @author: ncleju
nikcleju@21: """
nikcleju@21: 
nikcleju@21: #from numpy import *
nikcleju@21: #from scipy import *
nikcleju@27: import numpy
nikcleju@27: import numpy.linalg
nikcleju@21: import scipy as sp
nikcleju@27: #import scipy.stats #from scipy import stats
nikcleju@27: #import scipy.linalg #from scipy import lnalg
nikcleju@21: import math
nikcleju@21: 
nikcleju@21: from numpy.random import RandomState
nikcleju@21: rng = RandomState()
nikcleju@21: 
nikcleju@21: def Generate_Analysis_Operator(d, p):
nikcleju@21:   # generate random tight frame with equal column norms
nikcleju@21:   if p == d:
nikcleju@21:     T = rng.randn(d,d);
nikcleju@27:     [Omega, discard] = numpy.qr(T);
nikcleju@21:   else:
nikcleju@21:     Omega = rng.randn(p, d);
nikcleju@27:     T = numpy.zeros((p, d));
nikcleju@21:     tol = 1e-8;
nikcleju@21:     max_j = 200;
nikcleju@21:     j = 1;
nikcleju@27:     while (sum(sum(abs(T-Omega))) > numpy.dot(tol,numpy.dot(p,d)) and j < max_j):
nikcleju@21:         j = j + 1;
nikcleju@21:         T = Omega;
nikcleju@27:         [U, S, Vh] = numpy.linalg.svd(Omega);
nikcleju@21:         V = Vh.T
nikcleju@21:         #Omega = U * [eye(d); zeros(p-d,d)] * V';
nikcleju@27:         Omega2 = numpy.dot(numpy.dot(U, numpy.concatenate((numpy.eye(d), numpy.zeros((p-d,d))))), V.transpose())
nikcleju@21:         #Omega = diag(1./sqrt(diag(Omega*Omega')))*Omega;
nikcleju@27:         Omega = numpy.dot(numpy.diag(1.0 / numpy.sqrt(numpy.diag(numpy.dot(Omega2,Omega2.transpose())))), Omega2)
nikcleju@21:     #end
nikcleju@21:     ##disp(j);
nikcleju@21:     #end
nikcleju@21:   return Omega
nikcleju@21: 
nikcleju@21: 
nikcleju@21: def Generate_Data_Known_Omega(Omega, d,p,m,k,noiselevel, numvectors, normstr):
nikcleju@21:   #function [x0,y,M,LambdaMat] = Generate_Data_Known_Omega(Omega, d,p,m,k,noiselevel, numvectors, normstr)
nikcleju@21:   
nikcleju@21:   # Building an analysis problem, which includes the ingredients: 
nikcleju@21:   #   - Omega - the analysis operator of size p*d
nikcleju@21:   #   - M is anunderdetermined measurement matrix of size m*d (m<d)
nikcleju@21:   #   - x0 is a vector of length d that satisfies ||Omega*x0||=p-k
nikcleju@21:   #   - Lambda is the true location of these k zeros in Omega*x0
nikcleju@21:   #   - a measurement vector y0=Mx0 is computed
nikcleju@21:   #   - noise contaminated measurement vector y is obtained by
nikcleju@21:   #     y = y0 + n where n is an additive gaussian noise with norm(n,2)/norm(y0,2) = noiselevel
nikcleju@21:   # Added by Nic:
nikcleju@21:   #   - Omega = analysis operator
nikcleju@21:   #   - normstr: if 'l0', generate l0 sparse vector (unchanged). If 'l1',
nikcleju@21:   #   generate a vector of Laplacian random variables (gamma) and
nikcleju@21:   #   pseudoinvert to find x
nikcleju@21: 
nikcleju@21:   # Omega is known as input parameter
nikcleju@21:   #Omega=Generate_Analysis_Operator(d, p);
nikcleju@21:   # Omega = randn(p,d);
nikcleju@21:   # for i = 1:size(Omega,1)
nikcleju@21:   #     Omega(i,:) = Omega(i,:) / norm(Omega(i,:));
nikcleju@21:   # end
nikcleju@21:   
nikcleju@21:   #Init
nikcleju@27:   LambdaMat = numpy.zeros((k,numvectors))
nikcleju@27:   x0 = numpy.zeros((d,numvectors))
nikcleju@27:   y = numpy.zeros((m,numvectors))
nikcleju@27:   realnoise = numpy.zeros((m,numvectors))
nikcleju@21:   
nikcleju@21:   M = rng.randn(m,d);
nikcleju@21:   
nikcleju@21:   #for i=1:numvectors
nikcleju@21:   for i in range(0,numvectors):
nikcleju@21:     
nikcleju@21:     # Generate signals
nikcleju@21:     #if strcmp(normstr,'l0')
nikcleju@21:     if normstr == 'l0':
nikcleju@21:         # Unchanged
nikcleju@21:         
nikcleju@21:         #Lambda=randperm(p); 
nikcleju@21:         Lambda = rng.permutation(int(p));
nikcleju@27:         Lambda = numpy.sort(Lambda[0:k]); 
nikcleju@21:         LambdaMat[:,i] = Lambda; # store for output
nikcleju@21:         
nikcleju@21:         # The signal is drawn at random from the null-space defined by the rows 
nikcleju@21:         # of the matreix Omega(Lambda,:)
nikcleju@27:         [U,D,Vh] = numpy.linalg.svd(Omega[Lambda,:]);
nikcleju@21:         V = Vh.T
nikcleju@21:         NullSpace = V[:,k:];
nikcleju@27:         #print numpy.dot(NullSpace, rng.randn(d-k,1)).shape
nikcleju@21:         #print x0[:,i].shape
nikcleju@27:         x0[:,i] = numpy.squeeze(numpy.dot(NullSpace, rng.randn(d-k,1)));
nikcleju@21:         # Nic: add orthogonality noise
nikcleju@21:         #     orthonoiseSNRdb = 6;
nikcleju@21:         #     n = randn(p,1);
nikcleju@21:         #     #x0(:,i) = x0(:,i) + n / norm(n)^2 * norm(x0(:,i))^2 / 10^(orthonoiseSNRdb/10);
nikcleju@21:         #     n = n / norm(n)^2 * norm(Omega * x0(:,i))^2 / 10^(orthonoiseSNRdb/10);
nikcleju@21:         #     x0(:,i) = pinv(Omega) * (Omega * x0(:,i) + n);
nikcleju@21:         
nikcleju@21:     #elseif strcmp(normstr, 'l1')
nikcleju@21:     elif normstr == 'l1':
nikcleju@21:         print('Nic says: not implemented yet')
nikcleju@21:         raise Exception('Nic says: not implemented yet')
nikcleju@21:         #gamma = laprnd(p,1,0,1);
nikcleju@21:         #x0(:,i) = Omega \ gamma;
nikcleju@21:     else:
nikcleju@21:         #error('normstr must be l0 or l1!');
nikcleju@21:         print('Nic says: not implemented yet')
nikcleju@21:         raise Exception('Nic says: not implemented yet')
nikcleju@21:     #end
nikcleju@21:     
nikcleju@21:     # Acquire measurements
nikcleju@27:     y[:,i]  = numpy.dot(M, x0[:,i])
nikcleju@21: 
nikcleju@21:     # Add noise
nikcleju@27:     t_norm = numpy.linalg.norm(y[:,i],2);
nikcleju@27:     n = numpy.squeeze(rng.randn(m, 1));
nikcleju@27:     y[:,i] = y[:,i] + noiselevel * t_norm * n / numpy.linalg.norm(n, 2);
nikcleju@27:     realnoise[:,i] = noiselevel * t_norm * n / numpy.linalg.norm(n, 2)
nikcleju@21:   #end
nikcleju@21: 
nikcleju@21:   return x0,y,M,LambdaMat,realnoise
nikcleju@21: 
nikcleju@21: #####################
nikcleju@21: 
nikcleju@21: #function [xhat, arepr, lagmult] = ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, ilagmult, params)
nikcleju@21: def ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, ilagmult, params):
nikcleju@21:   
nikcleju@21:     #
nikcleju@21:     # This function aims to compute
nikcleju@21:     #    xhat = argmin || Omega(Lambdahat, :) * x ||_2   subject to  || y - M*x ||_2 <= epsilon.
nikcleju@21:     # arepr is the analysis representation corresponding to Lambdahat, i.e.,
nikcleju@21:     #    arepr = Omega(Lambdahat, :) * xhat.
nikcleju@21:     # The function also returns the lagrange multiplier in the process used to compute xhat.
nikcleju@21:     #
nikcleju@21:     # Inputs:
nikcleju@21:     #    y : observation/measurements of an unknown vector x0. It is equal to M*x0 + noise.
nikcleju@21:     #    M : Measurement matrix
nikcleju@21:     #    MH : M', the conjugate transpose of M
nikcleju@21:     #    Omega : analysis operator
nikcleju@21:         #    OmegaH : Omega', the conjugate transpose of Omega. Also, synthesis operator.
nikcleju@21:     #    Lambdahat : an index set indicating some rows of Omega.
nikcleju@21:     #    xinit : initial estimate that will be used for the conjugate gradient algorithm.
nikcleju@21:     #    ilagmult : initial lagrange multiplier to be used in
nikcleju@21:     #    params : parameters
nikcleju@21:     #        params.noise_level : this corresponds to epsilon above.
nikcleju@21:     #        params.max_inner_iteration : `maximum' number of iterations in conjugate gradient method.
nikcleju@21:     #        params.l2_accurary : the l2 accuracy parameter used in conjugate gradient method
nikcleju@21:     #        params.l2solver : if the value is 'pseudoinverse', then direct matrix computation (not conjugate gradient method) is used. Otherwise, conjugate gradient method is used.
nikcleju@21:     #
nikcleju@21:     
nikcleju@21:     #d = length(xinit)
nikcleju@21:     d = xinit.size
nikcleju@21:     lagmultmax = 1e5;
nikcleju@21:     lagmultmin = 1e-4;
nikcleju@21:     lagmultfactor = 2.0;
nikcleju@21:     accuracy_adjustment_exponent = 4/5.;
nikcleju@21:     lagmult = max(min(ilagmult, lagmultmax), lagmultmin);
nikcleju@21:     was_infeasible = 0;
nikcleju@21:     was_feasible = 0;
nikcleju@21:     
nikcleju@21:     #######################################################################
nikcleju@21:     ## Computation done using direct matrix computation from matlab. (no conjugate gradient method.)
nikcleju@21:     #######################################################################
nikcleju@21:     #if strcmp(params.l2solver, 'pseudoinverse')
nikcleju@21:     if params['l2solver'] == 'pseudoinverse':
nikcleju@21:     #if strcmp(class(M), 'double') && strcmp(class(Omega), 'double')
nikcleju@21:       if M.dtype == 'float64' and Omega.dtype == 'double':
nikcleju@21:         while True:
nikcleju@21:             alpha = math.sqrt(lagmult);
nikcleju@27:             xhat = numpy.linalg.lstsq(numpy.concatenate((M, alpha*Omega[Lambdahat,:])), numpy.concatenate((y, numpy.zeros(Lambdahat.size))))[0]
nikcleju@27:             temp = numpy.linalg.norm(y - numpy.dot(M,xhat), 2);
nikcleju@21:             #disp(['fidelity error=', num2str(temp), ' lagmult=', num2str(lagmult)]);
nikcleju@21:             if temp <= params['noise_level']:
nikcleju@21:                 was_feasible = True;
nikcleju@21:                 if was_infeasible:
nikcleju@21:                     break;
nikcleju@21:                 else:
nikcleju@21:                     lagmult = lagmult*lagmultfactor;
nikcleju@21:             elif temp > params['noise_level']:
nikcleju@21:                 was_infeasible = True;
nikcleju@21:                 if was_feasible:
nikcleju@21:                     xhat = xprev.copy();
nikcleju@21:                     break;
nikcleju@21:                 lagmult = lagmult/lagmultfactor;
nikcleju@21:             if lagmult < lagmultmin or lagmult > lagmultmax:
nikcleju@21:                 break;
nikcleju@21:             xprev = xhat.copy();
nikcleju@27:         arepr = numpy.dot(Omega[Lambdahat, :], xhat);
nikcleju@21:         return xhat,arepr,lagmult;
nikcleju@21: 
nikcleju@21: 
nikcleju@21:     ########################################################################
nikcleju@21:     ## Computation using conjugate gradient method.
nikcleju@21:     ########################################################################
nikcleju@21:     #if strcmp(class(MH),'function_handle') 
nikcleju@21:     if hasattr(MH, '__call__'):
nikcleju@21:         b = MH(y);
nikcleju@21:     else:
nikcleju@27:         b = numpy.dot(MH, y);
nikcleju@21:     
nikcleju@27:     norm_b = numpy.linalg.norm(b, 2);
nikcleju@21:     xhat = xinit.copy();
nikcleju@21:     xprev = xinit.copy();
nikcleju@21:     residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b;
nikcleju@21:     direction = -residual;
nikcleju@21:     iter = 0;
nikcleju@21:     
nikcleju@21:     while iter < params.max_inner_iteration:
nikcleju@21:         iter = iter + 1;
nikcleju@27:         alpha = numpy.linalg.norm(residual,2)**2 / numpy.dot(direction.T, TheHermitianMatrix(direction, M, MH, Omega, OmegaH, Lambdahat, lagmult));
nikcleju@21:         xhat = xhat + alpha*direction;
nikcleju@21:         prev_residual = residual.copy();
nikcleju@21:         residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b;
nikcleju@27:         beta = numpy.linalg.norm(residual,2)**2 / numpy.linalg.norm(prev_residual,2)**2;
nikcleju@21:         direction = -residual + beta*direction;
nikcleju@21:         
nikcleju@27:         if numpy.linalg.norm(residual,2)/norm_b < params['l2_accuracy']*(lagmult**(accuracy_adjustment_exponent)) or iter == params['max_inner_iteration']:
nikcleju@21:             #if strcmp(class(M), 'function_handle')
nikcleju@21:             if hasattr(M, '__call__'):
nikcleju@27:                 temp = numpy.linalg.norm(y-M(xhat), 2);
nikcleju@21:             else:
nikcleju@27:                 temp = numpy.linalg.norm(y-numpy.dot(M,xhat), 2);
nikcleju@21:             
nikcleju@21:             #if strcmp(class(Omega), 'function_handle')
nikcleju@21:             if hasattr(Omega, '__call__'):
nikcleju@21:                 u = Omega(xhat);
nikcleju@27:                 u = math.sqrt(lagmult)*numpy.linalg.norm(u(Lambdahat), 2);
nikcleju@21:             else:
nikcleju@27:                 u = math.sqrt(lagmult)*numpy.linalg.norm(Omega[Lambdahat,:]*xhat, 2);
nikcleju@21:             
nikcleju@21:             
nikcleju@21:             #disp(['residual=', num2str(norm(residual,2)), ' norm_b=', num2str(norm_b), ' omegapart=', num2str(u), ' fidelity error=', num2str(temp), ' lagmult=', num2str(lagmult), ' iter=', num2str(iter)]);
nikcleju@21:             
nikcleju@21:             if temp <= params['noise_level']:
nikcleju@21:                 was_feasible = True;
nikcleju@21:                 if was_infeasible:
nikcleju@21:                     break;
nikcleju@21:                 else:
nikcleju@21:                     lagmult = lagmultfactor*lagmult;
nikcleju@21:                     residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b;
nikcleju@21:                     direction = -residual;
nikcleju@21:                     iter = 0;
nikcleju@21:             elif temp > params['noise_level']:
nikcleju@21:                 lagmult = lagmult/lagmultfactor;
nikcleju@21:                 if was_feasible:
nikcleju@21:                     xhat = xprev.copy();
nikcleju@21:                     break;
nikcleju@21:                 was_infeasible = True;
nikcleju@21:                 residual = TheHermitianMatrix(xhat, M, MH, Omega, OmegaH, Lambdahat, lagmult) - b;
nikcleju@21:                 direction = -residual;
nikcleju@21:                 iter = 0;
nikcleju@21:             if lagmult > lagmultmax or lagmult < lagmultmin:
nikcleju@21:                 break;
nikcleju@21:             xprev = xhat.copy();
nikcleju@21:         #elseif norm(xprev-xhat)/norm(xhat) < 1e-2
nikcleju@21:         #    disp(['rel_change=', num2str(norm(xprev-xhat)/norm(xhat))]);
nikcleju@21:         #    if strcmp(class(M), 'function_handle')
nikcleju@21:         #        temp = norm(y-M(xhat), 2);
nikcleju@21:         #    else
nikcleju@21:         #        temp = norm(y-M*xhat, 2);
nikcleju@21:         #    end
nikcleju@21:     #
nikcleju@21:     #        if temp > 1.2*params.noise_level
nikcleju@21:     #            was_infeasible = 1;
nikcleju@21:     #            lagmult = lagmult/lagmultfactor;
nikcleju@21:     #            xprev = xhat;
nikcleju@21:     #        end
nikcleju@21:     
nikcleju@21:     #disp(['fidelity_error=', num2str(temp)]);
nikcleju@21:     print 'fidelity_error=',temp
nikcleju@21:     #if iter == params['max_inner_iteration']:
nikcleju@21:         #disp('max_inner_iteration reached. l2_accuracy not achieved.');
nikcleju@21:     
nikcleju@21:     ##
nikcleju@21:     # Compute analysis representation for xhat
nikcleju@21:     ##
nikcleju@21:     #if strcmp(class(Omega),'function_handle') 
nikcleju@21:     if hasattr(Omega, '__call__'):
nikcleju@21:         temp = Omega(xhat);
nikcleju@21:         arepr = temp(Lambdahat);
nikcleju@21:     else:    ## here Omega is assumed to be a matrix
nikcleju@27:         arepr = numpy.dot(Omega[Lambdahat, :], xhat);
nikcleju@21:     
nikcleju@21:     return xhat,arepr,lagmult
nikcleju@21: 
nikcleju@21: 
nikcleju@21: ##
nikcleju@21: # This function computes (M'*M + lm*Omega(L,:)'*Omega(L,:)) * x.
nikcleju@21: ##
nikcleju@21: #function w = TheHermitianMatrix(x, M, MH, Omega, OmegaH, L, lm)
nikcleju@21: def TheHermitianMatrix(x, M, MH, Omega, OmegaH, L, lm):
nikcleju@21:     #if strcmp(class(M), 'function_handle')
nikcleju@21:     if hasattr(M, '__call__'):
nikcleju@21:         w = MH(M(x));
nikcleju@21:     else:    ## M and MH are matrices
nikcleju@27:         w = numpy.dot(numpy.dot(MH, M), x);
nikcleju@21:     
nikcleju@21:     if hasattr(Omega, '__call__'):
nikcleju@21:         v = Omega(x);
nikcleju@27:         vt = numpy.zeros(v.size);
nikcleju@21:         vt[L] = v[L].copy();
nikcleju@21:         w = w + lm*OmegaH(vt);
nikcleju@21:     else:    ## Omega is assumed to be a matrix and OmegaH is its conjugate transpose
nikcleju@27:         w = w + lm*numpy.dot(numpy.dot(OmegaH[:, L],Omega[L, :]),x);
nikcleju@21:     
nikcleju@21:     return w
nikcleju@21: 
nikcleju@21: def GAP(y, M, MH, Omega, OmegaH, params, xinit):
nikcleju@21:   #function [xhat, Lambdahat] = GAP(y, M, MH, Omega, OmegaH, params, xinit)
nikcleju@21:   
nikcleju@21:   ##
nikcleju@21:   # [xhat, Lambdahat] = GAP(y, M, MH, Omega, OmegaH, params, xinit)
nikcleju@21:   #
nikcleju@21:   # Greedy Analysis Pursuit Algorithm
nikcleju@21:   # This aims to find an approximate (sometimes exact) solution of
nikcleju@21:   #    xhat = argmin || Omega * x ||_0   subject to   || y - M * x ||_2 <= epsilon.
nikcleju@21:   #
nikcleju@21:   # Outputs:
nikcleju@21:   #   xhat : estimate of the target cosparse vector x0.
nikcleju@21:   #   Lambdahat : estimate of the cosupport of x0.
nikcleju@21:   #
nikcleju@21:   # Inputs:
nikcleju@21:   #   y : observation/measurement vector of a target cosparse solution x0,
nikcleju@21:   #       given by relation  y = M * x0 + noise.
nikcleju@21:   #   M : measurement matrix. This should be given either as a matrix or as a function handle
nikcleju@21:   #       which implements linear transformation.
nikcleju@21:   #   MH : conjugate transpose of M. 
nikcleju@21:   #   Omega : analysis operator. Like M, this should be given either as a matrix or as a function
nikcleju@21:   #           handle which implements linear transformation.
nikcleju@21:   #   OmegaH : conjugate transpose of OmegaH.
nikcleju@21:   #   params : parameters that govern the behavior of the algorithm (mostly).
nikcleju@21:   #      params.num_iteration : GAP performs this number of iterations.
nikcleju@21:   #      params.greedy_level : determines how many rows of Omega GAP eliminates at each iteration.
nikcleju@21:   #                            if the value is < 1, then the rows to be eliminated are determined by
nikcleju@21:   #                                j : |omega_j * xhat| > greedy_level * max_i |omega_i * xhat|.
nikcleju@21:   #                            if the value is >= 1, then greedy_level is the number of rows to be
nikcleju@21:   #                            eliminated at each iteration.
nikcleju@21:   #      params.stopping_coefficient_size : when the maximum analysis coefficient is smaller than
nikcleju@21:   #                                         this, GAP terminates.
nikcleju@21:   #      params.l2solver : legitimate values are 'pseudoinverse' or 'cg'. determines which method
nikcleju@21:   #                        is used to compute
nikcleju@21:   #                        argmin || Omega_Lambdahat * x ||_2   subject to  || y - M * x ||_2 <= epsilon.
nikcleju@21:   #      params.l2_accuracy : when l2solver is 'cg', this determines how accurately the above 
nikcleju@21:   #                           problem is solved.
nikcleju@21:   #      params.noise_level : this corresponds to epsilon above.
nikcleju@21:   #   xinit : initial estimate of x0 that GAP will start with. can be zeros(d, 1).
nikcleju@21:   #
nikcleju@21:   # Examples:
nikcleju@21:   #
nikcleju@21:   # Not particularly interesting:
nikcleju@21:   # >> d = 100; p = 110; m = 60; 
nikcleju@21:   # >> M = randn(m, d);
nikcleju@21:   # >> Omega = randn(p, d);
nikcleju@21:   # >> y = M * x0 + noise;
nikcleju@21:   # >> params.num_iteration = 40;
nikcleju@21:   # >> params.greedy_level = 0.9;
nikcleju@21:   # >> params.stopping_coefficient_size = 1e-4;
nikcleju@21:   # >> params.l2solver = 'pseudoinverse';
nikcleju@21:   # >> [xhat, Lambdahat] = GAP(y, M, M', Omega, Omega', params, zeros(d, 1));
nikcleju@21:   #
nikcleju@21:   # Assuming that FourierSampling.m, FourierSamplingH.m, FDAnalysis.m, etc. exist:
nikcleju@21:   # >> n = 128;
nikcleju@21:   # >> M = @(t) FourierSampling(t, n);
nikcleju@21:   # >> MH = @(u) FourierSamplingH(u, n);
nikcleju@21:   # >> Omega = @(t) FDAnalysis(t, n);
nikcleju@21:   # >> OmegaH = @(u) FDSynthesis(t, n);
nikcleju@21:   # >> params.num_iteration = 1000;
nikcleju@21:   # >> params.greedy_level = 50;
nikcleju@21:   # >> params.stopping_coefficient_size = 1e-5;
nikcleju@21:   # >> params.l2solver = 'cg';   # in fact, 'pseudoinverse' does not even make sense.
nikcleju@21:   # >> [xhat, Lambdahat] = GAP(y, M, MH, Omega, OmegaH, params, zeros(d, 1));
nikcleju@21:   #
nikcleju@21:   # Above: FourierSampling and FourierSamplingH are conjugate transpose of each other.
nikcleju@21:   #        FDAnalysis and FDSynthesis are conjugate transpose of each other.
nikcleju@21:   #        These routines are problem specific and need to be implemented by the user.
nikcleju@21:   
nikcleju@21:   #d = length(xinit(:));
nikcleju@21:   d = xinit.size
nikcleju@21:   
nikcleju@21:   #if strcmp(class(Omega), 'function_handle')
nikcleju@21:   #    p = length(Omega(zeros(d,1)));
nikcleju@21:   #else    ## Omega is a matrix
nikcleju@21:   #    p = size(Omega, 1);
nikcleju@21:   #end
nikcleju@21:   if hasattr(Omega, '__call__'):
nikcleju@27:       p = Omega(numpy.zeros((d,1)))
nikcleju@21:   else:
nikcleju@21:       p = Omega.shape[0]
nikcleju@21:   
nikcleju@21:   
nikcleju@21:   iter = 0
nikcleju@21:   lagmult = 1e-4
nikcleju@21:   #Lambdahat = 1:p;
nikcleju@27:   Lambdahat = numpy.arange(p)
nikcleju@21:   #while iter < params.num_iteration
nikcleju@21:   while iter < params["num_iteration"]:
nikcleju@21:       iter = iter + 1
nikcleju@21:       #[xhat, analysis_repr, lagmult] = ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, lagmult, params);
nikcleju@21:       xhat,analysis_repr,lagmult = ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, lagmult, params)
nikcleju@21:       #[to_be_removed, maxcoef] = FindRowsToRemove(analysis_repr, params.greedy_level);
nikcleju@27:       to_be_removed,maxcoef = FindRowsToRemove(analysis_repr, params["greedy_level"])
nikcleju@21:       #disp(['** maxcoef=', num2str(maxcoef), ' target=', num2str(params.stopping_coefficient_size), ' rows_remaining=', num2str(length(Lambdahat)), ' lagmult=', num2str(lagmult)]);
nikcleju@21:       #print '** maxcoef=',maxcoef,' target=',params['stopping_coefficient_size'],' rows_remaining=',Lambdahat.size,' lagmult=',lagmult
nikcleju@21:       if check_stopping_criteria(xhat, xinit, maxcoef, lagmult, Lambdahat, params):
nikcleju@21:           break
nikcleju@21: 
nikcleju@21:       xinit = xhat.copy()
nikcleju@21:       #Lambdahat[to_be_removed] = []
nikcleju@27:       Lambdahat = numpy.delete(Lambdahat.squeeze(),to_be_removed)
nikcleju@21:   
nikcleju@21:       #n = sqrt(d);
nikcleju@21:       #figure(9);
nikcleju@21:       #RR = zeros(2*n, n-1);
nikcleju@21:       #RR(Lambdahat) = 1;
nikcleju@21:       #XD = ones(n, n);
nikcleju@21:       #XD(:, 2:end) = XD(:, 2:end) .* RR(1:n, :);
nikcleju@21:       #XD(:, 1:(end-1)) = XD(:, 1:(end-1)) .* RR(1:n, :);
nikcleju@21:       #XD(2:end, :) = XD(2:end, :) .* RR((n+1):(2*n), :)';
nikcleju@21:       #XD(1:(end-1), :) = XD(1:(end-1), :) .* RR((n+1):(2*n), :)';
nikcleju@21:       #XD = FD2DiagnosisPlot(n, Lambdahat);
nikcleju@21:       #imshow(XD);
nikcleju@21:       #figure(10);
nikcleju@21:       #imshow(reshape(real(xhat), n, n));
nikcleju@21:   
nikcleju@21:   #return;
nikcleju@27:   return xhat,Lambdahat
nikcleju@21:  
nikcleju@21: def FindRowsToRemove(analysis_repr, greedy_level):
nikcleju@21: #function [to_be_removed, maxcoef] = FindRowsToRemove(analysis_repr, greedy_level)
nikcleju@21: 
nikcleju@21:     #abscoef = abs(analysis_repr(:));
nikcleju@27:     abscoef = numpy.abs(analysis_repr)
nikcleju@21:     #n = length(abscoef);
nikcleju@21:     n = abscoef.size
nikcleju@21:     #maxcoef = max(abscoef);
nikcleju@21:     maxcoef = abscoef.max()
nikcleju@21:     if greedy_level >= 1:
nikcleju@21:         #qq = quantile(abscoef, 1.0-greedy_level/n);
nikcleju@21:         qq = sp.stats.mstats.mquantile(abscoef, 1.0-greedy_level/n, 0.5, 0.5)        
nikcleju@21:     else:
nikcleju@21:         qq = maxcoef*greedy_level
nikcleju@21: 
nikcleju@21:     #to_be_removed = find(abscoef >= qq);
nikcleju@27:     # [0] needed because nonzero() returns a tuple of arrays!
nikcleju@27:     to_be_removed = numpy.nonzero(abscoef >= qq)[0]
nikcleju@21:     #return;
nikcleju@27:     return to_be_removed,maxcoef
nikcleju@21: 
nikcleju@21: def check_stopping_criteria(xhat, xinit, maxcoef, lagmult, Lambdahat, params):
nikcleju@21: #function r = check_stopping_criteria(xhat, xinit, maxcoef, lagmult, Lambdahat, params)
nikcleju@21: 
nikcleju@21:     #if isfield(params, 'stopping_coefficient_size') && maxcoef < params.stopping_coefficient_size
nikcleju@21:     if ('stopping_coefficient_size' in params) and maxcoef < params['stopping_coefficient_size']:
nikcleju@21:         return 1
nikcleju@21: 
nikcleju@21:     #if isfield(params, 'stopping_lagrange_multiplier_size') && lagmult > params.stopping_lagrange_multiplier_size
nikcleju@21:     if ('stopping_lagrange_multiplier_size' in params) and lagmult > params['stopping_lagrange_multiplier_size']:
nikcleju@21:         return 1
nikcleju@21: 
nikcleju@21:     #if isfield(params, 'stopping_relative_solution_change') && norm(xhat-xinit)/norm(xhat) < params.stopping_relative_solution_change
nikcleju@27:     if ('stopping_relative_solution_change' in params) and numpy.linalg.norm(xhat-xinit)/numpy.linalg.norm(xhat) < params['stopping_relative_solution_change']:
nikcleju@21:         return 1
nikcleju@21: 
nikcleju@21:     #if isfield(params, 'stopping_cosparsity') && length(Lambdahat) < params.stopping_cosparsity
nikcleju@21:     if ('stopping_cosparsity' in params) and Lambdahat.size() < params['stopping_cosparsity']:
nikcleju@21:         return 1
nikcleju@21:     
nikcleju@21:     return 0