1

I'm profiling some code, and cProfile reports that almost all the time is spent in <string>:1(<lambda>). What does that mean?
Here's the output of cProfile:

10182 function calls in 191.365 seconds
Ordered by: standard name
ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   727  187.331    0.258  187.331    0.258 <string>:1(<lambda>)
     1    0.000    0.000  191.365  191.365 <string>:1(<module>)
   727    0.791    0.001    0.837    0.001 function_base.py:822(gradient)
   727    0.003    0.000    0.031    0.000 numeric.py:65(zeros_like)
     1    0.000    0.000  191.365  191.365 ode.py:316(integrate)
     1    2.447    2.447  191.365  191.365 ode.py:685(run)
  1454    0.376    0.000    1.213    0.001 useful.py:117(q_powers_opr)
   727    0.296    0.000  188.918    0.260 useful.py:95(dWave)
  1454    0.001    0.000    0.001    0.000 {len}
   727    0.001    0.000    0.001    0.000 {method 'append' of 'list' objects}
   727    0.013    0.000    0.013    0.000 {method 'astype' of 'numpy.ndarray' objects}
     1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
   727    0.024    0.000    0.024    0.000 {method 'fill' of 'numpy.ndarray'objects}
   727    0.004    0.000    0.004    0.000 {numpy.core.multiarray.empty_like} 
   727    0.001    0.000    0.001    0.000 {range}
   727    0.078    0.000    1.291    0.002 {sum}

Here's the profiling code:

import useful
import numpy as np
import sympy as sp
import cProfile
u=sp.Symbol('u')
q=sp.Symbol('q')
x_0=sp.Symbol('x_0')
p_0=sp.Symbol('p_0')
w=sp.Symbol('w')
L_u=sp.Lambda(u,u**4/4-5*u**2/2)
V_q=sp.Lambda(q,10*q)
sys=useful.system(L_u)
pAxis=np.linspace(-60+0j,60,10000)
pWave=useful.numpify(useful.gaussian_wavefunction.subs({x_0:5,p_0:0,w:1}))(pAxis)
pWave_t=sys.time_evolve(pWave,pAxis,V_q)
print "Starting Integration"
cProfile.run("pWave2=pWave_t.integrate(0.2)")

And here's useful.py, which actually does everything:

import numpy as np
import sympy as sp
import functools
from scipy import integrate
from itertools import *
u=sp.Symbol('u')
q=sp.Symbol('q')
p=sp.Symbol('p')
w=sp.Symbol('w')
x_0=sp.Symbol('x_0')
p_0=sp.Symbol('p_0')
hStep=sp.Lambda(p,(1+sp.sign(p))/2)
gaussian_wavefunction=sp.Lambda(p,(2*sp.pi*w**2)**(-sp.Rational(1,4))*\
    sp.exp(-sp.I*x_0*p)*sp.exp(-(p-p_0)**2/(4*w**2)))
