function assigns = Mus_V_DJ_cuts_insertions_dinuc_ntbias_model_assignments(model, seq, V, D, J, genV, genD, genJ, assigns, probability_threshold_factor, no_error, ignore_deep_error, do_smoothing, force_all_alleles, READ_LENGTH_CORRECTION )
%#codegen
% Function that calculates probabilities of various assignments for a given
% sequence.
%disp('Aaaaaa')
nseq = seq2nseq(seq)';
% probability_threshold_factor=.001; % assignments with probability < this factor* highest probability so far, are ignored.
log_probability_threshold_factor = log(probability_threshold_factor);
log_probability_hopeless_threshold_factor = log_probability_threshold_factor + log(10^(-6));
max_skips = 300000;
in=0; % assignment number
%max_proba_index = 0;

skips = 0; % number of assignments that have been skipped over for being low probability.
deep_error_limit = 50;
max_J_depth = model.max_J_deletions;
max_V_depth = length(seq) + model.max_V_deletions;
J_max_error = 0;
D_max_error = 0;

assume_palindrome_negative_deletions = true;
log_max_model_p_nt_DJ  = log(max(model.RnucleotideDJ_per_nucleotideDJ_3prime,[],2));
log_max_model_p_nt_VD  = log(max(model.RnucleotideVD_per_nucleotideVD_5prime,[],2));
log_max_model_p_nt = log(max([model.RnucleotideDJ_per_nucleotideDJ_3prime(:); model.RnucleotideVD_per_nucleotideVD_5prime(:)]));

log_RnucleotideDJ_per_nucleotideDJ_3prime = log(model.RnucleotideDJ_per_nucleotideDJ_3prime(:))';
log_RnucleotideVD_per_nucleotideVD_5prime = log(model.RnucleotideVD_per_nucleotideVD_5prime(:))';

if model.Rerror_per_sequenced_nucleotide > 0
    log_Rerror_per_sequenced_nucleotide_divided_by_3 = log(model.Rerror_per_sequenced_nucleotide/3);
else
    log_Rerror_per_sequenced_nucleotide_divided_by_3 = -1000;
end

L_err = length(seq) - model.high_error_region;
log_proba_Rerror_normalization = -L_err*log(1+model.Rerror_per_sequenced_nucleotide);

log_max_model_pcutV_given_gene = log(max(model.PcutV_given_V,[],1));
log_max_model_pcutJ_given_gene = log(max(model.PcutJ_given_J,[],1));
log_max_model_pcutD_given_gene = log(max(reshape(model.PcutDlcutDr_given_D, size(model.PcutDlcutDr_given_D,1)^2, size(model.PcutDlcutDr_given_D,3) )) );

log_max_model_pinsVD = log(max(model.PinsVD));
log_max_model_pinsDJ = log(max(model.PinsDJ));
log_max_model_pins = log_max_model_pinsVD + log_max_model_pinsDJ;

% if do_smoothing is false, only one assignment (hopefully the most likely)
% is considered for the sequence.
if do_smoothing && ~no_error
    % the loop over deletions starts from 'negative' deletions.
    nd_start = -model.negative_excess_deletions_max; % usually -3.
    np_start_from_max = false;
else
    % don't consider negative deletions
    nd_start = 0;
    % use max possible palindrome
    np_start_from_max = true;
    
end

% highest probability so far
log_highest_probability = -1000;
best_D_align_length = 0;
% Main loop over assignments


max_nerrorsv_1 = 0;
% Loop over V alleles
for v= 1:length(V.alleles_all)
    
    % if too many skips, quit
    if skips > max_skips
        assigns.n_assignments = in;
        assigns.skips = skips;
        return;
    end
    
    v_a = V.alleles_all(v); % allele number
    v_g = genV(v_a).gene_index; % gene number
    
    if ignore_deep_error % this is sort of a hack,I only use it for the overlapping seq calculation.
        high_error_region = max(model.high_error_region, V.align_length(v) - model.max_excess_V_deletions - deep_error_limit);
    else
        high_error_region = model.high_error_region;
    end
    
    if v==1
        v_err_pos = V.error_positions(v,:); % error positions in V alignment
        deletable_errors = sum(v_err_pos > (V.align_length(v) - model.max_excess_V_deletions));
        min_nerrorsv = V.n_errors(v) - deletable_errors - sum(v_err_pos <= high_error_region & v_err_pos > 0);
        max_nerrorsv_1 = V.n_errors(v) - sum(v_err_pos <= high_error_region & v_err_pos > 0);
        
        if no_error && min_nerrorsv > 0
            continue;
        end
        
        
    else
        % Check if this V has too many more errors than the best V
        v_err_pos = V.error_positions(v,:); % error positions in V alignment
        deletable_errors = sum(v_err_pos > (V.align_length(v) - model.max_excess_V_deletions));
        min_nerrorsv = V.n_errors(v) - deletable_errors - sum(v_err_pos <= high_error_region & v_err_pos > 0);
        max_nerrorsv = V.n_errors(v) - sum(v_err_pos <= high_error_region & v_err_pos > 0);
        if no_error && min_nerrorsv > 0
            continue;
        end
        
        log_min_diff_perrv = (max_nerrorsv - max_nerrorsv_1)*log_Rerror_per_sequenced_nucleotide_divided_by_3;
        if log_min_diff_perrv < log_probability_threshold_factor
            %disp('skip 108')
            skips = skips + 1;
            continue;
        end
        
    end
    
    log_max_pcutV_loop_v = log_max_model_pcutV_given_gene(v_g);
    
    log_highest_probability_GIVEN_current_V_allele = -1000; % Highest probability so far GIVEN current D allele AND outer loop variables {V, J alleles}
    
    
    v_break_out = false;
    n_assignments_v_gene = 0;
    % Loop over J alleles
    for j = 1:length(J.alleles_all)
        
        if skips > max_skips
            assigns.n_assignments = in;
            assigns.skips = skips;
            return;
        end
        
        if  J.align_length(j) < model.min_J_align_length
            continue;
        end
        
        j_a = J.alleles_all(j); % allele number
        j_g = genJ(j_a).gene_index; % gene number
        
        niVD_DJ0 = J.align_position(1,j) - V.align_length(v) -1;
        log_max_pcutJ_loop_j = log_max_model_pcutJ_given_gene(j_g);
        
        
        log_highest_probability_GIVEN_current_J_allele = -1000; % Highest probability so far GIVEN current D allele AND outer loop variables {V, J alleles}
        n_assignments_j_gene = 0;
        j_break_out = false;
        % Loop over D alleles
        for d_i=1:length(D.allele_order)
            
            if skips > max_skips
                assigns.n_assignments = in;
                assigns.skips = skips;
                return;
            end
            
            d = D.allele_order(d_i);
            %             d
            % Skip if D is short compared to best alignment
            %             if d_i > 1 && D.align_length(1,d) <= (D.align_length(1,D.allele_order(1)) - 5)
            %                 break;
            %             end
            
            d_a = d; % allele number
            d_g = genD(d_a).gene_index; % gene number
            l_d_seq = length(genD(d_a).seq);
            
            n_assignments_d_gene = 0;
            d_break_out = false;
            
            % Base Probability of gene choices, multiplied by conditional probabilities for allele choices, given genes.
            % J has no distinguishable alleles (differences are beyond sequencing primer ends)
            probabase = model.PV(v_g)*model.PDJ(d_g,j_g)*model.PVallele_given_gene(genV(v_a).allele, v_g)*model.PDallele_given_gene(genD(d_a).allele, d_g);
            
            p_max_Dl = max(flat(D.p_region_max_length_left(:,1,d)));
            p_max_Dr = max(flat(D.p_region_max_length_right(:,1,d)));
            p_max_J = (J.p_region_max_length(j,1));
            % if not zero deletions, no pals.
            if J.min_deletions(j)>2
                p_max_J=0;
            end
            % Upper bound on insertions nt bias factor
            niVD_DJ_min = max(0, niVD_DJ0 - D.align_length(1,d) - V.p_region_max_length(v,1) - p_max_J - p_max_Dl - p_max_Dr );
            log_p_max_nt_VD_DJ_d_allele = (niVD_DJ_min)*log_max_model_p_nt;
            
            
            log_max_pcutD_loop_d = log_max_model_pcutD_given_gene(d_g);
            
            log_max_pcutVDJ_loop_d = log_max_pcutV_loop_v + log_max_pcutJ_loop_j + log_max_pcutD_loop_d;
                        
            
            if probabase==0 || ( (log(probabase) + log_max_pcutVDJ_loop_d + log_max_model_pins + log_p_max_nt_VD_DJ_d_allele) < (log_highest_probability + log_probability_threshold_factor) )
                %disp('skip 185')
                
                skips = skips + 1;
                continue;
            end
            
            log_probabase = log(probabase);
            

            log_highest_probability_GIVEN_current_D_allele = -1000; % Highest probability so far GIVEN current D allele AND outer loop variables {V, J alleles}
            % Loop over V deletions.
            % ndV1 is the deviation from the no. of deletions implied by the alignment (min_deletions).
            % Lower bound is nd_start (usually -3).
            % Upper bound is set by maximum deletions allowed and by minimum match length required.
            ndV_max = min([model.max_excess_V_deletions, V.align_length(v) - model.min_V_length,model.max_V_deletions - V.min_deletions(v)]);
            %ndV_max
            for ndV1=nd_start:ndV_max
                if skips > max_skips
                    assigns.n_assignments = in;
                    assigns.skips = skips;
                    return;
                end
                
                % Actual number of deletions
                ndV = V.min_deletions(v) + ndV1;

                % This is a violation of a constraint.
                % Happens for eg. when min_deletions is 0 and I consider 'negative' deletions.
                if ndV < 0 || ndV > model.max_V_deletions
                    % Go to next iteration of loop.
                    continue;
                end

                % Calculate number of errors in V section, by leaving out errors that are in deleted region for this assignment.
                % Also leaves out errors in left most 'high_error_region'.
                
                if ignore_deep_error
                    high_error_region = max(model.high_error_region, V.align_length(v) - ndV1 - deep_error_limit);
                else
                    high_error_region = model.high_error_region;
                end
                
                v_err_pos = V.error_positions(v,:); % error positions in V alignment
                v_err_excess_pos = V.excess_error_positions(v,:); % error positions in extended V alignment for 'negative' deletions
                if ndV1 < 0
                    % Case of 'negative' deletions
                    v_ex_errs_i = v_err_excess_pos <= (V.align_length(v) - ndV1) & v_err_excess_pos > 0 ;
                    v_ex_errs = sum(v_ex_errs_i);
                    if v_ex_errs > 1, % If the nts beyond the 1st -ve nt don't match, just skip.
                        continue;
                    end
                    nerrorsv = V.n_errors(v) + v_ex_errs - sum(v_err_pos <= high_error_region & v_err_pos > 0);
                else
                    % Case of positive deletions
                    v_ex_errs=0;
                    v_ex_errs_i = [];
                    nerrorsv = V.n_errors(v) - sum(v_err_pos > (V.align_length(v) - ndV1)) - sum(v_err_pos <= high_error_region & v_err_pos > 0);
                end
                
                if no_error && nerrorsv > 0
                    continue;
                end
                
                log_perrv = nerrorsv*log_Rerror_per_sequenced_nucleotide_divided_by_3;
                
                
                if ndV==0 || ~assume_palindrome_negative_deletions
                    npV_max = min(model.max_palindrome, V.p_region_max_length(v, 1+ndV));
                    npV_potential_max = npV_max;
                else
                    npV_max = 0;
                    npV_potential_max = min(model.max_palindrome, V.p_region_max_length(v, 1+ndV));
                end
                
                
                %Upper bound on insertions nt bias factor
                
                niVD_DJ_min = max(0, niVD_DJ0 + ndV1 - D.align_length(1,d) - npV_max - p_max_J - p_max_Dl - p_max_Dr);
