computeConformedMultiInstanceKernel.m

function [conformedMultiInstanceKernel, rawConformedMultiInstanceKernel, instancesTransformationKernel, clusterCentres] = computeConformedMultiInstanceKernel(instanceStarts, instanceEnds, instanceWideKernel, trainIndices, allSeqsAsBags, nClusters, sig, Y, computationVersion, debugLevel, debugMsgLocation)

% COMPUTECONFORMEDMULTIINSTANCEKERNEL
% Computes the conformally transformed multi-instance kernel

% INPUT PARAMS
% Param 'instanceStarts' (vector)
% Indices where instances (segments) of any sequence begin
%
% Param 'instanceEnds' (vector)
% Indices where instances (segments) of any sequence end
% 
% Param 'instanceWideKernel' (matrix)
% The instance-wide kernel matrix of dimension #Instances x #Instances
%
% Param 'trainIndices' (vector)
% Indices of the training sequences
%
% Param 'allSeqsAsBags' (cell array)
% Collection of all segments per sequence
% 
% Param 'nClusters'
% number of clusters for k-means
%
% Param 'sig' 
% Gaussian bandwidth (sigma) value for the Gaussian transformation
%
% Param 'Y' (vector)
% Vector of labels
%
% Param 'computationVersion' ('AccumArray'/'Looping')
% Flag specifying what approach to use for computing the kernel
%
% Param 'debugLevel' (0/1/2)
%
% Param 'debugMsgLocation' (1/fileID)
%
%
% OUTPUT PARAMS
% Param 'conformedMultiInstanceKernel' (matrix)
% Normalized conformally transformed MI kernel
% 
% Param 'rawConformedMultiInstanceKernel' (matrix)
% Unnormalized conformally transformed MI kernel
%
% Param 'instancesTransformationKernel' (matrix)
% Transformation of the training instances
%
% Param 'clusterCentres' (matrix)
% Cluster centres obtained upon k-means
%
% ADDITIONAL NOTES
% -- Stratification
%    . If stratification is required, see below, set stratified = 1
% -- Replicates for k-means
%    . currently set to 1
%    . can be increased
%
% Author: snikumbh@mpi-inf.mpg.de

 
% 1.1 Obtain expansion points
%
% -- Convert allSeqsAsBags in to a matrix X of dimensions nInstances x FVlength 
% 
% -- spread out the instances, get the nInstances
% -- trainIndices are trainingIndices
% Edit: We earlier expected trainIndices to be a row vector. Now it can be either 
% a row or a column vector.

if nargin < 10
    % default
    computationVersion = 'Looping';
end
nBags = length(allSeqsAsBags);
% nBags, only as many bags as training examples
idxInstances = zeros(1,nBags);
stratified = 0;
% The following two lines also used for stratified representation.
if stratified
    posInd = size(find(Y > 0),2);
    negInd = size(find(Y < 0),2);
end

if strcmp(computationVersion, 'AccumArray')
    instanceIDVector = zeros(1,instanceEnds(end));
end


% Maintain a count of the instances arising from training indices/bags
% This is used later to define the number of rows for XtoPassFrom
currentTotal = 0; 
for i=1:nBags
    % #instances at each index/bag
    idxInstances(i) = size(allSeqsAsBags{i}, 2);
    if strcmp(computationVersion, 'AccumArray')
        instanceIDVector(instanceStarts(i):instanceEnds(i)) = repmat(i, [1 size(allSeqsAsBags{i}, 2)]);
    end
    % The following is useful to have stratified repressentation, 
    % we don't need this now.
    if stratified
        if i <= posInd
            nInstancesPos = nInstancesPos + size(allSeqsAsBags{i}, 2);
        end 
    end
    if any(i==trainIndices) == 1 %check, if a trainingIndex, get instances
        currentTotal = currentTotal+idxInstances(i);
    end 
    % useful to count instances only in training bags
end
forIndicesInX = cumsum(idxInstances);
nInstances = sum(idxInstances);

if strcmp(computationVersion, 'AccumArray')
    assert(size(instanceIDVector,2) == nInstances);
end

