Back to the main page.
Bug 3058 - error using cross-correlations ft_spike_xcorr
| Status | CLOSED FIXED |
| Reported | 2016-02-01 11:45:00 +0100 |
| Modified | 2016-05-05 20:31:15 +0200 |
| Product: | FieldTrip |
| Component: | spike |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Windows |
| Importance: | P5 normal |
| Assigned to: | Jan-Mathijs Schoffelen |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Danial Arabali - 2016-02-01 11:45:53 +0100
Hello, I recently started to learn Fieldtrip toolbox functions (v. 20151228) to use for my spike-LPF analysis. Thank you for providing the toolbox and your supports! While trying out the spike tutorial provided by the http://www.fieldtriptoolbox.org/tutorial/spike, I had a problem running part of the code related to cross-correlations: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% spike = ft_read_spike('p029_sort_final_01.nex'); cfg = []; cfg.spikechannel = {'sig001U_wf', 'sig002U_wf', 'sig003U_wf', 'sig004U_wf'}; % select only the MUA spike = ft_spike_select(cfg, spike); cfg = []; cfg.dataset = 'p029_sort_final_01.nex'; cfg.trialfun = 'trialfun_stimon'; cfg = ft_definetrial(cfg); cfg.timestampspersecond = spike.hdr.FileHeader.Frequency; % 40000 spikeTrials = ft_spike_maketrials(cfg,spike); cfg = []; cfg.maxlag = 0.2; % maximum 200 ms cfg.binsize = 0.001; % bins of 1 ms cfg.outputunit = 'proportion'; % make unit area cfg.latency = [-2.5 0]; cfg.vartriallen = 'no'; % do not allow variable trial lengths cfg.method = 'xcorr'; % compute the normal cross-correlogram Xc = ft_spike_xcorr(cfg,spikeTrials); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% I got the following error: Subscript indices must either be real positive integers or logicals. Error in ft_spike_xcorr (line 350) stat.label = spike.label(chansel); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The "spikeTrials.label" is a 1*4 cellarray, but the function "ft_channelcombination" inside the "ft_spike_xcorr" requires "cfg.channelcmb" as a two-column array. So there is a discrepancy which I don't understand how we need to provide the input for the "ft_spike_xcorr" or more specifically "spikeTrials.label"! Or maybe there is something else that I don't understand because I am naive to the toolbox. I would appreciate to help me figuring out the issue. Thank you.
valentina - 2016-02-22 12:39:34 +0100
(In reply to Danial Arabali from comment #0) Hi Danial, I'm stuck exactly at your point with the same error. Today I tried to do some detective work because I think the first error is at the line 138: chansel = unique(cmbindx(:)); % get the unique channels because it takes also the value 0 as index for channels. So i corrected like: chansel = unique(cmbindx(:)); % get the unique channels chansel = chansel(chansel~=0); % revome zeros and the script works fine to the end, but the output is different: infact I don't get the instruction at line 133 for k=1:size(cfg.channelcmb,1) cmbindx(k,1) = strmatch(cfg.channelcmb(k,1), spike.label, 'exact'); cmbindx(k,2) = strmatch(cfg.channelcmb(k,2), spike.label, 'exact'); end Why does it is set to search for only two channels if we specified 4 channels? (cfg.spikechannel = {'sig001U_wf', 'sig002U_wf', 'sig003U_wf', 'sig004U_wf'}; % select only the MUA) Maybe the instruction at this point should be something like this for k=1:size(cfg.channelcmb,1) for c=1:length(spike.label) cmbindx(k,c) = strmatch(cfg.channelcmb(k,c), spike.label, 'exact'); end end Anyway, I tried this solution and the output is of the right size now (xcorr: [4x4x400 double] instead of xcorr: [2x2x400 double]), but I don't think the content is right. The last test I did was to select only the channels 3 and 4 cfg.spikechannel = { 'sig003U_wf', 'sig004U_wf'}; % select only the MUA and run the original script with only the first modified instruction chansel = unique(cmbindx(:)); % get the unique channels chansel = chansel(chansel~=0); % revome zeros Then I followed the second part of the tutorial % compute the shuffled correlogram cfg.method = 'shiftpredictor'; % compute the shift predictor Xshuff = ft_spike_xcorr(cfg,spikeTrials); iCmb = 3; jCmb = 4; figure xcSmoothed = conv(squeeze(Xc.xcorr(iCmb,jCmb,:)),ones(1,5)./5,'same'); % do some smoothing hd = plot(Xc.time(3:end-2),xcSmoothed(3:end-2),'k'); % leave out borders (because of smoothing) hold on xcSmoothed = conv(squeeze(Xshuff.shiftpredictor(iCmb,jCmb,:)),ones(1,5)./5,'same'); plot(Xc.time(3:end-2),xcSmoothed(3:end-2),'r') hold on xlabel('delay') ylabel('proportion of coincidences') title([Xc.label{iCmb} Xc.label{jCmb}]) axis tight And I just changed iCmb = 3; jCmb = 4; in iCmb = 1; jCmb = 2; And I got exactly the same figure output of the tutorial. I hope this can be a little help for the work Valentina
Danial Arabali - 2016-02-22 14:30:08 +0100
(In reply to valentina from comment #1) Hallo and thanks Valentina, that's right. I don't understand as well how the following loop should look like, but i seems that it is not correct like this! line 133 @ ft_spike_xcorr for k=1:size(cfg.channelcmb,1) cmbindx(k,1) = strmatch(cfg.channelcmb(k,1), spike.label, 'exact'); cmbindx(k,2) = strmatch(cfg.channelcmb(k,2), spike.label, 'exact'); end I would prefer your approach with two for loops, But I stuck in a few lines back and more specifically, what should we expect as an output of the ft_channelcombination cfg.channelcmb = ft_channelcombination(cfg.channelcmb, spike.label,true); it gives a vector of all channel combinations?! (of size(1,20)) e.g. for the currect example in tutorial, if I say ft_channelcombination(cfg.channelcmb, spike.label(1:2),true) I will get: {'sig001U_wf' 'sig002U_wf' 'sig002U_wf' 'sig001U_wf' 'sig001U_wf' 'sig002U_wf'} which I don't know how it should be read?! By the way, with correcting the code according to your suggestions, I got the Xc and Xshuff without erros but if I plot the results with iCmb = 3; jCmb = 4; I got a blank figure due to the NAN vector of xcSmoothed and with iCmb = 1; jCmb = 2; I got a different plot as in the tutorial with different pick value! Cheers, Danial
Danial Arabali - 2016-02-22 14:40:42 +0100
(In reply to Danial Arabali from comment #2) (In reply to valentina from comment #1) So, if I specify only two channels: cfg.spikechannel = {'sig003U_wf', 'sig004U_wf'} and let the ft_spike_xcorr with both corrections: 1) two for loops to get cmbindx(k,c); 2) chansel = chansel(chansel~=0); the code runs, but for plotting and getting the same figure as in the tutorial I specified: iCmb = 2; jCmb = 1; The pick of the graph is at (0.0005, 0.003195)
valentina - 2016-02-22 16:04:28 +0100
Created attachment 776 cross correlation 1 vs 2
valentina - 2016-02-22 17:13:54 +0100
(In reply to Danial Arabali from comment #3) (In reply to valentina from comment #4) Yes, exactly Danial. I got the same results of the tutorial as you described in comment #3 Instead, selecting all the 4 channels {'sig001U_wf', 'sig002U_wf', 'sig003U_wf', 'sig004U_wf'} and using both 1) two for loops to get cmbindx(k,c); 2) chansel = chansel(chansel~=0); I got the same NaN vector xcSmoothed as yours for iCmb = 3; jCmb = 4; and the figure that I attached in comment #4 for iCmb = 1; jCmb = 2;
valentina - 2016-02-22 17:26:30 +0100
(In reply to valentina from comment #5) (In reply to Danial Arabali from comment #3) (In reply to valentina from comment #4) THANK YOU, DANIAL! I solved it thanks your suggestion of looking at the output of cfg.channelcmb = ft_channelcombination(cfg.channelcmb, spike.label,true); As it's specified in the comment at the beginning of the script ft_channelcombination: "You should specify channel combinations as a two-column cell array" So I just changed the last line (170) of ft_channelcombination script from collect = [datachannel(indx(:,1)) datachannel(indx(:,2))]; to collect = [datachannel(indx(:,1))' datachannel(indx(:,2))']; to have an output with 2 columns and on each row all the possible combinations of selected channels. So, following the original tutorial, the output will be cfg.channelcmb: 'sig001U_wf' 'sig001U_wf' 'sig002U_wf' 'sig001U_wf' 'sig003U_wf' 'sig001U_wf' 'sig004U_wf' 'sig001U_wf' 'sig002U_wf' 'sig002U_wf' 'sig003U_wf' 'sig002U_wf' 'sig004U_wf' 'sig002U_wf' 'sig003U_wf' 'sig003U_wf' 'sig004U_wf' 'sig003U_wf' 'sig004U_wf' 'sig004U_wf' I deleted the changes I made beforen ft_spike_xcorr 1) two for loops to get cmbindx(k,c); 2) chansel = chansel(chansel~=0); because at this point they were useless! Summarizing, the only thing you have to change is the line 170 of ft_channelcombination script in collect = [datachannel(indx(:,1))' datachannel(indx(:,2))'];
Jan-Mathijs Schoffelen - 2016-02-22 17:32:55 +0100
Hi both, thanks for taking this up. Is it then fair to summarize the issue as a row-vector versus column vector problem? If so, would it also work if you do cfg.spikechannel = {'sig001U_wf', 'sig002U_wf', 'sig003U_wf', 'sig004U_wf'}'; instead of cfg.spikechannel = {'sig001U_wf', 'sig002U_wf', 'sig003U_wf', 'sig004U_wf'}; (and then not adjust ft_channelcombination?)
valentina - 2016-02-22 18:06:53 +0100
(In reply to Jan-Mathijs Schoffelen from comment #7) Yes, it's definitley a a row-vector versus column vector problem. Do you mean: cfg.spikechannel = {'sig001U_wf'; 'sig002U_wf'; 'sig003U_wf'; 'sig004U_wf'}'; instead of cfg.spikechannel = {'sig001U_wf', 'sig002U_wf', 'sig003U_wf', 'sig004U_wf'};? If it's so, it works fine (without modify ft_channelcombination) but it's not enough to change it when you select the channels at the beginning of the tutorial (at this point: % read in the data, select the channels and define the trials spike = ft_read_spike('p029_sort_final_01.nex'); cfg = []; cfg.spikechannel = {'sig001U_wf', 'sig002U_wf', 'sig003U_wf', 'sig004U_wf'}; % select only the MUA)) I had to redefine cfg.spikechannel = {'sig001U_wf'; 'sig002U_wf'; 'sig003U_wf'; 'sig004U_wf'}'; before performing Xc = ft_spike_xcorr(cfg,spikeTrials); i.e.: spikeTrials.label={'sig001U_wf'; 'sig002U_wf'; 'sig003U_wf'; 'sig004U_wf'}; cfg = []; cfg.maxlag = 0.2; % maximum 200 ms cfg.binsize = 0.001; % bins of 1 ms cfg.outputunit = 'proportion'; % make unit area cfg.latency = [-2.5 0]; cfg.vartriallen = 'no'; % do not allow variable trial lengths cfg.method = 'xcorr'; % compute the normal cross-correlogram Xc = ft_spike_xcorr(cfg,spikeTrials); Valentina
Danial Arabali - 2016-02-22 18:42:26 +0100
(In reply to valentina from comment #8) (In reply to Jan-Mathijs Schoffelen from comment #7) Hi both, thanks! It seems that as Valentina explained we just have to change the line 170 in ft_channelcombination. I also checked and it works fine. Thank you. Best
Jan-Mathijs Schoffelen - 2016-02-22 20:53:23 +0100
I still don't fully agree with the proposed change. The cause seems to be that spike.label is a row-vector. This is against convention, because data.label is typically a column-vector. If we add a transposition in ft_channelcombination it will fix it for this instance, but it will break everything else. I think ft_read_spike should be changed such that it returns a column-vector in the label field. Better still, is to ensure that the input into ft_channelcombination is a column-vector. Even better still, is not to rely in ft_channelcombination on the input to be a column vector, but impose this. The very last solution is the most generic, but would require a change at the beginning of the function. Do you agree with this reasoning?
valentina - 2016-02-23 10:53:41 +0100
(In reply to Jan-Mathijs Schoffelen from comment #10) Hi Jan-Mathijs, this was my very first approach to FieldTrip, so I don't know how changing ft_channelcombination could affect other processes, but I think it surely will. So I agree with you that the best way is to impose a column vector to the label changing ft_read_spike (at least in the part for nex and plx formats, I don't know the others filetype) at line 160 and 169 I just tried to change spike.label{chan} = deblank(hdr.VarHeader(i).Name); with spike.label{chan,1} = deblank(hdr.VarHeader(i).Name); and the same at line 204 spike.label{i,1} = deblank(hdr.ChannelHeader(i).Name); and for the spike cross correlation tutorial everything works fine without changing anything else. Do you think it's enough or should we change something else?
Jan-Mathijs Schoffelen - 2016-02-23 11:04:57 +0100
Hi Valentina, In principle I would say indeed that your proposed change is fine. Yet I read in the the help of ft_read_spike that the label-field is designed to be 1xN. Although this is not according to convention of the other types of data, I am not sure whether anywhere else in the code the label-field is expected/assumed to be a row-vector. My proposal would be to explicitly make datachannel a column vector in ft_channelcombination. This will probably be the most robust.
Jan-Mathijs Schoffelen - 2016-02-23 11:23:03 +0100
I have changed the code, and it should now work again. The current version of the code can be grabbed from github.com/fieldtrip/fieldtrip For more information about the changes, see: https://github.com/fieldtrip/fieldtrip/pull/103/files Thanks for looking into this!
valentina - 2016-02-23 11:45:13 +0100
(In reply to Jan-Mathijs Schoffelen from comment #13) Thanks Jan-Mathijs! I didn't notice in the help of ft_read spike that the label-field is designed to be 1xN, so I agree that your way is the most correct one since also in ft_channelselection, line 138, it's required that datachannel is a column vector % ensure that both inputs are column vectors datachannel = datachannel(:); Thank you again for your time Best, Valentina