peak_fit.py 16 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
__authors__ = ["D. Naudet"]
__date__ = "01/06/2016"
__license__ = "MIT"

33
import logging
34 35 36
import time
import ctypes
import multiprocessing as mp
Damien Naudet's avatar
Damien Naudet committed
37
from threading import Thread
38
import multiprocessing.sharedctypes as mp_sharedctypes
39
from silx.math.fit import snip1d
40 41 42

import numpy as np

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

52 53 54 55

_logger = logging.getLogger(__name__)


56 57
disp_times = False

Damien Naudet's avatar
WIP  
Damien Naudet committed
58

59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
def background_estimation(mode, data):
    """Estimates a background of mode kind for data

    :param BackgroundTypes mode: The kind of background to compute
    :param numpy.ndarray data: The data array to process
    :return: The estimated background as an array with same shape as input
    :rtype: numpy.ndarray
    :raises ValueError: In case mode is not known
    """
    # Background subtraction
    if mode == BackgroundTypes.CONSTANT:
        # Shift data so that smallest value is 0
        return np.ones_like(data) * np.nanmin(data)

    elif mode == BackgroundTypes.LINEAR:
        # Simple linear background
        return np.linspace(data[0], data[-1], num=len(data), endpoint=True)

77 78 79 80
    elif mode == BackgroundTypes.SNIP:
        # Using snip background
        return snip1d(data, snip_width=len(data))

81 82 83 84 85 86 87
    elif mode == BackgroundTypes.NONE:
        return np.zeros_like(data)

    else:
        raise ValueError("Unsupported background mode")


Damien Naudet's avatar
Damien Naudet committed
88
class PeakFitter(Thread):
Damien Naudet's avatar
Damien Naudet committed
89
    """
Thomas Vincent's avatar
Thomas Vincent committed
90
    :param str qspace_f: path to the HDF5 file containing the qspace cubes
Damien Naudet's avatar
Damien Naudet committed
91

Thomas Vincent's avatar
Thomas Vincent committed
92
    :param FitTypes fit_type:
Damien Naudet's avatar
Damien Naudet committed
93 94 95 96

    :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
97
    :type indices: *optional* `array_like`
Damien Naudet's avatar
Damien Naudet committed
98

Thomas Vincent's avatar
Thomas Vincent committed
99 100
    :param Union[int,None] n_proc:
        Number of process to use. If None, the config value is used.
101 102

    :param BackgroundTypes background: The background subtraction to perform
Damien Naudet's avatar
Damien Naudet committed
103
    """
Damien Naudet's avatar
Damien Naudet committed
104

Damien Naudet's avatar
Damien Naudet committed
105 106 107 108
    READY, RUNNING, DONE, ERROR, CANCELED = __STATUSES = range(5)

    def __init__(self,
                 qspace_f,
Damien Naudet's avatar
Damien Naudet committed
109
                 fit_type=FitTypes.GAUSSIAN,
110
                 n_peaks=1,
Damien Naudet's avatar
Damien Naudet committed
111 112
                 indices=None,
                 n_proc=None,
113
                 roi_indices=None,
Thomas Vincent's avatar
Thomas Vincent committed
114
                 background=BackgroundTypes.NONE):
Damien Naudet's avatar
Damien Naudet committed
115 116 117 118 119
        super(PeakFitter, self).__init__()

        self.__results = None
        self.__thread = None
        self.__progress = 0
Damien Naudet's avatar
Damien Naudet committed
120
        self.__callback = None
Damien Naudet's avatar
Damien Naudet committed
121 122 123 124 125 126 127

        self.__status = self.READY

        self.__indices = None

        self.__qspace_f = qspace_f
        self.__fit_type = fit_type
128
        self.__background = background
129
        self.__n_peaks = n_peaks
Damien Naudet's avatar
Damien Naudet committed
130 131 132 133

        if n_proc:
            self.__n_proc = n_proc
        else:
Thomas Vincent's avatar
Thomas Vincent committed
134
            n_proc = self.__n_proc = config.DEFAULT_PROCESS_NUMBER
Damien Naudet's avatar
Damien Naudet committed
135 136 137 138 139

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

        if roi_indices is not None:
Damien Naudet's avatar
Damien Naudet committed
140
            self.__roi_indices = np.array(roi_indices[:])
Damien Naudet's avatar
Damien Naudet committed
141 142 143 144 145
        else:
            self.__roi_indices = None

        if fit_type not in FitTypes.ALLOWED:
            self.__set_status(self.ERROR)
146 147 148 149 150
            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))
151

Damien Naudet's avatar
Damien Naudet committed
152 153 154 155
        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
156

Damien Naudet's avatar
Damien Naudet committed
157
                n_points = qdata_shape[0]
158

Damien Naudet's avatar
Damien Naudet committed
159
                if indices is None:
Damien Naudet's avatar
Damien Naudet committed
160
                    indices = list(range(n_points))
Damien Naudet's avatar
Damien Naudet committed
161 162 163 164 165
                else:
                    indices = indices[:]
        except IOError:
            self.__set_status(self.ERROR)
            raise
Damien Naudet's avatar
Damien Naudet committed
166

Damien Naudet's avatar
Damien Naudet committed
167
        self.__indices = np.array(indices)
Damien Naudet's avatar
Damien Naudet committed
168 169 170 171 172 173 174

    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
175 176
    results = property(lambda self: self.__results)

Damien Naudet's avatar
Damien Naudet committed
177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
    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
208
        n_peaks = self.__n_peaks
Damien Naudet's avatar
Damien Naudet committed
209 210 211 212 213

        t_total = time.time()

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

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

Damien Naudet's avatar
Damien Naudet committed
217 218
        try:
            with QSpaceH5.QSpaceH5(qspace_f) as qspace_h5:
Damien Naudet's avatar
Damien Naudet committed
219
                with qspace_h5:
Damien Naudet's avatar
Damien Naudet committed
220 221 222 223 224 225 226 227 228 229 230 231
                    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
232
        if fit_type == FitTypes.GAUSSIAN:
Damien Naudet's avatar
Damien Naudet committed
233
            fit_class = GaussianFitter
234 235 236 237
            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
238
            fit_class = CentroidFitter
Damien Naudet's avatar
Damien Naudet committed
239
            shared_results = CentroidResults(n_points=n_indices)
240 241 242 243 244
        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
245 246
        else:
            raise RuntimeError('Unknown Fit Type')
Damien Naudet's avatar
Damien Naudet committed
247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287

        # 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
288
                                 fit_class,
Damien Naudet's avatar
Damien Naudet committed
289 290 291
                                 (n_indices, 9),
                                 idx_queue,
                                 qspace_f,
292 293
                                 read_lock,
                                 self.__background))
Damien Naudet's avatar
Damien Naudet committed
294 295

        if disp_times:
Thomas Vincent's avatar
Thomas Vincent committed
296
            class MyTimes(object):
Damien Naudet's avatar
Damien Naudet committed
297 298 299 300 301 302 303 304 305 306 307 308 309
                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
310
            res_times = MyTimes()
Damien Naudet's avatar
Damien Naudet committed
311 312 313 314 315 316 317 318
            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
319 320 321
            res = pool.apply_async(_fit_process,
                                   args=arg_list,
                                   callback=callback)
Damien Naudet's avatar
Damien Naudet committed
322 323 324
            res_list.append(res)

        # sending the image indices
325 326
        for i_fit, i_cube in enumerate(indices):
            idx_queue.put([i_fit, i_cube])
Damien Naudet's avatar
Damien Naudet committed
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:
337 338 339 340 341
            _logger.info('Total : {0}.'.format(t_total))
            _logger.info('Read {0}'.format(res_times.t_read))
            _logger.info('Mask {0}'.format(res_times.t_mask))
            _logger.info('Fit {0}'.format(res_times.t_fit))
            _logger.info('Write {0}'.format(res_times.t_write))
Damien Naudet's avatar
Damien Naudet committed
342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359

        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
360 361
                                                 q_z=q_z,
                                                 background_mode=self.__background)
Damien Naudet's avatar
Damien Naudet committed
362 363 364 365 366

        self.__results = fit_results

        self.__set_status(self.DONE)

Damien Naudet's avatar
Damien Naudet committed
367 368 369
        if self.__callback:
            self.__callback()

Damien Naudet's avatar
Damien Naudet committed
370
        return fit_results
Damien Naudet's avatar
WIP  
Damien Naudet committed
371 372


Damien Naudet's avatar
Damien Naudet committed
373
def _init_thread(shared_res_,
Damien Naudet's avatar
Damien Naudet committed
374
                 shared_prog_,
Damien Naudet's avatar
Damien Naudet committed
375
                 fit_class_,
Damien Naudet's avatar
Damien Naudet committed
376 377 378
                 result_shape_,
                 idx_queue_,
                 qspace_f_,
379 380
                 read_lock_,
                 background_):
Damien Naudet's avatar
Damien Naudet committed
381
    global shared_res, \
Damien Naudet's avatar
Damien Naudet committed
382
        shared_progress, \
Damien Naudet's avatar
Damien Naudet committed
383
        fit_class, \
Damien Naudet's avatar
Damien Naudet committed
384 385 386
        result_shape, \
        idx_queue, \
        qspace_f, \
387 388
        read_lock, \
        fit_background
Damien Naudet's avatar
Damien Naudet committed
389 390

    shared_res = shared_res_
Damien Naudet's avatar
Damien Naudet committed
391
    shared_progress = shared_prog_
Damien Naudet's avatar
Damien Naudet committed
392
    fit_class = fit_class_
Damien Naudet's avatar
Damien Naudet committed
393 394 395 396
    result_shape = result_shape_
    idx_queue = idx_queue_
    qspace_f = qspace_f_
    read_lock = read_lock_
397
    fit_background = background_
Damien Naudet's avatar
Damien Naudet committed
398

Damien Naudet's avatar
WIP  
Damien Naudet committed
399

Damien Naudet's avatar
Damien Naudet committed
400
def _fit_process(th_idx, roiIndices=None):
Damien Naudet's avatar
Damien Naudet committed
401
    print('Thread {0} started.'.format(th_idx))
Damien Naudet's avatar
Damien Naudet committed
402 403 404
    try:
        t_read = 0.
        t_fit = 0.
405
        t_mask = 0.
Damien Naudet's avatar
Damien Naudet committed
406

407
        l_shared_res = shared_res.local_copy()
Damien Naudet's avatar
Damien Naudet committed
408
        progress = np.frombuffer(shared_progress, dtype='int32')
409

Damien Naudet's avatar
Damien Naudet committed
410
        qspace_h5 = QSpaceH5.QSpaceH5(qspace_f)
Damien Naudet's avatar
Damien Naudet committed
411

Damien Naudet's avatar
Damien Naudet committed
412
        # Put this in the main thread
Damien Naudet's avatar
Damien Naudet committed
413 414 415 416 417 418 419
        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
420
        with qspace_h5:
Damien Naudet's avatar
Damien Naudet committed
421 422 423 424 425 426 427
            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
428

429
            if roiIndices is not None:
Damien Naudet's avatar
Damien Naudet committed
430 431 432 433 434
                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
435 436 437
            mask = np.where(histo > 0)
            weights = histo[mask]

Damien Naudet's avatar
Damien Naudet committed
438 439 440
        # read_lock.release()
        read_cube = np.ascontiguousarray(np.zeros(q_shape[1:]),
                                         dtype=q_dtype)
441

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

Damien Naudet's avatar
Damien Naudet committed
444 445
        while True:
            # TODO : timeout
Thomas Vincent's avatar
Thomas Vincent committed
446
            data = idx_queue.get()
447

Thomas Vincent's avatar
Thomas Vincent committed
448
            if data is None:
Damien Naudet's avatar
Damien Naudet committed
449
                break
450

Thomas Vincent's avatar
Thomas Vincent committed
451
            i_fit, i_cube = data
452 453

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

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

Damien Naudet's avatar
Damien Naudet committed
458
            t0 = time.time()
Damien Naudet's avatar
Damien Naudet committed
459
            with qspace_h5.qspace_dset_ctx() as dset:
Damien Naudet's avatar
Damien Naudet committed
460
                dset.read_direct(read_cube,
Damien Naudet's avatar
Damien Naudet committed
461 462
                                 source_sel=np.s_[i_cube],
                                 dest_sel=None)
463
            t_read += time.time() - t0
Damien Naudet's avatar
Damien Naudet committed
464

465
            if roiIndices is not None:
Damien Naudet's avatar
Damien Naudet committed
466 467 468 469
                cube = read_cube[xSlice, ySlice, zSlice]
            else:
                cube = read_cube

470
            t0 = time.time()
Damien Naudet's avatar
Damien Naudet committed
471
            cube[mask] = cube[mask] / weights
472 473
            t_mask = time.time() - t0

474 475 476
            t0 = time.time()

            z_sum = cube.sum(axis=0).sum(axis=0)
Damien Naudet's avatar
Damien Naudet committed
477 478 479
            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
480

Thomas Vincent's avatar
Thomas Vincent committed
481
            if fit_background != BackgroundTypes.NONE:
482 483 484
                # Background subtraction
                for array in (z_sum, y_sum, x_sum):
                    array -= background_estimation(fit_background, array)
485

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

Damien Naudet's avatar
Damien Naudet committed
488 489
            t_fit += time.time() - t0

Damien Naudet's avatar
Damien Naudet committed
490 491 492 493
            t0 = time.time()

            t_write = time.time() - t0

Damien Naudet's avatar
Damien Naudet committed
494
    except Exception as ex:
Damien Naudet's avatar
Damien Naudet committed
495
        print('EX', ex)
Damien Naudet's avatar
Damien Naudet committed
496

497
    times = (t_read, t_mask, t_fit, t_write)
Damien Naudet's avatar
Damien Naudet committed
498 499
    if disp_times:
        print('Thread {0} done ({1}).'.format(th_idx, times))
Damien Naudet's avatar
Damien Naudet committed
500
    return times