Upload files to 'NCM-errors'

NCM-errors/NCM-errors-analysis-code/ncm-error-description

@@ -0,0 +1,5 @@
+This folder contains figure scripts for reproducing figure pieces from NCM error chapter of dissertation
+
+There is no code for figure 2
+
+Note that file paths will need to be changed! 

NCM-errors/NCM-errors-analysis-code/ncm_errors_error_detection_pipeline.m

@@ -0,0 +1,218 @@
+% it is very time-consuming to perform the analyses that determine whether
+% neurons are responsive to perturbations of different types (as this
+% involves over 50 distinct rank-sum calculations per neuron).
+%
+% to save time in analysis, we run this code one time after running the
+% prelim pipeline on each bird. the result is large multi-bird data
+% structures summarizing error-responsiveness data
+
+the_good_birds = {'D:\error_expmts\pipeline_outputs\pipeline_7635_210723_LH_NCM_g0.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_9007_210722_LH_NCM_g0.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_8975_210610_LH_NCM_g1.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_9047_210722_LH_NCM_g0.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_9049_210724_LH_NCM_g0.mat'
+    };
+
+pipeline_output_save_target = 'D:\error_expmts\multibird_summary_ncm_errors.mat';
+
+big_aud_bcs = [];
+big_rare_dev_iti = [];
+big_rare_dev_std = [];
+big_fam_late_iti = [];
+big_fam_late_std = [];
+
+big_rdi_mat = [];
+big_fli_mat = [];
+big_std_mat = [];
+big_rdi_ids =[];
+big_fli_ids = [];
+big_std_ids = [];
+big_resp_mat = []; % we'll use this to calculate indices for errors (each row is neuron x pert combo)
+brm_mean_norm = [];
+
+big_adapt_resp_raw = [];
+big_adapt_resp_bcs = [];
+
+big_hab_adapt_raw = [];
+big_odd_adapt_raw = [];
+
+big_aud_resp_traj = [];
+all_ids = 0;
+total_ids = [];
+for big_bird_num = 1:length(the_good_birds)
+    clearvars -except big_bird_num the_good_birds big_aud_bcs ...
+        big_rare_dev_iti big_fam_late_iti all_ids total_ids...
+        big_rare_dev_std big_fam_late_std big_adapt_resp_bcs...
+        big_adapt_resp_raw big_rdi_mat big_fli_mat big_std_mat...
+        big_hab_adapt_raw big_odd_adapt_raw big_aud_resp_traj big_resp_mat...
+        big_rdi_ids big_fli_ids big_std_ids pipeline_output_save_target ...
+        brm_mean_norm
+    close all
+    
+    this_bird = the_good_birds{big_bird_num};
+    load(this_bird)
+        
+    disp(['starting bird ' num2str(big_bird_num)])
+    
+    
+    %% sound feature sensitivity
+    
+    iti_resp_aud = iti_resp(1:5,good_unit_ids);
+    iti_std = std(iti_resp_aud, [], 1);
+    iti_mean = mean(iti_resp_aud, 1);
+    
+    iti_ceil_aud = iti_mean + 1.96*iti_std;
+    iti_floor_aud = iti_mean - 1.96*iti_std;
+    
+    aud_resp_bc = zeros(length(unique_sylls),length(good_unit_ids));
+    adapt_resp_bc = zeros(length(unique_sylls),length(good_unit_ids));
+    for sIdx = 1:length(unique_sylls)
+        this_syll = unique_sylls{sIdx};
+        
+        syll_resp = reshape(all_trial_resp,size(all_trial_resp,1)*size(all_trial_resp,2),size(all_trial_resp,3));
+        syll_resp = syll_resp(:, good_unit_ids);
+        
+        these_trials_x = find(strcmp(syllableInfo(1:2375,3),this_syll));
+        these_trials = these_trials_x(1:5);
+        these_final_trials = these_trials_x(end-4:end);
+        
+        if sIdx >= 19
+            if strcmp(stim_date, '210610')
+                hab_trials_x = find(strcmp(syllableInfo(12476:14475,3),this_syll));
+                odd_trials_x = find(strcmp(syllableInfo(2376:12475,3),this_syll));
+                last_pre_hab = 12475;
+                last_pre_odd = 2375;
+            elseif strcmp(stim_date, '210613')
+                hab_trials_x = find(strcmp(syllableInfo(12436:14435,3),this_syll));
+                odd_trials_x = find(strcmp(syllableInfo(2376:12435,3),this_syll));
+                last_pre_hab = 12435;
+                last_pre_odd = 2375;
+            end
+            hab_trials = hab_trials_x(1:5);
+            hab_final_trials = hab_trials_x(end-4:end);
+            
+            odd_trials = odd_trials_x(1:5);
+            
+            
+            final_hab = syll_resp(last_pre_hab+hab_final_trials,:);
+            mean_hab_final = mean(final_hab,1);
+            these_hab = syll_resp(last_pre_hab+hab_trials,:);
+            mean_hab_init = mean(these_hab,1);
+            
+            odd_init = syll_resp(last_pre_odd+odd_trials,:);
+            mean_odd_init = mean(odd_init,1);
+            
+            hab_adapt_raw(sIdx,:) = (mean_hab_init - mean_hab_final) ./ (mean_hab_init + mean_hab_final);
+            odd_adapt_raw(sIdx,:) = (mean_odd_init - mean_hab_final) ./ (mean_odd_init + mean_hab_final);
+        end
+        
+        final_resp = syll_resp(these_final_trials,:);
+        mean_final = mean(final_resp,1);
+        ceil_final = mean_final + 1.96*std(final_resp,[],1);
+        floor_final = mean_final - 1.96*std(final_resp,[],1);
+        
+        
+        these_resp = syll_resp(these_trials,:);
+        mean_resp = mean(these_resp, 1);
+        this_std_resp = std(these_resp, [], 1);
+        ceil_resp = mean_resp + 1.96*this_std_resp;
+        floor_resp = mean_resp - 1.96*this_std_resp;
+        
+        this_resp_cond = (floor_resp > iti_ceil_aud) | (ceil_resp < iti_floor_aud);
+        
+        aud_resp_bc(sIdx,:) = this_resp_cond;
+        
+        adapt_resp_raw(sIdx,:) = (mean_resp - mean_final) ./ (mean_resp + mean_final);
+        adapt_resp_bc(sIdx,:) = (floor_resp > ceil_final); % 1 if end and start don't overlap
+        
+        if sIdx >= 19
+            big_aud_resp_traj = [big_aud_resp_traj; ...
+                mean_resp' mean_final' mean_odd_init' mean_hab_init' mean_hab_final'];
+        end
+            
+    end
+    
+    big_hab_adapt_raw = [big_hab_adapt_raw hab_adapt_raw];
+    big_odd_adapt_raw = [big_odd_adapt_raw odd_adapt_raw];
+
+
+    
+    big_aud_bcs = [big_aud_bcs aud_resp_bc];
+    
+    big_adapt_resp_raw = [big_adapt_resp_raw adapt_resp_raw];
+    big_adapt_resp_bcs = [big_adapt_resp_bcs adapt_resp_bc];
+    
+    big_rare_dev_iti = [big_rare_dev_iti rare_dev_bc_iti];
+    big_rare_dev_std = [big_rare_dev_std rare_dev_bc_std];
+    
+    big_fam_late_iti = [big_fam_late_iti fam_late_bc_iti];
+    big_fam_late_std = [big_fam_late_std fam_late_bc_std];
+    
+    
+    std_aud_resp = find(sum(aud_resp_bc(19:23,:),1) > 0);
+    
+    %% big matrices of error responses
+    
+    sig_thresh = .05/18; % bonferroni-corrected sig val for detection
+       
+    for unitIdx = 1:length(good_unit_ids)
+        disp(['pert ids: unit ' num2str(unitIdx) ' of ' num2str(length(good_unit_ids)) '.'])
+        this_id = good_unit_ids(unitIdx);
+        this_prepert_iti = squeeze(arbt_prepert_iti(:,:,this_id));
+        this_pert_iti = squeeze(arbt_pert_iti(:,:,this_id));
+        this_postpert_iti = squeeze(arbt_postpert_iti(:,:,this_id));
+        
+        
+        for pIdx = 1:18
+            mw_iti = ranksum(this_prepert_iti(pIdx,:), this_pert_iti(pIdx,:), 'method', 'exact', 'tail', 'left');
+            this_erumd = mw_iti < sig_thresh; % 
+            
+            mw_iti = ranksum(this_prepert_iti(pIdx,1:10), this_prepert_iti(pIdx,end-9:end), 'method', 'exact');
+            this_intra_er = mw_iti < sig_thresh;
+            
+            mw_iti = ranksum(this_postpert_iti(pIdx,:), this_pert_iti(pIdx,:), 'method', 'exact', 'tail', 'left');
+            this_ltumd = mw_iti < sig_thresh;
+            
+            mw_iti = ranksum(this_postpert_iti(pIdx,1:10), this_postpert_iti(pIdx,end-9:end), 'method', 'exact');
+            this_intra_lt = mw_iti < sig_thresh;
+            
+            mw_iti = ranksum(this_prepert_iti(pIdx,:), this_postpert_iti(pIdx,:), 'method', 'exact', 'tail', 'left');
+            this_erult = mw_iti < sig_thresh;
+            
+            mw_iti = ranksum(this_pert_iti(pIdx,:), this_postpert_iti(pIdx,:), 'method', 'exact', 'tail', 'left');
+            this_mdult = mw_iti < sig_thresh;
+            
+            mw_iti = ranksum(this_pert_iti(pIdx,1:10), this_pert_iti(pIdx,end-9:end), 'method', 'exact');
+            this_intra_md = mw_iti < sig_thresh;
+            
+            
+            this_resp_patt = [this_prepert_iti(pIdx,:) this_pert_iti(pIdx,:) ...
+                this_postpert_iti(pIdx,:)];
+            
+            this_resp_patt = this_resp_patt ./ max(abs(this_resp_patt));
+            trp_mean_norm = this_resp_patt ./ mean(this_prepert_iti(pIdx,:));
+            
+            all_ids = all_ids + 1;
+    
+            if (this_mdult && ~this_intra_md && this_erult && ~this_intra_er)
+                big_fli_mat = [big_fli_mat; this_resp_patt];
+                big_fli_ids = [big_fli_ids; all_ids];
+            elseif (this_erumd && ~this_intra_er && this_ltumd && ~this_intra_lt)
+                big_rdi_mat = [big_rdi_mat; this_resp_patt];
+                big_rdi_ids = [big_rdi_ids; all_ids];
+            elseif (~this_intra_er && ~this_intra_md && ~this_intra_lt)
+                big_std_mat = [big_std_mat; this_resp_patt];
+                big_std_ids = [big_std_ids; all_ids];
+            end
+
+            big_resp_mat = [big_resp_mat; this_resp_patt];
+            total_ids = [total_ids; all_ids];
+            
+            brm_mean_norm = [brm_mean_norm; trp_mean_norm];
+
+        end
+    end
+
+end
+
+save(pipeline_output_save_target)

NCM-errors/NCM-errors-analysis-code/ncm_errors_fig1_ssa_recap.m

@@ -0,0 +1,162 @@
+% 1. load prelim bird data and plot summary of ssa indices
+% 2. for specific single unit, generate motif-wise SSA psth
+
+% notes:
+%   for each neuron, for each song, ssa index is 
+%   (avg motif response during first bout - avg motif response during last bout) / (avg during first bout)
+%   then, for each neuron, we take the max ssa_index across songs (because
+%   different neurons habituate differently to different songs)
+
+%% prep and set variables
+
+clear
+close all
+prelim_ssa_bird = 'D:\ssa_expmts\pipelines\pipeline_9059_210811_LH_NCM_g0.mat';
+
+psth_unit_id = 25;
+psth_unit_type = 'su';
+
+
+%% ssa index summary    
+
+load(prelim_ssa_bird)
+
+big_ssa_mat = [];
+
+p1_resp = all_motif_resp(1:500,:);
+init_resps = p1_resp([1:5 126:130 251:255 376:380],:);
+max_resps = max(init_resps,[],1);
+
+% eliminate units with no initial resp to any song, to avoid infs
+% NB: this misses units that potentiate, but we're not looking at those
+% for this basic confirmatory analysis
+valid_ids = find(max_resps > 0);
+p1_resp_valid = p1_resp(:,valid_ids);
+max_resps_valid = max_resps(valid_ids);
+
+% normalize responses to strongest initial response
+p1_norm = p1_resp_valid ./ max_resps_valid;
+
+% get responses for each unit for each of 4 std songs
+s1p1_resp = p1_resp(1:125,valid_ids);
+s2p1_resp = p1_resp(126:250,valid_ids);
+s3p1_resp = p1_resp(251:375,valid_ids);
+s4p1_resp = p1_resp(376:500,valid_ids);
+
+
+for unitIdx = 1:length(valid_ids)
+    this_ssa_vec = [];
+    for songIdx = 1:4 % at least one will be nonzero because of valid_units
+        eval(['this_resp = s' num2str(songIdx) 'p1_resp(:,unitIdx);'])
+        this_init = mean(this_resp(1:5));
+        this_end = mean(this_resp(121:125));
+        test_ssa = (this_init-this_end)/this_init;
+        if ~isinf(test_ssa)
+            this_ssa_vec = [this_ssa_vec; test_ssa];
+        else
+            this_ssa_vec = [this_ssa_vec; 0];
+        end
+    end
+    big_ssa_mat = [big_ssa_mat; max(this_ssa_vec)]; %max(this_ssa_vec) or not
+end
+
+
+h=figure; 
+histogram(max(big_ssa_mat,0),20)
+xlabel('SSA Index')
+ylabel('Number of Neurons')
+title('Distribution of SSA Index')
+xlim([-0.2 1.2])
+
+print('figure_pieces/fig1_ssaIdx', '-dsvg', '-r300')
+saveas(h, 'figure_pieces/fig1_ssaIdx.fig')
+
+
+%% plot psth
+
+% unoptimized and counter-intuitive, sorry. 
+
+% transition times for the stim
+s1s2 = 125;
+s2s3 = 250;
+s3s4 = 375;
+s4s1b = 500;
+s1bs2b = 550;
+s2bs3b = 600;
+s3bs4b = 650;
+std_end = 700;
+
+if strcmp(psth_unit_type, 'su')
+    this_row = all_motif_resp(1:std_end,psth_unit_id);
+else
+    this_row = all_motif_resp_mu(1:std_end,psth_unit_id);
+end
+
+h2=figure(1006);
+subplot(4,2,[1 2 3 4])
+plot(1:s1s2, this_row(1:s1s2), 'r-')
+hold on;
+plot(s1s2+1:s1s2+50, this_row(s4s1b+1:s1bs2b), 'r.')
+newstart = s1s2+51;
+next_piece = this_row(s1s2+1:s2s3);
+newend = newstart+length(next_piece)-1;
+plot(newstart:newend, this_row(s1s2+1:s2s3), 'b-')
+
+newstart = newend+1;
+next_piece = this_row(s1bs2b+1:s2bs3b);
+newend = newstart+length(next_piece)-1;
+plot(newstart:newend, next_piece, 'b.')
+
+newstart = newend+1;
+next_piece = this_row(s2s3+1:s3s4);
+newend = newstart+length(next_piece)-1;
+plot(newstart:newend, next_piece, 'g-')
+
+newstart = newend+1;
+next_piece = this_row(s2bs3b+1:s3bs4b);
+newend = newstart+length(next_piece)-1;
+plot(newstart:newend, next_piece, 'g.')
+
+newstart = newend+1;
+next_piece = this_row(s3s4+1:s4s1b);
+newend = newstart+length(next_piece)-1;
+plot(newstart:newend, next_piece, 'k-')
+
+newstart = newend+1;
+next_piece = this_row(s3bs4b+1:std_end);
+newend = newstart+length(next_piece)-1;
+plot(newstart:newend, next_piece, 'k.')
+
+hold off;
+ylabel('resp strength (hz)')
+title(['Std Motif Responses - ' num2str(this_id)])
+xlim([0 std_end])
+
+subplot(4,2,5)
+plot(this_row(1:25), 'r*-')
+xlim([1 25])
+title('First Responses to Std 1')
+ylim([min(0, min(this_row)) max(this_row)+10])
+
+subplot(4,2,6)
+plot(this_row(s1s2+1:s1s2+25), 'b*-')
+xlim([1 25])
+title('First Responses to Std 2')
+ylim([min(0, min(this_row)) max(this_row)+10])
+
+subplot(4,2,7)
+plot(this_row(s2s3+1:s2s3+25), 'g*-')
+xlim([1 25])
+title('First Responses to Std 3')
+ylim([min(0, min(this_row)) max(this_row)+10])
+
+subplot(4,2,8)
+plot(this_row(s3s4+1:s3s4+25), 'k*-')
+xlim([1 25])
+title('First Responses to Std 4, p1')
+ylim([min(0, min(this_row)) max(this_row)+10])
+
+print('figure_pieces/fig1_ssaExample', '-dsvg', '-r300')
+saveas(h2, 'figure_pieces/fig1_ssaExample.fig')
+
+disp('all done!')

NCM-errors/NCM-errors-analysis-code/ncm_errors_fig3_learned_error_summary.m

@@ -0,0 +1,246 @@
+% 1. load error detection pipeline output
+% 2. plot learned error summary imagesc
+% 3. across all birds, generate learned error response index histogram
+% 4. for given bird,unit_type,generate example raster (probably need to
+%   re-load individual bird pipeline output for this)
+
+
+%% prep and set parameters
+clear
+close all
+error_detection_pipeline = 'D:\error_expmts\multibird_summary_ncm_errors.mat';
+
+example_info = [1, 390, 15; 3, 458, 11]; % bird, unit, perturbation. as many as you want
+
+load(error_detection_pipeline)
+
+rng(23) % just affects the sorting of summary fig
+%% plot learned error imagesc
+% extremely straightforward - all the heavy lifting is done in the error
+% detection script
+h=figure(4);
+
+imagesc(flipud(big_fli_mat(randperm(length(big_fli_mat(:,1))),:)'))
+title('Familiar Pert Summary')
+ylabel('Trial')
+xlabel('Neuron x Pert Type Combo')
+
+print('figure_pieces/fig3_learnedSummary', '-dsvg', '-r300')
+saveas(h, 'figure_pieces/fig3_learnedSummary.fig')
+
+%% plot learned error response indices - color red for significant ids
+
+a1_total = mean(big_resp_mat(:, 1:25),2);
+odd_total = mean(big_resp_mat(:,26:50),2);
+a2_total = mean(big_resp_mat(:,51:75),2);
+
+err_idx = a2_total - ((a1_total +  odd_total)/2);
+[cts, edges] = histcounts(err_idx);
+
+
+h2=figure(20);
+histogram(err_idx, 'BinEdges', edges)
+hold on;
+histogram(err_idx(big_fli_ids),'BinEdges',edges)
+hold off;
+title('Learned Error Indices - All Interactions')
+
+print('figure_pieces/fig3_learnedIdxSummary', '-dsvg', '-r300')
+saveas(h2, 'figure_pieces/fig3_learnedIdxSummary.fig')
+
+
+% now we want 'max per neuron'
+% find max pert response per neuron, given that typically neuron responds
+% to just 1
+num_neurons = length(err_idx)/18;
+ei2 = zeros(num_neurons,1);
+
+sig_nrn_ids = [];
+for nIdx = 1:num_neurons
+    start_idx = ((nIdx-1)*18)+1;
+    end_idx = start_idx + 17;
+    
+    sig_nrn_bool = ~isempty(intersect(find(big_fli_ids >= start_idx), find(big_fli_ids <= end_idx)));
+    if sig_nrn_bool
+        sig_nrn_ids = [sig_nrn_ids; nIdx];
+    end
+    
+    ei2(nIdx) = max(err_idx(start_idx:end_idx));
+end
+
+[cts2, edges2] = histcounts(ei2);
+h3=figure(15);
+histogram(ei2, 'BinEdges', edges2)
+hold on;
+histogram(ei2(sig_nrn_ids),'BinEdges',edges2)
+hold off;
+title('Learned Error Indices - Max per Neuron')
+
+print('figure_pieces/fig3_learnedIdxSummary_mpn', '-dsvg', '-r300')
+saveas(h3, 'figure_pieces/fig3_learnedIdxSummary_mpn.fig')
+
+
+%% plot raster examples
+% nb: this script is more involved
+
+for ex_idx = 1:length(example_info(:,1))
+    this_bird = example_info(ex_idx,1);
+    load(the_good_birds{this_bird});
+    
+    %     pClass = {'_W', '_I'};
+    %     pType = {'S', 'H', 'N'};
+    
+    % window is slightly before motif to slightly after
+    
+    motif_1_length_secs = motif_length;
+    timeWindow_1 = [-0.1 motif_1_length_secs+0.225]; % see much after because we perturb final syllable
+    
+    pids = exp_info.pertSyllIDs_1;
+    
+    pert_time = exp_info.sub_inserts_1; % % same for both!
+    bStart = (cell2mat(syllableInfo(pids(1),1)) - cell2mat(syllableInfo(1,1)));
+    dStart = (cell2mat(syllableInfo(pids(2),1)) - cell2mat(syllableInfo(1,1)));
+    fStart = (cell2mat(syllableInfo(pids(3),1)) - cell2mat(syllableInfo(1,1)));
+    pert_time_1 = pert_time + [bStart dStart fStart];
+    pert_time_2 = 0 + [bStart dStart fStart];
+    
+    bLen = cell2mat(syllableInfo(pids(1),4)) - cell2mat(syllableInfo(pids(1),1));
+    dLen = cell2mat(syllableInfo(pids(2),4)) - cell2mat(syllableInfo(pids(2),1));
+    fLen = cell2mat(syllableInfo(pids(3),4)) - cell2mat(syllableInfo(pids(3),1));
+    
+    
+    
+    if strcmp(stim_date, '210610')
+        bl_1 = 1;
+        pert = 476;
+        hab = 2496;
+        bl_2 = 2896;
+    elseif strcmp(stim_date, '210613')
+        bl_1 = 1;
+        pert = 476;
+        hab = 2488;
+        bl_2 = 2888;
+    end
+    
+    target_pert = pert_types{example_info(ex_idx,3)};
+    target_syll = target_pert(2);
+    target_type = target_pert(4);
+    target_class = target_pert(5);
+    
+    switch target_syll
+        case 'b'
+            pert_time_2 = pert_time_2(1);
+            syll_lens_1 = bLen;
+        case 'd' 
+            pert_time_2 = pert_time_2(2);
+            syll_lens_1 = dLen;
+        case 'f'
+            pert_time_2 = pert_time_2(3);
+            syll_lens_1 = fLen;
+    end
+    
+    
+    
+    
+    relevant_syll = syllableInfo; % this is perturb phase
+    motif_starts = cell2mat(motifInfo(:,1));
+    syllTypes = syllableInfo(:,3);
+    mTypes = motifInfo(:,2);
+    unique_types = unique(mTypes);
+    
+    % we now need to build structure with start times of each relevant
+    % trial
+    
+    % what are relevant trials?
+    relevant_types = {'@@@@@'}; % this just needs to be something absent from strings so that first cell is empty
+    for uIdx = 1:length(unique_types)
+        thisU = unique_types{uIdx};
+        if contains(thisU, target_pert) % so is okay if target type doesn't exist...
+            relevant_types{end+1} = thisU;
+        end
+    end
+    
+    
+    start_times = cell(length(relevant_types),1);
+    
+    % first prehab
+    for mIdx = bl_1:pert-1
+        prevStart = motif_starts(19*(floor((mIdx-1)/19)) + 1);
+        
+        thisStart = motif_starts(mIdx);
+        thisType = mTypes{mIdx};
+        typeNum = [];
+        relIdx = 1;
+        while isempty(typeNum) && (relIdx <= length(relevant_types))
+            if contains(thisType, relevant_types{relIdx})
+                typeNum = relIdx;
+            end
+            relIdx = relIdx+1;
+        end
+        if ~isempty(typeNum)
+            start_times{typeNum} = [start_times{typeNum}; thisStart];
+            if typeNum == 2 % only add std once
+                start_times{1} = [start_times{1}; prevStart]; % add previous standard motif
+            end
+        end
+    end
+    
+    
+    % pert phase
+    for mIdx = pert:hab-1
+        if mIdx-pert > 1
+            prevStart = motif_starts(mIdx-1);
+        end
+        thisStart = motif_starts(mIdx);
+        thisType = mTypes{mIdx};
+        typeNum = [];
+        relIdx = 1;
+        while isempty(typeNum) && (relIdx <= length(relevant_types))
+            if contains(thisType, relevant_types{relIdx})
+                typeNum = relIdx;
+            end
+            relIdx = relIdx+1;
+        end
+        if ~isempty(typeNum)
+            start_times{typeNum} = [start_times{typeNum}; thisStart];
+            start_times{1} = [start_times{1}; prevStart]; % add previous standard motif
+        end
+    end
+    
+    
+    % second prehab
+    for mIdx = bl_2:length(motifTimingData(:,1))
+        prevStart = motif_starts(19*(floor((mIdx-1)/19)) + 1);
+        
+        thisStart = motif_starts(mIdx);
+        thisType = mTypes{mIdx};
+        typeNum = [];
+        relIdx = 1;
+        while isempty(typeNum) && (relIdx <= length(relevant_types))
+            if contains(thisType, relevant_types{relIdx})
+                typeNum = relIdx;
+            end
+            relIdx = relIdx+1;
+        end
+        if ~isempty(typeNum)
+            start_times{typeNum} = [start_times{typeNum}; thisStart];
+            if typeNum == 2 % only add std once
+                start_times{1} = [start_times{1}; prevStart]; % add previous standard motif
+            end
+        end
+    end
+    
+    
+    
+    this_id = (example_info(ex_idx,2)); % not clearly grabbing correct unit? or else wrong pert
+    theseSpikes = sp.st(sp.clu==goodUnits(this_id));
+
+    perturb_response_plot_psth_vF(theseSpikes, start_times, this_id, timeWindow_1, micData, target_type, pert_time_2, syll_lens_1, 1, ex_idx)
+    
+    h5 = gcf;
+    print(['figure_pieces/fig3_learnedExample' '_' num2str(ex_idx)], '-dsvg', '-r300')
+    saveas(h5, ['figure_pieces/fig3_learnedExample'  '_' num2str(ex_idx) '.fig'])
+end
+
+
+disp('i think im done')

NCM-errors/NCM-errors-analysis-code/ncm_errors_fig4_oddball_error_summary.m

@@ -0,0 +1,244 @@
+% 1. load error detection pipeline output
+% 2. plot learned error summary imagesc
+% 3. across all birds, generate learned error response index histogram
+% 4. for given bird,unit_type,generate example raster (probably need to
+%   re-load individual bird pipeline output for this)
+
+
+%% prep and set parameters
+clear
+close all
+error_detection_pipeline = 'D:\error_expmts\multibird_summary_ncm_errors.mat';
+
+example_info = [2, 11, 11; 1, 280, 6]; % bird, unit, perturbation. as many as you want
+
+load(error_detection_pipeline)
+
+rng(23) % just affects the sorting of summary fig
+%% plot learned error imagesc
+% extremely straightforward - all the heavy lifting is done in the error
+% detection script
+figure(8);
+imagesc(flipud(big_rdi_mat(randperm(length(big_rdi_mat(:,1))),:)'))
+title('Rare Dev. Summary')
+ylabel('Trial')
+xlabel('Neuron x Pert Type Combo')
+
+h5 = gcf;
+print('figure_pieces/fig4_oddballSummary', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig4_oddballSummary' '.fig'])
+
+
+%% plot learned error response indices - color red for significant ids
+
+a1_total = mean(big_resp_mat(:, 1:25),2);
+odd_total = mean(big_resp_mat(:,26:50),2);
+a2_total = mean(big_resp_mat(:,51:75),2);
+
+odd_idx = odd_total - ((a1_total + a2_total)/2);
+[cts, edges] = histcounts(odd_idx);
+
+
+figure(15);
+histogram(odd_idx, 'BinEdges', edges)
+hold on;
+histogram(odd_idx(big_rdi_ids),'BinEdges',edges)
+hold off;
+title('Oddball Error Indices - All Interactions')
+h5 = gcf;
+print(['figure_pieces/fig4_oddballIdx'], '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig4_oddballIdx' '.fig'])
+
+% now we want 'max per neuron'
+% find max pert response per neuron, given that typically neuron responds
+% to just 1
+num_neurons = length(odd_idx)/18;
+oi2 = zeros(num_neurons,1);
+
+sig_nrn_ids = [];
+for nIdx = 1:num_neurons
+    start_idx = ((nIdx-1)*18)+1;
+    end_idx = start_idx + 17;
+    
+    sig_nrn_bool = ~isempty(intersect(find(big_rdi_ids >= start_idx), find(big_rdi_ids <= end_idx)));
+    if sig_nrn_bool
+        sig_nrn_ids = [sig_nrn_ids; nIdx];
+    end
+    
+    oi2(nIdx) = max(odd_idx(start_idx:end_idx));
+end
+
+[cts2, edges2] = histcounts(oi2);
+figure(7);
+histogram(oi2, 'BinEdges', edges2)
+hold on;
+histogram(oi2(sig_nrn_ids),'BinEdges',edges2)
+hold off;
+title('Oddball Error Indices - Max per Neuron')
+
+h5 = gcf;
+print(['figure_pieces/fig4_oddballIdx_mpn'], '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig4_oddballIdx_mpn' '.fig'])
+
+
+
+%% plot raster examples
+% nb: this script is more involved
+
+for ex_idx = 1:length(example_info(:,1))
+    this_bird = example_info(ex_idx,1);
+    load(the_good_birds{this_bird});
+    
+    %     pClass = {'_W', '_I'};
+    %     pType = {'S', 'H', 'N'};
+    
+    % window is slightly before motif to slightly after
+    
+    motif_1_length_secs = motif_length;
+    timeWindow_1 = [-0.1 motif_1_length_secs+0.225]; % see much after because we perturb final syllable
+    
+    pids = exp_info.pertSyllIDs_1;
+    
+    pert_time = exp_info.sub_inserts_1; % % same for both!
+    bStart = (cell2mat(syllableInfo(pids(1),1)) - cell2mat(syllableInfo(1,1)));
+    dStart = (cell2mat(syllableInfo(pids(2),1)) - cell2mat(syllableInfo(1,1)));
+    fStart = (cell2mat(syllableInfo(pids(3),1)) - cell2mat(syllableInfo(1,1)));
+    pert_time_1 = pert_time + [bStart dStart fStart];
+    pert_time_2 = 0 + [bStart dStart fStart];
+    
+    bLen = cell2mat(syllableInfo(pids(1),4)) - cell2mat(syllableInfo(pids(1),1));
+    dLen = cell2mat(syllableInfo(pids(2),4)) - cell2mat(syllableInfo(pids(2),1));
+    fLen = cell2mat(syllableInfo(pids(3),4)) - cell2mat(syllableInfo(pids(3),1));
+    
+    
+    
+    if strcmp(stim_date, '210610')
+        bl_1 = 1;
+        pert = 476;
+        hab = 2496;
+        bl_2 = 2896;
+    elseif strcmp(stim_date, '210613')
+        bl_1 = 1;
+        pert = 476;
+        hab = 2488;
+        bl_2 = 2888;
+    end
+    
+    target_pert = pert_types{example_info(ex_idx,3)};
+    target_syll = target_pert(2);
+    target_type = target_pert(4);
+    target_class = target_pert(5);
+    
+    switch target_syll
+        case 'b'
+            pert_time_2 = pert_time_2(1);
+            syll_lens_1 = bLen;
+        case 'd' 
+            pert_time_2 = pert_time_2(2);
+            syll_lens_1 = dLen;
+        case 'f'
+            pert_time_2 = pert_time_2(3);
+            syll_lens_1 = fLen;
+    end
+    
+    
+    
+    
+    relevant_syll = syllableInfo; % this is perturb phase
+    motif_starts = cell2mat(motifInfo(:,1));
+    syllTypes = syllableInfo(:,3);
+    mTypes = motifInfo(:,2);
+    unique_types = unique(mTypes);
+    
+    % we now need to build structure with start times of each relevant
+    % trial
+    
+    % what are relevant trials?
+    relevant_types = {'@@@@@'}; % this just needs to be something absent from strings so that first cell is empty
+    for uIdx = 1:length(unique_types)
+        thisU = unique_types{uIdx};
+        if contains(thisU, target_pert) % so is okay if target type doesn't exist...
+            relevant_types{end+1} = thisU;
+        end
+    end
+    
+    
+    start_times = cell(length(relevant_types),1);
+    
+    % first prehab
+    for mIdx = bl_1:pert-1
+        prevStart = motif_starts(19*(floor((mIdx-1)/19)) + 1);
+        
+        thisStart = motif_starts(mIdx);
+        thisType = mTypes{mIdx};
+        typeNum = [];
+        relIdx = 1;
+        while isempty(typeNum) && (relIdx <= length(relevant_types))
+            if contains(thisType, relevant_types{relIdx})
+                typeNum = relIdx;
+            end
+            relIdx = relIdx+1;
+        end
+        if ~isempty(typeNum)
+            start_times{typeNum} = [start_times{typeNum}; thisStart];
+            if typeNum == 2 % only add std once
+                start_times{1} = [start_times{1}; prevStart]; % add previous standard motif
+            end
+        end
+    end
+    
+    
+    % pert phase
+    for mIdx = pert:hab-1
+        if mIdx-pert > 1
+            prevStart = motif_starts(mIdx-1);
+        end
+        thisStart = motif_starts(mIdx);
+        thisType = mTypes{mIdx};
+        typeNum = [];
+        relIdx = 1;
+        while isempty(typeNum) && (relIdx <= length(relevant_types))
+            if contains(thisType, relevant_types{relIdx})
+                typeNum = relIdx;
+            end
+            relIdx = relIdx+1;
+        end
+        if ~isempty(typeNum)
+            start_times{typeNum} = [start_times{typeNum}; thisStart];
+            start_times{1} = [start_times{1}; prevStart]; % add previous standard motif
+        end
+    end
+    
+    
+    % second prehab
+    for mIdx = bl_2:length(motifTimingData(:,1))
+        prevStart = motif_starts(19*(floor((mIdx-1)/19)) + 1);
+        
+        thisStart = motif_starts(mIdx);
+        thisType = mTypes{mIdx};
+        typeNum = [];
+        relIdx = 1;
+        while isempty(typeNum) && (relIdx <= length(relevant_types))
+            if contains(thisType, relevant_types{relIdx})
+                typeNum = relIdx;
+            end
+            relIdx = relIdx+1;
+        end
+        if ~isempty(typeNum)
+            start_times{typeNum} = [start_times{typeNum}; thisStart];
+            if typeNum == 2 % only add std once
+                start_times{1} = [start_times{1}; prevStart]; % add previous standard motif
+            end
+        end
+    end
+    
+    this_id = (example_info(ex_idx,2)); % not clearly grabbing correct unit? or else wrong pert
+    theseSpikes = sp.st(sp.clu==goodUnits(this_id));
+    perturb_response_plot_psth_vF(theseSpikes, start_times, this_id, timeWindow_1, micData, target_type, pert_time_2, syll_lens_1, 1, ex_idx)
+    
+    h5 = gcf;
+    print(['figure_pieces/fig4_oddballExample' '_' num2str(ex_idx)], '-dsvg', '-r300')
+    saveas(h5, ['figure_pieces/fig4_oddballExample'  '_' num2str(ex_idx) '.fig'])
+end
+
+disp('i think im done')

NCM-errors/NCM-errors-analysis-code/ncm_errors_fig5_error_specific_coding.m

@@ -0,0 +1,396 @@
+% 1. load multibird summary data
+% 2. plot distributions of # of error responses per neuron
+% 3. plot likelihood of sensitivity for syll types and pert types
+% 4. plot indices (effect mag) for different types (syll and pert)
+
+% NB: missing stat tests, which would be nice because some trends seem
+% significant - probably just binomial?
+% could possibly show that mag distributions are not significantly
+% different
+
+clear
+close all
+error_detection_pipeline = 'D:\error_expmts\multibird_summary_ncm_errors.mat';
+
+rng(23) % important for simulating number of pert responses per neuron -NOT USED
+%% setup 
+
+load(error_detection_pipeline)
+
+% these permutations will be useful
+% sort by ptype and then by intra v whole syll pert
+pType_perm = [1 7 13 2 8 14 3 9 15 4 10 16 5 11 17 6 12 18];
+
+% sort by ptype and then by H v N v S
+pKind_perm = [1 7 13 4 10 16 2 8 14 5 11 17 3 9 15 6 12 18];
+
+% get magnitude for each pert, learned and oddball
+a1_total = mean(big_resp_mat(:, 1:25),2);
+odd_total = mean(big_resp_mat(:,26:50),2);
+a2_total = mean(big_resp_mat(:,51:75),2);
+
+err_idx = a2_total - ((a1_total +  odd_total)/2);
+odd_idx = odd_total - ((a1_total + a2_total)/2);
+
+ei_sig = err_idx(big_fli_ids);
+oi_sig = odd_idx(big_rdi_ids);
+
+%% get sim numbers - what if error responses uniformly distributed? - NOT USED!
+nrns_per_bird = [173; 106; 120; 87; 85];
+bStart = [1; 174; 280; 400; 487];
+bEnd = [173; 279; 399; 486; 571];
+
+total_nrns = sum(nrns_per_bird);
+
+brdi_blank = rand(100,size(big_rare_dev_iti,1), size(big_rare_dev_iti,2));
+bfli_blank = rand(100,size(big_fam_late_iti,1), size(big_fam_late_iti,2));
+
+brdi_sim = brdi_blank;
+bfli_sim = bfli_blank;
+for bird_idx = 1:5
+    these_ids = bStart(bird_idx):bEnd(bird_idx);
+    brdi_sim(:,:,these_ids) = brdi_blank(:,:,these_ids) < ...
+        sum(sum(big_rare_dev_iti(:,these_ids)) / ...
+        numel(big_rare_dev_iti(:,these_ids)));
+    
+    bfli_sim(:,:,these_ids) = bfli_blank(:,:,these_ids) < ...
+        sum(sum(big_fam_late_iti(:,these_ids)) / ...
+        numel(big_fam_late_iti(:,these_ids)));
+end
+
+% brdi_sim = brdi_blank < (sum(sum(big_rare_dev_iti))/numel(big_rare_dev_iti));
+% bfli_sim = bfli_blank < (sum(sum(big_fam_late_iti))/numel(big_fam_late_iti));
+
+sim_edges = -0.5:1:18.5;
+
+rpn_rd = zeros(100,19);
+rpn_fl = zeros(100,19);
+for sim_idx = 1:100
+    this_brdi_sim = squeeze(brdi_sim(sim_idx,:,:));
+    this_bfli_sim = squeeze(bfli_sim(sim_idx,:,:));
+    
+    this_rpn_rd = sum(this_brdi_sim,1);
+    this_rpn_fl = sum(this_bfli_sim,1);
+    
+    this_hc_rd = histcounts(this_rpn_rd,sim_edges);
+    this_hc_fl = histcounts(this_rpn_fl,sim_edges);
+    
+    rpn_rd(sim_idx,:) = this_hc_rd;
+    rpn_fl(sim_idx,:) = this_hc_fl;
+end
+%% distributions of # of learned error responses 
+
+num_neurons = length(big_resp_mat(:, 1))/18;
+
+
+resp_per_nrn = [];
+sig_nrn_ids = [];
+for nIdx = 1:num_neurons
+    start_idx = ((nIdx-1)*18)+1;
+    end_idx = start_idx + 17;
+    
+    num_resp = length(intersect(find(big_fli_ids >= start_idx), find(big_fli_ids <= end_idx)));
+    
+    sig_nrn_bool = ~isempty(intersect(find(big_fli_ids >= start_idx), find(big_fli_ids <= end_idx)));
+    if sig_nrn_bool
+        sig_nrn_ids = [sig_nrn_ids; nIdx];
+    end
+    
+    resp_per_nrn = [resp_per_nrn; num_resp];
+end
+
+% sim results
+mean_sim = mean(rpn_fl,1); % average across sims
+sd_sim = std(rpn_fl,0,1); % stds
+
+figure(81);
+histogram(resp_per_nrn(sig_nrn_ids))
+% hold on;
+% plot(mean_sim(2:19));
+% hold off;
+title('Learned Error Responses per Neuron')
+
+h5 = gcf;
+print('figure_pieces/fig5_errorsPerNeuron', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig5_errorsPerNeuron' '.fig'])
+
+%% distributions of # of oddball error responses 
+
+num_neurons = length(big_resp_mat(:, 1))/18;
+
+resp_per_nrn = [];
+sig_nrn_ids = [];
+for nIdx = 1:num_neurons
+    start_idx = ((nIdx-1)*18)+1;
+    end_idx = start_idx + 17;
+    
+    num_resp = length(intersect(find(big_rdi_ids >= start_idx), find(big_rdi_ids <= end_idx)));
+    
+    sig_nrn_bool = ~isempty(intersect(find(big_rdi_ids >= start_idx), find(big_rdi_ids <= end_idx)));
+    if sig_nrn_bool
+        sig_nrn_ids = [sig_nrn_ids; nIdx];
+    end
+    
+    resp_per_nrn = [resp_per_nrn; num_resp];
+end
+
+% sim results
+mean_sim = mean(rpn_rd,1); % average across sims
+sd_sim = std(rpn_rd,0,1); % stds
+
+figure(82);
+histogram(resp_per_nrn(sig_nrn_ids))
+title('Oddball Error Responses per Neuron')
+% hold on;
+% plot(mean_sim(2:19));
+% hold off;
+
+h5 = gcf;
+print('figure_pieces/fig5_oddErrorsPerNeuron', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig5_oddErrorsPerNeuron' '.fig'])
+
+
+%% distributions of sensitivity by syllable type of learned error responses 
+
+bfli_pTypes = sum(big_fam_late_iti,2);
+bfli_neurons = big_fam_late_iti(:,familiar_late_iti);
+
+pfli_bars = [sum(bfli_pTypes(1:6)), sum(bfli_pTypes(7:12)), sum(bfli_pTypes(13:18))];
+figure(83);
+subplot(311)
+bar(pfli_bars)
+title('Fam. Late - Syllable Type Sorting (b, d, f)')
+
+subplot(312)
+bar([sum(bfli_pTypes(pType_perm(1:9))), sum(bfli_pTypes(pType_perm(10:18)))]);
+title('Fam. Late - Pert. Length Sorting (intra v whole)')
+
+subplot(313)
+bar([sum(bfli_pTypes(pKind_perm(1:6))), sum(bfli_pTypes(pKind_perm(7:12))),...
+    sum(bfli_pTypes(pKind_perm(13:18)))])
+title('Fam. Late - Pert. Type Sorting (HS, N, Sil)')
+
+h5 = gcf;
+print('figure_pieces/fig5_flErrorTypes', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig5_flErrorTypes' '.fig'])
+
+
+%% distributions of sensitivity by syllable type of oddball error responses 
+
+brdi_pTypes = sum(big_rare_dev_iti,2);
+brdi_neurons = big_rare_dev_iti(:,rare_dev_iti);
+
+prdi_bars = ([sum(brdi_pTypes(1:6)), sum(brdi_pTypes(7:12)), sum(brdi_pTypes(13:18))]);
+figure(84);
+subplot(311)
+bar(prdi_bars)
+title('Rare Dev. - Syllable Type Sorting (b, d, f)')
+
+subplot(312)
+bar([sum(brdi_pTypes(pType_perm(1:9))), sum(brdi_pTypes(pType_perm(10:18)))]);
+title('Rare Dev. - Pert. Length Sorting (intra v whole)')
+
+subplot(313)
+bar([sum(brdi_pTypes(pKind_perm(1:6))), sum(brdi_pTypes(pKind_perm(7:12))),...
+    sum(brdi_pTypes(pKind_perm(13:18)))])
+title('Rare Dev. - Pert. Type Sorting (HS, N, Sil)')
+
+h5 = gcf;
+print('figure_pieces/fig5_rdErrorTypes', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig5_rdErrorTypes' '.fig'])
+
+
+%% get effect magnitudes for learned and oddball error responses
+% first get effect magnitudes by perturbation
+
+ei_mags_by_p = {};
+oi_mags_by_p = {};
+
+fl_idx = 0;
+rd_idx = 0;
+for p_idx = 1:18
+    this_row_fl = big_fam_late_iti(p_idx,:);
+    this_row_rd = big_rare_dev_iti(p_idx,:);
+    
+    num_fls = length(find(this_row_fl==1));
+    fl_new = fl_idx + num_fls;
+    if num_fls > 0
+        ei_mags_by_p{p_idx} = ei_sig(fl_idx+1:fl_new);
+    end
+    fl_idx = fl_new;
+    
+    num_rds = length(find(this_row_rd==1));
+    rd_new = rd_idx + num_rds;
+    if num_rds > 0
+        oi_mags_by_p{p_idx} = oi_sig(rd_idx+1:rd_new);
+    end
+    rd_idx = rd_new;
+end
+    
+%% plot learned error effect magnitudes (syllable)
+
+% syll type
+syll_rows = [1:6 ; 7:12 ; 13:18];
+fl_mag_syll_dist = cell(3,1);
+for syll_idx = 1:3
+    for row_idx = syll_rows(syll_idx,:)
+        fl_mag_syll_dist{syll_idx} =...
+            [fl_mag_syll_dist{syll_idx}; ei_mags_by_p{row_idx}];
+    end
+end
+
+
+edges = [0:0.05:1];
+for syll_idx = 1:3
+    figure(86);
+    subplot(3,1,syll_idx)
+    histogram(fl_mag_syll_dist{syll_idx},edges)
+    if syll_idx == 1
+        title('Fam. Late - Syllable Type Sorting (b, d, f)')
+    end
+end
+
+h5 = gcf;
+print('figure_pieces/fig5_flErrorMags', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig5_flErrorMags' '.fig'])
+
+
+%% plot learned error effect magnitudes (error type (intra v whole))
+
+% pert type
+pType_rows = [1 7 13 2 8 14 3 9 15; 4 10 16 5 11 17 6 12 18];
+fl_mag_pType_dist = cell(2,1);
+for pType_idx = 1:2
+    for row_idx = pType_rows(pType_idx,:)
+        fl_mag_pType_dist{pType_idx} =...
+            [fl_mag_pType_dist{pType_idx}; ei_mags_by_p{row_idx}];
+    end
+end
+
+
+edges = [0:0.05:1];
+for pType_idx = 1:2
+    figure(87);
+    subplot(2,1,pType_idx)
+    histogram(fl_mag_pType_dist{pType_idx},edges)
+    if pType_idx == 1
+        title('Fam. Late - Pert. Type Sorting (intra v whole)')
+    end
+end
+
+h5 = gcf;
+print('figure_pieces/fig5_flErrorMags_type', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig5_flErrorMags_type' '.fig'])
+
+
+%% plot learned error effect magnitudes (pert kind - harm stack v noise v silence)
+
+% sort by ptype and then by H v N v S
+pKind_rows = [1 7 13 4 10 16; 2 8 14 5 11 17; 3 9 15 6 12 18];
+fl_mag_pKind_dist = cell(3,1);
+for pKind_idx = 1:3
+    for row_idx = pKind_rows(pKind_idx,:)
+        fl_mag_pKind_dist{pKind_idx} =...
+            [fl_mag_pKind_dist{pKind_idx}; ei_mags_by_p{row_idx}];
+    end
+end
+
+
+edges = [0:0.05:1];
+for pKind_idx = 1:3
+    figure(88);
+    subplot(3,1,pKind_idx)
+    histogram(fl_mag_pKind_dist{pKind_idx},edges)
+    if pKind_idx == 1
+        title('Fam. Late - Pert. Kind Sorting (H v N v S)')
+    end
+end
+
+h5 = gcf;
+print('figure_pieces/fig5_flErrorMags_kind', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig5_flErrorMags_kind' '.fig'])
+
+%% plot oddball error effect magnitudes (syllable)
+
+% syll type
+syll_rows = [1:6 ; 7:12 ; 13:18];
+rd_mag_syll_dist = cell(3,1);
+for syll_idx = 1:3
+    for row_idx = syll_rows(syll_idx,:)
+        rd_mag_syll_dist{syll_idx} =...
+            [rd_mag_syll_dist{syll_idx}; oi_mags_by_p{row_idx}];
+    end
+end
+
+
+edges = [0:0.05:1];
+for syll_idx = 1:3
+    figure(89);
+    subplot(3,1,syll_idx)
+    histogram(rd_mag_syll_dist{syll_idx},edges)
+    if syll_idx == 1
+        title('Rare Dev. - Syllable Type Sorting (b, d, f)')
+    end
+end
+
+h5 = gcf;
+print('figure_pieces/fig5_rdErrorMags', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig5_rdErrorMags' '.fig'])
+
+
+%% plot oddball error effect magnitudes (error type (intra v whole))
+
+% pert type
+pType_rows = [1 7 13 2 8 14 3 9 15; 4 10 16 5 11 17 6 12 18];
+rd_mag_pType_dist = cell(2,1);
+for pType_idx = 1:2
+    for row_idx = pType_rows(pType_idx,:)
+        rd_mag_pType_dist{pType_idx} =...
+            [rd_mag_pType_dist{pType_idx}; oi_mags_by_p{row_idx}];
+    end
+end
+
+
+edges = [0:0.05:1];
+for pType_idx = 1:2
+    figure(90);
+    subplot(2,1,pType_idx)
+    histogram(rd_mag_pType_dist{pType_idx},edges)
+    if pType_idx == 1
+        title('Rare Dev. - Pert. Type Sorting (intra v whole)')
+    end
+end
+
+h5 = gcf;
+print('figure_pieces/fig5_rdErrorMags_type', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig5_rdErrorMags_type' '.fig'])
+
+%% plot oddball error effect magnitudes (pert kind - h stck v noise v silence)
+
+% sort by ptype and then by H v N v S
+pKind_rows = [1 7 13 4 10 16; 2 8 14 5 11 17; 3 9 15 6 12 18];
+rd_mag_pKind_dist = cell(3,1);
+for pKind_idx = 1:3
+    for row_idx = pKind_rows(pKind_idx,:)
+        rd_mag_pKind_dist{pKind_idx} =...
+            [rd_mag_pKind_dist{pKind_idx}; oi_mags_by_p{row_idx}];
+    end
+end
+
+
+edges = [0:0.05:1];
+for pKind_idx = 1:3
+    figure(91);
+    subplot(3,1,pKind_idx)
+    histogram(rd_mag_pKind_dist{pKind_idx},edges)
+    if pKind_idx == 1
+        title('Rare Dev. - Pert. Kind Sorting (H v N v S)')
+    end
+end
+
+h5 = gcf;
+print('figure_pieces/fig5_rdErrorMags_kind', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig5_rdErrorMags_kind' '.fig'])
+
+%%
+disp('all done')

NCM-errors/NCM-errors-analysis-code/ncm_errors_fig6_fs_and_rs.m

@@ -0,0 +1,237 @@
+% 1. load waveform data (possibly from separate pipeline)
+% 2. perform 1d separation of fs v rs neurons using end slope
+% 3. distribution of error sensitivity, fs v rs (learned & oddball)
+% 4. oddball index, fs v rs
+% 5. learned error index, fs v rs
+%
+% message is: "this distinction doesn't really matter," which further
+% supports central claim that NCM is so densely interconnected that
+% everything is multiplexed, pointing towards a dense model for error
+% detection, rather than a sparse one
+
+% at the very least, this shows that circuit-breaking approaches used when
+% studying mammalian systems won't really apply to NCM
+
+%% setup
+clear
+close all
+
+the_good_birds = {'D:\error_expmts\pipeline_outputs\pipeline_7635_210723_LH_NCM_g0.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_9007_210722_LH_NCM_g0.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_8975_210610_LH_NCM_g1.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_9047_210722_LH_NCM_g0.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_9049_210724_LH_NCM_g0.mat'
+    };
+
+error_detection_pipeline = 'D:\error_expmts\multibird_summary_ncm_errors.mat';
+
+rng(23) % just affects the sorting of summary fig
+
+%% get waveform data across birds and process it somewhat
+big_end_slopes = [];
+big_wfs = [];
+
+for bird_id = 1:length(the_good_birds)
+    this_bird = the_good_birds{bird_id};
+    load(this_bird)
+    load(['D:/' bird_id '_wfData_' rec_date '.mat'])
+    
+    end_slopes = [];
+    scaled_wfs = [];
+    for wf_idx = 1:size(best_wfs,1) % only for good units - no problem
+        twf = best_wfs(wf_idx,:);
+        [~,this_trough] = min(twf(26:end-25)); % should be in middle...
+        this_trough = this_trough + 25;
+        this_chopped = twf(this_trough-25:this_trough+25);
+        dmWF = this_chopped - mean(this_chopped(1:10));
+        dmWF = dmWF ./ abs(min(dmWF));
+        dwf = diff(dmWF);
+        
+        scaled_wfs = [scaled_wfs; dmWF];
+        
+        this_end = (dmWF(31) - dmWF(26)) / 5; 
+        end_slopes = [end_slopes; this_end];
+    end
+    
+    % then we can look at patterns outside of loop
+    big_wfs = [big_wfs; scaled_wfs];
+    big_end_slopes = [big_end_slopes; end_slopes];
+end
+
+% also load error pipeline output
+load(error_detection_pipeline)
+
+disp('data is loaded!')
+
+%% perform 1d separation
+
+wf_types = kmeans(big_end_slopes,2);
+
+fs_ids = find(wf_types==2); % larger slopes
+rs_ids = find(wf_types==1); % smaller slopes
+
+edges = 0:0.006125:.5;
+figure(12);
+subplot(211)
+histogram(big_end_slopes,edges)
+title('End Slopes of Waveform')
+
+subplot(212)
+histogram(big_end_slopes(rs_ids),edges)
+hold on;
+histogram(big_end_slopes(fs_ids),edges)
+hold off;
+title('K-means Clustering')
+legend('Putative RS Cells', 'Putative FS Cells', 'Location', 'Best')
+
+h5 = gcf;
+print('figure_pieces/fig6_fs_rs_split', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig6_fs_rs_split' '.fig'])
+
+%% processing for error types, etc
+
+% get magnitude for each pert, learned and oddball
+a1_total = mean(big_resp_mat(:, 1:25),2);
+odd_total = mean(big_resp_mat(:,26:50),2);
+a2_total = mean(big_resp_mat(:,51:75),2);
+
+err_idx = a2_total - ((a1_total +  odd_total)/2);
+odd_idx = odd_total - ((a1_total + a2_total)/2);
+
+ei_sig = err_idx(big_fli_ids);
+oi_sig = odd_idx(big_rdi_ids);
+
+num_neurons = length(big_resp_mat(:, 1))/18;
+
+resp_per_nrn_fl = [];
+sig_nrn_ids_fl = [];
+resp_per_nrn_rd = [];
+sig_nrn_ids_rd = [];
+for nIdx = 1:num_neurons
+    start_idx = ((nIdx-1)*18)+1;
+    end_idx = start_idx + 17;
+    
+    % fl processing
+    num_resp_fl = length(intersect(find(big_fli_ids >= start_idx), find(big_fli_ids <= end_idx)));
+    sig_nrn_bool = ~isempty(intersect(find(big_fli_ids >= start_idx), find(big_fli_ids <= end_idx)));
+    if sig_nrn_bool
+        sig_nrn_ids_fl = [sig_nrn_ids_fl; nIdx];
+    end
+    resp_per_nrn_fl = [resp_per_nrn_fl; num_resp_fl];
+    ei2(nIdx) = max(err_idx(start_idx:end_idx));
+    
+    % rd processing
+    num_resp_rd = length(intersect(find(big_rdi_ids >= start_idx), find(big_rdi_ids <= end_idx)));
+    sig_nrn_bool = ~isempty(intersect(find(big_rdi_ids >= start_idx), find(big_rdi_ids <= end_idx)));
+    if sig_nrn_bool
+        sig_nrn_ids_rd = [sig_nrn_ids_rd; nIdx];
+    end
+    resp_per_nrn_rd = [resp_per_nrn_rd; num_resp_rd];
+    oi2(nIdx) = max(odd_idx(start_idx:end_idx));
+    
+end
+
+%% probability of (learned, oddball) responses given (fs, rs)
+
+% what is probability that cell has error given that it is fs v rs?
+% learned errors
+fs_learned_prop = length(intersect(sig_nrn_ids_fl, fs_ids))/length(fs_ids);
+rs_learned_prop = length(intersect(sig_nrn_ids_fl, rs_ids))/length(rs_ids);
+total_learned_prop = length(sig_nrn_ids_fl) / num_neurons;
+
+% oddball errors
+fs_oddball_prop = length(intersect(sig_nrn_ids_rd, fs_ids))/length(fs_ids);
+rs_oddball_prop = length(intersect(sig_nrn_ids_rd, rs_ids))/length(rs_ids);
+total_oddball_prop = length(sig_nrn_ids_rd) / num_neurons;
+
+
+% what is probability that cell is fs v rs given that is has error?
+learned_fs_prop = length(intersect(sig_nrn_ids_fl, fs_ids))/length(sig_nrn_ids_fl);
+learned_rs_prop = length(intersect(sig_nrn_ids_fl, rs_ids))/length(sig_nrn_ids_fl);
+total_rs_prop = length(rs_ids) / num_neurons;
+
+% oddball errors
+oddball_fs_prop = length(intersect(sig_nrn_ids_rd, fs_ids))/length(sig_nrn_ids_rd);
+oddball_rs_prop = length(intersect(sig_nrn_ids_rd, rs_ids))/length(sig_nrn_ids_rd);
+total_fs_prop = length(fs_ids) / num_neurons;
+
+
+
+
+figure(14);
+subplot(2,2,3)
+bar([total_learned_prop fs_learned_prop rs_learned_prop])
+title('Prop. of Neurons with Learned Errors (total, fs, rs)')
+
+subplot(2,2,4)
+bar([total_oddball_prop fs_oddball_prop rs_oddball_prop])
+title('Prop. of Neurons with Oddball Errors (total, fs, rs)')
+
+
+subplot(2,2,2)
+bar([total_fs_prop learned_fs_prop oddball_fs_prop])
+title('Prop. of Neurons that are FS (total, learned, error)')
+
+subplot(2,2,1)
+bar([total_rs_prop learned_rs_prop oddball_rs_prop])
+title('Prop. of Neurons that are RS (total, learned, error)')
+
+h5 = gcf;
+print('figure_pieces/fig6_fs_rs_errors', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig6_fs_rs_errors' '.fig'])
+
+%% response magnitudes, fs v rs
+
+% first get sig ids for fs, rs
+sig_err_rs = intersect(sig_nrn_ids_fl,rs_ids);
+sig_odd_rs = intersect(sig_nrn_ids_rd,rs_ids);
+sig_err_fs = intersect(sig_nrn_ids_fl,fs_ids);
+sig_odd_fs = intersect(sig_nrn_ids_rd,fs_ids);
+
+% learned error indices
+[cts, edges] = histcounts(ei2(rs_ids));
+[cts2, edges2] = histcounts(ei2(fs_ids));
+figure(7);
+subplot(211)
+histogram(ei2(rs_ids), 'BinEdges', edges)
+hold on;
+histogram(ei2(sig_err_rs),'BinEdges',edges)
+hold off;
+title('Learned Error Indices - Max per Neuron - RS')
+
+subplot(212)
+histogram(ei2(fs_ids), 'BinEdges', edges2)
+hold on;
+histogram(ei2(sig_err_fs),'BinEdges',edges2)
+hold off;
+title('Learned Error Indices - Max per Neuron - FS')
+
+h5 = gcf;
+print('figure_pieces/fig6_fs_rs_flMags', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig6_fs_rs_flMags' '.fig'])
+
+
+% oddball error indices
+[cts, edges] = histcounts(oi2(rs_ids));
+[cts2, edges2] = histcounts(oi2(fs_ids));
+figure(8);
+subplot(211)
+histogram(oi2(rs_ids), 'BinEdges', edges)
+hold on;
+histogram(oi2(sig_odd_rs),'BinEdges',edges)
+hold off;
+title('Oddball Error Indices - Max per Neuron - RS')
+
+subplot(212)
+histogram(oi2(fs_ids), 'BinEdges', edges2)
+hold on;
+histogram(oi2(sig_odd_fs),'BinEdges',edges2)
+hold off;
+title('Oddball Error Indices - Max per Neuron - FS')
+
+h5 = gcf;
+print('figure_pieces/fig6_fs_rs_rdMags', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig6_fs_rs_rdMags' '.fig'])
+
+%%
+disp('All done')

NCM-errors/NCM-errors-analysis-code/ncm_errors_fig7_habituation_and_errors.m

@@ -0,0 +1,210 @@
+% 1. load big error pipeline output
+% 2. cartoon of complementary habituation and error emergence curves v actual 
+% 3. Scatter plots of error v habituation idxs, oddball v habituation idxs
+% 4. network-wide habituation during each of 4 phases
+%
+% nb: habituation index can be less than 1 bc big_resp_mat values are
+% between -1 and +1 (since all responses are corrected by iti response rate)
+%
+% point is to illustrate insufficiency of simple matched model of error
+% emergence, while maintaining the possibility that the two processes are
+% related, due to habituation over the course of the experiment (and of
+% course it's also possible that habituation triggers some un-observable
+% long-term process that only completes after habituation phase).
+
+% add 1 more plot - imagesc of activity during the first period, then
+% during the second, third and fourth periods (standards only)
+
+%% setup
+clear
+close all
+
+the_good_birds = {'D:\error_expmts\pipeline_outputs\pipeline_7635_210723_LH_NCM_g0.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_9007_210722_LH_NCM_g0.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_8975_210610_LH_NCM_g1.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_9047_210722_LH_NCM_g0.mat',...
+    'D:\error_expmts\pipeline_outputs\pipeline_9049_210724_LH_NCM_g0.mat'
+    };
+
+error_detection_pipeline = 'D:\error_expmts\multibird_summary_ncm_errors.mat';
+
+rng(23) % just affects the sorting of summary fig
+load(error_detection_pipeline)
+
+%% model cartoon (nb: these don't necessarily reflect real data)
+x = 0:100;
+decay_model = exp(-x/16);
+decay_actual = exp(-x/4);
+growth_model = 1 ./ (1 + exp(-0.2*(x-25)));
+growth_actual = zeros(101,1);
+growth_actual(52:101) = 1;
+
+figure(1); 
+subplot(211)
+plot(x, decay_model./max(decay_model), 'linewidth', 4); 
+hold on; 
+plot(x, growth_model, 'linewidth', 4); hold off;
+xlabel('Time')
+ylabel('Response Strength')
+legend('SSA', 'Error Sensitivity', 'Location', 'Best')
+title('Matched Model of Error Emergence')
+xticks([])
+
+figure(1);
+subplot(212)
+plot(x, decay_actual./max(decay_actual), 'linewidth', 4); 
+hold on; 
+plot(x, growth_actual, 'linewidth', 4); hold off;
+xlabel('Time')
+ylabel('Response Strength')
+legend('SSA', 'Error Sensitivity', 'Location', 'Best')
+title('Experimental Data')
+xticks([])
+
+h5 = gcf;
+print('figure_pieces/fig7_timescale_cartoon', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig7_timescale_cartoon' '.fig'])
+
+%% pre-processing to id significant oddball and error responses
+% get magnitude for each pert, learned and oddball
+a1_start = mean(big_resp_mat(:,1:5),2);
+a1_end = mean(big_resp_mat(:,21:25),2);
+a2_end = mean(big_resp_mat(:,71:75),2);
+
+a1_total = mean(big_resp_mat(:, 1:25),2);
+odd_total = mean(big_resp_mat(:,26:50),2);
+a2_total = mean(big_resp_mat(:,51:75),2);
+
+hab_idx = a1_start - a2_end;
+odd_idx = odd_total - ((a1_total + a2_total)/2);
+err_idx = a2_total - ((a1_total +  odd_total)/2);
+
+ei_sig = err_idx(big_fli_ids);
+oi_sig = odd_idx(big_rdi_ids);
+
+num_neurons = length(big_resp_mat(:, 1))/18;
+
+resp_per_nrn_fl = [];
+sig_nrn_ids_fl = [];
+resp_per_nrn_rd = [];
+sig_nrn_ids_rd = [];
+ei2 = []; oi2 = []; hi2 = [];
+for nIdx = 1:num_neurons
+    start_idx = ((nIdx-1)*18)+1;
+    end_idx = start_idx + 17;
+    
+    % fl processing
+    num_resp_fl = length(intersect(find(big_fli_ids >= start_idx), find(big_fli_ids <= end_idx)));
+    sig_nrn_bool = ~isempty(intersect(find(big_fli_ids >= start_idx), find(big_fli_ids <= end_idx)));
+    if sig_nrn_bool
+        sig_nrn_ids_fl = [sig_nrn_ids_fl; nIdx];
+    end
+    resp_per_nrn_fl = [resp_per_nrn_fl; num_resp_fl];
+    ei2(nIdx) = max(err_idx(start_idx:end_idx));
+    
+    % rd processing
+    num_resp_rd = length(intersect(find(big_rdi_ids >= start_idx), find(big_rdi_ids <= end_idx)));
+    sig_nrn_bool = ~isempty(intersect(find(big_rdi_ids >= start_idx), find(big_rdi_ids <= end_idx)));
+    if sig_nrn_bool
+        sig_nrn_ids_rd = [sig_nrn_ids_rd; nIdx];
+    end
+    resp_per_nrn_rd = [resp_per_nrn_rd; num_resp_rd];
+    oi2(nIdx) = max(odd_idx(start_idx:end_idx)); 
+    
+    hi2(nIdx) = max(hab_idx(start_idx:end_idx));
+end
+
+%% scatter plots of learned and oddball error idx v habituation 
+
+figure(5)
+subplot(2,1,1)
+plot(hi2, ei2, 'k.', 'MarkerSize', 12)
+hold on;
+plot(hi2(sig_nrn_ids_fl), ei2(sig_nrn_ids_fl), 'r.', 'MarkerSize', 12)
+hold off;
+xlabel('Hab. Idx.')
+ylabel('Err. Idx.')
+axis([-0.5 1.5 -0.5 0.75])
+title('Learned Error Idx versus Habituation Index - max per neuron')
+
+subplot(2,1,2)
+plot(hi2, oi2, 'k.', 'MarkerSize', 12)
+hold on;
+plot(hi2(sig_nrn_ids_rd), oi2(sig_nrn_ids_rd), 'r.', 'MarkerSize', 12)
+hold off;
+xlabel('Hab. Idx.')
+ylabel('Odd. Idx.')
+axis([-0.5 1.5 -0.5 0.75])
+title('Oddball Error Idx versus Habituation Index - max per neuron')
+
+
+h5 = gcf;
+print('figure_pieces/fig7_hab_scatters', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig7_hab_scatters' '.fig'])
+%% habituation pre-processing - need to get all standard responses
+% load individual bird
+
+std_ids_bl1 = 1:19:475;
+std_ids_pert = find(strcmp(motifInfo(476:2495,2),'std_A')); %motifInfo always same
+std_ids_hab = 2496:2895;
+std_ids_bl2 = 2896:19:3370;
+
+big_std_bl1 = [];
+big_std_pert = [];
+big_std_hab = [];
+big_std_bl2 = [];
+for bird_id = 1:length(the_good_birds)
+    this_bird = the_good_birds{bird_id};
+    load(this_bird)
+    
+    % are motif_resp elements normalized to the ITI at all? NO (fixable if issue)
+    big_std_bl1 = [big_std_bl1; motif_resp(std_ids_bl1,good_unit_ids)'];
+    big_std_pert = [big_std_pert; motif_resp(std_ids_pert,good_unit_ids)'];
+    big_std_hab = [big_std_hab; motif_resp(std_ids_hab,good_unit_ids)'];
+    big_std_bl2 = [big_std_bl2; motif_resp(std_ids_bl2,good_unit_ids)'];
+end
+
+
+%% habituation across all phases
+
+% for every good neuron across all datasets, we want the response to every
+% single standard motif. we'll then calculate a habituation index for bl1,
+% pert, hab, and bl2 phases by calculating 
+% (mean_init - mean_final) ./ (mean_init + mean_final)
+
+hab_bl1 = (mean(big_std_bl1(:,1:5),2) - mean(big_std_bl1(:,end-4:end),2)) ./ ...
+    (mean(big_std_bl1(:,1:5),2) + mean(big_std_bl1(:,end-4:end),2));
+
+hab_pert = (mean(big_std_pert(:,1:5),2) - mean(big_std_pert(:,end-4:end),2)) ./ ...
+    (mean(big_std_pert(:,1:5),2) + mean(big_std_pert(:,end-4:end),2));
+
+hab_hab = (mean(big_std_hab(:,1:5),2) - mean(big_std_hab(:,end-4:end),2)) ./ ...
+    (mean(big_std_hab(:,1:5),2) + mean(big_std_hab(:,end-4:end),2));
+
+hab_bl2 = (mean(big_std_bl2(:,1:5),2) - mean(big_std_bl2(:,end-4:end),2)) ./ ...
+    (mean(big_std_bl2(:,1:5),2) + mean(big_std_bl2(:,end-4:end),2));
+
+
+figure(245);
+subplot(141)
+histogram(hab_bl1)
+title('Baseline Phase 1')
+
+subplot(142)
+histogram(hab_pert)
+title('Pert. Phase')
+
+subplot(143)
+histogram(hab_hab)
+title('Habituation Phase')
+
+subplot(144)
+histogram(hab_bl2)
+title('Baseline Phase 2')
+
+h5 = gcf;
+print('figure_pieces/fig7_hab_across_phases', '-dsvg', '-r300')
+saveas(h5, ['figure_pieces/fig7_hab_across_phases' '.fig'])
+
+%%
+disp('all done')

NCM-errors/NCM-errors-analysis-code/ncm_errors_prelim_pipeline.m

@@ -0,0 +1,524 @@
+% run from C:\Users\Fee/ Lab\Documents\MATLAB\NCM/ Project\useful/
+% scripts\custom\
+
+addpath(genpath('neuropix_pipeline'))
+
+close all;
+clear;
+bird_id = '9049';
+rec_date = '210724';
+stim_date = '210610';
+dataset_name = '9049_210724_LH_NCM_g0';
+
+bin_name = [dataset_name '_t0.nidq.bin'];
+path = ['D:\\' dataset_name];
+
+DemoReadSGLXData(bin_name, path, rec_date)
+
+%% we don't have imec if we didn't cluster...
+myKsDir = [path '\' dataset_name '_imec0'];
+
+% [spikeTimes, spikeAmps, spikeDepths, spikeSites] = ksDriftmap(myKsDir);
+% disp('done with dir')
+
+thisBinName = [dataset_name '_t0.nidq.bin'];
+thisPath = ['D:\\' dataset_name];
+
+[sp, b, Audio_eventTimes, IMEC_eventTimes] = correct_IMEC_Spiketimes(thisPath, thisBinName);
+
+save([thisPath '\syncData_' rec_date '.mat'], 'sp', 'b', 'Audio_eventTimes', 'IMEC_eventTimes')
+
+disp('done with sync!')
+save('sync_checkpoint.mat')
+
+%%
+spwf = sp;
+cids = sp.cids(find(sp.cgs==2));
+
+best_wfs = zeros(length(cids),82);
+
+myKsDir = [path '\' dataset_name '_imec0'];
+for neuron = 1:length(cids)
+    gwfparams.dataDir = myKsDir;    % KiloSort/Phy output folder
+    apD = dir(fullfile(myKsDir, '*ap*.bin')); % AP band file from spikeGLX specifically
+    gwfparams.fileName = apD(1).name;         % .dat file containing the raw
+    gwfparams.dataType = 'int16';            % Data type of .dat file (this should be BP filtered)
+    gwfparams.nCh = 385;                      % Number of channels that were streamed to disk in .dat file
+    gwfparams.wfWin = [-40 41];              % Number of samples before and after spiketime to include in waveform
+    gwfparams.nWf = 100;                    % Number of waveforms per unit to pull out
+    gwfparams.spikeTimes = ceil(spwf.st(spwf.clu==cids(neuron))*30000); % Vector of cluster spike times (in samples) same length as .spikeClusters
+    gwfparams.spikeClusters = spwf.clu(spwf.clu==cids(neuron));
+    
+    wf = getWaveForms(gwfparams);
+    
+    theseWFs = squeeze(wf.waveFormsMean);
+    wfRange = range(theseWFs,2);
+    [~, maxIdx] = max(wfRange);
+    
+    best_wf = theseWFs(maxIdx,:);
+    
+    best_wfs(neuron,:) = best_wf;
+end
+
+save([thisPath '\wfData_' rec_date '.mat'], 'best_wfs')
+
+disp('done with waveforms')
+
+%% timeline info
+% motif times start here
+[micData, fs] = audioread([path '\micData_' rec_date '.wav']);
+
+load(['workspace_' stim_date '.mat'], 'stim_timeline', 'timeline')
+
+
+save_info = [path '/timingData_' rec_date '.mat'];
+
+%theoretically, this varies... (and it got us into trouble)...
+minThresh = 0.015;
+
+timeOffset = find(micData>minThresh,1);
+
+motifStarts = stim_timeline.motifs{2,1};
+motifStarts = motifStarts + (timeOffset/fs);
+
+x_MotifStarts = motifStarts;
+
+motifStarts_samples = round(motifStarts*fs);
+motif_lengths = diff(motifStarts_samples);
+
+mStarts_2 = zeros(1,length(motifStarts));
+mStarts_2(1) = motifStarts_samples(1);
+
+precess = [];
+
+for mIdx = 2:length(motifStarts)
+    prev_m_start = mStarts_2(mIdx-1); 
+    prev_m_end = prev_m_start + motif_lengths(mIdx-1);
+    windowTime1 = 0.02;
+    windowTime2 = 0.02;
+    
+    end_check = prev_m_end - round(windowTime1*fs); % try 10ms step back - equal to shortest ISI (in case overshot)
+    fr_samp = [end_check prev_m_end+round(windowTime2*fs)]; % also go forward just in case undershot (can be arbitrarily forward bc only first will be found)
+    next_m_start = find(micData(fr_samp(1):fr_samp(2)) > minThresh,1);
+%     figure(3); plot(micData(fr_samp(1):fr_samp(2))); pause;
+    mStarts_2(mIdx) = end_check + next_m_start;
+    
+    precess = [precess; prev_m_end-mStarts_2(mIdx)];
+end
+    
+mStarts_samples_2 = mStarts_2;
+mStarts_2 = mStarts_samples_2 / fs;
+
+% good for checking!
+for n = 1:length(mStarts_2)
+    figure(1); 
+    if n < length(mStarts_2)
+        plot(micData((mStarts_samples_2(n)):mStarts_samples_2(n+1)))
+    else
+        plot(micData((mStarts_samples_2(n)):mStarts_samples_2(n)+fs))
+    end
+    title(num2str(n))
+    pause(0.05);%(0.05);
+end
+
+% ALRIGHT, mStarts_samples_2 and mStarts_2 look good! - 
+% maybe now we could do the other analyses
+
+timingData = cell(length(mStarts_2), 2);
+timingData(:,1) = stim_timeline.motifs{1};
+timingData(:,2) = mat2cell(mStarts_2', ones(1,length(mStarts_2)),1);
+
+motifTimingData = timingData; % this is what needs to be saved!
+
+save(save_info, 'motifTimingData');
+
+disp('ALL DONE WITH TIMING DATA!')
+
+%% load checkpoint
+load([path '\timingData_' rec_date '.mat']);
+load([path '\syncData_' rec_date '.mat']);
+load(['workspace_' stim_date '.mat'], 'timeline', 'exp_info');
+
+fs = 40000; % update fs to that of microphone
+
+
+
+[micData, ~] = audioread([path '\micData_' rec_date '.wav']);
+
+
+% song 6-10-21
+if strcmp(stim_date, '210610')
+    syllable_lengths = [2397; 5777; 3735; 4824; 5491] / 40000;
+elseif strcmp(stim_date, '210613')
+    syllable_lengths = [3954; 3482; 6598; 6327; 4992] / 40000;
+end
+
+syllable_offsets = [0; syllable_lengths(1)+.02; sum(syllable_lengths(1:2))+.04; ...
+    sum(syllable_lengths(1:3))+.06; sum(syllable_lengths(1:4))+.08];
+
+motif_length = sum(syllable_lengths)+(4*0.02);
+
+
+motifInfo = cell(length(motifTimingData(:,1)),3);
+syllableInfo = cell(5*length(motifTimingData(:,1)),4);
+for mIdx = 1:length(motifTimingData(:,1))
+    motifInfo(mIdx,1) = motifTimingData(mIdx,2);
+    motifInfo(mIdx,2) = motifTimingData(mIdx,1);
+    motifInfo(mIdx,3) = {cell2mat(motifInfo(mIdx,1))+motif_length};
+    
+    mType = motifTimingData{mIdx,1};
+    
+    if strcmp(stim_date, '210610')
+        sTypes = {'std_a', 'std_b', 'std_c', 'std_d', 'std_f'};
+    elseif strcmp(stim_date, '210613')
+        sTypes = {'std_a', 'std_b', 'std_c', 'std_d', 'std_h'};
+    end
+    
+    if ~strcmp(mType, 'std_A')
+        pSyll = mType(2);
+        switch pSyll
+            case 'b'
+                sTypes{2} = mType;
+            case 'c'
+                sTypes{3} = mType;
+            case 'd'
+                sTypes{4} = mType;
+            case 'f'
+                sTypes{5} = mType;
+            case 'h'
+                sTypes{5} = mType;
+        end
+    end
+    
+    for sub_idx = 1:5
+        this_syll_idx = (mIdx*5)-5+sub_idx;
+        
+        syllableInfo(this_syll_idx,1) = ...
+            {cell2mat(motifTimingData(mIdx,2)) + syllable_offsets(sub_idx)};
+        syllableInfo{this_syll_idx,2} = mType;
+        syllableInfo{this_syll_idx,3} = sTypes{sub_idx};
+        syllableInfo(this_syll_idx,4) = ...
+            {cell2mat(syllableInfo(this_syll_idx,1))...
+            + syllable_lengths(sub_idx)};
+    end
+end
+
+save('motif_checkpoint.mat')
+disp('Ready for action!')
+
+%% check unit quality etc before plotting!
+%% unit quality screen - sometimes ksDriftmap runs out of memory??? no idea what to do in that case...
+% i could probably go all the way back to kilosort and disable some
+% channels? ksdriftmap memory requirement depends on number of spikes
+valid_contacts = [1 240];
+
+[spikeTimes, spikeAmps, spikeDepths, spikeSites] = ksDriftmap(myKsDir);
+disp('drift info here and ready')
+
+valid_spikes = find(spikeSites<=240);
+
+spikeSites = spikeSites(valid_spikes);
+spikeDepths = spikeDepths(valid_spikes);
+spikeAmps = spikeAmps(valid_spikes);
+spikeTimes = spikeTimes(valid_spikes);
+
+sp.st = sp.st(valid_spikes);
+sp.spikeTemplates = sp.spikeTemplates(valid_spikes);
+sp.clu = sp.clu(valid_spikes);
+
+%% now, create response matrices
+goodUnits = sp.cids(find(sp.cgs==2));
+good_unit_ids = 1:length(goodUnits);
+
+resp_win = [0 0.2]; % resp win extends 200ms past syllable end. NOT NECESSARILY THE RIGHT THING TO DO
+iti_win = [-0.5 0]; % should always be at least .75 between motifs, but we want to allow up to 250ms for motif response i think...
+
+motif_types = motifInfo(:,2);
+inter_motif_intervals = cell2mat(motifInfo(2:end,1)) - cell2mat(motifInfo(1:end-1,3));
+bout_start_motifs = [1; find(inter_motif_intervals > 0.5)+1]; % first one always new 'bout'
+motif_starts = cell2mat(motifInfo(:,1));
+syllable_types = syllableInfo(:,3);
+unique_sylls = unique(syllable_types);
+
+% should always be the same for stim 6/10
+if strcmp(stim_date, '210610')
+    pert_transition = 2376;
+    pert_end = 12475;
+    post_transition = 14476;
+elseif strcmp(stim_date, '210613')
+    pert_transition = 2376;
+    pert_end = 12435;
+    post_transition = 14436;
+end
+
+% initialize matrix
+all_trial_resp = zeros(length(sTypes), size(motifInfo,1), length(good_unit_ids));
+iti_resp = zeros(length(bout_start_motifs), length(good_unit_ids));
+motif_resp = zeros(size(motifInfo,1), length(good_unit_ids));
+std_motif_sylls = find(strcmp(syllableInfo(:,2), 'std_A'));
+
+pert_types = unique_sylls(1:18);
+
+avg_resp_by_type_expmt = zeros(length(unique_sylls), length(good_unit_ids));
+
+all_resp_by_type_prepert = zeros(length(pert_types), 25, length(good_unit_ids)); 
+all_resp_by_type_pert = zeros(length(pert_types), 25, length(good_unit_ids));
+all_resp_by_type_postpert = zeros(length(pert_types), 25, length(good_unit_ids));
+
+arbt_prepert_iti = zeros(length(pert_types),25,length(good_unit_ids));
+arbt_pert_iti = arbt_prepert_iti;
+arbt_postpert_iti = arbt_pert_iti;
+
+
+for unitIdx = 1:length(good_unit_ids)
+    disp(['resp matrices for unit ' num2str(unitIdx) ' of ' num2str(length(good_unit_ids)) '.'])
+    this_id = good_unit_ids(unitIdx);
+    theseSpikes = sp.st(sp.clu==goodUnits(this_id));
+    
+    % create massive response matrix! 
+    %(syllable archetype x motif num x neuron)
+    for syllIdx = 1:size(syllableInfo,1)
+        this_motif_num = ceil(syllIdx/length(sTypes));
+        this_syll_num = mod(syllIdx, length(sTypes));
+        if (this_syll_num == 0)
+            this_syll_num = 5;
+        end
+        this_win = cell2mat(syllableInfo(syllIdx,[1 4])) + resp_win;
+        syll_spikes = find((theseSpikes > this_win(1)) & (theseSpikes < this_win(2)));
+        resp_fr = length(syll_spikes) / (this_win(2) - this_win(1));
+        
+        all_trial_resp(this_syll_num, this_motif_num, unitIdx) = resp_fr;
+    end
+    
+    % for each inter-bout period (when motifs are separated by at least
+    % 750ms), find the firing rate during that period (within iti_win)
+    for itiIdx = 1:length(bout_start_motifs)
+        this_motif = bout_start_motifs(itiIdx);
+        this_win = motif_starts(this_motif) + iti_win;
+        iti_spikes = find((theseSpikes > this_win(1)) & (theseSpikes < this_win(2)));
+        resp_fr = length(iti_spikes) / (this_win(2) - this_win(1));
+        
+        iti_resp(itiIdx, unitIdx) = resp_fr;
+    end
+    
+    % for each motif, find the response to that motif (again allowing a
+    % resp_win response window)
+    for motifIdx = 1:size(motifInfo,1)
+        this_motif_win = cell2mat(motifInfo(motifIdx, [1 3])) + resp_win;
+        motif_spikes = find((theseSpikes > this_motif_win(1)) & (theseSpikes < this_motif_win(2)));
+        resp_fr = length(motif_spikes) / (this_motif_win(2) - this_motif_win(1));
+        
+        motif_resp(motifIdx, unitIdx) = resp_fr;
+    end
+    
+    %for each syllable type, find average response to that syllable simply
+    %by finding response to that syllable each time it occurs and averaging
+    %the firing rates
+    for typeIdx = 1:length(unique_sylls)
+        this_type = unique_sylls(typeIdx);
+        these_sylls = find(strcmp(syllable_types, this_type));
+        
+        temp_resp = [];
+        
+        for syllIdx = 1:length(these_sylls)
+            this_win = cell2mat(syllableInfo(these_sylls(syllIdx),[1 4])) + resp_win;
+            syll_spikes = find((theseSpikes > this_win(1)) & (theseSpikes < this_win(2)));
+            resp_fr = length(syll_spikes) / (this_win(2) - this_win(1));
+        
+            temp_resp = [temp_resp; resp_fr];
+        end
+        avg_resp_by_type_expmt(typeIdx, unitIdx) = mean(temp_resp);
+    end
+    
+    
+    % for each p type, for each occurrence of that ptype, find the previous
+    % standard trial and normalize syllable response by that. average 
+    % across these, separated by phase (prepert, pert, postpert)
+    % also do the same with iti normalization
+    for ptypeIdx = 1:length(pert_types)
+        this_p_type = pert_types{ptypeIdx};
+        % to id ITI, find motif number of pert, and then find iti
+        % "bout_start_motifs" entry that has preceding iti
+        
+        
+        these_stds = intersect(std_motif_sylls, find(strcmp(syllableInfo(:,3),['std_' this_p_type(2)])));
+        these_perts = find(strcmp(syllableInfo(:,3),this_p_type));
+        these_motif_nums = ceil(these_perts/length(sTypes));
+        key_itis = zeros(length(these_motif_nums),1);
+        for this_motif_idx = 1:length(these_motif_nums)
+            this_this = these_motif_nums(this_motif_idx);
+            this_list = bout_start_motifs - this_this;
+            this_list(this_list>0) = -999999;
+            [~,key_iti] = max(this_list);
+            key_itis(this_motif_idx) = key_iti;
+        end
+        these_itis = iti_resp(key_itis);
+        
+        temp_resp = [];
+        temp_resp_iti = [];
+        for pIdx = 1:length(these_perts)
+            this_p = these_perts(pIdx);
+            this_std_cand = find(these_stds < this_p);
+            this_std = this_std_cand(end);
+            
+            this_p_win = cell2mat(syllableInfo(this_p,[1 4])) + resp_win;
+            syll_spikes = find((theseSpikes > this_p_win(1)) & (theseSpikes < this_p_win(2)));
+            resp_fr_p = length(syll_spikes) / (this_p_win(2) - this_p_win(1));
+            
+            this_bl_win = cell2mat(syllableInfo(this_std,[1 4])) + resp_win;
+            syll_spikes = find((theseSpikes > this_bl_win(1)) & (theseSpikes < this_bl_win(2)));
+            resp_fr_bl = length(syll_spikes) / (this_bl_win(2) - this_bl_win(1));
+            
+            resp_fr_iti = these_itis(pIdx);
+            
+            temp_resp = [temp_resp; resp_fr_p - resp_fr_bl];
+            temp_resp_iti = [temp_resp_iti; resp_fr_p - resp_fr_iti];
+        end
+        
+        this_prepert = find(these_perts < pert_transition);
+        this_pert = find(these_perts >= pert_transition & these_perts < pert_end);
+        this_postpert = find(these_perts >= post_transition);
+        
+        all_resp_by_type_prepert(ptypeIdx, :, unitIdx) = temp_resp(this_prepert);
+        all_resp_by_type_pert(ptypeIdx, :, unitIdx) = temp_resp(this_pert);
+        all_resp_by_type_postpert(ptypeIdx, :, unitIdx) = temp_resp(this_postpert);
+        
+        arbt_prepert_iti(ptypeIdx,:,unitIdx) = temp_resp_iti(this_prepert);
+        arbt_pert_iti(ptypeIdx,:,unitIdx) = temp_resp_iti(this_pert);
+        arbt_postpert_iti(ptypeIdx,:,unitIdx) = temp_resp_iti(this_postpert);
+    end        
+end
+
+
+%% eliminate low fr
+fr_thresh = 1/60; % 1 spike per minute
+
+% exclude low fr - dangerous because neuron may ONLY spike during pert
+low_fr = [];
+
+goodUnits = sp.cids(find(sp.cgs==2)); % 1 for MUs, 2 for isolated
+for unitIdx = 1:length(goodUnits)
+    theseSpikes = sp.st(sp.clu==goodUnits(unitIdx));
+    this_fr = length(theseSpikes) / (length(micData)/fs);
+    if this_fr < fr_thresh
+        low_fr = [low_fr; unitIdx];
+    end
+end
+
+disp(['found ' num2str(length(low_fr)) ' low fr units.'])
+
+%% eliminate units without auditory responses
+% really what we need is average response for each syllable type, then keep
+% doing what we're doing (except instead of 2*sd, use 2.87*sd for
+% bonferroni correction - pvalue of responsiveness is then .05/23)
+
+% technically, to claim that a unit is definitely auditory, we'd want 2.87
+% sds because of bonferroni. But really we just want to eliminate units
+% that are clearly non-auditory, so we'll use a single std. dev
+
+% we just can't make claims about nonexistence of certain units, but we
+% wouldn't be able to do that anyway
+
+mean_iti = mean(iti_resp,1);
+sd_iti = std(iti_resp,1); 
+
+iti_floor = mean_iti - sd_iti; % bonferroni corrected bounds
+iti_ceil = mean_iti + sd_iti;
+
+% just looking for syllable type that drives greatest response
+max_resp = max(avg_resp_by_type_expmt, [], 1);
+min_resp = min(avg_resp_by_type_expmt, [], 1);
+
+non_aud = find(iti_floor < min_resp & iti_ceil > max_resp);
+
+disp(['found  ' num2str(length(non_aud)) ' non-aud'])
+
+%% eliminate units with significant, long-lasting fr changes
+drift_thresh= 10;
+
+high_fr_drift = [];
+
+win_size = 100;
+step_size = 5;
+moving_avg = [];
+
+for this_start = 1:step_size:size(iti_resp,1)-(win_size-1)
+    this_end = this_start + (win_size-1);
+    
+    this_row = mean(iti_resp(this_start:this_end,:),1);
+    moving_avg = [moving_avg; this_row];
+end
+
+mv_avg_denom = moving_avg;
+mv_avg_denom(mv_avg_denom == 0) = 0.1;
+
+for this_col = 1:size(moving_avg,2)
+    if ((max(moving_avg(:,this_col))/min(moving_avg(:,this_col)))/mean(moving_avg(:,this_col))) > drift_thresh
+        high_fr_drift = [high_fr_drift; this_col];
+    end
+end
+
+disp(['found ' num2str(length(high_fr_drift)) ' with fr drift'])
+    
+% 1. use motif_resp
+% 2. take moving average of that (10 trials, step of 1 trial)
+% 3. fr_drift = diff of that moving average
+% 4. eliminate neurons where fr_drift exceeds 5 (maybe 10) hz (experiment)
+% 5. update all_trial_resp, iti_resp and good_unit_ids accordingly
+
+%% eliminate units that drift significantly in space
+sd_win_size = 10*fs; % 10 second window
+sd_step_size = 2*fs; % 1 second 
+spatial_drift_thresh = 30;
+
+high_spatial_drift = [];
+
+% deffo_bad = load('def_bad_210119b.mat');
+tape_length = round(max(round(spikeTimes*fs))+(5*fs), 5, 'significant');
+win_starts = 1:sd_step_size:(tape_length-sd_win_size)+1;
+win_ends = win_starts + sd_win_size - 1;
+depth_traces = zeros(length(goodUnits),length(win_starts));
+density_traces = depth_traces;
+ctr_samples = (sd_win_size/2):sd_step_size:(tape_length-(sd_win_size/2));
+
+depth_means = zeros(length(goodUnits),1);
+depth_sds = depth_means;
+
+for unitIdx = 1:length(goodUnits)
+    theseDepths = spikeDepths(sp.clu==goodUnits(unitIdx));
+    theseSpikes = spikeTimes(sp.clu==goodUnits(unitIdx));
+    
+    depth_means(unitIdx) = mean(theseDepths);
+    depth_sds(unitIdx) = std(theseDepths);
+    
+    for winIdx = 1:length(win_starts)
+        this_win_start = win_starts(winIdx)/fs;
+        this_win_end  = win_ends(winIdx)/fs;
+        these_spike_ids = find(theseSpikes>this_win_start & theseSpikes<this_win_end);
+        these_sds = theseDepths(these_spike_ids);
+        depth_traces(unitIdx,winIdx) = mean(these_sds);
+        density_traces(unitIdx,winIdx) = length(these_spike_ids)/(this_win_end-this_win_start);
+    end
+end
+
+total_drift = range(depth_traces,2);
+high_spatial_drift = find(total_drift>=spatial_drift_thresh);
+
+disp(['found ' num2str(length(high_spatial_drift)) ' unstable units'])
+
+% 1. signal_loc is average depth of unit across 10 sec window with 2 sec
+% steps
+% 2. eliminate neurons where range of signal_loc exceeds 30u (experiment)
+% 3. update all_trial_resp, iti_resp and good_unit_ids accordingly
+
+%% CHECK IN - HOW MANY NEURONS REMAIN? ARE THEY REASONABLE?
+
+bad_units = union(union(low_fr, non_aud), union(high_fr_drift, high_spatial_drift));
+good_unit_ids = setdiff(1:length(goodUnits), bad_units);
+
+disp(['there are ' num2str(length(good_unit_ids)) ' good single units'])
+
+
+%% 
+save(['D:\\pipeline_outputs\pipeline_' dataset_name '.mat'])
+disp('pipeline complete!')