
function filter_somabc_xt(fname,varargin)
% determine which cells are barcoded



if ~exist('fname','var')
    fname='filt_neurons-bc-somas.mat';
end
% parse thresholds
if ~isempty(varargin)
    varargin=reshape(varargin,2,[]);
    for i=1:size(varargin,2)
        switch lower(varargin{1,i})
            case 'complexity'
                complexity_thresh=varargin{2,i};
            case 'score'
                score_thresh=varargin{2,i};
            case 'signal'
                sig_thresh=varargin{2,i};
            case 'classifier'
                Mdl=varargin{2,i};
        end
    end
end


% default thresholds. Change the default score and sig threshold to NOT do
% thresholding using these.
if ~exist('complexity_thresh','var')
    complexity_thresh=-0.9;
end
if ~exist('score_thresh','var')
    score_thresh=0;
end
if ~exist('sig_thresh','var')
    sig_thresh=0;
end
if ~exist('classifier','var')
    Mdl=[];
end


%load dataset
S=load(fname,'filt_neurons');
filt_neurons=S.filt_neurons;

%calculate complexity of barcodes
bc_complexity=zeros(numel(filt_neurons.id),1);
for i=1:numel(filt_neurons.id)
    bc_complexity(i)=log10(prod(lingseqcomplexity(char(48+filt_neurons.soma_bc(i,:)))));
end


score1s=sort(filt_neurons.soma_bc_score,2,'ascend');
sig1m=sort(squeeze(max(filt_neurons.soma_bc_sig,[],2)),2,'ascend');



pass_complexity=bc_complexity>=complexity_thresh;
% high_score=median(filt_neurons.soma_bc_score,2)>=score_thresh; %median score
% high_sig=median(max(filt_neurons.soma_bc_sig,[],2),3)>=sig_thresh; %median signal of max channel
base_num=size(score1s,2);
high_score=score1s(:,min(base_num,3))>=score_thresh; % thresholding using the 3rd lowest score
high_sig=sig1m(:,min(base_num,3))>=sig_thresh; % thresholding using the 3rd lowest intensity value

% check with trained logistic regression model. The model takes 3rd lowest
% score/4096 and 3rd lowest intensity value to predict which cells are good. Can
% be combined with the straight score and sig thresholding, but not sure
% that's helpful in most cases.

if ~isempty(Mdl)
    X=[sig1m(:,3)/4096, ...
        score1s(:,3)];
    pass_mdl=predict(Mdl, X);
else
    pass_mdl=ones(size(high_sig));
end

filt_neurons.is_barcoded=pass_complexity&high_score&high_sig&pass_mdl>0;


save('filt_neurons-bc-somas.mat','filt_neurons','-v7.3');

end

