g++ error (involving sum I believe)

[nbecker]

The following code results in a g++ error
(and the longest error message in history!)
https://gist.github.com/nbecker/37393e16d5ae8044e37e80b040c4e548

[serge-sans-paille]

14147 lines of error :-) mmmmh, tasty!

[nbecker]

Did I set the record? Is there a prize??

[serge-sans-paille]

:-) That's not a metric I'm keeping track of :-)

there's two constructs pythran does not support (yet) in tour code:

• indexing with a list and an integer (as in Lq[i,Nij])
• removing from a list view (as in Nij.remove(j)) # actually pythran considers this as a view but Python does not

Thanks for providing such a delicate test case, food for thoughts :-)

[serge-sans-paille]

The following compiles correctly (but is the result correct?) with an upt-o-date master

import numpy as np

#pythran export BinaryProduct(int[:,:], int[:])
def BinaryProduct(X,Y):
        
    """ Binary Matrices or Matrix-vector product in Z/2Z. Works with scipy.sparse.csr_matrix matrices X,Y too."""
 
    A = X.dot(Y)
    
    return A%2

#pythran export InCode(int[:,:], int[:])
def InCode(H,x):

    """ Computes Binary Product of H and x. If product is null, x is in the code.
        Returns appartenance boolean. 
    """
        
    return  (BinaryProduct(H,x)==0).all()

#pythran export Decoding_logBP(int[:,:], int list list, int list list, float[:,:], float[:], int)
def Decoding_logBP(H,Bits,Nodes,Lq,Lc,max_iter):
    """ Decoding function using Belief Propagation algorithm (logarithmic version)
    IMPORTANT: if H is large (n>1000), H should be scipy.sparse.csr_matrix object to speed up calculations
    (highly recommanded. )
    -----------------------------------
    Parameters:
    
    H: 2D-array (OR scipy.sparse.csr_matrix object) Parity check matrix, shape = (m,n) 
    y: n-vector recieved after transmission in the channel. (In general, returned 
    by Coding Function)
    Signal-Noise Ratio: SNR = 10log(1/variance) in decibels of the AWGN used in coding.
    
    max_iter: (default = 1) max iterations of the main loop. Increase if decoding is not error-free.
     """
        
    m,n=H.shape

    if not len(Lc)==n:
        raise ValueError('La taille de y doit correspondre au nombre de colonnes de H')

    if m>=n:
        raise ValueError('H doit avoir plus de colonnes que de lignes')
    
    
    # var = 10**(-SNR/10)

    # ### ETAPE 0: initialisation 
    
    # Lc = 2*y/var

    # Lq=np.zeros(shape=(m,n))

    Lr = np.zeros(shape=(m,n))
    
    count=0
    
    prod=np.prod
    tanh = np.tanh
    log = np.log
    
    Lq += Lc
    
    #Bits,Nodes = BitsAndNodes(H)
    while(True):

        count+=1 #Compteur qui empêche la boucle d'être infinie .. 

        #### ETAPE 1 : Horizontale
        for i in range(m):
            Ni = Bits[i]
            for j in Ni:
                Nij = list(Ni)
                Nij.remove(j)
                X = prod(tanh(0.5*np.array(Lq[i])[Nij]))
                num = 1 + X
                denom = 1 - X
                # if num == 0: 
                #     Lr[i,j] = -1
                # elif denom  == 0:
                #     Lr[i,j] =  1
                # else: 
                #     Lr[i,j] = log(num/denom)
                Lr[i,j] = log(num) - log(denom)
        Lr = np.clip (Lr, -100, 100)

        #### ETAPE 2 : Verticale
        for j in range(n):
            Mj = Nodes[j]
            
            for i in Mj:
                Mji = list(Mj)
                Mji.remove(i)

                Lq[i,j] = Lc[j]+sum(np.array(Lr[Mji])[:,j])
        
 
        #### LLR a posteriori:
        # L_posteriori = np.zeros(n)
        # for j in range(n):
        #     Mj = Nodes[j]

        #     L_posteriori[j] = Lc[j] + sum(Lr[Mj,j])
        extrinsic = np.empty(n)
        for j in range(n):
            Mj = Nodes[j]

            extrinsic[j] = sum(np.array(Lr[Mj])[:,j])

        L_posteriori = extrinsic + Lc
        #x = np.array(L_posteriori <= 0).astype(int)
        x = np.array(extrinsic <= 0).astype(int)
        product = InCode(H,x)
        #print (count, product)

        if product or count >= max_iter:
            break
    # print (count)
    return np.array(L_posteriori <= 0).astype(int), Lq - Lc, extrinsic, product

