iFit: Mathematical operators

1. A reminder on axes, monitors and signal errors
2. Unary operators
   1. Fourier transform
   2. Wavelet transform
   3. Derivatives and variable change (Jacobian)
3. Binary operators
   1. Binary operators between objects of different type (Data sets, Models, ...)
   2. Convolution/correlation
   3. Comparisons
4. Statistics, fitting and peak searching, segmentation (k-means)
   1. Statistics
   2. Peak searching and model fitting
   3. Data segmentation (partitioning/clustering, k-means)
5. Projection, integration and sum
6. Cutting, Merging and combining data sets
   1. Cutting and selecting portions of a data set
   2. The difference between addition and combination
   3. Gluing an object content inside an other one
7. Appending and slicing data sets
   1. Graphical user interface to explore data sets and slices
   2. Graphical user interface to select regions-of-interest (ROI)
8. Creating arrays of objects
9. Interpolation, smoothing, intersection and union
   1. Intersection and union
   2. Interpolation
   3. Event data sets
   4. Re-binning and resizing
   5. Smoothing
   6. Manipulating object dimensions

A reminder on axes, monitors and signal errors

>> a = iData([ ifitpath 'Data/ILL_IN6.dat' ]);
>> a{0}; % or getaxis(a, 0) to get the Signal (0-rank)
>> a{1}; % or getaxis(a, 1) to get the 1-st rank axis, that is 'x' for vector-type data, and 'y' for 2D+
>> a{2}; % or getaxis(a, 2) to get the 2-nd rank axis, that is 'x' for 2D+ data sets

• All mathematical methods propagate the Error together with the Signal, assuming Error distributions are Gaussian. The sections below detail the way Errors are propagated.
• Binary operators use the Monitor value to weight the respective ratio when performing the operation. The Monitor value is also affected by Mathematical operations that renormalize the Signal. The Monitor is usually a value larger than 1. When set to 0 (or empty), no weighting is used.
  
• The Axes are used in binary operators to determine intersection of Signal ranges. A rebinning of objects is automatically performed when the Axes sampling do not match.
  

Unary operators

>> a = iData([ ifitpath 'Data/ILL_IN6.dat' ]);
>> b = [ log(a) floor(a) sqrt(a) ]
b =  array [1  3] iData object:

Index     [Tag]     [Dimension]                           [Title] [Last command]          [Label]
    1  tp362134     [340  1024] 'File ILL_IN6.dat ILL Data (n...' tp362134=log(tp36213...  
    2  tp362140     [340  1024] 'File ILL_IN6.dat ILL Data (n...' tp362140=floor(tp362...  
    3  tp362145     [340  1024] 'File ILL_IN6.dat ILL Data (n...' tp362145=sqrt(tp3621...  

Fourier transform

>> t=linspace(0,1,1000);
>> a = iData(t,0.7*sin(2*pi*50*t)+sin(2*pi*120*t)+0.05*randn(size(t)));
>> c=fft(a); d=ifft(c); subplot([ a abs(c) d ]);

Wavelet transform

>> a=iData([ ifitpath 'Data/Diff_BananaPSD_1314088587.th_y');
>> w=cwt(a, 2);	% the cwt along the axis rank 2 (scattering angle)
>> plot(a); slice(log(abs(w).^2))	% plot the initial data set and its cwt power spectrum

Derivatives and variable change (Jacobian)

>> a=iData(peaks);
>> g=gradient(a); subplot([ a g(1) g(2) ]);

Binary operators

• s1, e1 and m1 are the Signal, Error and Monitor of the 1st object, resp.
• s2, e2 and m2 are the Signal, Error and Monitor of the 2nd object, resp.
• s3, e3 and m3 are the Signal, Error and Monitor of the resulting object, resp.
  

>> a = iData([ ifitpath 'Data/ILL_IN6.dat' ]);
>> b = iData([ ifitpath 'Data/ILL_IN20.dat' ]); 
>> set(a,'Monitor', mean(b.Monitor)/10);
>> c = a+b;
>> subplot([ log(a) b log(c) ] ,'tight')

>> a = iData([ ifitpath 'Data/ILL_IN6.dat' ]);
>> c = a+[ ifitpath 'Data/ILL_IN20.dat' ];

Binary operators between objects of different type (Data sets, Models, ...)

Convolution/correlation

>> c = conv(a,b);		% returned convoluted object with size which is size(a)+size(b)+1
>> c = conv(a,b, 'same');	% returned convoluted object with size which is size(a)
>> c = conv(a,b, 'valid');	% returned convoluted object with size which is size(a)-size(b)+1
>> c = conv(a,b, 'pad');	% pads 'a' with starting/ending values to minimize border effects
>> c = conv(a,b, 'center');	% centers 'b' so that convolution does not shift 'a' axes
>> c = conv(a,b, 'norm');	% normalizes 'b' so that convolution does not change 'a' integral
>> c = conv(a,b, 'background');	% subtracts minimal value in 'b' so that convolution does not change 'a' integral
>> c = conv(a,b, 'deconv');	% deconvolution, but very sensitive to noise (use with caution)

>> c = conv(a,b, 'same pad background center norm');
>> c = convn(a,b);	% same as above in a shorter call

>> a = iData([ ifitpath 'Data/ILL_IN6.dat' ]);
>> c = convn(a,a);

>> c = convn(a, 3)	% convolution of an iData with a normalized Gaussian of width 3

>> global b;			% holds the convolution function as a single vector/matrix or iData object
>> h = convn(gauss, 'double(b)');	% creates a fit function with convolution
>> h.Constraint= 'global b'; 	% declares b as an external global variable
>> p = fits(a, h);

>> c = interp(b, a);

Comparisons

• greater: gt >
• greater or equal: ge >=
• lower: lt <
• lower or equal: le <=
• equal: eq ==
• not equal: ne ~=

Statistics, fitting and peak searching, segmentation (k-means)

Statistics

>> a = iData([ ifitpath 'Data/sv1850.scn' ]);
>> [w,x]=std(a)
w =
    0.0036
x =
    1.0007
>> [ min(a) max(a) median(a) mean(a) ]
ans =

            0         7387          119       1630.7

>> fits(a)
   'Amplitude'    'Centre'    'HalfWidth'    'Background'
      0.67857       1.0008    0.0034641   0.00024067

>> a = iData(peaks);
>> [w, x] = std(a, -(1:ndims(a))); % we subtract the minimum value (negative dimension spec.)

Peak searching and model fitting

>> a = iData([ ifitpath 'Data/MCA.dat' ]);
>> setalias(a,'Signal','Data.Counts(:,2)');
>> [half_width, center, amplitude, baseline]=peaks(a);
>> plot(a); hold on; 
>> plot(center,amplitude,'ro'); set(gca, 'yscale', 'log')

Data segmentation (partitioning/clustering, k-means)

>> a = iData(peaks);
>> b = kmeans(a, 3);	% request segmentation in 3 groups
>> plot(a, b+5);	% overlay data set and segmentation result

Projection, integration and sum

>> a = iData([ ifitpath 'Data/ILL_IN6.dat' ]);
>> xlabel(a, 'Time channel'); % 2nd axis
>> ylabel(a, 'Angle channel');% 1st axis
>> subplot([ log(a) log(camproj(a)) log(sum(a)) ],'tight');

>> subplot([ (a) cumsum(a) ] ,'tight')

>> a = iData([ ifitpath 'Data/*.cbf' ]);
>> b = camproj(a,'radial'); % returns a distribution to be transformed into an histogram
>> c = cart2sph(a);
>> subplot(log(a), hist(b,500), log(c))

>> b = iData(flow);		% is 3D
>> resize(b, [25 50 1]);	% makes a 2D object from a 3D one

Cutting, Merging and combining data sets

Cutting and selecting portions of a data set

>> a = iData([ ifitpath 'Data/ILL_IN6.dat' ]);
>> b = xlim(a, [100 200]);					% selects the object with X axis ranging from 100 to 200
>> b = a(:, find(  100 <= a{2} & a{2} <= 200 ));				% same as above: X is rank 2(columns)
>> c = a(find(50 <= a{1} & a{1} <= 150),find(100 <= a{2} & a{2} <= 200));	% cuts along X(rank2) and Y axes(rank1, rows)

The difference between addition and combination

• combination: (S1+S2) over monitor (M1+M2)
• addition (sum): (M1+M2)*(S1/M1+S2/M2) over monitor(M1+M2)

>> a = iData([ ifitpath 'Data/ILL_IN6*.dat' ]);
>> a(1)=setalias(a(1),'Monitor', 1);
>> a(2)=setalias(a(2),'Monitor', 10);
>> b=combine(a)             % combine/merge
>> c=a(1)+a(2);             % addition
>> subplot([ log(a) log(b) log(c) ], 'tight');

Gluing an object content inside an other one

>> b(1) = a;

Appending and slicing data sets

>> x=-pi:0.01:pi; a=iData(x,x); 
>> a.Error=0;                         % replace default Error=sqrt(Signal) by no-error.
>> b=sin(a); c=cos(a); d=exp(-a.*a);  % create new objects by applying operator on the initial linear one
>> e=cat(1, [a b c d ]);              % catenate along all X axes, result is a single vector
>> f=copyobj(e);
>> rmaxis(f,1);                       % reset object axis to the default
>> g=cat(2, [a b c d]);               % append objects side-by-side onto new dimension
>> subplot([e f g],'tight')           % plot all
>> h=dog(2, g);                       % back to [a b c d] with 'dog'

>> a = iData(peaks);
>> b = a(10, :);		% row 10
>> c = a([ 5 7 ], 1:2:end);	% rows 5 and 7, one every 2 columns

Graphical user interface to explore data sets and slices

Graphical user interface to select regions-of-interest (ROI)

Creating arrays of objects

>> a = zeros(iData, [5 3]);    % creates a 5x3 empty iData array
>> a = iData(peaks);
>> b = zeros(a, 5, 3);         % creates a 5x3 constant array of 'peaks' (2D surface)

>> a = iData(peaks);
>> b = linspace(a, cos(a), 5);
>> c = logspace(a, sin(a), 5);
>> subplot([b c],[2 5])

Linspace (first row) and logspace (2nd row)

Interpolation, smoothing, intersection and union

Intersection and union

>> a = iData(peaks);
>> b = copyobj(a);
>> a{1} = a{1}+10; a{2} = a{2}+10;  % shift along X and Y
>> a.Signal=a.Signal+5;             % shift along Z
>> [ai,bi]=intersect(a,b);
>> [au,bu]=union(a,b);
>> subplot(1,3,1); plot([a  b ],'transparent tight')
>> subplot(1,3,2); plot([ai bi],'transparent tight')
>> subplot(1,3,3); plot([au bu],'transparent tight')

Interpolation

>> a = iData(peaks(10))+2
>> b = interp(a,2)	% doubled binning, but we rather recommend the 'resize' method (see below)
>> c = interp(a,1:.25:15,3:.25:12)
>> subplot([a b c],[1 3],'mesh tight')

Event data sets

It is possible to create an event data set from a column-type object (that is a matrix which columns are axes/data) by means of the 'event' method (which creates a cloud of points) or the 'hist' method (which creates an histogram out of an event list).

For event based data sets, the hist operator transforms the event object into an histogrammed data set. The histogram axes can be given as vectors (where event will be accumulated), or a number of bins within the [min max] of the distributions. The hist operator is much faster than the interpolation one, interp for event based data.

>> a=iData([ ifitpath 'Data/Monitor_GV*']);
>> b=hist(a);
>> surf(b, 'median');

If the number of events is limited compared to the number of bins in the histogram, some areas of the histogram can remain empty when filling it with the events. It is then possible to 'fill' the empty slots with a fast neighboring algorithm, using a discrete cosine transform.

>> a=iData([ ifitpath 'Data/Monitor_GV*']);
>> b=hist(a); 
>> c=fill(b);
>> subplot(a,b,c);

Re-binning and resizing

To compress an object and reduce its binning, that is the number of elements on each dimension, you can use the reducevolume method. Without argument, it reduces the number of elements to less than 1e6. Additional arguments indicate the compression factor on each dimension. This method is automatically used when plotting large objects (except if the plot options contains a 'full' or 'whole' word).

In a similar manner, the resize method is a very fast yet accurate algorithm (discrete cosine transform) to make a fast rebinning, that is change the number of elements on each dimension/rank, but preserve the data set meaning.

>> a = iData(peaks(40))+2
>> subplot(a, resize(a, [20 20]), resize(a, 40, 40)) % resize operation

This allows for instance to quickly reduce the dimensionality of a data set:

>> b = iData(flow);
>> resize(b, [25 50 1]);	% makes a 2D object from a 3D one

Reference: Narasimha M. et al, On the computation of the discrete cosine transform, IEEE Trans Comm, 26, 6, 1978, pp 934-936.

Smoothing

To explicitly smooth a data set, you can use the 'smooth' method, which makes use of a discrete cosine transform filter spline smoothing (shape preserving):

>> a = iData(peaks(40).*(1+0.1*randn(40)));
>> b = smooth(a);
>> subplot(a,b);

Reference: Garcia D, Robust smoothing of gridded data in one and higher dimensions with missing values. Computational Statistics & Data Analysis, 2010.

The smoothing can also use a Savitsky-Golay moving average, but may lead to significant sharp features modifications (be cautious).

>> b = smooth(a, 'sgolay');
>> b = smooth(a, span, 'sgolay', degree); % Savitsky Golay degree and span (length) specified

 Last, you can use a nD Gaussian blur

>> b = smooth(a, 'gauss',1);		% Gaussian smoothing, all width=1
>> b = smooth(a, 'gauss',[ 1 10 ]);	% idem, with specific width, here for a 2D smoothing

Last, in the case a data set contains NaN (Not-a-Number) values (that is missing values), you can fill these values with a fast interpolation on neighboring elements using the fill method (based on a discrete cosine transform).

Manipulating object dimensions

It is also possible to manipulate an object Signal in the same way as a matrix with methods reshape (re-order elements in the object), permute (permute ranks, that is a generalized transposition) and circshift (shift ranks/dimensions), which work just as their matrix equivalent.

Last, the sort operator sorts object axes, in ascending or descending order. The corresponding Signal, Error, and Monitor values are also re-ordered accordingly.

E. Farhi - iFit/Math operators -  Nov. 27, 2018 2.0.2 - back to Main iFit Page