peak_fit.py 15.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
#!/usr/bin/python
# coding: utf8
# /*##########################################################################
#
# Copyright (c) 2015-2016 European Synchrotron Radiation Facility
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
# ###########################################################################*/

Damien Naudet's avatar
WIP  
Damien Naudet committed
27 28
from __future__ import absolute_import

29 30 31 32 33 34 35
__authors__ = ["D. Naudet"]
__date__ = "01/06/2016"
__license__ = "MIT"

import time
import ctypes
import multiprocessing as mp
Damien Naudet's avatar
Damien Naudet committed
36
from threading import Thread
37 38 39 40
import multiprocessing.sharedctypes as mp_sharedctypes

import numpy as np

Damien Naudet's avatar
Damien Naudet committed
41
# from silx.math import curve_fit
Thomas Vincent's avatar
Thomas Vincent committed
42
from ... import config
Damien Naudet's avatar
Damien Naudet committed
43 44 45
from ...io import QSpaceH5
from ...fit import (GaussianFitter, GaussianResults,
                    CentroidFitter, CentroidResults,
46
                    SilxFitter, SilxResults)
Damien Naudet's avatar
Damien Naudet committed
47
from .sharedresults import FitTypes
48 49 50

disp_times = False

Damien Naudet's avatar
WIP  
Damien Naudet committed
51

52
class BackgroundTypes(object):
53 54 55 56 57 58
    """Enum of background subtraction types:

    - NONE: No background
    - CONSTANT: Remove constant (= min of the data) background
    - LINEAR: Remove linear background using line from first to last data
    """
Thomas Vincent's avatar
Thomas Vincent committed
59 60 61 62
    NONE = "None"
    CONSTANT = "Constant"
    LINEAR = "Linear"
    ALLOWED = NONE, CONSTANT, LINEAR
63 64


Damien Naudet's avatar
Damien Naudet committed
65
class PeakFitter(Thread):
Damien Naudet's avatar
Damien Naudet committed
66
    """
Thomas Vincent's avatar
Thomas Vincent committed
67
    :param str qspace_f: path to the HDF5 file containing the qspace cubes
Damien Naudet's avatar
Damien Naudet committed
68

Thomas Vincent's avatar
Thomas Vincent committed
69
    :param FitTypes fit_type:
Damien Naudet's avatar
Damien Naudet committed
70 71 72 73

    :param indices: indices of the cubes (in the input HDF5 dataset) for which
        the qx/qy/qz peaks coordinates will be computed. E.g : if the array
        [1, 2, 3] is provided, only those cubes will be fitted.
Thomas Vincent's avatar
Thomas Vincent committed
74
    :type indices: *optional* `array_like`
Damien Naudet's avatar
Damien Naudet committed
75

Thomas Vincent's avatar
Thomas Vincent committed
76 77
    :param Union[int,None] n_proc:
        Number of process to use. If None, the config value is used.
78 79

    :param BackgroundTypes background: The background subtraction to perform
Damien Naudet's avatar
Damien Naudet committed
80
    """
Damien Naudet's avatar
Damien Naudet committed
81

Damien Naudet's avatar
Damien Naudet committed
82 83 84 85
    READY, RUNNING, DONE, ERROR, CANCELED = __STATUSES = range(5)

    def __init__(self,
                 qspace_f,
Damien Naudet's avatar
Damien Naudet committed
86
                 fit_type=FitTypes.GAUSSIAN,
87
                 n_peaks=1,
Damien Naudet's avatar
Damien Naudet committed
88 89
                 indices=None,
                 n_proc=None,
90
                 roi_indices=None,
Thomas Vincent's avatar
Thomas Vincent committed
91
                 background=BackgroundTypes.NONE):
Damien Naudet's avatar
Damien Naudet committed
92 93 94 95 96
        super(PeakFitter, self).__init__()

        self.__results = None
        self.__thread = None
        self.__progress = 0
Damien Naudet's avatar
Damien Naudet committed
97
        self.__callback = None
Damien Naudet's avatar
Damien Naudet committed
98 99 100 101 102 103 104

        self.__status = self.READY

        self.__indices = None

        self.__qspace_f = qspace_f
        self.__fit_type = fit_type
105
        self.__background = background
106
        self.__n_peaks = n_peaks
Damien Naudet's avatar
Damien Naudet committed
107 108 109 110

        if n_proc:
            self.__n_proc = n_proc
        else:
Thomas Vincent's avatar
Thomas Vincent committed
111
            n_proc = self.__n_proc = config.DEFAULT_PROCESS_NUMBER
