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 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

Damien Naudet's avatar
Damien Naudet committed
52
class PeakFitter(Thread):
Damien Naudet's avatar
Damien Naudet committed
53 54 55 56 57 58 59 60 61 62 63 64
    """
    :param qspace_f: path to the HDF5 file containing the qspace cubes
    :type data_h5f: `str`

    :param fit_type:
    :type img_indices: *optional*

    :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.
    :type img_indices: *optional* `array_like`

Thomas Vincent's avatar
Thomas Vincent committed
65 66
    :param Union[int,None] n_proc:
        Number of process to use. If None, the config value is used.
Damien Naudet's avatar
Damien Naudet committed
67
    """
Damien Naudet's avatar
Damien Naudet committed
68

Damien Naudet's avatar
Damien Naudet committed
69 70 71 72
    READY, RUNNING, DONE, ERROR, CANCELED = __STATUSES = range(5)

    def __init__(self,
                 qspace_f,
Damien Naudet's avatar
Damien Naudet committed
73
                 fit_type=FitTypes.GAUSSIAN,
74
                 n_peaks=1,
Damien Naudet's avatar
Damien Naudet committed
75 76 77 78 79 80 81 82
                 indices=None,
                 n_proc=None,
                 roi_indices=None):
        super(PeakFitter, self).__init__()

        self.__results = None
        self.__thread = None
        self.__progress = 0
Damien Naudet's avatar
Damien Naudet committed
83
        self.__callback = None
Damien Naudet's avatar
Damien Naudet committed
84 85 86 87 88 89 90

        self.__status = self.READY

        self.__indices = None

        self.__qspace_f = qspace_f
        self.__fit_type = fit_type
91
        self.__n_peaks = n_peaks
Damien Naudet's avatar
Damien Naudet committed
92 93 94 95

        if n_proc:
            self.__n_proc = n_proc
        else:
Thomas Vincent's avatar
Thomas Vincent committed
96
            n_proc = self.__n_proc = config.DEFAULT_PROCESS_NUMBER
Damien Naudet's avatar
Damien Naudet committed
97 98 99 100 101

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

        if roi_indices is not None:
Damien Naudet's avatar
Damien Naudet committed
102
            self.__roi_indices = np.array(roi_indices[:])
Damien Naudet's avatar
Damien Naudet committed
103 104 105 106 107 108
        else:
            self.__roi_indices = None

        if fit_type not in FitTypes.ALLOWED:
            self.__set_status(self.ERROR)
            raise ValueError('Unknown fit type : {0}'.format(fit_type))
109

Damien Naudet's avatar
Damien Naudet committed
110 111 112 113
        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
114

Damien Naudet's avatar
Damien Naudet committed
115
                n_points = qdata_shape[0]
116

Damien Naudet's avatar
Damien Naudet committed
117
                if indices is None:
Damien Naudet's avatar
Damien Naudet committed
118
                    indices = list(range(n_points))
Damien Naudet's avatar
Damien Naudet committed
119 120 121 122 123
                else:
                    indices = indices[:]
        except IOError:
            self.__set_status(self.ERROR)
            raise
Damien Naudet's avatar
Damien Naudet committed
124

Damien Naudet's avatar
Damien Naudet committed
125
        self.__indices = np.array(indices)
Damien Naudet's avatar
Damien Naudet committed
126 127 128 129 130 131 132

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

Damien Naudet's avatar
Damien Naudet committed
135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
    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
166
        n_peaks = self.__n_peaks
Damien Naudet's avatar
Damien Naudet committed
167 168 169 170 171

        t_total = time.time()

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

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

Damien Naudet's avatar
Damien Naudet committed
175 176
        try:
            with QSpaceH5.QSpaceH5(qspace_f) as qspace_h5:
Damien Naudet's avatar
Damien Naudet committed
177
                with qspace_h5:
Damien Naudet's avatar
Damien Naudet committed
178 179 180 181 182 183 184 185 186 187 188 189
                    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
190
        if fit_type == FitTypes.GAUSSIAN:
Damien Naudet's avatar
Damien Naudet committed
191
            fit_class = GaussianFitter
192 193 194 195
            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
196
            fit_class = CentroidFitter
Damien Naudet's avatar
Damien Naudet committed
197
            shared_results = CentroidResults(n_points=n_indices)
198 199 200 201 202
        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)
Damien Naudet's avatar
Damien Naudet committed
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243

        # 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
244
                                 fit_class,
Damien Naudet's avatar
Damien Naudet committed
245 246 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
                                 (n_indices, 9),
                                 idx_queue,
                                 qspace_f,
                                 read_lock))

        if disp_times:
            class myTimes(object):
                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_

            res_times = myTimes()
            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
274 275 276
            res = pool.apply_async(_fit_process,
                                   args=arg_list,
                                   callback=callback)
Damien Naudet's avatar
Damien Naudet committed
277 278 279
            res_list.append(res)

        # sending the image indices
280 281
        for i_fit, i_cube in enumerate(indices):
            idx_queue.put([i_fit, i_cube])
Damien Naudet's avatar
Damien Naudet committed
282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320

        # 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,
                                                 q_z=q_z)

        self.__results = fit_results

        self.__set_status(self.DONE)

Damien Naudet's avatar
Damien Naudet committed
321 322 323
        if self.__callback:
            self.__callback()

Damien Naudet's avatar
Damien Naudet committed
324
        return fit_results
Damien Naudet's avatar
WIP  
Damien Naudet committed
325 326


