[WIP] Another attempt to fix the windows CI

.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 --no-progress python --version 3.8
+          choco install --no-progress GnuWin
+          python -m pip install --upgrade pip
           python -m pip install -e .[test]
           python -m pip install tox
 

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

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")

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 */

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

gala/potential/potential/builtin/builtin_potentials.c

@@ -1,6 +1,7 @@
 #include <math.h>
 #include <string.h>
 #include "extra_compile_macros.h"
+#include "builtin_potentials.h"
 
 #if USE_GSL == 1
 #include <gsl/gsl_sf_gamma.h>

gala/potential/potential/builtin/builtin_potentials.h

@@ -1,3 +1,9 @@
+#include "extra_compile_macros.h"
+
+#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);
@@ -37,14 +43,16 @@ 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);
-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);

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] = <energyfunc>(powerlawcutoff_value)
             self.cpotential.density[0] = <densityfunc>(powerlawcutoff_density)
             self.cpotential.gradient[0] = <gradientfunc>(powerlawcutoff_gradient)

gala/potential/potential/src/cpotential.c

@@ -1,3 +1,4 @@
+#include <stdlib.h>
 #include <math.h>
 #include "cpotential.h"
 
@@ -33,23 +34,28 @@ 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];
     }
 
     // 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,15 +83,18 @@ 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;
 }
 
 
 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);
 }
 
 

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])