From b93d6552169d3642ffc5ee7c07bd6aac5b18efcc Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Thu, 2 Jun 2022 20:17:59 -0400 Subject: [PATCH 1/4] Add pyproject.toml to have numpy available during installation --- pyproject.toml | 2 + setup.py | 183 +++++++++++++++++++++++-------------------------- 2 files changed, 88 insertions(+), 97 deletions(-) create mode 100644 pyproject.toml diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..a497ed1 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,2 @@ +[build-system] +requires = ["setuptools","wheel","numpy"] \ No newline at end of file diff --git a/setup.py b/setup.py index 35f9312..4156862 100644 --- a/setup.py +++ b/setup.py @@ -6,26 +6,16 @@ import tempfile import subprocess import shutil +import numpy from setuptools import setup, Extension from setuptools.command.build_ext import build_ext -# can we build recfile? +# To build recfile packages = ["esutil"] ext_modules = [] -try: - import numpy - - include_dirs = [numpy.get_include()] - include_dirs += ["esutil/include"] - have_numpy = True -except ImportError: - have_numpy = False - ext_modules = [] - include_dirs = [] - - stdout.write("Numpy not found: Not building C extensions\n") - time.sleep(5) +include_dirs = [numpy.get_include()] +include_dirs += ["esutil/include"] extra_compile_args = [] extra_link_args = [] @@ -172,88 +162,87 @@ def build_extensions(self): # Make the extensions to be built # -if have_numpy: - # recfile - include_dirs += ["esutil/recfile"] - recfile_sources = [ - "esutil/recfile/records.cpp", - "esutil/recfile/records_wrap.cpp", - ] - recfile_module = Extension( - "esutil.recfile._records", - extra_compile_args=extra_compile_args, - extra_link_args=extra_link_args, - sources=recfile_sources, - include_dirs=include_dirs, - ) - ext_modules.append(recfile_module) - packages.append("esutil.recfile") - - # cosmology package - cosmo_sources = glob("esutil/cosmology/*.c") - cosmo_module = Extension( - "esutil.cosmology._cosmolib", - extra_compile_args=extra_compile_args, - extra_link_args=extra_link_args, - sources=cosmo_sources, - include_dirs=include_dirs, - ) - ext_modules.append(cosmo_module) - packages.append("esutil.cosmology") - - # HTM - include_dirs += ["esutil/htm", "esutil/htm/htm_src"] - htm_sources = glob("esutil/htm/htm_src/*.cpp") - htm_sources += ["esutil/htm/htmc.cc", "esutil/htm/htmc_wrap.cc"] - htm_module = Extension( - "esutil.htm._htmc", - extra_compile_args=extra_compile_args, - extra_link_args=extra_link_args, - sources=htm_sources, - include_dirs=include_dirs, - ) - - ext_modules.append(htm_module) - packages.append("esutil.htm") - - # stat package - # include_dirs += ['esutil/stat'] - # chist_sources = ['chist.cc','chist_wrap.cc'] - # chist_sources = ['esutil/stat/'+s for s in chist_sources] - chist_sources = ["esutil/stat/chist_pywrap.c"] - chist_module = Extension( - "esutil.stat._chist", - extra_compile_args=extra_compile_args, - extra_link_args=extra_link_args, - sources=chist_sources, - include_dirs=include_dirs, - ) - ext_modules.append(chist_module) - stat_util_sources = ["_stat_util.c"] - stat_util_sources = ["esutil/stat/" + s for s in stat_util_sources] - stat_util_module = Extension( - "esutil.stat._stat_util", - extra_compile_args=extra_compile_args, - extra_link_args=extra_link_args, - sources=stat_util_sources, - include_dirs=include_dirs, - ) - ext_modules.append(stat_util_module) - packages.append("esutil.stat") - - # integrate package - - # cgauleg_sources = glob('esutil/integrate/*.cc') - cgauleg_sources = glob("esutil/integrate/cgauleg_pywrap.c") - cgauleg_module = Extension( - "esutil.integrate._cgauleg", - extra_compile_args=extra_compile_args, - extra_link_args=extra_link_args, - sources=cgauleg_sources, - include_dirs=include_dirs, - ) - ext_modules.append(cgauleg_module) - packages.append("esutil.integrate") +# recfile +include_dirs += ["esutil/recfile"] +recfile_sources = [ + "esutil/recfile/records.cpp", + "esutil/recfile/records_wrap.cpp", +] +recfile_module = Extension( + "esutil.recfile._records", + extra_compile_args=extra_compile_args, + extra_link_args=extra_link_args, + sources=recfile_sources, + include_dirs=include_dirs, +) +ext_modules.append(recfile_module) +packages.append("esutil.recfile") + +# cosmology package +cosmo_sources = glob("esutil/cosmology/*.c") +cosmo_module = Extension( + "esutil.cosmology._cosmolib", + extra_compile_args=extra_compile_args, + extra_link_args=extra_link_args, + sources=cosmo_sources, + include_dirs=include_dirs, +) +ext_modules.append(cosmo_module) +packages.append("esutil.cosmology") + +# HTM +include_dirs += ["esutil/htm", "esutil/htm/htm_src"] +htm_sources = glob("esutil/htm/htm_src/*.cpp") +htm_sources += ["esutil/htm/htmc.cc", "esutil/htm/htmc_wrap.cc"] +htm_module = Extension( + "esutil.htm._htmc", + extra_compile_args=extra_compile_args, + extra_link_args=extra_link_args, + sources=htm_sources, + include_dirs=include_dirs, +) + +ext_modules.append(htm_module) +packages.append("esutil.htm") + +# stat package +# include_dirs += ['esutil/stat'] +# chist_sources = ['chist.cc','chist_wrap.cc'] +# chist_sources = ['esutil/stat/'+s for s in chist_sources] +chist_sources = ["esutil/stat/chist_pywrap.c"] +chist_module = Extension( + "esutil.stat._chist", + extra_compile_args=extra_compile_args, + extra_link_args=extra_link_args, + sources=chist_sources, + include_dirs=include_dirs, +) +ext_modules.append(chist_module) +stat_util_sources = ["_stat_util.c"] +stat_util_sources = ["esutil/stat/" + s for s in stat_util_sources] +stat_util_module = Extension( + "esutil.stat._stat_util", + extra_compile_args=extra_compile_args, + extra_link_args=extra_link_args, + sources=stat_util_sources, + include_dirs=include_dirs, +) +ext_modules.append(stat_util_module) +packages.append("esutil.stat") + +# integrate package + +# cgauleg_sources = glob('esutil/integrate/*.cc') +cgauleg_sources = glob("esutil/integrate/cgauleg_pywrap.c") +cgauleg_module = Extension( + "esutil.integrate._cgauleg", + extra_compile_args=extra_compile_args, + extra_link_args=extra_link_args, + sources=cgauleg_sources, + include_dirs=include_dirs, +) +ext_modules.append(cgauleg_module) +packages.append("esutil.integrate") long_description = """ A python package including a wide variety of utilities, focused primarily on @@ -291,7 +280,7 @@ def build_extensions(self): packages=packages, cmdclass={"build_ext": MyBuilder}, ext_modules=ext_modules, - install_requires=['numpy'], + install_requires=['numpy','scipy'], ) # If we get to here, then all was fine. Go ahead and delete the files in the From 1550b028b56757463ca0b85a2b3f18bacbf78236 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Thu, 2 Jun 2022 20:18:23 -0400 Subject: [PATCH 2/4] Simplify workflow: don't use conda, just let pip install dependencies --- .github/workflows/test.yaml | 24 ++++-------------------- 1 file changed, 4 insertions(+), 20 deletions(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 0638f02..0df652a 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -2,8 +2,6 @@ name: tests on: push: - branches: - - master pull_request: null jobs: @@ -18,29 +16,15 @@ jobs: steps: - uses: actions/checkout@v2 - - uses: conda-incubator/setup-miniconda@v2 - with: - python-version: ${{ matrix.pyver }} - channels: conda-forge - channel-priority: strict - show-channel-urls: true - - - name: configure conda and install code - shell: bash -l {0} - run: | - conda config --set always_yes yes - conda install -q mamba - - mamba install -q flake8 pytest scipy - - pip install --no-deps -e . + - name: Install code + run: pip install -e . - name: lint - shell: bash -l {0} run: | + pip install flake8 flake8 esutil - name: test - shell: bash -l {0} run: | + pip install pytest pytest -vv esutil From 140cf42586946563124a738e0024fca431394c30 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Mon, 9 Jan 2023 10:09:34 -0500 Subject: [PATCH 3/4] Setup python in tests --- .github/workflows/test.yaml | 5 + esutil/cosmology/cosmolib.cpp | 218 +++++++ esutil/cosmology/cosmolib_pywrap.cpp | 831 +++++++++++++++++++++++++++ esutil/integrate/cgauleg_pywrap.cpp | 141 +++++ esutil/stat/_stat_util.cpp | 122 ++++ esutil/stat/chist_pywrap.cpp | 148 +++++ tmp/tmpd_h00cl2.cpp | 7 + tmp/tmpv19mw1dh.exe | Bin 0 -> 35036 bytes tmp/tmpvz1p6dh5.os | Bin 0 -> 3088 bytes 9 files changed, 1472 insertions(+) create mode 100644 esutil/cosmology/cosmolib.cpp create mode 100644 esutil/cosmology/cosmolib_pywrap.cpp create mode 100644 esutil/integrate/cgauleg_pywrap.cpp create mode 100644 esutil/stat/_stat_util.cpp create mode 100644 esutil/stat/chist_pywrap.cpp create mode 100644 tmp/tmpd_h00cl2.cpp create mode 100755 tmp/tmpv19mw1dh.exe create mode 100644 tmp/tmpvz1p6dh5.os diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 0df652a..7a4a34e 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -16,6 +16,11 @@ jobs: steps: - uses: actions/checkout@v2 + - name: Set up Python ${{ matrix.pyver }} + uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.pyver }} + - name: Install code run: pip install -e . diff --git a/esutil/cosmology/cosmolib.cpp b/esutil/cosmology/cosmolib.cpp new file mode 100644 index 0000000..5432bf4 --- /dev/null +++ b/esutil/cosmology/cosmolib.cpp @@ -0,0 +1,218 @@ +#include +#include +#include "cosmolib.h" + + +struct cosmo* cosmo_new( + double DH, + int flat, + double omega_m, + double omega_l, + double omega_k) { + + struct cosmo* c; + c=(struct cosmo* ) calloc(1,sizeof(struct cosmo)); + + if (c == NULL) { + return NULL; + } + + c->DH=DH; + c->flat=flat; + c->omega_m=omega_m; + c->omega_l=omega_l; + c->omega_k=omega_k; + + c->tcfac = 0; + if (c->flat != 1) { + if (c->omega_k > 0) { + c->tcfac = sqrt(c->omega_k)/c->DH; + } else { + c->tcfac = sqrt(-c->omega_k)/c->DH; + } + } + + gauleg(-1.0,1.0, NPTS, c->x, c->w); + gauleg(-1.0,1.0, VNPTS, c->vx, c->vw); + + return c; +} + + +/* comoving distance in Mpc */ +double Dc(struct cosmo* c, double zmin, double zmax) { + return c->DH*ez_inverse_integral(c, zmin, zmax); +} + + +// transverse comoving distance +double Dm(struct cosmo* c, double zmin, double zmax) { + + double d; + + d = Dc(c, zmin, zmax); + + if (!c->flat) { + if (c->omega_k > 0) { + d= sinh(d*c->tcfac)/c->tcfac; + } else { + d= sin(d*c->tcfac)/c->tcfac; + } + } + return d; +} + + + +// angular diameter distances +double Da(struct cosmo* c, double zmin, double zmax) { + double d; + d = Dm(c, zmin, zmax); + d /= (1.+zmax); + return d; +} + + + + +// luminosity distances +double Dl(struct cosmo* c, double zmin, double zmax) { + double d; + d = Dm(c, zmin, zmax); + d *= (1.+zmax); + return d; +} + +// comoving volume element +double dV(struct cosmo* c, double z) { + double da, ezinv, oneplusz; + double dv; + + oneplusz = 1.+z; + + da = Da(c, 0.0, z); + ezinv = ez_inverse(c, z); + dv = c->DH*da*da*ezinv*oneplusz*oneplusz; + + return dv; +} + +// comoving volume between zmin and zmax +double V(struct cosmo* c, double zmin, double zmax) { + int i; + double f1,f2,z; + double dv; + double v=0; + + f1 = (zmax-zmin)/2.; + f2 = (zmax+zmin)/2.; + + for (i=0; ivx[i]*f1 + f2; + dv = dV(c, z); + v += f1*dv*c->vw[i]; + } + + return v*4.*M_PI; + +} + + +// inverse critical density for lensing +double scinv(struct cosmo* c, double zl, double zs) { + double dl, ds, dls; + + if (zs <= zl) { + return 0.0; + } + + dl = Da(c, 0.0, zl); + ds = Da(c, 0.0, zs); + dls = Da(c, zl, zs); + return dls*dl/ds*FOUR_PI_G_OVER_C_SQUARED; +} + + + +double ez_inverse(struct cosmo* c, double z) { + double oneplusz, oneplusz2; + double ezi; + + oneplusz = 1.+z; + if (c->flat) { + ezi = c->omega_m*oneplusz*oneplusz*oneplusz + c->omega_l; + } else { + oneplusz2 = oneplusz*oneplusz; + ezi = c->omega_m*oneplusz2*oneplusz + c->omega_k*oneplusz2 + c->omega_l; + } + ezi = sqrt(1.0/ezi); + return ezi; +} + + +double ez_inverse_integral(struct cosmo* c, double zmin, double zmax) { + int i; + double f1, f2, z, ezinv_int=0, ezinv; + + f1 = (zmax-zmin)/2.; + f2 = (zmax+zmin)/2.; + + ezinv_int = 0.0; + + for (i=0;ix[i]*f1 + f2; + + ezinv = ez_inverse(c, z); + ezinv_int += f1*ezinv*c->w[i]; + } + + return ezinv_int; + +} + +void gauleg(double x1, double x2, int npts, double* x, double* w) { + int i, j, m; + double xm, xl, z1, z, p1, p2, p3, pp=0, EPS, abszdiff; + + EPS = 4.e-11; + + m = (npts + 1)/2; + + xm = (x1 + x2)/2.0; + xl = (x2 - x1)/2.0; + z1 = 0.0; + + for (i=1; i<= m; ++i) + { + + z=cos( M_PI*(i-0.25)/(npts+.5) ); + + abszdiff = fabs(z-z1); + + while (abszdiff > EPS) + { + p1 = 1.0; + p2 = 0.0; + for (j=1; j <= npts;++j) + { + p3 = p2; + p2 = p1; + p1 = ( (2.0*j - 1.0)*z*p2 - (j-1.0)*p3 )/j; + } + pp = npts*(z*p1 - p2)/(z*z -1.); + z1=z; + z=z1 - p1/pp; + + abszdiff = fabs(z-z1); + + } + + x[i-1] = xm - xl*z; + x[npts+1-i-1] = xm + xl*z; + w[i-1] = 2.0*xl/( (1.-z*z)*pp*pp ); + w[npts+1-i-1] = w[i-1]; + + + } + +} diff --git a/esutil/cosmology/cosmolib_pywrap.cpp b/esutil/cosmology/cosmolib_pywrap.cpp new file mode 100644 index 0000000..1311297 --- /dev/null +++ b/esutil/cosmology/cosmolib_pywrap.cpp @@ -0,0 +1,831 @@ +/* + + This is a python class definition, wrapping the cosmological distance + calculations in cosmolib.c. The "struct cosmo" is the underlying + "class" and the functions are the methods. + + These wrappers are minimal. Scalars are converted as needed, but there is no + conversion of the input types to arrays for the "vec" vectorized versions of + the functions. It is the responsibility of the python wrapper in cosmology.py + to take care of that. + + I think this is the right compromise: it is messier to write the C versions + of these type checks, but is fairly trivial to write the python checks and + conversions. + + I also could have generated this with SWIG. But I'm experienced enough in + creating these classes that it is actually less work to write it explicitly + than mess around with SWIG complications. And this file is a factor of + ten smaller than the corresponding SWIG wrapper. The size of the SWIG wrapper + is dominated by all the type conversions, which I do in the python wrapper. + + April 2011 + Erin Sheldon, Brookhaven National Laboratory + + */ + +#include +#include "cosmolib.h" +#include + +struct PyCosmoObject { + PyObject_HEAD + struct cosmo* cosmo; +}; + + + +static void +PyCosmoObject_dealloc(struct PyCosmoObject* self) +{ + free(self->cosmo); + +#if PY_MAJOR_VERSION >= 3 + // introduced in python 2.6 + Py_TYPE(self)->tp_free((PyObject*)self); +#else + // old way, removed in python 3 + self->ob_type->tp_free((PyObject*)self); +#endif + + +} + + +static int +PyCosmoObject_init(struct PyCosmoObject* self, PyObject *args, PyObject *kwds) +{ + double DH; + int flat; + double omega_m, omega_l, omega_k; + + free(self->cosmo); + + if (!PyArg_ParseTuple(args, + (char*)"diddd", + &DH, &flat, &omega_m, &omega_l, &omega_k)) { + printf("failed to Parse init"); + return -1; + } + + self->cosmo = cosmo_new(DH, flat, omega_m, omega_l, omega_k); + if (self->cosmo == NULL) { + PyErr_SetString(PyExc_MemoryError, "Failed to allocate struct cosmo"); + return -1; + } + return 0; +} + +static PyObject * +PyCosmoObject_repr(struct PyCosmoObject* self) { +#if PY_MAJOR_VERSION >= 3 + const char* code="y"; +#else + const char* code="s"; +#endif + + char repr[255]; + if (self->cosmo != NULL) { + sprintf(repr, "flat: %d\n" + "DH: %f\n" + "omega_m: %f\n" + "omega_l: %f\n" + "omega_k: %f", + self->cosmo->flat, + self->cosmo->DH, + self->cosmo->omega_m, + self->cosmo->omega_l, + self->cosmo->omega_k); + return Py_BuildValue(code, repr); + } else { + return Py_BuildValue(code, ""); + } +} + +static PyObject* PyCosmoObject_DH(struct PyCosmoObject* self) { + return PyFloat_FromDouble(self->cosmo->DH); +} +static PyObject* PyCosmoObject_flat(struct PyCosmoObject* self) { + return PyLong_FromLong((long) self->cosmo->flat); +} +static PyObject* PyCosmoObject_omega_m(struct PyCosmoObject* self) { + return PyFloat_FromDouble(self->cosmo->omega_m); +} +static PyObject* PyCosmoObject_omega_l(struct PyCosmoObject* self) { + return PyFloat_FromDouble(self->cosmo->omega_l); +} +static PyObject* PyCosmoObject_omega_k(struct PyCosmoObject* self) { + return PyFloat_FromDouble(self->cosmo->omega_k); +} + +/* + The wrapper methods and vectorizations. + + For the array inputs, the caller is responsible for making sure the input is + an array, contiguous, of the right data type. That is much more easily + done in the python wrapper. +*/ + + +static PyObject* +PyCosmoObject_ez_inverse(struct PyCosmoObject* self, PyObject* args) { + double z; + double ezinv; + + if (!PyArg_ParseTuple(args, (char*)"d", &z)) { + return NULL; + } + + ezinv = ez_inverse(self->cosmo, z); + return PyFloat_FromDouble(ezinv); +} +static PyObject* +PyCosmoObject_ez_inverse_vec(struct PyCosmoObject* self, PyObject* args) { + PyObject* zObj=NULL, *resObj=NULL;; + double *z, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"O", &zObj)) { + return NULL; + } + + n = PyArray_SIZE(zObj); + z = (double* )PyArray_DATA(zObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, z[i]); + } + + return resObj; + +} + + + +static PyObject* +PyCosmoObject_ez_inverse_integral(struct PyCosmoObject* self, PyObject* args) { + double zmin, zmax; + double ezinv_int; + + if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { + return NULL; + } + + ezinv_int = ez_inverse_integral(self->cosmo, zmin, zmax); + return PyFloat_FromDouble(ezinv_int); +} + + + + +// comoving distance and vectorizations +static PyObject* +PyCosmoObject_Dc(struct PyCosmoObject* self, PyObject* args) { + double zmin, zmax; + double d; + + if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { + return NULL; + } + + d = Dc(self->cosmo, zmin, zmax); + return PyFloat_FromDouble(d); + +} + + +static PyObject* +PyCosmoObject_Dc_vec1(struct PyCosmoObject* self, PyObject* args) { + PyObject* zminObj=NULL, *resObj=NULL;; + double *zmin, zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"Od", &zminObj, &zmax)) { + return NULL; + } + + n = PyArray_SIZE(zminObj); + zmin = (double* )PyArray_DATA(zminObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo->DH*ez_inverse_integral(self->cosmo, zmin[i], zmax); + } + + return resObj; + +} + +static PyObject* +PyCosmoObject_Dc_vec2(struct PyCosmoObject* self, PyObject* args) { + PyObject* zmaxObj=NULL, *resObj=NULL;; + double zmin, *zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"dO", &zmin, &zmaxObj)) { + return NULL; + } + + n = PyArray_SIZE(zmaxObj); + zmax = (double* )PyArray_DATA(zmaxObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo->DH*ez_inverse_integral(self->cosmo, zmin, zmax[i]); + } + + return resObj; +} + +static PyObject* +PyCosmoObject_Dc_2vec(struct PyCosmoObject* self, PyObject* args) { + PyObject* zmaxObj, *zminObj=NULL, *resObj=NULL; + double *zmin, *zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"OO", &zminObj, &zmaxObj)) { + return NULL; + } + + n = PyArray_SIZE(zminObj); + zmin = (double* )PyArray_DATA(zminObj); + zmax = (double* )PyArray_DATA(zmaxObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo->DH*ez_inverse_integral(self->cosmo, zmin[i], zmax[i]); + } + + return resObj; +} + +// transverse comoving distance and vectorizations +static PyObject* +PyCosmoObject_Dm(struct PyCosmoObject* self, PyObject* args) { + double zmin, zmax; + double d; + + if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { + return NULL; + } + + d = Dm(self->cosmo, zmin, zmax); + return PyFloat_FromDouble(d); + +} + +static PyObject* +PyCosmoObject_Dm_vec1(struct PyCosmoObject* self, PyObject* args) { + PyObject* zminObj=NULL, *resObj=NULL;; + double *zmin, zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"Od", &zminObj, &zmax)) { + return NULL; + } + + n = PyArray_SIZE(zminObj); + zmin = (double* )PyArray_DATA(zminObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zmin[i], zmax); + } + + return resObj; + +} + +static PyObject* +PyCosmoObject_Dm_vec2(struct PyCosmoObject* self, PyObject* args) { + PyObject* zmaxObj=NULL, *resObj=NULL;; + double zmin, *zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"dO", &zmin, &zmaxObj)) { + return NULL; + } + + n = PyArray_SIZE(zmaxObj); + zmax = (double* )PyArray_DATA(zmaxObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zmin, zmax[i]); + } + + return resObj; +} + +static PyObject* +PyCosmoObject_Dm_2vec(struct PyCosmoObject* self, PyObject* args) { + PyObject* zmaxObj, *zminObj=NULL, *resObj=NULL; + double *zmin, *zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"OO", &zminObj, &zmaxObj)) { + return NULL; + } + + n = PyArray_SIZE(zminObj); + zmin = (double* )PyArray_DATA(zminObj); + zmax = (double* )PyArray_DATA(zmaxObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zmin[i], zmax[i]); + } + + return resObj; +} + + +// Angular diameter distance +static PyObject* +PyCosmoObject_Da(struct PyCosmoObject* self, PyObject* args) { + double zmin, zmax; + double d; + + if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { + return NULL; + } + + d = Da(self->cosmo, zmin, zmax); + return PyFloat_FromDouble(d); + +} + +static PyObject* +PyCosmoObject_Da_vec1(struct PyCosmoObject* self, PyObject* args) { + PyObject* zminObj=NULL, *resObj=NULL;; + double *zmin, zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"Od", &zminObj, &zmax)) { + return NULL; + } + + n = PyArray_SIZE(zminObj); + zmin = (double* )PyArray_DATA(zminObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zmin[i], zmax); + } + + return resObj; + +} + +static PyObject* +PyCosmoObject_Da_vec2(struct PyCosmoObject* self, PyObject* args) { + PyObject* zmaxObj=NULL, *resObj=NULL;; + double zmin, *zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"dO", &zmin, &zmaxObj)) { + return NULL; + } + + n = PyArray_SIZE(zmaxObj); + zmax = (double* )PyArray_DATA(zmaxObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zmin, zmax[i]); + } + + return resObj; +} + +static PyObject* +PyCosmoObject_Da_2vec(struct PyCosmoObject* self, PyObject* args) { + PyObject* zmaxObj, *zminObj=NULL, *resObj=NULL; + double *zmin, *zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"OO", &zminObj, &zmaxObj)) { + return NULL; + } + + n = PyArray_SIZE(zminObj); + zmin = (double* )PyArray_DATA(zminObj); + zmax = (double* )PyArray_DATA(zmaxObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zmin[i], zmax[i]); + } + + return resObj; +} + + +// luminosity distance +static PyObject* +PyCosmoObject_Dl(struct PyCosmoObject* self, PyObject* args) { + double zmin, zmax; + double d; + + if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { + return NULL; + } + + d = Dl(self->cosmo, zmin, zmax); + return PyFloat_FromDouble(d); + +} + +static PyObject* +PyCosmoObject_Dl_vec1(struct PyCosmoObject* self, PyObject* args) { + PyObject* zminObj=NULL, *resObj=NULL;; + double *zmin, zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"Od", &zminObj, &zmax)) { + return NULL; + } + + n = PyArray_SIZE(zminObj); + zmin = (double* )PyArray_DATA(zminObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zmin[i], zmax); + } + + return resObj; + +} + +static PyObject* +PyCosmoObject_Dl_vec2(struct PyCosmoObject* self, PyObject* args) { + PyObject* zmaxObj=NULL, *resObj=NULL;; + double zmin, *zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"dO", &zmin, &zmaxObj)) { + return NULL; + } + + n = PyArray_SIZE(zmaxObj); + zmax = (double* )PyArray_DATA(zmaxObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zmin, zmax[i]); + } + + return resObj; +} + +static PyObject* +PyCosmoObject_Dl_2vec(struct PyCosmoObject* self, PyObject* args) { + PyObject* zmaxObj, *zminObj=NULL, *resObj=NULL; + double *zmin, *zmax, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"OO", &zminObj, &zmaxObj)) { + return NULL; + } + + n = PyArray_SIZE(zminObj); + zmin = (double* )PyArray_DATA(zminObj); + zmax = (double* )PyArray_DATA(zmaxObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zmin[i], zmax[i]); + } + + return resObj; +} + +// Comoving volume element and vectorization +static PyObject* +PyCosmoObject_dV(struct PyCosmoObject* self, PyObject* args) { + double z; + double dv; + + if (!PyArg_ParseTuple(args, (char*)"d", &z)) { + return NULL; + } + + dv = dV(self->cosmo, z); + return PyFloat_FromDouble(dv); + +} + +static PyObject* +PyCosmoObject_dV_vec(struct PyCosmoObject* self, PyObject* args) { + PyObject* zObj=NULL, *resObj=NULL;; + double *z, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"O", &zObj)) { + return NULL; + } + + n = PyArray_SIZE(zObj); + z = (double* )PyArray_DATA(zObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, z[i]); + } + + return resObj; + +} + +// Comoving volume between zmin and zmax +static PyObject* +PyCosmoObject_V(struct PyCosmoObject* self, PyObject* args) { + double zmin, zmax; + double v; + + if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { + return NULL; + } + + v = V(self->cosmo, zmin, zmax); + return PyFloat_FromDouble(v); + +} + + + +// Inverse critical density +static PyObject* +PyCosmoObject_scinv(struct PyCosmoObject* self, PyObject* args) { + double zl, zs; + double d; + + if (!PyArg_ParseTuple(args, (char*)"dd", &zl, &zs)) { + return NULL; + } + + d = scinv(self->cosmo, zl, zs); + return PyFloat_FromDouble(d); + +} + +static PyObject* +PyCosmoObject_scinv_vec1(struct PyCosmoObject* self, PyObject* args) { + PyObject* zlObj=NULL, *resObj=NULL;; + double *zl, zs, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"Od", &zlObj, &zs)) { + return NULL; + } + + n = PyArray_SIZE(zlObj); + zl = (double* )PyArray_DATA(zlObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zl[i], zs); + } + + return resObj; + +} + +static PyObject* +PyCosmoObject_scinv_vec2(struct PyCosmoObject* self, PyObject* args) { + PyObject* zsObj=NULL, *resObj=NULL;; + double zl, *zs, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"dO", &zl, &zsObj)) { + return NULL; + } + + n = PyArray_SIZE(zsObj); + zs = (double* )PyArray_DATA(zsObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zl, zs[i]); + } + + return resObj; +} + +static PyObject* +PyCosmoObject_scinv_2vec(struct PyCosmoObject* self, PyObject* args) { + PyObject* zsObj, *zlObj=NULL, *resObj=NULL; + double *zl, *zs, *res; + npy_intp n, i; + + if (!PyArg_ParseTuple(args, (char*)"OO", &zlObj, &zsObj)) { + return NULL; + } + + n = PyArray_SIZE(zlObj); + zl = (double* )PyArray_DATA(zlObj); + zs = (double* )PyArray_DATA(zsObj); + + resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); + res = (double* )PyArray_DATA(resObj); + + for (i=0; icosmo, zl[i], zs[i]); + } + + return resObj; +} + + + + +static PyMethodDef PyCosmoObject_methods[] = { + {"DH", (PyCFunction)PyCosmoObject_DH, METH_VARARGS, "DH\n\nGet the Hubble distance"}, + {"flat", (PyCFunction)PyCosmoObject_flat, METH_VARARGS, "flat\n\nReturn if universe if flat"}, + {"omega_m", (PyCFunction)PyCosmoObject_omega_m, METH_VARARGS, "omega_m\n\nGet omega matter"}, + {"omega_l", (PyCFunction)PyCosmoObject_omega_l, METH_VARARGS, "omega_m\n\nGet omega lambda"}, + {"omega_k", (PyCFunction)PyCosmoObject_omega_k, METH_VARARGS, "omega_m\n\nGet omega curvature"}, + {"ez_inverse", (PyCFunction)PyCosmoObject_ez_inverse, METH_VARARGS, "ez_inverse(z)\n\nGet 1/E(z)"}, + {"ez_inverse_vec", (PyCFunction)PyCosmoObject_ez_inverse_vec, METH_VARARGS, "ez_inverse_vec(z)\n\nGet 1/E(z) for z an array"}, + {"ez_inverse_integral", (PyCFunction)PyCosmoObject_ez_inverse_integral, METH_VARARGS, "ez_inverse_integral(zmin, zmax)\n\nGet integral of 1/E(z) from zmin to zmax"}, + {"Dc", (PyCFunction)PyCosmoObject_Dc, METH_VARARGS, "Dc(zmin,zmax)\n\nComoving distance between zmin and zmax"}, + {"Dc_vec1", (PyCFunction)PyCosmoObject_Dc_vec1, METH_VARARGS, "Dc_vec1(zmin,zmax)\n\nComoving distance between zmin(array) and zmax"}, + {"Dc_vec2", (PyCFunction)PyCosmoObject_Dc_vec2, METH_VARARGS, "Dc_vec2(zmin,zmax)\n\nComoving distance between zmin and zmax(array)"}, + {"Dc_2vec", (PyCFunction)PyCosmoObject_Dc_2vec, METH_VARARGS, "Dc_2vec(zmin,zmax)\n\nComoving distance between zmin and zmax both arrays"}, + {"Dm", (PyCFunction)PyCosmoObject_Dm, METH_VARARGS, "Dm(zmin,zmax)\n\nTransverse comoving distance between zmin and zmax"}, + {"Dm_vec1", (PyCFunction)PyCosmoObject_Dm_vec1, METH_VARARGS, "Dm_vec1(zmin,zmax)\n\nTransverse Comoving distance between zmin(array) and zmax"}, + {"Dm_vec2", (PyCFunction)PyCosmoObject_Dm_vec2, METH_VARARGS, "Dm_vec2(zmin,zmax)\n\nTransverse Comoving distance between zmin and zmax(array)"}, + {"Dm_2vec", (PyCFunction)PyCosmoObject_Dm_2vec, METH_VARARGS, "Dm_2vec(zmin,zmax)\n\nTransverse Comoving distance between zmin and zmax both arrays"}, + {"Da", (PyCFunction)PyCosmoObject_Da, METH_VARARGS, "Da(zmin,zmax)\n\nAngular diameter distance distance between zmin and zmax"}, + {"Da_vec1", (PyCFunction)PyCosmoObject_Da_vec1, METH_VARARGS, "Da_vec1(zmin,zmax)\n\nAngular diameter distance distance between zmin(array) and zmax"}, + {"Da_vec2", (PyCFunction)PyCosmoObject_Da_vec2, METH_VARARGS, "Da_vec2(zmin,zmax)\n\nAngular diameter distance distance between zmin and zmax(array)"}, + {"Da_2vec", (PyCFunction)PyCosmoObject_Da_2vec, METH_VARARGS, "Da_2vec(zmin,zmax)\n\nAngular diameter distance distance between zmin and zmax both arrays"}, + {"Dl", (PyCFunction)PyCosmoObject_Dl, METH_VARARGS, "Dl(zmin,zmax)\n\nLuminosity distance distance between zmin and zmax"}, + {"Dl_vec1", (PyCFunction)PyCosmoObject_Dl_vec1, METH_VARARGS, "Dl_vec1(zmin,zmax)\n\nLuminosity distance distance between zmin(array) and zmax"}, + {"Dl_vec2", (PyCFunction)PyCosmoObject_Dl_vec2, METH_VARARGS, "Dl_vec2(zmin,zmax)\n\nLuminosity distance distance between zmin and zmax(array)"}, + {"Dl_2vec", (PyCFunction)PyCosmoObject_Dl_2vec, METH_VARARGS, "Dl_2vec(zmin,zmax)\n\nLuminosity distance distance between zmin and zmax both arrays"}, + {"dV", (PyCFunction)PyCosmoObject_dV, METH_VARARGS, "dV(z)\n\nComoving volume element at redshift z"}, + {"dV_vec", (PyCFunction)PyCosmoObject_dV_vec, METH_VARARGS, "dV(z)\n\nComoving volume element at redshift z(array)"}, + {"V", (PyCFunction)PyCosmoObject_V, METH_VARARGS, "V(z)\n\nComoving volume between zmin and zmax"}, + {"scinv", (PyCFunction)PyCosmoObject_scinv, METH_VARARGS, "scinv(zl,zs)\n\nInverse critical density distance between zl and zs"}, + {"scinv_vec1", (PyCFunction)PyCosmoObject_scinv_vec1, METH_VARARGS, "scinv_vec1(zl,zs)\n\nInverse critical density distance between zl(array) and zs"}, + {"scinv_vec2", (PyCFunction)PyCosmoObject_scinv_vec2, METH_VARARGS, "scinv_vec2(zl,zs)\n\nInverse critical density distance between zl and zs(array)"}, + {"scinv_2vec", (PyCFunction)PyCosmoObject_scinv_2vec, METH_VARARGS, "scinv_2vec(zl,zs)\n\nInverse critical density distance between zl and zs both arrays"}, + + {NULL} /* Sentinel */ +}; + + + + + + +static PyTypeObject PyCosmoType = { +#if PY_MAJOR_VERSION >= 3 + PyVarObject_HEAD_INIT(NULL, 0) +#else + PyObject_HEAD_INIT(NULL) + 0, /*ob_size*/ +#endif + "_cosmolib.cosmo", /*tp_name*/ + sizeof(struct PyCosmoObject), /*tp_basicsize*/ + 0, /*tp_itemsize*/ + (destructor)PyCosmoObject_dealloc, /*tp_dealloc*/ + 0, /*tp_print*/ + 0, /*tp_getattr*/ + 0, /*tp_setattr*/ + 0, /*tp_compare*/ + //0, /*tp_repr*/ + (reprfunc)PyCosmoObject_repr, /*tp_repr*/ + 0, /*tp_as_number*/ + 0, /*tp_as_sequence*/ + 0, /*tp_as_mapping*/ + 0, /*tp_hash */ + 0, /*tp_call*/ + 0, /*tp_str*/ + 0, /*tp_getattro*/ + 0, /*tp_setattro*/ + 0, /*tp_as_buffer*/ + Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/ + "Cosmology Class", /* tp_doc */ + 0, /* tp_traverse */ + 0, /* tp_clear */ + 0, /* tp_richcompare */ + 0, /* tp_weaklistoffset */ + 0, /* tp_iter */ + 0, /* tp_iternext */ + PyCosmoObject_methods, /* tp_methods */ + 0, /* tp_members */ + 0, /* tp_getset */ + 0, /* tp_base */ + 0, /* tp_dict */ + 0, /* tp_descr_get */ + 0, /* tp_descr_set */ + 0, /* tp_dictoffset */ + //0, /* tp_init */ + (initproc)PyCosmoObject_init, /* tp_init */ + 0, /* tp_alloc */ + //PyCosmoObject_new, /* tp_new */ + PyType_GenericNew, /* tp_new */ +}; + + +static PyMethodDef cosmotype_methods[] = { + {NULL} /* Sentinel */ +}; + + +#if PY_MAJOR_VERSION >= 3 + static struct PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "_cosmolib", /* m_name */ + "Define cosmo type and methods ", /* m_doc */ + -1, /* m_size */ + cosmotype_methods, /* m_methods */ + NULL, /* m_reload */ + NULL, /* m_traverse */ + NULL, /* m_clear */ + NULL, /* m_free */ + }; +#endif + + +#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */ +#define PyMODINIT_FUNC void +#endif +PyMODINIT_FUNC +#if PY_MAJOR_VERSION >= 3 +PyInit__cosmolib(void) +#else +init_cosmolib(void) +#endif +{ + PyObject* m; + + PyCosmoType.tp_new = PyType_GenericNew; + +#if PY_MAJOR_VERSION >= 3 + if (PyType_Ready(&PyCosmoType) < 0) { + return NULL; + } + m = PyModule_Create(&moduledef); + if (m==NULL) { + return NULL; + } + +#else + + if (PyType_Ready(&PyCosmoType) < 0) + return; + + m = Py_InitModule3("_cosmolib", cosmotype_methods, "Define cosmo type and methods."); + + if (m==NULL) { + return; + } +#endif + + Py_INCREF(&PyCosmoType); + PyModule_AddObject(m, "cosmo", (PyObject *)&PyCosmoType); + + import_array(); + + +#if PY_MAJOR_VERSION >= 3 + return m; +#endif +} diff --git a/esutil/integrate/cgauleg_pywrap.cpp b/esutil/integrate/cgauleg_pywrap.cpp new file mode 100644 index 0000000..f00abf6 --- /dev/null +++ b/esutil/integrate/cgauleg_pywrap.cpp @@ -0,0 +1,141 @@ +#include +#include "numpy/arrayobject.h" + +static PyObject* PyCGauleg_cgauleg(PyObject* self, PyObject* args) { + + double x1=0, x2=0; + long npts_long=0; + + npy_intp npts = 0; + PyObject* xarray = NULL; + PyObject* warray = NULL; + double *x=NULL, *w=NULL; + + int i, j, m; + double xm, xl, z1, z, p1, p2, p3, pp=0, pi, EPS, abszdiff; + + PyObject* output_tuple = NULL; + + if (!PyArg_ParseTuple(args, (char*)"ddl", &x1, &x2, &npts_long)) { + return NULL; + } + + npts = npts_long; + + xarray = PyArray_ZEROS( + 1, + &npts, + NPY_FLOAT64, + 0 + ); + warray = PyArray_ZEROS( + 1, + &npts, + NPY_FLOAT64, + 0 + ); + + x = (double* ) PyArray_DATA(xarray); + w = (double* ) PyArray_DATA(warray); + + EPS = 4.e-11; + pi = 3.141592653589793; + + m = (npts + 1)/2; + + xm = (x1 + x2)/2.0; + xl = (x2 - x1)/2.0; + z1 = 0.0; + + for (i=1; i<= m; ++i) + { + + z=cos( pi*(i-0.25)/(npts+.5) ); + + abszdiff = fabs(z-z1); + + while (abszdiff > EPS) + { + p1 = 1.0; + p2 = 0.0; + for (j=1; j <= npts;++j) + { + p3 = p2; + p2 = p1; + p1 = ( (2.0*j - 1.0)*z*p2 - (j-1.0)*p3 )/j; + } + pp = npts*(z*p1 - p2)/(z*z -1.); + z1=z; + z=z1 - p1/pp; + + abszdiff = fabs(z-z1); + + } + + x[i-1] = xm - xl*z; + x[npts+1-i-1] = xm + xl*z; + w[i-1] = 2.0*xl/( (1.-z*z)*pp*pp ); + w[npts+1-i-1] = w[i-1]; + + + } + + + output_tuple = PyTuple_New(2); + PyTuple_SetItem(output_tuple, 0, xarray); + PyTuple_SetItem(output_tuple, 1, warray); + return output_tuple; +} + +static PyMethodDef cgauleg_methods[] = { + {"cgauleg", (PyCFunction)PyCGauleg_cgauleg, METH_VARARGS, "run gauleg"}, + {NULL} /* Sentinel */ +}; + +#if PY_MAJOR_VERSION >= 3 + static struct PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "_cgauleg", /* m_name */ + "Define c version of gauleg", /* m_doc */ + -1, /* m_size */ + cgauleg_methods, /* m_methods */ + NULL, /* m_reload */ + NULL, /* m_traverse */ + NULL, /* m_clear */ + NULL, /* m_free */ + }; +#endif + +#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */ +#define PyMODINIT_FUNC void +#endif +PyMODINIT_FUNC +#if PY_MAJOR_VERSION >= 3 +PyInit__cgauleg(void) +#else +init_cgauleg(void) +#endif +{ + PyObject* m; + +#if PY_MAJOR_VERSION >= 3 + m = PyModule_Create(&moduledef); + if (m==NULL) { + return NULL; + } + +#else + + m = Py_InitModule3("_cgauleg", cgauleg_methods, "Define c version of gauleg."); + + if (m==NULL) { + return; + } +#endif + + import_array(); + +#if PY_MAJOR_VERSION >= 3 + return m; +#endif +} diff --git a/esutil/stat/_stat_util.cpp b/esutil/stat/_stat_util.cpp new file mode 100644 index 0000000..127affb --- /dev/null +++ b/esutil/stat/_stat_util.cpp @@ -0,0 +1,122 @@ +#include +#include + +static PyObject * +PyStatUtil_random_sample(PyObject *self, PyObject *args) +{ + PyObject* randind_obj=NULL; + npy_intp *randind=NULL; + long int nmax=0, nrand=0, seed=0; + npy_intp dims[1]; + long int i=0, ntoselect=0, ntocheck=0; + double prob=0; + if (!PyArg_ParseTuple(args, (char*)"lll", &nmax, &nrand, &seed)) { + return NULL; + } + + if (nmax <= 0 || nrand <= 0) { + PyErr_Format(PyExc_ValueError,"nmax/nrand must be >= 0, got %ld/%ld", nmax, nrand); + return NULL; + } + if (nrand > nmax) { + PyErr_Format(PyExc_ValueError,"nrand must be <= nmax, got %ld/%ld", nmax, nrand); + return NULL; + } + + + srand48(seed); + + dims[0] = nrand; + randind_obj = PyArray_SimpleNew(1, dims, NPY_INTP); + randind = PyArray_DATA(randind_obj); + + ntoselect=nrand; + ntocheck=nmax; + + ntocheck=nmax; + for (i=0; i= 3 + static struct PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "_stat_util", /* m_name */ + "Defines the some gmix fit methods", /* m_doc */ + -1, /* m_size */ + stat_util_module_methods, /* m_methods */ + NULL, /* m_reload */ + NULL, /* m_traverse */ + NULL, /* m_clear */ + NULL, /* m_free */ + }; +#endif + +#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */ +#define PyMODINIT_FUNC void +#endif +PyMODINIT_FUNC +#if PY_MAJOR_VERSION >= 3 +PyInit__stat_util(void) +#else +init_stat_util(void) +#endif +{ + PyObject* m; + + + //PyGMixEMObjectType.tp_new = PyType_GenericNew; + +#if PY_MAJOR_VERSION >= 3 + /* + if (PyType_Ready(&PyGMixEMObjectType) < 0) { + return NULL; + } + */ + m = PyModule_Create(&moduledef); + if (m==NULL) { + return NULL; + } + +#else + /* + if (PyType_Ready(&PyGMixEMObjectType) < 0) { + return; + } + */ + m = Py_InitModule3("_stat_util", stat_util_module_methods, + "This module gmix fit related routines.\n"); + if (m==NULL) { + return; + } +#endif + + /* + Py_INCREF(&PyGMixEMObjectType); + PyModule_AddObject(m, "GMixEM", (PyObject *)&PyGMixEMObjectType); + */ + + import_array(); +#if PY_MAJOR_VERSION >= 3 + return m; +#endif +} diff --git a/esutil/stat/chist_pywrap.cpp b/esutil/stat/chist_pywrap.cpp new file mode 100644 index 0000000..23c647a --- /dev/null +++ b/esutil/stat/chist_pywrap.cpp @@ -0,0 +1,148 @@ +#include +#include "numpy/arrayobject.h" + +/* + + assume hist and rev are contiguous, which we can + guarantee in the caller +*/ + +static PyObject* PyCHist_chist(PyObject* self, PyObject* args) { + + PyObject* data_pyobj=NULL; + double datamin=0; + PyObject* sort_pyobj=NULL; + double binsize=0; + PyObject* hist_pyobj=NULL; + PyObject* rev_pyobj=NULL; + + npy_int64 *hist=NULL, *rev=NULL; + + int dorev=0; + npy_intp nbin = 0, ndata=0, nrev=0; + npy_int64 + i=0, + binnum_old = 0, + offset = 0, data_index = 0, binnum=0, tbin = 0; + double thisdata=0; + + if (!PyArg_ParseTuple(args, (char*)"OdOdOO", + &data_pyobj, + &datamin, + &sort_pyobj, + &binsize, + &hist_pyobj, + &rev_pyobj)) { + return NULL; + } + + if (rev_pyobj != Py_None) { + dorev=1; + rev = (npy_int64 *) PyArray_DATA(rev_pyobj); + nrev = PyArray_SIZE(rev_pyobj); + } + + ndata = PyArray_SIZE(sort_pyobj); + nbin = PyArray_SIZE(hist_pyobj); + + hist=(npy_int64 *) PyArray_DATA(hist_pyobj); + + // this is my reverse engineering of the IDL reverse + // indices + binnum_old = -1; + + for (i=0; i= 0 && binnum < nbin) { + // Should we upate the reverse indices? + if (dorev && (binnum > binnum_old) ) { + tbin = binnum_old + 1; + while (tbin <= binnum) { + rev[tbin] = offset; + tbin++; + } + } + // Update the histogram + hist[binnum] = hist[binnum] + 1; + binnum_old = binnum; + } + } + + tbin = binnum_old + 1; + while (tbin <= nbin) { + if (dorev) { + rev[tbin] = nrev; + } + tbin++; + } + + Py_RETURN_NONE; + +} + +static PyMethodDef chist_methods[] = { + {"chist", (PyCFunction)PyCHist_chist, METH_VARARGS, "histogrammer"}, + {NULL} /* Sentinel */ +}; + +#if PY_MAJOR_VERSION >= 3 + static struct PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "_chist", /* m_name */ + "Define c version of histogrammer", /* m_doc */ + -1, /* m_size */ + chist_methods, /* m_methods */ + NULL, /* m_reload */ + NULL, /* m_traverse */ + NULL, /* m_clear */ + NULL, /* m_free */ + }; +#endif + +#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */ +#define PyMODINIT_FUNC void +#endif +PyMODINIT_FUNC +#if PY_MAJOR_VERSION >= 3 +PyInit__chist(void) +#else +init_chist(void) +#endif +{ + PyObject* m; + +#if PY_MAJOR_VERSION >= 3 + m = PyModule_Create(&moduledef); + if (m==NULL) { + return NULL; + } + +#else + + m = Py_InitModule3("_chist", chist_methods, "Define c version of histogrammer."); + + if (m==NULL) { + return; + } +#endif + + import_array(); + +#if PY_MAJOR_VERSION >= 3 + return m; +#endif +} diff --git a/tmp/tmpd_h00cl2.cpp b/tmp/tmpd_h00cl2.cpp new file mode 100644 index 0000000..111fe4f --- /dev/null +++ b/tmp/tmpd_h00cl2.cpp @@ -0,0 +1,7 @@ + + #include + int main() { + std::cout<<"Hello World\n"; + return 0; + } + \ No newline at end of file diff --git a/tmp/tmpv19mw1dh.exe b/tmp/tmpv19mw1dh.exe new file mode 100755 index 0000000000000000000000000000000000000000..5b02def6a668c6e93d7e5143cbe936e8c4f21acd GIT binary patch literal 35036 zcmeI5e{56N702&;&jAxs>J~_#Elq&7>nI~}LIYX5#t_UB(E&@+k*{6Tz{({vC8k!tekAsZfj-^k;pcsj2rFE@zJ7Oe;XdN zs$Q{Mke%(%Y5T+6%_L=u-!Jt}zaBND;cL$IcVD&D?rFO{+1dW=;?Yb>f82-#!!346 zuD^us?>4(WV`m)X`j)aiye=3H$Y6L)1dVe2?Y8|zZ3krM`pEdEoZjpAAyO_{(6~U> zEL+-Oy$Tq0UpLyC$z%1mEinelgp7fl##IZE^6|Q* zOT3HfGHWh@LP{Rr8jrPNO;gN83mPA_YXXQ0+iu!ojT=Q`3M2}<2>I%5 zLZneR9{O?!rOf&zzckJ}sPHB!y;X?f8`Ry6GK%xqi!U8+f3EC@-&(oq)Qs=7zg>>9 z5{NZstU|58)UD-zO`U~0bNt9D(K?+EAzw6hEl`Z%DWbGJ9xDw6n@JkF;)dQfXW^Vc z7c4YFx(D<2AgBJNPsH_DytFye(N)?O3NfC(@GCcp%k025#WOn?b60Vco%m;e)C0!)Aj zFaajO1egF5U;<2l2{3{G5`m%hzf2BkqNlh-`(yfk(X_clxJcikB;P=1K;LPeLB)}%7fSVs0g8dTmr>WhNAzt~C+~C?Ch8y5ocYC@#oA`kyEm^` zOiXJIv{y4W`ftj%5qmf6?3f#2e;e|*k&z-|wG%qU(~Wsexc~O>%$Fu(%M}qoGI=d zMC@8^!E0A$D#`R%0e@4$b_#qJrBWonnqt7gYKoWqZs=Mb@DrZ#uy*Una5deGU686b1A$CZjJjV zh(US}6Gh_pw|alXk>aramT1&(G2) zWa)Qi>Gx*oGqUstvh-P5`kXAiJWH?4(!ZFcFSPWLA1qKRA1pYeDo+>eQ6>qCJX!{4%kye#{QkBfMGBTpvVIHLKR|?Ll;+$e1cCFPz zp-9P!NGue{r?%7qXKu~B#Ln+S4))@9yoXK?&MM<<{Nu16g>1CXxs$yM<#aCnBUH2^ zr?bQb!U2EOTd1?~*o^=#ZFF3CC6rtVYoM-#*%G#O4%wBk8n_ZzKV~~R2TkqA+PTxt zv8It+33C$2=^QoLU|{CfId@mWI_vF9U@m4k-aDm)`Q!%^U;<2l2`~XBzyz286JP>N zfC(@GCcp%k025#WOn?b60Vco%m;e)C0!)AjFaajO1egF5U;<2l2`~XBzyz286JP>N zfC(@GCcp%k025#WOn?b60Vco%m;e)C0!)AjFoFML0wtJxEo2vje&c-f(gNt#Jpee# z51kVI{-a zj|&ODfk;w5Vg1N`y&i9G(_w#o>5>McLP~dKDB|~p^s4r_F4y?{y5X%~+E6CTg8^@L zv-!bDTsHgSdS#1lpo_P|thIc&D^L^k%Shaa>Ato)|I!BOcKcg>F=@nnK_g!0_j*9`feApk-eKGGdG@k~e%Ftu@1$dvKn|uPHCOsD2W>lQ-)r#;*b}fc=w1s-zqM_; z-+~fMaXWm{ofwqpJ_?w$B+W$Y6~?e_3hI_^<8GDf?0q?Qe~x`2$5yOhqyCF>?BzLj zD97HIV{gx~cjegobL>}h?6)i%Pb<-V9duU*aiKdp=zb2mi-YdrpgTC|yQXjcPRK;a zBuF9TM*Zkc54z9eR{ba2K9GqLEDKr({Ap*lu>a*xn>J)?pDLd#WwhOpw0rve*iNu5 zt!>xCe(Y&C+*WD7mbp7{qiigbUT?!ZxqL~(TuHl6CV;yV0ek^-i^JZDST}P_eyi*~ zy}5mj3}Q4sBND62c(1g)H_Vs#qhXBhie3aYx3`cl+=9(56q2Jmqxf=7jnq5+deo4H zuQ{ZzPB>dah{I z)AtA8KQ`m+@c?s(&<9v|L(v_u(ro>R%o$+h!yWzWC<0iWgo={@V4%zW~!w BjOPFV literal 0 HcmV?d00001 diff --git a/tmp/tmpvz1p6dh5.os b/tmp/tmpvz1p6dh5.os new file mode 100644 index 0000000000000000000000000000000000000000..94d4a425c0cc1c1d464da36310bdd6c1b90091d7 GIT binary patch literal 3088 zcmb7GU2GIp6uvXFZMRU{Vz++?X$d|MbxF4^rF}3gw5t#ZHr+L#0=LuIVY@oJ+wM-G zq$U)2X#x)|@<5OTAK4WXF&fjv#5VY#;Q=KEYhuFIkXY03Kt$-0)cT$IS*JrF@n-LQ zbI!fzJAbox=GK$@|5Pw0FbuksU^a-jG-F~kVon&VhWONFK{8dHWP@SLWw@v5AteUq zx;mX6QLd4TiPd^XP^%P4}x ztPq_v#e2?`jGm4sW0uCvH)_T8g0=oR70NN(+~7+khGa$OBgy0OWR&xLJ&x^c!s6?< zlr!KXl@~>$x?!xMY=xo$(`*GgCShW1#jM5*!ibQJXkRbn4Z!MhTtK{*beJ@VpCW>fQ0O^)mrHV3lbdd~mv{j0FUdMHUn`m6z%a8#9JrArBpKcoCO(BVzL)CSk{&CGugd z&7Eh{XT;jG-n#L=Es}Rl(%ybvTiC^>*5Ci|d|p~_ z#8`mu-O|o49&Z-1`7#E+M%1sRMR|csSc!@MeKLP+M9O~9ICi13WLq(MB$p?!|iW^XF?S6 zLX@%ZM9+^R`(DJ~_-E0R7lk>x#UpHvl@(zfn}awfRNNEUJ>jysv|H7*RP%vUT8pl+ z_Y%Hl3iOv=H7q@g!n}&H;u&P0kCT?iao2jkbv`u!YZ>$VfH+tD73Nvu_WO!B6iFoH2gl$GgO@TQ~%)f!q%y-d-T$^9O-10;?_A3-m!=Kvp=j z{NDg+{SzSS<-Z8T`K8gH19%vR5yKhnr*pp&NN1oI2o46SeRM_wUxzK7mm1ax4<7^~ z^wM77dYA@**jL(V>02y)v!w?Z>t$lK&ntQ>yy6AnA`Bjbq!MV4o4Bb&t?c*5544VH z*Z-DF{@Rl5`d=s6VgD(~=(7F(xC&oQ!VI(fco+VvAlrQ(x7Mq{LE|=Nl_Nfi* w z{7yNH{Z4sneia+RzhkQHzqxY-iIhTCMWuJ1%^R9Sqxa{#^ z@ore}($#b#o|JX<7%8G!uYvo+dK>3~&Uh-r2jq;}1}UTJxQ`CRbuQyIucVThcvMaE zw5kuLQ^!N5;a--qj+2~MtG!!heF})ymxQ7ojV;AtF%qgR5=}bAI`ZlD)Rv~7Ef|T)|Ek9 zrBdaJ3QbPNAZi-E8Lc*87mMvn40 z=w&^X-eYIcVP+BTiz*6$BbG_^b!R>R;DLM&(Ck zMb)X=Exa`zHMf~h$TmekF@*kv485i&hT0fiTiKi4-Hwl8O>H&wfv5BVIOVa?Q9h)m TGpVGk#q|^XSc|cnjyU!&e8`$D literal 0 HcmV?d00001 From 8ba24de0afa0e15eff98ace38857f2fd3ad6a09d Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Mon, 9 Jan 2023 10:13:12 -0500 Subject: [PATCH 4/4] Remove accidentally included files --- esutil/cosmology/cosmolib.cpp | 218 ------- esutil/cosmology/cosmolib_pywrap.cpp | 831 --------------------------- esutil/integrate/cgauleg_pywrap.cpp | 141 ----- esutil/stat/_stat_util.cpp | 122 ---- esutil/stat/chist_pywrap.cpp | 148 ----- tmp/tmpd_h00cl2.cpp | 7 - tmp/tmpv19mw1dh.exe | Bin 35036 -> 0 bytes tmp/tmpvz1p6dh5.os | Bin 3088 -> 0 bytes 8 files changed, 1467 deletions(-) delete mode 100644 esutil/cosmology/cosmolib.cpp delete mode 100644 esutil/cosmology/cosmolib_pywrap.cpp delete mode 100644 esutil/integrate/cgauleg_pywrap.cpp delete mode 100644 esutil/stat/_stat_util.cpp delete mode 100644 esutil/stat/chist_pywrap.cpp delete mode 100644 tmp/tmpd_h00cl2.cpp delete mode 100755 tmp/tmpv19mw1dh.exe delete mode 100644 tmp/tmpvz1p6dh5.os diff --git a/esutil/cosmology/cosmolib.cpp b/esutil/cosmology/cosmolib.cpp deleted file mode 100644 index 5432bf4..0000000 --- a/esutil/cosmology/cosmolib.cpp +++ /dev/null @@ -1,218 +0,0 @@ -#include -#include -#include "cosmolib.h" - - -struct cosmo* cosmo_new( - double DH, - int flat, - double omega_m, - double omega_l, - double omega_k) { - - struct cosmo* c; - c=(struct cosmo* ) calloc(1,sizeof(struct cosmo)); - - if (c == NULL) { - return NULL; - } - - c->DH=DH; - c->flat=flat; - c->omega_m=omega_m; - c->omega_l=omega_l; - c->omega_k=omega_k; - - c->tcfac = 0; - if (c->flat != 1) { - if (c->omega_k > 0) { - c->tcfac = sqrt(c->omega_k)/c->DH; - } else { - c->tcfac = sqrt(-c->omega_k)/c->DH; - } - } - - gauleg(-1.0,1.0, NPTS, c->x, c->w); - gauleg(-1.0,1.0, VNPTS, c->vx, c->vw); - - return c; -} - - -/* comoving distance in Mpc */ -double Dc(struct cosmo* c, double zmin, double zmax) { - return c->DH*ez_inverse_integral(c, zmin, zmax); -} - - -// transverse comoving distance -double Dm(struct cosmo* c, double zmin, double zmax) { - - double d; - - d = Dc(c, zmin, zmax); - - if (!c->flat) { - if (c->omega_k > 0) { - d= sinh(d*c->tcfac)/c->tcfac; - } else { - d= sin(d*c->tcfac)/c->tcfac; - } - } - return d; -} - - - -// angular diameter distances -double Da(struct cosmo* c, double zmin, double zmax) { - double d; - d = Dm(c, zmin, zmax); - d /= (1.+zmax); - return d; -} - - - - -// luminosity distances -double Dl(struct cosmo* c, double zmin, double zmax) { - double d; - d = Dm(c, zmin, zmax); - d *= (1.+zmax); - return d; -} - -// comoving volume element -double dV(struct cosmo* c, double z) { - double da, ezinv, oneplusz; - double dv; - - oneplusz = 1.+z; - - da = Da(c, 0.0, z); - ezinv = ez_inverse(c, z); - dv = c->DH*da*da*ezinv*oneplusz*oneplusz; - - return dv; -} - -// comoving volume between zmin and zmax -double V(struct cosmo* c, double zmin, double zmax) { - int i; - double f1,f2,z; - double dv; - double v=0; - - f1 = (zmax-zmin)/2.; - f2 = (zmax+zmin)/2.; - - for (i=0; ivx[i]*f1 + f2; - dv = dV(c, z); - v += f1*dv*c->vw[i]; - } - - return v*4.*M_PI; - -} - - -// inverse critical density for lensing -double scinv(struct cosmo* c, double zl, double zs) { - double dl, ds, dls; - - if (zs <= zl) { - return 0.0; - } - - dl = Da(c, 0.0, zl); - ds = Da(c, 0.0, zs); - dls = Da(c, zl, zs); - return dls*dl/ds*FOUR_PI_G_OVER_C_SQUARED; -} - - - -double ez_inverse(struct cosmo* c, double z) { - double oneplusz, oneplusz2; - double ezi; - - oneplusz = 1.+z; - if (c->flat) { - ezi = c->omega_m*oneplusz*oneplusz*oneplusz + c->omega_l; - } else { - oneplusz2 = oneplusz*oneplusz; - ezi = c->omega_m*oneplusz2*oneplusz + c->omega_k*oneplusz2 + c->omega_l; - } - ezi = sqrt(1.0/ezi); - return ezi; -} - - -double ez_inverse_integral(struct cosmo* c, double zmin, double zmax) { - int i; - double f1, f2, z, ezinv_int=0, ezinv; - - f1 = (zmax-zmin)/2.; - f2 = (zmax+zmin)/2.; - - ezinv_int = 0.0; - - for (i=0;ix[i]*f1 + f2; - - ezinv = ez_inverse(c, z); - ezinv_int += f1*ezinv*c->w[i]; - } - - return ezinv_int; - -} - -void gauleg(double x1, double x2, int npts, double* x, double* w) { - int i, j, m; - double xm, xl, z1, z, p1, p2, p3, pp=0, EPS, abszdiff; - - EPS = 4.e-11; - - m = (npts + 1)/2; - - xm = (x1 + x2)/2.0; - xl = (x2 - x1)/2.0; - z1 = 0.0; - - for (i=1; i<= m; ++i) - { - - z=cos( M_PI*(i-0.25)/(npts+.5) ); - - abszdiff = fabs(z-z1); - - while (abszdiff > EPS) - { - p1 = 1.0; - p2 = 0.0; - for (j=1; j <= npts;++j) - { - p3 = p2; - p2 = p1; - p1 = ( (2.0*j - 1.0)*z*p2 - (j-1.0)*p3 )/j; - } - pp = npts*(z*p1 - p2)/(z*z -1.); - z1=z; - z=z1 - p1/pp; - - abszdiff = fabs(z-z1); - - } - - x[i-1] = xm - xl*z; - x[npts+1-i-1] = xm + xl*z; - w[i-1] = 2.0*xl/( (1.-z*z)*pp*pp ); - w[npts+1-i-1] = w[i-1]; - - - } - -} diff --git a/esutil/cosmology/cosmolib_pywrap.cpp b/esutil/cosmology/cosmolib_pywrap.cpp deleted file mode 100644 index 1311297..0000000 --- a/esutil/cosmology/cosmolib_pywrap.cpp +++ /dev/null @@ -1,831 +0,0 @@ -/* - - This is a python class definition, wrapping the cosmological distance - calculations in cosmolib.c. The "struct cosmo" is the underlying - "class" and the functions are the methods. - - These wrappers are minimal. Scalars are converted as needed, but there is no - conversion of the input types to arrays for the "vec" vectorized versions of - the functions. It is the responsibility of the python wrapper in cosmology.py - to take care of that. - - I think this is the right compromise: it is messier to write the C versions - of these type checks, but is fairly trivial to write the python checks and - conversions. - - I also could have generated this with SWIG. But I'm experienced enough in - creating these classes that it is actually less work to write it explicitly - than mess around with SWIG complications. And this file is a factor of - ten smaller than the corresponding SWIG wrapper. The size of the SWIG wrapper - is dominated by all the type conversions, which I do in the python wrapper. - - April 2011 - Erin Sheldon, Brookhaven National Laboratory - - */ - -#include -#include "cosmolib.h" -#include - -struct PyCosmoObject { - PyObject_HEAD - struct cosmo* cosmo; -}; - - - -static void -PyCosmoObject_dealloc(struct PyCosmoObject* self) -{ - free(self->cosmo); - -#if PY_MAJOR_VERSION >= 3 - // introduced in python 2.6 - Py_TYPE(self)->tp_free((PyObject*)self); -#else - // old way, removed in python 3 - self->ob_type->tp_free((PyObject*)self); -#endif - - -} - - -static int -PyCosmoObject_init(struct PyCosmoObject* self, PyObject *args, PyObject *kwds) -{ - double DH; - int flat; - double omega_m, omega_l, omega_k; - - free(self->cosmo); - - if (!PyArg_ParseTuple(args, - (char*)"diddd", - &DH, &flat, &omega_m, &omega_l, &omega_k)) { - printf("failed to Parse init"); - return -1; - } - - self->cosmo = cosmo_new(DH, flat, omega_m, omega_l, omega_k); - if (self->cosmo == NULL) { - PyErr_SetString(PyExc_MemoryError, "Failed to allocate struct cosmo"); - return -1; - } - return 0; -} - -static PyObject * -PyCosmoObject_repr(struct PyCosmoObject* self) { -#if PY_MAJOR_VERSION >= 3 - const char* code="y"; -#else - const char* code="s"; -#endif - - char repr[255]; - if (self->cosmo != NULL) { - sprintf(repr, "flat: %d\n" - "DH: %f\n" - "omega_m: %f\n" - "omega_l: %f\n" - "omega_k: %f", - self->cosmo->flat, - self->cosmo->DH, - self->cosmo->omega_m, - self->cosmo->omega_l, - self->cosmo->omega_k); - return Py_BuildValue(code, repr); - } else { - return Py_BuildValue(code, ""); - } -} - -static PyObject* PyCosmoObject_DH(struct PyCosmoObject* self) { - return PyFloat_FromDouble(self->cosmo->DH); -} -static PyObject* PyCosmoObject_flat(struct PyCosmoObject* self) { - return PyLong_FromLong((long) self->cosmo->flat); -} -static PyObject* PyCosmoObject_omega_m(struct PyCosmoObject* self) { - return PyFloat_FromDouble(self->cosmo->omega_m); -} -static PyObject* PyCosmoObject_omega_l(struct PyCosmoObject* self) { - return PyFloat_FromDouble(self->cosmo->omega_l); -} -static PyObject* PyCosmoObject_omega_k(struct PyCosmoObject* self) { - return PyFloat_FromDouble(self->cosmo->omega_k); -} - -/* - The wrapper methods and vectorizations. - - For the array inputs, the caller is responsible for making sure the input is - an array, contiguous, of the right data type. That is much more easily - done in the python wrapper. -*/ - - -static PyObject* -PyCosmoObject_ez_inverse(struct PyCosmoObject* self, PyObject* args) { - double z; - double ezinv; - - if (!PyArg_ParseTuple(args, (char*)"d", &z)) { - return NULL; - } - - ezinv = ez_inverse(self->cosmo, z); - return PyFloat_FromDouble(ezinv); -} -static PyObject* -PyCosmoObject_ez_inverse_vec(struct PyCosmoObject* self, PyObject* args) { - PyObject* zObj=NULL, *resObj=NULL;; - double *z, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"O", &zObj)) { - return NULL; - } - - n = PyArray_SIZE(zObj); - z = (double* )PyArray_DATA(zObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, z[i]); - } - - return resObj; - -} - - - -static PyObject* -PyCosmoObject_ez_inverse_integral(struct PyCosmoObject* self, PyObject* args) { - double zmin, zmax; - double ezinv_int; - - if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { - return NULL; - } - - ezinv_int = ez_inverse_integral(self->cosmo, zmin, zmax); - return PyFloat_FromDouble(ezinv_int); -} - - - - -// comoving distance and vectorizations -static PyObject* -PyCosmoObject_Dc(struct PyCosmoObject* self, PyObject* args) { - double zmin, zmax; - double d; - - if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { - return NULL; - } - - d = Dc(self->cosmo, zmin, zmax); - return PyFloat_FromDouble(d); - -} - - -static PyObject* -PyCosmoObject_Dc_vec1(struct PyCosmoObject* self, PyObject* args) { - PyObject* zminObj=NULL, *resObj=NULL;; - double *zmin, zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"Od", &zminObj, &zmax)) { - return NULL; - } - - n = PyArray_SIZE(zminObj); - zmin = (double* )PyArray_DATA(zminObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo->DH*ez_inverse_integral(self->cosmo, zmin[i], zmax); - } - - return resObj; - -} - -static PyObject* -PyCosmoObject_Dc_vec2(struct PyCosmoObject* self, PyObject* args) { - PyObject* zmaxObj=NULL, *resObj=NULL;; - double zmin, *zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"dO", &zmin, &zmaxObj)) { - return NULL; - } - - n = PyArray_SIZE(zmaxObj); - zmax = (double* )PyArray_DATA(zmaxObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo->DH*ez_inverse_integral(self->cosmo, zmin, zmax[i]); - } - - return resObj; -} - -static PyObject* -PyCosmoObject_Dc_2vec(struct PyCosmoObject* self, PyObject* args) { - PyObject* zmaxObj, *zminObj=NULL, *resObj=NULL; - double *zmin, *zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"OO", &zminObj, &zmaxObj)) { - return NULL; - } - - n = PyArray_SIZE(zminObj); - zmin = (double* )PyArray_DATA(zminObj); - zmax = (double* )PyArray_DATA(zmaxObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo->DH*ez_inverse_integral(self->cosmo, zmin[i], zmax[i]); - } - - return resObj; -} - -// transverse comoving distance and vectorizations -static PyObject* -PyCosmoObject_Dm(struct PyCosmoObject* self, PyObject* args) { - double zmin, zmax; - double d; - - if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { - return NULL; - } - - d = Dm(self->cosmo, zmin, zmax); - return PyFloat_FromDouble(d); - -} - -static PyObject* -PyCosmoObject_Dm_vec1(struct PyCosmoObject* self, PyObject* args) { - PyObject* zminObj=NULL, *resObj=NULL;; - double *zmin, zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"Od", &zminObj, &zmax)) { - return NULL; - } - - n = PyArray_SIZE(zminObj); - zmin = (double* )PyArray_DATA(zminObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zmin[i], zmax); - } - - return resObj; - -} - -static PyObject* -PyCosmoObject_Dm_vec2(struct PyCosmoObject* self, PyObject* args) { - PyObject* zmaxObj=NULL, *resObj=NULL;; - double zmin, *zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"dO", &zmin, &zmaxObj)) { - return NULL; - } - - n = PyArray_SIZE(zmaxObj); - zmax = (double* )PyArray_DATA(zmaxObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zmin, zmax[i]); - } - - return resObj; -} - -static PyObject* -PyCosmoObject_Dm_2vec(struct PyCosmoObject* self, PyObject* args) { - PyObject* zmaxObj, *zminObj=NULL, *resObj=NULL; - double *zmin, *zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"OO", &zminObj, &zmaxObj)) { - return NULL; - } - - n = PyArray_SIZE(zminObj); - zmin = (double* )PyArray_DATA(zminObj); - zmax = (double* )PyArray_DATA(zmaxObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zmin[i], zmax[i]); - } - - return resObj; -} - - -// Angular diameter distance -static PyObject* -PyCosmoObject_Da(struct PyCosmoObject* self, PyObject* args) { - double zmin, zmax; - double d; - - if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { - return NULL; - } - - d = Da(self->cosmo, zmin, zmax); - return PyFloat_FromDouble(d); - -} - -static PyObject* -PyCosmoObject_Da_vec1(struct PyCosmoObject* self, PyObject* args) { - PyObject* zminObj=NULL, *resObj=NULL;; - double *zmin, zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"Od", &zminObj, &zmax)) { - return NULL; - } - - n = PyArray_SIZE(zminObj); - zmin = (double* )PyArray_DATA(zminObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zmin[i], zmax); - } - - return resObj; - -} - -static PyObject* -PyCosmoObject_Da_vec2(struct PyCosmoObject* self, PyObject* args) { - PyObject* zmaxObj=NULL, *resObj=NULL;; - double zmin, *zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"dO", &zmin, &zmaxObj)) { - return NULL; - } - - n = PyArray_SIZE(zmaxObj); - zmax = (double* )PyArray_DATA(zmaxObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zmin, zmax[i]); - } - - return resObj; -} - -static PyObject* -PyCosmoObject_Da_2vec(struct PyCosmoObject* self, PyObject* args) { - PyObject* zmaxObj, *zminObj=NULL, *resObj=NULL; - double *zmin, *zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"OO", &zminObj, &zmaxObj)) { - return NULL; - } - - n = PyArray_SIZE(zminObj); - zmin = (double* )PyArray_DATA(zminObj); - zmax = (double* )PyArray_DATA(zmaxObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zmin[i], zmax[i]); - } - - return resObj; -} - - -// luminosity distance -static PyObject* -PyCosmoObject_Dl(struct PyCosmoObject* self, PyObject* args) { - double zmin, zmax; - double d; - - if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { - return NULL; - } - - d = Dl(self->cosmo, zmin, zmax); - return PyFloat_FromDouble(d); - -} - -static PyObject* -PyCosmoObject_Dl_vec1(struct PyCosmoObject* self, PyObject* args) { - PyObject* zminObj=NULL, *resObj=NULL;; - double *zmin, zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"Od", &zminObj, &zmax)) { - return NULL; - } - - n = PyArray_SIZE(zminObj); - zmin = (double* )PyArray_DATA(zminObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zmin[i], zmax); - } - - return resObj; - -} - -static PyObject* -PyCosmoObject_Dl_vec2(struct PyCosmoObject* self, PyObject* args) { - PyObject* zmaxObj=NULL, *resObj=NULL;; - double zmin, *zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"dO", &zmin, &zmaxObj)) { - return NULL; - } - - n = PyArray_SIZE(zmaxObj); - zmax = (double* )PyArray_DATA(zmaxObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zmin, zmax[i]); - } - - return resObj; -} - -static PyObject* -PyCosmoObject_Dl_2vec(struct PyCosmoObject* self, PyObject* args) { - PyObject* zmaxObj, *zminObj=NULL, *resObj=NULL; - double *zmin, *zmax, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"OO", &zminObj, &zmaxObj)) { - return NULL; - } - - n = PyArray_SIZE(zminObj); - zmin = (double* )PyArray_DATA(zminObj); - zmax = (double* )PyArray_DATA(zmaxObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zmin[i], zmax[i]); - } - - return resObj; -} - -// Comoving volume element and vectorization -static PyObject* -PyCosmoObject_dV(struct PyCosmoObject* self, PyObject* args) { - double z; - double dv; - - if (!PyArg_ParseTuple(args, (char*)"d", &z)) { - return NULL; - } - - dv = dV(self->cosmo, z); - return PyFloat_FromDouble(dv); - -} - -static PyObject* -PyCosmoObject_dV_vec(struct PyCosmoObject* self, PyObject* args) { - PyObject* zObj=NULL, *resObj=NULL;; - double *z, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"O", &zObj)) { - return NULL; - } - - n = PyArray_SIZE(zObj); - z = (double* )PyArray_DATA(zObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, z[i]); - } - - return resObj; - -} - -// Comoving volume between zmin and zmax -static PyObject* -PyCosmoObject_V(struct PyCosmoObject* self, PyObject* args) { - double zmin, zmax; - double v; - - if (!PyArg_ParseTuple(args, (char*)"dd", &zmin, &zmax)) { - return NULL; - } - - v = V(self->cosmo, zmin, zmax); - return PyFloat_FromDouble(v); - -} - - - -// Inverse critical density -static PyObject* -PyCosmoObject_scinv(struct PyCosmoObject* self, PyObject* args) { - double zl, zs; - double d; - - if (!PyArg_ParseTuple(args, (char*)"dd", &zl, &zs)) { - return NULL; - } - - d = scinv(self->cosmo, zl, zs); - return PyFloat_FromDouble(d); - -} - -static PyObject* -PyCosmoObject_scinv_vec1(struct PyCosmoObject* self, PyObject* args) { - PyObject* zlObj=NULL, *resObj=NULL;; - double *zl, zs, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"Od", &zlObj, &zs)) { - return NULL; - } - - n = PyArray_SIZE(zlObj); - zl = (double* )PyArray_DATA(zlObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zl[i], zs); - } - - return resObj; - -} - -static PyObject* -PyCosmoObject_scinv_vec2(struct PyCosmoObject* self, PyObject* args) { - PyObject* zsObj=NULL, *resObj=NULL;; - double zl, *zs, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"dO", &zl, &zsObj)) { - return NULL; - } - - n = PyArray_SIZE(zsObj); - zs = (double* )PyArray_DATA(zsObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zl, zs[i]); - } - - return resObj; -} - -static PyObject* -PyCosmoObject_scinv_2vec(struct PyCosmoObject* self, PyObject* args) { - PyObject* zsObj, *zlObj=NULL, *resObj=NULL; - double *zl, *zs, *res; - npy_intp n, i; - - if (!PyArg_ParseTuple(args, (char*)"OO", &zlObj, &zsObj)) { - return NULL; - } - - n = PyArray_SIZE(zlObj); - zl = (double* )PyArray_DATA(zlObj); - zs = (double* )PyArray_DATA(zsObj); - - resObj = PyArray_ZEROS(1, &n, NPY_FLOAT64, 0); - res = (double* )PyArray_DATA(resObj); - - for (i=0; icosmo, zl[i], zs[i]); - } - - return resObj; -} - - - - -static PyMethodDef PyCosmoObject_methods[] = { - {"DH", (PyCFunction)PyCosmoObject_DH, METH_VARARGS, "DH\n\nGet the Hubble distance"}, - {"flat", (PyCFunction)PyCosmoObject_flat, METH_VARARGS, "flat\n\nReturn if universe if flat"}, - {"omega_m", (PyCFunction)PyCosmoObject_omega_m, METH_VARARGS, "omega_m\n\nGet omega matter"}, - {"omega_l", (PyCFunction)PyCosmoObject_omega_l, METH_VARARGS, "omega_m\n\nGet omega lambda"}, - {"omega_k", (PyCFunction)PyCosmoObject_omega_k, METH_VARARGS, "omega_m\n\nGet omega curvature"}, - {"ez_inverse", (PyCFunction)PyCosmoObject_ez_inverse, METH_VARARGS, "ez_inverse(z)\n\nGet 1/E(z)"}, - {"ez_inverse_vec", (PyCFunction)PyCosmoObject_ez_inverse_vec, METH_VARARGS, "ez_inverse_vec(z)\n\nGet 1/E(z) for z an array"}, - {"ez_inverse_integral", (PyCFunction)PyCosmoObject_ez_inverse_integral, METH_VARARGS, "ez_inverse_integral(zmin, zmax)\n\nGet integral of 1/E(z) from zmin to zmax"}, - {"Dc", (PyCFunction)PyCosmoObject_Dc, METH_VARARGS, "Dc(zmin,zmax)\n\nComoving distance between zmin and zmax"}, - {"Dc_vec1", (PyCFunction)PyCosmoObject_Dc_vec1, METH_VARARGS, "Dc_vec1(zmin,zmax)\n\nComoving distance between zmin(array) and zmax"}, - {"Dc_vec2", (PyCFunction)PyCosmoObject_Dc_vec2, METH_VARARGS, "Dc_vec2(zmin,zmax)\n\nComoving distance between zmin and zmax(array)"}, - {"Dc_2vec", (PyCFunction)PyCosmoObject_Dc_2vec, METH_VARARGS, "Dc_2vec(zmin,zmax)\n\nComoving distance between zmin and zmax both arrays"}, - {"Dm", (PyCFunction)PyCosmoObject_Dm, METH_VARARGS, "Dm(zmin,zmax)\n\nTransverse comoving distance between zmin and zmax"}, - {"Dm_vec1", (PyCFunction)PyCosmoObject_Dm_vec1, METH_VARARGS, "Dm_vec1(zmin,zmax)\n\nTransverse Comoving distance between zmin(array) and zmax"}, - {"Dm_vec2", (PyCFunction)PyCosmoObject_Dm_vec2, METH_VARARGS, "Dm_vec2(zmin,zmax)\n\nTransverse Comoving distance between zmin and zmax(array)"}, - {"Dm_2vec", (PyCFunction)PyCosmoObject_Dm_2vec, METH_VARARGS, "Dm_2vec(zmin,zmax)\n\nTransverse Comoving distance between zmin and zmax both arrays"}, - {"Da", (PyCFunction)PyCosmoObject_Da, METH_VARARGS, "Da(zmin,zmax)\n\nAngular diameter distance distance between zmin and zmax"}, - {"Da_vec1", (PyCFunction)PyCosmoObject_Da_vec1, METH_VARARGS, "Da_vec1(zmin,zmax)\n\nAngular diameter distance distance between zmin(array) and zmax"}, - {"Da_vec2", (PyCFunction)PyCosmoObject_Da_vec2, METH_VARARGS, "Da_vec2(zmin,zmax)\n\nAngular diameter distance distance between zmin and zmax(array)"}, - {"Da_2vec", (PyCFunction)PyCosmoObject_Da_2vec, METH_VARARGS, "Da_2vec(zmin,zmax)\n\nAngular diameter distance distance between zmin and zmax both arrays"}, - {"Dl", (PyCFunction)PyCosmoObject_Dl, METH_VARARGS, "Dl(zmin,zmax)\n\nLuminosity distance distance between zmin and zmax"}, - {"Dl_vec1", (PyCFunction)PyCosmoObject_Dl_vec1, METH_VARARGS, "Dl_vec1(zmin,zmax)\n\nLuminosity distance distance between zmin(array) and zmax"}, - {"Dl_vec2", (PyCFunction)PyCosmoObject_Dl_vec2, METH_VARARGS, "Dl_vec2(zmin,zmax)\n\nLuminosity distance distance between zmin and zmax(array)"}, - {"Dl_2vec", (PyCFunction)PyCosmoObject_Dl_2vec, METH_VARARGS, "Dl_2vec(zmin,zmax)\n\nLuminosity distance distance between zmin and zmax both arrays"}, - {"dV", (PyCFunction)PyCosmoObject_dV, METH_VARARGS, "dV(z)\n\nComoving volume element at redshift z"}, - {"dV_vec", (PyCFunction)PyCosmoObject_dV_vec, METH_VARARGS, "dV(z)\n\nComoving volume element at redshift z(array)"}, - {"V", (PyCFunction)PyCosmoObject_V, METH_VARARGS, "V(z)\n\nComoving volume between zmin and zmax"}, - {"scinv", (PyCFunction)PyCosmoObject_scinv, METH_VARARGS, "scinv(zl,zs)\n\nInverse critical density distance between zl and zs"}, - {"scinv_vec1", (PyCFunction)PyCosmoObject_scinv_vec1, METH_VARARGS, "scinv_vec1(zl,zs)\n\nInverse critical density distance between zl(array) and zs"}, - {"scinv_vec2", (PyCFunction)PyCosmoObject_scinv_vec2, METH_VARARGS, "scinv_vec2(zl,zs)\n\nInverse critical density distance between zl and zs(array)"}, - {"scinv_2vec", (PyCFunction)PyCosmoObject_scinv_2vec, METH_VARARGS, "scinv_2vec(zl,zs)\n\nInverse critical density distance between zl and zs both arrays"}, - - {NULL} /* Sentinel */ -}; - - - - - - -static PyTypeObject PyCosmoType = { -#if PY_MAJOR_VERSION >= 3 - PyVarObject_HEAD_INIT(NULL, 0) -#else - PyObject_HEAD_INIT(NULL) - 0, /*ob_size*/ -#endif - "_cosmolib.cosmo", /*tp_name*/ - sizeof(struct PyCosmoObject), /*tp_basicsize*/ - 0, /*tp_itemsize*/ - (destructor)PyCosmoObject_dealloc, /*tp_dealloc*/ - 0, /*tp_print*/ - 0, /*tp_getattr*/ - 0, /*tp_setattr*/ - 0, /*tp_compare*/ - //0, /*tp_repr*/ - (reprfunc)PyCosmoObject_repr, /*tp_repr*/ - 0, /*tp_as_number*/ - 0, /*tp_as_sequence*/ - 0, /*tp_as_mapping*/ - 0, /*tp_hash */ - 0, /*tp_call*/ - 0, /*tp_str*/ - 0, /*tp_getattro*/ - 0, /*tp_setattro*/ - 0, /*tp_as_buffer*/ - Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/ - "Cosmology Class", /* tp_doc */ - 0, /* tp_traverse */ - 0, /* tp_clear */ - 0, /* tp_richcompare */ - 0, /* tp_weaklistoffset */ - 0, /* tp_iter */ - 0, /* tp_iternext */ - PyCosmoObject_methods, /* tp_methods */ - 0, /* tp_members */ - 0, /* tp_getset */ - 0, /* tp_base */ - 0, /* tp_dict */ - 0, /* tp_descr_get */ - 0, /* tp_descr_set */ - 0, /* tp_dictoffset */ - //0, /* tp_init */ - (initproc)PyCosmoObject_init, /* tp_init */ - 0, /* tp_alloc */ - //PyCosmoObject_new, /* tp_new */ - PyType_GenericNew, /* tp_new */ -}; - - -static PyMethodDef cosmotype_methods[] = { - {NULL} /* Sentinel */ -}; - - -#if PY_MAJOR_VERSION >= 3 - static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, - "_cosmolib", /* m_name */ - "Define cosmo type and methods ", /* m_doc */ - -1, /* m_size */ - cosmotype_methods, /* m_methods */ - NULL, /* m_reload */ - NULL, /* m_traverse */ - NULL, /* m_clear */ - NULL, /* m_free */ - }; -#endif - - -#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */ -#define PyMODINIT_FUNC void -#endif -PyMODINIT_FUNC -#if PY_MAJOR_VERSION >= 3 -PyInit__cosmolib(void) -#else -init_cosmolib(void) -#endif -{ - PyObject* m; - - PyCosmoType.tp_new = PyType_GenericNew; - -#if PY_MAJOR_VERSION >= 3 - if (PyType_Ready(&PyCosmoType) < 0) { - return NULL; - } - m = PyModule_Create(&moduledef); - if (m==NULL) { - return NULL; - } - -#else - - if (PyType_Ready(&PyCosmoType) < 0) - return; - - m = Py_InitModule3("_cosmolib", cosmotype_methods, "Define cosmo type and methods."); - - if (m==NULL) { - return; - } -#endif - - Py_INCREF(&PyCosmoType); - PyModule_AddObject(m, "cosmo", (PyObject *)&PyCosmoType); - - import_array(); - - -#if PY_MAJOR_VERSION >= 3 - return m; -#endif -} diff --git a/esutil/integrate/cgauleg_pywrap.cpp b/esutil/integrate/cgauleg_pywrap.cpp deleted file mode 100644 index f00abf6..0000000 --- a/esutil/integrate/cgauleg_pywrap.cpp +++ /dev/null @@ -1,141 +0,0 @@ -#include -#include "numpy/arrayobject.h" - -static PyObject* PyCGauleg_cgauleg(PyObject* self, PyObject* args) { - - double x1=0, x2=0; - long npts_long=0; - - npy_intp npts = 0; - PyObject* xarray = NULL; - PyObject* warray = NULL; - double *x=NULL, *w=NULL; - - int i, j, m; - double xm, xl, z1, z, p1, p2, p3, pp=0, pi, EPS, abszdiff; - - PyObject* output_tuple = NULL; - - if (!PyArg_ParseTuple(args, (char*)"ddl", &x1, &x2, &npts_long)) { - return NULL; - } - - npts = npts_long; - - xarray = PyArray_ZEROS( - 1, - &npts, - NPY_FLOAT64, - 0 - ); - warray = PyArray_ZEROS( - 1, - &npts, - NPY_FLOAT64, - 0 - ); - - x = (double* ) PyArray_DATA(xarray); - w = (double* ) PyArray_DATA(warray); - - EPS = 4.e-11; - pi = 3.141592653589793; - - m = (npts + 1)/2; - - xm = (x1 + x2)/2.0; - xl = (x2 - x1)/2.0; - z1 = 0.0; - - for (i=1; i<= m; ++i) - { - - z=cos( pi*(i-0.25)/(npts+.5) ); - - abszdiff = fabs(z-z1); - - while (abszdiff > EPS) - { - p1 = 1.0; - p2 = 0.0; - for (j=1; j <= npts;++j) - { - p3 = p2; - p2 = p1; - p1 = ( (2.0*j - 1.0)*z*p2 - (j-1.0)*p3 )/j; - } - pp = npts*(z*p1 - p2)/(z*z -1.); - z1=z; - z=z1 - p1/pp; - - abszdiff = fabs(z-z1); - - } - - x[i-1] = xm - xl*z; - x[npts+1-i-1] = xm + xl*z; - w[i-1] = 2.0*xl/( (1.-z*z)*pp*pp ); - w[npts+1-i-1] = w[i-1]; - - - } - - - output_tuple = PyTuple_New(2); - PyTuple_SetItem(output_tuple, 0, xarray); - PyTuple_SetItem(output_tuple, 1, warray); - return output_tuple; -} - -static PyMethodDef cgauleg_methods[] = { - {"cgauleg", (PyCFunction)PyCGauleg_cgauleg, METH_VARARGS, "run gauleg"}, - {NULL} /* Sentinel */ -}; - -#if PY_MAJOR_VERSION >= 3 - static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, - "_cgauleg", /* m_name */ - "Define c version of gauleg", /* m_doc */ - -1, /* m_size */ - cgauleg_methods, /* m_methods */ - NULL, /* m_reload */ - NULL, /* m_traverse */ - NULL, /* m_clear */ - NULL, /* m_free */ - }; -#endif - -#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */ -#define PyMODINIT_FUNC void -#endif -PyMODINIT_FUNC -#if PY_MAJOR_VERSION >= 3 -PyInit__cgauleg(void) -#else -init_cgauleg(void) -#endif -{ - PyObject* m; - -#if PY_MAJOR_VERSION >= 3 - m = PyModule_Create(&moduledef); - if (m==NULL) { - return NULL; - } - -#else - - m = Py_InitModule3("_cgauleg", cgauleg_methods, "Define c version of gauleg."); - - if (m==NULL) { - return; - } -#endif - - import_array(); - -#if PY_MAJOR_VERSION >= 3 - return m; -#endif -} diff --git a/esutil/stat/_stat_util.cpp b/esutil/stat/_stat_util.cpp deleted file mode 100644 index 127affb..0000000 --- a/esutil/stat/_stat_util.cpp +++ /dev/null @@ -1,122 +0,0 @@ -#include -#include - -static PyObject * -PyStatUtil_random_sample(PyObject *self, PyObject *args) -{ - PyObject* randind_obj=NULL; - npy_intp *randind=NULL; - long int nmax=0, nrand=0, seed=0; - npy_intp dims[1]; - long int i=0, ntoselect=0, ntocheck=0; - double prob=0; - if (!PyArg_ParseTuple(args, (char*)"lll", &nmax, &nrand, &seed)) { - return NULL; - } - - if (nmax <= 0 || nrand <= 0) { - PyErr_Format(PyExc_ValueError,"nmax/nrand must be >= 0, got %ld/%ld", nmax, nrand); - return NULL; - } - if (nrand > nmax) { - PyErr_Format(PyExc_ValueError,"nrand must be <= nmax, got %ld/%ld", nmax, nrand); - return NULL; - } - - - srand48(seed); - - dims[0] = nrand; - randind_obj = PyArray_SimpleNew(1, dims, NPY_INTP); - randind = PyArray_DATA(randind_obj); - - ntoselect=nrand; - ntocheck=nmax; - - ntocheck=nmax; - for (i=0; i= 3 - static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, - "_stat_util", /* m_name */ - "Defines the some gmix fit methods", /* m_doc */ - -1, /* m_size */ - stat_util_module_methods, /* m_methods */ - NULL, /* m_reload */ - NULL, /* m_traverse */ - NULL, /* m_clear */ - NULL, /* m_free */ - }; -#endif - -#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */ -#define PyMODINIT_FUNC void -#endif -PyMODINIT_FUNC -#if PY_MAJOR_VERSION >= 3 -PyInit__stat_util(void) -#else -init_stat_util(void) -#endif -{ - PyObject* m; - - - //PyGMixEMObjectType.tp_new = PyType_GenericNew; - -#if PY_MAJOR_VERSION >= 3 - /* - if (PyType_Ready(&PyGMixEMObjectType) < 0) { - return NULL; - } - */ - m = PyModule_Create(&moduledef); - if (m==NULL) { - return NULL; - } - -#else - /* - if (PyType_Ready(&PyGMixEMObjectType) < 0) { - return; - } - */ - m = Py_InitModule3("_stat_util", stat_util_module_methods, - "This module gmix fit related routines.\n"); - if (m==NULL) { - return; - } -#endif - - /* - Py_INCREF(&PyGMixEMObjectType); - PyModule_AddObject(m, "GMixEM", (PyObject *)&PyGMixEMObjectType); - */ - - import_array(); -#if PY_MAJOR_VERSION >= 3 - return m; -#endif -} diff --git a/esutil/stat/chist_pywrap.cpp b/esutil/stat/chist_pywrap.cpp deleted file mode 100644 index 23c647a..0000000 --- a/esutil/stat/chist_pywrap.cpp +++ /dev/null @@ -1,148 +0,0 @@ -#include -#include "numpy/arrayobject.h" - -/* - - assume hist and rev are contiguous, which we can - guarantee in the caller -*/ - -static PyObject* PyCHist_chist(PyObject* self, PyObject* args) { - - PyObject* data_pyobj=NULL; - double datamin=0; - PyObject* sort_pyobj=NULL; - double binsize=0; - PyObject* hist_pyobj=NULL; - PyObject* rev_pyobj=NULL; - - npy_int64 *hist=NULL, *rev=NULL; - - int dorev=0; - npy_intp nbin = 0, ndata=0, nrev=0; - npy_int64 - i=0, - binnum_old = 0, - offset = 0, data_index = 0, binnum=0, tbin = 0; - double thisdata=0; - - if (!PyArg_ParseTuple(args, (char*)"OdOdOO", - &data_pyobj, - &datamin, - &sort_pyobj, - &binsize, - &hist_pyobj, - &rev_pyobj)) { - return NULL; - } - - if (rev_pyobj != Py_None) { - dorev=1; - rev = (npy_int64 *) PyArray_DATA(rev_pyobj); - nrev = PyArray_SIZE(rev_pyobj); - } - - ndata = PyArray_SIZE(sort_pyobj); - nbin = PyArray_SIZE(hist_pyobj); - - hist=(npy_int64 *) PyArray_DATA(hist_pyobj); - - // this is my reverse engineering of the IDL reverse - // indices - binnum_old = -1; - - for (i=0; i= 0 && binnum < nbin) { - // Should we upate the reverse indices? - if (dorev && (binnum > binnum_old) ) { - tbin = binnum_old + 1; - while (tbin <= binnum) { - rev[tbin] = offset; - tbin++; - } - } - // Update the histogram - hist[binnum] = hist[binnum] + 1; - binnum_old = binnum; - } - } - - tbin = binnum_old + 1; - while (tbin <= nbin) { - if (dorev) { - rev[tbin] = nrev; - } - tbin++; - } - - Py_RETURN_NONE; - -} - -static PyMethodDef chist_methods[] = { - {"chist", (PyCFunction)PyCHist_chist, METH_VARARGS, "histogrammer"}, - {NULL} /* Sentinel */ -}; - -#if PY_MAJOR_VERSION >= 3 - static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, - "_chist", /* m_name */ - "Define c version of histogrammer", /* m_doc */ - -1, /* m_size */ - chist_methods, /* m_methods */ - NULL, /* m_reload */ - NULL, /* m_traverse */ - NULL, /* m_clear */ - NULL, /* m_free */ - }; -#endif - -#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */ -#define PyMODINIT_FUNC void -#endif -PyMODINIT_FUNC -#if PY_MAJOR_VERSION >= 3 -PyInit__chist(void) -#else -init_chist(void) -#endif -{ - PyObject* m; - -#if PY_MAJOR_VERSION >= 3 - m = PyModule_Create(&moduledef); - if (m==NULL) { - return NULL; - } - -#else - - m = Py_InitModule3("_chist", chist_methods, "Define c version of histogrammer."); - - if (m==NULL) { - return; - } -#endif - - import_array(); - -#if PY_MAJOR_VERSION >= 3 - return m; -#endif -} diff --git a/tmp/tmpd_h00cl2.cpp b/tmp/tmpd_h00cl2.cpp deleted file mode 100644 index 111fe4f..0000000 --- a/tmp/tmpd_h00cl2.cpp +++ /dev/null @@ -1,7 +0,0 @@ - - #include - int main() { - std::cout<<"Hello World\n"; - return 0; - } - \ No newline at end of file diff --git a/tmp/tmpv19mw1dh.exe b/tmp/tmpv19mw1dh.exe deleted file mode 100755 index 5b02def6a668c6e93d7e5143cbe936e8c4f21acd..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 35036 zcmeI5e{56N702&;&jAxs>J~_#Elq&7>nI~}LIYX5#t_UB(E&@+k*{6Tz{({vC8k!tekAsZfj-^k;pcsj2rFE@zJ7Oe;XdN zs$Q{Mke%(%Y5T+6%_L=u-!Jt}zaBND;cL$IcVD&D?rFO{+1dW=;?Yb>f82-#!!346 zuD^us?>4(WV`m)X`j)aiye=3H$Y6L)1dVe2?Y8|zZ3krM`pEdEoZjpAAyO_{(6~U> zEL+-Oy$Tq0UpLyC$z%1mEinelgp7fl##IZE^6|Q* zOT3HfGHWh@LP{Rr8jrPNO;gN83mPA_YXXQ0+iu!ojT=Q`3M2}<2>I%5 zLZneR9{O?!rOf&zzckJ}sPHB!y;X?f8`Ry6GK%xqi!U8+f3EC@-&(oq)Qs=7zg>>9 z5{NZstU|58)UD-zO`U~0bNt9D(K?+EAzw6hEl`Z%DWbGJ9xDw6n@JkF;)dQfXW^Vc z7c4YFx(D<2AgBJNPsH_DytFye(N)?O3NfC(@GCcp%k025#WOn?b60Vco%m;e)C0!)Aj zFaajO1egF5U;<2l2{3{G5`m%hzf2BkqNlh-`(yfk(X_clxJcikB;P=1K;LPeLB)}%7fSVs0g8dTmr>WhNAzt~C+~C?Ch8y5ocYC@#oA`kyEm^` zOiXJIv{y4W`ftj%5qmf6?3f#2e;e|*k&z-|wG%qU(~Wsexc~O>%$Fu(%M}qoGI=d zMC@8^!E0A$D#`R%0e@4$b_#qJrBWonnqt7gYKoWqZs=Mb@DrZ#uy*Una5deGU686b1A$CZjJjV zh(US}6Gh_pw|alXk>aramT1&(G2) zWa)Qi>Gx*oGqUstvh-P5`kXAiJWH?4(!ZFcFSPWLA1qKRA1pYeDo+>eQ6>qCJX!{4%kye#{QkBfMGBTpvVIHLKR|?Ll;+$e1cCFPz zp-9P!NGue{r?%7qXKu~B#Ln+S4))@9yoXK?&MM<<{Nu16g>1CXxs$yM<#aCnBUH2^ zr?bQb!U2EOTd1?~*o^=#ZFF3CC6rtVYoM-#*%G#O4%wBk8n_ZzKV~~R2TkqA+PTxt zv8It+33C$2=^QoLU|{CfId@mWI_vF9U@m4k-aDm)`Q!%^U;<2l2`~XBzyz286JP>N zfC(@GCcp%k025#WOn?b60Vco%m;e)C0!)AjFaajO1egF5U;<2l2`~XBzyz286JP>N zfC(@GCcp%k025#WOn?b60Vco%m;e)C0!)AjFoFML0wtJxEo2vje&c-f(gNt#Jpee# z51kVI{-a zj|&ODfk;w5Vg1N`y&i9G(_w#o>5>McLP~dKDB|~p^s4r_F4y?{y5X%~+E6CTg8^@L zv-!bDTsHgSdS#1lpo_P|thIc&D^L^k%Shaa>Ato)|I!BOcKcg>F=@nnK_g!0_j*9`feApk-eKGGdG@k~e%Ftu@1$dvKn|uPHCOsD2W>lQ-)r#;*b}fc=w1s-zqM_; z-+~fMaXWm{ofwqpJ_?w$B+W$Y6~?e_3hI_^<8GDf?0q?Qe~x`2$5yOhqyCF>?BzLj zD97HIV{gx~cjegobL>}h?6)i%Pb<-V9duU*aiKdp=zb2mi-YdrpgTC|yQXjcPRK;a zBuF9TM*Zkc54z9eR{ba2K9GqLEDKr({Ap*lu>a*xn>J)?pDLd#WwhOpw0rve*iNu5 zt!>xCe(Y&C+*WD7mbp7{qiigbUT?!ZxqL~(TuHl6CV;yV0ek^-i^JZDST}P_eyi*~ zy}5mj3}Q4sBND62c(1g)H_Vs#qhXBhie3aYx3`cl+=9(56q2Jmqxf=7jnq5+deo4H zuQ{ZzPB>dah{I z)AtA8KQ`m+@c?s(&<9v|L(v_u(ro>R%o$+h!yWzWC<0iWgo={@V4%zW~!w BjOPFV diff --git a/tmp/tmpvz1p6dh5.os b/tmp/tmpvz1p6dh5.os deleted file mode 100644 index 94d4a425c0cc1c1d464da36310bdd6c1b90091d7..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 3088 zcmb7GU2GIp6uvXFZMRU{Vz++?X$d|MbxF4^rF}3gw5t#ZHr+L#0=LuIVY@oJ+wM-G zq$U)2X#x)|@<5OTAK4WXF&fjv#5VY#;Q=KEYhuFIkXY03Kt$-0)cT$IS*JrF@n-LQ zbI!fzJAbox=GK$@|5Pw0FbuksU^a-jG-F~kVon&VhWONFK{8dHWP@SLWw@v5AteUq zx;mX6QLd4TiPd^XP^%P4}x ztPq_v#e2?`jGm4sW0uCvH)_T8g0=oR70NN(+~7+khGa$OBgy0OWR&xLJ&x^c!s6?< zlr!KXl@~>$x?!xMY=xo$(`*GgCShW1#jM5*!ibQJXkRbn4Z!MhTtK{*beJ@VpCW>fQ0O^)mrHV3lbdd~mv{j0FUdMHUn`m6z%a8#9JrArBpKcoCO(BVzL)CSk{&CGugd z&7Eh{XT;jG-n#L=Es}Rl(%ybvTiC^>*5Ci|d|p~_ z#8`mu-O|o49&Z-1`7#E+M%1sRMR|csSc!@MeKLP+M9O~9ICi13WLq(MB$p?!|iW^XF?S6 zLX@%ZM9+^R`(DJ~_-E0R7lk>x#UpHvl@(zfn}awfRNNEUJ>jysv|H7*RP%vUT8pl+ z_Y%Hl3iOv=H7q@g!n}&H;u&P0kCT?iao2jkbv`u!YZ>$VfH+tD73Nvu_WO!B6iFoH2gl$GgO@TQ~%)f!q%y-d-T$^9O-10;?_A3-m!=Kvp=j z{NDg+{SzSS<-Z8T`K8gH19%vR5yKhnr*pp&NN1oI2o46SeRM_wUxzK7mm1ax4<7^~ z^wM77dYA@**jL(V>02y)v!w?Z>t$lK&ntQ>yy6AnA`Bjbq!MV4o4Bb&t?c*5544VH z*Z-DF{@Rl5`d=s6VgD(~=(7F(xC&oQ!VI(fco+VvAlrQ(x7Mq{LE|=Nl_Nfi* w z{7yNH{Z4sneia+RzhkQHzqxY-iIhTCMWuJ1%^R9Sqxa{#^ z@ore}($#b#o|JX<7%8G!uYvo+dK>3~&Uh-r2jq;}1}UTJxQ`CRbuQyIucVThcvMaE zw5kuLQ^!N5;a--qj+2~MtG!!heF})ymxQ7ojV;AtF%qgR5=}bAI`ZlD)Rv~7Ef|T)|Ek9 zrBdaJ3QbPNAZi-E8Lc*87mMvn40 z=w&^X-eYIcVP+BTiz*6$BbG_^b!R>R;DLM&(Ck zMb)X=Exa`zHMf~h$TmekF@*kv485i&hT0fiTiKi4-Hwl8O>H&wfv5BVIOVa?Q9h)m TGpVGk#q|^XSc|cnjyU!&e8`$D