Module sigma

Functions

def is_prime(a)

Expand source code

def is_prime(a):
    if a<2:
        return False
    for i in range(2,int(a**0.5)+1) :
        if a%i==0:
            return False
    return True

def sigma_new(N, M, I)

computation of sigma using algorithm 3.4 Input: N (integer): length of the vector x M (integer): sparsity of x, len(I) I (list): the support indices of x

Output: optsig (integer): best prime among K primes < N/2 dsigma (integer): minimal distance contained in sigmaD.

Expand source code

def sigma_new(N,M,I):
    """ computation of sigma using algorithm 3.4
    Input:  N (integer): length of the vector x
            M (integer): sparsity of x, len(I)
            I (list):    the support indices of x 

    Output: optsig (integer): best prime among K primes < N/2
            dsigma (integer): minimal distance contained in sigmaD.
    """
    if M==1 or M==0:
        return 1, N-1

    # sigmas = possible primes for sigma
    sigmas=[]  
    i=int(N/2)
    K=max([k for k in range(1,M) if k*np.log2(k)<M ])
    while len(sigmas)<K and i>1:
        if is_prime(i):
            sigmas.append(i)
        i=i-1
    
    Dsigma=0
    dsigma=min([(I[(i+1)%M]-I[i])%N for i in range(M)])
    optsig=1
    for s in sigmas:
        
        #compute sigma*I sorted
        sigI=[(s*i)%N for i in I]
        sigI=sorted(sigI)
        
        #compute sigmaD (=distances between neighboring values in sigI)
        sigD=[(sigI[(i+1)%M]-sigI[i])%N for i in range(M)]
        sigmaD=[sigD[M-2]]+[sigD[M-1]]+sigD
        ks=np.argmin(sigmaD[1:M+1])+1
        
        #update Dsigma and dsigma
        Dsigma_act=max([np.abs((1/np.sin(np.pi*sigmaD[ks]/N))+np.abs(1/np.sin(np.pi*sigmaD[ks-1]/N))), np.abs((1/np.sin(np.pi*sigmaD[ks]/N))+np.abs(1/np.sin(np.pi*sigmaD[ks+1]/N)))])
        
        if Dsigma_act>Dsigma:
            optsig=s
            mindis=min(sigmaD)
            if mindis>0:
                Dsigma=Dsigma_act
                dsigma=mindis
            
        elif Dsigma_act==Dsigma:
            if abs(sum(np.exp(-2*np.pi*1j*s*i/N) for i in I))<abs(sum(np.exp(-2*np.pi*1j*optsig*i/N) for i in I)):
                optsig=s
                mindis=min(sigmaD)
                if mindis>0:
                    Dsigma=Dsigma_act
                    dsigma=mindis
    return optsig, dsigma

def sigma_opt(N, M, I)

computation of sigma using algorithm 4.5 from Plonka et al. (2018) Input: N (integer): length of the vector x M (integer): sparsity of x, len(I) I (list): the support indices of x

Output: optsig (integer): best prime among M primes < N/2 dsigma (integer): minimal distance contained in sigmaD.

Expand source code

def sigma_opt(N,M,I):
    """ computation of sigma using algorithm 4.5 from Plonka et al. (2018)
    Input:  N (integer): length of the vector x
            M (integer): sparsity of x, len(I)
            I (list):    the support indices of x 

    Output: optsig (integer): best prime among M primes < N/2
            dsigma (integer): minimal distance contained in sigmaD.
    """
    
    if M==1 or M==0:
        return 1, N-1

    # sigmas = possible primes for sigma
    sigmas=[]  
    i=int(N/2)
    while len(sigmas)<M and i>1:
        if is_prime(i):
            sigmas.append(i)
        i=i-1
   
    dsigma=min([(I[(i+1)%M]-I[i])%N for i in range(M)])
    optsig=1
    for s in sigmas:
        
        #compute sigma*I sorted
        sigI=[(s*i)%N for i in I]
        sigI=sorted(sigI)
        
        #compute sigmaD (=distances between neighboring values in sigI)
        sigmaD=[(sigI[(i+1)%M]-sigI[i])%N for i in range(M)]
       
        #update dsigma (=max(min(sigmaD)))
        if min(sigmaD)>dsigma:
            optsig=s
            dsigma=min(sigmaD)
        elif min(sigmaD)==dsigma:
            if abs(sum(np.exp(-2*np.pi*1j*s*i/N) for i in I))<abs(sum(np.exp(-2*np.pi*1j*optsig*i/N) for i in I)):
                optsig=s
                dsigma=min(sigmaD)
                
    return optsig, dsigma