%                 niVD_DJ_min
%                 niVD_DJ0
%                 ndV1
%                 npV_max
%                 p_max_J
%                 p_max_Dl
%                 p_max_Dr
%                 D.align_length(1,d)
                
                
                log_p_max_nt_VD_DJ_dV = (niVD_DJ_min)*log_max_model_p_nt;
                
                log_highest_probability_GIVEN_current_V_deletions = -1000; % Highest probability so far GIVEN current V deletions AND outer loop variables {V, J, D alleles}                

                test_proba = log_probabase + log_perrv + log_max_pcutVDJ_loop_d + log_max_model_pins + log_p_max_nt_VD_DJ_dV;
%                 if(test_proba > log_highest_probability_GIVEN_current_D_allele), log_highest_probability_GIVEN_current_D_allele = test_proba; end;
%                 if(test_proba > log_highest_probability_GIVEN_current_V_allele), log_highest_probability_GIVEN_current_V_allele = test_proba; end;
%                 if(test_proba > log_highest_probability_GIVEN_current_J_allele), log_highest_probability_GIVEN_current_J_allele = test_proba; end;
%                 if(test_proba > log_highest_probability_GIVEN_current_V_deletions), log_highest_probability_GIVEN_current_V_deletions = test_proba; end;
                
                if test_proba < log_probability_threshold_factor + log_highest_probability
                    %disp('skip 272')
                    skips = skips + 1;
                    continue;
                end
                
                V_align_length = V.align_length(v) - ndV1 - READ_LENGTH_CORRECTION;
                

                
                
                
                % Loop over J right deletions.
                % ndJ1 is the deviation from the no. of deletions implied by the alignment (min_deletions).
                % Lower bound is nd_start (usually -3).
                % Upper bound is set by maximum deletions allowed and by minimum match length required
                ndJ_max = min([model.max_excess_J_deletions,J.align_length(j)-model.min_J_assign_length,model.max_J_deletions - J.min_deletions(j)]);
                ndJ_max;
                for ndJ1=nd_start:ndJ_max
                    if skips > max_skips
                        assigns.n_assignments = in;
                        assigns.skips = skips;
                        return;
                    end
                    
                    
                    
                    % actual number of J deletions
                    ndJ = J.min_deletions(j) + ndJ1;
                    % This is a violation of a constraint.
                    % Happens for eg. when min_deletions is 0 and I consider 'negative' deletions.
                    if ndJ < 0 || ndJ > model.max_J_deletions
                        % Go to next iteration of loop.
                        continue;
                    end
                    
                    
                    niVD_DJ_total_with_ps = niVD_DJ0 + ndV1 + ndJ1;
                    % Violation of constraint. V and J overlap in this assignment. So move on until this is not true.
                    if niVD_DJ_total_with_ps < 0
                        continue;
                    end
                    
                    
                    % Calculate number of errors in J section, by leaving out errors that are in deleted region for this assignment.
                    j_err_pos = J.error_positions(j,:); % error positions in J alignment
                    j_err_excess_pos = J.excess_error_positions(j,:); % error positions in extended J alignment for 'negative' deletions
                    if ndJ1 < 0
                        % Case of 'negative' deletions
                        j_ex_errs_i = j_err_excess_pos >= (J.align_position(1,j) + ndJ1);
                        j_ex_errs = sum(j_ex_errs_i);
                        if j_ex_errs > 1
                            continue;
                        end
                        nerrorsj = J.n_errors(j) + j_ex_errs;
                    else
                        % Case of positive deletions
                        j_ex_errs = 0;
                        j_ex_errs_i = [];
                        nerrorsj = sum( j_err_pos >= (J.align_position(1,j) + ndJ1));
                    end
                    
                    if no_error && nerrorsj > 0
                        continue;
                    end
                    
                    
                    
                    if nerrorsj > J_max_error || nerrorsj/(J.align_length(j)-ndJ1) > 0.3
                        continue;
                    end
                    
                    
                    log_perrj = nerrorsj*log_Rerror_per_sequenced_nucleotide_divided_by_3;
                    
                    %                     % Upper bound on insertions nt bias factor
                    niVD_DJ_min = max(0, niVD_DJ0 + ndV1 + ndJ1 - D.align_length(1,d) - V.p_region_max_length(v,1 + ndV) - J.p_region_max_length(j,1 + ndJ) - p_max_Dl - p_max_Dr);
                    log_p_max_nt_VD_DJ_dJ = (niVD_DJ_min)*log_max_model_p_nt;
                    
                    log_highest_probability_GIVEN_current_J_deletions = -1000;                   
                    
                    test_proba = log_probabase + log_perrv + log_perrj + log_max_pcutVDJ_loop_d + log_max_model_pins + log_p_max_nt_VD_DJ_dJ;
%                     if(test_proba > log_highest_probability_GIVEN_current_D_allele), log_highest_probability_GIVEN_current_D_allele = test_proba; end;
%                     if(test_proba > log_highest_probability_GIVEN_current_V_allele), log_highest_probability_GIVEN_current_V_allele = test_proba; end;
%                     if(test_proba > log_highest_probability_GIVEN_current_J_allele), log_highest_probability_GIVEN_current_J_allele = test_proba; end;
%                     if(test_proba > log_highest_probability_GIVEN_current_V_deletions), log_highest_probability_GIVEN_current_V_deletions = test_proba; end;
%                     if(test_proba > log_highest_probability_GIVEN_current_J_deletions), log_highest_probability_GIVEN_current_J_deletions = test_proba; end;
                    

%                     ndJ1
                    if test_proba < log_probability_threshold_factor + log_highest_probability
                        %disp('skip 360')
                        continue;
                    end
                    
                    % Maximum possible half-lengths of palindromes, given current deletions values
                    if ndJ==0 || ~assume_palindrome_negative_deletions
                        npJ_max = min(model.max_palindrome, J.p_region_max_length(j, 1+ndJ));
                        npJ_potential_max = npJ_max;
                    else
                        npJ_max = 0;
                        npJ_potential_max = min(model.max_palindrome, J.p_region_max_length(j, 1+ndJ));
                    end
                    

                    
                    if np_start_from_max
                        np_step = -1;
                        npV_start = npV_max;
                        npV_end = 0;
                        npJ_start = npJ_max;
                        npJ_end = 0;
                        
                    else
                        np_step = 1;
                        npV_start = 0;
                        npV_end = npV_max;
                        npJ_start = 0;
                        npJ_end = npJ_max;
                    end
                    
                    J_align_length = J.align_length(j) - ndJ1;
                    
                    
                    % Loop over half-length of V palindrome
                    for npV=npV_start:np_step:npV_end
                        if skips > max_skips
                            assigns.n_assignments = in;
                            assigns.skips = skips;
                            return;
                        end
                        
                        
                        V_end = V.align_length(v) - ndV1 + npV;
                        
                        % V with palindrome overlaps J.
                        if V_end >= (J.align_position(1,j) + ndJ1)
                            continue;
                        end
                        
                        % Cut variable
                        ncutV = ndV - npV;
                        
                        PcutV = model.PcutV_given_V( -model.min_V_cut +  1 + ncutV, v_g);
                        
                        if PcutV==0
                            %disp('skip 415')
                            skips = skips + 1;
                            continue;
                        end
                        
                        log_PcutV = log(PcutV);
                        log_max_pcutVDJ_loop_pV = log_PcutV + log_max_pcutD_loop_d + log_max_pcutJ_loop_j;
                        
                        %                     % Upper bound on insertions nt bias factor
                        niVD_DJ_min = max(0, niVD_DJ0 + ndV1 + ndJ1 - D.align_length(1,d) - npV - J.p_region_max_length(j,1 + ndJ) - p_max_Dl - p_max_Dr);
                        log_p_max_nt_VD_DJ_pV = (niVD_DJ_min)*log_max_model_p_nt;
                        
                        test_proba = log_probabase + log_perrv + log_perrj + log_max_pcutVDJ_loop_pV + log_max_model_pins + log_p_max_nt_VD_DJ_pV;