[nbecker]

Thanks. Compiles using current released version, but output is not correct :( Also compiles with replacing e.g., Lr[Mj,j] -> Lr[Mj][j]. Shouldn't that work?

…

On Wed, Mar 28, 2018 at 3:55 PM serge-sans-paille ***@***.***> wrote: The following compiles correctly (but is the result correct?) with an upt-o-date master import numpy as np #pythran export BinaryProduct(int[:,:], int[:]) def BinaryProduct(X,Y): """ Binary Matrices or Matrix-vector product in Z/2Z. Works with scipy.sparse.csr_matrix matrices X,Y too.""" A = X.dot(Y) return A%2 #pythran export InCode(int[:,:], int[:]) def InCode(H,x): """ Computes Binary Product of H and x. If product is null, x is in the code. Returns appartenance boolean. """ return (BinaryProduct(H,x)==0).all() #pythran export Decoding_logBP(int[:,:], int list list, int list list, float[:,:], float[:], int) def Decoding_logBP(H,Bits,Nodes,Lq,Lc,max_iter): """ Decoding function using Belief Propagation algorithm (logarithmic version) IMPORTANT: if H is large (n>1000), H should be scipy.sparse.csr_matrix object to speed up calculations (highly recommanded. ) ----------------------------------- Parameters: H: 2D-array (OR scipy.sparse.csr_matrix object) Parity check matrix, shape = (m,n) y: n-vector recieved after transmission in the channel. (In general, returned by Coding Function) Signal-Noise Ratio: SNR = 10log(1/variance) in decibels of the AWGN used in coding. max_iter: (default = 1) max iterations of the main loop. Increase if decoding is not error-free. """ m,n=H.shape if not len(Lc)==n: raise ValueError('La taille de y doit correspondre au nombre de colonnes de H') if m>=n: raise ValueError('H doit avoir plus de colonnes que de lignes') # var = 10**(-SNR/10) # ### ETAPE 0: initialisation # Lc = 2*y/var # Lq=np.zeros(shape=(m,n)) Lr = np.zeros(shape=(m,n)) count=0 prod=np.prod tanh = np.tanh log = np.log Lq += Lc #Bits,Nodes = BitsAndNodes(H) while(True): count+=1 #Compteur qui empêche la boucle d'être infinie .. #### ETAPE 1 : Horizontale for i in range(m): Ni = Bits[i] for j in Ni: Nij = list(Ni) Nij.remove(j) X = prod(tanh(0.5*np.array(Lq[i])[Nij])) num = 1 + X denom = 1 - X # if num == 0: # Lr[i,j] = -1 # elif denom == 0: # Lr[i,j] = 1 # else: # Lr[i,j] = log(num/denom) Lr[i,j] = log(num) - log(denom) Lr = np.clip (Lr, -100, 100) #### ETAPE 2 : Verticale for j in range(n): Mj = Nodes[j] for i in Mj: Mji = list(Mj) Mji.remove(i) Lq[i,j] = Lc[j]+sum(np.array(Lr[Mji])[:,j]) #### LLR a posteriori: # L_posteriori = np.zeros(n) # for j in range(n): # Mj = Nodes[j] # L_posteriori[j] = Lc[j] + sum(Lr[Mj,j]) extrinsic = np.empty(n) for j in range(n): Mj = Nodes[j] extrinsic[j] = sum(np.array(Lr[Mj])[:,j]) L_posteriori = extrinsic + Lc #x = np.array(L_posteriori <= 0).astype(int) x = np.array(extrinsic <= 0).astype(int) product = InCode(H,x) #print (count, product) if product or count >= max_iter: break # print (count) return np.array(L_posteriori <= 0).astype(int), Lq - Lc, extrinsic, product — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#833 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAHK0A-9tTtXMA4JqvY4PwsTL118WqeCks5ti-rHgaJpZM4S-32J> .

[nbecker]

Good news, I've got it working. The time spent in my test program in the decoder function dropped from 95% in the best python version to 89% in the pythran version.
The working code is here https://gist.github.com/nbecker/a1085d485a1e63d6126820a255302d3e

Along the way, I found that if I mistakenly wrote e.g.,
extrinsic[j] = sum(np.array(Lr[Mj][j]))
instead of:
extrinsic[j] = sum(np.array(Lr[Mj][:,j]))
that the code segfaults at runtime. I don't suppose pythran could catch that?

[nbecker]

I made a little faster by removing common subexpressions for the arrays:
https://gist.github.com/nbecker/472edc3bdf1365a7e41db1ec1b50e2ce

I don't know why I need to use e.g.,
Lrj = np.array(Lr[:,j])
...
Lq[i,j] = Lc[j]+sum(np.array(Lrj[Mji]))
Don't know why np.array(...) is needed here in either case, but doesn't compile without

[serge-sans-paille]

@nbecker I improved your code while improving pythran too, check the HEAD of fix/ldpc-decoder branch and the source code in pythran/tests/cases/ldpc-decoder.py and tell me about the speedups :-)

