[SciPy-user] Help with array functions

[eric jones]

So, try two...

I don't know of a pure Python solution either.  My next option is 
usually playing around in weave.  Here is a weave snippet that 
implements the algorithm in C.  It runs, but just like the last one is 
untested.  Given the track record here, a decent test is advisable... :-) 

Here is the output from my cursory check to see if the code would work:

C:\temp>wv.py
weave results:  [ 0. 0.8 0.96 0.992 0.9984  0.99968  0.999936  
0.9999872  0.99999744  1.]
python results: [ 0. 0.8 0.96 0.992 0.9984  0.99968  0.999936  
0.9999872  0.99999744  1.]
C time: 0.00707567096418
Python time: 0.316835262617
weave vs. Python speedup: 44.7781226997

This is using gcc as the compiler in weave and 2048 samples in a 
signal.  If you drop to 256, the improvement is cut in half.

Now the code:

# test.py
from scipy import zeros, ones, Float64, arange
from scipy import stats
import weave

def python_algorithm(gain, complexout, attack, peak, level, samples):
    # gain is changed *inplace*
    i = arange(1,samples)
    for x in i:
        if abs(complexout[x-1]) > 0:
           gain[x] = ( (gain[x-1] * (1.0 - attack)) + ((attack * peak) / 
abs(complexout[x-1])) )       
           if gain[x] > level:
               gain[x] = level
        else:
           gain[x] = 1

def weave_algorithm(gain, complexout, attack, peak, level, samples):
    # gain is changed *inplace*
    code = """
           for (int x=1; x < samples; x++)
           {
               if (abs(complexout[x-1]) > 0.0)
               {
                  gain[x] = ( (gain[x-1] * (1.0 - attack)) +
                              ((attack * peak) / abs(complexout[x-1])) 
);       
                  if (gain[x] > level)
                      gain[x] = level;
               }
               else          
                  gain[x] = 1;
           }
           """
   
    weave.inline(code,
                 ['gain', 'complexout', 'samples', 'attack', 'peak', 
'level'],
                 compiler='gcc')
   

# some default parameters...
samples = 2048
gain = zeros(samples,typecode=Float64)
complexout = stats.choice([-1.0,0.0,1.0],size=samples)

attack = .8
peak = 1.0
level = 1.0


# equivalence tests
weave_algorithm(gain, complexout, attack, peak, level, samples)
print 'weave results:', gain[:10]

gain = zeros(samples,typecode=Float64)
python_algorithm(gain, complexout, attack, peak, level, samples)
print 'python results:', gain[:10]


# speed tests
import time

N = 100

t1 = time.clock()
for it in range(N):
    weave_algorithm(gain, complexout, attack, peak, level, samples)
t2 = time.clock()
ct = t2 - t1
print 'C time:', ct

t1 = time.clock()
for it in range(N):
    python_algorithm(gain, complexout, attack, peak, level, samples)
t2 = time.clock()
pt = t2-t1
print 'Python time:', pt

ratio = pt/ct
print 'weave vs. Python speedup:', ratio


eric jones wrote:

> Yikes.  Your right!
>
> eric
>
> Perry Greenfield wrote:
>
>> Eric Jones wrote:
>>  
>>
>>> Bob.Cowdery at CGI-Europe.com wrote:
>>>
>>> The general ideas of uisng Numeric in Python is to avoid loops, even 
>>> at the expense of sometimes doing extra computations.  In this case, 
>>> the signal processing calculations can be done for every entry in 
>>> the array (except the first).  The conditions checked for in the 
>>> "if" statements can then be filtered out in subsequent array 
>>> processing steps.
>>>
>>>
>>> Here is some (totally untested) code that should give you some ideas:
>>>
>>>    gain = self.m_procGain
>>>    abs_complexout = abs(self.m_complexout)
>>>
>>>    # Do singal processing calculation for all elements in array
>>>    gain[1:] = ( (gain[:-1] * (1 - attack)) +
>>>                 ((attack * peak) / abs_complexout[:-1]) )
>>>
>>>   
>>
>> Does this do what is needed? I think the point he was making
>> was that gain[i] depends on gain[i-1] which in turn depends
>> on gain[i-2] and so on. The expression above doesn't capture that
>> since from his code he doesn't show any initial value for gain
>> (implicitly the gain[0] must exist). So using an array expression
>> which uses a shifted gain array in its input won't work.
>> This is the sort of problem that doesn't fit the array model
>> very well unless there are some primitive array operations
>> that happen to do the iterative thing needed (sometimes
>> the accumulate or convolve functions can be used for this
>> thing, but I can't see any obvious existing function. Seems
>> like a candidate for a C extension but maybe I'm missing something.
>>
>> Perry
>>  
>>
>>
>> _______________________________________________
>> SciPy-user mailing list
>> SciPy-user at scipy.net
>> http://www.scipy.net/mailman/listinfo/scipy-user
>>  
>>
>
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user