Update from development repo

.github/workflows/build.yml

@@ -9,11 +9,11 @@ jobs:
 
     strategy:
       matrix:
-        python: ["3.12", "3.11", "3.10", "3.9"]
+        python: ["3.12", "3.11", "3.10"]
 
     steps:
-      - uses: actions/checkout@v1
-      - uses: actions/setup-python@v1
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python }}
 
@@ -30,19 +30,18 @@ jobs:
 
   test-macos:
     name: "Run tests on MacOS"
-    runs-on: macos-12
+    runs-on: macos-latest
     env:
-      # LDFLAGS: "-ld64" # For MacOS 13 and above (XCode CLT 15 and above.)
-      CC: gcc-12
-      CXX: gcc-12
-      FC: gfortran-12
+      CC: gcc-14
+      CXX: gcc-14
+      FC: gfortran-14
       DUCC0_NUM_THREADS: 2
 
     steps:
-      - uses: actions/checkout@v1
-      - uses: actions/setup-python@v1
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
         with:
-          python-version: "3.10"
+          python-version: "3.11"
 
       - name: Install Dependencies (MacOS)
         run: |
@@ -91,7 +90,7 @@ jobs:
     strategy:
       matrix:
         # macos-13 is an intel runner, macos-14 is apple silicon
-        os: [macos-14]
+        os: [macos-latest]
 
     steps:
       - uses: actions/checkout@v4
@@ -113,9 +112,9 @@ jobs:
     name: Build source distribution
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
 
-      - uses: actions/setup-python@v2
+      - uses: actions/setup-python@v5
         name: Install Python
         with:
           python-version: '3.10'

Makefile

@@ -76,3 +76,17 @@ dist: clean ## builds source and wheel package
 
 install: clean ## install the package to the active Python's site-packages
 	pip install .
