nikcleju@15: # -*- coding: utf-8 -*-
nikcleju@15: """
nikcleju@15: Created on Sat Nov 05 21:29:09 2011
nikcleju@15: 
nikcleju@15: @author: Nic
nikcleju@15: """
nikcleju@15: 
nikcleju@15: # -*- coding: utf-8 -*-
nikcleju@15: """
nikcleju@15: Created on Sat Nov 05 18:39:54 2011
nikcleju@15: 
nikcleju@15: @author: Nic
nikcleju@15: """
nikcleju@15: 
nikcleju@15: import numpy as np
nikcleju@15: 
nikcleju@15: #function s=SL0(A, x, sigma_min, sigma_decrease_factor, mu_0, L, A_pinv, true_s)
nikcleju@15: def SL0_approx(A, x, eps, sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, A_pinv=None, true_s=None):
nikcleju@15:   
nikcleju@15:   if A_pinv is None:
nikcleju@15:     A_pinv = np.linalg.pinv(A)
nikcleju@15:   
nikcleju@15:   if true_s is not None:
nikcleju@15:       ShowProgress = True
nikcleju@15:   else:
nikcleju@15:       ShowProgress = False
nikcleju@15:   
nikcleju@15:   # Initialization
nikcleju@15:   #s = A\x;
nikcleju@15:   s = np.dot(A_pinv,x)
nikcleju@15:   sigma = 2.0 * np.abs(s).max()
nikcleju@15:   
nikcleju@15:   # Main Loop
nikcleju@15:   while sigma>sigma_min:
nikcleju@15:       for i in np.arange(L):
nikcleju@15:           delta = OurDelta(s,sigma)
nikcleju@15:           s = s - mu_0*delta
nikcleju@15:           # At this point, s no longer exactly satisfies x = A*s
nikcleju@15:           # The original SL0 algorithm projects s onto {s | x = As} with
nikcleju@15:           # s = s - np.dot(A_pinv,(np.dot(A,s)-x))   # Projection
nikcleju@15:           # We want to project s onto {s | |x-As| < eps}
nikcleju@15:           # We move onto the direction -A_pinv*(A*s-x), but only with a
nikcleju@15:           # smaller step:
nikcleju@15:           direction = np.dot(A_pinv,(np.dot(A,s)-x))
nikcleju@15:           if (np.linalg.norm(np.dot(A,direction)) >= eps):
nikcleju@15:             s = s - (1.0 - eps/np.linalg.norm(np.dot(A,direction))) * direction
nikcleju@15: 
nikcleju@15:           #assert(np.linalg.norm(x - np.dot(A,s)) < eps + 1e-6)          
nikcleju@15:       
nikcleju@15:       if ShowProgress:
nikcleju@15:           #fprintf('     sigma=#f, SNR=#f\n',sigma,estimate_SNR(s,true_s))
nikcleju@15:           string = '     sigma=%f, SNR=%f\n' % sigma,estimate_SNR(s,true_s)
nikcleju@15:           print string
nikcleju@15:       
nikcleju@15:       sigma = sigma * sigma_decrease_factor
nikcleju@15:   
nikcleju@15:   return s
nikcleju@15:   
nikcleju@15:   
nikcleju@15: ####################################################################
nikcleju@15: #function delta=OurDelta(s,sigma)
nikcleju@15: def OurDelta(s,sigma):
nikcleju@15:   
nikcleju@15:   return s * np.exp( (-np.abs(s)**2) / sigma**2)
nikcleju@15:   
nikcleju@15: ####################################################################
nikcleju@15: #function SNR=estimate_SNR(estim_s,true_s)
nikcleju@15: def estimate_SNR(estim_s, true_s):
nikcleju@15:   
nikcleju@15:   err = true_s - estim_s
nikcleju@15:   return 10*np.log10((true_s**2).sum()/(err**2).sum())
nikcleju@15: