NLPAR GPU (#56)

* GPU enabled version of NLPAR
* PyEBSDIndex will now default to using discrete graphics if both integrated and discrete GPUs are found.

.github/dependabot.yml

@@ -5,7 +5,10 @@
 
 version: 2
 updates:
-  - package-ecosystem: "" # See documentation for possible values
+  - package-ecosystem: "github-actions" # See documentation for possible values
     directory: "/" # Location of package manifests
+    target-branch: "develop"
     schedule:
       interval: "weekly"
+    open-pull-requests-limit: 2
+

.github/workflows/tests.yml

@@ -49,10 +49,10 @@ jobs:
             DEPENDENCIES: matplotlib==3.3 numba==0.52 numpy==1.19 ray[default]==1.13
             LABEL: -oldest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python-version }}
 

License

@@ -18,5 +18,5 @@ works bear some notice that they are derived from it, and any modified versions
 some notice that they have been modified.
 
 Author: David Rowenhorst; 
-The US Naval Research Laboratory Date: 21 Aug 2020
+The US Naval Research Laboratory Date: 22 May 2024
 

pyebsdindex/__init__.py

@@ -10,10 +10,18 @@
 __version__ = "0.2.2dev"
 
 
-# Try to import only once
+# Try to import only once - also will perform check that at least one GPU is found.
 try:
+    _pyopencl_installed = False
     import pyopencl
-    _pyopencl_installed = True
+    from pyebsdindex.opencl import openclparam
+    testcl = openclparam.OpenClParam()
+    try:
+        gpu = testcl.get_gpu()
+        if len(gpu) > 0:
+            _pyopencl_installed = True
+    except:
+        raise ImportError('pyopencl could not find GPU')
 except ImportError:
     _pyopencl_installed = False
 

pyebsdindex/_ebsd_index_single.py

@@ -18,7 +18,14 @@
 # some notice that they have been modified.
 #
 # Author: David Rowenhorst;
