python - 在cython中声明一个类的实例

Question

我对声明类实例的方式有疑问。我有两个类(class)，即 cosmology 和 NFWHalo。我需要在另一个中使用 cosmology 类中的一些方法。

cdef extern from "gsl/gsl_math.h":
    ctypedef struct gsl_function:
        double (* function) (double x, void * params)
        void * params

cdef extern from "gsl/gsl_integration.h":
    ctypedef struct gsl_integration_workspace
    gsl_integration_workspace *  gsl_integration_workspace_alloc(size_t n)
    void  gsl_integration_workspace_free(gsl_integration_workspace * w)
    int  gsl_integration_qags(const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double *result, double *abserr)
cdef double func_callback(double x, void* params): 
     return (<cosmology>params).__angKernel(x) 

std_G=4.3e-9 # Newton's const   in Mpc (km/s)^2 M_sol^{-1}
v_c = 299792.458 #km/s
cdef class cosmology(object):
    cdef public: 
         double omega_m
         double omega_l
         double omega_r
         double omega_c
         double h
         double w
         double G
         double v_c

    def __init__(self, omega_m = 0.3,  omega_l = 0.7,  h = 0.7, w = -1, omega_r = 0., G = std_G):

        self.omega_m = omega_m
        self.omega_l = omega_l
        self.omega_r = omega_r
        self.omega_c = (1. - omega_m - omega_l)
        self.h = h
        self.w = w
        self.G = G
        self.v_c = v_c
    def __copy__(self):

        return cosmology(omega_m = self.omega_m, omega_l = self.omega_l, h = self.h, w = self.w, omega_r = self.omega_r, G = self.G)

    cpdef double a(self, double z):
        return 1./(1.+z)

    cpdef double E(self, double a):

        return (self.omega_r*a**(-4) + self.omega_m*a**(-3) + self.omega_c*a**(-2) + self.omega_l)**0.5
    @cython.boundscheck(False)
    @cython.wraparound(False)
    @cython.nonecheck(False)
    cdef double __angKernel(self, double x):
         return self.E(x**-1)**-1

    cpdef double Da(self, double z, double z_ref=0):
          if z < 0:
             raise ValueError("Redshift z must not be negative")
          if z < z_ref:
             raise ValueError("Redshift z must not be smaller than the reference redshift")
          cdef gsl_integration_workspace* w =gsl_integration_workspace_alloc(1000)

          cdef gsl_function F
          F.function = &func_callback
          F.params = <void*>self    
          cdef double result = 3, error = 5
          cdef double d, err, rk
          gsl_integration_qags(&F, z_ref+1, z+1, 0, 1e-7, 1000, w, &result, &error)
          d, err = result, error 
          # check for curvature
          rk = (fabs(self.omega_c))**0.5
          if (rk*d > 0.01):
             if self.omega_c > 0:
                d = sinh(rk*d)/rk
             if self.omega_c < 0:
                d = sin(rk*d)/rk    
          gsl_integration_workspace_free(w) 
          return d/(1.+z)

NFWHalo 类

