<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Einasto Profile
In this notebook we visually test the Einasto profile (and do some timing etc.)
Step1:  Density Profile
First, to check for blatant errors, we run with full vectors for both $r$ and $m$
Step2:  Now plot versus $r$
Step3:  Now plot versus $m$
Step4:  Fourier Transform
First plot against $k$
Step5:  Now plot against $m$
Step6:  We may have to be a bit wary of the high-mass, high-k tail, but other than that we should be okay. Finally, to make sure things run properly, try full matrix
Step7:  Perhaps it's better to pre-cache results.
Attempt to Cache Results
<ASSISTANT_TASK:>
Python Code:
%pylab inline
from halomod import HaloModel
from scipy.interpolate import InterpolatedUnivariateSpline as spline
hm = HaloModel(profile_model="Einasto")
Explanation: Einasto Profile
In this notebook we visually test the Einasto profile (and do some timing etc.)
End of explanation
_ = hm.profile.rho(hm.r,hm.m)
Explanation: Density Profile
First, to check for blatant errors, we run with full vectors for both $r$ and $m$:
End of explanation
hm.update(profile_model="Einasto")
plot(hm.r, hm.profile.rho(hm.r,1e12),label="m=12",color="b")
plot(hm.r, hm.profile.rho(hm.r,1e14),label="m=14",color="r")
plot(hm.r, hm.profile.rho(hm.r,1e16),label="m=16",color="g")
hm.update(profile_model="NFW")
plot(hm.r, hm.profile.rho(hm.r,1e12),label="m=12",ls="--",color="b")
plot(hm.r, hm.profile.rho(hm.r,1e14),label="m=14",ls="--",color="r")
plot(hm.r, hm.profile.rho(hm.r,1e16),label="m=16",ls="--",color="g")
legend(loc=0)
xscale('log')
yscale('log')
show()
Explanation: Now plot versus $r$:
End of explanation
hm.update(profile_model="Einasto")
plot(hm.m, hm.profile.rho(0.01,hm.m),label="r=0.01",color="b")
plot(hm.m, hm.profile.rho(0.1,hm.m),label="r=0.1",color="r")
plot(hm.m, hm.profile.rho(1.0,hm.m),label="r=1",color="g")
hm.update(profile_model="NFW")
plot(hm.m, hm.profile.rho(0.01,hm.m),ls="--",color="b")
plot(hm.m, hm.profile.rho(0.1,hm.m),ls="--",color="r")
plot(hm.m, hm.profile.rho(1.0,hm.m),ls="--",color="g")
legend(loc=0)
xscale('log')
yscale('log')
show()
Explanation: Now plot versus $m$:
End of explanation
hm.update(profile_model="Einasto")
#plot(hm.k, hm.profile.u(hm.k,1e12),label="m=12",color="b")
#print hm.profile.u(hm.k,1e12)
#plot(hm.k, hm.profile.u(hm.k,1e14),label="m=14",color="r")
#plot(hm.k, hm.profile.u(hm.k,1e16),label="m=16",color="g")
plot(hm.k, hm.profile.u(hm.k,1e12),label="m=12",color="b")
#print hm.profile.u(hm.k,1e12)
plot(hm.k, hm.profile.u(hm.k,1e14),label="m=14",color="r")
plot(hm.k, hm.profile.u(hm.k,1e16),label="m=16",color="g")
hm.update(profile_model="NFW")
plot(hm.k, hm.profile.u(hm.k,1e12),label="m=12",ls="--",color="b")
plot(hm.k, hm.profile.u(hm.k,1e14),label="m=14",ls="--",color="r")
plot(hm.k, hm.profile.u(hm.k,1e16),label="m=16",ls="--",color="g")
hm.update(profile_model="Einasto")
legend(loc=0)
xscale('log')
yscale('log')
show()
Explanation: Fourier Transform
First plot against $k$:
End of explanation
hm.update(profile_model="Einasto")
plot(hm.m, hm.profile.u(0.01,hm.m),label="k=0.01",color="b")
plot(hm.m, hm.profile.u(5,hm.m),label="k=5",color="r")
plot(hm.m, hm.profile.u(1000,hm.m),label="k=1000",color="g")
hm.update(profile_model="NFW")
plot(hm.m, hm.profile.u(0.01,hm.m),ls="--",color="b")
plot(hm.m, hm.profile.u(5,hm.m),ls="--",color="r")
plot(hm.m, hm.profile.u(1000,hm.m),ls="--",color="g")
legend(loc=0)
xscale('log')
yscale('log')
show()
Explanation: Now plot against $m$:
End of explanation
hm.update(profile_model="Einasto")
%timeit hm.profile.u(hm.k,hm.m)
hm.update(profile_model="NFW")
%timeit hm.profile.u(hm.k,hm.m)
Explanation: We may have to be a bit wary of the high-mass, high-k tail, but other than that we should be okay. Finally, to make sure things run properly, try full matrix:
End of explanation
def f(x,a=0.18):
    return np.exp((-2/a)*(x**a-1))
def _p(K, c):
    minsteps = 1000
    res = np.zeros((len(K),len(c)))
    for ik, kappa in enumerate(K):
        smallest_period = np.pi / kappa
        dx = smallest_period / 8
        nsteps = max(int(np.ceil(c.max() / dx)),minsteps)
        x, dx = np.linspace(0, c.max(), nsteps, retstep=True)
        
        spl = spline(x, x*f(x)*np.sin(kappa*x)/kappa)
        intg = spl.antiderivative()
        res[ik,:] = intg(c) - intg(0)
    
    return np.clip(res,0,None)
K = np.logspace(-4,4,500)
c = np.logspace(0,2,1000)
pk = _p(K,c)
#plot(K,pk)
#xscale('log')
#yscale('log')
np.savez("uKc_einasto.npz",pk=pk,K=K,c=c)
from scipy.interpolate import RectBivariateSpline
def _newp(K,c):
    data = np.load("uKc_einasto.npz")
    pk = data['pk']
    _k = data['K']
    _c = data['c']
    
    c = np.atleast_1d(c)
    if np.isscalar(K):
        K = np.atleast_2d(K)
    if K.ndim < 2:
        if len(K)!=len(c):
            K = np.atleast_2d(K).T # should be len(rs) x len(k)
        else:
            K = np.atleast_2d(K)
    pk[pk<=0] = 1e-8
    
    spl = RectBivariateSpline(np.log(_k),np.log(_c),np.log(pk))
    cc = np.repeat(c,K.shape[0])
    return np.exp(hm.profile._reduce(spl.ev(np.log(K.flatten()),np.log(cc)).reshape(K.shape)))
    
c,K = hm.profile._get_k_variables(hm.k,hm.m)
%timeit _newp(K,c)
plot(np.logspace(-4,4,500),pk[:,0]/ _newp(np.logspace(-4,4,500),1))
plot(np.logspace(-4,4,500),_p(np.logspace(-4,4,500),np.atleast_1d(1)).flatten()/_newp(np.logspace(-4,4,500),1))
plot(np.logspace(-4,4,500),pk[:,0]/ _p(np.logspace(-4,4,500),np.atleast_1d(1)).flatten())
plot()
xscale('log')
Explanation: Perhaps it's better to pre-cache results.
Attempt to Cache Results
End of explanation
<END_TASK>