peak_fit.py 14.2 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

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

51 52 53 54

_logger = logging.getLogger(__name__)


55 56
disp_times = False

Damien Naudet's avatar
WIP  
Damien Naudet committed
57

58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
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)

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

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

    else:
        raise ValueError("Unsupported background mode")


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

Thomas Vincent's avatar
Thomas Vincent committed
91
    :param FitTypes fit_type:
Damien Naudet's avatar
Damien Naudet committed
92 93

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

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

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

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

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

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

        self.__status = self.READY

        self.__indices = None

        self.__qspace_f = qspace_f
        self.__fit_type = fit_type
127
        self.__background = background
128
        self.__n_peaks = n_peaks
Damien Naudet's avatar
Damien Naudet committed
129

130
        self.__n_proc = n_proc if n_proc else config.DEFAULT_PROCESS_NUMBER
Damien Naudet's avatar
Damien Naudet committed
131
        self.__shared_progress = mp_sharedctypes.RawArray(ctypes.c_int32,
132
                                                          self.__n_proc)
Damien Naudet's avatar
Damien Naudet committed
133 134

        if roi_indices is not None:
Damien Naudet's avatar
Damien Naudet committed
135
            self.__roi_indices = np.array(roi_indices[:])
Damien Naudet's avatar
Damien Naudet committed
136 137 138 139 140
        else:
            self.__roi_indices = None

        if fit_type not in FitTypes.ALLOWED:
            self.__set_status(self.ERROR)
141 142 143 144 145
            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))
146

Damien Naudet's avatar
Damien Naudet committed
147 148 149 150
        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
151

Damien Naudet's avatar
Damien Naudet committed
152
                n_points = qdata_shape[0]
153

Damien Naudet's avatar
Damien Naudet committed
154
                if indices is None:
Damien Naudet's avatar
Damien Naudet committed
155
                    indices = list(range(n_points))
Damien Naudet's avatar
Damien Naudet committed
156 157 158 159 160
                else:
                    indices = indices[:]
        except IOError:
            self.__set_status(self.ERROR)
            raise
Damien Naudet's avatar
Damien Naudet committed
161

Damien Naudet's avatar
Damien Naudet committed
162
        self.__indices = np.array(indices)
Damien Naudet's avatar
Damien Naudet committed
163 164 165 166 167 168 169

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

172
    def peak_fit(self, blocking=True, callback=None):
Damien Naudet's avatar
Damien Naudet committed
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194
        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)

        t_total = time.time()

195
        progress = np.frombuffer(self.__shared_progress, dtype='int32')
Damien Naudet's avatar
Damien Naudet committed
196
        progress[:] = 0
197

198
        n_indices = len(self.__indices)
Damien Naudet's avatar
Damien Naudet committed
199

Damien Naudet's avatar
Damien Naudet committed
200
        try:
201
            with QSpaceH5.QSpaceH5(self.__qspace_f) as qspace_h5:
Damien Naudet's avatar
Damien Naudet committed
202
                with qspace_h5:
203 204
                    x_pos = qspace_h5.sample_x[self.__indices]
                    y_pos = qspace_h5.sample_y[self.__indices]
Damien Naudet's avatar
Damien Naudet committed
205 206 207 208
        except IOError:
            self.__set_status(self.ERROR)
            raise

209
        if self.__fit_type == FitTypes.GAUSSIAN:
Damien Naudet's avatar
Damien Naudet committed
210
            fit_class = GaussianFitter
211
            shared_results = GaussianResults(n_points=n_indices,
212 213
                                             n_peaks=max(1, self.__n_peaks))
        elif self.__fit_type == FitTypes.CENTROID:
Damien Naudet's avatar
Damien Naudet committed
214
            fit_class = CentroidFitter
Damien Naudet's avatar
Damien Naudet committed
215
            shared_results = CentroidResults(n_points=n_indices)
216
        elif self.__fit_type == FitTypes.SILX:
217 218
            fit_class = SilxFitter
            shared_results = SilxResults(n_points=n_indices,
219
                                         n_peaks=max(1, self.__n_peaks))
Thomas Vincent's avatar
Thomas Vincent committed
220 221
        else:
            raise RuntimeError('Unknown Fit Type')
Damien Naudet's avatar
Damien Naudet committed
222 223 224 225 226 227

        manager = mp.Manager()

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

228
        pool = mp.Pool(self.__n_proc,
Damien Naudet's avatar
Damien Naudet committed
229 230
                       initializer=_init_thread,
                       initargs=(shared_results,
231
                                 self.__shared_progress,
Damien Naudet's avatar
Damien Naudet committed
232
                                 fit_class,
Damien Naudet's avatar
Damien Naudet committed
233 234
                                 (n_indices, 9),
                                 idx_queue,
Thomas Vincent's avatar
Thomas Vincent committed
235
                                 self.__qspace_f,
236 237
                                 read_lock,
                                 self.__background))
