path: root/PyMca5/PyMcaPhysics/xas/XASStackBatch.py

#/*##########################################################################
#
# The PyMca X-Ray Fluorescence Toolkit
#
# Copyright (c) 2004-2015 European Synchrotron Radiation Facility
#
# This file is part of the PyMca X-ray Fluorescence Toolkit developed at
# the ESRF by the Software group.
#
# 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.
#
#############################################################################*/
__author__ = "V.A. Sole - ESRF Data Analysis"
__contact__ = "sole@esrf.fr"
__license__ = "MIT"
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
__doc__ = """
Module to process a stack of absorption spectra.
"""
import os
import numpy
import h5py
import posixpath
import logging
from PyMca5.PyMca import XASClass
from PyMca5.PyMcaIO import ConfigDict
import time


_logger = logging.getLogger(__name__)


class XASStackBatch(object):
    def __init__(self, analyzer=None):
        if analyzer is None:
            analyzer = XASClass.XASClass()
        self._analyzer = analyzer

    def setConfiguration(self, configuration):
        if "XASParameters" in configuration:
            self._analyzer.setConfiguration(configuration["XASParameters"])
        else:
            self._analyzer.setConfiguration(configuration)

    def setConfigurationFile(self, ffile):
        if not os.path.exists(ffile):
            raise IOError("File <%s> does not exists" % ffile)
        configuration = ConfigDict.ConfigDict()
        configuration.read(ffile)
        self.setConfiguration(configuration)

    def processMultipleSpectra(self, x, y,
                               xmin=None,
                               xmax=None,
                               configuration=None,
                               ysum=None,
                               weight=None,
                               mask=None,
                               directory=None,
                               name=None,
                               entry=None):
        """
        This method performs the actual work.

        :param x: 1D array containing the x axis (usually the channels) of the spectra.
        :param y: 3D array containing the spectra as [nrows, ncolumns, nchannels]
        :param weight: 0 Means no weight, 1 Use an average weight, 2 Individual weights (slow)
        :return: A dictionnary with the results as keys.
        """

        t0 = time.time()
        if configuration is not None:
            self._analyzer.setConfiguration(configuration)

        # read the current configuration
        config = self._analyzer.getConfiguration()

        #
        if weight is None:
            # dictated by the current configuration
            pass
        else:
            _logger.warning("WARNING: weight not handled yet")
        weightPolicy = 0 # no weight
        #weightPolicy = 1 # use average weight from the sum spectrum
        #weightPolicy = 2 # individual pixel weights (slow)
        if hasattr(x, "value"):
            # hdf5 dataset
            x = x.value

        if hasattr(y, "info") and hasattr(y, "data"):
            data = y.data
            mcaIndex = y.info.get("McaIndex", -1)
        else:
            data = y
            mcaIndex = -1

        if len(data.shape) != 3:
            txt = "For the time being only three dimensional arrays supported"
            raise IndexError(txt)
        if mcaIndex not in [-1, 2]:
            txt = "For the time being only mca arrays supported"
            raise IndexError(txt)
        firstSpectrum = None
        if ysum is not None:
            firstSpectrum = ysum
        if weightPolicy:
            # if the cumulated spectrum is present it should be better
            nRows = data.shape[0]
            nColumns = data.shape[1]
            nPixels =  nRows * nColumns
            if ysum is not None:
                firstSpectrum = ysum
            elif weightPolicy == 1:
                # we need to calculate the sum spectrum to derive the uncertainties
                totalSpectra = data.shape[0] * data.shape[1]
                jStep = min(5000, data.shape[1])
                ysum = numpy.zeros((data.shape[mcaIndex],), numpy.float)
                for i in range(0, data.shape[0]):
                    if i == 0:
                        chunk = numpy.zeros((data.shape[0], jStep), numpy.float)
                    jStart = 0
                    while jStart < data.shape[1]:
                        jEnd = min(jStart + jStep, data.shape[1])
                        ysum += data[i, jStart:jEnd, :].sum(axis=0, dtype=numpy.float)
                        jStart = jEnd
                firstSpectrum = ysum
            else:
                firstSpectrum = data[0, :, :].sum(axis=0, dtype=numpy.float)

        if firstSpectrum is None:
            firstSpectrum = data[0, 0, :]
        # TODO: Check if only one X and it is well behaved in order to
        # avoid unnecessary calculation on each spectrum
        self._analyzer.setSpectrum(x, firstSpectrum)
        # initialize the output arrays
        ddict = self._analyzer.processSpectrum()

        # initialize the arrays from the first results
        entry0 = "PyMcaResults"
        usedEnergy = ddict["Energy"]
        usedMu = ddict["Mu"]
        normalizedIdx = (ddict["NormalizedEnergy"] >= ddict["NormalizedPlotMin"]) & \
              (ddict["NormalizedEnergy"] <= ddict["NormalizedPlotMax"])
        normalizedSpectrumX = ddict["NormalizedEnergy"][normalizedIdx]
        normalizedSpectrumY = ddict["NormalizedMu"][normalizedIdx]
        exafsIdx = (ddict["EXAFSKValues"] >= ddict["KMin"]) & \
                   (ddict["EXAFSKValues"] <= ddict["KMax"])
        exafsSpectrumX = ddict["EXAFSKValues"][exafsIdx]
        exafsSpectrumY = ddict["EXAFSNormalized"][exafsIdx]
        xFT = ddict["FT"]["FTRadius"]
        yFT = ddict["FT"]["FTIntensity"]

        if directory is None:
            directory = os.getcwd()
        if name is None:
            name = "XAS_Result"
        fname = os.path.join(directory, name)
        if entry is None:
            entry = posixpath.join("xas_analysis")
        else:
            entry = posixpath.join(entry, "xas_analysis")
        if not fname.endswith(".h5"):
            fname = fname + ".h5"
        out = h5py.File(fname, "w")
        e0Path = posixpath.join(entry, "edge")
        jumpPath = posixpath.join(entry, "jump")
        spectrumXPath = posixpath.join(entry, "spectrum", "energy")
        spectrumYPath = posixpath.join(entry, "spectrum", "mu")
        normalizedXPath = posixpath.join(entry, "normalized", "energy")
        normalizedYPath = posixpath.join(entry, "normalized", "mu")
        exafsXPath = posixpath.join(entry, "exafs", "k")
        exafsYPath = posixpath.join(entry, "exafs", "signal")
        ftXPath = posixpath.join(entry, "FT", "Radius")
        ftYPath = posixpath.join(entry, "FT", "Intensity")
        ftImaginaryPath = posixpath.join(entry, "FT", "Imaginary")

        iXMin = 0
        iXMax = data.shape[-1] - 1
        e0 = out.require_dataset(e0Path,
                                 shape=data.shape[:-1],
                                 dtype=numpy.float32,
                                 chunks=None,
                                 compression=None)
        jump = out.require_dataset(jumpPath,
                                   shape=data.shape[:-1],
                                   dtype=numpy.float32,
                                   chunks=None,
                                   compression=None)
        shape = list(data.shape[:-1]) + [usedEnergy.size]
        spectrumX = out.require_dataset(spectrumXPath,
                                   shape=[usedEnergy.size],
                                   dtype=numpy.float32,
                                   chunks=None,
                                   compression=None)
        spectrumY = out.require_dataset(spectrumYPath,
                                   shape=shape,
                                   dtype=numpy.float32,
                                   chunks=None,
                                   compression=None)
        shape = list(data.shape[:-1]) + [normalizedSpectrumX.size]
        normalizedX = out.require_dataset(normalizedXPath,
                                   shape=[normalizedSpectrumX.size],
                                   dtype=numpy.float32,
                                   chunks=None,
                                   compression=None)
        normalizedY = out.require_dataset(normalizedYPath,
                                   shape=shape,
                                   dtype=numpy.float32,
                                   chunks=None,
                                   compression=None)
        shape = list(data.shape[:-1]) + [exafsSpectrumX.size]
        exafsX = out.require_dataset(exafsXPath,
                                     shape=[exafsSpectrumX.size],
                                     dtype=numpy.float32,
                                     chunks=None,
                                     compression=None)
        exafsY = out.require_dataset(exafsYPath,
                                     shape=shape,
                                     dtype=numpy.float32,
                                     chunks=None,
                                     compression=None)
        shape = list(data.shape[:-1]) + [xFT.size]
        ftX = out.require_dataset(ftXPath,
                                     shape=[xFT.size],
                                     dtype=numpy.float32,
                                     chunks=None,
                                     compression=None)
        ftY = out.require_dataset(ftYPath,
                                     shape=shape,
                                     dtype=numpy.float32,
                                     chunks=None,
                                     compression=None)
        ftImaginary = out.require_dataset(ftImaginaryPath,
                                     shape=shape,
                                     dtype=numpy.float32,
                                     chunks=None,
                                     compression=None)
        spectrumX[:] = ddict["Energy"]
        normalizedX[:] = ddict["NormalizedEnergy"][normalizedIdx]
        exafsX[:] = ddict["EXAFSKValues"][exafsIdx]
        ftX[:] = ddict["FT"]["FTRadius"]

        t0 = time.time()
        totalSpectra = data.shape[0] * data.shape[1]
        jStep = min(100, data.shape[1])
        if weightPolicy == 2:
            SVD = False
            sigma_b = None
        elif weightPolicy == 1:
            # the +1 is to prevent misbehavior due to weights less than 1.0
            sigma_b = 1 + numpy.sqrt(dummySpectrum)/nPixels
            SVD = True
        else:
            SVD = True
            sigma_b = None
        last_svd = None
        #for i in range(10):
        for i in range(0, data.shape[0]):
            #print(i)
            #chunks of nColumns spectra
            if i == 0:
                chunk = numpy.zeros((jStep,
                                     iXMax-iXMin+1),
                                     numpy.float)
            jStart = 0
            j = 0
            while jStart < data.shape[1]:
                jEnd = min(jStart + jStep, data.shape[1])
                #chunk[:,:(jEnd - jStart)] = data[i, jStart:jEnd, iXMin:iXMax+1].T
                spectra  = data[i, jStart:jEnd, iXMin:iXMax+1]
                nSpectra = spectra.shape[0]
                for spectrumNumber in range(nSpectra):
                    if mask is not None:
                        if mask[i, j] == 0:
                            continue
                    self._analyzer.setSpectrum(x, spectra[spectrumNumber])
                    ddict = self._analyzer.processSpectrum()
                    spectrumY[i, j] = ddict["Mu"]
                    e0[i, j] = ddict["Edge"]
                    jump[i, j] = ddict["Jump"]
                    #normalizedX[i, j] = ddict["NormalizedEnergy"][normalizedIdx]
                    normalizedY[i, j] = ddict["NormalizedMu"][normalizedIdx]
                    #exafsX[i, j] = ddict["EXAFSKValues"][exafsIdx]
                    exafsY[i, j] = ddict["EXAFSNormalized"][exafsIdx]
                    #ftX[i, j] = ddict["FT"]["FTRadius"]
                    ftY[i, j] = ddict["FT"]["FTIntensity"]
                    ftImaginary[i, j] = ddict["FT"]["FTImaginary"]
                    j +=1
                jStart = jEnd
        outputDict = {}
        outputDict["names"] = ["Jump", "Edge"]
        output = numpy.zeros((2, e0.shape[0], e0.shape[1]), dtype = e0.dtype)
        output[0, :] = jump.value
        output[1, :] = e0.value
        outputDict["images"] = output
        out.flush()
        out.close()

        t = time.time() - t0
        _logger.debug("First fit elapsed = %f", t)
        _logger.debug("Spectra per second = %f", data.shape[0]*data.shape[1]/float(t))
        t0 = time.time()
        return outputDict

if __name__ == "__main__":
    _logger.setLevel(logging.DEBUG)
    analyzer = XASClass.XASClass()
    instance = XASStackBatch(analyzer=analyzer)
    configurationFile = "test.ini"
    dataFile = h5py.File("testdata.h5", "r")
    for entry in dataFile:
        data = dataFile[entry]["data"]
        energy = dataFile[entry]["energy"]
        break
    instance.setConfigurationFile(configurationFile)
    instance.processMultipleSpectra(energy, data)
    dataFile.close()