Damien Naudet's avatar
Damien Naudet committed
327
def _init_thread(shared_res_,
Damien Naudet's avatar
Damien Naudet committed
328
                 shared_prog_,
Damien Naudet's avatar
Damien Naudet committed
329
                 fit_class_,
Damien Naudet's avatar
Damien Naudet committed
330 331 332 333
                 result_shape_,
                 idx_queue_,
                 qspace_f_,
                 read_lock_):
Damien Naudet's avatar
Damien Naudet committed
334
    global shared_res, \
Damien Naudet's avatar
Damien Naudet committed
335
        shared_progress, \
Damien Naudet's avatar
Damien Naudet committed
336
        fit_class, \
Damien Naudet's avatar
Damien Naudet committed
337 338 339 340 341 342
        result_shape, \
        idx_queue, \
        qspace_f, \
        read_lock

    shared_res = shared_res_
Damien Naudet's avatar
Damien Naudet committed
343
    shared_progress = shared_prog_
Damien Naudet's avatar
Damien Naudet committed
344
    fit_class = fit_class_
Damien Naudet's avatar
Damien Naudet committed
345 346 347 348
    result_shape = result_shape_
    idx_queue = idx_queue_
    qspace_f = qspace_f_
    read_lock = read_lock_
Damien Naudet's avatar
Damien Naudet committed
349

Damien Naudet's avatar
WIP  
Damien Naudet committed
350

Damien Naudet's avatar
Damien Naudet committed
351
def _fit_process(th_idx, roiIndices=None):
Damien Naudet's avatar
Damien Naudet committed
352
    print('Thread {0} started.'.format(th_idx))
Damien Naudet's avatar
Damien Naudet committed
353 354 355
    try:
        t_read = 0.
        t_fit = 0.
356
        t_mask = 0.
Damien Naudet's avatar
Damien Naudet committed
357

358
        l_shared_res = shared_res.local_copy()
Damien Naudet's avatar
Damien Naudet committed
359
        progress = np.frombuffer(shared_progress, dtype='int32')
360

Damien Naudet's avatar
Damien Naudet committed
361
        qspace_h5 = QSpaceH5.QSpaceH5(qspace_f)
Damien Naudet's avatar
Damien Naudet committed
362

Damien Naudet's avatar
Damien Naudet committed
363
        # Put this in the main thread
Damien Naudet's avatar
Damien Naudet committed
364 365 366 367 368 369 370
        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
371
        with qspace_h5:
Damien Naudet's avatar
Damien Naudet committed
372 373 374 375 376 377 378
            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
379

380
            if roiIndices is not None:
Damien Naudet's avatar
Damien Naudet committed
381 382 383 384 385
                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
386 387 388
            mask = np.where(histo > 0)
            weights = histo[mask]

Damien Naudet's avatar
Damien Naudet committed
389 390 391
        # read_lock.release()
        read_cube = np.ascontiguousarray(np.zeros(q_shape[1:]),
                                         dtype=q_dtype)
392

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

Damien Naudet's avatar
Damien Naudet committed
395 396
        while True:
            # TODO : timeout
397
            next = idx_queue.get()
398

399
            if next is None:
Damien Naudet's avatar
Damien Naudet committed
400
                break
401

402 403 404
            i_fit, i_cube = next

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

406
            if i_cube % 100 == 0:
Damien Naudet's avatar
Damien Naudet committed
407
                print(
408
                'Processing cube {0}/{1}.'.format(i_fit, result_shape[0]))
409

Damien Naudet's avatar
Damien Naudet committed
410
            t0 = time.time()
Damien Naudet's avatar
Damien Naudet committed
411
            with qspace_h5.qspace_dset_ctx() as dset:
Damien Naudet's avatar
Damien Naudet committed
412
                dset.read_direct(read_cube,
Damien Naudet's avatar
Damien Naudet committed
413 414
                                 source_sel=np.s_[i_cube],
                                 dest_sel=None)
415
            t_read += time.time() - t0
Damien Naudet's avatar
Damien Naudet committed
416

417
            if roiIndices is not None:
Damien Naudet's avatar
Damien Naudet committed
418 419 420 421
                cube = read_cube[xSlice, ySlice, zSlice]
            else:
                cube = read_cube

422
            t0 = time.time()
Damien Naudet's avatar
Damien Naudet committed
423
            cube[mask] = cube[mask] / weights
424 425
            t_mask = time.time() - t0

426 427 428
            t0 = time.time()

            z_sum = cube.sum(axis=0).sum(axis=0)
Damien Naudet's avatar
Damien Naudet committed
429 430 431
            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
432

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

Damien Naudet's avatar
Damien Naudet committed
435 436
            t_fit += time.time() - t0

Damien Naudet's avatar
Damien Naudet committed
437 438 439 440
            t0 = time.time()

            t_write = time.time() - t0

Damien Naudet's avatar
Damien Naudet committed
441
    except Exception as ex:
Damien Naudet's avatar
Damien Naudet committed
442
        print('EX', ex)
Damien Naudet's avatar
Damien Naudet committed
443

444
    times = (t_read, t_mask, t_fit, t_write)
Damien Naudet's avatar
Damien Naudet committed
445 446
    if disp_times:
        print('Thread {0} done ({1}).'.format(th_idx, times))
Damien Naudet's avatar
Damien Naudet committed
447
    return times
Damien Naudet's avatar
Damien Naudet committed
448 449 450 451


if __name__ == '__main__':
    pass