
function conn=match_bc2neurons_sindbis_mr1(filt_neurons,bc,varargin)
% match bc barcodes to filt_neurons.
% Optional inputs:
%           mismatch: number of mismatches allowed, default 1.
%           filterslice: slices that will be filtered out, relative to soma
%               slice number. For example, [-1 0] means filtering out the
%               previous slice and the soma slice. Default is [] (no
%               filtering).
%           allowduplicates: whether to allow bc rolonies to match to
%               multiple cells. Default is 0. This affect filt_mat, not mat.
%           output: filename to which conn is saved. should be a .mat
%               filename-compatible string/char array. Default is conn.mat.
% Output conn with the following fields:
%           neuron_uid: corresponds to filt_neurons.uid and rows of
%           matrices.
%           bc_uid: corresponds to bc.uid and columns of matrices.
%           match: number of matches for valid matching barcodes. Equals
%           barcode length - NM tag in bowtie output. This include matching
%           barcodes that are non-best matches, but withint the mismatch
%           tolerance.
%           mat: sparse logical matrix that indicates which neuron each
%           barcode matches to. Allows a barcode to match to multiple if
%           they are equally best matches.
%           filt_mat: sparse logical matrix that indicates barcode matching
%           after filtering by slice numbers and by duplicates. If
%           allowduplicates == 0 and filterslice == [], then this is the
%           same as mat.

%
%



% conn.neuron_uid;
% conn.bc_uid;
% conn.mat;
% conn.filt_slice_dist;

% check if uid is present for both neurons and bc
if ~isfield(filt_neurons,'uid')||~isfield(bc,'uid')
    error('uid is absent for filt_neurons or bc. Add uid first before matching barcodes.\n');
end
%%
% default thresholds
mismatch_tolerance=1;
filterslice=[];
allow_duplicates=0;
fname='conn.mat';
rad=300; % in um
hough_dist_thresh=70; % in um
pixel_size=0.65; % um/pxl
houghthresh=0.3; % relative to max H in hough transform
nhoodsize=[3266,101]; % first number is in um, second in deg. This is converted to [rescaled pixel, deg]
fillgap=rad/1.5; % in um, need to convert to pxls
minlength=fillgap*2; % in um, need to convert to pxls
coor='CCF';
numpeaks=20;
dens_disk_rad=65; % in um
uniq_slice=unique(filt_neurons.slice);
suspicious_slides=dictionary(uniq_slice,{uniq_slice}); % by default, only remove center slice.
stitching_tolerance_pxl= 0; % in pixels, default is not filtering for this


%%
% parse thresholds
if ~isempty(varargin)
    varargin=reshape(varargin,2,[]);
    for i=1:size(varargin,2)
        switch lower(varargin{1,i})
            case 'mismatch'
                mismatch_tolerance=varargin{2,i};
                assert(mismatch_tolerance==round(mismatch_tolerance)&&mismatch_tolerance>=0, ...
                    'mismatch should be non-negative integer. Abort.\n');
            case 'filterslice'
                filterslice=varargin{2,i};
                assert(sum(round(filterslice)~=filterslice)==0, ...
                    'filterslice should be integer array or scalar. Abort.\n');
            case 'allowduplicates'
                allow_duplicates=varargin{2,i};
                assert(isnumeric(allow_duplicates)||islogical(allow_duplicates), ...
                    'allow_duplicates should be 0 or 1. Abort.\n');
            case 'output'
                fname=varargin{2,i};
            case 'sigma'
                rad=varargin{2,i};
            case 'dist2linethreshold'
                hough_dist_thresh=varargin{2,i};
            case 'pixel_size'
                pixel_size=varargin{2,i};
            case 'coordinates'
                coor=varargin{2,i};
            case 'houghthreshold'
                houghthresh=varargin{2,i};
            case 'nhoodsize'
                nhoodsize=varargin{2,i};
            case 'fillgap'
                fillgap=varargin{2,i};
            case 'minlength'
                minlength=varargin{2,i};
            case 'numpeaks'
                numpeaks=varargin{2,i};
            case 'diskradius'
                dens_disk_rad=varargin{2,i};
            case 'type'
                experiment_type=varargin{2,i};
            case 'contaminatedslices'
                suspicious_slides=varargin{2,i};% this needs to be a dictionary
            case 'stitchingtolerance'
                stitching_tolerance_pxl=varargin{2,i};

        end
    end
