demixed_SCA/dsca_cv.m at master · CCNHuntLab/demixed_SCA · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
%%% history
%%%   - 2020/10/22 y.takagi - initially created with modifying Dmtry Kobak's dPCA program
%%%     see also: https://github.com/machenslab/dPCA
function [meanError] = dsca_cv(Xfull, ...
    Xtrial, Yfull, Ytrial, numOfTrialsX, numOfTrialsY, varargin)


% default input parameters
options = struct('numComps',       25,                  ...
                 'lambda',         1e-10, ...
                 'numRep',         10,                  ...
                 'combinedParams', [],                  ...
                 'method',         'naive',          ...
                 'simultaneous',   false);

% read input parameters
optionNames = fieldnames(options);
if mod(length(varargin),2) == 1
	error('Please provide propertyName/propertyValue pairs')
end
for pair = reshape(varargin,2,[])    % pair is {propName; propValue}
	if any(strcmp(pair{1}, optionNames))
        options.(pair{1}) = pair{2};
    else
        error('%s is not a recognized parameter name', pair{1})
	end
end

if min(numOfTrialsX(:)) <= 0 || min(numOfTrialsY(:)) <= 0
    error('dSCA:tooFewTrials0','Some neurons seem to have no trials in some condition(s)')
elseif min(numOfTrialsX(:)) == 1 || min(numOfTrialsY(:)) == 1
    error('dSCA:tooFewTrials1','Cannot perform cross-validation')
end

Xsum = bsxfun(@times, Xfull, numOfTrialsX);
Ysum = bsxfun(@times, Yfull, numOfTrialsY);

for rep = 1:options.numRep
    if options.simultaneous
        [Xtest, Ytest, ~, ~] = dsca_getTestTrials_noave(Xtrial, Ytrial, min(numOfTrialsX), ...
            'simultaneous', options.simultaneous);
    else
        [Xtest, ~] = dpca_getTestTrials(Xtrial, numOfTrialsX, ...
            'simultaneous', options.simultaneous);
        [Ytest, ~] = dpca_getTestTrials(Ytrial, numOfTrialsY, ...
            'simultaneous', options.simultaneous);
    end
    Xtrain = bsxfun(@times, Xsum - Xtest, 1./(numOfTrialsX-1));
    Ytrain = bsxfun(@times, Ysum - Ytest, 1./(numOfTrialsY-1));

    XtestCen = bsxfun(@minus, Xtest, mean(Xtest(:,:),2));
    YtestCen = bsxfun(@minus, Ytest, mean(Ytest(:,:),2));
    YtestMargs = dpca_marginalize(YtestCen, 'combinedParams', options.combinedParams, ...
                    'ifFlat', 'yes');
    for i=1:length(YtestMargs)
        margTestVar(i) = sum(YtestMargs{i}(:).^2);
    end

    if strcmp(options.method, 'naive')
        margVar_toNormalize = margTestVar;
    end


    [W,V,whichMarg] = dsca(Xtrain, Ytrain, options.numComps, ...
        'combinedParams', options.combinedParams, ...
        'lambda',options.lambda);

    for i=1:length(YtestMargs)
        recError = 0;

        if strcmp(options.method, 'naive')
            recError = sum(sum((YtestMargs{i} - V(:,whichMarg==i)*W(:,whichMarg==i)'*XtestCen(:,:)).^2));
        end

        errorsMarg(i, rep) = recError/margVar_toNormalize(i);
    end
end

meanError = mean(errorsMarg,2);