SpeechSpikeAnalysis/PoissonSurpriseBruteTrial.m at master · Brain-Modulation-Lab/SpeechSpikeAnalysis · 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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

%%%% combination method
figdir = '/Users/brainmodulationlab/Dropbox (Brain Modulation Lab)/Speech Project/MerData/figures_DBS2000';


for idx=1:length(DATA)
    Stim = DATA(idx).Trial;
    clear BB PP BBbase PPbase;
    N0=4;
    fs = 1000;

    DDbase = arrayfun(@(x) DATA(idx).D((x-fs+1):x), DATA(idx).SpeechOnset.base_samples(DATA(idx).trange), ...
        'uniformoutput', false);

    %r = fs*nnz(cell2mat(DDbase))/length(cell2mat(DDbase));

    for trial=1:length(Stim.DD)
        %r = Stim.basemeanIFR;
        r = fs*nnz(Stim.DD{trial})/length(Stim.DD{trial});

        clear B P

        %% check for long trials and truncate to SpeechEnd + 1s
        D = Stim.DD{trial};

        if length(D)>(Stim.idxEventTimes{trial}(end)+fs)
            fprintf('Trial segment shortered from %4.2f secs to %4.2f secs.\n', length(D)/fs, ...
                (Stim.idxEventTimes{trial}(end)+fs)/fs)
            D = D(1:(Stim.idxEventTimes{trial}(end)+fs));
        end

        D = D(Stim.idxEventTimes{trial}(2):end);  %% Start detecting at command stim (Cue)

        if nnz(DDbase{trial}) > N0;

            fprintf('Unit %d, trial %d baseline: %d spikes...', idx, trial, nnz(DDbase{trial}));
            [Sb, Sp, comb, ~] = PoissonBrute(DDbase{trial}, fs, r, N0);
            fprintf(' done.\n');

            [Smax, Sind] = max(Sb);
            Sind = Sind(1);
            ii = comb{Sind};
            %B.ii = setdiff(ii, ii(end))';
            B.ii = ii';
            B.S = double(Smax);

            [Smax, Sind] = max(Sp);
            Sind = Sind(1);
            ii = comb{Sind};
            %P.ii = setdiff(ii, ii(end))';
            P.ii = ii';
            P.S = double(Smax);

        else
            B.ii = [];
            B.S = 0;
            P.ii = [];
            P.S = 0;
        end
        BBbase(trial) = B;
        PPbase(trial) = P;

        if nnz(D) > N0;

            fprintf('Unit %d, trial %d: %d spikes...', idx, trial, nnz(D));
            [Sb, Sp, comb, t] = PoissonBrute(D, fs, r, N0);
            fprintf(' done.\n');

            [Smax, Sind] = max(Sb);
            Sind = Sind(1);
            ii = comb{Sind};
            %B.ii = setdiff(ii, ii(end))';
            B.ii = ii';
            B.t = t(B.ii);
            [~,tmax] = max(DATA(idx).IFR(Stim.start_samples(DATA(idx).trange(trial)) + Stim.idxEventTimes{trial}(2) +...
            (round(fs*B.t(1)):round(fs*B.t(end)))));
            B.tmax = B.t(1) + tmax/fs;
            B.S = double(Smax);

            [Smax, Sind] = max(Sp);
            Sind = Sind(1);
            ii = comb{Sind};
            %P.ii = setdiff(ii, ii(end))';
            P.ii = ii';
            P.t = t(P.ii);
            m = max(DATA(idx).ISITS(Stim.start_samples(DATA(idx).trange(trial)) + Stim.idxEventTimes{trial}(2) +...
                (round(fs*P.t(1)):round(fs*P.t(end)))));
            Ptmax = bwconncomp(DATA(idx).ISITS(Stim.start_samples(DATA(idx).trange(trial)) + Stim.idxEventTimes{trial}(2) +...
                (round(fs*P.t(1)):round(fs*P.t(end))))==m);
            [~,Pidx] = max(cellfun(@length, Ptmax.PixelIdxList));
            tmax = floor(mean(Ptmax.PixelIdxList{Pidx}));
            P.tmax = P.t(1) + tmax/fs;
            P.S = double(Smax);

        else
            B.ii = [];
            B.t = [];
            B.tmax = [];
            B.S = 0;
            P.ii = [];
            P.t = [];
            P.tmax = [];
            P.S = 0;
        end
        BB(trial) = B;
        PP(trial) = P;
    end

%     figtitle = sprintf('Cue2: Rec %d %s, %5.3f, %s, Unit x, (%s)', ...
%         idx, DATA(idx).SubjectID, DATA(idx).depth, DATA(idx).elec, DATA(idx).grade);
%
%     RT = DATA(idx).RT(DATA(idx).trange);
%     [~,order] = sort(RT);
%
%     plot_rasterSpeechTaskBurstsOrder( [figdir,'/TrialBaselinedBurstMarked'], figtitle, ...
%         DATA(idx).Cue2, DATA(idx).trange, BB, quantile([BBbase(:).S],0.75), [], order, RT, 'b');
%
%     figtitle = sprintf('Cue2: Rec %d %s, %5.3f, %s, Unit x, (%s)', ...
%         idx, DATA(idx).SubjectID, DATA(idx).depth, DATA(idx).elec, DATA(idx).grade);
%
%     RT = DATA(idx).RT(DATA(idx).trange);
%     [~,order] = sort(RT);
%
%     plot_rasterSpeechTaskBurstsOrder( [figdir,'/TrialBaselinedPauseMarked'], figtitle, ...
%         DATA(idx).Cue2, DATA(idx).trange, PP, quantile([PPbase(:).S],0.75), [], order, RT, 'b');

    DATA(idx).Trial.BBbase = BBbase;
    DATA(idx).Trial.PPbase = PPbase;
    DATA(idx).Trial.BB = BB;
    DATA(idx).Trial.PP = PP;
end