end

%%

rescale=17/pixel_size; % this is scaled to ~30 um per pixel.
rad_pxl=rad/pixel_size/rescale; % rad in rescaled pxls
h=fspecial('gaussian',2*ceil(2*rad_pxl)+1,rad_pxl); % make the gaussian filter used in imgaussfilt.
conf_thresh=max(h,[],'all');%i.e., the intensity expected from a single rol at the same location.
fillgap_pxl=fillgap/pixel_size/rescale; % in rescaled pxls
minlength_pxl=minlength/pixel_size/rescale; % in rescaled pxls
nhoodsize_pxl=[nhoodsize(1)/pixel_size/rescale, nhoodsize(2)]; % in rescaled pxls and deg
nhoodsize_pxl(1)=min(round(nhoodsize_pxl(1)/2)*2+1,round(sqrt(sum(max(filt_neurons.pos).^2))*pixel_size/rescale/5)*2+1);% in case the slice is very small to begin with
dens_disk_rad_pxl=dens_disk_rad/pixel_size/rescale; % in rescaled pxls

prc_thresh=95; % filter out the top pixels when figuring out the background level variations in slices. Currently not tunable


%% make unified soma sequences
% only match to neurons that are is_barcoded
% if soma_bc_hd exists (not all 5s), substitute soma_bc with soma_bc_hd
% then convert sequences to 'GTAC'

somabc=filt_neurons.soma_bc;
hd=filt_neurons.soma_bc_hd(:,1)~=5; % 5 is N which is used when soma_bc_hd is not called.
somabc(hd,:)=filt_neurons.soma_bc_hd(hd,:); %use hd if available

dic='GTACN';
somabcC=dic(somabc);
bcC=dic(bc.seq);

%% write the sequences to fasta files and run bowtie2. Need wsl and bowtie2 toolbox extension installed
% write fasta files
tic
fasta_folder='fasta';
mkdir(fasta_folder);
myfastawrite(fullfile(fasta_folder,'neurons.fa'), ...
    filt_neurons.uid(filt_neurons.is_barcoded), ...
    somabcC(filt_neurons.is_barcoded,:));
myfastawrite(fullfile(fasta_folder,'bc.fa'), ...
    bc.uid, ...
    bcC);
t1=toc;
fprintf('Writing fasta files took %.2g secs.\n',t1);


% make bowtie index
tic
ref_folder='bowtie_ref';
mkdir(ref_folder);
bowtie2build(fullfile(fasta_folder,'neurons.fa'),fullfile(ref_folder,'ref'));
t2=toc;
fprintf('Building bowtie indices took %.2g secs.\n',t2);

% run bowtie

numcores=feature('numcores');
[~,mem]=memory;
maxcore=min(round(numcores*2), floor(mem.PhysicalMemory.Total/2000000000));

alignOptions=Bowtie2AlignOptions('AlignReverseComplementStrand',0, ...
    'NoGapPosition',1, ...
    'ReadGapCosts',[10 3], ... % not alowing gaps
    'RefGapCosts',[10,3], ... % not alowing gaps
    'NumAlignments','All', ...
    'NumThreads',maxcore, ...
    'NumReseeding',3, ...
    'MismatchPenalty',[4 2], ... % this may need to change if more tolernace is needed than 2
    'SeedLength',8, ...
    'ExcludeUnaligned',1, ...
    'NumSeedMismatches',0);
tic

output_dir='bowtie_output';
mkdir(output_dir);
bowtie2(fullfile(ref_folder,'ref'), ...
    fullfile(fasta_folder,'bc.fa'),'', ...
    fullfile(output_dir,'results_bowtie.txt'), ...
    alignOptions);
t4=toc;
fprintf('bowtie matching took %.2g secs.\n',t4);
% read bowtie results
fprintf('Reading bowtie results ...');
tic
%alignment=samread(fullfile(output_dir,'results_bowtie.txt'));
alignment=myread_bowtie_block(fullfile(output_dir,'results_bowtie.txt'));
t5=toc;
fprintf('Done, took %.2g secs.\n', t5);

%% filter by mismatch tolerance
fprintf('Filtering out barcodes with more than %u mismatches ...',mismatch_tolerance);




