gin commit from chris-mh16

Deleted files: 1

CODE/figfunctions/figure_neurotrace.asv

@@ -1,479 +0,0 @@
-function figure_neurotrace(fld) 
-
-fprintf('\n======================================================\n');
-fprintf('-- Creating Figure neurotrace --\n');
-fprintf('======================================================\n');
-
-%% settings
-green = [23, 105, 13]./255;
-red = [201, 0, 34]./255;
-cols = [green; 0.3 0.3 0.3; red; 1.00,0.84,0.00];
-
-snrthres = 2.5;
-smoothfact = 15; mindays = 3;
-falpha = 1; % set to 1 for converting to illustrator
-
-do_cleantraces = true; % if true we throw out artefact traces
-% citerion: avg emplitude in first 100ms after stimulus onset should be 
-% at least 3 times the standard deviation in the preceding 100 ms
-do_latency = true;
-
-%% load RTs
-if exist(fullfile(fld.basedir, 'results','figure_behavior','RTs.mat'),'file')
-    load(fullfile(fld.basedir, 'results','figure_behavior','RTs.mat'),'RTs');  
-else
-    error('RTs missing, run figure_behavior.m first');
-end
-
-%% plots
-% load traces
-monkeys = {'M1','M2'};
-wm = [];
-
-%% plot
-f1 = figure;
-set(f1,'Position',[500 500 1600 800])
-% individual animals
-for mi = 1:2
-    monkey = monkeys{mi};
-    savedir = fullfile(fld.basedir,'results','figure_neurotrace');
-    load(fullfile(savedir, [monkey '_averages_snr' num2str(snrthres) '_mindays' num2str(mindays) '.mat']),...
-        'traces','tracesLUT','grandtraces','grandtracesLUT','env_t');
-
-    traces_ALL{mi} = traces; %#ok<*AGROW,*NASGU> 
-    tracesLUT_ALL{mi} = tracesLUT;
-    grandtraces_ALL{mi} = grandtraces;
-    grandtracesLUT_ALL{mi} = grandtracesLUT;
-    env_t_ALL{mi} = env_t;
-    
-    figure(f1); subplot(2,3,mi+1); hold on;
-    mns = []; sems = []; trcs = cell(0);
-    for cond = 1:3
-                    
-        if do_cleantraces
-            incl = grandtracesLUT(:,2)==cond;
-            get_traces = grandtraces(incl,:);
-            sel_LUT = grandtracesLUT(incl,:);
-            incl_ALL{mi,cond} = incl;
-
-            prewin = env_t > -0.1 & env_t < 0;
-            postwin = env_t > 0 & env_t < 0.1;
-            sd_pre = std(get_traces(:,prewin),0,2);
-            m_post = mean(get_traces(:,postwin),2);
-            snr_sel = m_post > 3*sd_pre;
-
-            incl = incl(snr_sel');
-            get_traces = get_traces(snr_sel',:);
-            newLUT = sel_LUT(snr_sel',:);
-
-            % now average sites over sessions
-            get_traces_ = get_traces; % backup
-            get_traces = [];
-            for rs = unique(newLUT(:,1))'
-                if sum(newLUT(:,1)==rs)>1
-                    get_traces = [get_traces;...
-                        mean(get_traces_(newLUT(:,1)==rs,:))];
-                else
-                    get_traces = [get_traces; get_traces_(newLUT(:,1)==rs,:)];
-                end
-            end
-            recsites{mi,cond}=unique(newLUT(:,1));
-        else 
-            incl = tracesLUT(:,2)==cond;
-            get_traces = traces(incl,:);
-            sel_LUT = tracesLUT(incl,:);
-            incl_ALL{mi,cond} = incl;
-
-            newLUT = sel_LUT;
-            recsites{mi,cond}=unique(newLUT(:,1));
-        end
-
-        % calculate the mean across channels        
-        s_traces = [];
-        for t = 1:size(get_traces,1)
-            s_traces(t,:) = smooth(get_traces(t,:),smoothfact);
-        end
-        s_traces_ALL{mi,cond} = s_traces;
-
-        mn = mean(s_traces);
-        mns(cond,:) = mn; % for calculating the difference between conditions
-        trcs{cond} = s_traces; % for calculating the modulation onset
-        sem = std(s_traces)./sqrt(sum(incl));
-        sems(cond,:) = sem;
-
-        t2 = [env_t*1e3, fliplr(env_t*1e3)];
-        inBetween = [mn+sem, fliplr(mn-sem)];
-        fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
-            'FaceAlpha',falpha); hold all
-        plot(env_t*1e3,mn,'Color',cols(cond,:));
-        
-        RT=RTs{mi}(cond);
-        plot([RT RT],[0 1],'--','Color',cols(cond,:),'LineWidth',2)
-        
-        xlim([0 0.25*1e3]); ylim([0 1]);
-    end
-    title(monkeys{mi}); xlabel('time(ms)'); ylabel('MUA'); 
-
-
-    % only include sites with all conditions 
-    if size(recsites{mi,1},1) <= size(recsites{mi,2},1) && ...
-            size(recsites{mi,1},1) <= size(recsites{mi,3},1)
-        rsidx = [1 2 3];
-    elseif size(recsites{mi,2},1) < size(recsites{mi,1},1) && ...
-            size(recsites{mi,2},1) < size(recsites{mi,3},1)
-        rsidx = [2 1 3];
-    elseif size(recsites{mi,3},1) < size(recsites{mi,1},1) && ...
-            size(recsites{mi,3},1) < size(recsites{mi,2},1)
-        rsidx = [3 1 2];
-    end
-    rsinc{mi,rsidx(1)} = ismember(recsites{mi,rsidx(1)},recsites{mi,rsidx(1)});
-    rsinc{mi,rsidx(2)} = ismember(recsites{mi,rsidx(2)},recsites{mi,rsidx(1)});
-    rsinc{mi,rsidx(3)} = ismember(recsites{mi,rsidx(3)},recsites{mi,rsidx(1)});
-        
-    inRF = {'target','','distractor'};
-    for cond = [1 3 4]
-        figure(f1)
-        subplot(2,3,mi+4);
-        if cond < 4
-            df = trcs{cond}(rsinc{mi,cond},:) - ...
-                trcs{2}(rsinc{mi,2},:); % difference for each channel separately
-            mn = mean(df);
-            csem = std(trcs{cond}(rsinc{mi,cond},:) - ...
-                trcs{2}(rsinc{mi,2},:))./...
-                sqrt(size(trcs{cond}(rsinc{mi,cond},:),1));
-        else
-            df = trcs{1}(rsinc{mi,1},:) - ...
-                trcs{3}(rsinc{mi,3},:);
-            mn = mean(df);
-            csem = std(trcs{1}(rsinc{mi,1},:) - ...
-                trcs{3}(rsinc{mi,3},:))./...
-                sqrt(size(trcs{1}(rsinc{mi,1},:),1));
-        end
-            
-        inBetween = [mn+csem, fliplr(mn-csem)];
-        fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
-            'FaceAlpha',falpha); hold all
-        plot(env_t*1e3,mn,'Color',cols(cond,:));
-        xlim([0 0.25*1e3]); ylim([-0.15 0.25]);
-        
-        % let's calculate the latency of modulation per channel
-        if do_latency
-            do_latfig=0;
-            if mi==1
-                until = find(env_t<0.25,1,'last');
-            else
-                until = find(env_t<0.21,1,'last');
-            end
-            lat = []; coeff = []; rs = []; converge = [];
-            for i = 1:size(df,1) % loop channels
-                cur_resp = df(i,1:until)';
-                if cond==3
-                    % the latency function cannot fit downward modulation
-                    cur_resp = cur_resp*-1;
-                end
-                [lat(i),~] = latencyfit_jp(cur_resp,env_t(1:until)'*1000,do_latfig);
-            end
-
-            % latency on average
-            cur_resp = mn(1:until)';
-            if cond==3
-                % the latency function cannot fit downward modulation
-                cur_resp = cur_resp*-1;
-            end
-            lt = latencyfit_jp(smooth(cur_resp,20),env_t(1:until)'*1000,do_latfig);
-
-            if cond < 4
-                rng('default') % for reproducibility
-                nsamp = 100;
-                btstrp_LT = bootstrp(nsamp,@latencyfit_jp,smooth(cur_resp,20),env_t(1:until)'*1000,0);
-                fprintf([monkey ' ' inRF{cond} ' - Bootstrapped latency [' num2str(nsamp) ' samples] =====\n'])
-                fprintf(['mean: ' num2str(mean(btstrp_LT,'omitnan')) ', std:' num2str(std(btstrp_LT,'omitnan')) '\n'])
-                BLT{mi,cond}=btstrp_LT;
-            end
-
-            if ~isnan(lt) % if lt = nan, latencyfit_jp gives no figure, so we shouldn't make a title either
-                if cond < 4
-                    title([monkey ', fit on average response, ' inRF{cond} ' in RF']);
-                else
-                    title([monkey ', fit on average response modulation']);
-                end
-            end
-
-            if cond < 4
-                disp([monkey ', latency of modulation onset when ' inRF{cond} ' is in RF: ' ...
-                    num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
-                    'ms as estimated on the average response)']);
-            else
-                disp([monkey ', latency of targ-distr modulation: ' ...
-                    num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
-                    'ms as estimated on the average response)']);
-            end
-            %disp(['individual channel estimates: ' num2str(lat)]);
-
-            ModLatency{mi,cond} = lt;
-            figure(f1);
-            plot([lt lt],[-0.2 0.25],'--','Color',cols(cond,:),'LineWidth',2)
-            if isnan(lt)
-                plot([round(mean(lat,'omitnan')) round(mean(lat,'omitnan'))],[-0.2 0.25],...
-                    ':','Color',cols(cond,:),'LineWidth',2)
-            end
-        end
-    end
-    figure(f1)
-    plot([0 .25*1e3],[0 0],'k');
-    title(monkeys{mi}); xlabel('time(ms)'); ylabel('Modulation(MUA)'); 
-
-    % save the extra cleaned traces
-    save(fullfile(savedir,'Traces_SNR_trialbased.mat'),...
-        's_traces_ALL','recsites','rsinc');
-
-    %% statistics
-    window_means = [];
-    win_idx = find(env_t>0.15 & env_t<0.20);
-
-    for cond = 1:3
-        if do_cleantraces
-            incl = grandtracesLUT(:,2)==cond;
-            get_traces = grandtraces(incl,win_idx); %#ok<*FNDSB> 
-        else
-            incl = tracesLUT(:,2)==cond;
-            get_traces = traces(incl,win_idx);
-        end
-        incl = tracesLUT(:,2)==cond;
-
-        % calculate the average in a window
-        s_traces = [];
-        for t = 1:size(get_traces,1)
-            s_traces(t,:) = smooth(get_traces(t,:),smoothfact);
-        end
-        window_means(:,cond) = mean(s_traces,2);
-    end
-    anova1(window_means);
-    
-    wm = [wm; window_means];
-    wm2{mi} = window_means;
-end
-
-figure(f1); subplot(2,3,1); hold on;
-% pooled animals
-for cond = 1:3
-    traces = [s_traces_ALL{1,cond};s_traces_ALL{2,cond}];
-
-    % calculate the mean across channels
-    % get_traces = traces(incl,:);
-    get_traces = traces;
-    s_traces = [];
-    for t = 1:size(get_traces,1)
-        s_traces(t,:) = smooth(get_traces(t,:),smoothfact);
-    end
-
-    mn = mean(s_traces);
-    mns(cond,:) = mn; % for calculating the difference between conditions
-    trcs{cond} = s_traces; % for calculating the modulation onset
-    sem = std(s_traces)./sqrt(sum(size(s_traces,1)));
-    sems(cond,:) = sem;
-    
-    t2 = [env_t*1e3, fliplr(env_t*1e3)];
-    inBetween = [mn+sem, fliplr(mn-sem)];
-    fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
-        'FaceAlpha',falpha); hold all
-    
-    RT=RTs{3}(cond);
-    plot([RT RT],[0 1],'--','Color',cols(cond,:),'LineWidth',2)
-    plot(env_t*1e3,mn,'Color',cols(cond,:));
-    xlim([0 0.25*1e3]); ylim([0 1]);
-end
-title('BOTH animals pooled'); xlabel('time(ms)'); ylabel('MUA'); 
-
-inRF = {'target','','distractor'};
-
-for cond=1:3
-    rsinc{3,cond} = [rsinc{1,cond};rsinc{2,cond}];
-end
-
-for cond = [1 3 4]
-    figure(f1)
-    subplot(2,3,4);
-
-    if cond < 4
-        df = trcs{cond}(rsinc{3,cond},:) - ...
-            trcs{2}(rsinc{3,2},:); % difference for each channel separately
-        mn = mean(df);
-        csem = std(trcs{cond}(rsinc{3,cond},:) - ...
-            trcs{2}(rsinc{3,2},:))./...
-            sqrt(size(trcs{cond}(rsinc{3,cond},:),1));
-    else
-        df = trcs{1}(rsinc{3,1},:) - ...
-            trcs{3}(rsinc{3,3},:);
-        mn = mean(df);
-        csem = std(trcs{1}(rsinc{3,1},:) - ...
-            trcs{3}(rsinc{3,3},:))./...
-            sqrt(size(trcs{1}(rsinc{3,1},:),1));
-    end
-    
-    inBetween = [mn+csem, fliplr(mn-csem)];
-    fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
-        'FaceAlpha',falpha); hold all
-    plot(env_t*1e3,mn,'Color',cols(cond,:));
-    xlim([0 0.25*1e3]); ylim([-0.15 0.25]);
-    
-    % let's calculate the latency of modulation per channel
-    if do_latency
-        do_latfig=0;
-        if mi==1
-            until = find(env_t<0.25,1,'last');
-        else
-            until = find(env_t<0.21,1,'last');
-        end
-        lat = []; coeff = []; rs = []; converge = [];
-        for i = 1:size(df,1) % loop channels
-            cur_resp = df(i,1:until)';
-            if cond==3
-                % the latency function cannot fit downward modulation
-                cur_resp = cur_resp*-1;
-            end
-            [lat(i),~] = latencyfit_jp(cur_resp,env_t(1:until)'*1000,do_latfig);
-        end
-        
-        %collect latencty
-        clat{cond} = lat;
-
-        % latency on average
-        do_latfig=1;
-        cur_resp = mn(1:until)';
-        if cond==3
-            % the latency function cannot fit downward modulation
-            cur_resp = cur_resp*-1;
-        end
-        lt = latencyfit_jp(smooth(cur_resp,20),env_t(1:until)'*1000,do_latfig);
-
-        if cond < 4
-            rng('default') % for reproducibility
-            nsamp = 100;
-            btstrp_LT = bootstrp(nsamp,@latencyfit_jp,smooth(cur_resp,20),env_t(1:until)'*1000,0);
-            fprintf([inRF{cond} ' - Bootstrapped latency [' num2str(nsamp) ' samples] =====\n'])
-            fprintf(['mean: ' num2str(mean(btstrp_LT,'omitnan')) ', std:' num2str(std(btstrp_LT,'omitnan')) '\n'])
-            BLT{3,cond}=btstrp_LT;
-        end
-
-        if ~isnan(lt) % if lt = nan, latencyfit_jp gives no figure, so we shouldn't make a title either
-            if cond < 4
-                title(['BOTH monkeys, fit on average response, ' inRF{cond} ' in RF']);
-            else
-                title('BOTH monkeys, fit on average response modulation');
-            end
-        end
-        if cond < 4
-            disp(['BOTH monkeys, latency of modulation onset when ' inRF{cond} ' is in RF: ' ...
-                num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
-                'ms as estimated on the average response)']);
-        else
-            disp(['BOTH monkeys, latency of targ-distr modulation: ' ...
-                num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
-                'ms as estimated on the average response)']);
-        end
-        %disp(['individual channel estimates: ' num2str(lat)]);
-
-        figure(f1)
-        plot([lt lt],[-0.2 0.25],'--','Color',cols(cond,:),'LineWidth',2)
-        if isnan(lt)
-            plot([round(mean(lat,'omitnan')) round(mean(lat,'omitnan'))],[-0.2 0.25],...
-                ':','Color',cols(cond,:),'LineWidth',2)
-        end
-
-    end
-end
-figure(f1)
-plot([0 .25*1e3],[0 0],'k');
-title('ALL animals pooled'); xlabel('time(ms)'); ylabel('Modulation (MUA)'); 
-
-%% do a bootstrapped latency analysis to get some stats
-rng('default') % for reproducibility
-nsamp = 1000;
-btstrp_LAT = bootstrp(nsamp,@nanmean,[clat{1}' clat{3}']);
-[h,p,ci,stats]=ttest(btstrp_LAT(:,1)-btstrp_LAT(:,2));
-fprintf('== BTSTRP comparison of latencies ===\n');
-fprintf(['Mean difference in latency (T-SD): ' ...
-    num2str(mean(btstrp_LAT(:,1)-btstrp_LAT(:,2))) ', SEM: ' ...
-    num2str(std(btstrp_LAT(:,1)-btstrp_LAT(:,2))./sqrt(size(btstrp_LAT,1))) ...
-    'ms\n']);
-fprintf([num2str(sum((btstrp_LAT(:,1)-btstrp_LAT(:,2))<0)) '/' num2str(nsamp) ...
-    ' T before SD, p = ' num2str(sum((btstrp_LAT(:,1)-btstrp_LAT(:,2))>0)/nsamp) ...
-    ' that H0[T not before SD] is rejected \n']);
-fprintf(['ttest t(' num2str(stats.df) ') = ' num2str(stats.tstat) ...
-    ', p = ' num2str(p) '\n']);
-
-%% use the bootstrapped latency estimates for stats
-fprintf('== Comparison of bootstrapped T-SD latencies ===\n');
-for mi=1:3
-    [h,p,ci,stats]=ttest(BLT{mi,1}-BLT{mi,3});
-    nn = sum(~isnan());
-    if mi<3
-        fprintf(['-- Monkey ' num2str(mi) ' --\n']);
-    else
-        fprintf('-- BOTH Monkeys --\n');
-    end
-    fprintf(['Mean difference in latency (T-SD): ' ...
-        num2str(mean(BLT{mi,1}-BLT{mi,3},'omitnan')) ', SEM: ' ...
-        num2str(std(BLT{mi,1}-BLT{mi,3},'omitnan')./sqrt(sum(~isnan(BLT{mi,1})))) ...
-        'ms\n']);
-    fprintf([num2str(sum((BLT{mi,1}-BLT{mi,3})<0)) '/' num2str(100) ...
-        ' T before SD, p = ' num2str(sum((BLT{mi,1}-BLT{mi,3})>0)/100) ...
-        ' that H0[T not before SD] is rejected \n']);
-    fprintf(['ttest t(' num2str(stats.df) ') = ' num2str(stats.tstat) ...
-        ', p = ' num2str(p) '\n']);
-end
-
-%% statistics
-% do a simple paired t-test for each monkey
-fprintf('=============================\n');
-fprintf(' Paired T-test \n');
-fprintf('=============================\n');
-clear h p stats
-for mi = 1:2
-    mns = wm2{mi};
-    p = []; vp=[]; ci=1;
-    for cond = [1 3]
-        [vh,vp(end+1)]=vartest2(mns(:,cond),mns(:,2));
-        [h,p(end+1),~,stats{mi,ci}] = ttest(mns(:,cond), mns(:,2));
-        ci=ci+1;
-    end
-    [vh,vp(end+1)]=vartest2(mns(:,1),mns(:,3));
-    [h,p(end+1),~,stats{mi,3}] = ttest(mns(:,1), mns(:,3));
-        
-    disp([monkeys{mi} ', target response increased: p=' num2str(p(1)) ...
-        ' t = ' num2str(stats{mi,1}.tstat) ', df = ' num2str(stats{mi,1}.df) ... 
-        ' (unequal var p = ' num2str(vp(1)) ')']);
-    disp([monkeys{mi} ', distractor response decreased: p=' num2str(p(2)) ...
-        ' t = ' num2str(stats{mi,2}.tstat) ', df = ' num2str(stats{mi,2}.df) ... 
-        ' (unequal var p = ' num2str(vp(2)) ')']);
-    disp([monkeys{mi} ', targ vs distractor: p=' num2str(p(3)) ...
-        ' t = ' num2str(stats{mi,3}.tstat) ', df = ' num2str(stats{mi,3}.df) ... 
-        ' (unequal var p = ' num2str(vp(3)) ')']);
-end
-
-%Both together
-mns = [wm2{1}; wm2{2}];
-p = [];vp=[]; ci=1;
-for cond = [1 3]
-    [vh,vp(end+1)]=vartest2(mns(:,cond),mns(:,2));
-    [h,p(end+1),~,stats{3,ci}] = ttest(mns(:,cond), mns(:,2));
-    ci=ci+1;
-end
-[vh,vp(end+1)]=vartest2(mns(:,1),mns(:,3));
-[h,p(end+1),~,stats{3,3}] = ttest(mns(:,1), mns(:,3));
-
-disp(['BOTH, target response increased: p=' num2str(p(1)) ...
-    ' t = ' num2str(stats{3,1}.tstat) ', df = ' num2str(stats{3,1}.df) ... 
-    ' (unequal var p = ' num2str(vp(1)) ')']);
-disp(['BOTH, distractor response decreased: p=' num2str(p(2)) ...
-    ' t = ' num2str(stats{3,2}.tstat) ', df = ' num2str(stats{3,2}.df) ... 
-    ' (unequal var p = ' num2str(vp(2)) ')']);
-disp(['BOTH, targ vs distractor: p=' num2str(p(3)) ...
-    ' t = ' num2str(stats{3,3}.tstat) ', df = ' num2str(stats{3,3}.df) ... 
-    ' (unequal var p = ' num2str(vp(3)) ')']);
-
-
-%% save figure
-figure(f1)
-saveas(gcf,fullfile(savedir,'figure_neurotrace.svg'));