logMessages(debugMsgLocation, sprintf('\n--- currentTotal for setting a zeros XtoPassFrom: %d\n', currentTotal), debugLevel);
XtoPassFrom = zeros( currentTotal, size(allSeqsAsBags{1},1) );
% Make sure that only instances from training bags are taken
% trainIndices store the right set of indices as passed to it
currentTotal = 0;
for i=1:length(trainIndices)%indices of relevant bags/ training bags
    XtoPassFrom( currentTotal+1 : currentTotal + idxInstances(trainIndices(i)) , : ) = transpose(allSeqsAsBags{ trainIndices(i) });
    currentTotal = currentTotal + idxInstances( trainIndices(i) );
end
logMessages(debugMsgLocation, sprintf('--- allSeqsAsBags is now n x p matrix, ready for kmeans\n'), debugLevel);
%
% 1.2. Apply kmeans 
%      -- we are using Matlab's algorithm for kmeans.
%      -- using the buckshot hueristic means:
%         - randomly sample sqrt(nClusters * nInstances) 
%           data points from nInstances
%         - this will give nClusters

% If sending complete XtoPassFrom for k-means clustering set XtoPass = XtoPassFrom;
% Replicates
nofRep = 1;
% Randomly sample without replacement
% XtoPassFrom contains only the relevant instances, hence freely select any instances from it
nSampled = round(sqrt(nClusters * currentTotal));
logMessages(debugMsgLocation, sprintf('--- for kmeans, Total: %d, nSampled: %d, nClusters: %d\n', currentTotal, nSampled, nClusters), debugLevel);
XtoPass = XtoPassFrom( transpose(randsample( currentTotal, round(sqrt(nClusters * currentTotal)) )), :);
%% Without the buckshot heuristic, use
% XtoPass = XtoPassFram
clear('XtoPassFrom');
[clusterCentres, matrixOfDistancesFromCentresForPassedVectors] = getExpansionPoints(XtoPass, nClusters, nofRep, debugLevel, debugMsgLocation);
clear('XtoPass');
% clusterCentres is a k-by-p matrix
% matrixOfDistancesFromCentresForPassedVectors is n-by-k matrix

%
% 2. As many kernels as expansion points
%

% Use the fact that only comparisons with the clusterCentre 
% are to be computed anew. Other values can be fetched from the instance-wide kernel.
% transformationKernel stores all products of the instances with the cluster centers
% transformationKernel is Gaussian RBF
tic;
logMessages(debugMsgLocation, sprintf('--- applying transformation...'), debugLevel);
[instancesTransformationKernel] = applyTransformation(allSeqsAsBags, clusterCentres, sig, debugLevel, debugMsgLocation);
logMessages(debugMsgLocation, sprintf('done in %.3f  seconds\n', toc), debugLevel);
% instancesTransformationKernel is cell array for number of bags
% instancesTransformationKernel{eachBag} is nInstanceInBag-by-nClusters

% Initialize kernels
for i=1:nClusters
    conformedMultiInstanceKernel{i} = zeros(nBags);
    rawConformedMultiInstanceKernel{i} = zeros(nBags);
end