%                         if(test_proba > log_highest_probability_GIVEN_current_D_allele), log_highest_probability_GIVEN_current_D_allele = test_proba; end;
%                         if(test_proba > log_highest_probability_GIVEN_current_V_allele), log_highest_probability_GIVEN_current_V_allele = test_proba; end;
%                         if(test_proba > log_highest_probability_GIVEN_current_J_allele), log_highest_probability_GIVEN_current_J_allele = test_proba; end;
%                         if(test_proba > log_highest_probability_GIVEN_current_V_deletions), log_highest_probability_GIVEN_current_V_deletions = test_proba; end;
%                         if(test_proba > log_highest_probability_GIVEN_current_J_deletions), log_highest_probability_GIVEN_current_J_deletions = test_proba; end;
                       
                        if test_proba < log_probability_threshold_factor + log_highest_probability
                            %disp('skip 435')
                            skips = skips + 1;
                            continue;
                        end
                        
                        % Loop over half-length of J palindrome
                        for npJ= npJ_start: np_step: npJ_end
                            if skips > max_skips
                                assigns.n_assignments = in;
                                assigns.skips = skips;
                                return;
                            end
                            
                            
                            J_start = J.align_position(1,j) + ndJ1 - npJ;
                            % V and J overlap
                            if V_end >= J_start
                                continue;
                            end
                            
                            
                            % Cut variable
                            ncutJ = ndJ - npJ;
                            PcutJ = model.PcutJ_given_J( -model.min_J_cut + 1 + ncutJ, j_g);
                            
                            if PcutJ==0
                                %disp('skip 461')
                                skips = skips + 1;
                                continue;
                            end
                            
                            log_PcutJ = log(PcutJ);
                            log_max_pcutVDJ_loop_pJ = log_PcutV + log_PcutJ + log_max_pcutD_loop_d;
                            
                            
                            n_D_aligns = D.n_alignments(d);
                            start_n_D_aligns = 1;
                            % find first completely valid alignment
                            completely_valid_na = false;
                            c_na=0;
                            while ~completely_valid_na && c_na < n_D_aligns
                                c_na = c_na + 1;
                                completely_valid_na = D.align_position_left(c_na,d) > V_end && D.align_position_right(c_na,d) < J_start;
                            end
                            
                            
                            if completely_valid_na
                                first_valid_length = D.align_length(c_na,d);
                            else
                                % find first partly valid alignment
                                partly_valid_na = false;
                                p_na=0;
                                while ~partly_valid_na && p_na < n_D_aligns
                                    p_na = p_na + 1;
                                    partly_valid_na = D.align_position_left(p_na,d) < J_start && D.align_position_right(p_na,d) > V_end;
                                end
                                
                                if partly_valid_na
                                    first_valid_length = min(J_start-1,D.align_position_right(p_na,d)) - max(V_end+1,D.align_position_left(p_na,d)) + 1;
                                    start_n_D_aligns = p_na;
                                else
                                    first_valid_length = 0;
                                    % do only ZERO D. No even partly valid alignment.
                                    start_n_D_aligns = n_D_aligns + 1;
                                end
                                
                            end
                            
                            best_D_align_length = max(best_D_align_length,first_valid_length);
                            do_zeroD = best_D_align_length <=8;
                            
                            
                            if do_zeroD
                                end_n_D_aligns = n_D_aligns + 1;
                            else
                                end_n_D_aligns = n_D_aligns;
                            end
                            
                            % NDN length
                            niVD_DJ_total = J_start - V_end - 1;
                            
                            % Upper bound on nt bias factor
                            niVD_DJ_min = max(0,niVD_DJ_total - first_valid_length - D.p_region_max_length_left(start_n_D_aligns,D.deletions_left(start_n_D_aligns,d)+1,d) - D.p_region_max_length_right(start_n_D_aligns,D.deletions_right(start_n_D_aligns,d)+1,d) );
                            log_p_max_nt_VD_DJ_pJ = niVD_DJ_min*log_max_model_p_nt;
                            
                            test_proba = log_probabase + log_perrv + log_perrj + log_max_pcutVDJ_loop_pJ + log_max_model_pins + log_p_max_nt_VD_DJ_pJ;
%                             if(test_proba > log_highest_probability_GIVEN_current_D_allele), log_highest_probability_GIVEN_current_D_allele = test_proba; end;
%                             if(test_proba > log_highest_probability_GIVEN_current_V_allele), log_highest_probability_GIVEN_current_V_allele = test_proba; end;
%                             if(test_proba > log_highest_probability_GIVEN_current_J_allele), log_highest_probability_GIVEN_current_J_allele = test_proba; end;
%                             if(test_proba > log_highest_probability_GIVEN_current_V_deletions), log_highest_probability_GIVEN_current_V_deletions = test_proba; end;
%                             if(test_proba > log_highest_probability_GIVEN_current_J_deletions), log_highest_probability_GIVEN_current_J_deletions = test_proba; end;
                            if test_proba < log_probability_threshold_factor + log_highest_probability
                                %disp('skip 528')
                                skips = skips + 1;
                                continue;
                            end
                            
                            
                           
                            for na = start_n_D_aligns:end_n_D_aligns;
                                %na

                                zeroD = na > n_D_aligns;
                                if ~zeroD && D.align_length(na,d) < (best_D_align_length - 8)
                                    break;
                                end

                                % If the whole D alignment lies within the V or J
                                if ~zeroD && (D.align_position_right(na,d) <= V_end || D.align_position_left(na,d) >= J_start)
                                    continue;
                                end
                                
                                
                                
                                if ~zeroD
                                    ndDl_min = D.deletions_left(na,d) + nd_start;
                                    % Can go all the way until 1 nt left. Not to zero.
                                    ndDl_max = D.deletions_left(na,d) + min([model.max_excess_D_deletions, model.max_D_deletions - D.deletions_left(na,d), D.align_length(na,d)-1]);
                                    
                                    d_errs_internal_max = D.n_errors(na,d);
                                    
                                else
                                    ndDl_min = 0;
                                    ndDl_max = min(model.max_D_deletions,length(genD(d_a).seq));
                                    d_errs_internal_max = 0;
                                end
                                
                                
                                
                                for ndDl=ndDl_min:ndDl_max
                                    
                                    
                                    if ndDl < 0 || ndDl > model.max_D_deletions
                                        continue;
                                    end
                                    
                                    
                                    if ~zeroD
                                        ndDl1 = ndDl - D.deletions_left(na,d);
                                        npDl_potential_max = min(model.max_palindrome, D.p_region_max_length_left(na, 1+ndDl,d));
                                    else
                                        ndDl1 = 0;
                                        npDl_potential_max = 0;
                                    end
                                    
                                    if ndDl==0
                                        npDl_max = npDl_potential_max;
                                    else
                                        npDl_max = 0;
                                    end
                                    
                                    % errors
                                    if ndDl1 < 0
                                        d_err_excess_pos_left = D.excess_error_positions_left(na,:,d); % error positions in extended J alignment for 'negative' deletions
                                        % Case of 'negative' deletions
                                        d_ex_errs_left_i = d_err_excess_pos_left >= D.align_position_left(na,d) + ndDl1;
                                        d_ex_errs_left = sum(d_ex_errs_left_i);
                                        if d_ex_errs_left > 1
                                            continue;
                                        end
                                    else
                                        % Case of positive deletions
                                        d_ex_errs_left = 0;
                                        d_err_excess_pos_left = [];
                                        d_ex_errs_left_i = [];
                                    end
                                    
                                    
                                    
                                    if ~zeroD
                                        ndDr_min = D.deletions_right(na,d) + nd_start;
                                        ndDr_max = D.deletions_right(na,d) + min([model.max_excess_D_deletions,model.max_D_deletions - D.deletions_right(na,d),(D.align_length(na,d) - ndDl1 -1 )]);
                                        
                                        %                                         niVD_DJ_min = max(0,niVD_DJ_total - D.align_length(na,d) + ndDl1 - npDl_max - p_max_Dr_na);
                                        
                                        insVD_nseq_min_loopdDl = nseq( V_end + 1 - READ_LENGTH_CORRECTION :  D.align_position_left(na,d) + ndDl1 - npDl_max - 1 - READ_LENGTH_CORRECTION);
                                        if length(insVD_nseq_min_loopdDl) > model.max_insertions
                                            continue;
                                        end
                                        min_mononucleotideVD = [sum(insVD_nseq_min_loopdDl==1); sum(insVD_nseq_min_loopdDl==2); sum(insVD_nseq_min_loopdDl==3); sum(insVD_nseq_min_loopdDl==4)];
                                        
                                        log_p_max_nt_VD_loop_Dl = min_mononucleotideVD(1)*log_max_model_p_nt_VD(1) + min_mononucleotideVD(2)*log_max_model_p_nt_VD(2) + min_mononucleotideVD(3)*log_max_model_p_nt_VD(3) + min_mononucleotideVD(4)*log_max_model_p_nt_VD(4);
                                        
                                    else
                                        % For zero D case, ndDr is determined by ndDl since we are doing a convolution.
                                        ndDr_min = length(genD(d_a).seq) - ndDl;
                                        ndDr_max = ndDr_min;
                                        
                                        log_p_max_nt_VD_loop_Dl = 0;
                                        
                                        
                                        % for compilation purposes
                                        insVD_nseq_min_loopdDl = [];
                                    end
                                    
                                    log_highest_probability_GIVEN_current_Dl_deletions = -1000;
                                    
                                    
                                    
                                    for ndDr = ndDr_min:ndDr_max
                                        
                                        
                                        if ndDr < 0 || ndDr > model.max_D_deletions
                                            continue;
                                        end
                                        
                                        if ~zeroD
                                            npDr_potential_max = min(model.max_palindrome,  D.p_region_max_length_right(na, 1 + ndDr, d));
                                        else
                                            npDr_potential_max = 0;
                                        end
                                        
                                        if ndDr==0
                                            npDr_max = npDr_potential_max;
                                        else
                                            npDr_max = 0;
                                        end
                                        
                                        
                                        
                                        if ~zeroD
                                            ndDr1 = ndDr - D.deletions_right(na,d);
                                            D_align_length = D.align_length(na,d) - ndDl1 - ndDr1;
                                            assert(D_align_length > 0);
                                            
                                            insDJ_nseq_min_loopdDr = nseq( D.align_position_right(na,d) - ndDr1 + npDr_max + 1 - READ_LENGTH_CORRECTION : J_start - 1 - READ_LENGTH_CORRECTION);
                                            if length(insDJ_nseq_min_loopdDr) > model.max_insertions
                                                continue;
                                            end
                                            min_mononucleotideDJ = [sum(insDJ_nseq_min_loopdDr==1); sum(insDJ_nseq_min_loopdDr==2); sum(insDJ_nseq_min_loopdDr==3); sum(insDJ_nseq_min_loopdDr==4)];
                                            log_p_max_nt_DJ_loop_Dr = min_mononucleotideDJ(1)*log_max_model_p_nt_DJ(1) + min_mononucleotideDJ(2)*log_max_model_p_nt_DJ(2) + min_mononucleotideDJ(3)*log_max_model_p_nt_DJ(3) + min_mononucleotideDJ(4)*log_max_model_p_nt_DJ(4);
                                            
                                            log_p_max_nt_VD_DJ_na_loop_Dr = log_p_max_nt_DJ_loop_Dr + log_p_max_nt_VD_loop_Dl;
                                            
                                            log_max_pins_loop_dDr = log(model.PinsDJ(1 + length(insDJ_nseq_min_loopdDr))*model.PinsVD(1 + length(insVD_nseq_min_loopdDl)));
                                            
                                            
                                            % Calculate number of errors in J section, by leaving out errors that are in deleted region for this assignment.
                                            d_err_pos = D.error_positions(na,:,d); % error positions in J alignment
                                            d_errs_i = (d_err_pos >= D.align_position_left(na,d) + ndDl1) & d_err_pos <=(D.align_position_right(na,d) -ndDr1);
                                            d_errs = sum(d_errs_i);
                                            
                                        else
                                            
                                            niVD_DJ_min = niVD_DJ_total;
                                            %niVD_DJ_total
                                            log_p_max_nt_VD_DJ_na_loop_Dr = niVD_DJ_min*log_max_model_p_nt;
                                            
                                            log_max_pins_loop_dDr = log_max_model_pins;
                                            
                                            ndDr1 = 0;
                                            D_align_length = 0;
                                            d_errs = 0;
                                            d_errs_i = [];
                                            d_err_pos = [];
                                        end
                                        
                                        
                                        
                                        % errors
                                        if ndDr1 < 0
                                            d_err_excess_pos_right = D.excess_error_positions_right(na,:,d); % error positions in extended J alignment for 'negative' deletions
                                            % Case of 'negative' deletions
                                            d_ex_errs_right_i = d_err_excess_pos_right > 0 & d_err_excess_pos_right <= (D.align_position_right(na,d) - ndDr1);
                                            d_ex_errs_right = sum(d_ex_errs_right_i);
                                            if d_ex_errs_right > 1
                                                continue;
                                            end
                                        else
                                            % Case of positive deletions
                                            d_ex_errs_right = 0;
                                            d_err_excess_pos_right = [];
                                            d_ex_errs_right_i = [];
                                        end
                                        
                                        nerrorsd = d_errs + d_ex_errs_left + d_ex_errs_right;
                                        
                                        if nerrorsd > D_max_error
                                            continue;
                                        end
                                        
                                        
                                        log_highest_probability_GIVEN_current_Dr_deletions = -1000;
                                        
                                        
                                        nerrors = nerrorsv + nerrorsj + nerrorsd;
                                        % Cost for errors is (error_rate/3)^(n_errors) because we know the mistaken nucleotides.
                                        if nerrors==0
                                            log_perr = 0;
                                        else
                                            log_perr = nerrors*log_Rerror_per_sequenced_nucleotide_divided_by_3;
                                        end
                                        
                                        test_proba = log_probabase + log_perr + log_max_pcutVDJ_loop_pJ + log_max_pins_loop_dDr + log_p_max_nt_VD_DJ_na_loop_Dr;