Damien Naudet's avatar
Damien Naudet committed
112 113 114 115 116

        self.__shared_progress = mp_sharedctypes.RawArray(ctypes.c_int32,
                                                          n_proc)

        if roi_indices is not None:
Damien Naudet's avatar
Damien Naudet committed
117
            self.__roi_indices = np.array(roi_indices[:])
Damien Naudet's avatar
Damien Naudet committed
118 119 120 121 122
        else:
            self.__roi_indices = None

        if fit_type not in FitTypes.ALLOWED:
            self.__set_status(self.ERROR)
123 124 125 126 127
            raise ValueError('Unknown fit type: {0}'.format(fit_type))

        if background not in BackgroundTypes.ALLOWED:
            self.__set_status(self.ERROR)
            raise ValueError('Unknown background type: {}'.format(background))
128

Damien Naudet's avatar
Damien Naudet committed
129 130 131 132
        try:
            with QSpaceH5.QSpaceH5(qspace_f) as qspace_h5:
                with qspace_h5.qspace_dset_ctx() as dset:
                    qdata_shape = dset.shape
Damien Naudet's avatar
Damien Naudet committed
133

Damien Naudet's avatar
Damien Naudet committed
134
                n_points = qdata_shape[0]
135

Damien Naudet's avatar
Damien Naudet committed
136
                if indices is None:
Damien Naudet's avatar
Damien Naudet committed
137
                    indices = list(range(n_points))
Damien Naudet's avatar
Damien Naudet committed
138 139 140 141 142
                else:
                    indices = indices[:]
        except IOError:
            self.__set_status(self.ERROR)
            raise
Damien Naudet's avatar
Damien Naudet committed
143

Damien Naudet's avatar
Damien Naudet committed
144
        self.__indices = np.array(indices)
Damien Naudet's avatar
Damien Naudet committed
145 146 147 148 149 150 151

    def __set_status(self, status):
        assert status in self.__STATUSES
        self.__status = status

    status = property(lambda self: self.__status)

Damien Naudet's avatar
Damien Naudet committed
152 153
    results = property(lambda self: self.__results)

Damien Naudet's avatar
Damien Naudet committed
154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
    def peak_fit(self,
                 blocking=True,
                 callback=None):

        if self.__thread and self.__thread.is_alive():
            raise RuntimeError('A fit is already in progress.')

        self.__results = None

        if blocking:
            return self.__peak_fit()
        else:
            thread = self.__thread = Thread(target=self.__peak_fit)
            self.__callback = callback
            thread.start()

    def progress(self):
        return (100.0 *
                np.frombuffer(self.__shared_progress, dtype='int32').max() /
                (len(self.__indices) - 1))

    def __peak_fit(self):

        self.__set_status(self.RUNNING)

        qspace_f = self.__qspace_f
        fit_type = self.__fit_type
        indices = self.__indices
        n_proc = self.__n_proc
        roi_indices = self.__roi_indices
        shared_progress = self.__shared_progress
185
        n_peaks = self.__n_peaks
Damien Naudet's avatar
Damien Naudet committed
186 187 188 189 190

        t_total = time.time()

        progress = np.frombuffer(shared_progress, dtype='int32')
        progress[:] = 0
191 192

        n_indices = len(indices)
Damien Naudet's avatar
Damien Naudet committed
193

Damien Naudet's avatar
Damien Naudet committed
194 195
        try:
            with QSpaceH5.QSpaceH5(qspace_f) as qspace_h5:
Damien Naudet's avatar
Damien Naudet committed
196
                with qspace_h5:
Damien Naudet's avatar
Damien Naudet committed
197 198 199 200 201 202 203 204 205 206 207 208
                    x_pos = qspace_h5.sample_x[indices]
                    y_pos = qspace_h5.sample_y[indices]
        except IOError:
            self.__set_status(self.ERROR)
            raise

            # shared_res = mp_sharedctypes.RawArray(ctypes.c_double, n_indices * 9)
            # # TODO : find something better
            # shared_success = mp_sharedctypes.RawArray(ctypes.c_bool, n_indices)
            # # success = np.ndarray((n_indices,), dtype=np.bool)
            # # success[:] = True

Damien Naudet's avatar
Damien Naudet committed
209
        if fit_type == FitTypes.GAUSSIAN:
Damien Naudet's avatar
Damien Naudet committed
210
            fit_class = GaussianFitter
211 212 213 214
            n_peaks = n_peaks if n_peaks >= 1 else 1
            shared_results = GaussianResults(n_points=n_indices,
                                             n_peaks=n_peaks)
        elif fit_type == FitTypes.CENTROID:
