From c2c76f626235e2a2cbe3c55e1709b7badcbef7df Mon Sep 17 00:00:00 2001
From: Jonathan Wright <wright@esrf.fr>
Date: Thu, 21 Mar 2024 16:56:52 +0100
Subject: [PATCH] limit threads and do not overwrite geometry

---
 ImageD11/sinograms/point_by_point.py | 34 ++++++++++++++++++----------
 setup.py                             |  1 +
 2 files changed, 23 insertions(+), 12 deletions(-)

diff --git a/ImageD11/sinograms/point_by_point.py b/ImageD11/sinograms/point_by_point.py
index 115863fc..4d3feb9a 100755
--- a/ImageD11/sinograms/point_by_point.py
+++ b/ImageD11/sinograms/point_by_point.py
@@ -4,11 +4,14 @@
 # Try to build the point-by-point mapping code ...
 from __future__ import print_function, division
 
-import multiprocessing
-# method = multiprocessing.set_start_method("forkserver")
-from ImageD11 import cImageD11  # NO!!! patches to forkserver
-
+from ImageD11 import cImageD11
+cImageD11.check_multiprocessing(patch=True) # forkserver
+try:
+    import threadpoolctl
+except ImportError:
+    threadpoolctl = None
 import sys, os
+import multiprocessing
 from multiprocessing.managers import SharedMemoryManager
 import time, random
 import numpy as np
@@ -63,6 +66,8 @@ def dtymask(i, j, cosomega, sinomega, dtyi, step, y0):
     return dtyi == dty_icalc
 
 
+
+
 def idxpoint(
     i,
     j,
@@ -194,6 +199,8 @@ def proxy(args):
 
 def initializer(parfile, symmetry, loglevel=3):
     global ucglobal, symglobal, parglobal
+    if threadpoolctl is not None:
+        threadpoolctl.threadpool_limits(limits=1)
     parglobal = parameters.read_par_file(parfile)
     ucglobal = unitcell.unitcell_from_parameters(parglobal)
     symglobal = sym_u.getgroup(symmetry)()
@@ -255,7 +262,8 @@ def setpeaks(self, colf):
         """
         # Load the peaks
         colf.parameters.loadparameters(self.parfile)
-        colf.updateGeometry()  # for ds
+        if 'ds' not in colf.titles:
+            colf.updateGeometry()  # for ds
         #
         uc = unitcell.unitcell_from_parameters(colf.parameters)
         uc.makerings(colf.ds.max())
@@ -269,9 +277,11 @@ def setpeaks(self, colf):
         hmax = 0
         for i in range(len(uc.ringds)):
             ds = uc.ringds[i]
-            rm = abs(colf.ds - ds) < self.ds_tol
             hkls = uc.ringhkls[ds]
-            if i in self.foridx and rm.sum() > 0:
+            if i in self.foridx:
+                rm = abs(colf.ds - ds) < self.ds_tol
+                if rm.sum() == 0:
+                    continue
                 icut = np.max(colf.sum_intensity[rm]) * self.ifrac
                 rm = rm & (colf.sum_intensity > icut)
                 isel |= rm
@@ -287,7 +297,7 @@ def setpeaks(self, colf):
                 )
                 sel |= rm
             else:
-                print(i, "%.4f" % (ds), hkls[-1], len(hkls), rm.sum(), "skip")
+                print(i, "%.4f" % (ds), hkls[-1], len(hkls), "skipped")
             for hkl in hkls:
                 hmax = max(np.abs(hkl).max(), hmax)
         if self.forgen is None:
@@ -318,13 +328,13 @@ def setpeaks(self, colf):
         )
         self.hmax = hmax
 
-    def iplot(self):
+    def iplot(self, skip=1):
         import pylab as pl
 
         f, ax = pl.subplots(2, 1)
-        ax[0].plot(self.colf.tth, self.colf.sum_intensity, ",")
-        ax[0].plot(self.icolf.tth, self.icolf.sum_intensity, ",")
-        ax[0].set(yscale="log", ylabel="sum intensity", xlabel="tth")
+        ax[0].plot(self.colf.ds[::skip], self.colf.sum_intensity[::skip], ",")
+        ax[0].plot(self.icolf.ds[::skip], self.icolf.sum_intensity[::skip], ",")
+        ax[0].set(yscale="log", ylabel="sum intensity", xlabel="d-star")
         histo = self.dset.sinohist(
             omega=self.icolf.omega, dty=self.icolf.dty, weights=self.icolf.sum_intensity
         )
diff --git a/setup.py b/setup.py
index 8df47dc0..70f440fc 100644
--- a/setup.py
+++ b/setup.py
@@ -166,6 +166,7 @@ def build_extension(self, ext):
     "fast_histogram",
     "scikit-image",
     "tqdm",
+    'threadpoolctl ; python_version >= "3" ',
 ]