[SciPy-User] [SciPy-user] Speedup with integrate.romberg() and vector arguments

Fri Apr 29 02:44:38 EDT 2011

Thank you for your tips. A fixed grid wouldn't be a real option.

But by writing the romberg() loop your way, the time for one curve_fit()
call just takes 3.8 s instead of 15.7 s. 

def phi_romb_vec(q, a, alpha, beta, r_0, r=r): # phi(q,...,r)
    y = array([integrate.romberg(dy_dr, r[0], r[-1], # dy(r,...)
         args=(qi, A, a, alpha, beta, r_0), vec_func=True) for qi in q])
    return y        

Interestingly, quad() "vectorized" doesn't change a thing: 9.24 s --> 9.22
s.

I am new to this and still have to figure out what you do with
np.concatenate(). It didn't seem to work that way.

Regards,
Arne

josef.pktd wrote:
> 
> On Wed, Apr 27, 2011 at 9:14 AM,  <josef.pktd at gmail.com> wrote:
>> On Wed, Apr 27, 2011 at 4:50 AM,  <arsbbr at gmx.net> wrote:
>>> Hi,
>>>
>>> I need to integrate a function with an [1xn] array argument q.
>>> For this I used integrate.quad(), but then I have to loop through every
>>> element in q.
>>>
>>> In http://www.astrobetter.com/interpolation-and-integration-in-python/ I
>>> found a hint how to speed up the integration.
>>>
>>> Unfortunately I can't seem to get romberg() to work:
>>>
>>> q is a [1 x n] numpy array, so I went from
>>>
>>> def phi_quad(q, a, alpha, beta, r_0, r=r): # phi(q,...,r)
>>>    y = zeros(r.size)
>>>    for k in range(r.size):
>>>        y[k] = integrate.quad(dy_dr, r[0], r[-1], # dy(r,...)
>>>            args=(q[k], A, a, alpha, beta, r_0))[0]
>>>    return y
>>>
>>> to
>>>
>>> def phi_romb_vec(q, a, alpha, beta, r_0, r=r): # phi(q,...,r)
>>>    y = integrate.romberg(dy_dr, r[0], r[-1], # dy(r,...)
>>>         args=(q, A, a, alpha, beta, r_0), vec_func=True)
>>>    return y
>>>
>>> which raises an error, if q is an array. The other arguments are all
>>> floats. I attached a working example of the problem. Is there any way to
>>> speed up the integration?
>>
>> I never tried this, but from the blog and the documentation, this is
>> my interpretation
>>
>> It looks like romberg takes a vector argument but handles only a
>> single integration problem (scalar output).
>>
>> in the function to be integrated, romberg uses a vector of r to speed
>> up the integration
>> dy_dr(r, q, A, a, alpha, beta, r_0)
>> however, q and the other parameters need to be scalars
>>
>> So you still need to loop over q, doing one integration problem at a
>> time.
>>
>> It might be possible to enhance romberg to proper broadcasting, but I
>> think it would require that the number of iterations, number of
>> subdivisions is the same for all integration problems.
> 
> this actually seems to work with minimal changes
> if np.all(err < tol) or np.all(err < rtol*abs(result)):
> 
> but needs some broadcasting checks.
> 
>>>> np.concatenate([phi_romb_vec(qi, a, alpha, beta, r_0, r=r) for qi in
>>>> q])
> array([ 0.00065848,  0.00127215,  0.00216926,  0.00339121,  0.00497182,
>         0.00693667,  0.0093027 ,  0.01207799,  0.01526184,  0.01884502,
>         0.02281024,  0.02713285,  0.03178164,  0.03671976,  0.04190578,
>         0.0472947 ,  0.05283912,  0.05849023,  0.06419887,  0.06991646,
>         0.07559587,  0.08119216,  0.08666319,  0.09197017,  0.09707798,
>         0.1019555 ,  0.1065757 ,  0.11091574,  0.11495692,  0.11868458,
>         0.12208791,  0.12515975,  0.12789636,  0.13029711,  0.13236424,
>         0.13410255,  0.1355191 ,  0.13662298,  0.137425  ,  0.13793743,
>         0.13817382,  0.13814868,  0.13787734,  0.13737571,  0.13666008,
>         0.13574698,  0.13465298,  0.13339455,  0.13198792,  0.13044894,
>         0.128793  ,  0.12703493,  0.12518888,  0.12326832,  0.12128595,
>         0.11925369,  0.11718264])
>>>> phi_romb_vec(q[None,:], a, alpha, beta, r_0, r=r)
> array([[ 0.00065848,  0.00127215,  0.00216926,  0.00339121,  0.00497182,
>          0.00693667,  0.0093027 ,  0.01207799,  0.01526184,  0.01884502,
>          0.02281024,  0.02713285,  0.03178164,  0.03671976,  0.04190578,
>          0.0472947 ,  0.05283912,  0.05849023,  0.06419887,  0.06991646,
>          0.07559587,  0.08119216,  0.08666319,  0.09197017,  0.09707798,
>          0.1019555 ,  0.1065757 ,  0.11091574,  0.11495692,  0.11868458,
>          0.1220879 ,  0.12515975,  0.12789636,  0.13029711,  0.13236424,
>          0.13410255,  0.1355191 ,  0.13662298,  0.137425  ,  0.13793743,
>          0.13817382,  0.13814868,  0.13787734,  0.13737571,  0.13666008,
>          0.13574698,  0.13465298,  0.13339455,  0.13198792,  0.13044894,
>          0.128793  ,  0.12703493,  0.12518888,  0.12326832,  0.12128595,
>          0.11925369,  0.11718264]])
> 
> I haven't managed the call to curve_fit with the vectorized romberg
> 
> Josef
> 
>>
>> If you are willing to work with a fixed grid, then integrate.romb
>> might be useful, since from the docs it takes an axis argument, so a
>> broadcasted (r, q) grid might work.
>>
>> Josef
>>
>>
>>>
>>> Best regards,
>>> Arne
>>> --
>>> NEU: FreePhone - kostenlos mobil telefonieren und surfen!
>>> Jetzt informieren: http://www.gmx.net/de/go/freephone
>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>>>
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
> 
> 

-- 
View this message in context: http://old.nabble.com/Speedup-with-integrate.romberg%28%29-and-vector-arguments-tp31485594p31502972.html
Sent from the Scipy-User mailing list archive at Nabble.com.