function [seq_v_ok, V] = match_Vs(seq, genV, V_minimum_alignment_length, V_maximum_deletion, negative_excess_deletions_max, V_allowed_errors, error_cost)
%#codegen
%This function goes through each genomic V, finds best alignments, sorts these alignments by score and then returns the ones that pass certain filters.

align_length = zeros(length(genV),1); % the length of alignment for each allele
align_position = zeros(2,length(genV)); % V.align_position(1,v) contains the first matched nt in seq, so it's always 1 ; V.align_position(2,:) the first matched nt in genV.seq{v}
n_errors=zeros(length(genV),1); % number of mismatches starting from position 1.
error_positions = zeros(length(genV),V_allowed_errors); % positions of these mismatches on the sequence.

% number of deletions from genomic V for this alignment.
% It's called minimum because the number of deletions could be greater with the matched nucleotides randomly inserted.
% Of course, the number of deletions could also be smaller because of
% sequencing error. These would be 'negative' deletions.
min_deletions=zeros(length(genV),1);


v_large_deletion_flag=zeros(length(genV),1); % Flag if min_deletions is too large.

% Loop through alleles.
for v=1:length(genV)
    % Get highest scoring alignment for this allele, with acceptable number
    % of errors
    [align_length(v), align_position(:,v), n_errors(v), error_positions(v,:)] = align_with_constraints_fast_left(seq, genV(v).seq, V_allowed_errors, error_cost,[length(genV(v).seq):-1:1]');
    % Calculate the number of deletions implied by the alignment.
    min_deletions(v) = length(genV(v).seq)-(align_position(2,v) + align_length(v)-1);

    % Flag too many deletions
    if(min_deletions(v) > V_maximum_deletion)
        v_large_deletion_flag(v) = 1;
        %continue;
    end

end
% Check!
assert(all(min_deletions >=0));

% sort the genomic Vs in descending order of 'score' and ascending order of min_deletions
scores  = align_length - error_cost*n_errors;
S = [-scores, min_deletions];
[~,order]=sortrows(S);

% Set a score threshold for alleles.
min_score = max(scores) - 3*V_minimum_alignment_length;
max_errors = n_errors(order(1)) + 3;

% Subset of alleles that have high enough score, long enough alignment, and
% not too many deletions.
v_list = find( scores(order) >= min_score & align_length(order) >= V_minimum_alignment_length & n_errors(order) <= max_errors & v_large_deletion_flag(order)==0 );

if isempty(v_list)
    %disp('No reasonable V match for sequence.')
    %disp(seq);
    ok_order = [];
    seq_v_ok = false; % Sequence not ok!
else
    ok_order = order(v_list);
    seq_v_ok = true; % Sequence ok!
end

if seq_v_ok

    % Store all the alignment information in V for acceptable alleles
    V.align_length = align_length(ok_order);
    V.align_position =  align_position(:,ok_order);
    V.min_deletions = min_deletions(ok_order);
    V.n_errors = n_errors(ok_order);
    V.error_positions = error_positions(ok_order,:);
    V.p_region_max_length = zeros(numel(ok_order),V_maximum_deletion + 1);
    V.excess_error_positions = zeros(numel(ok_order),negative_excess_deletions_max);
    V.alleles_all = ok_order;

    % Find gene indices of the acceptable alleles.
    %%%codegen compatible version of:
    %%%genV_ok_gene_inds = [genV(ok_order).gene_index];
    genV_ok = genV(ok_order);
    genV_ok_gene_inds = zeros(1,numel(ok_order));
    for g=1:numel(ok_order)
        genV_ok_gene_inds(g) = genV_ok(g).gene_index;
    end

    % Get list of distinct genes from acceptable alleles.
    [~,Vg] = unique(genV_ok_gene_inds,'first' ); % Pick the first allele listed for each gene.
    Vg = sort(Vg); % sort it because unique returns in ascending order of ok_order. This puts it back in order of best match.
    % Only one (best) allele from each gene
    V.alleles_from_distinct_genes = ok_order(Vg');

    %Now go through and find palindromic nucleotides for various number of deletions,
    %as well as error positions for 'negative' deletions.
    for v=1:numel(ok_order)
        % Loop over number of deletions (nd is actual number of genomic deletions).
        % Lower bound is maximum of 0 and min_deletions - negative_excess_deletions_max (usually 3).
        % Upper bound is minimum of V_maximum_deletion (usually 12) and the never actually possible case of the whole alignment being random insertions
        for nd= (max(0,V.min_deletions(v) - negative_excess_deletions_max)):min(V_maximum_deletion, V.align_length(v) + V.min_deletions(v))

            % For each value of deletions, find longest half-palindrome
            % from the implied end of the gene sequence
            p=0;
            still_palindrome=true;
            % Upper bound on p is sequence length left on the right beyond end of V match (accounting for deletions) OR the never actually possible case of the length of the V match (accounting for deletions)
            while still_palindrome && p < min( length(seq) - V.align_length(v) + nd - V.min_deletions(v) , V.align_length(v) - (nd - V.min_deletions(v)) ) % Upper bound is because of some constraints
                still_palindrome = genV(ok_order(v)).seq(V.align_position(2,v) + V.align_length(v) -1 - p - (nd - V.min_deletions(v)) ) == dnacomplement(seq(V.align_length(v) + p - (nd - V.min_deletions(v)) + 1) );
                if still_palindrome
                    p = p + 1;
                end
            end

            % Store length of longest half-palindrome for each value of
            % deletions.
            V.p_region_max_length(v,1+nd) = p;
        end

        % Now for the 'negative' deletions, store where the mismatches are.
        % This is so we can count number of errors for these possibilities.

        % Number of 'negative' deletions to consider:
        n_excess = min([negative_excess_deletions_max, V.min_deletions(v), length(seq) - V.align_length(v) ]);

        % Find mismatches between sequence and genomic V at those positions.
        excess_err_pos = V.align_length(v) + find((genV(ok_order(v)).seq((V.align_position(2,v)+V.align_length(v)): (V.align_position(2,v)+V.align_length(v) + n_excess-1)  )~=seq( (V.align_length(v)+1): (V.align_length(v) + n_excess) )));
        V.excess_error_positions(v,1:length(excess_err_pos)) = excess_err_pos;
    end

else

    % Since sequence is not ok, return empty values.
    V.align_length = [];
    V.align_position =  [];
    V.min_deletions = [];
    V.n_errors = [];
    V.error_positions = [];
    V.p_region_max_length = [];
    V.excess_error_positions = [];
    V.alleles_all = [];
    V.alleles_from_distinct_genes = [];


end


end