% convert NM tags, not needed when using myread_bowtie
%tags={alignment.Tags};
%nm=cellfun(@(x) x.NM,tags); % number of edits
%xm=cellfun(@(x) x.XM,tags); % number of substitutions. With large gap penalty, this hsould be the same as nm

% apply mismatch threshold
%keep=nm<=mismatch_tolerance; % version for samread
keep=[alignment.NM]<=mismatch_tolerance; % version for myread_bowtie

alignment=alignment(keep);






fprintf('Done.\n');

%% build connectivity matrix without other filtering (sparse logical)

% 
% tags={alignment.Tags};
% nm1=cellfun(@(x) x.NM,tags); %
% match=size(filt_neurons.soma_bc,2)-nm1; %version for samread
match=size(filt_neurons.soma_bc,2)-[alignment.NM]; % version for myread_bowtie


[~,bc_idx]=ismember({alignment.QueryName},bc.uid);
[~,neuron_idx]=ismember({alignment.ReferenceName},filt_neurons.uid);
% dont' do any other filtering, but only keep the best (possibly tied)
% matches

conn=struct;
conn.neuron_uid=filt_neurons.uid;
conn.bc_uid=bc.uid;

conn.match=sparse(neuron_idx, ...
    bc_idx, ...
    double(match), ...
    numel(filt_neurons.uid), ...
    numel(bc.uid));


% remove non-best matches
fprintf('Filtering out non-best matches ...');
% tic
is_primary=[alignment(:).Flag]'==0;
keep_2nd=ones(numel(is_primary),1,'logical');
primary_match=match(is_primary);
secondary_match=match(~is_primary);
primary_QueryName={alignment(is_primary).QueryName};
[~,I]=ismember({alignment(~is_primary).QueryName},primary_QueryName); %I is the index of primary query  for every secondary query.
keep_2nd(~is_primary)=secondary_match==primary_match(I); % for non-primary matches, only keep when match scores are the same as primary match scores

%make indices
mat_row=neuron_idx(keep_2nd);
mat_col=bc_idx(keep_2nd);

%check if filtering overlapping rolonies
if stitching_tolerance_pxl>0
    [mat_row,mat_col]=filter_overlapping_bc(mat_row,mat_col, ...
        bc.slice, ...
        bc.pos, ...
        bc.fov, ...
        stitching_tolerance_pxl, ...
        numcores);
end

% make sparse matrix
conn.mat=sparse(mat_row, ...
    mat_col, ...
    true, ...
    numel(filt_neurons.uid), ...
    numel(bc.uid)); % connectiivty matrix with allowed mismatches, and only best matches
clear alignment
% toc
fprintf('Done.\n');

% % this step seems slow, but optimize in the future
% tic
% conn.mat(conn.match>0&conn.match~=max(conn.match,[],1))=0;% remove all connections that are not highest scored for that barcode\
% toc