Damien Naudet's avatar
Damien Naudet committed
215
            fit_class = CentroidFitter
Damien Naudet's avatar
Damien Naudet committed
216
            shared_results = CentroidResults(n_points=n_indices)
217 218 219 220 221
        elif fit_type == FitTypes.SILX:
            fit_class = SilxFitter
            n_peaks = n_peaks if n_peaks >= 1 else 1
            shared_results = SilxResults(n_points=n_indices,
                                         n_peaks=n_peaks)
Thomas Vincent's avatar
Thomas Vincent committed
222 223
        else:
            raise RuntimeError('Unknown Fit Type')
Damien Naudet's avatar
Damien Naudet committed
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264

        # with h5py.File(qspace_f, 'r') as qspace_h5:
        #
        #     q_x = qspace_h5['bins_edges/x'][:]
        #     q_y = qspace_h5['bins_edges/y'][:]
        #     q_z = qspace_h5['bins_edges/z'][:]
        #     qdata = qspace_h5['data/qspace']
        #
        #     n_points = qdata.shape[0]
        #
        #     if indices is None:
        #         indices = range(n_points)
        #
        #     n_indices = len(indices)
        #
        #     x_pos = qspace_h5['geom/x'][indices]
        #     y_pos = qspace_h5['geom/y'][indices]
        #
        #     shared_res = mp_sharedctypes.RawArray(ctypes.c_double, n_indices * 9)
        #     # TODO : find something better
        #     shared_success = mp_sharedctypes.RawArray(ctypes.c_bool, n_indices)
        #     success = np.ndarray((n_indices,), dtype=np.bool)
        #     success[:] = True
        #
        #     # this has to be done otherwise h5py complains about not being
        #     # able to open compressed datasets from other processes
        #     del qdata

        # results = np.ndarray((n_indices, 11), dtype=np.double)
        # results[:, 0] = x_pos
        # results[:, 1] = y_pos

        manager = mp.Manager()

        read_lock = manager.Lock()
        idx_queue = manager.Queue()

        pool = mp.Pool(n_proc,
                       initializer=_init_thread,
                       initargs=(shared_results,
                                 shared_progress,
Damien Naudet's avatar
Damien Naudet committed
265
                                 fit_class,
Damien Naudet's avatar
Damien Naudet committed
266 267 268
                                 (n_indices, 9),
                                 idx_queue,
                                 qspace_f,
269 270
                                 read_lock,
                                 self.__background))
Damien Naudet's avatar
Damien Naudet committed
271 272

        if disp_times:
Thomas Vincent's avatar
Thomas Vincent committed
273
            class MyTimes(object):
Damien Naudet's avatar
Damien Naudet committed
274 275 276 277 278 279 280 281 282 283 284 285 286
                def __init__(self):
                    self.t_read = 0.
                    self.t_mask = 0.
                    self.t_fit = 0.
                    self.t_write = 0.

                def update(self, arg):
                    (t_read_, t_mask_, t_fit_, t_write_) = arg
                    self.t_read += t_read_
                    self.t_mask += t_mask_
                    self.t_fit += t_fit_
                    self.t_write += t_write_

Thomas Vincent's avatar
Thomas Vincent committed
287
            res_times = MyTimes()
Damien Naudet's avatar
Damien Naudet committed
288 289 290 291 292 293 294 295
            callback = res_times.update
        else:
            callback = None

        # creating the processes
        res_list = []
        for th_idx in range(n_proc):
            arg_list = (th_idx, roi_indices)
Damien Naudet's avatar
Damien Naudet committed
296 297 298
            res = pool.apply_async(_fit_process,
                                   args=arg_list,
                                   callback=callback)
Damien Naudet's avatar
Damien Naudet committed
299 300 301
            res_list.append(res)

        # sending the image indices
302 303
        for i_fit, i_cube in enumerate(indices):
            idx_queue.put([i_fit, i_cube])