class system:
    def __init__(self,K_u):
        """
        This function accepts a definition of kinetic energy in terms of
        velocity (u). Lagrangeans with cross terms between u and q are at the
        moment outside the scope of this program; this may change.
        """
        #This file will use the convention {function}_{input variables}
        #I can't figure out how to allow parameterization of the lagrangean,
        #   so I won't for now.
        self.p_u=sp.Lambda(u,sp.diff(K_u(u),u))
        u_crit=sp.solve(sp.diff(self.p_u(u),u))
        self.p_crit=map(self.p_u,u_crit)
        n_p_crit=sorted([complex(sp.N(p_crit))
            for p_crit in self.p_crit],key=lambda x: x.real)
        #print n_p_crit
        self.H_u=sp.Lambda(u,sp.expand(self.p_u(u)*u-K_u(u)))
        self.u_p=[sp.Lambda(p,u_p_i) for u_p_i in sp.solve(self.p_u(u)-p,u)]
        #This will fail horribly if there is only one critical point;
        #   really you shouldn't be using it in that case, though, so..
        test_points=[2*n_p_crit[0]-n_p_crit[1]] #Left of the leftmost point
        test_points+=[(p_1+p_2)/2 for (p_1,p_2) in
            izip(n_p_crit[:-1],n_p_crit[1:])] #average of every pair of points
        test_points+=[2*n_p_crit[-1]-n_p_crit[-2]]#Right of the rightmost point
        def slopeSign(f):
            #Given a symbolic function f, returns a numerical function over the
            #   same domain, which gives the sign of the slope of f.
            def out(x):
                deriv=sp.diff(f(p),p).subs(p,x)#Note this is still sympified
                return int(np.sign(complex(deriv).real))
            return out
        #Note that the below is a generator. A streight list will not close
        #   arround u_p_i, and they'll all reference the last one. Annoying,
        #   although a generator is better for this job anyway.
        n_orientation_p=(slopeSign(u_p_i) for u_p_i in self.u_p)
        orientation=[[n_orientation_p_i(point)
            for point in test_points] for n_orientation_p_i in n_orientation_p]
        #The following determines the validity of u(p) in each region for each
        #   velocity functuion
        valid=[[-0.0001<complex(u_p_i(point)).imag<0.0001
            for point in test_points] for u_p_i in self.u_p ]
        step=[sp.Lambda(p,hStep(n_p_crit[0]-p))]
        step+=[sp.Lambda(p,hStep(p-p_1)*hStep(p_2-p))
            for (p_1,p_2) in izip(n_p_crit[:-1],n_p_crit[1:])]
        step+=[sp.Lambda(p,hStep(p-n_p_crit[-1]))]
        #The effective factors are functions of p, valued at -1,0,or 1, that
        #   are multiplied by the various possible functions of p based on
        #   functions of u (see H_eff_p for example)
        def signedSteps(sign,nonzero):
            for (step_i,sign_i,nonzero_i) in izip(step,sign,nonzero):
                yield step_i(p)*sign_i*int(nonzero_i)
        def eff_factor(orientation_i,valid_i):
            return sp.Lambda(p,sum(signedSteps(orientation_i,valid_i)))
        self.eff_factors_p=[eff_factor(orientation_i,valid_i)
                for (orientation_i,valid_i) in izip(orientation,valid)]
        self.H_eff_p=self.effective_function(self.H_u)
        #a system will have the following:
        #p_u,H_u,u_p,H_p,eff_factors_p,H_eff_p
    def effective_function(self,f_u):
        f_p=[sp.Lambda(p,sp.expand(f_u(u_p_i(p)))) for u_p_i in self.u_p]
        #The above is a list of functions in terms of p.
        return sp.Lambda(p,sum(f_p_i(p)*eff_factors_p_i(p)\
            for (f_p_i,eff_factors_p_i) in izip(f_p,self.eff_factors_p)))
    def time_evolve(self,pWave,pAxis,V_q=0):
        H_eff_p_n=numpify(self.H_eff_p)
        if V_q==0:
            return pWave*np.exp(-(0+1j)*t*H_eff_p_n(pAxis))
        else:
            try:
                #[-2::-1] gives the coeffs in order q^1,q^2,q^3...
                V_q_coeffs=V_q(q).as_poly().all_coeffs()[-2::-1]
                V_q_coeffs=map(complex,V_q_coeffs)#De-sympifies the coeffs
            except TypeError, te:#Not everything is a polynomial...
                print "Nonpolynomial potentials are not supported"
                raise
            if (np.abs(np.diff(pAxis,2))<0.0001).all():
                spacing=pAxis[1]-pAxis[0]
            else:
                spacing=np.gradient(pAxis)
            print spacing
            def dWave(t,pWave,pAxis):
                dWave=-(0+1j)*H_eff_p_n(pAxis)*pWave
                dWave+=-(0+1j)*sum(q_powers_opr(pWave,V_q_coeffs,spacing))
                """
            def dWave(pWave,H_eff_p_n,V_q_coeffs,pAxis,spacing):
                dWave=np.empty_like(pWave)
                dWave[1::2]=-H_eff_p_n(pAxis)*pWave[::2]
                dWave[::2]=H_eff_p_n(pAxis)*pWave[1::2]
                vWave=sum(q_powers_opr(pWave,V_q_coeffs,spacing))
                dWave[1::2]=-vWave[::2]
                dWave[::2]=vWave[1::2]
                del vWave
                """
                return dWave
            ode=integrate.ode(dWave)
            ode.set_integrator('zvode')
            ode.set_f_params(pAxis)
            ode.set_initial_value(pWave)
            #ode.set_f_params(H_eff_p_n,V_q_coeffs,pAxis,spacing)
            return ode
            #integrate.odeint(dWave,pWave,[0,t],
            #    args=(H_eff_p_n,V_q_coeffs,spacing))