Damien Naudet's avatar
Damien Naudet committed
238 239

        if disp_times:
Thomas Vincent's avatar
Thomas Vincent committed
240
            class MyTimes(object):
Damien Naudet's avatar
Damien Naudet committed
241 242 243 244 245 246 247 248 249 250 251 252 253
                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
254
            res_times = MyTimes()
Damien Naudet's avatar
Damien Naudet committed
255 256 257 258 259 260
            callback = res_times.update
        else:
            callback = None

        # creating the processes
        res_list = []
261 262
        for th_idx in range(self.__n_proc):
            arg_list = (th_idx, self.__roi_indices)
Damien Naudet's avatar
Damien Naudet committed
263 264 265
            res = pool.apply_async(_fit_process,
                                   args=arg_list,
                                   callback=callback)
Damien Naudet's avatar
Damien Naudet committed
266 267 268
            res_list.append(res)

        # sending the image indices
269
        for i_fit, i_cube in enumerate(self.__indices):
270
            idx_queue.put([i_fit, i_cube])
Damien Naudet's avatar
Damien Naudet committed
271 272

        # sending the None value to let the threads know that they should return
273
        for th_idx in range(self.__n_proc):
Damien Naudet's avatar
Damien Naudet committed
274 275 276 277 278 279 280
            idx_queue.put(None)

        pool.close()
        pool.join()

        t_total = time.time() - t_total
        if disp_times:
281 282 283 284 285
            _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
286

Thomas Vincent's avatar
Thomas Vincent committed
287
        with QSpaceH5.QSpaceH5(self.__qspace_f) as qspace_h5:
288
            q_dim0, q_dim1, q_dim2 = qspace_h5.qspace_dimension_values
Damien Naudet's avatar
Damien Naudet committed
289

290 291 292 293
        if self.__roi_indices is not None:
            q_dim0 = q_dim0[self.__roi_indices[0][0]:self.__roi_indices[0][1]]
            q_dim1 = q_dim1[self.__roi_indices[1][0]:self.__roi_indices[1][1]]
            q_dim2 = q_dim2[self.__roi_indices[2][0]:self.__roi_indices[2][1]]
Damien Naudet's avatar
Damien Naudet committed
294 295 296

        fit_results = shared_results.fit_results(sample_x=x_pos,
                                                 sample_y=y_pos,
297 298 299
                                                 q_x=q_dim0,
                                                 q_y=q_dim1,
                                                 q_z=q_dim2,
Thomas Vincent's avatar
Thomas Vincent committed
300
                                                 background_mode=self.__background)
Damien Naudet's avatar
Damien Naudet committed
301 302 303 304 305

        self.__results = fit_results

        self.__set_status(self.DONE)

Damien Naudet's avatar
Damien Naudet committed
306 307 308
        if self.__callback:
            self.__callback()

Damien Naudet's avatar
Damien Naudet committed
309
        return fit_results
Damien Naudet's avatar
WIP  
Damien Naudet committed
310 311


Damien Naudet's avatar
Damien Naudet committed
312
def _init_thread(shared_res_,
Damien Naudet's avatar
Damien Naudet committed
313
                 shared_prog_,
Damien Naudet's avatar
Damien Naudet committed
314
                 fit_class_,
Damien Naudet's avatar
Damien Naudet committed
315 316 317
                 result_shape_,
                 idx_queue_,
                 qspace_f_,
318 319
                 read_lock_,
                 background_):
Damien Naudet's avatar
Damien Naudet committed
320
    global shared_res, \
Damien Naudet's avatar
Damien Naudet committed
321
        shared_progress, \
Damien Naudet's avatar
Damien Naudet committed
322
        fit_class, \
Damien Naudet's avatar
Damien Naudet committed
323 324 325
        result_shape, \
        idx_queue, \
        qspace_f, \
326 327
        read_lock, \
        fit_background
Damien Naudet's avatar
Damien Naudet committed
328 329

    shared_res = shared_res_
Damien Naudet's avatar
Damien Naudet committed
330
    shared_progress = shared_prog_
Damien Naudet's avatar
Damien Naudet committed
331
    fit_class = fit_class_
Damien Naudet's avatar
Damien Naudet committed
332 333 334 335
    result_shape = result_shape_
    idx_queue = idx_queue_
    qspace_f = qspace_f_
    read_lock = read_lock_
336
    fit_background = background_
Damien Naudet's avatar
Damien Naudet committed
337

Damien Naudet's avatar
WIP  
Damien Naudet committed
338

Damien Naudet's avatar
Damien Naudet committed
339
def _fit_process(th_idx, roiIndices=None):
Damien Naudet's avatar
Damien Naudet committed
340
    print('Thread {0} started.'.format(th_idx))
Damien Naudet's avatar
Damien Naudet committed
341 342 343
    try:
        t_read = 0.
        t_fit = 0.