%                                         if(test_proba > log_highest_probability_GIVEN_current_D_allele), log_highest_probability_GIVEN_current_D_allele = test_proba; end;
%                                         if(test_proba > log_highest_probability_GIVEN_current_V_allele), log_highest_probability_GIVEN_current_V_allele = test_proba; end;
%                                         if(test_proba > log_highest_probability_GIVEN_current_J_allele), log_highest_probability_GIVEN_current_J_allele = test_proba; end;
%                                         if(test_proba > log_highest_probability_GIVEN_current_V_deletions), log_highest_probability_GIVEN_current_V_deletions = test_proba; end;
%                                         if(test_proba > log_highest_probability_GIVEN_current_J_deletions), log_highest_probability_GIVEN_current_J_deletions = test_proba; end;
%                                         if(test_proba > log_highest_probability_GIVEN_current_Dl_deletions), log_highest_probability_GIVEN_current_Dl_deletions = test_proba; end;
%                                         if(test_proba > log_highest_probability_GIVEN_current_Dr_deletions), log_highest_probability_GIVEN_current_Dr_deletions = test_proba; end;
                                        %if (~do_zeroD && (test_proba < log_probability_threshold_factor + log_highest_probability)) || (test_proba < log_probability_hopeless_threshold_factor + log_highest_probability)
                                        %if (test_proba < log_probability_threshold_factor + log_highest_probability)
                                        
                                        if ((test_proba < log_probability_threshold_factor + log_highest_probability)) || (test_proba < log_probability_hopeless_threshold_factor + log_highest_probability)
                                            %disp('skip 737')
%                                             test_proba
%                                             log_probabase
%                                             log_perr
%                                             log_max_pcutVDJ_loop_pJ
%                                             log_max_pins_loop_dDr
%                                             log_p_max_nt_VD_DJ_na_loop_Dr
%                                             log_highest_probability
%                                             niVD_DJ_min
%                                             ndJ1
%                                             ndV1
                                            skips = skips + 1;
                                            continue;
                                        end
                                        
                                        
                                        
                                        if np_start_from_max
                                            np_step = -1;
                                            npDl_start = npDl_max;
                                            npDl_end = 0;
                                            npDr_start = npDr_max;
                                            npDr_end = 0;
                                        else
                                            np_step = 1;
                                            npDl_start = 0;
                                            npDl_end = npDl_max;
                                            npDr_start = 0;
                                            npDr_end = npDr_max;
                                            
                                        end
                                        
                                        
                                        

                                        
                                        % Loop over half-length of D left palindrome
                                        for npDl=npDl_start: np_step : npDl_end
                                            if skips > max_skips
                                                assigns.n_assignments = in;
                                                assigns.skips = skips;
                                                return;
                                            end
                                            
                                            
                                            if ~zeroD
                                                D_left_start = D.align_position_left(na,d) + ndDl1 - npDl;
                                                % V with palindrome overlaps D.
                                                if V_end >= D_left_start
                                                    continue;
                                                end
                                                
                                                insVD_nseq_min_looppDl = nseq( V_end + 1 - READ_LENGTH_CORRECTION :  D.align_position_left(na,d) + ndDl1 - npDl - 1 - READ_LENGTH_CORRECTION);
                                                if length(insVD_nseq_min_looppDl) > model.max_insertions
                                                    continue;
                                                end
                                                min_mononucleotideVD = [sum(insVD_nseq_min_looppDl==1); sum(insVD_nseq_min_looppDl==2); sum(insVD_nseq_min_looppDl==3); sum(insVD_nseq_min_looppDl==4)];
                                                
                                                log_p_max_nt_VD_loop_pDl = min_mononucleotideVD(1)*log_max_model_p_nt_VD(1) + min_mononucleotideVD(2)*log_max_model_p_nt_VD(2) + min_mononucleotideVD(3)*log_max_model_p_nt_VD(3) + min_mononucleotideVD(4)*log_max_model_p_nt_VD(4);
                                                
                                                
                                                
                                            else
                                                log_p_max_nt_VD_loop_pDl = 0;
                                                % For compilation purposes, this needs to be assigned a value.
                                                D_left_start = -1;
                                                
                                            end
                                            
                                            
                                            % Loop over half-length of D right palindrome
                                            for npDr= npDr_start: np_step : npDr_end
                                                if skips > max_skips
                                                    assigns.n_assignments = in;
                                                    assigns.skips = skips;
                                                    return;
                                                end
                                                
                                                
                                                
                                                % if deletions is > 0 and palindromes > 0, skip.
                                                if ~zeroD
                                                    D_right_end = D.align_position_right(na,d) -ndDr1 + npDr;
                                                    % J with palindrome overlaps D.
                                                    if J_start <= D_right_end
                                                        continue;
                                                    end
                                                    
                                                    insDJ_nseq_min_looppDr = nseq( D.align_position_right(na,d) - ndDr1 + npDr + 1 - READ_LENGTH_CORRECTION : J_start - 1 - READ_LENGTH_CORRECTION);
                                                    if length(insDJ_nseq_min_looppDr) > model.max_insertions
                                                        continue;
                                                    end
                                                    min_mononucleotideDJ = [sum(insDJ_nseq_min_looppDr==1); sum(insDJ_nseq_min_looppDr==2); sum(insDJ_nseq_min_looppDr==3); sum(insDJ_nseq_min_looppDr==4)];
                                                    log_p_max_nt_DJ_loop_pDr = min_mononucleotideDJ(1)*log_max_model_p_nt_DJ(1) + min_mononucleotideDJ(2)*log_max_model_p_nt_DJ(2) + min_mononucleotideDJ(3)*log_max_model_p_nt_DJ(3) + min_mononucleotideDJ(4)*log_max_model_p_nt_DJ(4);
                                                    
                                                    
                                                    log_p_max_nt_VD_DJ_na_loop_pDr = log_p_max_nt_VD_loop_pDl + log_p_max_nt_DJ_loop_pDr;
                                                    
                                                    
                                                else
                                                    niVD_DJ_min = niVD_DJ_total;
                                                    log_p_max_nt_VD_DJ_na_loop_pDr = niVD_DJ_min*log_max_model_p_nt;
                                                    
                                                    % For compilation purposes, this needs to be assigned a value.
                                                    D_right_end = -1;
                                                end
                                                
                                                
                                                if zeroD
                                                    % For zero D, we convolute over VD and DJ insertions.
                                                    if niVD_DJ_total > 2*model.max_insertions
                                                        continue;
                                                    end
                                                    
                                                    if do_smoothing
                                                        niVD_min = max(0, niVD_DJ_total - model.max_insertions);
                                                        niVD_max = min(model.max_insertions, niVD_DJ_total);
                                                        
                                                        if niVD_min > model.max_insertions || (niVD_DJ_total - niVD_max) > model.max_insertions
                                                            continue;
                                                        end
                                                        
                                                    else
                                                        % For strawman, we just split in two.
                                                        niVD_min = floor(niVD_DJ_total/2);
                                                        niVD_max = niVD_min;
                                                        if niVD_min > model.max_insertions || (niVD_DJ_total - niVD_max) > model.max_insertions
                                                            continue;
                                                        end
                                                        
                                                    end
                                                    
                                                    
                                                    log_max_pins_loop_pDr = log(max( model.PinsVD(1+niVD_min:1+niVD_max).*model.PinsDJ((1+niVD_DJ_total-niVD_min):-1:(1+niVD_DJ_total-niVD_max)) ));
                                                    
                                                    % For compilation purposes, this needs to be assigned a value.
                                                    niVD_nonzero = -1;
                                                    niDJ_nonzero = -1;
                                                else
                                                    % For nonzero D, ndDr is fixed, so are niVD and niDJ.
                                                    niVD_nonzero = D_left_start - V_end - 1;
                                                    niDJ_nonzero = J_start - D_right_end - 1;
                                                    niVD_min = niVD_nonzero;
                                                    niVD_max = niVD_nonzero;
                                                    
                                                    if niVD_nonzero > model.max_insertions || niDJ_nonzero > model.max_insertions
                                                        continue;
                                                    end
                                                    
                                                    log_max_pins_loop_pDr = log(model.PinsDJ(1 + niDJ_nonzero)*model.PinsVD(1+niVD_nonzero));
                                                end
                                                
                                                % cut variables.
                                                ncutDl = ndDl - npDl;
                                                ncutDr = ndDr - npDr;
                                                % D deletions probability
                                                PcutDlDr = model.PcutDlcutDr_given_D( -model.min_D_cut + 1 + ncutDl, -model.min_D_cut + 1 + ncutDr, d_g);
                                                
                                                if PcutDlDr==0
                                                    %disp('skip 887')
                                                    skips = skips + 1;
                                                    continue;
                                                end
                                                
                                                
                                                log_PcutVDJ = log(PcutDlDr) + log_PcutV + log_PcutJ;
                                                
                                                test_proba = log_probabase + log_perr + log_PcutVDJ + log_max_pins_loop_pDr + log_p_max_nt_VD_DJ_na_loop_pDr;