%% filter out barcodes by slice and by duplicates.
% note that filtering by duplicates should be done after dupliates are
% removed in filt_neurons, especially those due to technical reasons (i.e.,
% a neuron imaged twice in overlaping FOVs, or neurons sliced in half.

% check slices
% tic
if ~isempty(filterslice)
    fprintf('Filtering barcodes based on slice numbers.\n');
    [neuron_idx1, bc_idx1]=find(conn.mat);
    neuron_slice=filt_neurons.slice(neuron_idx1);
    bc_slice=bc.slice(bc_idx1);
    remove_slice=zeros(numel(neuron_idx1),1);
    for m=1:numel(filterslice)
        remove_slice(bc_slice==neuron_slice+filterslice(m))=1;
    end
    mat_row=neuron_idx1(~remove_slice);
    mat_col=bc_idx1(~remove_slice);

    %check if filtering overlapping rolonies
    if stitching_tolerance_pxl>0
        [mat_row,mat_col]=filter_overlapping_bc(mat_row,mat_col, ...
            bc.slice, ...
            bc.pos, ...
            bc.fov, ...
            stitching_tolerance_pxl, ...
            numcores);
    end
    
    %make matrix
    conn.filt_mat=sparse(mat_row, ...
        mat_col, ...
        true, ...
        numel(filt_neurons.uid), ...
        numel(bc.uid)); % this is the connectiivty matrix with slice filtering applied
else
    fprintf('Skipping filtering based on slice numbers.\n');
    conn.filt_mat=conn.mat;
end
% toc

% filter out duplicates if applicable
if allow_duplicates==0
    conn.filt_mat(:,sum(conn.mat,1)>1)=0;
end

% filter out barcodes that are likely duplicates from neighboring FOVS that
% were imaged twice.




save(fname,'conn','-v7.3','-nocompression');

%% filter for axonal barseq

if ismember(lower(experiment_type),"axonalbarseq")
    fprintf('Filtering for axonal BARseq ...');
    h1=tic;
    assert(isfield(filt_neurons,coor)&&isfield(bc,coor),...
        sprintf('%s does not exist in either filt_neurons or bc.',coor));
    
    rol_num=sum(conn.mat,2);
    idx=find(rol_num>0);% all neurons with matching bacodes

    filt_neurons_slice=filt_neurons.slice(idx); % for slicing in parfor, sorted by idx already
    conn_mat=conn.mat(idx,:); % for slicing in parfor, sorted by idx already
    bc_coor=bc.(coor);
    %bc_coor(bc_coor<1)=1; % fix negative coordinates. Maybe these should be filtered out.
    bc_slice=bc.slice;
    save_pairs={};% used to record which neurons to save
    
    %numel(idx)
    parfor(n=1:numel(idx),round(maxcore/2))
        %tic
        soma_slice=filt_neurons_slice(n);
        %[~,I]=ismember(soma_slice,uniq_slice);
        bcid=find(conn_mat(n,:)&min(bc_coor,[],2)'>0); % this is actually index, not id. Filtering out any barcodes with negative coordinates. This shouldn't be in the brain in any case
        trusted_bcid=bcid(~ismember(bc_slice(bcid),suspicious_slides{soma_slice}));
        %suspicious_bcid=bcid(ismember(bc.slice(bcid),suspicious_slides{soma_slice}));


        % method 1: define a confidence map using barcodes from trusted slices
        conf_map=full(sparse(ceil(bc_coor(trusted_bcid,2)/rescale), ...
            ceil(bc_coor(trusted_bcid,1)/rescale), ...
            1, ...
            ceil(max(bc_coor(:,2))/rescale), ...
            ceil(max(bc_coor(:,1))/rescale)));
        conf_map_blur=imgaussfilt((conf_map), rad_pxl);
        conf_score=conf_map_blur( ...
            sub2ind(size(conf_map_blur), ...
            ceil(bc_coor(bcid,2)/rescale), ...
            ceil(bc_coor(bcid,1)/rescale)));
        %t11=toc;

        % method 2: find rolonies that look like lines in rolony enriched
        % areas
        % 2.1: hough line transform
        rol_map=full(sparse(ceil(bc_coor(bcid,2)/rescale), ...
            ceil(bc_coor(bcid,1)/rescale), ...
            1, ...
            ceil(max(bc_coor(:,2))/rescale), ...
            ceil(max(bc_coor(:,1))/rescale)));
        rol_map3=imclose(imdilate(rol_map>0, ...
            strel('disk',round(dens_disk_rad_pxl))), ...
            strel('disk',round(dens_disk_rad_pxl*1.5))); % make rolonies bigger, then connect them
        [H,T,R]=hough(rol_map3); % hough transform, then find peaks and fine lines
        P=houghpeaks(H,numpeaks, ...
            'threshold',ceil(houghthresh*max(H(:))), ...
            'NHoodSize',nhoodsize_pxl);
        hlines = houghlines(rol_map3,T,R,P, ...
            'FillGap',fillgap_pxl, ...
            'MinLength',minlength_pxl);
        %t12=toc;

        % 2.2: calculate distance between each rol and each hough line.
        dist2lines=zeros(numel(bcid),numel(hlines));
        for nnn=1:numel(hlines)
            dist2lines(:,nnn)=dist2linesegment(bc_coor(bcid,:)/rescale, ...
                hlines(nnn).point1, ...
                hlines(nnn).point2);
        end
        min_dist=min(dist2lines,[],2);

        % 2.3: make a density map of all rolonies of that bc
        dens_map=full(sparse(ceil(bc_coor(bcid,2)/rescale), ...
            ceil(bc_coor(bcid,1)/rescale), ...
            1, ...
            ceil(max(bc_coor(:,2))/rescale), ...
            ceil(max(bc_coor(:,1))/rescale)));
        dens_map_blur=imgaussfilt(dens_map, rad_pxl);
        thresh_val1=prctile(dens_map_blur(:),prc_thresh); % only pixels below this value in dens_map_blur are used to find background variations.
        thresh_val=std(dens_map_blur(dens_map_blur<thresh_val1))*10+median(dens_map_blur(:)); % background variations * 10 + median of variations. This is higher in "noisy" images.
        dens_score=dens_map_blur( ...
            sub2ind(size(dens_map_blur), ...
            ceil(bc_coor(bcid,2)/rescale), ...
            ceil(bc_coor(bcid,1)/rescale)));
        %t13=toc;

        % 2.4: save points if it's a trusted slice, if has trusted bc
        % nearby (2.1), or [is close to a hough line (2.2) AND in a
        % dense area (based on all bc) relative to the rest of the slice (2.3)].

        save_point=ismember(bcid',trusted_bcid) | ... % is in one of the non-problematic slices
            conf_score>conf_thresh | ... % or in an area with concentration of confident rolonies
            min_dist<=hough_dist_thresh/rescale&dens_score>thresh_val; % or on a line and is concentrated
        if sum(save_point)>0
            save_pairs{n,1}=[repmat(idx(n),sum(save_point),1),bcid(save_point)'];
        else
            save_pairs{n,1}=[];
        end

        %t14=toc;
        % if rem(n,10)==0
        %     fprintf('cell %u took %.2g step1, %.2g step2.1, %.2g step2.3, %.2g complete.\n', ...
        %         n, t11, t12, t13, t14);
        % end
    end
    save_pairs1=cell2mat(save_pairs);
    mat_row=save_pairs1(:,1);
    mat_col=save_pairs1(:,2);

    %check if filtering overlapping rolonies
    if stitching_tolerance_pxl>0
        [mat_row,mat_col]=filter_overlapping_bc(mat_row,mat_col, ...
            bc.slice, ...
            bc.pos, ...
            bc.fov, ...
            stitching_tolerance_pxl, ...
            numcores);
    end

    %make matrix
    conn.barseq_mat=sparse(mat_row, ...
        mat_col, ...
        true, ...
        size(conn.mat,1), ...
        size(conn.mat,2));
    t10=toc(h1);
    fprintf('Done, took %u secs.\n',t10);


end

%%
save(fname,'conn','-v7.3','-nocompression');
fprintf('Saved connectivity to %s.\n', fname);

end




% from bowtie2 manual:
%
% AS:i:<N>
% Alignment score. Can be negative. Can be greater than 0 in --local mode
% (but not in --end-to-end mode). Only present if SAM record is for an
% aligned read.
%
% XS:i:<N>
% Alignment score for the best-scoring alignment found other than the
% alignment reported. Can be negative. Can be greater than 0 in --local
% mode (but not in --end-to-end mode). Only present if the SAM record is
% for an aligned read and more than one alignment was found for the read.
% Note that, when the read is part of a concordantly-aligned pair, this
% score could be greater than AS:i.
%
% YS:i:<N>
% Alignment score for opposite mate in the paired-end alignment. Only
% present if the SAM record is for a read that aligned as part of a
% paired-end alignment.
%
% XN:i:<N>
% The number of ambiguous bases in the reference covering this alignment.
% Only present if SAM record is for an aligned read.
%
% XM:i:<N>
% The number of mismatches in the alignment. Only present if SAM record
% is for an aligned read.
%
% XO:i:<N>
% The number of gap opens, for both read and reference gaps, in the
% alignment. Only present if SAM record is for an aligned read.
%
% XG:i:<N>
% The number of gap extensions, for both read and reference gaps, in the
% alignment. Only present if SAM record is for an aligned read.
%
% NM:i:<N>
% The edit distance; that is, the minimal number of one-nucleotide edits
% (substitutions, insertions and deletions) needed to transform the read
% string into the reference string. Only present if SAM record is for an
% aligned read.
%
% YF:Z:<S>
% String indicating reason why the read was filtered out. See also:
% Filtering. Only appears for reads that were filtered out.
%
% YT:Z:<S>
% Value of UU indicates the read was not part of a pair. Value of CP
% indicates the read was part of a pair and the pair aligned concordantly.
% Value of DP indicates the read was part of a pair and the pair aligned
% discordantly. Value of UP indicates the read was part of a pair but the
% pair failed to aligned either concordantly or discordantly.
%
% MD:Z:<S>
% A string representation of the mismatched reference bases in the
% alignment. See SAM Tags format specification for details. Only present
% if SAM record is for an aligned read.






    function [neuron_idx,bc_idx]=filter_overlapping_bc(neuron_idx,bc_idx,bc_slice,bc_pos,bc_fov,stitching_tolerance_pxl,maxcore)
        %mat=conn.mat;
        %stitching_tolerance_pxl=5;
        cell_idx=unique(neuron_idx);
        %rol_num=sum(mat,2);
        %cell_idx=find(rol_num>0);

        to_remove=cell(1,numel(cell_idx));
        %mat_sorted=mat(cell_idx,:);
        %processed=zeros(numel(rol_num),1);

        %tic
        parfor(nnn=1:numel(cell_idx),maxcore)
            %if processed(nnn)==0
            %bcid=find(mat_sorted(nnn,:));
            bcid=bc_idx(neuron_idx==cell_idx(nnn))';
            if size(bcid,2)<size(bcid,1)
                bcid=permute(bcid,[2,1]);
            end
            uniq_slice=unique(bc_slice(bcid));
            to_remove{nnn}=[];
            for mmm=1:numel(uniq_slice)
                bcid_inslice=bcid(bc_slice(bcid)==uniq_slice(mmm));

                % calculate pdist for every barcode, .* dist in FOV
                d=pdist(bc_pos(bcid_inslice,:));
                diff_fov=pdist(double(bc_fov(bcid_inslice)),'hamming');
                d_idx=tril(true(numel(bcid_inslice),numel(bcid_inslice)),-1);
                [row,~]=ind2sub(size(d_idx),find(d_idx));
                to_remove{nnn}=[to_remove{nnn},bcid_inslice(row(diff_fov&d<stitching_tolerance_pxl))]; % always remove the latter item in pdist output.
            end
            %processed(mat_sorted(:,bcid(1)))=1;
            %end

        end
        remove_flag=ismember(bc_idx,cell2mat(to_remove));
        neuron_idx=neuron_idx(~remove_flag);
        bc_idx=bc_idx(~remove_flag);
        %mat(:,unique(cell2mat(to_remove)))=0;
        %toc
    end



    function myfastawrite(fname,bcid,bcC)
        %% custom way of writing fasta, 10x faster than built in for large arrays.
        % make parts
        init=repmat('>',size(bcid,1),1);
        nl=repmat(newline,size(bcid,1),1);

        % combine them
        joint=init+string(bcid)+nl+bcC+nl+nl;

        % convert to a single char
        joint2=char(joint);
        joint3=reshape(joint2',1,[]);

        %write to file
        if isfile(fname)
            warning('%s exists, overwriting.\n',fname)
        end

        fid=fopen(fname,'w');
        fprintf(fid,'%s',joint3);
        fclose(fid);
        %toc

    end



    function d = dist2linesegment(pt, v1, v2)
        % calculate distance of pt to a line segment defined by endpoints v1/v2.
        % pt can be a MxN matrix representing M points. v1 and v2 must be single
        % points
        if size(pt,2)==2
            pt(:,3)=0;
            v1(:,3)=0;
            v2(:,3)=0;
        end


        a = repmat(v1 - v2,size(pt,1),1);
        b = pt - v2;
        c = pt - v1;

        d1=sqrt(sum(b.^2,2));
        d2=sqrt(sum(c.^2,2));

        ab_isacute=(sum(a.^2,2)+d1.^2)>d2.^2;
        ac_isacute=(sum(a.^2,2)+d2.^2)>d1.^2;

        d0=sqrt(sum(cross(a,b,2).^2,2))./sqrt(sum(a.^2,2));

        % if both angles are acute, then use d0. If not, then use shorter distance
        % to one of the two ends.
        d=(ab_isacute&ac_isacute).*d0+(~ab_isacute|~ac_isacute).*min(d1,d2);

    end