344
        t_mask = 0.
Damien Naudet's avatar
Damien Naudet committed
345

346
        l_shared_res = shared_res.local_copy()
Damien Naudet's avatar
Damien Naudet committed
347
        progress = np.frombuffer(shared_progress, dtype='int32')
348

Damien Naudet's avatar
Damien Naudet committed
349
        qspace_h5 = QSpaceH5.QSpaceH5(qspace_f)
Damien Naudet's avatar
Damien Naudet committed
350

Damien Naudet's avatar
Damien Naudet committed
351
        # Put this in the main thread
Damien Naudet's avatar
Damien Naudet committed
352
        if roiIndices is not None:
353 354 355
            dim0Slice = slice(roiIndices[0][0], roiIndices[0][1], 1)
            dim1Slice = slice(roiIndices[1][0], roiIndices[1][1], 1)
            dim2Slice = slice(roiIndices[2][0], roiIndices[2][1], 1)
Damien Naudet's avatar
Damien Naudet committed
356 357 358

        # TODO : timeout to check if it has been canceled
        # read_lock.acquire()
Damien Naudet's avatar
Damien Naudet committed
359
        with qspace_h5:
360 361
            q_dim0, q_dim1, q_dim2 = qspace_h5.qspace_dimension_values

Damien Naudet's avatar
Damien Naudet committed
362 363 364 365
            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
366

367
            if roiIndices is not None:
368 369 370 371
                q_dim0 = q_dim0[dim0Slice]
                q_dim1 = q_dim1[dim1Slice]
                q_dim2 = q_dim2[dim2Slice]
                histo = histo[dim0Slice, dim1Slice, dim2Slice]
Damien Naudet's avatar
Damien Naudet committed
372

Damien Naudet's avatar
Damien Naudet committed
373 374 375
            mask = np.where(histo > 0)
            weights = histo[mask]

Damien Naudet's avatar
Damien Naudet committed
376 377 378
        # read_lock.release()
        read_cube = np.ascontiguousarray(np.zeros(q_shape[1:]),
                                         dtype=q_dtype)
379

380
        fitter = fit_class(q_dim0, q_dim1, q_dim2, l_shared_res)
Damien Naudet's avatar
Damien Naudet committed
381

Damien Naudet's avatar
Damien Naudet committed
382 383
        while True:
            # TODO : timeout
Thomas Vincent's avatar
Thomas Vincent committed
384
            data = idx_queue.get()
385

Thomas Vincent's avatar
Thomas Vincent committed
386
            if data is None:
Damien Naudet's avatar
Damien Naudet committed
387
                break
388

Thomas Vincent's avatar
Thomas Vincent committed
389
            i_fit, i_cube = data
390 391

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

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

Damien Naudet's avatar
Damien Naudet committed
396
            t0 = time.time()
Damien Naudet's avatar
Damien Naudet committed
397
            with qspace_h5.qspace_dset_ctx() as dset:
Damien Naudet's avatar
Damien Naudet committed
398
                dset.read_direct(read_cube,
Damien Naudet's avatar
Damien Naudet committed
399 400
                                 source_sel=np.s_[i_cube],
                                 dest_sel=None)
401
            t_read += time.time() - t0
Damien Naudet's avatar
Damien Naudet committed
402

403
            if roiIndices is not None:
404
                cube = read_cube[dim0Slice, dim1Slice, dim2Slice]
Damien Naudet's avatar
Damien Naudet committed
405 406 407
            else:
                cube = read_cube

408
            t0 = time.time()
Damien Naudet's avatar
Damien Naudet committed
409
            cube[mask] = cube[mask] / weights
410 411
            t_mask = time.time() - t0

412 413 414
            t0 = time.time()

            z_sum = cube.sum(axis=0).sum(axis=0)
Damien Naudet's avatar
Damien Naudet committed
415 416 417
            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
418

Thomas Vincent's avatar
Thomas Vincent committed
419
            if fit_background != BackgroundTypes.NONE:
420 421 422
                # Background subtraction
                for array in (z_sum, y_sum, x_sum):
                    array -= background_estimation(fit_background, array)
423

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

Damien Naudet's avatar
Damien Naudet committed
426 427
            t_fit += time.time() - t0

Damien Naudet's avatar
Damien Naudet committed
428 429 430 431
            t0 = time.time()

            t_write = time.time() - t0

Damien Naudet's avatar
Damien Naudet committed
432
    except Exception as ex:
Damien Naudet's avatar
Damien Naudet committed
433
        print('EX', ex)
Damien Naudet's avatar
Damien Naudet committed
434

435
    times = (t_read, t_mask, t_fit, t_write)
Damien Naudet's avatar
Damien Naudet committed
436 437
    if disp_times:
        print('Thread {0} done ({1}).'.format(th_idx, times))
Damien Naudet's avatar
Damien Naudet committed
438
    return times