function tensorHc = convertMeasureTensorSe(tensorFull, indexVector, ...
    mapping, dim)
% convert the tensor of mixed volumes for application to hC-vectors
%
% THIS FUNCTION IS NO USER FUNCTION.

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.

% substitute nan with zero
tensorFull(isnan(tensorFull)) = 0;

switch dim
    case 1
        % mapping:   nC x nC
        % tensorFull: 1 x nS
        % tensor:     1 x nC
        % tensorHc:   1 x nC
        tensor = tensorFull(indexVector);
        
        % meas = sum(i): mixV(i) * lam(i)
        % meas = sum(i): mixV(i) * (sum(q): mapping(i, q) * hc(q))
        % meas = sum(q): hc(j) * (sum(i): mapping(i, q) * mixV(i))
        % meas = sum(q): hc(q) * tensor(q)
        tensorHc = tensor * mapping;
        
    case 2
        % mapping:    nC x nC
        % tensorFull: nS x nS
        % tensor:     nC x nC
        % tensorHc:   nC x nC
        tensor = tensorFull(indexVector, indexVector);
        nCol = size(mapping, 2);
        
        % meas = sum(i,j): mixV(i,j) * lam(i) * lam(j)
        % meas = sum(i,j): mixV(i,j) * (sum(q): mapping(i, q) * hc(q))
        %                            * (sum(r): mapping(j, r) * hc(r))
        % meas = sum(q,r): hc(q) * hc(r)
        %      * sum(i,j): mapping(i, q) * mixV(i,j) *  mapping(j, r)
        % meas = sum(q,r): hc(q) * hc(r) * tensor(q,r)
        tensorHc = mapping' * tensor * mapping;
        
        % to make the hC-tensors comparable, we must see that the
        %   application is for i~=j the effective value is:
        %      tensor(i,j) + tensor(j,i).
        %   so that tensors with j<i are forced to 0 with the following
        filterLower = logical(tril(ones(nCol)) - eye(nCol));
        tensorHc = tensorHc + (tensorHc.*filterLower)';
        tensorHc(filterLower) = 0;
        
    case 3
        % mapping:    nC x nC
        % tensorFull: nS x nS x nS
        % tensor:     nC x nC x nC
        % tensorHc:   nC x nC x nC
        tensor = tensorFull(indexVector, indexVector, indexVector);
        nCol = size(mapping, 2);
        nRow = size(mapping, 1);
        
        % meas = sum(i,j,k): mixV(i,j,k) * lam(i) * lam(j) * lam(k)
        % meas = sum(i,j,k): mixV(i,j,k) * (sum(q): mapping(i, q) * hc(q))
        %                                * (sum(r): mapping(j, r) * hc(r))
        %                                * (sum(s): mapping(k, s) * hc(s))
        % meas = sum(q,r,s): hc(q) * hc(r) * hc(s)
        %      * sum(k) mapping(k,s) * (
        %          sum(i,j): mapping(i, q) * mixV(i,j,k) *  mapping(j, s)
        % meas = sum(q,r,s): hc(q) * hc(r) * hc(s) * tensor(q,r,s)
        % 
        % where the inner most sum can be calculated as above, and is
        %   summed up in a for loop
       
        tensorHc = zeros([nCol, nCol, nCol]);
        for iResult = 1:nCol
            
            tensorHcTmp = zeros([nCol, nCol, nCol]);
            
            for iSum = 1:nRow
                tensorHcTmp(:,:,iSum) = mapping(iSum, iResult) * ...
                              (mapping' * tensor(:,:,iSum) * mapping);
            end
            
            tensorHc(:,:,iResult) = sum(tensorHcTmp, 3);
        end
        
        % to make the hC-tensors comparable, we must see that the
        %   effective coefficient for (i=j=k) is:
        %      V(i,j,k)
        %   and for (i = j, j < k) is:
        %      V(i,j,k) + V(i,k,j)
        %   and for (i < j, j = k) is:
        %      V(i,j,k) + V(j,i,k)
        %   and for (i < j < k) is:
        %      V(i,j,k) + V(i,k,j) + V(j,i,k) + V(j,k,i) + V(k,i,j) +
        %      V(k,j,i)
        
        iTensor = ones([nCol, nCol, nCol]);
        jTensor = iTensor;
        kTensor = iTensor;
        for i = 1:nCol
            iTensor(i,:,:) = iTensor(i,:,:) * i;
            jTensor(:,i,:) = jTensor(:,i,:) * i;
            kTensor(:,:,i) = kTensor(:,:,i) * i;
        end
        
        % We assign the values that should not be zero!
        %   i = j = k (nothing to do)
        %   i < j < k
        filter = (iTensor < jTensor) & (jTensor < kTensor);
        tensorHc = tensorHc ...
            + permute(tensorHc.*permute(filter, [1 3 2]), [1 3 2]) ...
            + permute(tensorHc.*permute(filter, [2 1 3]), [2 1 3]) ...
            + permute(tensorHc.*permute(filter, [3 1 2]), [2 3 1]) ...
            + permute(tensorHc.*permute(filter, [2 3 1]), [3 1 2]) ...
            + permute(tensorHc.*permute(filter, [3 2 1]), [3 2 1]);
                              
        %   i = j, j < k
        filter = (iTensor == jTensor) & (jTensor < kTensor);
        tensorHc = tensorHc + ...
            + permute(tensorHc.*permute(filter, [1 3 2]), [1 3 2]) ...
            + permute(tensorHc.*permute(filter, [3 1 2]), [2 3 1]);
        
        % j = k, i < j
        filter = (iTensor < jTensor) & (jTensor == kTensor);
        tensorHc = tensorHc + ...
            + permute(tensorHc.*permute(filter, [2 3 1]), [3 1 2]) ...
            + permute(tensorHc.*permute(filter, [2 1 3]), [2 1 3]);
        
        % And we delete the values that are not permitted
        % i > j or j > k
        tensorHc((iTensor > jTensor) | (jTensor > kTensor)) = 0;
        
end

end