def q_powers_opr(wave,coeffs,spacing):
    #nextWave=np.empty_like(wave)
    for c in coeffs:
        yield c*(0+1j)*np.gradient(wave)/spacing
        """
        nextWave[1::2]=c*np.gradient(wave[::2])/spacing
        nextWave[::2]=-c*np.gradient(wave[1::2])/spacing
        wave=nextWave
        """
def numpify(func):
    try:
        return map(numpify,func)
    except TypeError, te:#If there is only one function
        return sp.lambdify(func.variables,func(*func.variables),"numpy")
def niceIfft(pAxis,pWave):
    pRange=pAxis[-1]-pAxis[0]
    p0=pAxis[0]
    length=np.size(pAxis)
    xRange=2*np.pi*length/pRange
    xAxis=np.linspace(-xRange/2,xRange/2,length)
    iAxis=np.arange(length)
    xWave=(1/np.sqrt(2*np.pi)*pRange*np.fft.ifftshift(np.fft.ifft(pWave))*
        np.exp(2*np.pi*(0+1j)*(-p0)/pRange*iAxis))
    return (xAxis,xWave)
7
  • read the next few lines in the report: almost all time is spent in def dWave(t,pWave,pAxis)
    – jfs
    Commented Jun 13, 2012 at 16:19
  • True, but it doesn't explain what part of dWave is slow. Do I just need to profile dWave independently to get the details of what's happening inside it? Commented Jun 13, 2012 at 16:30
  • After profiling dWave, I get similarly useless results. All the time in dWave is spent in <string>:1(<lambda>). Commented Jun 13, 2012 at 16:41
  • cProfile shows times for functions e.g., there is the result for sum()
    – jfs
    Commented Jun 13, 2012 at 16:43
  • I found where dWave is spending most of its time: It's on H_eff_p_n(pAxis). When I profile just that line, the profiler gives the same <string>:1(<lambda>). I'm not posting this as an answer, because I would still be very interested to find out how I could have found this if dWave was longer (without binary searching through the code to isolate out the slow part). Commented Jun 13, 2012 at 16:49

1 Answer 1

1

I think the problem is your numpify function. It's used to create H_eff_p_n, which is the only function called by dWave that isn't explicitly named in the output. I don't know how lambdify works, but it must, in some way or another, create a lambda function by execing a string. And of course, a lambda is anonymous!

>>> (lambda x: x).func_name
'<lambda>'

I just saw your comment above; a simple workaround is to simply use def. As in...

H_eff_p_n_lambda=numpify(self.H_eff_p)
def H_eff_p_n(*args, **kwargs):
    return H_eff_p_n_lambda(*args, **kwargs)

This won't perform as well, but it will at least make it clear (via cumtime) which lambda is being slow. I thought modifying the __name__ or func_name attribute might help, but it didn't work for me.

Looking more closely at the source of profile (which I assume does roughly the same thing as cProfile here), it seems that the function name that profile uses comes not from func.func_name, but rather from func.func_code.co_name, which is read only! So giving a lambda a name for profiling purposes might be hard.

1
  • Yep, that's right. Turns out I can precompute that value and it cuts down execution time by a factor of 60. Commented Jun 13, 2012 at 17:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.