added paths for chain_summary.m and chain_stats.m, plus more comments

Example/Example_run.m

@@ -14,17 +14,17 @@
 % Input data
 input.fs = 30000; %acquisition rate
 input.ts = 82; %wf time samples
-input.l2_weights = 1500;
-input.threshold = 10;
-input.dim_mask = logical([1,1,1,0,0,0,0,0,0,1]);
-input.dim_mask_physical = logical([1,1,1,0,0,0,0,0,0,0]);
-input.dim_mask_wf = logical([0,0,0,0,0,0,0,0,0,1]);
+input.threshold = 10;    % distance threshold for calling matches as real
+input.l2_weights = 1500; % weight for waveform difference, see comments in Pipeline\EMD_matlab\weighted_gdf_nt
+input.dim_mask = logical([1,1,1,0,0,0,0,0,0,1]); % used to calculate full distance metric
+input.dim_mask_physical = logical([1,1,1,0,0,0,0,0,0,0]); % used to calculate 'position only' part of distance metric
+input.dim_mask_wf = logical([0,0,0,0,0,0,0,0,0,1]); % used to calculate 'waveform only' part of distance metric
 input.chan_pos_name = 'channel_positions.npy';
 input.wf_name = 'ksproc_mean_waveforms.npy';
 input.KSLabel_name = 'cluster_KSLabel.tsv';
 input.validation = 0;
 
-numData = 2; 
+numData = 5; 
 
 
 
@@ -49,13 +49,13 @@
     all_input(id) = load(fullfile(input.input_path,['result',num2str(id),num2str(id+1)], "Input.mat")); 
     all_output(id) = load(fullfile(input.input_path,['result',num2str(id),num2str(id+1)],'Output.mat'));
 end
-[chain_all,z_loc,len] = chain_summary(all_input,all_output,numData);
+[chain_all,z_loc,len] = chain_summary(all_input,all_output,numData,input.input_path);
 
 
 
 %----------Plot chains of interest (waveform, firing rate, original location, drift-corrected location, L2)----------
 full_chain = chain_all(len == numData,:); %find chains with length across all datasets
-[L2_weight,fr_all,fr_change,x_loc_all,z_loc_all] = chain_stats(all_input,all_output,full_chain,numData);
+[L2_weight,fr_all,fr_change,x_loc_all,z_loc_all] = chain_stats(all_input,all_output,full_chain,numData,input.input_path);
 
 numChain = size(full_chain,1);
 ichain = 1; %which chain to plot, please enter a number between 1 and numChain as input, NEED CHANGE  

Pipeline/NT_main.m

@@ -1,7 +1,8 @@
 % This is the main function of the algorithm
 % Input: kilosort cluster label, channel map, mean waveforms, 
 % Output: Unit match assignment
-% For more comparisons, users need to write their own loops
+% For more comparisons, users need to write their own loops; see
+% Example_run.m
 function NT_main(input,chan_pos,mwf1,mwf2)
 
 % Estimate location 

Pipeline/chain/chain_stats.m

@@ -1,5 +1,5 @@
 % calculate L2, FR, and locations of all chains
-function [L2_weight,fr_all,fr_change,x_loc,z_loc] = chain_stats(all_input,all_output,full_chain,numData)
+function [L2_weight,fr_all,fr_change,x_loc,z_loc] = chain_stats(all_input,all_output,full_chain,numData,output_path)
 
 for ichain = 1:size(full_chain,1)
     for id = 1:numData-1
@@ -68,5 +68,5 @@
     end
 end
 
-save('C:\Users\labadmin\Desktop\Neuron Tracking Pipeline\User version\chain_stats.mat','full_chain','L2_weight','fr_all','fr_change','z_loc','x_loc')
+save(fullfile(output_path,'chain_stats.mat'),'full_chain','L2_weight','fr_all','fr_change','z_loc','x_loc')
 end

Pipeline/chain_summary.m

@@ -1,7 +1,7 @@
 % This is the second function of the algorithm
 % Input: Unit assignment from multiple datasets 
 % Output: Summary of chains within distance threshold
-function [chain_all,z_loc,len] = chain_summary(all_input,all_output,numData)
+function [chain_all,z_loc,len] = chain_summary(all_input,all_output,numData,output_path)
 
 % Summarize chains
 numChain = 200; %set to a number large number
@@ -51,5 +51,5 @@
     len = sum(chain_all ~= 0,2); %chain lengths
 end
 
-save('C:\Users\labadmin\Desktop\Neuron Tracking Pipeline\User version\chain_summary.mat', "chain_all", "z_loc", "len");
+save(fullfile(output_path,'chain_summary.mat'), "chain_all", "z_loc", "len");
 end

Pipeline/create_input/wf_metric/wave_metrics.m

@@ -1,4 +1,8 @@
 function meas_out = wave_metrics(mw, chan_pos, input)
+% Calculate an array of simple waveform metrics from the mean waveforms for
+% each unit.
+% 1D waveform characteristics based on allen ecephys pipeline
+% Position fitting based on Boussard 2D fit to 1/R
 
 [nUnit, nChan, nt] = size(mw);
 pp_all = squeeze(max(mw,[],3)-min(mw,[],3));
@@ -106,20 +110,21 @@
     % 2D waveform metrics    
     pp_unit = squeeze(pp_all(i,:))';
 
-% Julien Boussard - style fit of peak-to-peak voltage vs position
-% if background sub pp_unit > 60 uV, attempt a fit of the
-% background subtracted pp_all
-if max(squeeze(pp_all(i,:))) > 60
-    fitvals = fit_loc(i, pp_all, chan_pos);
-    fitX = fitvals(1);
-    fitZ = fitvals(2);
-    fitY = fitvals(3);
-else
-    fitX = chan_pos(pk_chan(i),1);
-    fitZ = chan_pos(pk_chan(i),2);
-    fitY = -1;      % a marker for no fit
-end
+    % Julien Boussard - style fit of peak-to-peak voltage vs position
+    % if background sub pp_unit > 60 uV, attempt a fit of the
+    % background subtracted pp_all
+    if max(squeeze(pp_all(i,:))) > 60
+        fitvals = fit_loc(i, pp_all, chan_pos);
+        fitX = fitvals(1);
+        fitZ = fitvals(2);
+        fitY = fitvals(3);
+    else
+        fitX = chan_pos(pk_chan(i),1);
+        fitZ = chan_pos(pk_chan(i),2);
+        fitY = -1;      % a marker for no fit
+    end
 
+    % calculate spread of waveforms in z
     sp_thresh = 0.2*max(pp_unit);
     chan_above = pp_unit > sp_thresh;
     zmax = max(chan_pos(chan_above,2));