I am new to community and please pardon me if I didn't provide information as intended.
This code is supposed to be creating a custom model which will be used with lmfit for curve fitting purposes. hc(q,radius,pd) is the function. It is one of the simplest functions that I will be using but even this function is taking quite a bit in python compared to Matlab. q values which are created for the sake of simplicity will be replaced by the experimental values and it will be used to fit the intensity from 1D SAXS results.
The difference between the vectorization and for loop is almost double. Although the documentation suggests to use for Loop, it gives slower results.
I am trying to learn Python, coming from Matlab. I have a very simple code for starting purposes:
from numpy import vectorize
from scipy import integrate
from scipy.special import j1
from math import sqrt, exp, pi, log
import matplotlib.pyplot as plt
import numpy as np
from numpy import empty
def plot_hc(radius, pd):
q = np.linspace(0.008, 1.0, num=500)
y = hc(q, radius, pd)
plt.loglog(q, y)
plt.show()
def hc_formfactor(q, radius):
y = (1.0 / q) * (radius * j1(q * radius))
y = y ** 2
return y
def g_distribution(z, radius, pd):
return (1 / (sqrt(2 * pi) * pd)) * exp(
-((z - radius) / (sqrt(
2) * pd)) ** 2)
def ln_distribution(z, radius, pd):
return (1 / (sqrt(2 * pi) * pd * z / radius)) * exp(
-(log(z / radius) / (sqrt(2) * pd)) ** 2)
# Dist=1(for G_Distribution)
# Dist=2(for LN Distribution)
Dist = 1
@vectorize
def hc(x, radius, pd):
global d
if Dist == 1:
nmpts = 4
va = radius - nmpts * pd
vb = radius + nmpts * pd
if va < 0:
va = 0
d = integrate.quad(lambda z: g_distribution(z, radius, pd), va, vb)
elif Dist == 2:
nmpts = 4
va = radius - nmpts * pd
vb = radius + nmpts * pd
if va < 0:
va = 0
d = integrate.quad(lambda z: ln_distribution(z, radius, pd), va, vb)
else:
d = 1
def fun(z, x, radius, pd):
if Dist == 1:
return hc_formfactor(x, z) * g_distribution(z, radius, pd)
elif Dist == 2:
return hc_formfactor(x, z) * ln_distribution(z, radius, pd)
else:
return hc_formfactor(x, z)
y = integrate.quad(lambda z: fun(z, x, radius, pd), va, vb)[0]
return y/d[0]
if __name__ == '__main__':
plot_hc(radius=40, pd=0.5)
As suggested in the documentation, I should use for loop, but that reduced the speed even more. The code using for loop is as follows:
from numpy import vectorize
from scipy import integrate
from scipy.special import j1
from math import sqrt, exp, pi, log
import matplotlib.pyplot as plt
import numpy as np
from numpy import empty
def plot_hc(radius, pd):
q = np.linspace(0.008, 1.0, num=500)
y = hc(q, radius, pd)
plt.loglog(q, y)
plt.show()
def hc_formfactor(q, radius):
y = (1.0 / q) * (radius * j1(q * radius))
y = y ** 2
return y
def g_distribution(z, radius, pd):
return (1 / (sqrt(2 * pi) * pd)) * exp(
-((z - radius) / (sqrt(
2) * pd)) ** 2)
def ln_distribution(z, radius, pd):
return (1 / (sqrt(2 * pi) * pd * z / radius)) * exp(
-(log(z / radius) / (sqrt(2) * pd)) ** 2)
# Dist=1(for G_Distribution)
# Dist=2(for LN Distribution)
Dist = 1
def hc(q, radius, pd):
if Dist == 1:
nmpts = 4
va = radius - nmpts * pd
vb = radius + nmpts * pd
if va < 0:
va = 0
d = integrate.quad(lambda z: g_distribution(z, radius, pd), va,vb)
elif Dist == 2:
nmpts = 4
va = radius - nmpts * pd
vb = radius + nmpts * pd
if va < 0:
va = 0
d = integrate.quad(lambda z: ln_distribution(z, radius, pd), va, vb)
else:
d = 1
def fun(z, q, radius, pd):
if Dist == 1:
return hc_formfactor(q, z) * g_distribution(z, radius, pd)
elif Dist == 2:
return hc_formfactor(q, z) * ln_distribution(z, radius, pd)
else:
return hc_formfactor(q, z)
y = empty([len(q)])
for n in range(len(q)):
y[n] = integrate.quad(lambda z: fun(z, q[n], radius, pd), va, vb)[0]
return y / d[0]
if __name__ == '__main__':
plot_hc(radius=40, pd=0.5)
I don't understand, what should I do to optimise the code? If I run the same program for the same values in Matlab it is very fast. I don't know what mistake I did here. Also, some suggested to use jit from numba to speed up the integration, but I am not sure how to implement it. Please help :)