Damien Naudet's avatar
Damien Naudet committed
304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336

        # sending the None value to let the threads know that they should return
        for th_idx in range(n_proc):
            idx_queue.put(None)

        pool.close()
        pool.join()

        t_total = time.time() - t_total
        if disp_times:
            print('Total : {0}.'.format(t_total))
            print('Read {0}'.format(res_times.t_read))
            print('Mask {0}'.format(res_times.t_mask))
            print('Fit {0}'.format(res_times.t_fit))
            print('Write {0}'.format(res_times.t_write))

        with QSpaceH5.QSpaceH5(qspace_f) as qspace_h5:
            q_x = qspace_h5.qx
            q_y = qspace_h5.qy
            q_z = qspace_h5.qz

        if roi_indices is not None:
            xSlice = slice(roi_indices[0][0], roi_indices[0][1], 1)
            ySlice = slice(roi_indices[1][0], roi_indices[1][1], 1)
            zSlice = slice(roi_indices[2][0], roi_indices[2][1], 1)
            q_x = q_x[xSlice]
            q_y = q_y[ySlice]
            q_z = q_z[zSlice]

        fit_results = shared_results.fit_results(sample_x=x_pos,
                                                 sample_y=y_pos,
                                                 q_x=q_x,
                                                 q_y=q_y,
Thomas Vincent's avatar
Thomas Vincent committed
337 338
                                                 q_z=q_z,
                                                 background_mode=self.__background)
Damien Naudet's avatar
Damien Naudet committed
339 340 341 342 343

        self.__results = fit_results

        self.__set_status(self.DONE)

Damien Naudet's avatar
Damien Naudet committed
344 345 346
        if self.__callback:
            self.__callback()

Damien Naudet's avatar
Damien Naudet committed
347
        return fit_results
Damien Naudet's avatar
WIP  
Damien Naudet committed
348 349


Damien Naudet's avatar
Damien Naudet committed
350
def _init_thread(shared_res_,
Damien Naudet's avatar
Damien Naudet committed
351
                 shared_prog_,
Damien Naudet's avatar
Damien Naudet committed
352
                 fit_class_,
Damien Naudet's avatar
Damien Naudet committed
353 354 355
                 result_shape_,
                 idx_queue_,
                 qspace_f_,
356 357
                 read_lock_,
                 background_):
Damien Naudet's avatar
Damien Naudet committed
358
    global shared_res, \
Damien Naudet's avatar
Damien Naudet committed
359
        shared_progress, \
Damien Naudet's avatar
Damien Naudet committed
360
        fit_class, \
Damien Naudet's avatar
Damien Naudet committed
361 362 363
        result_shape, \
        idx_queue, \
        qspace_f, \
364 365
        read_lock, \
        fit_background
Damien Naudet's avatar
Damien Naudet committed
366 367

    shared_res = shared_res_
Damien Naudet's avatar
Damien Naudet committed
368
    shared_progress = shared_prog_
Damien Naudet's avatar
Damien Naudet committed
369
    fit_class = fit_class_
Damien Naudet's avatar
Damien Naudet committed
370 371 372 373
    result_shape = result_shape_
    idx_queue = idx_queue_
    qspace_f = qspace_f_
    read_lock = read_lock_
374
    fit_background = background_
Damien Naudet's avatar
Damien Naudet committed
375

Damien Naudet's avatar
WIP  
Damien Naudet committed
376

Damien Naudet's avatar
Damien Naudet committed
377
def _fit_process(th_idx, roiIndices=None):
Damien Naudet's avatar
Damien Naudet committed
378
    print('Thread {0} started.'.format(th_idx))
Damien Naudet's avatar
Damien Naudet committed
379 380 381
    try:
        t_read = 0.
        t_fit = 0.
382
        t_mask = 0.
Damien Naudet's avatar
Damien Naudet committed
383

384
        l_shared_res = shared_res.local_copy()
Damien Naudet's avatar
Damien Naudet committed
385
        progress = np.frombuffer(shared_progress, dtype='int32')
386

Damien Naudet's avatar
Damien Naudet committed
387
        qspace_h5 = QSpaceH5.QSpaceH5(qspace_f)
Damien Naudet's avatar
Damien Naudet committed
388

Damien Naudet's avatar
Damien Naudet committed
389
        # Put this in the main thread
Damien Naudet's avatar
Damien Naudet committed
390 391 392 393 394 395 396
        if roiIndices is not None:
            xSlice = slice(roiIndices[0][0], roiIndices[0][1], 1)
            ySlice = slice(roiIndices[1][0], roiIndices[1][1], 1)
            zSlice = slice(roiIndices[2][0], roiIndices[2][1], 1)

        # TODO : timeout to check if it has been canceled
        # read_lock.acquire()
Damien Naudet's avatar
Damien Naudet committed
397
        with qspace_h5:
Damien Naudet's avatar
Damien Naudet committed
398 399 400 401 402 403 404
            q_x = qspace_h5.qx
            q_y = qspace_h5.qy
            q_z = qspace_h5.qz
            with qspace_h5.qspace_dset_ctx() as dset:
                q_shape = dset.shape
                q_dtype = dset.dtype
            histo = qspace_h5.histo