cdef class NFWHalo(cosmology):
    cdef object cosmo
    cdef public double M, c, z
    cdef public double rs, rs_arcsec
    cdef int SIZE
    cdef double[::1] ks,zs, halo_pos

    cdef char* path
    @cython.boundscheck(False)
    @cython.cdivision(True)
    @cython.wraparound(False)
    @cython.nonecheck(False)
    def __cinit__(self, M, c, z, halo_pos, filename=None, zs=None, cosmo = None, omega_m=None, omega_l=None, h=None):
        if omega_m is None: 
           omega_m=0.3
        if omega_l is None:
           omega_l=1-omega_m
        if h is None:                
           h=1.
        if cosmo is None:
           self.cosmo = cosmology(omega_m=omega_m, omega_l=omega_l, h=h)
        else:
           self.cosmo = copy.copy(cosmo)
           print "dimensionless Hubble constant: ", self.cosmo.h
           print "Omega Lambda:", self.cosmo.omega_l
           print "Omega Matter: ", self.cosmo.omega_m

        if filename is None:
           raise ValueError("Could not find a path to the file which contains the table of angular diameter distances")
        self.M = M
        self.c = c
        self.z = z
        self.halo_pos = halo_pos
        self.path=filename

        if zs is None:
           raise ValueError("You must give an array which contains the steps where the redshift probability distribution are computed!")
        self.zs=zs
        # calculate scale radius

        cdef double a = self.cosmo.a(self.z)

        cdef double R200 = 1.63e-5/(1+self.z) * (self.M * self.__omega(a)/self.__omega(1))**0.3333 # in Mpc/h
        self.rs = R200/self.c
        cdef double dl=  self.cosmo.Da(self.z)*3000 #in Mpc/h
        cdef double scale = self.rs / dl
        cdef double arcsec2rad = 1./206265
        self.rs_arcsec = scale/arcsec2rad
        print self.rs_arcsec

    def __repr__(self):
        c = self.__class__.__name__
        return "%s(angular diameter distances info at %s)" % (c, self.path)
    cpdef double __omega(self, double a):
          return self.cosmo.omega_m/(self.cosmo.E(a)**2 * a**3)
    cpdef double __ks(self, double z_s):
          """Lensing strength of halo as function of source redshift.
          """
          cdef double rho_c,Sigma_c, d0, a, ez, rho_s 
          # critical density and surface density
          rho_c = 2.7722e11
          Sigma_c = 5.5444e14
          a  = self.cosmo.a(self.z)
          ez = self.cosmo.E(a)
          d0 = 200./3 * self.c**3/(log(1+self.c) - (1.*self.c)/(1+self.c))
          rho_s = rho_c * ez**2 *d0

          cdef double dl, k_s
          dl  = self.cosmo.Da(z_s, self.z) * self.cosmo.Da(self.z) / self.cosmo.Da(z_s)
          k_s = dl * self.rs * rho_s / Sigma_c
          return k_s

    @cython.cdivision(True)    
    @cython.boundscheck(False)
    @cython.wraparound(False)
    cdef void get_ks(self):
         cdef FILE* handle
         cdef Py_ssize_t i, nz
         nz = len(self.zs)
         # allocate number * sizeof(double) bytes of memory
         cdef double[::1] array = np.empty((nz,))
         if not hasattr(self, 'ks'):  # does self.ks not exist?
            try:
               ## first, check for existing file, see if we can load in self.ks
               handle = fopen(self.path, "r")
               if handle == NULL:
                  raise ValueError("cannot open file {}".format(self.path))

               for i in range(nz-1,-1,-1):
                   fscanf(handle,"%f", &array[i])
               fclose(handle)
               self.ks= array

            except IOError:

               self.ks = self.calculate_ks() 

    cdef double[::1] calculate_ks(self):
         cdef Py_ssize_t i, nz
         nz  = len(self.zs)

         cdef double[::1] k_s = np.empty((nz,))
         for i from nz > i >= 0:         
             k_s[i]= self.__ks(self.zs[i])
         #write the calculated k_s in a file
         cdef FILE* handle=<FILE *>fopen(self.path,"wb")

         fwrite(<void*>&k_s[0], sizeof(double), nz, handle)
         fclose(handle)
         return k_s

我当前的代码已编译，但是当我在 python 中import 并尝试使用它时，它会引发此错误消息:更新

import numpy as np
import WLUtilities
M=4.7779e14
c=3.57
halo_z=0.2577079
halo_pos=np.array([1274.252,1439.9])
path='output.cat'
omega_m=0.3
omega_l=1-omega_m
shear_z=np.array([  0.0e+00,   1.0e-03,   2.0e-03])
nfw = WLUtilities.NFWHalo(M, c, halo_z, halo_pos, path, shear_z, omega_m, omega_l) 

 File "WLUtilities.pyx", line 286, in WLUtilities.NFWHalo.__cinit__ (WLUtilities.c:6406)
    print "dimensionless Hubble constant: ", self.cosmo.h
AttributeError: 'float' object has no attribute 'h'

我不确定我是否正确声明了我的类的实例变量。我想知道如何在 __cinit__ 方法中使用 None，同时在声明实例时它不会引发错误消息？有什么建议吗？？

Answer 1

您正在传递 omega_m，它是一个 float ，作为参数，并在您的函数中将其视为宇宙学对象。这就是引发此错误的原因。

你可能想添加

self.cosmo = cosmology(omega_m=omega_m, omega_l=omega_l, h=h)

或者类似这样的东西，在 if 语句中。

编辑。让我试着更好地解释一下

当您创建 NFWHalo 对象时，您调用

nfw = WLUtilities.NFWHalo(M, c, halo_z, halo_pos, path, shear_z, omega_m, omega_l) 

请注意，您的参数与定义中的参数不匹配，它们是:

 def __cinit__(self, M, c, z, halo_pos, filename=None, zs=None, cosmo = None, omega_m=None, omega_l=None, h=None)

你可能打算改变这一行

nfw = WLUtilities.NFWHalo(M, c, halo_z, halo_pos, path, shear_z, omega_m, omega_l) 

对此

nfw = WLUtilities.NFWHalo(M, c, halo_z, halo_pos, path, shear_z, None, omega_m, omega_l) 

您对许多参数和各自的默认值感到困惑。

[再次]编辑:

为了减少你对此的困惑，而且有很多争论，你为什么不显式地传递你的参数呢？就像这样:

nfw = WLUtilities.NFWHalo(M=M, c=c, z=halo_z, halo_pos = halo_pos, filename = path, zs=shear_z, omega_m = omega_m, omega_l = omega_l)