[nbecker]

looking at the test case, I guess we still have to write
G(Lq[i][Nij]
instead of
G(Lq[i,Nij])
correct?

[serge-sans-paille]

correct!

[nbecker]

Here's the latest version
https://gist.github.com/nbecker/c00202dc72a343a27923a1ad408c9e90
In 2 places, I factored out common subexpressions from the inner for loops.

The code compiles (yay!) with latest pythran from git (yay!) and speed has improved.

[serge-sans-paille]

Can you provide me with relevant arguments so that I can run benchmarks on my own? I tried that ``` python -m timeit -s 'import nbb, numpy as np; x = np.arange(200).reshape((10,20)); xp = np.arange(200.).reshape((10,20)); y = [[i/2]*20 for i in range(20)]; z = np.arange(20.)' 'nbb.Decoding_logBP(x, y, y, xp, z, 2)' ```

[nbecker]

Hopefully, this has everything needed to try running ldpc_nr.py, which is what I used for benchmarking

ldpc.zip

It should print
True
True
0
10 times and exit.

No, on review I see more dependencies. Let me try to clean that up.

[nbecker]

OK, try this one. You should only also need to install pyldpc, but that is available via pip (pypi)
ldpc2.zip

[nbecker]

I see the function phi0 translates to
return (-pythonic::math::functor::log{}(pythonic::math::functor::tanh{}((pythonic::operator_::div(x__, 2L)))));

The use of '2L' seems suspicious here, this is a floating divide, not integer. Does it matter?

[serge-sans-paille]

The use of '2L' seems suspicious here, this is a floating divide, not integer. Does it matter?

No it does not :-)

[serge-sans-paille]

Can you confirm #838 fixes these issues?

[nbecker]

The version I'm testing has tags default/fix/ldpc-decoder and includes "Various view improvements". I'm afraid I'm using mercurial not git (via hg-git plugin), so it's a little hard to compare versions but I think this is the version you're asking about. Still haven't become familiar with using git.

This one is working with the test code I sent. However if I comment out the Nij = list(Ni) and change it to Nij = Ni, and similarly for Mji = Mj, the the python process segfaults. I was checking if the list() was needed, although now thinking about it, it is needed in this version, but it shouldn't segfault.