1 function [c,U,V,c_perc,U_perc,V_perc] = tkcca(X,Y,tau,kappas)
   2 % [c, U, V, c_bstrp, U_bstrp, V_bstrp] = tkcca(X,Y,tau,kappas)
   3 %
   4 % 	temporal kernel CCA
   5 %	
   6 % INPUT
   7 %	X	a nDimensions-by-Samples data matrix
   8 %	Y	a mDimensions-by-Samples data matrix
   9 %	tau	the maximal time lag by which X is shifted with respect to Y
  10 %	kappas	the regularizers
  11 %	
  12 % OUTPUT
  13 %	c	the canonical correlogram, from -tau to tau
  14 %	U	the time resolved canonical variate for X
  15 %	V	the canonical variate for Y
  16 %	
  17 %	USAGE 
  18 %
  19 %>> [c,U,V] = tkcca(X,Y,10,.1)
  20 %		
  21 %	computes canonical correlogram c, a time resolved variate U and canonical 
  22 %	projection V for timelag -10 to 10 with fixed regularizer of .1
  23 %
  24 
  25 if nargin<4, kappas = [10.^-[0:7]]'; end
  26 
  27 % center data
  28 X = X - repmat(mean(X,2),1,size(X,2));
  29 Y = Y - repmat(mean(Y,2),1,size(Y,2));
  30 % embed one signal in time shifted copies of itself
  31 [X, timeidx, tauidx] = embed(X,tau);
  32 % compute the linear kernels
  33 kY = Y(:,timeidx)' * Y(:,timeidx);
  34 kX = X' * X;
  35 % find the right regularizer
  36 kappaOpt	= optimize_kappa(kX,kY,kappas);
  37 % compute kcca using the right regularizer
  38 [r, a, b]	= kcca(kX,kY,kappaOpt);
  39 % reconstruct the canonical variates and compute canonical correlogram
  40 [U, V, c]	= reconstruct(X,Y(:,timeidx),a,b,tauidx);
  41 
  42 
  43 function [eX, timeidx, tauidx] = embed(X,tau)
  44 	% embed the first signal in its temporal context
  45 	[D T] = size(X);
  46 	% in case tau is a scalar, make it a vector from -tau to tau
  47 	if length(tau)==1,tau = -tau:tau;end
  48 	startInd 	= abs(tau(1)) + 1;
  49 	stopInd		= T - abs(tau(end));
  50 	len			= stopInd - startInd + 1;
  51 	% create a column vector that contains the indices of the first segment
  52 	idx = repmat((startInd:stopInd)', 1, length(tau)) + repmat(tau, len, 1); 	
  53 	% create (linear) indices for the different dimensions
  54 	dim_offset = repmat( (0:D-1)*T, length(tau)*len, 1);
  55 	idx = repmat(idx(:), 1, D) + dim_offset;
  56 	% for the linear indices we need column-signals
  57 	X = X';
  58 	% get the data (D channels, segments are concatenated) and reshape it
  59 	eX = reshape(X(idx), len, length(tau)*D)';
  60 	tauidx = repmat(tau',D,1); 
  61 	timeidx = startInd:stopInd;
  62 
  63 function [r,a,b] = kcca(kX,kY,kappa)
  64 	% compute the dual coefficients 
  65 	n = size(kX,1);
  66 	options.disp = 0;
  67 	% force kernel symmetry
  68 	kX = (kX+kX')/2; kY = (kY+kY')/2;
  69 	% normalise the spectral norm of the matrices
  70 	kX = kX./max(eig(kX));
  71 	kY = kY./max(eig(kY));
  72 	% Generate LH
  73 	LH = [zeros(n) kX*kY';kY*kX' zeros(n)];
  74 	% generate RH with regularization ridge
  75 	RH = [kX*kX' zeros(n);zeros(n) kY*kY'] + eye(2*n)*kappa;
  76 	% make sure the matrices are symmetric
  77 	RH=(RH+RH')/2; LH=(LH+LH')/2;
  78 	% Compute the generalized eigenvectors
  79 	[Vs,r]=eigs(LH,RH,1,'LA',options);
  80 	a = Vs(1:n);
  81 	b = Vs(n+1:end);
  82 	
  83 function kappa = optimize_kappa(kX,kY,kappas,iterations)
  84 	if nargin<4, iterations=10; end
  85 	shcors = zeros(length(kappas),iterations);
  86 	skX = kX; skY = kY;
  87 	% try each regularizer
  88 	for iR = 1:length(kappas)
  89 		r(iR,1) = kcca(kX,kY,kappas(iR));
  90 		% for all iterations of the reshuffling procedure
  91 		for iS = 1:iterations
  92 			idx = randperm(size(kX,1));
  93 			skX = kX(idx,idx);		
  94 		    % do cca on shuffled data
  95      		shcors(iR,iS) = kcca(skX,skY,kappas(iR));
  96 		end
  97 	end
  98   	% pick that regularizer that maximizes the distance between 
  99   	% true correlations and shuffled data correlations
 100   	[val,pick]  = max(mean((repmat(r,1,(iterations))-shcors).^2,2));
 101   	fprintf('Picked kappa=[ %2.8f ]\ncorrelations  %0.2f (true) vs. %0.2f (shuffled)\n',...
 102   			kappas(pick),r(pick),mean(shcors(pick,:)))
 103   	kappa = kappas(pick);
 104 
 105 function [U,V,c] = reconstruct(eX,Y,a,b,tidx)
 106 	% the number of time lags
 107 	nTau = length(unique(tidx));
 108 	% the number of dimensions
 109 	D = size(eX,1)/nTau;
 110 	% the time-lag sorted indices
 111 	[sorted, sortInds] = sort(tidx);
 112 	% the zero lag canonical component
 113 	pY = ((Y * b)' * Y)';
 114 	% the time shifted canonical components
 115 	pX = repmat(eX * a, 1, size(eX, 2) ) .* eX;
 116 	if D>1
 117 		pX = squeeze(mean(reshape(pX(sortInds,:), [D,nTau, size(eX,2)] )));
 118 	end
 119 	% the correlations between the zero lag component of Y and the 
 120 	% time shifted components of X (i.e. the canonical correlogram)
 121 	c = corr(pY,pX');
 122 	% the time resolved variate
 123 	U = reshape(eX * a, nTau, D)';
 124 	% the other variate
 125 	V = Y * b;