Merge branch 'master' of https://github.com/ahay/src

Author: sfomel

.github/workflows/dockerimage.yml

@@ -13,4 +13,4 @@ jobs:
         name: madagascar/m8r
         username: ${{ secrets.DOCKER_USERNAME }}
         password: ${{ secrets.DOCKER_PASSWORD }}
-        dockerfile: admin/docker/2.0-dev-tex/Dockerfile
+        dockerfile: admin/docker/4.2-dev-tex/Dockerfile

admin/docker/4.2-dev-tex/Dockerfile

@@ -39,7 +39,7 @@ RUN apt-get update && apt-get install -y \
  && rm -rf /var/lib/apt/lists/*
 
 # install python packages
-RUN pip install numpy scipy
+RUN pip install --break-system-packages numpy scipy
 
 # get code from github
 RUN git clone https://github.com/ahay/src.git $HOME/RSFSRC
@@ -59,7 +59,6 @@ RUN apt-get update && apt-get install -y \
 		texlive-fonts-recommended \
 		texlive-bibtex-extra \
 		texlive-lang-english \
-		texlive-generic-extra \
 		biber \
         --no-install-recommends \
  && rm -rf /var/lib/apt/lists/*

api/python/mxarray.py

@@ -16,8 +16,10 @@
 ##   along with this program; if not, write to the Free Software
 ##   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
+import os
 import m8r
 import numpy as np
+import subprocess
 
 try:
     import xarray as xr
@@ -31,17 +33,15 @@
 
 def rsf_to_xarray(path, chunks = "auto"):
     """Convert an RSF file to an xarray DataArray with optimal Dask lazy loading."""
-    
     if xr is None:
         raise ImportError("xarray is required.")
 
     f = m8r.Input(path)
-    
     shape = f.shape()
     ndim = len(shape)
 
     dims = []
-    # units = []
+    units = []
     coords = {}
     for axis in range(1, ndim+1):
         n = f.int(f"n{axis}")
@@ -50,21 +50,35 @@ def rsf_to_xarray(path, chunks = "auto"):
         label = f.string(f"label{axis}")
         if label == None:
             label = f"dim{axis}"
-        # unit = f.string(f"unit{axis}")
+        unit = f.string(f"unit{axis}")
+        if unit == None:
+            unit = "??"
 
         coords[label] = np.arange(n) * d + o
         dims.append(label)
-        # units.append(unit)
+        units.append(unit)
 
     # data = np.asarray(f)
     binFile = f.string("in")
-    dtype = f.string("type")
-    if dtype is None or dtype == "float":
-        dtype = np.float32
-    elif dtype == "int":
-        dtype = np.int32
-    else:
-        raise ValueError(f"Unsupported data type: {dtype}")
+    label = f.string("label")
+    unit0 = f.string("unit")
+
+    rsf_type = getattr(f, 'type', None)
+    if not rsf_type:
+        fmt = f.string("data_format")
+        if fmt and "complex" in fmt: rsf_type = 'complex'
+        elif fmt and "int" in fmt: rsf_type = 'int'
+        else: rsf_type = 'float'
+    
+    dtype_map = {
+        'float':   np.float32,
+        'int':     np.int32,
+        'complex': np.complex64,
+        'uchar':   np.uint8,
+        'char':    np.int8
+    }
+    dtype = dtype_map.get(rsf_type, np.float32)
+    # -------------------------------
 
     mm = np.memmap(
         binFile,
@@ -80,33 +94,61 @@ def rsf_to_xarray(path, chunks = "auto"):
 
 
     # covert c order to python order
-    # data = data.reshape(shape_py)
     data = data.transpose(*reversed(range(data.ndim)))
 
     ds = xr.DataArray(
         data, 
         dims=dims, 
-        coords=coords
+        coords=coords,
     )
+    if label:
+        ds.attrs['long_name'] = label
+    if unit0:
+        ds.attrs['units'] = unit0
+        
+    for dim, unit in zip(dims, units):
+        ds.coords[dim].attrs['units'] = unit
 
     return ds
 
 def xarray_to_rsf(ds, outpath):
     """Convert an xarray Dataset to an RSF file."""
-    
     if xr is None:
         raise ImportError("xarray is required.")
 
     if not isinstance(ds, xr.DataArray):
         raise ValueError("Input must be an xarray DataArray.")
 
+    if np.issubdtype(ds.dtype, np.complexfloating):
+        rsf_type = 'complex'
+        out_dtype = np.complex64
+        fmt_str = "native_complex"
+    elif np.issubdtype(ds.dtype, np.integer):
+        rsf_type = 'int'
+        out_dtype = np.int32
+        fmt_str = "native_int"
+    else:
+        rsf_type = 'float' 
+        out_dtype = np.float32
+        fmt_str = "native_float"
+    # -------------------------------
+
     data = ds.values
     data = data.transpose(*reversed(range(data.ndim)))
 
     dims = ds.dims
+    out = m8r.Output(outpath)
+    label = ds.attrs.get('long_name', None)
+    unit0 = ds.attrs.get('units', None)
 
+    # Set Type
+    out.settype(rsf_type)
+    out.put("data_format", fmt_str)
 
-    out = m8r.Output(outpath)
+    if label:
+        out.put("label", label)
+    if unit0:
+        out.put("unit", unit0)
 
     for i, dim in enumerate(dims, start=1):
 
@@ -123,8 +165,9 @@ def xarray_to_rsf(ds, outpath):
         out.put(f"o{i}", o)
         out.put(f"d{i}", d)
         out.put(f"label{i}", str(dim))
-        
-    out.write(data.astype(np.float32))
+        out.put(f"unit{i}", str(ds.coords[dim].attrs.get('units', '??')))
+
+    out.write(data.astype(out_dtype))
     out.close()
 
 def rsf_to_xarrayds(path, chunks = "auto"):
@@ -133,4 +176,121 @@ def rsf_to_xarrayds(path, chunks = "auto"):
     da = rsf_to_xarray(path, chunks=chunks)
     ds = da.to_dataset(name="data")
 
-    return ds
+    return ds
+
+## Monkey patching m8r.Filter to handle xarray inputs/outputs
+
+def _patched_setcommand(self, kw, args=[]):
+    """
+    Patched version of Filter.setcommand to handle auxiliary xarrays.
+    Example: Filter('sfprog')(velocity=my_xarray)
+    """
+    if not hasattr(self, '_mx_aux_refs'):
+        self._mx_aux_refs = []
+
+    for key, val in list(kw.items()):
+        if isinstance(val, xr.DataArray):
+            tmp_name = m8r.Temp()
+            xarray_to_rsf(val, tmp_name)
+
+            f_obj = m8r.File(tmp_name, temp=True)
+            self._mx_aux_refs.append(f_obj)
+            
+            kw[key] = str(tmp_name)
+
+    new_args = []
+    for val in args:
+        if isinstance(val, xr.DataArray):
+            tmp_name = m8r.Temp()
+            xarray_to_rsf(val, tmp_name)
+            f_obj = m8r.File(tmp_name, temp=True)
+            self._mx_aux_refs.append(f_obj)
+            new_args.append(str(tmp_name))
+        else:
+            new_args.append(val)
+            
+    # Call the original function with strings/files only
+    original_func = getattr(self.__class__, '_original_setcommand_func', None)
+    if not original_func:
+        # Fallback if class attribute missing
+        real_module = m8r.wrapped if hasattr(m8r, 'wrapped') else m8r
+        original_func = real_module.Filter._original_setcommand_func
+        
+    return original_func(self, kw, new_args)
+
+
+def _patched_apply(self, *srcs):
+    """
+    Patched version of Filter.apply to handle xarray inputs/outputs.
+    """
+    
+    # Check if any input is an xarray
+    if not any(isinstance(s, xr.DataArray) for s in srcs):
+        # Fallback to original
+        if hasattr(self.__class__, '_original_apply_func'):
+            return self.__class__._original_apply_func(self, *srcs)
+        real_module = m8r.wrapped if hasattr(m8r, 'wrapped') else m8r
+        return real_module.Filter._original_apply_func(self, *srcs)
+
+    # Handle xarray Input
+    clean = []
+    rsf_inputs = []
+    try:
+        for s in srcs:
+            if isinstance(s, xr.DataArray):
+                tmp = m8r.Temp()
+                xarray_to_rsf(s, tmp)
+                rsf_inputs.append(str(tmp))
+                clean.extend([str(tmp), str(tmp)+'@'])
+            else:
+                rsf_inputs.append(str(s))
+
+        out_file = m8r.Temp()
+        
+        # Split command for correct pipe handling
+        first, pipe_char, rest = self.command.partition('|')
+        cmd_parts = [first]
+        
+        if len(rsf_inputs) > 0:
+            cmd_parts.append(f"< {rsf_inputs[0]}")
+            if len(rsf_inputs) > 1:
+                cmd_parts.extend(rsf_inputs[1:])
+        
+        if pipe_char:
+            cmd_parts.append(pipe_char)
+            cmd_parts.append(rest)
+            
+        cmd_parts.append(f"> {out_file}")
+        
+        full_cmd = " ".join(cmd_parts)
+        
+        result = subprocess.run(full_cmd, shell=True, stderr=subprocess.PIPE, stdout=subprocess.PIPE, universal_newlines=True)
+        
+        if result.returncode != 0:
+            raise RuntimeError(f"{full_cmd}\n{result.stderr}")
+
+        if self.plot: 
+            return m8r.Vplot(out_file, temp=True)
+            
+        res = rsf_to_xarray(out_file)
+        clean.extend([str(out_file), str(out_file)+'@'])
+        return res
+
+    finally:
+        for f in clean:
+            if os.path.exists(f): 
+                try: os.unlink(f)
+                except: pass
+                
+if hasattr(m8r, 'wrapped'):
+    real_m8r_module = m8r.wrapped
+else:
+    real_m8r_module = m8r
+
+if not hasattr(real_m8r_module.Filter, '_original_apply_func'):
+    real_m8r_module.Filter._original_apply_func = real_m8r_module.Filter.apply
+    real_m8r_module.Filter.apply = _patched_apply
+
+if not hasattr(real_m8r_module.Filter, '_original_setcommand_func'):
+    real_m8r_module.Filter._original_setcommand_func = real_m8r_module.Filter.setcommand
+    real_m8r_module.Filter.setcommand = _patched_setcommand