%                                                 if(test_proba > log_highest_probability_GIVEN_current_D_allele), log_highest_probability_GIVEN_current_D_allele = test_proba; end;
%                                                 if(test_proba > log_highest_probability_GIVEN_current_V_allele), log_highest_probability_GIVEN_current_V_allele = test_proba; end;
%                                                 if(test_proba > log_highest_probability_GIVEN_current_J_allele), log_highest_probability_GIVEN_current_J_allele = test_proba; end;
%                                                 if(test_proba > log_highest_probability_GIVEN_current_V_deletions), log_highest_probability_GIVEN_current_V_deletions = test_proba; end;
%                                                 if(test_proba > log_highest_probability_GIVEN_current_J_deletions), log_highest_probability_GIVEN_current_J_deletions = test_proba; end;
%                                                 if(test_proba > log_highest_probability_GIVEN_current_Dl_deletions), log_highest_probability_GIVEN_current_Dl_deletions = test_proba; end;
%                                                 if(test_proba > log_highest_probability_GIVEN_current_Dr_deletions), log_highest_probability_GIVEN_current_Dr_deletions = test_proba; end;
                                                
                                                if ((test_proba < log_probability_threshold_factor + log_highest_probability)) || (test_proba < log_probability_hopeless_threshold_factor + log_highest_probability)
%                                                     if zeroD==1, %disp('!!!!!');
%                                                         log_probabase
%                                                         log_max_pins_loop_pDr
%                                                         log_p_max_nt_VD_DJ_na_loop_pDr
%                                                         log_PcutVDJ
%                                                         log_highest_probability
%                                                         log_probability_threshold_factor
%                                                         ncutDl
%                                                         ncutDr
%                                                         d_g
%                                                         niVD_min
%                                                         niVD_max
%                                                         niVD_DJ_total
%                                                         niVD_DJ_min
%                                                         ncutV
%                                                         ncutJ
%                                                         
%                                                         
%                                                     end
                                                    %disp('skip 923')
                                                    skips = skips + 1;
                                                    
                                                    continue;
                                                end
                                                
                                                
                                                
                                                % Loop over niVD for zero D case.
                                                % For nonzero D, this is single instance loop.
                                                %niVD_max
                                                for niVD=niVD_min:niVD_max
                                                 %   niVD
                                                    if zeroD
                                                        niDJ = niVD_DJ_total - niVD;
                                                    else
                                                        niDJ = niDJ_nonzero;
                                                    end
                                                    
                                                    % Check other constraints
                                                    % For instance, insertions can become negative with 'negative' deletions and palindromes together.
                                                    if niDJ > model.max_insertions || niVD  > model.max_insertions
                                                        continue;
                                                    end
                                                    
                                                    % nucleotide frequencies in insertion sequences (not including palindromes).
                                                    insVD_nseq = nseq( (V_end + 1 - READ_LENGTH_CORRECTION) :  (V_end  + niVD - READ_LENGTH_CORRECTION)); %% check this.
                                                    insDJ_nseq = nseq(  (J_start - niDJ - READ_LENGTH_CORRECTION)  :  (J_start - 1 - READ_LENGTH_CORRECTION) );
                                                    
                                                    insVD_nseq_5prime_shift = nseq( (V_end - READ_LENGTH_CORRECTION ) :  (V_end + niVD - 1 - READ_LENGTH_CORRECTION) );
                                                    insDJ_nseq_3prime_shift = nseq( (J_start - niDJ + 1 - READ_LENGTH_CORRECTION ) :  (J_start - READ_LENGTH_CORRECTION) );
                                                    
                                                    mononucleotideVD = [sum(insVD_nseq==1); sum(insVD_nseq==2); sum(insVD_nseq==3); sum(insVD_nseq==4)];
                                                    mononucleotideDJ = [sum(insDJ_nseq==1); sum(insDJ_nseq==2); sum(insDJ_nseq==3); sum(insDJ_nseq==4)];
                                                    
                                                    nucleotideVD = zeros(4,4);
                                                    for nt_VD = 1:niVD
                                                        nucleotideVD(insVD_nseq(nt_VD), insVD_nseq_5prime_shift(nt_VD)) = nucleotideVD(insVD_nseq(nt_VD), insVD_nseq_5prime_shift(nt_VD)) + 1;
                                                    end
                                                    
                                                    nucleotideDJ = zeros(4,4);
                                                    for nt_DJ = 1:niDJ
                                                        nucleotideDJ(insDJ_nseq(nt_DJ), insDJ_nseq_3prime_shift(nt_DJ)) = nucleotideDJ(insDJ_nseq(nt_DJ), insDJ_nseq_3prime_shift(nt_DJ)) + 1;
                                                    end
                                                    
                                                    % Compute final probability of assignment by multiplying base probability by factors for
                                                    % palindromes, deletions and the number of errors.
                                                    % multiply probability by factor for insertions of given lengths and sequence specific factors for
                                                    % nucleotide bias
                                                    log_pins = log(model.PinsVD(1+niVD)*model.PinsDJ(1 + niDJ));
                                                    
                                                    log_pntbias_VD = log_RnucleotideVD_per_nucleotideVD_5prime*nucleotideVD(:);
                                                    log_pntbias_DJ = log_RnucleotideDJ_per_nucleotideDJ_3prime*nucleotideDJ(:);                                                    
                                                    log_proba = log_probabase + log_perr + log_PcutVDJ + log_pins + log_pntbias_DJ + log_pntbias_VD;
                                                    
                                                    
                                                    % Update all highest probability so far values, if necessary

                                                    
                                                    if(log_proba > log_highest_probability_GIVEN_current_D_allele), log_highest_probability_GIVEN_current_D_allele = log_proba; end;
                                                    if(log_proba > log_highest_probability_GIVEN_current_V_allele), log_highest_probability_GIVEN_current_V_allele = log_proba; end;
                                                    if(log_proba > log_highest_probability_GIVEN_current_J_allele), log_highest_probability_GIVEN_current_J_allele = log_proba; end;
                                                    if(log_proba > log_highest_probability_GIVEN_current_V_deletions), log_highest_probability_GIVEN_current_V_deletions = log_proba; end;
                                                    if(log_proba > log_highest_probability_GIVEN_current_J_deletions), log_highest_probability_GIVEN_current_J_deletions = log_proba; end;
                                                    if(log_proba > log_highest_probability_GIVEN_current_Dl_deletions), log_highest_probability_GIVEN_current_Dl_deletions = log_proba; end;
                                                    if(log_proba > log_highest_probability_GIVEN_current_Dr_deletions), log_highest_probability_GIVEN_current_Dr_deletions = log_proba; end;
                                                    if D_align_length > best_D_align_length, best_D_align_length = D_align_length; end
                                                    
                                                    % if D length isvery short (<3) I let it count low prob events because the total prob of D short
                                                    % is spread out over many choices of D deletions. So if I skip them,
                                                    %if (~do_zeroD && (log_proba < log_probability_threshold_factor + log_highest_probability)) || (log_proba < log_probability_hopeless_threshold_factor + log_highest_probability)
                                                    if ((log_proba < log_probability_threshold_factor + log_highest_probability)) || (log_proba < log_probability_hopeless_threshold_factor + log_highest_probability)
                                                        %disp('skip 996')
                                                        skips = skips + 1;
                                                        continue;
                                                    end

                                                    
                                                    nucleotideVD_5prime = repmat(sum(nucleotideVD,1),4,1);
                                                    nucleotideDJ_3prime = repmat(sum(nucleotideDJ,1),4,1);
                                                    
                                                    
                                                    VD_left_edge_dinucleotide = uint8([0,0]);
                                                    VD_right_edge_dinucleotide = uint8([0,0]);
                                                    if niVD > 0
                                                        VD_left_edge_dinucleotide = [nseq(V_end-READ_LENGTH_CORRECTION), nseq(V_end+1-READ_LENGTH_CORRECTION)];
                                                        if ~zeroD
                                                            VD_right_edge_dinucleotide = [nseq(D_left_start-1-READ_LENGTH_CORRECTION), nseq(D_left_start-READ_LENGTH_CORRECTION)];
                                                        end
                                                    end
                                                    
                                                    
                                                    DJ_left_edge_dinucleotide = uint8([0,0]);
                                                    DJ_right_edge_dinucleotide = uint8([0,0]);
                                                    if niDJ > 0
                                                        DJ_right_edge_dinucleotide = [nseq(J_start-1-READ_LENGTH_CORRECTION), nseq(J_start-READ_LENGTH_CORRECTION)];
                                                        if ~zeroD
                                                            DJ_left_edge_dinucleotide = [nseq(D_right_end-READ_LENGTH_CORRECTION), nseq(D_right_end+1-READ_LENGTH_CORRECTION)];
                                                        end
                                                    end
                                                    
                                                    
                                                    trinucleotideVD = zeros(4,4,4);
                                                    trinucleotideDJ = zeros(4,4,4);
                                                    
                                                    for nt_VD = 1:(niVD-2)
                                                        trinucleotideVD(insVD_nseq(nt_VD), insVD_nseq(nt_VD+1), insVD_nseq(nt_VD+2)) = trinucleotideVD(insVD_nseq(nt_VD), insVD_nseq(nt_VD+1), insVD_nseq(nt_VD+2)) + 1;
                                                    end
                                                    
                                                    for nt_DJ = 1:(niDJ-2)
                                                        trinucleotideDJ(insDJ_nseq(nt_DJ), insDJ_nseq(nt_DJ+1), insDJ_nseq(nt_DJ+2)) = trinucleotideDJ(insDJ_nseq(nt_DJ), insDJ_nseq(nt_DJ+1), insDJ_nseq(nt_DJ+2)) + 1;
                                                    end
                                                    
                                                    
                                                    
                                                    error_vs_position = zeros(model.read_length,1); % to store positions with errors for this assignment
                                                    
                                                    % Total nucleotides that are assigned to either V, D or J in this assignment
                                                    % (ignoring high error region to the left)
                                                    
                                                    if ~zeroD
                                                        genic_length = (V.align_length(v) + J.align_length(j) + npV + npJ -ndV1 -ndJ1 + D.align_length(na,d) + npDl + npDr - ndDl1 - ndDr1 - high_error_region - READ_LENGTH_CORRECTION);
                                                    else
                                                        genic_length = (V.align_length(v) + J.align_length(j) +npV + npJ - ndV1 - ndJ1 - high_error_region - READ_LENGTH_CORRECTION);
                                                    end
                                                    
                                                    % If genic length is too short, skip. This can happen with the deletions variables increasing.