Damien Naudet's avatar
Damien Naudet committed
405

406
            if roiIndices is not None:
Damien Naudet's avatar
Damien Naudet committed
407 408 409 410 411
                q_x = q_x[xSlice]
                q_y = q_y[ySlice]
                q_z = q_z[zSlice]
                histo = histo[xSlice, ySlice, zSlice]

Damien Naudet's avatar
Damien Naudet committed
412 413 414
            mask = np.where(histo > 0)
            weights = histo[mask]

Damien Naudet's avatar
Damien Naudet committed
415 416 417
        # read_lock.release()
        read_cube = np.ascontiguousarray(np.zeros(q_shape[1:]),
                                         dtype=q_dtype)
418

Damien Naudet's avatar
Damien Naudet committed
419
        fitter = fit_class(q_x, q_y, q_z, l_shared_res)
Damien Naudet's avatar
Damien Naudet committed
420

Damien Naudet's avatar
Damien Naudet committed
421 422
        while True:
            # TODO : timeout
Thomas Vincent's avatar
Thomas Vincent committed
423
            data = idx_queue.get()
424

Thomas Vincent's avatar
Thomas Vincent committed
425
            if data is None:
Damien Naudet's avatar
Damien Naudet committed
426
                break
427

Thomas Vincent's avatar
Thomas Vincent committed
428
            i_fit, i_cube = data
429 430

            progress[th_idx] = i_fit
Damien Naudet's avatar
Damien Naudet committed
431

432
            if i_cube % 100 == 0:
Thomas Vincent's avatar
Thomas Vincent committed
433
                print('Processing cube {0}/{1}.'.format(i_fit, result_shape[0]))
434

Damien Naudet's avatar
Damien Naudet committed
435
            t0 = time.time()
Damien Naudet's avatar
Damien Naudet committed
436
            with qspace_h5.qspace_dset_ctx() as dset:
Damien Naudet's avatar
Damien Naudet committed
437
                dset.read_direct(read_cube,
Damien Naudet's avatar
Damien Naudet committed
438 439
                                 source_sel=np.s_[i_cube],
                                 dest_sel=None)
440
            t_read += time.time() - t0
Damien Naudet's avatar
Damien Naudet committed
441

442
            if roiIndices is not None:
Damien Naudet's avatar
Damien Naudet committed
443 444 445 446
                cube = read_cube[xSlice, ySlice, zSlice]
            else:
                cube = read_cube

447
            t0 = time.time()
Damien Naudet's avatar
Damien Naudet committed
448
            cube[mask] = cube[mask] / weights
449 450
            t_mask = time.time() - t0

451 452 453
            t0 = time.time()

            z_sum = cube.sum(axis=0).sum(axis=0)
Damien Naudet's avatar
Damien Naudet committed
454 455 456
            cube_sum_z = cube.sum(axis=2)
            y_sum = cube_sum_z.sum(axis=0)
            x_sum = cube_sum_z.sum(axis=1)
Damien Naudet's avatar
Damien Naudet committed
457

458 459 460 461 462 463 464 465 466 467 468 469 470 471 472
            # Background subtraction
            if fit_background == BackgroundTypes.CONSTANT:
                # Shift data so that smallest value is 0
                for array in (z_sum, y_sum, x_sum):
                    array -= np.nanmin(array)

            elif fit_background == BackgroundTypes.LINEAR:
                # Simple linear background
                for array in (z_sum, y_sum, x_sum):
                    array -= np.linspace(
                        array[0], array[-1], num=len(array), endpoint=True)

            elif fit_background != BackgroundTypes.NONE:
                raise RuntimeError("Unsupported background subtraction")

473
            fitter.fit(i_fit, i_cube, x_sum, y_sum, z_sum)
Damien Naudet's avatar
Damien Naudet committed
474

Damien Naudet's avatar
Damien Naudet committed
475 476
            t_fit += time.time() - t0

Damien Naudet's avatar
Damien Naudet committed
477 478 479 480
            t0 = time.time()

            t_write = time.time() - t0

Damien Naudet's avatar
Damien Naudet committed
481
    except Exception as ex:
Damien Naudet's avatar
Damien Naudet committed
482
        print('EX', ex)
Damien Naudet's avatar
Damien Naudet committed
483

484
    times = (t_read, t_mask, t_fit, t_write)
Damien Naudet's avatar
Damien Naudet committed
485 486
    if disp_times:
        print('Thread {0} done ({1}).'.format(th_idx, times))
Damien Naudet's avatar
Damien Naudet committed
487
    return times