Fix for issue 132

sacc/covariance.py

@@ -1,4 +1,4 @@
-from astropy.table import Table
+from astropy.table import Table, Column
 import scipy.linalg
 import numpy as np
 import warnings
@@ -218,9 +218,9 @@ def to_table(self):
         table: astropy.table.Table instance
             Table that can be used to reconstruct the object.
         """
-        col_names = [f'col_{i}' for i in range(self.size)]
-        cols = [self.covmat[i] for i in range(self.size)]
-        table = Table(data=cols, names=col_names)
+        # Store as a single vector column ('row') to avoid FITS TFIELDS>999
+        # Each table row is one row of the covariance matrix
+        table = Table([Column(name='row', data=self.covmat)])
         table.meta['SIZE'] = self.size
         return table
 
@@ -239,7 +239,11 @@ def from_table(cls, table):
             Loaded covariance object
         """
         size = table.meta['SIZE']
-        covmat = np.array([table[f'col_{i}'] for i in range(size)])
+        # Support both legacy many-column format and new single-column format
+        if 'row' in table.colnames:
+            covmat = np.array(list(table['row']))
+        else:
+            covmat = np.array([table[f'col_{i}'] for i in range(size)])
         return cls(covmat)
 
 
@@ -375,8 +379,12 @@ def from_tables(cls, tables):
         for i in range(nblock):
             table = tables[f'block_{i}']
             block_size = table.meta['SACCBSZE']
-            cols = [table[f'block_col_{i}'] for i in range(block_size)]
-            blocks.append(np.array(cols))
+            if 'block_row' in table.colnames:
+                block = np.array(list(table['block_row']))
+            else:
+                cols = [table[f'block_col_{j}'] for j in range(block_size)]
+                block = np.array(cols)
+            blocks.append(block)
         return cls(blocks)
 
     def to_tables(self):
@@ -396,9 +404,8 @@ def to_tables(self):
         nblock = len(self.blocks)
         for j, block in enumerate(self.blocks):
             b = len(block)
-            col_names = [f'block_col_{i}' for i in range(b)]
-            cols = [block[i] for i in range(b)]
-            table = Table(data=cols, names=col_names)
+            # Use single vector column to minimize TFIELDS
+            table = Table([Column(name='block_row', data=block)])
             table.meta['SIZE'] = self.size
             table.meta['SACCBIDX'] = j
             table.meta['SACCBCNT'] = nblock

sacc/sacc.py

@@ -5,6 +5,10 @@
 
 from astropy.io import fits
 from astropy.table import Table
+
+# Module-level constants for Sacc file format versions
+SACCFVER = 2  # Current FITS version
+SACCHDF5VER = 1  # Current HDF5 version
 import numpy as np
 
 from .tracers import BaseTracer
@@ -940,10 +944,10 @@ def save_fits(self, filename, overwrite=False):
 
         # Create the actual fits object
         primary_header = fits.Header()
+        primary_header['SACCFVER'] = SACCFVER
         with warnings.catch_warnings():
             warnings.simplefilter("ignore", category=fits.verify.VerifyWarning)
-            hdus = [fits.PrimaryHDU(header=primary_header)] + \
-                    [fits.table_to_hdu(table) for table in tables]
+            hdus = [fits.PrimaryHDU(header=primary_header)] + [fits.table_to_hdu(table) for table in tables]
         hdu_list = fits.HDUList(hdus)
         io.astropy_buffered_fits_write(filename, hdu_list)
 
@@ -962,36 +966,34 @@ def load_fits(cls, filename):
         """
         cov = None
         metadata = None
+        fitsver = None
 
         with fits.open(filename, mode="readonly") as f:
             tables = []
-            for hdu in f:
+            for idx, hdu in enumerate(f):
                 if hdu.name.lower() == 'primary':
-                    # The primary table is not a data table,
-                    # but in older files it was used to store metadata
                     header = hdu.header
+                    fitsver = header.get('SACCFVER', None)
+                    if fitsver is None:
+                        fitsver = 1
+                    if fitsver > SACCFVER:
+                        raise RuntimeError(f"Unsupported SACC FITS version: {fitsver}")
                     if "NMETA" in header:
                         metadata = {}
-                        # Older format metadata is kept in the primary
-                        # header, with keys KEY0, VAL0, etc.
                         n_meta = header['NMETA']
                         for i in range(n_meta):
                             k = header[f'KEY{i}']
                             v = header[f'VAL{i}']
                             metadata[k] = v
                 elif hdu.name.lower() == 'covariance':
-                    # Legacy covariance - HDU will just be called covariance
-                    # instead of the full name given by BaseIO.
-                    # Note that this will also allow us to use multiple
-                    # covariances in future.
                     cov = BaseCovariance.from_hdu(hdu)
                 else:
                     tables.append(Table.read(hdu))
 
-        # add the metadata table, if we are in the legacy format
         if metadata is not None:
             tables.append(io.metadata_to_table(metadata))
 
+        # Pass version to from_tables if needed (future-proofing)
         return cls.from_tables(tables, cov=cov)
 
     def save_hdf5(self, filename, overwrite=False, compression='gzip', compression_opts=4):
@@ -1027,6 +1029,8 @@ def save_hdf5(self, filename, overwrite=False, compression='gzip', compression_o
                 table.meta['EXTNAME'] = extname
 
         with h5py.File(filename, 'w') as f:
+            # Write version dataset
+            f.create_dataset('sacc_hdf5_version', data=np.array([SACCHDF5VER], dtype='i4'))
             used_names = {}
             for table in tables:
                 # Build a meaningful dataset name
@@ -1087,9 +1091,19 @@ def load_hdf5(cls, filename):
         """
         import h5py
         recovered_tables = []
+        hdf5ver = None
         with h5py.File(filename, 'r') as f:
+            # Check version
+            if 'sacc_hdf5_version' in f:
+                hdf5ver = int(np.array(f['sacc_hdf5_version'])[0])
+            else:
+                hdf5ver = 1
+            if hdf5ver > SACCHDF5VER:
+                raise RuntimeError(f"Unsupported SACC HDF5 version: {hdf5ver}")
             # Read all datasets (not groups) in the order they appear
             for key in f.keys():
+                if key == 'sacc_hdf5_version':
+                    continue
                 item = f[key]
                 if isinstance(item, h5py.Dataset):
                     table = Table.read(f, path=key)

test/test_bug_132.py

@@ -0,0 +1,36 @@
+import os
+import tempfile
+
+import numpy as np
+
+import sacc
+
+
+def test_bug_132():
+    s = sacc.Sacc()
+
+    s.add_tracer(
+        "Map",
+        "tracer",
+        quantity="galaxy_density",
+        spin=0,
+        ell=np.arange(100),
+        beam=np.ones(100),
+    )
+
+    ndata = 1000
+    s.add_ell_cl("cl_00", "tracer", "tracer", np.arange(ndata), np.ones(ndata))
+    # Add covariance
+    s.add_covariance(np.eye(ndata))
+
+    # Write to SACC files using a temporary directory
+    with tempfile.TemporaryDirectory() as tmpdir:
+        fits_filename = os.path.join(tmpdir, "test.fits")
+        hdf5_filename = os.path.join(tmpdir, "test.hdf5")
+        s.save_fits(fits_filename, overwrite=True)
+        s.save_hdf5(hdf5_filename, overwrite=True)
+        # Read back in
+        s2 = sacc.Sacc.load_fits(fits_filename)
+        s3 = sacc.Sacc.load_hdf5(hdf5_filename)
+        assert s2 == s
+        assert s3 == s