logMessages(debugMsgLocation, sprintf('--- computing the conformedKernel...'), debugLevel);
if strcmp(computationVersion, 'AccumArray')
    logMessages(debugMsgLocation, sprintf('\nAccumarray version---\n'), debugLevel);tic;
    instanceSubs = combvec(instanceIDVector, instanceIDVector)';
    % instanceSubs = tempInstanceSubs(tempInstanceSubs(:,1) <= tempInstanceSubs(:,2),:);
    % instanceSubs are only the upper-triangular indices including the diagonal
    % clear('tempInstanceSubs');
    instanceWideKernelRepeated = repmat(instanceWideKernel, 1, 1, nClusters);
    allInstancesTransformations = cat(1,instancesTransformationKernel{:});
    for i=1:nClusters
        KernelAsVector = (allInstancesTransformations(:,i) * allInstancesTransformations(:,i)') .* instanceWideKernelRepeated(:,:,i);
        tempKernel = accumarray(instanceSubs, KernelAsVector(:));
        assert(size(tempKernel,1) == size(zeros(nBags),1));
        % tempKernel = triu(tempKernel, 1) + tempKernel'; %because instanceSubs stored indices for only the upper-triangular part of the matrix
        conformedMultiInstanceKernel{i} = normalizeKernel(tempKernel);
    end
    clear('instanceWideKernelRepeated');
    clear('instanceIDVector');
    clear('instanceSubs');
    clear('allInstancesTransformations');
    logMessages(debugMsgLocation, sprintf('done in %.3f  seconds', toc), debugLevel);
end
if strcmp(computationVersion, 'Looping')
    % One can convert the transformationKernels into 3D matrices
    logMessages(debugMsgLocation, sprintf('\n--- Looping-over-bags version...'), debugLevel);
    % Performing the reshaped version
    if issparse(allSeqsAsBags{1})
	logMessages(debugMsgLocation, sprintf('Doing the reshaped version..\n'), debugLevel);
        % Pre-allocate memory, doesn't help in speed
        % for i=1:nClusters
        %     rawConformedMultiInstanceKernel{i} = zeros(nBags);
        %     conformedMultiInstanceKernel{i} = zeros(nBags);
        % end
        tic;
        tempKernelCollection = zeros(nBags, nBags, nClusters);
        for b1=1:nBags
            b1r = reshape(instancesTransformationKernel{b1}, idxInstances(b1) , [], nClusters);
            for b2=b1:nBags
                b2r = reshape(instancesTransformationKernel{b2}, idxInstances(b2) , [], nClusters);
                b2rt = permute(b2r,[2,1,3]);
                b1b2_transformationProducts = bsxfun(@times, b1r , b2rt);
                tempToBeSummed = bsxfun(@times, b1b2_transformationProducts, instanceWideKernel(instanceStarts(b1):instanceEnds(b1), instanceStarts(b2):instanceEnds(b2)) );
		term = 1/(idxInstances(b1) * idxInstances(b2));
                tempKernelCollection(b1, b2, :) = sum(sum(tempToBeSummed,1));
            end
        end
        logMessages(debugMsgLocation, sprintf('Bag-wise kernel done in %.3f  seconds\n', toc), debugLevel);
        tic;
        for i=1:nClusters
            tempKernel = zeros(nBags);
            tempKernel = triu(tempKernelCollection(:,:,i), 1) + tempKernelCollection(:,:,i)';
            rawConformedMultiInstanceKernel{i} = tempKernel;
            conformedMultiInstanceKernel{i} = normalizeKernel(tempKernel);
        end
        logMessages(debugMsgLocation, sprintf('done in %.3f  seconds\n', toc), debugLevel);
    else % when not sparse
	% not tested either for correctness or efficiency. We recommend only using sparse vectors, thus the approach above.
        tic;
        for i=1:size(clusterCentres,1)
            tempKernel = zeros(nBags);
            for b1=1:nBags
    	    tempb1 = instancesTransformationKernel{b1}(:,i);
    	        for b2=b1:nBags
    	            tempKernel(b1,b2) = sum(sum((tempb1 * instancesTransformationKernel{b2}(:, i)') .* instanceWideKernel(instanceStarts(b1):instanceEnds(b1), instanceStarts(b2):instanceEnds(b2)) ));
    	        end
            end
            tempKernel = triu(tempKernel, 1) + tempKernel';
	    rawConformedMultiInstanceKernel{i} = tempKernel;
            conformedMultiInstanceKernel{i} = normalizeKernel(tempKernel);
        end
        logMessages(debugMsgLocation, sprintf('done in %.3f  seconds\n', toc));
    end
end
%%%%%

% check psd
for i=1:size(clusterCentres,1)
    [~,p] = chol(conformedMultiInstanceKernel{i});
    if p==0
    else
        logMessages(debugMsgLocation,sprintf('cluster %d Cholesky p: %.4f\n', i, p), debugLevel);
        % p~=0 means there is an eigen value < 0
	% print the min(eigen values) in the log
        logMessages(debugMsgLocation, sprintf('%e,', min(eig(conformedMultiInstanceKernel{i}))), debugLevel);
        logMessages(debugMsgLocation, sprintf('Check kernel w.r.t. clusterCentre %d out of %d, may be not psd!\n', i, size(clusterCentres, 1)), debugLevel);
    end
end

end % function ends