%                                                     if genic_length < 10
%                                                         skips = skips + 1; continue;
%                                                     end
                                                    
                                                    % Store all positions that have been assigned to either V, D or J.
                                                    % This is used to estimate error rate vs position by dividing totat errors in a position by coverage.
                                                    coverage = zeros(length(seq),1);
                                                    coverage( 1: (V.align_length(v) -ndV1 + npV - READ_LENGTH_CORRECTION )) = 1; % Set V gene segment to 1.
                                                    coverage( (J.align_position(1,j) + ndJ1 - npJ - READ_LENGTH_CORRECTION) : end) = 1; % Set J gene segment to 1.
                                                    if ~zeroD
                                                        coverage( (D.align_position_left(na,d) + ndDl1 -npDl - READ_LENGTH_CORRECTION) : (D.align_position_right(na,d) - ndDr1 + npDr - READ_LENGTH_CORRECTION) ) = 1; % Set D gene segment to 1.
                                                    end
                                                    
                                                    % store potential pal that is consistent with assignment.
                                                    npV_potential_max = min(npV_potential_max, niVD + npV);
                                                    npJ_potential_max = min(npJ_potential_max, niDJ + npJ);
                                                    
                                                    if ~zeroD
                                                        npDl_potential_max = min(npDl_potential_max, niVD + npDl);
                                                        npDr_potential_max = min(npDr_potential_max, niDJ + npDr);
                                                    else
                                                        npDl_potential_max = 0;
                                                        npDr_potential_max = 0;
                                                    end
                                                    
                                                    % Increment valid assignment number
                                                    in=in+1;
                                                    n_assignments_v_gene = n_assignments_v_gene + 1;
                                                    n_assignments_d_gene = n_assignments_d_gene + 1;
                                                    n_assignments_j_gene = n_assignments_j_gene + 1;
                                                    if(log_proba > log_highest_probability), log_highest_probability = log_proba; assigns.max_proba_index = in; end;                                                    
                                                                                                        
                                                    v_err_pos_rel = zeros(max_V_depth,1);
                                                    j_err_pos_rel = zeros(max_J_depth,1);
                                                    if ~no_error
                                                        
                                                        if nerrorsv > 0
                                                            if ndV1 < 0
                                                                
                                                                % 'Negative' V deletions case
                                                                if v_ex_errs > 0
                                                                    v_err_excess_pos_ok = v_err_excess_pos(v_ex_errs_i);
                                                                    error_vs_position(v_err_excess_pos_ok)=1; % errors in extended V alignment, from 'negative' deletions
                                                                end
                                                                v_err_pos_ok = v_err_pos(v_err_pos>0);
                                                                error_vs_position(v_err_pos_ok)=1; % errors in actual V alignment
                                                            else
                                                                v_err_pos_ok = v_err_pos(v_err_pos > 0 & v_err_pos <= (V.align_length(v) - ndV1) );
                                                                error_vs_position(v_err_pos_ok)=1;
                                                                
                                                            end
                                                            v_err_pos_rel(ndV + V_align_length:-1:1+ndV) = error_vs_position(1:V_align_length);
                                                            
                                                        end
                                                        
                                                        % Set positions of errors in J to 1, accounting for deletions
                                                        
                                                        if nerrorsj > 0
                                                            if ndJ1 < 0
                                                                
                                                                % 'Negative' J deletions case
                                                                if j_ex_errs > 0
                                                                    j_err_excess_pos_ok = j_err_excess_pos( j_ex_errs_i );
                                                                    error_vs_position(j_err_excess_pos_ok) = 1; % errors in extended J alignment, from 'negative' deletions
                                                                end
                                                                error_vs_position(j_err_pos(j_err_pos>0)) = 1; % errors in actual J alignment
                                                            else
                                                                error_vs_position( j_err_pos(j_err_pos >= (J.align_position(1,j)+ndJ1) ) ) = 1;
                                                            end
                                                            
                                                            j_err_pos_rel(ndJ +1 : ndJ + J_align_length) = error_vs_position(end - J_align_length + 1:end);
                                                            
                                                            
                                                        end
                                                        % Set positions of errors in D to 1. All of these come from 'negative' deletions since D alignments do not
                                                        % allow errors.
                                                        
                                                        if ~zeroD && nerrorsd > 0
                                                            
                                                            if ndDl1 < 0 && d_ex_errs_left > 0
                                                                error_vs_position(d_err_excess_pos_left(d_ex_errs_left_i)) = 1;
                                                            end
                                                            
                                                            if ndDr1 < 0 && d_ex_errs_right > 0
                                                                error_vs_position(d_err_excess_pos_right(d_ex_errs_right_i)) = 1;
                                                            end
                                                            
                                                            if d_errs > 0
                                                                error_vs_position(d_err_pos(d_errs_i)) = 1;
                                                            end
                                                        end
                                                        
                                                    end
                                                    
                                                    
                                                    % Store everything in assigns.
                                                    % Each variable in counter that begins with an 'nP' or 'nM' must be
                                                    % found in assigns without the 'nP' or 'nM'.
                                                    
                                                    %%% The following are model variables
                                                    %disp(['Done ' num2str(na)])
                                                    assigns.V(in,:) = v_g;
                                                    assigns.DJ(in,:)=[d_g, j_g]; % gene choices
                                                    assigns.Vallele_given_gene(in,:) = [genV(v_a).allele, v_g]; % allele choice given gene
                                                    assigns.Dallele_given_gene(in,:) = [genD(d_a).allele, d_g];
                                                    assigns.cutV_given_V(in,:) = [ 1 - model.min_V_cut + ncutV, v_g];
                                                    assigns.cutJ_given_J(in,:) = [1 - model.min_J_cut + ncutJ, j_g];
                                                    assigns.insVD(in,:)=1 + niVD; % insertions
                                                    assigns.insDJ(in,:)=1 + niDJ;

						    % The following idiocy is needed to avoid a bug in the matlab C compile it seems...
                                                    assigns.nucleotideVD(in,:,1) = nucleotideVD(:,1);
                                                    assigns.nucleotideVD(in,:,2) = nucleotideVD(:,2);
                                                    assigns.nucleotideVD(in,:,3) = nucleotideVD(:,3);
                                                    assigns.nucleotideVD(in,:,4) = nucleotideVD(:,4);
                                                    assigns.nucleotideVD_5prime(in,:,1) = nucleotideVD_5prime(:,1);
                                                    assigns.nucleotideVD_5prime(in,:,2) = nucleotideVD_5prime(:,2);
                                                    assigns.nucleotideVD_5prime(in,:,3) = nucleotideVD_5prime(:,3);
                                                    assigns.nucleotideVD_5prime(in,:,4) = nucleotideVD_5prime(:,4);
                                                    assigns.nucleotideDJ(in,:,1) = nucleotideDJ(:,1);
                                                    assigns.nucleotideDJ(in,:,2) = nucleotideDJ(:,2);
                                                    assigns.nucleotideDJ(in,:,3) = nucleotideDJ(:,3);
                                                    assigns.nucleotideDJ(in,:,4) = nucleotideDJ(:,4);
                                                    assigns.nucleotideDJ_3prime(in,:,1) = nucleotideDJ_3prime(:,1);
                                                    assigns.nucleotideDJ_3prime(in,:,2) = nucleotideDJ_3prime(:,2);
                                                    assigns.nucleotideDJ_3prime(in,:,3) = nucleotideDJ_3prime(:,3);
                                                    assigns.nucleotideDJ_3prime(in,:,4) = nucleotideDJ_3prime(:,4);
                                                    
                                                    assigns.error(in,:) = nerrors; % numerator for error rate estimate
                                                    assigns.sequenced_nucleotide(in,:) = genic_length; % denominator for error rate estimate
                                                    
                                                    
                                                    % For tracking
                                                    
                                                    assigns.v_alignment_index(in,:) = v;
                                                    assigns.j_alignment_index(in,:) = j;
                                                    assigns.d_alignment_index(in,:) = [na, d];
                                                    
                                                    
                                                    
                                                    % numbers of A,C,G and T in insertions: numerator for probabilities of A,C,G,T insertions
                                                    assigns.mononucleotideVD(in,:) = mononucleotideVD;
                                                    assigns.insertionVD(in,:) = [ niVD , niVD, niVD, niVD]; % denominator for probabilities of A,C,G,T insertions.
                                                    assigns.mononucleotideDJ(in,:) = mononucleotideDJ;
                                                    assigns.insertionDJ(in,:) = [niDJ , niDJ, niDJ, niDJ]; % denominator for probabilities of A,C,G,T insertions.
                                                    
                                                    assigns.VD_left_edge_dinucleotide(in,:) = VD_left_edge_dinucleotide;
                                                    assigns.VD_right_edge_dinucleotide(in,:) = VD_right_edge_dinucleotide;
                                                    
                                                    assigns.DJ_left_edge_dinucleotide(in,:) = DJ_left_edge_dinucleotide;
                                                    assigns.DJ_right_edge_dinucleotide(in,:) = DJ_right_edge_dinucleotide;
                                                    
                                                    assigns.trinucleotideVD(in,:,:,:) = trinucleotideVD;
                                                    assigns.trinucleotideDJ(in,:,:,:) = trinucleotideDJ;
                                                    
                                                    assigns.VV_err_pos(in,v_g,:) = v_err_pos_rel;
                                                    assigns.VV_align_length(in,:) = [v_g, 1 + V_align_length];
                                                    
                                                    assigns.JJ_err_pos(in,j_g,:) = j_err_pos_rel;
                                                    assigns.JJ_align_length(in,:) = [j_g, 1 + J_align_length];
                                                    
                                                    
                                                    assigns.pVmax_delV_V(in,:) = [1 + npV_potential_max, 1+ ndV, v_g];
                                                    assigns.pJmax_delJ_J(in,:) = [1 + npJ_potential_max, 1+ ndJ, j_g];
                                                    
                                                    
                                                    assigns.pDlmax_delDl_D(in,:) = [1 + npDl_potential_max, 1+ ndDl, d_g];
                                                    assigns.pDrmax_delDr_D(in,:) = [1 + npDr_potential_max, 1+ ndDr, d_g];
                                                    
                                                    assigns.zeroD(in,:) = zeroD;
                                                    assigns.V_align_length(in,:) = [1 + V_align_length];
                                                    assigns.D_align_length(in,:) = [1 + D_align_length];
                                                    assigns.J_align_length(in,:) = [1 + J_align_length];
                                                    assigns.VDJ(in,:)=[v_g, d_g, j_g]; % gene choices
                                                    assigns.pVdelV(in,:)=[ 1 + npV , 1 + ndV]; % palindromes and deletions
                                                    assigns.pJdelJ(in,:)=[ 1 + npJ , 1 + ndJ];
                                                    assigns.delVinsVD(in,:) = [1 + ndV, 1 + niVD];
                                                    assigns.delVinsDJ(in,:) = [1 + ndV, 1 + niDJ];
                                                    assigns.delVdelJ(in,:) = [1 + ndV, 1 + ndJ];
                                                    assigns.delJinsVD(in,:) = [1 + ndJ, 1 + niVD];
                                                    assigns.delJinsDJ(in,:) = [1 + ndJ, 1 + niDJ];
                                                    
                                                    assigns.insVDinsDJ(in,:) = [1 + niVD, 1 + niDJ];
                                                    
                                                    assigns.insDJ_D_align_length(in,:) = [1 + niDJ, 1 + D_align_length];
                                                    assigns.insVD_D_align_length(in,:) = [1 + niVD, 1 + D_align_length];
                                                    
                                                    assigns.insDJ_J_align_length(in,:) = [1 + niDJ, 1 + J_align_length];
                                                    assigns.insVD_J_align_length(in,:) = [1 + niVD, 1 + J_align_length];
                                                    
                                                    assigns.insDJ_V_align_length(in,:) = [1 + niDJ, 1 + V_align_length];
                                                    assigns.insVD_V_align_length(in,:) = [1 + niVD, 1 + V_align_length];
                                                    
                                                    assigns.Dallele_D_align_length(in,:) = [d_a, 1 + D_align_length];
                                                    
                                                    
                                                    
                                                    assigns.VdelV(in,:) = [v_g, 1+ ndV ] ;
                                                    assigns.JdelJ(in,:) =  [j_g, 1 + ndJ];
                                                    assigns.VinsVD(in,:) = [v_g, 1 + niVD];
                                                    assigns.DinsVD(in,:) = [d_g, 1 + niVD];
                                                    assigns.DinsDJ(in,:) = [d_g, 1 + niDJ];
                                                    assigns.JinsDJ(in,:) = [j_g, 1 + niDJ];
                                                    assigns.VdelJ(in,:) = [v_g, 1 + ndJ];
                                                    assigns.JdelV(in,:) = [j_g, 1 + ndV];
                                                    
                                                    assigns.DdelV(in,:) = [d_g, 1 + ndV];
                                                    assigns.DdelJ(in,:) = [d_g, 1 + ndJ];
                                                    assigns.VinsDJ(in,:) = [v_g, 1 + niDJ];
                                                    assigns.JinsVD(in,:) = [j_g, 1 + niVD];
                                                    
                                                    assigns.pVinsVD(in,:) = [1 + npV, 1 + niVD];
                                                    assigns.pVinsDJ(in,:) = [1 + npV, 1 + niDJ];
                                                    assigns.pVdelJ(in,:) = [1 + npV, 1 + ndJ];
                                                    assigns.VpV(in,:) = [v_g, 1 + npV];
                                                    assigns.JpV(in,:) = [j_g, 1 + npV];
                                                    assigns.DpV(in,:) = [d_g, 1 + npV];
                                                    
                                                    assigns.pJinsVD(in,:) = [1 + npJ, 1 + niVD];
                                                    assigns.pJinsDJ(in,:) = [1 + npJ, 1 + niDJ];
                                                    assigns.pJdelV(in,:) = [1 + npJ, 1 + ndV];
                                                    assigns.VpJ(in,:) = [v_g, 1 + npJ];
                                                    assigns.JpJ(in,:) = [j_g, 1 + npJ];
                                                    assigns.DpJ(in,:) = [d_g, 1 + npJ];
                                                    assigns.pVpJ(in,:)  = [1 + npV, 1+ npJ];
                                                    assigns.VcutV(in,:) = [ v_g, 1 - model.min_V_cut + ncutV];
                                                    assigns.JcutJ(in,:) = [ j_g, 1 - model.min_J_cut + ncutJ];
                                                    assigns.DcutV(in,:) = [ d_g, 1 - model.min_V_cut + ncutV];
                                                    assigns.DcutJ(in,:) = [ d_g, 1 - model.min_J_cut + ncutJ];
                                                    assigns.VcutJ(in,:) = [ v_g, 1 - model.min_J_cut + ncutJ];
                                                    assigns.JcutV(in,:) = [ j_g, 1 - model.min_V_cut + ncutV];
                                                    assigns.insVDcutV(in,:) = [ 1 + niVD, 1 - model.min_V_cut + ncutV];
                                                    assigns.insDJcutV(in,:) = [ 1 + niDJ, 1 - model.min_V_cut + ncutV];
                                                    assigns.insVDcutJ(in,:) = [ 1 + niVD, 1 - model.min_J_cut + ncutJ];
                                                    assigns.insDJcutJ(in,:) = [ 1 + niDJ, 1 - model.min_J_cut + ncutJ];
                                                    
                                                    assigns.cutDlcutDr_given_D(in,:) = [ 1 - model.min_D_cut + ncutDl, 1 - model.min_D_cut + ncutDr, d_g];
                                                    
                                                    % All possible pairwise joint distributions, just for tracking.
                                                    
                                                    assigns.pDldelDl(in,:)=[ 1 + npDl , 1 + ndDl];
                                                    assigns.pDrdelDr(in,:)=[ 1 + npDr , 1 + ndDr];
                                                    assigns.delVdelDl(in,:) = [1 + ndV, 1 + ndDl];
                                                    assigns.delVdelDr(in,:) = [1 + ndV, 1 + ndDr];
                                                    
                                                    assigns.delJdelDl(in,:) = [1 + ndJ, 1 + ndDl];
                                                    assigns.delJdelDr(in,:) = [1 + ndJ, 1 + ndDr];
                                                    assigns.delDlinsVD(in,:) = [1 + ndDl, 1 + niVD];
                                                    assigns.delDlinsDJ(in,:) = [1 + ndDl, 1 + niDJ];
                                                    assigns.delDldelDr(in,:) = [1 + ndDl, 1 + ndDr];
                                                    assigns.delDrinsVD(in,:) = [1 + ndDr, 1 + niVD];
                                                    assigns.delDrinsDJ(in,:) = [1 + ndDr, 1 + niDJ];
                                                    assigns.DdelDl(in,:) = [d_g, 1 + ndDl];
                                                    assigns.DdelDr(in,:) = [d_g, 1 + ndDr ];
                                                    assigns.VdelDl(in,:) = [v_g, 1 + ndDl];
                                                    assigns.VdelDr(in,:) = [v_g, 1 + ndDr];
                                                    assigns.JdelDl(in,:) = [j_g, 1 + ndDl];
                                                    assigns.JdelDr(in,:) = [j_g, 1 + ndDr];
                                                    assigns.pVdelDl(in,:) = [1 + npV, 1 + ndDl];
                                                    assigns.pVdelDr(in,:) = [1 + npV, 1 + ndDr];
                                                    assigns.pJdelDl(in,:) = [1 + npJ, 1 + ndDl];
                                                    assigns.pJdelDr(in,:) = [1 + npJ, 1 + ndDr];
                                                    assigns.pDlinsVD(in,:) = [1 + npDl, 1+niVD];
                                                    assigns.pDlinsDJ(in,:) = [1 + npDl, 1+niDJ];
                                                    assigns.pDldelV(in,:) = [1 + npDl, 1+ndV];
                                                    assigns.pDldelJ(in,:)  = [1 + npDl, 1+ndJ];
                                                    assigns.pDldelDr(in,:) = [1 + npDl, 1+ndDr];
                                                    assigns.VpDl(in,:) = [v_g, 1+ npDl];
                                                    assigns.JpDl(in,:) = [j_g, 1+ npDl];
                                                    assigns.DpDl(in,:) = [d_g, 1+ npDl];
                                                    
                                                    assigns.pDrinsVD(in,:) = [1 + npDr, 1 + niVD];
                                                    assigns.pDrinsDJ(in,:) = [1 + npDr, 1 + niDJ];
                                                    assigns.pDrdelV(in,:)=  [1 + npDr, 1 + ndV];
                                                    assigns.pDrdelJ(in,:)=  [1 + npDr, 1 + ndJ];
                                                    assigns.pDrdelDl(in,:) = [1 + npDr, 1 + ndDl];
                                                    assigns.VpDr(in,:) = [v_g, 1 + npDr];
                                                    assigns.JpDr(in,:) = [j_g, 1 + npDr];
                                                    assigns.DpDr(in,:) = [d_g, 1 + npDr];
                                                    
                                                    assigns.pVpDl(in,:)  = [1 + npV, 1+ npDl];
                                                    assigns.pVpDr(in,:)  = [1 + npV, 1+ npDr];
                                                    assigns.pDlpDr(in,:) = [1 + npDl, 1+ npDr];
                                                    assigns.pDlpJ(in,:) = [1 + npDl, 1+ npJ];
                                                    assigns.pDrpJ(in,:) = [1 + npDr, 1+ npJ];
                                                    
                                                    assigns.DcutDl(in,:) = [ d_g, 1 - model.min_D_cut + ncutDl];
                                                    assigns.DcutDr(in,:) = [ d_g, 1 - model.min_D_cut + ncutDr];
                                                    assigns.VcutDl(in,:) = [ v_g, 1 - model.min_D_cut + ncutDl];
                                                    assigns.VcutDr(in,:) = [ v_g, 1 - model.min_D_cut + ncutDr];
                                                    assigns.JcutDl(in,:) = [ j_g, 1 - model.min_D_cut + ncutDl];
                                                    assigns.JcutDr(in,:) = [ j_g, 1 - model.min_D_cut + ncutDr];
                                                    assigns.insVDcutDl(in,:) = [1 + niVD,  1 - model.min_D_cut + ncutDl];
                                                    assigns.insDJcutDl(in,:) = [1 + niDJ,  1 - model.min_D_cut + ncutDl];
                                                    assigns.insVDcutDr(in,:) = [1 + niVD,  1 - model.min_D_cut + ncutDr];
                                                    assigns.insDJcutDr(in,:) = [1 + niDJ,  1 - model.min_D_cut + ncutDr];
                                                    
                                                    assigns.error_vs_position(in,:) = error_vs_position;
                                                    assigns.coverage(in,:) = coverage;
                                                    
                                                    %%% probability of this assignment
                                                    assigns.proba(in)=exp(log_proba + log_proba_Rerror_normalization);
                                                    %if in==1 || in==10
                                                       %%%log_proba = log_probabase + log_perr + log_PcutVDJ + log_pins + log_pntbias_DJ + log_pntbias_VD
                                                       %in
                                                       %log_proba
                                                       %log_probabase
                                                       %log_perr
                                                       %log_PcutVDJ
                                                       %log_pins
                                                       %log_pntbias_DJ
                                                       %log_pntbias_VD
                                                       
                                                    %end
                                                    if nerrors==0
                                                        assigns.event_probability(in) = exp(log_proba);
                                                    else
                                                        assigns.event_probability(in) = 0;
                                                    end
                                                    
                                                    assigns.generation_probability = assigns.generation_probability + assigns.event_probability(in);
                                                    assigns.likelihood = assigns.likelihood + assigns.proba(in);
                                                    
                                                    % If maximum number of valid assignments is reached, stop.
                                                    % If we are not smoothing, stop after 1 assignment
                                                    % Usually, we don't reach maximum number of assignments because we break out
                                                    % of all the loops due to the various thresholds below.
                                                    if in == model.max_assignments || (~do_smoothing && in==1)
                                                        assigns.n_assignments = in;
                                                        assigns.skips = skips;
                                                        return;
                                                    end
                                                    
                                                    if force_all_alleles && n_assignments_v_gene > model.max_assignments/length(V.alleles_all)
                                                        v_break_out = true;
                                                        break;
                                                    end
                                                    if force_all_alleles && n_assignments_j_gene > model.max_assignments/length(J.alleles_all)
                                                        j_break_out = true;
                                                        break;
                                                    end
                                                    if force_all_alleles && n_assignments_d_gene > model.max_assignments/length(D.allele_order)
                                                        d_break_out = true;
                                                        break;
                                                    end
                                                    
                                                end
                                                
                                            end
                                        end
                                        if( ~zeroD && log_highest_probability_GIVEN_current_Dr_deletions> -900 &&  ( ( (ndDr1 > 2 && d_errs_internal_max==0) && (log_highest_probability_GIVEN_current_Dr_deletions < log_highest_probability + log_probability_threshold_factor)) || (ndDr1 > 0 && d_errs_internal_max==0 && log_highest_probability_GIVEN_current_Dr_deletions < log_highest_probability + log_probability_hopeless_threshold_factor)) ), break; end;
                                    end
                                    if ( ~zeroD && log_highest_probability_GIVEN_current_Dl_deletions > -900 && ((ndDl1 > 2 && d_errs_internal_max==0) && (log_highest_probability_GIVEN_current_Dl_deletions < log_highest_probability + log_probability_threshold_factor) ) || ( ndDl1 > 0 && d_errs_internal_max==0 && log_highest_probability_GIVEN_current_Dl_deletions < log_highest_probability + log_probability_hopeless_threshold_factor)), break; end;
                                end
                                
                                if v_break_out || j_break_out || d_break_out
                                    break;
                                end
                                
                            end
                            
                            if v_break_out || j_break_out || d_break_out
                                break;
                            end
                            
                        end
                        
                        if v_break_out || j_break_out || d_break_out
                            break;
                        end
                        
                    end
                    
                    if log_highest_probability_GIVEN_current_J_deletions > -900 && ( ((ndJ1 > 1) && nerrorsj==0 && (log_highest_probability_GIVEN_current_J_deletions < log_highest_probability + log_probability_threshold_factor) ) || (ndJ1 > 0 && (log_highest_probability_GIVEN_current_J_deletions < log_highest_probability + log_probability_hopeless_threshold_factor))) , break; end;
                    if v_break_out || j_break_out || d_break_out
                        break;
                    end
                    
                end
                
                if log_highest_probability_GIVEN_current_V_deletions > -900 && (( (ndV1 > 1) && nerrorsv==0 && (log_highest_probability_GIVEN_current_V_deletions < log_highest_probability + log_probability_threshold_factor) ) || (ndV1 > 0 && (log_highest_probability_GIVEN_current_V_deletions < log_highest_probability + log_probability_hopeless_threshold_factor))),  break; end;
                if v_break_out || j_break_out || d_break_out
                    break;
                end
                
            end
            if v_break_out || j_break_out
                break;
            end
            
            
        end
        
        if v_break_out
            break;
        end
        
        if log_highest_probability_GIVEN_current_J_allele > -900 && (log_highest_probability_GIVEN_current_J_allele < log_highest_probability + log_probability_threshold_factor), break; end;
        
    end
    
    if log_highest_probability_GIVEN_current_V_allele > -900 && (log_highest_probability_GIVEN_current_V_allele < log_highest_probability + log_probability_threshold_factor),  break; end;
    
end
% Store how many valid assignments were considered.
assigns.n_assignments = in;
assigns.skips = skips;
end
