diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml index 3fb48e5..2528b5b 100644 --- a/.github/workflows/publish.yml +++ b/.github/workflows/publish.yml @@ -29,7 +29,7 @@ jobs: # IMPORTANT: this permission is mandatory for trusted publishing: id-token: write steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - name: Set up Python uses: actions/setup-python@v6 diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index f6a5a8c..6d967d2 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -20,7 +20,7 @@ jobs: name: check manifest runs-on: ubuntu-latest steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - uses: actions/setup-python@v6 with: @@ -49,7 +49,7 @@ jobs: DEPENDENCIES: matplotlib==3.8 numba==0.56.4 pyopencl==2023.1.2 LABEL: -oldest steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v6 diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 2611d32..f0b6513 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -5,6 +5,39 @@ Changelog All notable changes to PyEBSDIndex will be documented in this file. The format is based on `Keep a Changelog `_. + +0.3.10.1 (2026-07-03) +================== +BUG/ERROR FIX RELEASE +--------------------- +There was an error when attempting to multi-process GPU calculations in NLPAR/NLSTEM +that caused new patterns to not be written to the NLPAR/NLSTEM files when using HDF5 as output. +With the GPU enabled calculations, non-parallel diskIO ends up being the rate limiting step, and +thus the multiprocess NLPAR is no longer the default method for discrete GPUs, and a single process +GPU is now the default. Overall safer, even if it does take longer to process files. + +The NLPAR calculations are still done in parallel, either with GPU or multithreaded CPU calculations. + +This should require no changes from the examples behaviors outlined in the tutorials. + +A thanks to Dan Pennachio for alerting me to this issue and useful example data. + +Added +----- + +Changed +------- + +Fixed +----- +- Properly send in the vendor keyword to gnomonic correction setup + (currently no correction is applied). +- Removed multi-process NLPAR as default for discrete GPUs as it mostly was fasteer due-to + slightly risky multi-process file writes, which are not compatible with standard HDF5. +- Fixed region-to-region minimum scaling for ``stem_scale=True`` in NLPAR, now using the file + global minimum. + + 0.3.10 (2026-05-11) ================== Added diff --git a/pyebsdindex/__version.py b/pyebsdindex/__version.py index aa4cd15..1385141 100644 --- a/pyebsdindex/__version.py +++ b/pyebsdindex/__version.py @@ -1 +1 @@ -__version__ = "0.3.10" +__version__ = "0.3.10.1" diff --git a/pyebsdindex/_ebsd_index_parallel.py b/pyebsdindex/_ebsd_index_parallel.py index 27f5de9..d26901d 100644 --- a/pyebsdindex/_ebsd_index_parallel.py +++ b/pyebsdindex/_ebsd_index_parallel.py @@ -48,6 +48,7 @@ from pyebsdindex import band_detect as band_detect os.environ["GRPC_VERBOSITY"] = "ERROR" +os.environ["RAY_ACCEL_ENV_VAR_OVERRIDE_ON_ZERO"] = "0" RAYIPADDRESS = '127.0.0.1' #RAYIPADDRESS = '0.0.0.0' @@ -493,7 +494,7 @@ def index_pats_distributed( ) ) gtaskindex.append(gjob) - + time.sleep(1) #gpu_launched += 1 gpuwrker_cycles = -1000 @@ -524,7 +525,8 @@ def index_pats_distributed( gpuwrker_cycles +=1 - donewrker, busy = ray.wait(gputask,num_returns=len(gputask), timeout=0.1) + #donewrker, busy = ray.wait(gputask,num_returns=len(gputask), timeout=0.1) + donewrker, busy = ray.wait(gputask, num_returns=min(1,len(gputask)), timeout=0.1) #print(len(donewrker), nret) #print() diff --git a/pyebsdindex/_ebsd_index_single.py b/pyebsdindex/_ebsd_index_single.py index 1aebc86..d6a7ca1 100644 --- a/pyebsdindex/_ebsd_index_single.py +++ b/pyebsdindex/_ebsd_index_single.py @@ -423,6 +423,7 @@ def __init__( ("nmatch", np.int32), ("matchattempts", np.int32, 4), ("totvotes", np.int32), + ("grainid", np.int64), ] ) @@ -560,7 +561,7 @@ def index_pats( except: pass - self.gnomonic.calccorrection(PCpat) + #self.gnomonic.calccorrection(PCpat) banddata, bandnorm = self._detectbands(pats, PCpat, xyloc=xyloc, clparams=clparams, verbose=verbose, chunksize=chunksize, gpu_id=gpuid) tic = timer() diff --git a/pyebsdindex/ebsd_pattern.py b/pyebsdindex/ebsd_pattern.py index 69430e1..eb8b210 100644 --- a/pyebsdindex/ebsd_pattern.py +++ b/pyebsdindex/ebsd_pattern.py @@ -218,6 +218,8 @@ def __init__(self,path, filetype=None): self.xStep = None # assumming square grid data, with constant step size self.yStep = None self.xyCenter = np.array([0.0, 0.0]) + self.datamax = None + self.datamin = None # This is the location of the center of the scan relative to center of SEM field-of-view self.hexflag = False self.filetype = filetype @@ -423,6 +425,12 @@ def set_scan_rc(self, rc=(0,0)): # helper function for pattern files that don't self.nRows = np.uint64(rc[0]) self.nPatterns = np.uint64(self.nCols * self.nRows) + def write_datamaxmin(self): # placeholder for non-HDF5 files + pass + + def read_datamaxmin(self): # placeholder for non-HDF5 files + pass + class UPFile(EBSDPatternFile): @@ -1547,6 +1555,37 @@ def pat_writer(self, pat2write, patStart, nPatToWrite, typewrite): patterndset[patStart:patStart+nPatToWrite, :, :] = pat2write[0:nPatToWrite,:,:] f.close() + def write_datamaxmin(self): + try: + f = h5py.File(Path(self.filepath).expanduser(), 'r+') + except: + print("File Not Found:", str(Path(self.filepath))) + return -1 + + patterndset = f[self.h5patdatpth] + + if self.datamax is not None: + patterndset.attrs["DataMax"] = np.float32(self.datamax) + if self.datamin is not None: + patterndset.attrs["DataMin"] = np.float32(self.datamin) + f.close() + + def read_datamaxmin(self): + try: + f = h5py.File(Path(self.filepath).expanduser(), 'r+') + + except: + print("File Not Found:", str(Path(self.filepath))) + return -1 + + patterndset = f[self.h5patdatpth] + + if 'DataMax' in patterndset.attrs: + self.datamax = np.float32(patterndset.attrs["DataMax"]) + if 'DataMin' in patterndset.attrs: + self.datamin = np.float32(patterndset.attrs["DataMin"]) + f.close() + def read_header(self, path=None): if path is not None: self.filepath = path @@ -1931,10 +1970,9 @@ class DM5(HDF5PatFile): def __init__(self, path=None): HDF5PatFile.__init__(self, path) self.vendor = 'GATAN' - # OXFORDOINA only attributes - self.filedatatype = None # np.uint8 - self.patternh5id = 'Data' # Could also be 'Raw Patterns' + self.filedatatype = None # np.uint8 + self.patternh5id = 'Data' if self.filepath is not None: self.get_data_paths() @@ -2012,7 +2050,8 @@ def read_header(self, path=None): #self.yStep = 1.0 #np.float32(headerpath['Y Step'][()][0]) self.xStep = ((f[self.h5patdatpth].parent)["Calibrations/Dimension/[2]"]).attrs['Scale'] #np.float32(headerpath['X Step'][()][0]) self.yStep = ((f[self.h5patdatpth].parent)["Calibrations/Dimension/[3]"]).attrs['Scale'] - + f.close() + self.read_datamaxmin() return 0 # note this function uses multiple returns def read_data(self, path=None, convertToFloat=False, patStartCount=[0, -1], returnArrayOnly=False): @@ -2215,4 +2254,5 @@ def write_data(self, newpatterns=None, patStartCount = [0,-1], writeHead=False, patterndset[int(rowstart):int(rowstart + nrowwrite), int(colstart):int(colstart + ncolwrite), :, :] = pat2write - f.close() \ No newline at end of file + f.close() + diff --git a/pyebsdindex/ebsdfile.py b/pyebsdindex/ebsdfile.py index 05f9c24..fe5bb13 100644 --- a/pyebsdindex/ebsdfile.py +++ b/pyebsdindex/ebsdfile.py @@ -279,7 +279,15 @@ def writeoh5(filename, indexer, data, f.create_dataset(datasetname + '/EBSD/Data/X Position', data=x) f.create_dataset(datasetname + '/EBSD/Data/Y Position', data=y) f.create_dataset(datasetname + '/EBSD/Data/Valid', data=np.zeros(npoints, dtype=np.int8)) - f.create_dataset(datasetname + '/EBSD/Data/SEM Signal', data=np.zeros(npoints, dtype=np.int32)) + semsig = (data[-1]['pq']).copy() + semsig -= semsig.min() + semsig *= float(int(2**8 - 1))/semsig.max() + + semsig = (np.round(semsig)-int(2**7-1)).astype(np.int8) + + f.create_dataset(datasetname + '/EBSD/Data/SEM Signal', data=semsig) + f.create_dataset(datasetname + '/EBSD/Data/Grain ID', data=data[-1]['grainid']) + if version == '8.6': versiontxt = 'OIM Analysis 8.6.103 x64 [29 Sep 2022]' diff --git a/pyebsdindex/grainlabel.py b/pyebsdindex/grainlabel.py new file mode 100644 index 0000000..d510c82 --- /dev/null +++ b/pyebsdindex/grainlabel.py @@ -0,0 +1,567 @@ +""" +================================================================================ +EBSD 3D GRAIN LABELING MODULE +================================================================================ + +This module implements high-performance grain segmentation for 2D or 3D Electron +Backscatter Diffraction (EBSD) data using Numba JIT compilation. + +ALGORITHM OVERVIEW: + - Performs connected-component labeling on 3D voxel grids + - Groups neighboring voxels into "grains" based on: + * Same crystallographic phase + * Misorientation angle below a user-defined tolerance (typically 2–15°) + - Uses quaternion representation for crystal orientations + - Accounts for crystal symmetry via symmetry quaternion operations + - Supports both 2/6-NN (face-sharing) and 9/27-NN (full cube) connectivity + +================================================================================ +""" + +import numpy as np +import numba as nb + +from pyebsdindex.rotlib import quat_multiply + + +# ============================================================================ +# PUBLIC API: High-level entry points +# ============================================================================ + +def grainlabel(pyindxdata, indxer, mistol=3.0, + mask=None, labels=None, nnonly=False, dimensions=None): + """ + High-level grain labeling interface for PyEBSDindex results. + + Extracts quaternions, phase IDs, and symmetry operators from an indexer + object. + + Parameters + ---------- + pyindxdata : result array from pyebsdindex.index_pats() indexing results; + Labels will be made for each phase and the most-likely phase + identifier, pyindxdata[-1]. Each entry in the list requires + the fields: + - 'quat' : ndarray of shape (npts, 4) — quaternion orientations + - 'phase' : ndarray of shape (npts,) — phase ID per point + indxer : object + Indexer object with attributes: + - phaseLib : list of phase objects, each with .qsymops (symmetry quaternions) + - fID : file/scan metadata with .nRows, .nCols for the scan. + mistol : float, optional + Misorientation tolerance in degrees (default: 2.0) + mask : ndarray, optional + Boolean or binary mask; voxels where mask == 0 are pre-labeled as -2, and will + not be labeled. + labels : ndarray, optional + Pre-existing label array; if None, initialized to -1 (unlabeled) + nnonly : bool, optional + If True, use 4/6-NN (face-sharing); if False, use 9/27-NN (default: False) + dimensions : list or array, optional + [nRows, nCols] scan dimensions; extracted from indxer if not provided + + Returns + ------- + labels : a list of length nphases+1 of ndarray of shape (nRows, nCols) or (nSlices, nRows, nCols) for each phase + and for the pyindxdata[-1] (use the most likely phase identifier). + Integer grain labels; -1 = unlabeled, ≥0 = grain ID + + Notes + ----- + - Automatically iterates over all phases in phaseLib + - Reshapes 1D/2D data to 3D for uniform processing + """ + + qsym = [] + for p in indxer.phaseLib: + qsym.append(p.qsymops) + + nphases = len(qsym) + + if dimensions is None: + dimensions = [indxer.fID.nRows, indxer.fID.nCols] + else: + dimensions = dimensions + + + labels = [] + for i in range(nphases): + quats = pyindxdata[i]['quat'] + phaseID = (pyindxdata[i]['phase']).clip(None, 0) + l = phaselabel(quats, phaseID, mistol, [qsym[i]], dimensions, + labels=None, nnonly=nnonly, mask=mask) + labels.append(l) + + quats = pyindxdata[-1]['quat'] + phaseID = (pyindxdata[-1]['phase']) + l = phaselabel(quats, phaseID, mistol, qsym, dimensions, + labels=None, nnonly=nnonly, mask=mask) + labels.append(l) + return labels + +def appendgrainlabels(pyindxdata, labels): + nphase = len(pyindxdata)-1 + + for i in range(nphase): + pyindxdata[i]['grainid'][:] = labels[i].flatten() + + pyindxdata[-1]['grainid'][:] = labels[-1].flatten() + + return pyindxdata + + +def phaselabel(quats, phaseID, mistol, quatsym, dimensions, + labels=None, nnonly=False, mask=None): + """ + Core grain labeling function for 2D or 3D EBSD data. + + Reshapes input data to 3D, initializes label array, and iterates over + each crystallographic phase, calling label3D() for each. + + Parameters + ---------- + quats : ndarray of shape (npts, 4) + Quaternion orientations (w, x, y, z convention) + phaseID : ndarray of shape (npts,) + Phase ID per voxel (integer, typically 0–N_phases) + mistol : float + Misorientation tolerance in degrees + quatsym : list of ndarray or single ndarray + Symmetry quaternions per phase: + - If list: quatsym[i] has shape (nsym_i, 4) for phase i. + - If ndarray: single set of symmetries applied to all phases + dimensions : list of int + [nRows, nCols] for 2D data, or [nSlices, nRows, nCols] for 3D + labels : ndarray, optional + Pre-existing label array; if None, initialized to -1 + nnonly : bool, optional + If True, use 4/6-NN; if False, use 9/27-NN (default: False) + mask : ndarray, optional + Binary mask; voxels where mask == 0 are labeled -1 (masked out) + indicating no label should be applied to these points. Useful for + alternate non-indexed point thresholds. + + Returns + ------- + labels : ndarray + Grain labels, reshaped to match input dimensions + - -1 : unlabeled (no grain assigned) + - ≥0 : grain ID + + Notes + ----- + - Automatically pads 2D data to 3D (nslice=1) for uniform processing + - Iterates over phases in order; each phase gets its own label sequence + - Modifies labels in-place for efficiency + """ + ndim = len(dimensions) + + # ---- Pad dimensions to 3D ---- + dim = np.ones(3, dtype=np.int64) + dim[-ndim:] = np.array(dimensions, dtype=np.int64) + + # ---- Reshape quaternions to 3D ---- + shpq = quats.shape + q = quats.reshape(dim[0], dim[1], dim[2], 4) + phaseID = phaseID.reshape(dim[0], dim[1], dim[2]) + unqphase = np.unique(phaseID[phaseID >= 0]) + nphase = len(unqphase) + + + + # ---- Initialize or reshape label array ---- + if labels is None: + labels = -1 * np.ones((dim[0], dim[1], dim[2]), dtype=np.int64) + else: + labels = labels.reshape(dim[0], dim[1], dim[2]) + + # ---- Apply mask (if provided) ---- + if mask is not None: + mask = mask.reshape(dim[0], dim[1], dim[2]) + labels[mask == 0] = -2 + + # --- mark the unindexed points --- # + # --- this is likely not necessary as PhaseID <= -1 are not labeled ---# + labels[phaseID < 0] = -2 + + # ---- Normalize quatsym to list ---- + if isinstance(quatsym, np.ndarray): + qsym = [quatsym] + else: + qsym = quatsym + + + if len(qsym) != nphase: + raise ValueError('qsym must be a list of length equal to the number of phases') + + # ---- Iterate over phases ---- + phase = 0 + for qsym_i in qsym: + mxl = labels.max() + 1 + labels = label3D(q, phaseID, labels, mistol, qsym_i, phase, mxl, nnonly=nnonly) + phase += 1 + + # ---- Reshape back to original dimensions ---- + if ndim == 2: + labels = labels.reshape(dim[1], dim[2]) + labels = labels.clip(-1) + return labels + + +# ============================================================================ +# HELPER: Quaternion Misorientation with Symmetry +# ============================================================================ + +@nb.njit(fastmath=True, inline='always', boundscheck=False, error_model='numpy') +def misquatdot(q1_0, q1_1, q1_2, q1_3, + q2_0, q2_1, q2_2, q2_3, + qsym, tol): + """ + Tests if the misorientation between two quaternions is within tolerance, + accounting for crystal symmetry. + + Computes the minimum misorientation angle between q1 and q2 by testing + symmetry-equivalent orientations of q2. Returns 1 if any symmetry + variant passes the tolerance check; otherwise 0. This will exit early if + any symmetry variant passes the tolerance check. + + MATHEMATICAL BACKGROUND: + - Quaternion dot product: dot(q1, q2) = w1*w2 + x1*x2 + y1*y2 + z1*z2 + - Misorientation angle θ relates to dot product via: + cos(θ/2) = |dot(q1, q2)| (absolute value accounts for q ≡ -q) + - Tolerance is pre-converted to cos(θ_tol/2) by the caller + - Crystal symmetry: test q1 · (q2 * sym) for symmetry ops + + Parameters + ---------- + q1_0, q1_1, q1_2, q1_3 : float + Components of first quaternion (w, x, y, z) + q2_0, q2_1, q2_2, q2_3 : float + Components of second quaternion (w, x, y, z) + qsym : ndarray of shape (nsym, 4) + Symmetry quaternions; qsym[0] SHOULD be identity [1, 0, 0, 0] + tol : float + Tolerance as cos(θ_tol/2), where θ_tol is the misorientation angle limit + + Returns + ------- + int + 1 if misorientation is within tolerance (accounting for symmetry) + 0 otherwise + + Optimization Notes + ------------------ + - Takes scalar components (avoids array view overhead in Numba) + - Fast path: checks identity symmetry first (most grains match without rotation) + - Inline compilation: function is inlined at call sites + - fastmath=True: allows unsafe floating-point optimizations + + """ + + # ---- Fast path: identity symmetry (qsym[0] = [1, 0, 0, 0]) ---- + # Most neighboring voxels inside a grain match without needing symmetry rotation. + # This avoids the full loop for the common case. + dot0 = abs(q1_0 * q2_0 + q1_1 * q2_1 + q1_2 * q2_2 + q1_3 * q2_3) + #print(dot0) + if dot0 >= tol: + return 1 + + # ---- Full symmetry loop ---- + nsym = qsym.shape[0] + # simple iteration in case identity is NOT first entry -- for PyEBSDindex it always is. + # This should exit early, so should not cause a significant delay. + startstop = [[1, nsym], [0,1]] + for itrange in startstop: + for i in range(itrange[0], itrange[1]): + s0 = qsym[i, 0] + s1 = qsym[i, 1] + s2 = qsym[i, 2] + s3 = qsym[i, 3] + + # --- w,x,y,z = quat_multiply(qsym_i, q2) --- + w = s0 * q2_0 - (s1 * q2_1 + s2 * q2_2 + s3 * q2_3) + x = s0 * q2_1 + q2_0 * s1 + (s2 * q2_3 - s3 * q2_2) + y = s0 * q2_2 + q2_0 * s2 + (s3 * q2_1 - s1 * q2_3) + z = s0 * q2_3 + q2_0 * s3 + (s1 * q2_2 - s2 * q2_1) + + # Dot product with q1 + dot = abs(q1_0 * w + q1_1 * x + q1_2 * y + q1_3 * z) + if dot >= tol: + return 1 + + return 0 + + +# ============================================================================ +# MAIN ALGORITHM: 3D EBSD Grain Labeling +# ============================================================================ + +@nb.njit(parallel=False, fastmath=True, boundscheck=False, error_model='numpy') +def label3D(q, phaseID, labels, mistol, qsym, phase, initialID, nnonly=False): + """ + Numba-optimized 3D connected-component grain labeling for EBSD data. + + Implements a depth-first search (DFS) flood-fill algorithm using a + dynamically-growing 1D stack. For each unlabeled voxel of the target phase, + seeds a new grain and grows it by checking neighbors for: + 1. Same phase ID + 2. Misorientation angle below tolerance (accounting for symmetry) + + ALGORITHM PSEUDOCODE: + ┌──────────────────────────────────────────────────────────────────┐ + │ for each voxel (z, y, x) in 3D grid: │ + │ if phaseID[z,y,x] == phase AND labels[z,y,x] == -1: │ + │ labels[z,y,x] ← labeli │ + │ push (z,y,x) onto stack │ + │ while stack not empty: │ + │ pop (cz, cy, cx) │ + │ for each neighbor (nz, ny, nx): │ + │ if labels[nz,ny,nx] == -1 AND phaseID[nz,ny,nx] == phase:│ + │ if misorientation(q[cz,cy,cx], q[nz,ny,nx]) < tol: │ + │ labels[nz,ny,nx] ← labeli │ + │ push (nz, ny, nx) onto stack │ + │ labeli ← labeli + 1 │ + └──────────────────────────────────────────────────────────────────┘ + + Parameters + ---------- + q : ndarray of shape (nslice, nrow, ncol, 4), dtype float32 + Quaternion orientations (w, x, y, z) at each voxel + phaseID : ndarray of shape (nslice, nrow, ncol), dtype int32 + Crystallographic phase ID at each voxel + labels : ndarray of shape (nslice, nrow, ncol), dtype int64 + Output grain labels (modified in-place) + - -1 : unlabeled + - -2 : masked out (pre-set by caller) + - ≥0 : grain ID + mistol : float + Misorientation tolerance in degrees (e.g., 2.0, 5.0, 15.0) + qsym : ndarray of shape (nsym, 4), dtype float32 + Symmetry quaternions for the current phase + CRITICAL: qsym[0] MUST be the identity quaternion [1, 0, 0, 0] + phase : int32 + Phase ID to label in this call + initialID : uint64 + Starting label ID for this phase (typically max(labels) + 1) + nnonly : bool, optional + If True, use 6-NN (face-sharing neighbors only) + If False, use 27-NN (full 3×3×3 cube, default) + + Returns + ------- + labels : ndarray + Updated label array (same object, modified in-place) + + CONNECTIVITY MODES + ------------------ + 4/6-NN (nnonly=True): + Neighbors: (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) — face-sharing only + Use case: Stricter grain boundaries, fewer spurious merges + + 9/27-NN (nnonly=False, default): + Neighbors: all (dz, dy, dx) ∈ {-1, 0, 1}³ except (0, 0, 0) + Use case: Standard EBSD analysis, captures edge/corner contacts + """ + + nslice, nrow, ncol = labels.shape + nyx = nrow * ncol # Precomputed stride for flattening (z-major order) + + # ---- Convert tolerance from degrees to cosine space ---- + # θ_tol (degrees) → cos(θ_tol/2) for quaternion dot product comparison + tol = np.cos(mistol * np.pi / 180.0 / 2.0) + labeli = initialID + + # ---- Dynamic 1D stack (flattened indices) ---- + # Initialized with capacity = 2*ncol; grows as needed + INITIAL_CAPACITY = 2*ncol + pts2check = np.empty(INITIAL_CAPACITY, dtype=np.int64) + + # ---- Main loop: iterate over all voxels ---- + for z in range(nslice): + for y in range(nrow): + for x in range(ncol): + + # Skip if not target phase + if phaseID[z, y, x] != phase: + continue + # Skip if already labeled (or masked out) + if labels[z, y, x] != -1: + continue + + # ---- Seed new grain ---- + labels[z, y, x] = labeli + pts2check[0] = z * nyx + y * ncol + x + stack_size = 1 + + # ---- Flood fill (DFS) ---- + while stack_size > 0: + stack_size -= 1 + flat = pts2check[stack_size] + + # Decompose flat index back to (z, y, x) + cz = flat // nyx + rem = flat - cz * nyx + cy = rem // ncol + cx = rem - cy * ncol + + # Hoist current quaternion (loaded ONCE per pop) + # Avoids repeated array indexing in the neighbor loop + qc0 = q[cz, cy, cx, 0] + qc1 = q[cz, cy, cx, 1] + qc2 = q[cz, cy, cx, 2] + qc3 = q[cz, cy, cx, 3] + + if not nnonly: + # ============================= + # 27-Nearest Neighbors (3×3×3 cube) + # ============================= + for dz in range(-1, 2): + nz = cz + dz + if nz < 0 or nz >= nslice: + continue + for dy in range(-1, 2): + ny = cy + dy + if ny < 0 or ny >= nrow: + continue + for dx in range(-1, 2): + nx = cx + dx + if nx < 0 or nx >= ncol: + continue + # Skip if already labeled + if labels[nz, ny, nx] != -1: + continue + # Skip if different phase + if phaseID[nz, ny, nx] != phase: + continue + + # Test misorientation + if misquatdot( + q[nz, ny, nx, 0], q[nz, ny, nx, 1], + q[nz, ny, nx, 2], q[nz, ny, nx, 3], + qc0, qc1, qc2, qc3, + qsym, tol + ) > 0: + # Assign label and push to stack + labels[nz, ny, nx] = labeli + if stack_size >= pts2check.shape[0]: + pts2check = _grow_stack_1d(pts2check) + pts2check[stack_size] = nz * nyx + ny * ncol + nx + stack_size += 1 + else: + # ============================= + # 6-Nearest Neighbors (face-sharing only) + # ============================= + # Unrolled for performance: avoids loop overhead + + # -Z neighbor + nz = cz - 1 + if nz >= 0 and labels[nz, cy, cx] == -1 and phaseID[nz, cy, cx] == phase: + if misquatdot( + q[nz, cy, cx, 0], q[nz, cy, cx, 1], + q[nz, cy, cx, 2], q[nz, cy, cx, 3], + qc0, qc1, qc2, qc3, qsym, tol) > 0: + labels[nz, cy, cx] = labeli + if stack_size >= pts2check.shape[0]: + pts2check = _grow_stack_1d(pts2check) + pts2check[stack_size] = nz * nyx + cy * ncol + cx + stack_size += 1 + + # +Z neighbor + nz = cz + 1 + if nz < nslice and labels[nz, cy, cx] == -1 and phaseID[nz, cy, cx] == phase: + if misquatdot( + q[nz, cy, cx, 0], q[nz, cy, cx, 1], + q[nz, cy, cx, 2], q[nz, cy, cx, 3], + qc0, qc1, qc2, qc3, qsym, tol) > 0: + labels[nz, cy, cx] = labeli + if stack_size >= pts2check.shape[0]: + pts2check = _grow_stack_1d(pts2check) + pts2check[stack_size] = nz * nyx + cy * ncol + cx + stack_size += 1 + + # -Y neighbor + ny = cy - 1 + if ny >= 0 and labels[cz, ny, cx] == -1 and phaseID[cz, ny, cx] == phase: + if misquatdot( + q[cz, ny, cx, 0], q[cz, ny, cx, 1], + q[cz, ny, cx, 2], q[cz, ny, cx, 3], + qc0, qc1, qc2, qc3, qsym, tol) > 0: + labels[cz, ny, cx] = labeli + if stack_size >= pts2check.shape[0]: + pts2check = _grow_stack_1d(pts2check) + pts2check[stack_size] = cz * nyx + ny * ncol + cx + stack_size += 1 + + # +Y neighbor + ny = cy + 1 + if ny < nrow and labels[cz, ny, cx] == -1 and phaseID[cz, ny, cx] == phase: + if misquatdot( + q[cz, ny, cx, 0], q[cz, ny, cx, 1], + q[cz, ny, cx, 2], q[cz, ny, cx, 3], + qc0, qc1, qc2, qc3, qsym, tol) > 0: + labels[cz, ny, cx] = labeli + if stack_size >= pts2check.shape[0]: + pts2check = _grow_stack_1d(pts2check) + pts2check[stack_size] = cz * nyx + ny * ncol + cx + stack_size += 1 + + # -X neighbor + nx = cx - 1 + if nx >= 0 and labels[cz, cy, nx] == -1 and phaseID[cz, cy, nx] == phase: + if misquatdot( + q[cz, cy, nx, 0], q[cz, cy, nx, 1], + q[cz, cy, nx, 2], q[cz, cy, nx, 3], + qc0, qc1, qc2, qc3, qsym, tol) > 0: + labels[cz, cy, nx] = labeli + if stack_size >= pts2check.shape[0]: + pts2check = _grow_stack_1d(pts2check) + pts2check[stack_size] = cz * nyx + cy * ncol + nx + stack_size += 1 + + # +X neighbor + nx = cx + 1 + if nx < ncol and labels[cz, cy, nx] == -1 and phaseID[cz, cy, nx] == phase: + if misquatdot( + q[cz, cy, nx, 0], q[cz, cy, nx, 1], + q[cz, cy, nx, 2], q[cz, cy, nx, 3], + qc0, qc1, qc2, qc3, qsym, tol) > 0: + labels[cz, cy, nx] = labeli + if stack_size >= pts2check.shape[0]: + pts2check = _grow_stack_1d(pts2check) + pts2check[stack_size] = cz * nyx + cy * ncol + nx + stack_size += 1 + + # ---- Increment label for next grain ---- + labeli += 1 + + return labels + + +# ============================================================================ +# UTILITY: Dynamic Stack Management +# ============================================================================ + +@nb.njit(inline='always', boundscheck=False) +def _grow_stack_1d(stack): + """ + Double the capacity of a 1D int64 stack buffer and copy contents. + + Called when the stack is full and a new element needs to be pushed. + Allocates a new array with 2× capacity and copies all existing elements. + + Parameters + ---------- + stack : ndarray of shape (n,), dtype int64 + Current stack buffer + + Returns + ------- + new_stack : ndarray of shape (2*n,), dtype int64 + Enlarged stack with contents copied + + """ + old_cap = stack.shape[0] + new_stack = np.empty(old_cap * 2, dtype=np.int64) + for i in range(old_cap): + new_stack[i] = stack[i] + return new_stack \ No newline at end of file diff --git a/pyebsdindex/nlpar.py b/pyebsdindex/nlpar.py index 34c8e09..66ee300 100644 --- a/pyebsdindex/nlpar.py +++ b/pyebsdindex/nlpar.py @@ -60,7 +60,8 @@ if gpusharedmem: from pyebsdindex.opencl.nlpar_cl import NLPAR else: - from pyebsdindex.opencl.nlpar_clray import NLPAR + #from pyebsdindex.opencl.nlpar_clray import NLPAR + from pyebsdindex.opencl.nlpar_cl import NLPAR elif gpuisthere and not _ray_installed: from pyebsdindex.opencl.nlpar_cl import NLPAR diff --git a/pyebsdindex/nlpar_cpu.py b/pyebsdindex/nlpar_cpu.py index 5b57d1d..1717f78 100644 --- a/pyebsdindex/nlpar_cpu.py +++ b/pyebsdindex/nlpar_cpu.py @@ -239,6 +239,9 @@ def calcsigma_cpu(self,chunksize=0,nn=1,saturation_protect=None,automask=None, s ndone = 0 nchunks = int(chunks[1] * chunks[0]) + if stem_scale is True: # need to get the file min (and might as well get the max) + dmin, dmax = self._getdatamaxmin(chunks, patternfile) + for rowchunk in range(chunks[1]): rstart = chunks[3][rowchunk, 0] rend = chunks[3][rowchunk, 1] @@ -254,8 +257,8 @@ def calcsigma_cpu(self,chunksize=0,nn=1,saturation_protect=None,automask=None, s if stem_scale is True: #data = data - data.min() + 1 #data = np.log(data) - data = data - data.min() - data = np.sqrt(data) + data = data.astype(np.float32) - patternfile.datamin + data = np.sqrt(data).astype(np.float32) shp = data.shape data = data.reshape(data.shape[0], phw) @@ -264,7 +267,7 @@ def calcsigma_cpu(self,chunksize=0,nn=1,saturation_protect=None,automask=None, s np.array([0,nrowchunk ], dtype=np.uint64), np.array([0,ncolchunk],dtype=np.uint64), indices,saturation_protect) - + sigma[rstart:rend, cstart:cend] = np.minimum(sigma[rstart:rend, cstart:cend], sigchunk) # temp = (d2 > thresh).choose(dthresh, d2) n2[rstart:rend, cstart:cend,:] = np.select( [n2chunk > 0], [n2chunk], default=n2[rstart:rend, cstart:cend,:]) @@ -273,11 +276,12 @@ def calcsigma_cpu(self,chunksize=0,nn=1,saturation_protect=None,automask=None, s ndone += 1 if verbose >= 2: print("tiles complete: ", ndone, "/", nchunks, sep='', end='\r') - + print() + print(sigma.min(), sigma.max()) return sigma, d2, n2 def calcnlpar_cpu(self, chunksize=0, searchradius=None, lam = None, dthresh = None, - saturation_protect=None, automask=None, stem_scale = None, # see NLPAR __init__ for default values + saturation_protect= None, automask=None, stem_scale = None, # see NLPAR __init__ for default values filename=None, fileout=None, reset_sigma=False, backsub = False, rescale = False,verbose=2, diff_offset=None, **kwargs): @@ -392,7 +396,8 @@ def calcnlpar_cpu(self, chunksize=0, searchradius=None, lam = None, dthresh = No else: # not int, so no rescale. rescale = False - + if stem_scale is True: # need to get the file min (and might as well get the max) + dmin, dmax = self._getdatamaxmin(chunks, patternfile) nthreadpos = numba.get_num_threads() #numba.set_num_threads(18) @@ -442,7 +447,8 @@ def calcnlpar_cpu(self, chunksize=0, searchradius=None, lam = None, dthresh = No jobid += 1 jqueue.append(job) - + newdatamax = -np.inf + newdatamin = np.inf while len(jqueue) > 0: j = jqueue.pop(0) @@ -464,11 +470,11 @@ def calcnlpar_cpu(self, chunksize=0, searchradius=None, lam = None, dthresh = No data, xyloc = patternfile.read_data(patStartCount=[[ cstart, rstart], [ncolchunk, nrowchunk]], convertToFloat=True, returnArrayOnly=True) if stem_scale is True: - datamin = data.min() + datamin = patternfile.datamin # data = data - datamin + 1 # data = np.log(data) - data = data - datamin - data = np.sqrt(data) + data = data.astype(np.float32) - datamin + data = np.sqrt(data).astype(np.float32) shpdata = data.shape @@ -479,7 +485,8 @@ def calcnlpar_cpu(self, chunksize=0, searchradius=None, lam = None, dthresh = No if calcsigma is True: sigchunk = self.sigma_numba(data, 1, nrowchunk, ncolchunk, - [0,nrowchunk], [0,ncolchunk], + np.array([0,nrowchunk], dtype=np.int64), + np.array([0,ncolchunk],dtype=np.int64), indices, saturation_protect)[0] sigchunk = np.minimum(sigma[rstart:rend,cstart:cend], sigchunk) @@ -500,6 +507,7 @@ def calcnlpar_cpu(self, chunksize=0, searchradius=None, lam = None, dthresh = No dataout = dataout.reshape(nrowchunk, ncolchunk, -1) dataout = dataout[rstartcalc: rstartcalc + nrowcalc, cstartcalc:cstartcalc + ncolcalc, :] + if stem_scale is True: #dataout = np.exp(dataout) - 1 + datamin dataout = dataout**2 + datamin @@ -512,6 +520,8 @@ def calcnlpar_cpu(self, chunksize=0, searchradius=None, lam = None, dthresh = No temp *= np.float32(mxval) / temp.max() dataout[i, :, :] = temp + newdatamax = max(newdatamax, np.max(dataout)) + newdatamin = min(newdatamin, np.min(dataout)) patternfileout.write_data(newpatterns=dataout, patStartCount=[[np.int64(cstart + cstartcalc), np.int64(rstart + rstartcalc)], [ncolcalc, nrowcalc]], @@ -520,67 +530,12 @@ def calcnlpar_cpu(self, chunksize=0, searchradius=None, lam = None, dthresh = No if verbose >= 2: print("tiles complete: ", ndone, "/", nchunks, sep='', end='\r') - - # for j in range(0,nrows,chunksize): - # #print('Row start', j) - # if verbose >= 2: - # print("begin row: ", j, "/", nrows, sep='', end='\r') - # - # rowstartread = np.int64(max(0, j-sr)) - # rowend = min(j + chunksize+sr,nrows) - # - # if (rowend - rowstartread) < (2*sr+1): - # rowstartread = np.int64(max(0, rowend - (2*sr+1))) - # rowcountread = np.int64(rowend-rowstartread) - # data, xyloc = patternfile.read_data(patStartCount = [[0,rowstartread], [ncols,rowcountread]], - # convertToFloat=True,returnArrayOnly=True) - # - # shpdata = data.shape - # - # if backsub is True: - # data = self.backsub(data) - # - # - # data = data.reshape(shpdata[0], phw) - # - # rowstartcount = np.asarray([0,rowcountread],dtype=np.int64) - # if calcsigma is True: - # sigchunk, tmp = self.sigma_numba(data,1,rowcountread,ncols,rowstartcount,colstartcount,indices,saturation_protect) - # del tmp - # tmp = (sigma[rowstartread:rowend,:] < sigchunk).choose(sigchunk,sigma[rowstartread:rowend,:]) - # sigma[rowstartread:rowend,:] = tmp - # else: - # sigchunk = sigma[rowstartread:rowend,:] - # - # #dataout = data - # - # dataout = self.nlpar_nb(data,lam, sr, dthresh, sigchunk, - # rowcountread,ncols,indices,saturation_protect, diff_offset=diff_offset) - # - # - # dataout = dataout.reshape(rowcountread, ncols, phw) - # dataout = dataout[j-rowstartread:, :, : ] - # shpout = dataout.shape - # dataout = dataout.reshape(shpout[0]*shpout[1], pheight, pwidth) - # if rescale == True: - # for i in range(dataout.shape[0]): - # temp = dataout[i,:,:] - # temp -= temp.min() - # temp *= np.float32(mxval)/temp.max() - # dataout[i,:,:] = temp - # - # patternfileout.write_data(newpatterns=dataout,patStartCount = [[0,j], [ncols, shpout[0]]], - # flt2int='clip',scalevalue=1.0 ) - # #self.patternfileout.write_data(newpatterns=dataout,patStartCount=[j*ncols,shpout[0]*shpout[1]], - # # flt2int='clip',scalevalue=1.0 ) - # #return dataout - # #sigma[j:j+rowstartcount[1],:] += \ - # # sigchunk[rowstartcount[0]:rowstartcount[0]+rowstartcount[1],:] - - if verbose >= 2: print('', end='') + patternfileout.datamin = newdatamin + patternfileout.datamax = newdatamax + patternfileout.write_datamaxmin() numba.set_num_threads(nthreadpos) return str(patternfileout.filepath) @@ -990,17 +945,46 @@ def opt_lambda(self, **kwargs): # helper function def calcnlpar(self, **kwargs): # helper function return self.calcnlpar_cpu(**kwargs) + def _getdatamaxmin(self, chunks, patternfile): + patternfile.read_datamaxmin() + if (patternfile.datamax is None) or (patternfile.datamin is None): + patternfile.datamax = -np.inf + patternfile.datamin = np.inf + for rowchunk in range(chunks[1]): + rstart = chunks[3][rowchunk, 0] + rend = chunks[3][rowchunk, 1] + nrowchunk = rend - rstart + + for colchunk in range(chunks[0]): + cstart = chunks[2][colchunk, 0] + cend = chunks[2][colchunk, 1] + ncolchunk = cend - cstart + data, xyloc = patternfile.read_data(patStartCount=[[cstart, rstart], [ncolchunk, nrowchunk]], + convertToFloat=True, returnArrayOnly=True) + patternfile.datamax = max(patternfile.datamax, data.max()) + patternfile.datamin = min(patternfile.datamin, data.min()) + patternfile.write_datamaxmin() + dmax = patternfile.datamax + dmin = patternfile.datamin + return dmin, dmax + + def _calcchunks(self, patdim, ncol, nrow, target_bytes=2e9, col_overlap=0, row_overlap=0, col_offset=0, row_offset=0): col_overlap = min(col_overlap, ncol - 1) row_overlap = min(row_overlap, nrow - 1) + + mincolchunk = 2 if ncol >= nrow else 1 + minrowchunk = 2 if nrow > ncol else 1 + + byteperpat = patdim[-1] * patdim[-2] * 4 * 2 # assume a 4 byte float input and output array byteperdataset = byteperpat * ncol * nrow nchunks = int(np.ceil(byteperdataset / target_bytes)) - ncolchunks = (max(np.round(np.sqrt(nchunks * float(ncol) / nrow)), 1)) + ncolchunks = (max(np.round(np.sqrt(nchunks * float(ncol) / nrow)), mincolchunk)) colstep = max((ncol / ncolchunks), 1) ncolchunks = max(ncol / colstep, 1) colstepov = min(colstep + 2 * col_overlap, ncol) @@ -1008,7 +992,7 @@ def _calcchunks(self, patdim, ncol, nrow, target_bytes=2e9, col_overlap=0, row_o colstep = max(int(np.round(colstep)), 1) colstepov = min(colstep + 2 * col_overlap, ncol) - nrowchunks = max(np.ceil(nchunks / ncolchunks), 1) + nrowchunks = max(np.ceil(nchunks / ncolchunks), minrowchunk) rowstep = max((nrow / nrowchunks), 1) nrowchunks = max(nrow / rowstep, 1) rowstepov = min(rowstep + 2 * row_overlap, nrow) diff --git a/pyebsdindex/opencl/nlpar_cl.py b/pyebsdindex/opencl/nlpar_cl.py index a60b68d..a786a51 100644 --- a/pyebsdindex/opencl/nlpar_cl.py +++ b/pyebsdindex/opencl/nlpar_cl.py @@ -247,6 +247,10 @@ def calcsigma_cl(self,nn=1, countchunk = np.zeros((mxchunk, nnn), dtype=np.float32) ndone = 0 nchunks = int(chunks[1] * chunks[0]) + + if stem_scale is True: # need to get the file min (and might as well get the max) + dmin, dmax = self._getdatamaxmin(chunks, patternfile) + for rowchunk in range(chunks[1]): rstart = chunks[3][rowchunk, 0] rend = chunks[3][rowchunk, 1] @@ -262,7 +266,7 @@ def calcsigma_cl(self,nn=1, if stem_scale is True: #data = data - data.min() + 1 #data = np.log(data) - dmin = data.min() + dmin = patternfile.datamin data = data - dmin data = np.sqrt(data).astype(np.float32) @@ -478,6 +482,10 @@ class OpenCLClalcError(Exception): #print("target mem:", target_mem) chunks = self._calcchunks( [pwidth, pheight], ncols, nrows, target_bytes=target_mem, col_overlap=sr, row_overlap=sr) + + if stem_scale is True: # need to get the file min (and might as well get the max) + dmin, dmax = self._getdatamaxmin(chunks, patternfile) + #print(chunks[2], chunks[3]) if verbose >=1: print("lambda:", lam, "search radius:", sr, "dthresh:", dthresh) @@ -549,7 +557,8 @@ class OpenCLClalcError(Exception): jobid += 1 jqueue.append(job) - + newdatamax = -np.inf + newdatamin = np.inf while len(jqueue) > 0: j = jqueue.pop(0) j["nattempts"] += 1 @@ -569,6 +578,12 @@ class OpenCLClalcError(Exception): convertToFloat=False, returnArrayOnly=True) + + + if stem_scale is True: + data = data - patternfile.datamin + data = np.sqrt(data).astype(np.float32) + mxval0 = data.max() mnval0 = data.min() mxval = mxval0 @@ -576,9 +591,6 @@ class OpenCLClalcError(Exception): data -= mnval0 mxval = mxval - mnval0 - if stem_scale is True: - data = np.sqrt(data).astype(np.float32) - if saturation_protect == False: mxval += 1.0 else: @@ -646,12 +658,14 @@ class OpenCLClalcError(Exception): # nothing. It will attempt to reprocess the data 3 times before just writing out # whatever it has. if (mxval0 < np.float32(1.e-8)) or ( mxtest < 0.5 ) or (j["nattempts"] >= 3): - if stem_scale is True: - data = data ** 2 + if mnval0 < 0: data += mnval0 + if stem_scale is True: + data = (data ** 2) + patternfile.datamin + data = data.reshape(nrowcalc * ncolcalc, pheight, pwidth) if rescale == True: for i in range(data.shape[0]): @@ -660,7 +674,8 @@ class OpenCLClalcError(Exception): temp *= np.float32(mxval) / temp.max() data[i, :, :] = temp - + newdatamax = max(newdatamax, np.max(data)) + newdatamin = min(newdatamin, np.min(data)) patternfileout.write_data(newpatterns=data, patStartCount=[[np.int64(cstart + cstartcalc), np.int64(rstart + rstartcalc)], [ncolcalc, nrowcalc]], @@ -692,6 +707,9 @@ class OpenCLClalcError(Exception): print('', end='') queue.finish() queue = None + patternfileout.datamin = newdatamin + patternfileout.datamax = newdatamax + patternfileout.write_datamaxmin() return str(patternfileout.filepath) diff --git a/pyebsdindex/opencl/nlpar_clray.py b/pyebsdindex/opencl/nlpar_clray.py index 7f260a9..eda1d83 100644 --- a/pyebsdindex/opencl/nlpar_clray.py +++ b/pyebsdindex/opencl/nlpar_clray.py @@ -29,6 +29,7 @@ import os, sys, platform import logging +import time from timeit import default_timer as timer import numpy as np import pyopencl as cl @@ -45,6 +46,7 @@ # RAYIPADDRESS = '0.0.0.0' # the localhost address does not work on macOS when on a VPN os.environ["GRPC_VERBOSITY"] = "ERROR" +os.environ["RAY_ACCEL_ENV_VAR_OVERRIDE_ON_ZERO"] = "0" class NLPAR(nlpar_cl.NLPAR): @@ -133,7 +135,8 @@ def calcsigma_clray(self, nn=1, saturation_protect=None, automask=None, # stem_scale=stem_scale, # gpu_id=gpu_id, **kwargs) - target_mem = clparams.gpu[gpu_id].max_mem_alloc_size // 2 + #target_mem = clparams.gpu[gpu_id].max_mem_alloc_size // 2 + target_mem = min(clparams.queue.device.max_mem_alloc_size // 2, np.int64(64e9)) max_mem = clparams.gpu[gpu_id].global_mem_size * 0.5 if target_mem * ngpuwrker > max_mem: #print('revisemem:') @@ -159,6 +162,10 @@ def calcsigma_clray(self, nn=1, saturation_protect=None, automask=None, col_overlap=nn, row_overlap=nn) + if stem_scale is True: # need to get the file min (and might as well get the max) + dmin, dmax = self._getdatamaxmin(chunks, self.patternfile) + + jobqueue = [] for rowchunk in range(chunks[1]): @@ -183,7 +190,6 @@ def calcsigma_clray(self, nn=1, saturation_protect=None, automask=None, [cstart, cend, rstart, rend], \ [cstartcalc, cendcalc, rstartcalc, rendcalc])) - # wrker = NLPARGPUWorker(actorid=1, gpu_id=gpu_id, cudavis=cudavis) # job = jobqueue[0] # @@ -282,12 +288,14 @@ def _sigmachunkcalc_cl(self, data, calclim, clparams=None): if stem_scale == True: + # data = data - data.min() + 1 # data = np.log(data) - dmin = data.min() - data = data - dmin + dmin = self.patternfile.datamin + data = data.astype(np.float32) - dmin data = np.sqrt(data).astype(np.float32) + ctx = clparams.ctx prg = clparams.prg clkern = clparams.kernels @@ -378,7 +386,7 @@ def _sigmachunkcalc_cl(self, data, calclim, clparams=None): sigmachunk_gpu, count_local, dist_local, np.int64(nn)) - clparams.queue.finish() + cl.enqueue_copy(clparams.queue, distchunk, dist_local, is_blocking=False) cl.enqueue_copy(clparams.queue, countchunk, count_local, is_blocking=False) @@ -386,14 +394,12 @@ def _sigmachunkcalc_cl(self, data, calclim, clparams=None): #sigmachunk_gpu.release() - clparams.queue.finish() - dist_local.release() count_local.release() datapad_gpu.release() clparams.queue.flush() clparams.queue = None - #self.sigma = sigma + return sigmachunk, distchunk, countchunk @@ -526,7 +532,7 @@ def calcnlpar_clray(self, searchradius=None, lam = None, dthresh = None, # gpu_id= gpu_id, # diff_offset=diff_offset) - target_mem = clparams.gpu[gpu_id].max_mem_alloc_size//6 + target_mem = min(clparams.gpu[gpu_id].max_mem_alloc_size// 6, np.int64(4e9)) #clparams.gpu[gpu_id].max_mem_alloc_size//6 max_mem = clparams.gpu[gpu_id].global_mem_size*0.4 if target_mem*ngpuwrker > max_mem: target_mem = max_mem/ngpuwrker @@ -535,6 +541,9 @@ def calcnlpar_clray(self, searchradius=None, lam = None, dthresh = None, chunks = self._calcchunks([pwidth, pheight], ncols, nrows, target_bytes=target_mem, col_overlap=sr, row_overlap=sr) + if stem_scale is True: # need to get the file min (and might as well get the max) + dmin, dmax = self._getdatamaxmin(chunks, self.patternfile) + nnn = int((2 * sr + 1) ** 2) jobqueue = [] @@ -587,6 +596,9 @@ def calcnlpar_clray(self, searchradius=None, lam = None, dthresh = None, idlewrker.append(NLPARGPUWorker.options(num_cpus=float(0.99), num_gpus=ngpu_per_wrker).remote( actorid=w, gpu_id=gpu_id, cudavis=cudavis)) + newdatamax = -np.inf + newdatamin = np.inf + njobs = len(jobqueue) ndone = 0 while ndone < njobs: @@ -598,12 +610,20 @@ def calcnlpar_clray(self, searchradius=None, lam = None, dthresh = None, tasks.append(wrker.runnlpar_chunk.remote(job, nlparobj=nlpar_remote)) busywrker.append(wrker) if len(tasks) > 0: - donetasks, stillbusy = ray.wait(tasks, num_returns=len(busywrker), timeout=0.1) - + #donetasks, stillbusy = ray.wait(tasks, num_returns=len(busywrker), timeout=0.1) + donetasks, stillbusy = ray.wait(tasks, num_returns=1, timeout=0.1) for tsk in donetasks: indx = tasks.index(tsk) - message, job, newdata = ray.get(tsk) + message, gpujob, newdata = ray.get(tsk) if message == 'Done': + newdatamax = max(newdatamax, np.max(newdata)) + newdatamin = min(newdatamin, np.min(newdata)) + + self.patternfileout.write_data(newpatterns=newdata, + patStartCount=[[gpujob.cstart + gpujob.cstartcalc, + gpujob.rstart + gpujob.rstartcalc], + [gpujob.ncolcalc, gpujob.nrowcalc]], + flt2int='clip', scalevalue=1.0) idlewrker.append(busywrker.pop(indx)) tasks.remove(tsk) ndone += 1 @@ -615,6 +635,10 @@ def calcnlpar_clray(self, searchradius=None, lam = None, dthresh = None, if verbose >= 2: print('\n', end='') + self.patternfileout.datamin = newdatamin + self.patternfileout.datamax = newdatamax + self.patternfileout.write_datamaxmin() + return str(self.patternfileout.filepath) def _nlparchunkcalc_cl(self, data, calclim, clparams=None): @@ -660,10 +684,16 @@ def _nlparchunkcalc_cl(self, data, calclim, clparams=None): if stem_scale == True: # data = data - data.min() + 1 # data = np.log(data) - mndat = data.min() - data = data - mndat + data = data.astype(np.float32) - self.patternfile.datamin data = np.sqrt(data).astype(np.float32) + mxval0 = data.max() + mnval0 = data.min() + mxval = mxval0 + if mnval0 < 0: + data -= mnval0 + mxval = mxval - mnval0 + #print(chunks[2], chunks[3]) #print(lam, sr, dthresh) @@ -728,18 +758,23 @@ def _nlparchunkcalc_cl(self, data, calclim, clparams=None): np.float32(dthresh), np.float32(diff_offset)) - data = data.astype(np.float32) # prepare to receive data back from GPU - data.reshape(-1)[:] = 0.0 + + data = np.zeros(int(npadmx), dtype=np.float32) #data.astype(np.float32) # prepare to receive data back from GPU + #data.reshape(-1)[:] = 0.0 data = data.reshape(nrowchunk, ncolchunk, pheight, pwidth) + #print(data.min(), data.max()) cl.enqueue_copy(clparams.queue, data, datapadout_gpu, is_blocking=True) sigmachunk_gpu.release() + datapadout_gpu.release() clparams.queue.finish() + if mnval0 < 0: + data += mnval0 + if stem_scale == True: # data = data - data.min() + 1 # data = np.log(data) - data = data**2 - data += mndat + data = data**2 + self.patternfile.datamin if self.rescale == True: @@ -780,7 +815,7 @@ def __init__(self, actorid=0, gpu_id=None, cudavis = '0'): def runsigma_chunk(self,gpujob, nlparobj=None, **kwargs): if gpujob is None: #time.sleep(0.001) - return 'Bored', (None, None, None) + return 'Bored', None, None try: # print(type(self.openCLParams.ctx)) gpujob._starttime() @@ -793,8 +828,6 @@ def runsigma_chunk(self,gpujob, nlparobj=None, **kwargs): convertToFloat=False, returnArrayOnly=True) - - newdata = nlparobj._sigmachunkcalc_cl(data, gpujob, clparams=self.openCLParams, **kwargs) if self.openCLParams.queue is not None: @@ -814,7 +847,7 @@ def runnlpar_chunk(self, gpujob, nlparobj=None): if gpujob is None: #time.sleep(0.001) - return 'Bored', (None, None, None) + return 'Bored', None, None try: # print(type(self.openCLParams.ctx)) gpujob._starttime() @@ -822,28 +855,29 @@ def runnlpar_chunk(self, gpujob, nlparobj=None): #if self.openCLParams is not None: # self.openCLParams.get_queue() + data, xyloc = nlparobj.patternfile.read_data(patStartCount=[[gpujob.cstart, gpujob.rstart], [gpujob.ncolchunk, gpujob.nrowchunk]], convertToFloat=False, returnArrayOnly=True) - + #print(data.min(), data.max()) newdata = nlparobj._nlparchunkcalc_cl(data, gpujob, clparams=self.openCLParams) - + #print(newdata.min(), newdata.max()) if self.openCLParams.queue is not None: print("queue still here") self.openCLParams.queue.finish() self.openCLParams.queue = None - nlparobj.patternfileout.write_data(newpatterns=newdata, - patStartCount=[[gpujob.cstart + gpujob.cstartcalc, - gpujob.rstart + gpujob.rstartcalc], - [gpujob.ncolcalc, gpujob.nrowcalc]], - flt2int='clip', scalevalue=1.0) + # nlparobj.patternfileout.write_data(newpatterns=newdata, + # patStartCount=[[gpujob.cstart + gpujob.cstartcalc, + # gpujob.rstart + gpujob.rstartcalc], + # [gpujob.ncolcalc, gpujob.nrowcalc]], + # flt2int='clip', scalevalue=1.0) gpujob._endtime() - return 'Done', gpujob, None + return 'Done', gpujob, newdata except Exception as e: print(e) gpujob.rate = None