+
+
+# Symlink-based alternative setup below. Not intended for general use
+# Standard meson build
+SHELL=/bin/bash
+.PHONY: build inline
+inline: build
+	(shopt -s nullglob; cd pixell; rm -f *.so; ln -s ../build/*.so ../build/*.dylib .)
+build: build/build.ninja
+	(cd build; meson compile)
+build/build.ninja:
+	rm -rf build
+	mkdir build
+	meson setup build

cython/cmisc.pyx

@@ -31,7 +31,7 @@ def alm2cl(ainfo, alm, alm2=None):
 	# I used to flatten here to make looping simple, but that caused a copy to be made
 	# when combined with np.broadcast. So instead I will use manual raveling
 	pshape = alm.shape[:-1]
-	npre   = int(np.product(pshape))
+	npre   = int(np.prod(pshape))
 	cdef float[::1]  cl_single_sp, alm_single_sp1, alm_single_sp2
 	cdef double[::1] cl_single_dp, alm_single_dp1, alm_single_dp2
 	cdef int64_t[::1] mstart = np.ascontiguousarray(ainfo.mstart).view(np.int64)

cython/cmisc_core.c

@@ -3,6 +3,7 @@
 #include <stdio.h>
 #include <stdint.h>
 #include <string.h>
+#include <math.h>
 #include <omp.h>
 
 int min(int a, int b) { return a < b ? a : b; }
@@ -269,3 +270,54 @@ void transfer_alm_sp(int lmax1, int mmax1, int64_t * mstart1, float * alm1, int
 		}
 	}
 }
+
+// wcs acceleration
+#define DEG (M_PI/180)
+
+// I pass the wcs information as individual doubles to avoid having to construct
+// numpy arrays on the python side. All of these have the arguments in the same
+// order, regardless of which way they go, so they can be defined in a macro
+//
+// We only implement plain spherical coordinates here - the final coordinate
+// rotation is missing. For cylindrical coordinates this means that we only
+// support dec0 = 0 - the result will be wrong for other values
+#define wcsdef(name) \
+void name(int64_t n, double * restrict dec, double * restrict ra, \
+		double * restrict y, double * restrict x, \
+		double crval0, double crval1, double cdelt0, double cdelt1, \
+		double crpix0, double crpix1) { \
+	double ra0 = crval0*DEG, dec0 = crval1*DEG; \
+	double dra = cdelt0*DEG, ddec = cdelt1*DEG; \
+	double x0  = crpix0-1,   y0   = crpix1-1; \
+	_Pragma("omp parallel for") \
+	for(int64_t i = 0; i < n; i++) {
+#define wcsend } \
+}
+
+wcsdef(wcs_car_sky2pix)
+x[i] = (ra [i]-ra0 )/dra +x0;
+y[i] = (dec[i]-dec0)/ddec+y0; // dec0 should be zero
+wcsend
+
+wcsdef(wcs_car_pix2sky)
+ra [i] = (x[i]-x0)*dra +ra0;
+dec[i] = (y[i]-y0)*ddec+dec0; // dec0 should be zero
+wcsend
+
+wcsdef(wcs_cea_sky2pix)
+x[i] = (ra [i]-ra0 )/dra +x0;
+y[i] = sin(dec[i])/ddec  +y0;
+(void)dec0; // mark dec0 as explicitly unused
+wcsend
+
+wcsdef(wcs_cea_pix2sky)
+ra [i] = (x[i]-x0)*dra +ra0;
+dec[i] = asin((y[i]-y0)*ddec);
+(void)dec0; // mark dec0 as explicitly unused
+wcsend
+
+void rewind_inplace(int64_t n, double * vals, double period, double ref) {
+	_Pragma("omp parallel for")
+	for(int64_t i = 0; i < n; i++)
+		vals[i] = fmod((vals[i]+ref),period)-ref;
+}

cython/distances_core.c

@@ -15,7 +15,7 @@ int xoffs[8] = {-1, +1,  0,  0, +1, -1, +1, -1 };
 double wall_time() { struct timeval tv; gettimeofday(&tv,0); return tv.tv_sec + 1e-6*tv.tv_usec; }
 int max(int a, int b) { return a > b ? a : b; }
 int min(int a, int b) { return a < b ? a : b; }
-int compar_int(int * a, int * b) { return *a-*b; }
+int compar_int(const void * a, const void * b) { return *(int*)a-*(int*)b; }
 int wrap1(int a, int n) { return a < 0 ? a+n : a >= n ? a-n : a; }
 
 // The simple functions are too slow to serve as the basis for a distance transform.

fortran/interpol.F90

@@ -1,6 +1,6 @@
 module fortran
 
-	private :: map_border, calc_weights
+	!private :: map_border, calc_weights
 
 contains
 

meson.build

@@ -45,6 +45,8 @@ incdir_f2py = run_command(
 run_command('make', '-C', 'fortran', check: true)
 
 fortran_include = include_directories(incdir_numpy, incdir_f2py)
+add_project_arguments('-Wno-tabs', language : 'fortran')
+add_project_arguments('-Wno-conversion', language : 'fortran')
 
 fortran_sources = {
     'fortran/interpol_32.f90': '_interpol_32',
@@ -87,7 +89,7 @@ helper_sources = {
 linkables = []
 
 foreach source_name, module_name : helper_sources
-    linkables += shared_library(
+    linkables += static_library(
         module_name,
         source_name,
         install: true,
@@ -120,4 +122,4 @@ endforeach
 # The actual python install itself is left up to a helper build
 # script deifned in pixell/
 subdir('pixell')
-subdir('scripts')
+subdir('scripts')

pixell/aberration.py

@@ -116,7 +116,7 @@ def __init__(self, shape, wcs, dir=dir_equ, beta=beta, spin=[0,2],
 		nthread = int(utils.fallback(utils.getenv("OMP_NUM_THREADS",nthread),0))
 		# 1. Calculate the aberration field. These are tiny
 		alm_dpos = calc_boost_field(-beta, dir, nthread=nthread)
-		# 2. Evaluate these on our target geometry. Hardcoded float64 because of get_deflected_angles
+		# 2. Evaluate these on our target geometry.
 		deflect = enmap.zeros(alm_dpos.shape[:-1]+shape[-2:], wcs, coord_dtype)
 		curvedsky.alm2map(alm_dpos.astype(coord_ctype, copy=False), deflect, spin=1, nthread=nthread)
 		# 3. Calculate the offset angles.

pixell/bunch.py

@@ -1,5 +1,6 @@
 """My own version of bunch, since the standard one lacks tab completion
 and has trouble printing sometimes."""
+import os
 class Bunch:
 	def __init__(self, *args, **kwargs):
 		self._dict = {}
@@ -53,61 +54,102 @@ def __repr__(self):
 # Some simple I/O routines. These can't handle everything that could
 # be in a bunch, but they cover all my most common use cases.
 
-def read(fname, fmt="auto", group=None):
-	if fmt == "auto":
-		if is_hdf_path(fname): fmt = "hdf"
-		else: raise ValueError("Could not infer format for '%s'" % fname)
-	if fmt == "hdf": return read_hdf(fname, group=group)
+def read(fname, fmt="auto", group=None, gmode="dot"):
+	if fmt == "auto": fmt="hdf"
+	if fmt == "hdf": return read_hdf(fname, group=group, gmode=gmode)
 	else: raise ValueError("Unrecognized format '%s'" % fmt)
 
-def write(fname, bunch, fmt="auto", group=None):
-	if fmt == "auto":
-		if is_hdf_path(fname): fmt = "hdf"
-		else: raise ValueError("Could not infer format for '%s'" % fname)
-	if fmt == "hdf": write_hdf(fname, bunch, group=group)
+def write(fname, bunch, fmt="auto", group=None, gmode="dot"):
+	if fmt == "auto": fmt = "hdf"
+	if fmt == "hdf": write_hdf(fname, bunch, group=group, gmode=gmode)
 	else: raise ValueError("Unrecognized format '%s'" % fmt)
 
-def write_hdf(fname, bunch, group=None):
+def read_hdf(fname, group=None, gmode="dot"):
 	import h5py
-	fname, group = split_hdf_path(fname, group)
+	if group is None:
+		fname, group = split_hdf_path(fname, group, mode=gmode)
+	with h5py.File(fname, "r") as hfile:
+		if group: hfile = hfile[group]
+		return read_hdf_recursive(hfile)
+
+def read_hdf_recursive(hfile):
+	import h5py
+	if isinstance(hfile, h5py.Dataset):
+		return hfile[()]
+	else:
+		bunch = Bunch()
+		for key in hfile:
+			bunch[key] = read_hdf_recursive(hfile[key])
+		return bunch
+
+def write_hdf(fname, bunch, group=None, gmode="dot"):
+	import h5py
+	if group is None:
+		fname, group = split_hdf_path(fname, group, mode=gmode)
 	with h5py.File(fname, "w") as hfile:
 		if group: hfile = hfile.create_group(group)
-		for key in bunch:
+		write_hdf_recursive(hfile, bunch)
+
+def write_hdf_recursive(hfile, bunch):
+	for key in bunch:
+		if isinstance(bunch[key],Bunch):
+			hfile.create_group(key)
+			write_hdf_recursive(hfile[key], bunch[key])
+		else:
 			hfile[key] = bunch[key]
 
-def read_hdf(fname, group=None):
-	import h5py
-	bunch = Bunch()
-	fname, group = split_hdf_path(fname, group)
-	with h5py.File(fname, "r") as hfile:
-		if group: hfile = hfile[group]
-		for key in hfile:
-			bunch[key] = hfile[key][()]
-	return bunch
+#def make_safe(val):
+#	import numpy as np
+#	if isinstance(val, np.ndarray):
+#		try: return np.char.encode(val)
+#		except TypeError: pass
+#	elif isinstance(val, str):
+#		return val.encode()
+#	else:
+#		return val
 
 def is_hdf_path(fname):
 	"""Returns true if the fname would be recognized by split_hdf_path"""
-	for suf in [".hdf", ".h5"]:
-		name, _, group = fname.rpartition(suf)
-		if name and (not group or group[0] == "/"): return True
-	return False
+	return True
 
-def split_hdf_path(fname, subgroup=None):
+def split_hdf_path(fname, subgroup=None, mode="dot"):
 	"""Split an hdf path of the form path.hdf/group, where the group part is
 	optional, into the path and the group parts. If subgroup is specified, then
 	it will be appended to the group informaiton. returns fname, group. The
 	fname will be a string, and the group will be a string or None. Raises
-	a ValueError if the fname is not recognized as a hdf file."""
-	for suf in [".hdf", ".h5"]:
-		name, _, group = fname.rpartition(suf)
-		if not name: continue
-		name += suf
-		if not group: return name, subgroup
-		elif group[0] == "/":
-			group = group[1:]
-			if subgroup: group += "/" + subgroup
-			return name, group
-	raise ValueError("Not an hdf path")
+	a ValueError if unsuccessful.
+
+	mode controles how the split is done:
+	* "none": Don't split. fname is returned unmodified
+	* "dot": The last entry in the path given by filename
+	    containing a "." will be taken to be the real
+	    file name, the rest till be the hdf group path.
+	    For example, with a/b/c.d/e/f, a/b/c.d would be returned
+	    as the file name and e/f as the hdf group
+	* "exists": As dot, but based on whether a file with that
+	    name can be found on disk. Seemed like a good idea,
+	    except it doesn't work when writing a new file.
+	"""
+	toks = fname.split("/")
+	if mode == "dot":
+		# Find last entry with a dot i in it
+		for i, tok in reversed(list(enumerate(toks))):
+			if "." in tok: break
+		else: raise ValueError("Could not split hdf path using 'dot' method: no . found")
+	elif mode == "exists":
+		for i in reversed(list(range(len(toks)))):
+			cand = "/".join(toks[:i+1])
+			if os.path.isfile(cand): break
+		else: raise ValueError("Could not split hdf path using 'exists' method: no file found")
+	elif mode == "none":
+		i = len(toks)
+	else: raise ValueError("Unrecognized split mode '%s'" % (str(mode)))
+	# Return the result
+	fname = "/".join(toks[:i+1])
+	gtoks = toks[i+1:]
+	if subgroup: gtoks.append(subgroup)
+	group = "/".join(gtoks) if len(gtoks)>0 else None
+	return fname, group
 
 def concatenate(bunches):
 	"""Go from a list of bunches to a bunch of lists."""

pixell/colorize.py

@@ -169,7 +169,7 @@ def mpl_register(names=None):
 	if isinstance(names, basestring): names = [names]
 	for name in names:
 		cmap = to_mpl_colormap(name, schemes[name])
-		matplotlib.cm.register_cmap(name, cmap)
+		matplotlib.colormaps.register(cmap)
 
 def mpl_setdefault(name):
 	import matplotlib.pyplot

pixell/curvedsky.py

@@ -1238,7 +1238,7 @@ def analyse_geometry(shape, wcs, tol=1e-6):
 	# TODO: Pseudo-cylindrical projections can be handled with standard ducc synthesis,
 	# so ideally our check would be less stringent than this. Supporinting them requires
 	# more work, so will just do it with the general interface for now.
-	separable = wcsutils.is_cyl(wcs)
+	separable = wcsutils.is_separable(wcs)
 	divides   = utils.hasoff(360/np.abs(wcs.wcs.cdelt[0]), 0, tol=tol)
 	if not separable or not divides:
 		# Not cylindrical or ra does not evenly divide the sky