Implementing buildDescriptor from vlsvvtkreader as a C function

.github/workflows/test_python_turso.yml

@@ -31,7 +31,7 @@ jobs:
           timeout-minutes: 5
           run: |
             export TMPDIR=$RUNNER_TEMP
-            uv venv CI_env
+            uv venv -p 3.10 CI_env
             . CI_env/bin/activate
             uv pip install --editable ../analysator[${{ matrix.extras }}]
         - name: Trial imports

analysator/cpphelpers/cpphelpers.cpp

@@ -0,0 +1,170 @@
+#define PY_SSIZE_T_CLEAN
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+#include <Python.h>
+#include <cstdint>
+#include <iostream>
+#include <numpy/arrayobject.h>
+#include <numpy/ndarraytypes.h>
+#include <sstream>
+#include <stdint.h>
+#include <unordered_map>
+#include <vector>
+using namespace std;
+
+// C's % operator is "remainder" operator not modulus like Python's % (not sure
+// if there is implementation of this in stnadard library, but its not big so
+// here it is)
+constexpr int mod(int a, int b) noexcept { return ((a % b) + b) % b; }
+constexpr int CHILDS = 8;
+constexpr int floordiv(int a, int b) noexcept {
+  int q = a / b;
+  int r = a % b;
+  if ((r != 0) && ((r > 0) != (b > 0))) {
+    q--;
+  }
+  return q;
+}
+
+//Convert unordered_map into a Python dictionary
+static PyObject *convertToDict(unordered_map<int, uint64_t> &map) {
+  PyObject *dict = PyDict_New();
+  for (const auto &it : map) { //structured binding?
+    PyObject *key = PyLong_FromLongLong(it.first);
+    PyObject *val = PyLong_FromLongLong(it.second);
+    PyDict_SetItem(dict, key, val); //error code handling?
+    Py_DECREF(key);
+    Py_DECREF(val);
+  }
+  return dict;
+}
+
+//could be changed to span?
+static void children(const int cid, const int level, vector<int64_t> &cid_offsets,
+                     vector<int64_t> &xcells, vector<int64_t> &ycells,
+                     vector<int64_t> &zcells, vector<int64_t> &out,
+                     const vector<vector<int32_t>>& delta
+                     ) {
+  long cellid = cid - 1 - cid_offsets[level];
+  vector<int32_t> cellind(3, -1); //could be reused
+  cellind[0] = mod(cellid, (xcells[level])) * 2;
+  cellind[1] = mod(floordiv(cellid, xcells[level]), (ycells[level])) * 2;
+  cellind[2] = floordiv(cellid, xcells[level] * ycells[level]) * 2;
+  //children vector always size 8 (hopefully)
+  for (size_t i = 0; i < CHILDS; ++i) {
+    out[i] =
+        cid_offsets[level + 1] + (cellind[0] + delta[i][0]) +
+        xcells[level + 1] * (cellind[1] + delta[i][1]) +
+        (cellind[2] + delta[i][2]) * xcells[level + 1] * ycells[level + 1] + 1;
+  }
+
+  // return out;
+}
+
+//Convert python dictionary to unordered_map
+static int convertToUnordMap(PyObject *dict,
+                             unordered_map<int, uint64_t> &map) {
+  PyObject *key, *value;
+
+  Py_ssize_t pos = 0;
+
+  while (PyDict_Next(dict, &pos, &key, &value)) {
+
+    uint64_t val = PyLong_AsUnsignedLongLong(value);
+    long keyval = PyLong_AsLong(key);
+    map[keyval] = val;
+    if (PyErr_Occurred()) {
+      return 1;
+    }
+  }
+  return 0;
+}
+
+static PyObject *pyBuildDescriptor(PyObject *self, PyObject *args) {
+  unsigned int max_ref_level;
+  PyObject *fileindex_for_cellid;
+  unsigned int xc, yc, zc;
+
+  stringstream descr; //would vector with push_back be faster?
+  // O!|OO (1 required arg (PythonObject) with 2 optional (not sure why we need
+  // ! on the first one))
+  //  more args O!|O|i for integer
+  if (!PyArg_ParseTuple(args, "O|i|i|i|i", &fileindex_for_cellid, &xc, &yc, &zc,
+                        &max_ref_level)) {
+    return NULL;
+  }
+
+  vector<int64_t> xcells(max_ref_level + 1, 0);
+  vector<int64_t> ycells(max_ref_level + 1, 0);
+  vector<int64_t> zcells(max_ref_level + 1, 0);
+  //bitshift very cool
+  for (size_t r = 0; r < max_ref_level + 1; r++) {
+    xcells[r] = xc << r;
+    ycells[r] = yc << r;
+    zcells[r] = zc << r;
+  }
+  unordered_map<int, uint64_t> idxToFileIndex;
+  unordered_map<int, uint64_t> fileindex_for_cellid_map;
+  if (convertToUnordMap(fileindex_for_cellid, fileindex_for_cellid_map) != 0) {
+    cout << "error" << endl;
+  }
+  int idx = 0;
+  vector<vector<int>> subdivided(max_ref_level + 1);
+  vector<int64_t> cid_offsets(max_ref_level + 1, 0);
+  uint64_t isum = 0;
+  for (size_t p = 0; p < max_ref_level; p++) {
+    isum = isum + pow(2, 3 * p) * xc * yc * zc;
+    cid_offsets[p + 1] = isum;
+  }
+  for (size_t c = 1; c < xc * yc * zc + 1; c++) {
+    if (fileindex_for_cellid_map.find(c) != fileindex_for_cellid_map.end()) {
+      // Write
+      descr.put('.');
+      idxToFileIndex[idx] = fileindex_for_cellid_map[c];
+    } else {
+      descr.put('R');
+      subdivided[0].push_back(c);
+    }
+    idx = idx + 1;
+  }
+  static const vector<vector<int32_t>> delta = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0},
+                                   {0, 0, 1}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1}}; //array preferable?
+
+  vector<int64_t> childs(8, 0); //array preferabl?
+  descr.put('|');
+  for (size_t l = 1; l < max_ref_level + 1; l++) {
+    auto &vecptr = subdivided[l - 1];
+    auto &subd = subdivided[l];
+    for (int it : vecptr) {
+      children(it, l - 1, cid_offsets, xcells, ycells, zcells, childs, delta);
+      for (int64_t child : childs) {
+        auto it2 = fileindex_for_cellid_map.find(child);
+        if (it2 != fileindex_for_cellid_map.end()) {
+          descr.put('.');
+          idxToFileIndex[idx] = it2->second;
+        } else {
+          descr.put('R');
+          subd.push_back(child);
+        }
+        idx = idx + 1;
+      }
+    }
+    if (l < max_ref_level) {
+      descr.put('|');
+    }
+  }
+  PyObject *outdict = convertToDict(idxToFileIndex);
+  return Py_BuildValue("sO", descr.str().data(), outdict);
+}
+
+static PyMethodDef cpphelpers_methods[] = {
+    {"buildDescriptor", pyBuildDescriptor, METH_VARARGS, "Build descriptor"}, {NULL} /* Sentinel */
+};
+static struct PyModuleDef cpphelpers = {
+    .m_methods = cpphelpers_methods,
+};
+PyMODINIT_FUNC PyInit_cpphelpers(void)
+// create the module
+{
+  import_array();
+  return PyModuleDef_Init(&cpphelpers);
+}

