Search code examples
pythonmathsumintegralnotation

Need help making a sum calculation faster in python


Here is a piece of my code with the summation. It includes an integral. all the variables peppered in there have float or integer values.

import numpy as np
import math
import scipy.integrate
from mpmath import nsum, inf

sum1 = nsum(lambda m: (1/(math.factorial(m))*((ra/a1)**(2*m))*(1/(math.factorial((m+1)-1)))*(scipy.integrate.quad(lambda u:((u**((m+1)-1))*(math.exp(-u))) ,0,(d**2)/(4*(a1**2)),limit = 50))[0]),[0,100])

Is there any way I can make this faster? I need to be able to loop through this several hundred times. For the moment it's taking more than an hour. Here is the original complete equation. The part of my code that I pasted here is just the first sum (sigma) in the equation where m goes from 0 to infinity. I changed it to m from 0 to 100 in my code here temporarily.

original equation

Many thanks!

Update:

I've tried a different way. I think those who've commented were right and I was using scipy incorrectly before. I think I've fixed that. It is still slowish for my purposes, but now only takes 30min for several hundred repetitions. (I'm also guessing you can use either the "def" structure or "lambda" structure to pass through quad). I'd appreciate if you have any advice for making something like this faster, but I appreciate you guys helping lead to the answer.

for m in range(0,100):

    # def integrand(u,m):
    #     return (u**((m+1)-1))*(math.exp(-u)) 

    f = lambda u,m:(u**((m+1)-1))*(math.exp(-u)) 
    
    sumfunc =1/(math.factorial(m))*((ra/a1)**(2*m))*(1/(math.factorial((m+1)-1)))*scipy.integrate.quad(f,0,(d**2)/(4*(a1**2)),args=(m,))[0] 
    
    tosum.append(sumfunc)
    
sum1 = sum(tosum)

Solution

    1. mpmath knows gammainc, use it;
    2. npsum is smart, don't try to be smarter, just sum to infinity. It'll do some fancy computations to stop just when it needs to;
    3. if you care about precision, try to use mpmath types as much as possible

    So here's the code:

    from mpmath import *
    
    a1 = mpf('0.62')
    ra = mpf('0.99')
    
    def the_sum(d, a, **kwargs): # added `a` for reuse in the complete formula
        r_ = (ra / a) ** 2
        d_ = d ** 2 / (4 * a ** 2)
        return nsum(
            lambda m: r_ ** m / factorial(m) * gammainc(m + 1, 0, d_, regularized=True),
            [0, inf],
            **kwargs
        )
    

    and some benchmark

    %time [sum_OP(d) for d in range(5)]                                                                                                                                                                
    CPU times: user 684 ms, sys: 2.97 ms, total: 687 ms
    Wall time: 685 ms
    [mpf('0.0'),
     mpf('0.93752967676144794'),
     mpf('5.3589089548240878'),
     mpf('10.711751972234245'),
     mpf('12.600516476806067')]
    
    %time [the_sum(d, a1) for d in range(5)]                                                                                                                                                               
    CPU times: user 28.7 ms, sys: 1.3 ms, total: 30 ms
    Wall time: 28.8 ms
    [mpf('0.0'),
     mpf('0.93752967676144794'),
     mpf('5.3589089548240905'),
     mpf('10.711751972234246'),
     mpf('12.600516476806067')]
    

    If you really want to makes things faster, want to be smarter than nsum, and are confident it will work in this case (I do), you can tell it to use the direct method:

     %time [the_sum(d, a1, method='d') for d in range(5)]                                                                                                                                               
    CPU times: user 20.7 ms, sys: 613 µs, total: 21.3 ms
    Wall time: 20.9 ms
    [mpf('0.0'),
     mpf('0.93752967676144794'),
     mpf('5.3589089548240905'),
     mpf('10.711751972234246'),
     mpf('12.600516476806067')]
    

    Finally, as it seems you don't need any multiprecision, no need to use mpmath, and stick to some low precision version:

    a1 = 0.62
    ra = 0.99
    def sum_low_prec(d, a):
        return sum(
            (ra/a)**(2*m) / (math.factorial(m)) * scipy.special.gammainc(m + 1, (d**2)/(4*(a**2)))
             for m in range(30)
        )
    
    %time [sum_low_prec(d, a1) for d in range(5)]                                                                                                                                                     
    CPU times: user 871 µs, sys: 3 µs, total: 874 µs
    Wall time: 877 µs
    [0.0,
     0.9375296767614478,
     5.358908954824089,
     10.711751972234243,
     12.600516476806066]