-# The US Naval Research Laboratory Date: 21 Aug 2020
+# The US Naval Research Laboratory Date: 22 May 2024
+
+# For further information see:
+# David J. Rowenhorst, Patrick G. Callahan, Håkon W. Ånes. Fast Radon transforms for
+# high-precision EBSD orientation determination using PyEBSDIndex.
+# Journal of Applied Crystallography, 57(1):3–19, 2024.
+# DOI: 10.1107/S1600576723010221
+
 
 """Setup and handling of Radon indexing runs of EBSD patterns on a
 single thread.
@@ -412,6 +419,12 @@ def update_file(self, filename=None, patDim=np.array([120, 120], dtype=np.int32)
         if filename is None:
             self.filein = None
             self.bandDetectPlan.band_detect_setup(patDim=patDim)
+        elif isinstance(filename, ebsd_pattern.EBSDPatternFile):
+            self.fID = filename  
+            self.bandDetectPlan.band_detect_setup(
+                patDim=[self.fID.patternH, self.fID.patternW]
+            )       
+            self.filein = filename.filepath   
         else:
             self.filein = filename
             self.fID = ebsd_pattern.get_pattern_file_obj(self.filein)
@@ -617,6 +630,102 @@ def _detectbands(self, pats, PC, xyloc=None, clparams=None, verbose=0, chunksize
 
     def _indexbandsphase(self, banddata, bandnorm, verbose=0):
 
+#
+        rhomax = self.bandDetectPlan.rhoMax * (1-self.bandDetectPlan.rhoMaskFrac)
+        shpBandDat = banddata.shape
+        npoints = int(banddata.size/(shpBandDat[-1])+0.1)
+        nPhases = len(self.phaseLib)
+        q = np.zeros((nPhases, npoints, 4))
+        indxData = np.zeros((nPhases + 1, npoints), dtype=self.dataTemplate)
+        bandmatchindex = np.zeros((nPhases, npoints,shpBandDat[-1],2), dtype=np.int32)-100
+        banddataout = banddata.copy()
+
+        indxData["phase"] = -1
+        indxData["fit"] = 180.0
+        indxData["totvotes"] = 0
+        indxData["nmatch"] = 0
+
+
+        if self.phaseLib[0] is None:
+            return indxData, banddata
+
+        if self.nband_earlyexit is None:
+            earlyexit = -1
+            for ph in self.phaselist:
+                if hasattr(ph, 'nband_earlyexit'):
+                    earlyexit = max(earlyexit, ph.nband_earlyexit)
+            if earlyexit < 0:
+                earlyexit = shpBandDat[1]  # default to all the poles.
+            self.nband_earlyexit = earlyexit
+        else:
+            earlyexit = self.nband_earlyexit
+
+
+        adj_intensity = (-1 * np.abs(banddata["rho"]) * 0.5 / rhomax + 1) * banddata["max"]
+        adj_intensity *= ((banddata["theta"] > (2 * np.pi / 180)).astype(np.float32) + 0.5) / 2
+        adj_intensity *= ((banddata["theta"] < (178.0 * np.pi / 180)).astype(np.float32) + 0.5) / 2
+
+        adj_intensity *= banddata["valid"]
+
+        #print(adj_intensity.shape)
+
+
+        for j in range(len(self.phaseLib)):
+
+            indxData['pq'][j, :] = np.sum(banddata['max'] * banddata['valid'], axis=1)
+            p2do = np.ravel(np.nonzero(np.max(indxData["nmatch"], axis=0) < earlyexit)[0])
+
+            if p2do.size ==0:
+                break
+            (
+                avequat,
+                fit,
+                cm,
+                bandmatch,
+                nMatch,
+                matchAttempts,
+                totvotes,
+            ) = self.phaseLib[j].bandindex(
+            bandnorm[p2do, ...], band_intensity=adj_intensity[p2do, ...], band_widths=banddata["width"][p2do, ...], verbose=verbose)
+            whgood = np.nonzero(nMatch >= 3 )[0]
+            if whgood.size > 0:
+                whgood2 = p2do[whgood]
+                q[j, whgood2, :] = avequat[whgood, ...]
+                indxData["fit"][j, whgood2] = fit[whgood]
+                indxData["cm"][j, whgood2] = cm[whgood]
+                indxData["phase"][j, whgood2] = j
+                indxData["nmatch"][j, whgood2] = nMatch[whgood]
+                indxData["matchattempts"][j, whgood2] = matchAttempts[whgood, ...]
+                indxData["totvotes"][j, whgood2] = totvotes[whgood]
+                bandmatchindex[j, whgood2, ..., 1] = bandmatch[whgood, ...]
+
+
+
+
+        qref2detect = self._detector2refframe()
+        q = q.reshape(nPhases * npoints, 4)
+        q = rotlib.quat_multiply(q, qref2detect)
+        q = rotlib.quatnorm(q)
+        q = q.reshape(nPhases, npoints, 4)
+        indxData["quat"][0:nPhases, :, :] = q
+        indxData[-1, :] = indxData[0, :]
+        banddataout['band_match_index'][:,:,:] = bandmatchindex[0,:,:,:].squeeze()
+        if nPhases > 1:
+            for j in range(1, nPhases):
+                # indxData[-1, :] = np.where(
+                #    (indxData[j, :]["cm"] * indxData[j, :]["nmatch"])
+                #    > (indxData[j + 1, :]["cm"] * indxData[j + 1, :]["nmatch"]),
+                #    indxData[j, :],
+                #    indxData[j + 1, :],
+                phasetest = ((3.0 - indxData[j, :]["fit"]) * indxData[j, :]["nmatch"]) \
+                            > ((3.0 - indxData[-1, :]["fit"]) * indxData[-1, :]["nmatch"])
+                whbetter = np.nonzero(phasetest)
+                indxData[-1, whbetter] = indxData[j, whbetter]
+                banddataout['band_match_index'][whbetter,:] =  bandmatchindex[j,whbetter,:,:].squeeze()
+        return indxData, banddataout
+
+    def _indexbandsphase_old(self, banddata, bandnorm, verbose=0):
+
 #
         rhomax = self.bandDetectPlan.rhoMax * (1-self.bandDetectPlan.rhoMaskFrac)
         shpBandDat = banddata.shape
@@ -706,7 +815,6 @@ def _indexbandsphase(self, banddata, bandnorm, verbose=0):
                 indxData[-1, whbetter] = indxData[j, whbetter]
                 banddataout['band_match_index'][whbetter,:] =  bandmatchindex[j,whbetter,:,:].squeeze()
         return indxData, banddataout
-
     def _detector2refframe(self):
         ven = str.upper(self.vendor)
         if ven in ["EDAX", "EMSOFT", "KIKUCHIPY"]:

pyebsdindex/band_detect.py

@@ -1,29 +1,41 @@
-"""This software was developed by employees of the US Naval Research Laboratory (NRL), an
-agency of the Federal Government. Pursuant to title 17 section 105 of the United States
-Code, works of NRL employees are not subject to copyright protection, and this software
-is in the public domain. PyEBSDIndex is an experimental system. NRL assumes no
-responsibility whatsoever for its use by other parties, and makes no guarantees,
-expressed or implied, about its quality, reliability, or any other characteristic. We
-would appreciate acknowledgment if the software is used. To the extent that NRL may hold
-copyright in countries other than the United States, you are hereby granted the
-non-exclusive irrevocable and unconditional right to print, publish, prepare derivative
-works and distribute this software, in any medium, or authorize others to do so on your
-behalf, on a royalty-free basis throughout the world. You may improve, modify, and
-create derivative works of the software or any portion of the software, and you may copy
-and distribute such modifications or works. Modified works should carry a notice stating
-that you changed the software and should note the date and nature of any such change.
-Please explicitly acknowledge the US Naval Research Laboratory as the original source.
-This software can be redistributed and/or modified freely provided that any derivative
-works bear some notice that they are derived from it, and any modified versions bear
-some notice that they have been modified.
-
-Author: David Rowenhorst;
-The US Naval Research Laboratory Date: 21 Aug 2020"""
-
-from os import environ
-from pathlib import PurePath
+# This software was developed by employees of the US Naval Research Laboratory (NRL), an
+# agency of the Federal Government. Pursuant to title 17 section 105 of the United States
+# Code, works of NRL employees are not subject to copyright protection, and this software
+# is in the public domain. PyEBSDIndex is an experimental system. NRL assumes no
+# responsibility whatsoever for its use by other parties, and makes no guarantees,
+# expressed or implied, about its quality, reliability, or any other characteristic. We
+# would appreciate acknowledgment if the software is used. To the extent that NRL may hold
+# copyright in countries other than the United States, you are hereby granted the
+# non-exclusive irrevocable and unconditional right to print, publish, prepare derivative
+# works and distribute this software, in any medium, or authorize others to do so on your
+# behalf, on a royalty-free basis throughout the world. You may improve, modify, and
+# create derivative works of the software or any portion of the software, and you may copy
+# and distribute such modifications or works. Modified works should carry a notice stating
+# that you changed the software and should note the date and nature of any such change.
+# Please explicitly acknowledge the US Naval Research Laboratory as the original source.
+# This software can be redistributed and/or modified freely provided that any derivative
+# works bear some notice that they are derived from it, and any modified versions bear
+# some notice that they have been modified.
+#
+#
+# Author: David Rowenhorst;
+# The US Naval Research Laboratory Date: 22 May 2024
+#
+# For further information see:
+# David J. Rowenhorst, Patrick G. Callahan, Håkon W. Ånes. Fast Radon transforms for
+# high-precision EBSD orientation determination using PyEBSDIndex.
+# Journal of Applied Crystallography, 57(1):3–19, 2024.
+# DOI: 10.1107/S1600576723010221
+#
+#
+
+"""Class that reads in EBSD patterns, and finds the bands within those patterns using a Radon Transform."""
+
+#from os import environ
+import os
+from pathlib import PurePath, Path
 import platform
-import tempfile
+# import tempfile
 from timeit import default_timer as timer
 
 import matplotlib.pyplot as plt
@@ -37,10 +49,12 @@
 
 from pyebsdindex import radon_fast
 
-
-tempdir = PurePath("/tmp" if platform.system() == "Darwin" else tempfile.gettempdir())
-tempdir = tempdir.joinpath('numba')
-environ["NUMBA_CACHE_DIR"] = str(tempdir)
+tempdir = PurePath(Path.home())
+#tempdir = PurePath("/tmp" if platform.system() == "Darwin" else tempfile.gettempdir())
+#tempdir = tempdir.joinpath('numbacache')
+tempdir = tempdir.joinpath('.pyebsdindex').joinpath('numbacache')
+Path(tempdir).mkdir(parents=True, exist_ok=True)
+os.environ["NUMBA_CACHE_DIR"] = str(tempdir)+str(os.sep)
 
 RADEG = 180.0/np.pi
 
@@ -80,6 +94,7 @@ def __init__(
     self.EDAXIQ = False
     self.backgroundsub = None
     self.patternmask = None
+    self.useCPU = True
 
     self.dataType = np.dtype([('id', np.int32), ('max', np.float32), \
                     ('maxloc', np.float32, (2)), ('avemax', np.float32), ('aveloc', np.float32, (2)),\