analysator/cpphelpers/pyproject.toml

@@ -0,0 +1,8 @@
+[build-system]
+requires = ["setuptools>=61","numpy"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "cpphelpers"
+version = "1.0"
+

analysator/cpphelpers/setup.py

@@ -0,0 +1,8 @@
+from setuptools import setup, Extension
+import numpy
+setup(
+    name="cpphelpers",
+    version="1.0",
+    ext_modules=[Extension("cpphelpers", sources=["cpphelpers.cpp"])],
+    include_dirs=[numpy.get_include()]
+)

analysator/cpphelpers/test.py

@@ -0,0 +1,95 @@
+import cpphelpers as chp
+import analysator as pt
+import numpy as np
+import time
+def buildDescriptor(f):
+    now = time.time() 
+    f._VlsvReader__read_fileindex_for_cellid()
+    print("read",time.time()-now)
+    fileindex_for_cellid = f._VlsvReader__fileindex_for_cellid
+    xc = f._VlsvReader__xcells
+    yc = f._VlsvReader__ycells
+    zc = f._VlsvReader__zcells
+    max_ref_level = f.get_max_refinement_level()
+
+    cid_offsets = np.zeros(max_ref_level+1, dtype=np.int64)
+    isum = 0
+    for i in range(0,max_ref_level):
+        isum = isum + 2**(3*i) * xc * yc * zc
+        cid_offsets[i+1] = isum
+    xcells = np.zeros((max_ref_level+1), dtype=np.int64)
+    ycells = np.zeros((max_ref_level+1), dtype=np.int64)
+    zcells = np.zeros((max_ref_level+1), dtype=np.int64)
+    for r in range(max_ref_level+1):
+        xcells[r] = xc*2**(r)
+        ycells[r] = yc*2**(r)
+        zcells[r] = zc*2**(r)
+    delta = np.array([[0,0,0],[1,0,0],[0,1,0],[1,1,0],
+            [0,0,1],[1,0,1],[0,1,1],[1,1,1]],dtype=np.int32)
+
+
+    idxToFileIndex = {}
+    from io import StringIO
+    descr = StringIO()
+    print("Building descriptor")
+    subdivided = [[] for l in range(max_ref_level+1)]
+    idx = 0
+    for c in range(1,int(np.prod(f.get_spatial_mesh_size()))+1):
+        if c in fileindex_for_cellid.keys():
+            descr.write(".")
+            idxToFileIndex[idx] = fileindex_for_cellid[c]
+        else:
+            descr.write("R")
+            subdivided[0].append(c)
+        idx += 1
+    # @jit(nopython=True)
+    def children(cid, level):
+        # cellind = get_cell_indices(cid, level) # get_cellind here for compilation
+        cellids = cid - 1 - cid_offsets[level]
+        cellind = np.full(3, -1,dtype=np.int32)
+        cellind[0] = (cellids)%(xcells[level])
+        cellind[1] = ((cellids)//(xcells[level]))%(ycells[level])
+        cellind[2] = (cellids)//(xcells[level]*ycells[level])
+        cellind = cellind*2
+
+
+        out = np.zeros((8,),dtype=np.int64)
+        out[:] = cid_offsets[level+1] + (cellind[0] + delta[:,0]) + xcells[level+1]*(cellind[1] + delta[:,1]) + (cellind[2] + delta[:,2])*xcells[level+1]*ycells[level+1] + 1
+        return out           
+    descr.write("|")
+    for l in range(1,max_ref_level+1):
+        for c in subdivided[l-1]:
+            for child in children(c, l-1):
+                if child in fileindex_for_cellid.keys():
+                    descr.write(".")
+                    idxToFileIndex[idx] = fileindex_for_cellid[child]
+                else:
+                    descr.write("R")
+                    subdivided[l].append(child)
+                idx += 1
+        if l < max_ref_level:
+            descr.write("|")
+    return descr.getvalue(), idxToFileIndex
+
+file="/wrk-vakka/group/spacephysics/vlasiator/3D/FID/bulk1/bulk1.0001213.vlsv"
+f=pt.vlsvfile.VlsvReader(file)
+
+now=time.time()
+var=buildDescriptor(f)
+print(time.time()-now)
+
+now=time.time()
+
+f._VlsvReader__read_fileindex_for_cellid()
+print(time.time()-now)
+fileindex_for_cellid = f._VlsvReader__fileindex_for_cellid
+xc= f._VlsvReader__xcells
+yc= f._VlsvReader__ycells
+zc= f._VlsvReader__zcells
+max_ref_level = f.get_max_refinement_level()     
+#print(f.get_spatial_mesh_size()) just [xc,yc,zc]
+test=chp.test(fileindex_for_cellid,xc,yc,zc,max_ref_level)
+print(time.time()-now)
+print(test[0]==var[0])
+print(test[1]==var[1])
+

analysator/pyVlsv/vlsvvtkinterface.py

@@ -28,7 +28,7 @@
 import pickle
 from operator import itemgetter
 import numbers
-
+from warnings import deprecated
 import analysator as pt
 try:
    import vtk
@@ -643,7 +643,8 @@ def SetFileName(self, filename):
    def GetFileName(self):
       return self.__FileName
 
-   def buildDescriptor(self):
+   @deprecated("The 'buildDescriptor' should be used instead as it is faster, this may be removed in the future but is still here in case")
+   def buildDescriptorPython(self):
       f = self.__reader
       f._VlsvReader__read_fileindex_for_cellid()
       fileindex_for_cellid = f._VlsvReader__fileindex_for_cellid
@@ -715,6 +716,20 @@ def children(cid, level):
             descr.write("|")
 
       return descr.getvalue(), idxToFileIndex
+
+   def buildDescriptor(self):
+        import cpphelpers
+        f = self.__reader
+        f._VlsvReader__read_fileindex_for_cellid()
+        fileindex_for_cellid = f._VlsvReader__fileindex_for_cellid
+        xc= f._VlsvReader__xcells
+        yc= f._VlsvReader__ycells
+        zc= f._VlsvReader__zcells
+        max_ref_level = f.get_max_refinement_level()     
+
+        ret=cpphelpers.buildDescriptor(fileindex_for_cellid,xc,yc,zc,max_ref_level)
+        #returns descriptor as first value and second is the idxToFileIndex
+        return ret[0],ret[1]
 
 
    def getDescriptor(self, reinit=False):

pyproject.toml

@@ -1,5 +1,5 @@
 [build-system]
-requires = ["hatchling"]
+requires = ["hatchling","numpy"]
 build-backend = "hatchling.build"
 
 [project]
@@ -9,7 +9,7 @@ authors = [
     { name = "Markku Alho", email = "markku.alho@helsinki.fi"},
     { name = "Et al"},
 ]
-maintainters = [
+maintainers = [
     { name = "Markku Alho", email = "markku.alho@helsinki.fi"},
 ]
 description = "Tools for reading and analysing Vlasiator .vlsv output files."
@@ -28,8 +28,10 @@ dependencies = [
     "scikit-image",
     "vtk>=9.2",
     "yt",
-    "rtree"
+    "rtree",
+    "cpphelpers @ {root:uri}/analysator/cpphelpers/"
 ]
+
 [project.optional-dependencies]
 none = [
 ]
@@ -50,3 +52,8 @@ packages = ["./analysator","./pytools"]
 
 [tool.hatch.build.targets.sdist]
 packages = ["./analysator","./pytools"]
+
+[tool.hatch.metadata]
+allow-direct-references = true
+
+