.. algorithm:: .. summary:: .. relatedalgorithms:: .. properties:: Description ----------- The maximum entropy method (MEM) is used as a signal processing technique for reconstructing images from noisy data. It selects a single image from the many images which fit the data with the same value of the statistic, :math:`\chi^2`. The maximum entropy method selects from this *feasible set* of images, the one which has minimum information (maximum entropy). More specifically, the algorithm maximizes the entropy :math:`S\left(x\right)` subject to the constraint: .. math:: \chi^2 = \frac{1}{N'}\sum_m \frac{\left(d_m - d_m^c\right)^2}{\sigma_m^2} \leq \chi^2_{target} where :math:`d_m` are the experimental data, :math:`\sigma_m` the associated errors, :math:`d_m^c` the calculated or reconstructed data, :math:`\chi^2_{target}` is value of the input algorithm property ChiTargetOverN and :math:`N'` the number of measured data point. The image is a set of numbers :math:`\{x_0, x_1, \dots, x_{N-1}\}` related to the measured data via a 1D Fourier transform: .. math:: d_m = \frac{1}{N} \sum_{j=0}^{N-1} x_j e^{i 2\pi m j / N} where :math:`N` is the number of image points, which can be made different from :math:`N'` and is controlled by the input algorithm property ResolutionFactor, see examples further below. Note that even for real input data the reconstructed image can be complex, which means that both the real and imaginary parts will be taken into account for the calculations. This is the default behaviour, which can be changed by setting the input property *ComplexImage* to *False*. Note that the algorithm will fail to converge if the image is complex (i.e. the data does not satisfy Friedel's law) and this option is set to *False*. For this reason, it is recommended to use the default when no prior knowledge is available. The entropy is defined on the image :math:`\{x_j\}` as: .. math:: S = \sum_j \left[ \sqrt{x_j^2 + A^2} - x_j \sinh^{-1}\left(\frac{x_j}{A}\right) \right] or .. math:: S = -\sum_j x_j \left(\log(x_j/A)-1\right) where :math:`A` is a constant and the formula which is used depends on the input property *PositiveImage*: when it is set to *False* the first equation will be applied, whereas the latter expression will be used if this property is set to *True*. The sensitive of the reconstructed image to reconstructed image will vary depending on the data. In general a smaller value would preduce a sharper image. See section 4.7 in Ref. [1] for recommended strategy to selected :math:`A`. The implementation used to find the solution where the entropy is maximized subject to the constraint in :math:`\chi^2` follows the approach by Skilling & Bryan [2], which is an algorithm that is not explicitly using Lagrange multiplier. Instead, they construct a subspace from a set of *search directions* to approach the maximum entropy solution. Initially, the image :math:`x` is set to the flat background :math:`A` and the search directions are constructed using the gradients of :math:`S` and :math:`\chi^2`: .. math:: \mathbf{e}_1 = \left|x\right|\left(\nabla S\right) .. math:: \mathbf{e}_2 = \left|x\right|\left(\nabla \chi^2\right) where :math:`a\left(b\right)` stands for a componentwise multiplication between vectors :math:`a` and :math:`b`. The algorithm next uses a quadratic approximation to determine the increment :math:`\delta \mathbf{x}` that *moves* the image one step closer to the solution: .. math:: \mathbf{x} = \mathbf{x} + \delta \mathbf{x} If *DataLinearAdj* :math:`l_m` and *DataConstAdj* :math:`c_m` are set then the reconstructed data is adjusted to .. math:: d_m = \frac{l_m}{N} \sum_{j=0}^{N-1} x_j e^{i 2\pi m j / N} + c_m If *PerSpectrumImage* is set false, then chi squared is calculated for all spectra together to form a single image .. math:: \chi^2 = \frac{1}{N'}\sum_{m,k} \frac{\left(d_{m,k} - d_{m,k}^c\right)^2}{\sigma_{m,k}^2} \leq \chi^2_{target} then the reconstructed data are calculated from this single image and the adjustments are applied for each spectrum: .. math:: d_{m,k} = \frac{l_{m,k}}{N} \sum_{j=0}^{N-1} x_j e^{i 2\pi m j / N} + c_{m,k} Output Workspaces ----------------- There are four output workspaces: *ReconstructedImage* and *ReconstructedData* contain the image and calculated data respectively, and *EvolChi* and *EvolAngle* record the algorithm's evolution. *ReconstructedData* can be used to check that the calculated data approaches the experimental, measured data, and *ReconstructedImage* corresponds to its Fourier transform. If the input data are real they contain twice the original number of spectra, as the Fourier transform will be a complex signal in general. The real and imaginary parts are organized as follows (see Table 1): assuming your input workspace has :math:`M` spectra, the real part of the reconstructed spectrum :math:`s` corresponds to spectrum :math:`s` in the output workspaces *ReconstructedImage* and *ReconstructedData*, while the imaginary part can be found in spectrum :math:`s+M`. When the input data are complex, the input workspace is expected to have :math:`2M` spectra, where real and imaginary parts of input :math:`s`, are arranged in spectra :math:`(s, s+M)` respectively. In other words, *ReconstructedData* and *ReconstructedImage* will contain :math:`2M` spectra, where spectra :math:`(s, s+M)` correspond to the real and imaginary parts reconstructed from the input signal at :math:`(s, s+M)` in the input workspace (see Table 2 below). The workspaces *EvolChi* and *EvolAngle* record the evolution of :math:`\chi^2` and the angle (or non-parallelism) between :math:`\nabla S` and :math:`\nabla \chi^2`. Note that the algorithm runs until a solution is found, and that an image is considered to be a true maximum entropy solution when the following two conditions are met simultaneously: .. math:: \chi^2_{Target} - \chi^2 < \epsilon_1, \qquad \frac{1}{2} \left| \frac{\nabla S}{\left|\nabla S\right|} - \frac{\nabla \chi^2}{\left|\nabla \chi^2\right|} \right| < \epsilon_2 While one of these conditions is not satisfied the algorithm keeps running until it reaches the maximum number of iterations. When the maximum number of iterations is reached, the last reconstructed image and its corresponding calculated data are returned, even if they do not correspond to a true maximum entropy solution satisfying the conditions above. In this way, the user can check by inspecting the output workspaces whether the algorithm was evolving towards the correct solution or not. On the other hand, the user must always check the validity of the solution by inspecting *EvolChi* and *EvolAngle*, whose values will be set to zero once the true maximum entropy solution is found. .. table:: Table 1. Output workspaces for a real input workspace with M histograms and N bins, taking J iterations. +-------------------+------------------------------+----------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | Workspace | Number of histograms | Number of bins | Description | +===================+==============================+================+====================================================================================================================================================================================================================================================================================================================+ | EvolChi | M | J | For spectrum :math:`k` in the input workspace, evolution of :math:`\chi^2` until the solution is found | +-------------------+------------------------------+----------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | EvolAngle | M | J | For spectrum :math:`k` in the input workspace, evolution of the angle between :math:`\nabla S` and :math:`\nabla \chi^2`, until the solution is found | +-------------------+------------------------------+----------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | ReconstructedImage| 2M | N | For spectrum :math:`k` in the input workspace, the reconstructed image is stored in spectra :math:`k` (real part) and :math:`k+M` (imaginary part) | +-------------------+------------------------------+----------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ | ReconstructedData | 2M | N | For spectrum :math:`k` in the input workspace, the reconstructed data are stored in spectrum :math:`k` (real part) and :math:`k+M` (imaginary part). Note that although the input is real, the imaginary part is recorded for debugging purposes, it should be zero for all data points. | +-------------------+------------------------------+----------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ .. table:: Table 2. Output workspaces for a complex input workspace with 2M histograms and N bins, taking J iterations. +-------------------+------------------------------+----------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------+ | Workspace | Number of histograms | Number of bins | Description | +===================+==============================+================+==============================================================================================================================================================+ | EvolChi | M | J | For spectrum :math:`(k, k+M)` in the input workspace, evolution of :math:`\chi^2` until the solution is found | +-------------------+------------------------------+----------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------+ | EvolAngle | M | J | For spectrum :math:`(k, k+M)` in the input workspace, evolution of the angle between :math:`\nabla S` and :math:`\nabla \chi^2`, until the solution is found | +-------------------+------------------------------+----------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------+ | ReconstructedImage| 2M | :math:`N` | For spectrum :math:`(k, k+M)` in the input workspace, the reconstructed image is stored in spectra :math:`k` (real part) and :math:`k+M` (imaginary part) | +-------------------+------------------------------+----------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------+ | ReconstructedData | 2M | :math:`N` | For spectrum :math:`(k, k+M)` in the input workspace, the reconstructed data are stored in spectra :math:`k` (real part) and :math:`k+M` (imaginary part) | +-------------------+------------------------------+----------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------+ Usage ----- Some of the usage examples use usage data files. .. include:: ../usagedata-note.txt **Example - Reconstruct Fourier coefficients** In the example below, a workspace containing five Fourier coefficients is created and used as input to :ref:`algm-MaxEnt`. In the figure we show the original and reconstructed data (left), and the reconstructed image, i.e. Fourier transform (right). .. testcode:: ExFourierCoeffs # Create an empty workspace X = [] Y = [] E = [] N = 200 for i in range(0,N): x = ((i-N/2) *1./N) X.append(x) Y.append(0) E.append(0.001) # Fill in five Fourier coefficients # The input signal must be symmetric to get a real image Y[5] = Y[195] = 0.85 Y[10] = Y[190] = 0.85 Y[20] = Y[180] = 0.85 Y[12] = Y[188] = 0.90 Y[14] = Y[186] = 0.90 CreateWorkspace(OutputWorkspace='inputws',DataX=X,DataY=Y,DataE=E,NSpec=1) evolChi, evolAngle, image, data = MaxEnt(InputWorkspace='inputws', A=0.0001) print("First reconstructed coefficient: {:.3f}".format(data.readY(0)[5])) print("Second reconstructed coefficient: {:.3f}".format(data.readY(0)[10])) print("Third reconstructed coefficient: {:.3f}".format(data.readY(0)[20])) print("Fourth reconstructed coefficient: {:.3f}".format(data.readY(0)[12])) print("Fifth reconstructed coefficient: {:.3f}".format(data.readY(0)[14])) print("Number of iterations: "+str( len(evolAngle.readX(0)))) Output: .. testoutput:: ExFourierCoeffs :options: +ELLIPSIS +NORMALIZE_WHITESPACE First reconstructed coefficient: 0.847 Second reconstructed coefficient: 0.846 Third reconstructed coefficient: 0.846 Fourth reconstructed coefficient: 0.896 Fifth reconstructed coefficient: 0.896 Number of iterations: ... .. figure:: ../images/MaxEntFourierCoefficients.png :align: center **Example - Reconstruct a real muon dataset** In this example, :ref:`algm-MaxEnt` is run on a pre-analyzed muon dataset. The corresponding figure shows the original and reconstructed data (left), and the real part of the image obtained with :ref:`algm-MaxEnt` and :ref:`algm-FFT` (right). .. testcode:: ExMUSR00022725 Load(Filename=r'MUSR00022725.nxs', OutputWorkspace='MUSR00022725') CropWorkspace(InputWorkspace='MUSR00022725', OutputWorkspace='MUSR00022725', XMin=0.11, XMax=1.6, EndWorkspaceIndex=0) RemoveExpDecay(InputWorkspace='MUSR00022725', OutputWorkspace='MUSR00022725') Rebin(InputWorkspace='MUSR00022725', OutputWorkspace='MUSR00022725', Params='0.016') evolChi, evolAngle, image, data = MaxEnt(InputWorkspace='MUSR00022725', A=0.005) # Compare MaxEnt to FFT imageFFT = FFT(InputWorkspace='MUSR00022725') print("Image at {:.3f}: {:.3f}".format(image.readX(0)[44], image.readY(0)[44])) print("Image at {:.3f}: {:.3f}".format(image.readX(0)[46], image.readY(0)[46])) print("Image at {:.3f}: {:.3f}".format(image.readX(0)[48], image.readY(0)[48])) # check background is not nosiy def getInt(originalX,originalY,xMin,xMax): import numpy as np yData=[] xData=[] for j in range(0,len(originalX)): if originalX[j]>xMin and originalX[j]