From b860830e008de43fc7ec4eeb50e9e87d4e3a392a Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 09:01:14 -0400 Subject: [PATCH 01/14] allow static to static frame transformations --- .../frame/builtin/transformations.py | 31 ++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/gala/potential/frame/builtin/transformations.py b/gala/potential/frame/builtin/transformations.py index af8ceece7..0bd298246 100644 --- a/gala/potential/frame/builtin/transformations.py +++ b/gala/potential/frame/builtin/transformations.py @@ -3,7 +3,8 @@ import numpy as np # Gala -from ....dynamics import Orbit +from gala.dynamics import Orbit +from gala.units import DimensionlessUnitSystem __all__ = ['static_to_constantrotating', 'constantrotating_to_static'] @@ -145,3 +146,31 @@ def constantrotating_to_static(frame_r, frame_i, w, t=None): """ return _constantrotating_static_helper(frame_r=frame_r, frame_i=frame_i, w=w, t=t, sign=-1.) + + +def static_to_static(frame_r, frame_i, w, t=None): + """ + No-op transform + + Parameters + ---------- + frame_i : `~gala.potential.StaticFrame` + frame_r : `~gala.potential.ConstantRotatingFrame` + w : `~gala.dynamics.PhaseSpacePosition`, `~gala.dynamics.Orbit` + t : quantity_like (optional) + Required if input coordinates are just a phase-space position. + + Returns + ------- + pos : `~astropy.units.Quantity` + Position in static, inertial frame. + vel : `~astropy.units.Quantity` + Velocity in static, inertial frame. + """ + tmp = [isinstance(frame_r.units, DimensionlessUnitSystem), + isinstance(frame_i.units, DimensionlessUnitSystem)] + if not all(tmp) and any(tmp): + raise ValueError( + "StaticFrame to StaticFrame transformations are only allowed if " + "both unit systems are physical, or both are dimensionless.") + return w.pos.xyz, w.vel.d_xyz From 3487eb65c463adfde33f0959acee14133ff86aea Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 09:01:30 -0400 Subject: [PATCH 02/14] fix galpy tests --- gala/dynamics/core.py | 8 ++++---- gala/dynamics/tests/test_core.py | 12 +++++------- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/gala/dynamics/core.py b/gala/dynamics/core.py index 6cdb630b0..bdf6c5249 100644 --- a/gala/dynamics/core.py +++ b/gala/dynamics/core.py @@ -874,8 +874,8 @@ def to_galpy_orbit(self, ro=None, vo=None): galpy_orbit : `galpy.orbit.Orbit` """ - import galpy from galpy.orbit import Orbit + from galpy.util.config import __config__ as galpy_config if self.frame is not None: from ..potential import StaticFrame @@ -884,11 +884,11 @@ def to_galpy_orbit(self, ro=None, vo=None): w = self if ro is None: - ro = galpy.config.__config__.getfloat('normalization', 'ro') + ro = galpy_config.getfloat('normalization', 'ro') ro = ro * u.kpc if vo is None: - vo = galpy.config.__config__.getfloat('normalization', 'vo') + vo = galpy_config.getfloat('normalization', 'vo') vo = vo * u.km/u.s # PhaseSpacePosition or Orbit: @@ -904,7 +904,7 @@ def to_galpy_orbit(self, ro=None, vo=None): o = Orbit(np.array([R, vR, vT, z, vz, phi]).T, ro=ro, vo=vo) - if hasattr(w, 't'): + if hasattr(w, 't') and w.t is not None: o.t = w.t.to_value(ro / vo) return o diff --git a/gala/dynamics/tests/test_core.py b/gala/dynamics/tests/test_core.py index dc30ec753..aa728f5d2 100644 --- a/gala/dynamics/tests/test_core.py +++ b/gala/dynamics/tests/test_core.py @@ -393,9 +393,9 @@ def test_frame_transform(): @pytest.mark.parametrize('obj', [ PhaseSpacePosition([1, 2, 3.]*u.kpc, [1, 2, 3.]*u.km/u.s), PhaseSpacePosition([1, 2, 3.]*u.kpc, [1, 2, 3.]*u.km/u.s, - StaticFrame(galactic)), + StaticFrame(units=galactic)), PhaseSpacePosition([1, 2, 3.]*u.kpc, [1, 2, 3.]*u.km/u.s, - ConstantRotatingFrame([1., 0, 0]*u.rad/u.Myr, + ConstantRotatingFrame(Omega=[1., 0, 0]*u.rad/u.Myr, units=galactic)), ]) def test_io(tmpdir, obj): @@ -412,11 +412,9 @@ def test_io(tmpdir, obj): @pytest.mark.parametrize('obj', [ - PhaseSpacePosition([1,2,3.]*u.kpc, [1,2,3.]*u.km/u.s), - PhaseSpacePosition([1,2,3.]*u.kpc, [1,2,3.]*u.km/u.s, - StaticFrame(galactic)), - PhaseSpacePosition([1,2,3.]*u.kpc, [1,2,3.]*u.km/u.s, - ConstantRotatingFrame([1.,0,0]*u.rad/u.Myr, galactic)), + PhaseSpacePosition([1, 2, 3.]*u.kpc, [1, 2, 3.]*u.km/u.s), + PhaseSpacePosition([1, 2, 3.]*u.kpc, [1, 2, 3.]*u.km/u.s, + StaticFrame(units=galactic)), ]) @pytest.mark.skipif(not HAS_GALPY, reason="requires galpy to run this test") From d971e8b46ba8c834dbbb5555de33872538c39cfa Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Wed, 14 Oct 2020 22:27:44 -0400 Subject: [PATCH 03/14] a true hail mary --- .github/workflows/windows-tests.yml | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/.github/workflows/windows-tests.yml b/.github/workflows/windows-tests.yml index 9043398dc..fe4e2e176 100644 --- a/.github/workflows/windows-tests.yml +++ b/.github/workflows/windows-tests.yml @@ -3,10 +3,10 @@ name: Windows-tests on: push: branches: - - DISABLED + - main pull_request: branches: - - DISABLED + - main jobs: build: @@ -21,14 +21,12 @@ jobs: steps: - uses: actions/checkout@v2 - - uses: goanpeca/setup-miniconda@v1 - with: - auto-update-conda: true - python-version: ${{ matrix.python-version }} - name: Install dependencies run: | - conda install -c conda-forge -q gsl libpython + choco install python --version 3.8 + choco --yes install GnuWin + python -m pip install --upgrade pip python -m pip install -e .[test] python -m pip install tox From 76049b2db48225ad94865f3a047d641623fdc608 Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 08:08:55 -0400 Subject: [PATCH 04/14] sigh, use malloc instead of C99 variable len array --- gala/potential/potential/src/cpotential.c | 27 ++++++++++++++++++----- 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/gala/potential/potential/src/cpotential.c b/gala/potential/potential/src/cpotential.c index 4fe8e7669..0fb2acc85 100644 --- a/gala/potential/potential/src/cpotential.c +++ b/gala/potential/potential/src/cpotential.c @@ -1,3 +1,4 @@ +#include #include #include "cpotential.h" @@ -33,9 +34,11 @@ void apply_rotate(double *q_in, double *R, int n_dim, int transpose, void apply_shift_rotate(double *q_in, double *q0, double *R, int n_dim, int transpose, double *q_out) { - double tmp[n_dim]; + double *tmp; int j; + tmp = (double*) malloc(n_dim * sizeof(double)); + // Shift to the specified origin for (j=0; j < n_dim; j++) { tmp[j] = q_in[j] - q0[j]; @@ -43,13 +46,16 @@ void apply_shift_rotate(double *q_in, double *q0, double *R, int n_dim, // Apply rotation matrix apply_rotate(&tmp[0], R, n_dim, transpose, q_out); + + free(tmp); } double c_potential(CPotential *p, double t, double *qp) { double v = 0; int i, j; - double qp_trans[p->n_dim]; + double *qp_trans; + qp_trans = (double*) malloc(p->n_dim * sizeof(double)); for (i=0; i < p->n_components; i++) { for (j=0; j < p->n_dim; j++) @@ -58,6 +64,7 @@ double c_potential(CPotential *p, double t, double *qp) { &qp_trans[0]); v = v + (p->value)[i](t, (p->parameters)[i], &qp_trans[0], p->n_dim); } + free(qp_trans); return v; } @@ -66,7 +73,8 @@ double c_potential(CPotential *p, double t, double *qp) { double c_density(CPotential *p, double t, double *qp) { double v = 0; int i, j; - double qp_trans[p->n_dim]; + double *qp_trans; + qp_trans = (double*) malloc(p->n_dim * sizeof(double)); for (i=0; i < p->n_components; i++) { for (j=0; j < p->n_dim; j++) @@ -75,6 +83,7 @@ double c_density(CPotential *p, double t, double *qp) { &qp_trans[0]); v = v + (p->density)[i](t, (p->parameters)[i], &qp_trans[0], p->n_dim); } + free(qp_trans); return v; } @@ -82,8 +91,10 @@ double c_density(CPotential *p, double t, double *qp) { void c_gradient(CPotential *p, double t, double *qp, double *grad) { int i, j; - double qp_trans[p->n_dim]; - double tmp_grad[p->n_dim]; + double *qp_trans; + double *tmp_grad; + qp_trans = (double*) malloc(p->n_dim * sizeof(double)); + tmp_grad = (double*) malloc(p->n_dim * sizeof(double)); for (i=0; i < p->n_dim; i++) { grad[i] = 0.; @@ -104,12 +115,15 @@ void c_gradient(CPotential *p, double t, double *qp, double *grad) { apply_rotate(&tmp_grad[0], (p->R)[i], p->n_dim, 1, &grad[0]); } + free(qp_trans); + free(tmp_grad); } void c_hessian(CPotential *p, double t, double *qp, double *hess) { int i; - double qp_trans[p->n_dim]; + double *qp_trans; + qp_trans = (double*) malloc(p->n_dim * sizeof(double)); for (i=0; i < pow(p->n_dim,2); i++) { hess[i] = 0.; @@ -127,6 +141,7 @@ void c_hessian(CPotential *p, double t, double *qp, double *hess) { // - Hessian calculation for potentials with rotations are disabled } + free(qp_trans); } From 99389938c4759bc2f9762794f6a22efef70cc5f9 Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 08:15:00 -0400 Subject: [PATCH 05/14] no-progress --- .github/workflows/windows-tests.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/windows-tests.yml b/.github/workflows/windows-tests.yml index fe4e2e176..1deaedc1d 100644 --- a/.github/workflows/windows-tests.yml +++ b/.github/workflows/windows-tests.yml @@ -24,8 +24,8 @@ jobs: - name: Install dependencies run: | - choco install python --version 3.8 - choco --yes install GnuWin + choco install --no-progress python --version 3.8 + choco --yes install --no-progress GnuWin python -m pip install --upgrade pip python -m pip install -e .[test] python -m pip install tox From 27b3398e69aec50eeaf15c704f9d83f10b326c24 Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 08:20:54 -0400 Subject: [PATCH 06/14] remove more var len array stuff --- gala/integrate/cyintegrators/dopri/dop853.c | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gala/integrate/cyintegrators/dopri/dop853.c b/gala/integrate/cyintegrators/dopri/dop853.c index 5420303c7..099021114 100644 --- a/gala/integrate/cyintegrators/dopri/dop853.c +++ b/gala/integrate/cyintegrators/dopri/dop853.c @@ -1002,7 +1002,8 @@ void Fwrapper_direct_nbody (unsigned full_ndim, double t, double *w, double *f, // Note: only really works with a static frame! This should be enforced int i, j, k; unsigned ndim = full_ndim / norbits; // phase-space dimensionality - double f2[ndim/2]; + double *f2; + f2 = (double*) malloc(ndim/2 * sizeof(double)); for (i=0; i < norbits; i++) { // call gradient function @@ -1029,6 +1030,8 @@ void Fwrapper_direct_nbody (unsigned full_ndim, double t, double *w, double *f, } } + free(f2); + } /* Needed for Lyapunov */ From 57daa0de67e1773b4c3857e2b8600d82f198827b Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 08:44:30 -0400 Subject: [PATCH 07/14] define M_PI in headers --- gala/potential/potential/builtin/builtin_potentials.h | 4 ++++ gala/potential/scf/src/bfe_helper.h | 4 ++++ gala/potential/scf/src/coeff_helper.h | 4 ++++ 3 files changed, 12 insertions(+) diff --git a/gala/potential/potential/builtin/builtin_potentials.h b/gala/potential/potential/builtin/builtin_potentials.h index e5624db92..8230efb97 100644 --- a/gala/potential/potential/builtin/builtin_potentials.h +++ b/gala/potential/potential/builtin/builtin_potentials.h @@ -1,3 +1,7 @@ +#ifndef M_PI + #define M_PI 3.14159265358979323846 +#endif + extern double nan_density(double t, double *pars, double *q, int n_dim); extern double nan_value(double t, double *pars, double *q, int n_dim); extern void nan_gradient(double t, double *pars, double *q, int n_dim, double *grad); diff --git a/gala/potential/scf/src/bfe_helper.h b/gala/potential/scf/src/bfe_helper.h index c7d9a6817..da8b67253 100644 --- a/gala/potential/scf/src/bfe_helper.h +++ b/gala/potential/scf/src/bfe_helper.h @@ -8,3 +8,7 @@ extern double phi_nlm(double s, double phi, double X, int n, int l, int m); extern void sph_grad_phi_nlm(double s, double phi, double X, int n, int l, int m, int lmax, double *sphgrad); // #endif + +#ifndef M_PI + #define M_PI 3.14159265358979323846 +#endif \ No newline at end of file diff --git a/gala/potential/scf/src/coeff_helper.h b/gala/potential/scf/src/coeff_helper.h index bc55d5b0f..dd526273c 100644 --- a/gala/potential/scf/src/coeff_helper.h +++ b/gala/potential/scf/src/coeff_helper.h @@ -12,3 +12,7 @@ extern void c_STnlm_discrete(double *s, double *phi, double *X, double *m_k, int extern void c_STnlm_var_discrete(double *s, double *phi, double *X, double *m_k, int K, int n, int l, int m, double *ST_var); + +#ifndef M_PI + #define M_PI 3.14159265358979323846 +#endif \ No newline at end of file From 4f11c3b00689d52a3ecdcea935ee1524cbbe2e64 Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 09:25:52 -0400 Subject: [PATCH 08/14] include --- gala/potential/potential/builtin/builtin_potentials.c | 1 + 1 file changed, 1 insertion(+) diff --git a/gala/potential/potential/builtin/builtin_potentials.c b/gala/potential/potential/builtin/builtin_potentials.c index fb0b92659..ad7659557 100644 --- a/gala/potential/potential/builtin/builtin_potentials.c +++ b/gala/potential/potential/builtin/builtin_potentials.c @@ -1,6 +1,7 @@ #include #include #include "extra_compile_macros.h" +#include "builtin_potentials.h" #if USE_GSL == 1 #include From fdf59db273ba7e17791f6f34b8c034de00fa5ab5 Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 09:33:29 -0400 Subject: [PATCH 09/14] fix incorrect redeclaration in header --- gala/potential/potential/builtin/builtin_potentials.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gala/potential/potential/builtin/builtin_potentials.h b/gala/potential/potential/builtin/builtin_potentials.h index 8230efb97..9ec910e9f 100644 --- a/gala/potential/potential/builtin/builtin_potentials.h +++ b/gala/potential/potential/builtin/builtin_potentials.h @@ -48,7 +48,7 @@ extern void powerlawcutoff_hessian(double t, double *pars, double *q, int n_dim, extern double stone_value(double t, double *pars, double *q, int n_dim); extern void stone_gradient(double t, double *pars, double *q, int n_dim, double *grad); -extern void stone_density(double t, double *pars, double *q, int n_dim); +extern double stone_density(double t, double *pars, double *q, int n_dim); extern void stone_hessian(double t, double *pars, double *q, int n_dim, double *hess); extern double sphericalnfw_value(double t, double *pars, double *q, int n_dim); From 8f5cff395344ecb505774247c5b767d53c0c548f Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 09:53:09 -0400 Subject: [PATCH 10/14] only define powerlaw in header if using gsl --- gala/potential/potential/builtin/builtin_potentials.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/gala/potential/potential/builtin/builtin_potentials.h b/gala/potential/potential/builtin/builtin_potentials.h index 9ec910e9f..6a22b2f49 100644 --- a/gala/potential/potential/builtin/builtin_potentials.h +++ b/gala/potential/potential/builtin/builtin_potentials.h @@ -1,3 +1,5 @@ +#include "extra_compile_macros.h" + #ifndef M_PI #define M_PI 3.14159265358979323846 #endif @@ -41,10 +43,12 @@ extern void jaffe_gradient(double t, double *pars, double *q, int n_dim, double extern double jaffe_density(double t, double *pars, double *q, int n_dim); extern void jaffe_hessian(double t, double *pars, double *q, int n_dim, double *hess); +#if USE_GSL == 1 extern double powerlawcutoff_value(double t, double *pars, double *q, int n_dim); extern void powerlawcutoff_gradient(double t, double *pars, double *q, int n_dim, double *grad); extern double powerlawcutoff_density(double t, double *pars, double *q, int n_dim); extern void powerlawcutoff_hessian(double t, double *pars, double *q, int n_dim, double *hess); +#endif extern double stone_value(double t, double *pars, double *q, int n_dim); extern void stone_gradient(double t, double *pars, double *q, int n_dim, double *grad); From 9396f8b1f49ba641ab8f369c5a939896a90876fc Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 09:56:59 -0400 Subject: [PATCH 11/14] remove --yes --- .github/workflows/windows-tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/windows-tests.yml b/.github/workflows/windows-tests.yml index 1deaedc1d..2f927aae7 100644 --- a/.github/workflows/windows-tests.yml +++ b/.github/workflows/windows-tests.yml @@ -25,7 +25,7 @@ jobs: - name: Install dependencies run: | choco install --no-progress python --version 3.8 - choco --yes install --no-progress GnuWin + choco install --no-progress GnuWin python -m pip install --upgrade pip python -m pip install -e .[test] python -m pip install tox From 580cb12573f55a203ea31ba965f5f67bab37e38d Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 10:20:17 -0400 Subject: [PATCH 12/14] try a compile time def --- gala/potential/potential/builtin/cybuiltin.pyx | 16 ++++++++++------ setup.py | 4 ++++ 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/gala/potential/potential/builtin/cybuiltin.pyx b/gala/potential/potential/builtin/cybuiltin.pyx index fb50414a3..ad8534d6b 100644 --- a/gala/potential/potential/builtin/cybuiltin.pyx +++ b/gala/potential/potential/builtin/cybuiltin.pyx @@ -73,11 +73,6 @@ cdef extern from "potential/potential/builtin/builtin_potentials.h": double jaffe_density(double t, double *pars, double *q, int n_dim) nogil void jaffe_hessian(double t, double *pars, double *q, int n_dim, double *hess) nogil - double powerlawcutoff_value(double t, double *pars, double *q, int n_dim) nogil - void powerlawcutoff_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil - double powerlawcutoff_density(double t, double *pars, double *q, int n_dim) nogil - void powerlawcutoff_hessian(double t, double *pars, double *q, int n_dim, double *hess) nogil - double stone_value(double t, double *pars, double *q, int n_dim) nogil void stone_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil double stone_density(double t, double *pars, double *q, int n_dim) nogil @@ -119,6 +114,15 @@ cdef extern from "potential/potential/builtin/builtin_potentials.h": double longmuralibar_density(double t, double *pars, double *q, int n_dim) nogil void longmuralibar_hessian(double t, double *pars, double *q, int n_dim, double *hess) nogil + +IF USE_GSL_C == 1: + cdef extern from "potential/potential/builtin/builtin_potentials.h": + double powerlawcutoff_value(double t, double *pars, double *q, int n_dim) nogil + void powerlawcutoff_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil + double powerlawcutoff_density(double t, double *pars, double *q, int n_dim) nogil + void powerlawcutoff_hessian(double t, double *pars, double *q, int n_dim, double *hess) nogil + + __all__ = ['NullPotential', 'HenonHeilesPotential', # Misc. potentials 'KeplerPotential', 'HernquistPotential', 'IsochronePotential', 'PlummerPotential', 'JaffePotential', 'StonePotential', 'PowerLawCutoffPotential', # Spherical models @@ -442,7 +446,7 @@ cdef class PowerLawCutoffWrapper(CPotentialWrapper): np.ascontiguousarray(q0), np.ascontiguousarray(R)) - if USE_GSL == 1: + IF USE_GSL_C == 1: self.cpotential.value[0] = (powerlawcutoff_value) self.cpotential.density[0] = (powerlawcutoff_density) self.cpotential.gradient[0] = (powerlawcutoff_gradient) diff --git a/setup.py b/setup.py index 1c820d273..ad54939c0 100755 --- a/setup.py +++ b/setup.py @@ -168,10 +168,14 @@ with open(extra_compile_macros_file, 'w') as f: if gsl_version is not None: + GSL_ENABLED = True f.writelines(['#define USE_GSL 1']) else: + GSL_ENABLED = False f.writelines(['#define USE_GSL 0']) +for ext in extensions: + ext.cython_compile_time_env = {'USE_GSL_C': int(GSL_ENABLED)} setup(use_scm_version={'write_to': os.path.join('gala', 'version.py'), 'write_to_template': VERSION_TEMPLATE}, From b5a4c8bd90589bbe3306416838cbdcf4f2a73302 Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 12:38:30 -0400 Subject: [PATCH 13/14] use compile-time USE_GSL_C for windows compiler... --- gala/potential/scf/bfe.pyx | 37 ++++++++++++++------------- gala/potential/scf/computecoeff.pyx | 20 ++++++++------- gala/potential/scf/src/bfe.h | 4 +++ gala/potential/scf/src/bfe_helper.h | 7 ++--- gala/potential/scf/src/coeff_helper.h | 8 +++++- 5 files changed, 45 insertions(+), 31 deletions(-) diff --git a/gala/potential/scf/bfe.pyx b/gala/potential/scf/bfe.pyx index 88b4037dc..13d102d49 100644 --- a/gala/potential/scf/bfe.pyx +++ b/gala/potential/scf/bfe.pyx @@ -20,21 +20,22 @@ cimport cython cdef extern from "extra_compile_macros.h": int USE_GSL -cdef extern from "scf/src/bfe_helper.h": - double rho_nlm(double s, double phi, double X, int n, int l, int m) nogil - double phi_nlm(double s, double phi, double X, int n, int l, int m) nogil - double sph_grad_phi_nlm(double s, double phi, double X, int n, int l, int m, double *grad) nogil - -cdef extern from "scf/src/bfe.h": - void scf_density_helper(double *xyz, int K, double M, double r_s, - double *Snlm, double *Tnlm, - int nmax, int lmax, double *dens) nogil - void scf_potential_helper(double *xyz, int K, double G, double M, double r_s, - double *Snlm, double *Tnlm, - int nmax, int lmax, double *potv) nogil - void scf_gradient_helper(double *xyz, int K, double G, double M, double r_s, - double *Snlm, double *Tnlm, - int nmax, int lmax, double *grad) nogil +IF USE_GSL_C == 1: + cdef extern from "scf/src/bfe_helper.h": + double rho_nlm(double s, double phi, double X, int n, int l, int m) nogil + double phi_nlm(double s, double phi, double X, int n, int l, int m) nogil + double sph_grad_phi_nlm(double s, double phi, double X, int n, int l, int m, double *grad) nogil + + cdef extern from "scf/src/bfe.h": + void scf_density_helper(double *xyz, int K, double M, double r_s, + double *Snlm, double *Tnlm, + int nmax, int lmax, double *dens) nogil + void scf_potential_helper(double *xyz, int K, double G, double M, double r_s, + double *Snlm, double *Tnlm, + int nmax, int lmax, double *potv) nogil + void scf_gradient_helper(double *xyz, int K, double G, double M, double r_s, + double *Snlm, double *Tnlm, + int nmax, int lmax, double *grad) nogil __all__ = ['density', 'potential', 'gradient'] @@ -82,7 +83,7 @@ cpdef density(double[:,::1] xyz, int nmax = Snlm.shape[0]-1 int lmax = Snlm.shape[1]-1 - if USE_GSL == 1: + IF USE_GSL_C == 1: scf_density_helper(&xyz[0,0], ncoords, M, r_s, &Snlm[0,0,0], &Tnlm[0,0,0], nmax, lmax, &dens[0]) @@ -134,7 +135,7 @@ cpdef potential(double[:,::1] xyz, int nmax = Snlm.shape[0]-1 int lmax = Snlm.shape[1]-1 - if USE_GSL == 1: + IF USE_GSL_C == 1: scf_potential_helper(&xyz[0,0], ncoords, G, M, r_s, &Snlm[0,0,0], &Tnlm[0,0,0], nmax, lmax, &potv[0]) @@ -187,7 +188,7 @@ cpdef gradient(double[:,::1] xyz, int nmax = Snlm.shape[0]-1 int lmax = Snlm.shape[1]-1 - if USE_GSL == 1: + IF USE_GSL_C == 1: scf_gradient_helper(&xyz[0,0], ncoords, G, M, r_s, &Snlm[0,0,0], &Tnlm[0,0,0], nmax, lmax, &grad[0,0]) diff --git a/gala/potential/scf/computecoeff.pyx b/gala/potential/scf/computecoeff.pyx index 3c8c010d0..75e0af3de 100644 --- a/gala/potential/scf/computecoeff.pyx +++ b/gala/potential/scf/computecoeff.pyx @@ -20,11 +20,13 @@ cdef extern from "math.h": double cos(double x) nogil double sin(double x) nogil -cdef extern from "scf/src/coeff_helper.h": - double c_Snlm_integrand(double phi, double X, double xsi, double density, int n, int l, int m) - double c_Tnlm_integrand(double phi, double X, double xsi, double density, int n, int l, int m) - void c_STnlm_discrete(double *s, double *phi, double *X, double *m_k, int K, int n, int l, int m, double *ST) - void c_STnlm_var_discrete(double *s, double *phi, double *X, double *m_k, int K, int n, int l, int m, double *ST_var) +IF USE_GSL_C == 1: + cdef extern from "scf/src/coeff_helper.h": + double c_Snlm_integrand(double phi, double X, double xsi, double density, int n, int l, int m) + double c_Tnlm_integrand(double phi, double X, double xsi, double density, int n, int l, int m) + void c_STnlm_discrete(double *s, double *phi, double *X, double *m_k, int K, int n, int l, int m, double *ST) + void c_STnlm_var_discrete(double *s, double *phi, double *X, double *m_k, int K, int n, int l, int m, double *ST_var) + __all__ = ['Snlm_integrand', 'Tnlm_integrand'] @@ -40,7 +42,7 @@ cpdef Snlm_integrand(double phi, double X, double xsi, double z = r * X double val = 0. - if USE_GSL == 1: + IF USE_GSL_C == 1: val = c_Snlm_integrand(phi, X, xsi, density_func(x, y, z, *args) / M * r_s*r_s*r_s, n, l, m) @@ -58,7 +60,7 @@ cpdef Tnlm_integrand(double phi, double X, double xsi, double z = r * X double val = 0. - if USE_GSL == 1: + IF USE_GSL_C == 1: val = c_Tnlm_integrand(phi, X, xsi, density_func(x, y, z, *args) / M * r_s*r_s*r_s, n, l, m) @@ -71,7 +73,7 @@ cpdef STnlm_discrete(double[::1] s, double[::1] phi, double[::1] X, double[::1] ST = np.zeros(2) int K = s.size - if USE_GSL == 1: + IF USE_GSL_C == 1: c_STnlm_discrete(&s[0], &phi[0], &X[0], &m_k[0], K, n, l, m, &ST[0]) return ST @@ -83,7 +85,7 @@ cpdef STnlm_var_discrete(double[::1] s, double[::1] phi, double[::1] X, double[::1] ST_var = np.zeros(2) int K = s.size - if USE_GSL == 1: + IF USE_GSL_C == 1: c_STnlm_var_discrete(&s[0], &phi[0], &X[0], &m_k[0], K, n, l, m, &ST_var[0]) return ST_var diff --git a/gala/potential/scf/src/bfe.h b/gala/potential/scf/src/bfe.h index 10121b915..5f2f22ab7 100644 --- a/gala/potential/scf/src/bfe.h +++ b/gala/potential/scf/src/bfe.h @@ -1,3 +1,6 @@ +#include "extra_compile_macros.h" + +#if USE_GSL == 1 extern void scf_density_helper(double *xyz, int K, double M, double r_s, double *Snlm, double *Tnlm, int nmax, int lmax, double *dens); @@ -15,3 +18,4 @@ extern void scf_gradient_helper(double *xyz, int K, extern double scf_value(double t, double *pars, double *q, int n_dim); extern void scf_gradient(double t, double *pars, double *q, int n_dim, double *grad); extern double scf_density(double t, double *pars, double *q, int n_dim); +#endif diff --git a/gala/potential/scf/src/bfe_helper.h b/gala/potential/scf/src/bfe_helper.h index da8b67253..408f184b8 100644 --- a/gala/potential/scf/src/bfe_helper.h +++ b/gala/potential/scf/src/bfe_helper.h @@ -1,5 +1,6 @@ -// #ifndef _BFE_HELPER_ -// #define _BFE_HELPER_ +#include "extra_compile_macros.h" + +#if USE_GSL == 1 extern double rho_nl(double s, int n, int l); extern double rho_nlm(double s, double phi, double X, int n, int l, int m); @@ -7,7 +8,7 @@ extern double phi_nl(double s, int n, int l); extern double phi_nlm(double s, double phi, double X, int n, int l, int m); extern void sph_grad_phi_nlm(double s, double phi, double X, int n, int l, int m, int lmax, double *sphgrad); -// #endif +#endif #ifndef M_PI #define M_PI 3.14159265358979323846 diff --git a/gala/potential/scf/src/coeff_helper.h b/gala/potential/scf/src/coeff_helper.h index dd526273c..47dc7b4ec 100644 --- a/gala/potential/scf/src/coeff_helper.h +++ b/gala/potential/scf/src/coeff_helper.h @@ -1,3 +1,6 @@ +#include "extra_compile_macros.h" + +#if USE_GSL == 1 extern double STnlm_integrand_help(double phi, double X, double xsi, double density, int n, int l, int m); @@ -13,6 +16,9 @@ extern void c_STnlm_discrete(double *s, double *phi, double *X, double *m_k, int extern void c_STnlm_var_discrete(double *s, double *phi, double *X, double *m_k, int K, int n, int l, int m, double *ST_var); +#endif + + #ifndef M_PI #define M_PI 3.14159265358979323846 -#endif \ No newline at end of file +#endif From c63e8e66f983e0e5b9fb6354a31bc784d8626abe Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan Date: Thu, 15 Oct 2020 13:24:23 -0400 Subject: [PATCH 14/14] use compile-time gsl in bfe --- gala/potential/scf/bfe_class.pyx | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/gala/potential/scf/bfe_class.pyx b/gala/potential/scf/bfe_class.pyx index 5af97980b..eaf2d94ef 100644 --- a/gala/potential/scf/bfe_class.pyx +++ b/gala/potential/scf/bfe_class.pyx @@ -33,10 +33,11 @@ cdef extern from "src/funcdefs.h": ctypedef double (*energyfunc)(double t, double *pars, double *q, int n_dim) nogil ctypedef void (*gradientfunc)(double t, double *pars, double *q, int n_dim, double *grad) nogil -cdef extern from "scf/src/bfe.h": - double scf_value(double t, double *pars, double *q, int n_dim) nogil - double scf_density(double t, double *pars, double *q, int n_dim) nogil - void scf_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil +IF USE_GSL_C == 1: + cdef extern from "scf/src/bfe.h": + double scf_value(double t, double *pars, double *q, int n_dim) nogil + double scf_density(double t, double *pars, double *q, int n_dim) nogil + void scf_gradient(double t, double *pars, double *q, int n_dim, double *grad) nogil __all__ = ['SCFPotential'] @@ -46,7 +47,7 @@ cdef class SCFWrapper(CPotentialWrapper): self.init([G] + list(parameters), np.ascontiguousarray(q0), np.ascontiguousarray(R)) - if USE_GSL == 1: + IF USE_GSL_C == 1: self.cpotential.value[0] = (scf_value) self.cpotential.density[0] = (scf_density) self.cpotential.gradient[0] = (scf_gradient)