def compute_Pi_R(self, CDR3_seq, Pi_JinsDJ_given_D):
        """Compute Pi_R.
    
        This function returns the Pi array from the model factors of the D and J 
        genomic contributions, P(D, J)*P(delJ|J)P(delDl, delDr |D) and
        the DJ (N2) insertions,
        first_nt_bias_insDJ(n_1)PinsDJ(\ell_{DJ})\prod_{i=2}^{\ell_{DJ}}Rdj(n_i|n_{i-1}).
        This corresponds to \sum_D {D^{x_2}}_{x_3}{N^{x_3}}_{x_4}J(D)^{x_4}.
        
        For clarity in parsing the algorithm implementation, we include which 
        instance attributes are used in the method as 'parameters.'
    
        Parameters
        ----------
        CDR3_seq : str
            CDR3 sequence composed of 'amino acids' (single character symbols
            each corresponding to a collection of codons as given by codons_dict).
        Pi_JinsDJ_given_D : list
            List of (4, 3L) ndarrays corresponding to {N^{x_3}}_{x_4}J(D)^{x_4}.
            
        self.cutD_genomic_CDR3_segs : list of strings
            List of all the D genomic nucleotide sequences with the maximum number 
            of palindromic insertions appended on both ends.
        self.PD_given_J : ndarray
            Probability distribution of D conditioned on J, i.e. P(D|J).
        self.PD_nt_pos_vec : list of ndarrays
            For each D allele, format P(delDl, delDr|D) into the correct form 
            for a Pi array as if each position were the first in a codon.
        self.PD_2nd_nt_pos_per_aa_vec : list of dicts
            For each D allele, and each 'amino acid', format P(delDl, delDr|D) 
            for positions in the middle of a codon into the correct form for a 
            Pi array as if each position were the middle of a codon 
            corresponding to the 'amino acid'.
        self.min_delDl_given_DdelDr : list of lists
            minimum delDl for each delDr, D combination.
        self.max_delDl_given_DdelDr : list of lists
            maximum delDl for each delDr, D combination.
        self.PdelDldelDr_given_D : ndarray
            Joint probability distribution of the D deletions given the D allele,
            i.e. P(delDl, delDr |D)
        self.zeroD_given_D : list of floats
            The probability that a given D allele is fully deleted away.
        self.codons_dict : dict
            Dictionary, keyed by the allowed 'amino acid' symbols with the values 
            being lists of codons corresponding to the symbol.
        self.sub_codons_right : dict
            Dictionary of the 1 and 2 nucleotide suffixes (read from 5') for 
            each codon in an 'amino acid' grouping
    
        Returns
        -------
        Pi_L : ndarray
            (4, 3L) array corresponding to 
            \sum_D {D^{x_2}}_{x_3}{N^{x_3}}_{x_4}J(D)^{x_4}.
        
        """
        
        
        #Need to consider all D alignments from all possible positions and right deletions.
        nt2num = {'A': 0, 'C': 1, 'G': 2, 'T': 3}
        
        #n_aaseq = [aa_dict[aa] for aa in CDR3_seq]
        Pi_R = np.zeros((4, len(CDR3_seq)*3))
        min_pos = -len(CDR3_seq)*3
        
        num_dell_pos, num_delr_pos, num_D_genes = self.PdelDldelDr_given_D.shape
        
        for D_in, cutD_gen_seg  in enumerate(self.cutD_genomic_CDR3_segs):
            l_D_seg = len(cutD_gen_seg)
            
            #start position is first nt in a codon
            for init_pos in range(-1,-len(CDR3_seq)*3-1,-3):
                Pi_R[:, init_pos] += Pi_JinsDJ_given_D[D_in][:, init_pos]*self.zeroD_given_D[D_in]
                second_pos_dict = {'A': np.zeros(4), 'C': np.zeros(4), 'G': np.zeros(4), 'T': np.zeros(4)}
                codon_prefix_dict = {}
    
                for last_nt in 'ACGT':
                    for second_nt in 'ACGT':
                        codon_prefix_dict[last_nt + second_nt] = np.zeros(4)
                        #for first_nt in ['ACGT'[nt] for nt in range(4) if Pi_JinsDJ_given_D[D_in][nt, init_pos] > 0]:
                        for first_nt in 'ACGT':
                            if last_nt + second_nt + first_nt in self.codons_dict[CDR3_seq[init_pos/3]]: #possible allowed codon
                                second_pos_dict[second_nt][nt2num[last_nt]] += Pi_JinsDJ_given_D[D_in][nt2num[first_nt], init_pos] #base weight for middle pos nt
                                codon_prefix_dict[last_nt + second_nt][0] += Pi_JinsDJ_given_D[D_in][nt2num[first_nt], init_pos] #base weight for last pos nt
                for nt1 in 'ACGT':
                    if np.sum(second_pos_dict[nt1]) == 0:
                        second_pos_dict.pop(nt1, None)
                    for nt2 in 'ACGT':
                        if np.sum(codon_prefix_dict[nt1+nt2])== 0:
                            codon_prefix_dict.pop(nt1+nt2, None)
                            
    #            if len(second_pos_dict)> 0:
    #                print second_pos_dict
    #                return -1
                
                for delDr in range(num_delr_pos):
                    if self.min_delDl_given_DdelDr[D_in][delDr] == -1: # P(delDr | D) = 0 for this delDr --> move to next
                        continue
                    #Check if first nt from the D segment is okay
                    if cutD_gen_seg[l_D_seg - delDr - 1] in second_pos_dict.keys():
                        #The dell pos may be out of range of the PdelDldelDr_given_D -- check!
                        if l_D_seg - delDr - 1 <= self.max_delDl_given_DdelDr[D_in][delDr]:
                            Pi_R[:, init_pos - 1] += self.PdelDldelDr_given_D[l_D_seg - delDr - 1, delDr, D_in]*second_pos_dict[cutD_gen_seg[l_D_seg - delDr - 1]]
                    else:
                        continue #not okay, reject the alignment
                    
                    #Check if the second nt from the D segment is okay
                    if cutD_gen_seg[l_D_seg - delDr - 2:l_D_seg - delDr] in codon_prefix_dict.keys():
                        #The dell pos may be out of range of the PdelDldelDr_given_D -- check!
                        if l_D_seg - delDr - 2 <= self.max_delDl_given_DdelDr[D_in][delDr]:
                            Pi_R[0, init_pos - 2] += self.PdelDldelDr_given_D[l_D_seg - delDr - 2, delDr, D_in]*codon_prefix_dict[cutD_gen_seg[l_D_seg - delDr - 2:l_D_seg - delDr]][0]
                        base_prob = codon_prefix_dict[cutD_gen_seg[l_D_seg - delDr - 2:l_D_seg - delDr]][0]
                    else:
                        continue #no longer aligned, move to next delDr
                        
                    #Enter main loop
                    for pos in range(init_pos - 3, max(init_pos - l_D_seg + delDr, min_pos)-1, -1):
                        #note delDl = D_pos
                        D_pos = l_D_seg - delDr - 1 - ((init_pos - 1) - pos)
                        
                        #The dell pos may be out of range of the PdelDldelDr_given_D -- check!
                        if D_pos > self.max_delDl_given_DdelDr[D_in][delDr]:
                            current_PdelDldelDr = 0
                        else:
                            current_PdelDldelDr = self.PdelDldelDr_given_D[D_pos, delDr, D_in]
                        #Position is the first nt in codon
                        if pos%3 == 2:
                            #check alignment
                            if cutD_gen_seg[D_pos] in self.sub_codons_right[CDR3_seq[pos/3]]:
                                Pi_R[:, pos] += current_PdelDldelDr*base_prob*self.PD_nt_pos_vec[D_in][:, D_pos]
                            else:
                                break #no longer aligned -- exit loop
                        #Position is the second nt in codon
                        elif pos%3 == 1:
                            #check alignment
                            if cutD_gen_seg[D_pos:D_pos + 2] in self.sub_codons_right[CDR3_seq[pos/3]]:
                                Pi_R[:, pos] += current_PdelDldelDr*base_prob*self.PD_2nd_nt_pos_per_aa_vec[D_in][CDR3_seq[pos/3]][ :, D_pos]
                            else:
                                break #no longer aligned --- exit loop
                        #Position is the last nt in codon
                        else:
                            #check alignment
                            if cutD_gen_seg[D_pos:D_pos + 3] in self.codons_dict[CDR3_seq[pos/3]]:
                                Pi_R[0, pos] += current_PdelDldelDr*base_prob
                            else:
                                break #no longer aligned --- exit loop
                        
                    
                        
            #start position is second nt in a codon
            for init_pos in range(-2,-len(CDR3_seq)*3-1,-3):
                Pi_R[:, init_pos] += Pi_JinsDJ_given_D[D_in][:, init_pos]*self.zeroD_given_D[D_in]
                allowed_final_nts = ['ACGT'[nt] for nt in range(4) if Pi_JinsDJ_given_D[D_in][nt, init_pos] > 0]
                
                for delDr in range(num_delr_pos):
                    if self.min_delDl_given_DdelDr[D_in][delDr] == -1: # P(delDr | D) = 0 for this delDr --> move to next
                        continue
                    #check first nt of the D region (last in the codon)
                    if cutD_gen_seg[l_D_seg - delDr - 1] in allowed_final_nts: #first nt match
                        base_prob = Pi_JinsDJ_given_D[D_in][nt2num[cutD_gen_seg[l_D_seg - delDr - 1]], init_pos]
                        #The dell pos may be out of range of the PdelDldelDr_given_D -- check!
                        if l_D_seg - delDr - 1 <= self.max_delDl_given_DdelDr[D_in][delDr]:
                            Pi_R[0, init_pos-1] += self.PdelDldelDr_given_D[l_D_seg - delDr - 1, delDr, D_in]*base_prob
                    else:
                        continue #no alignment
                        
                    #Enter main loop
                    for pos in range(init_pos - 2, max(init_pos - l_D_seg + delDr, min_pos)-1, -1):
                        #note delDl = D_pos
                        D_pos = l_D_seg - delDr - 1 - ((init_pos - 1) - pos)
                        
                        #The dell pos may be out of range of the PdelDldelDr_given_D -- check!
                        if D_pos > self.max_delDl_given_DdelDr[D_in][delDr]:
                            current_PdelDldelDr = 0
                        else:
                            current_PdelDldelDr = self.PdelDldelDr_given_D[D_pos, delDr, D_in]
                        #Position is the first nt in codon
                        if pos%3 == 2:
                            #check alignment
                            if cutD_gen_seg[D_pos] in self.sub_codons_right[CDR3_seq[pos/3]]:
                                Pi_R[:, pos] += current_PdelDldelDr*base_prob*self.PD_nt_pos_vec[D_in][:, D_pos]
                            else:
                                break #no longer aligned -- exit loop
                        #Position is the second nt in codon
                        elif pos%3 == 1:
                            #check alignment
                            if cutD_gen_seg[D_pos:D_pos + 2] in self.sub_codons_right[CDR3_seq[pos/3]]:
                                Pi_R[:, pos] += current_PdelDldelDr*base_prob*self.PD_2nd_nt_pos_per_aa_vec[D_in][CDR3_seq[pos/3]][ :, D_pos]
                            else:
                                break #no longer aligned --- exit loop
                        #Position is the last nt in codon
                        else:
                            #check alignment
                            if cutD_gen_seg[D_pos:D_pos + 3] in self.codons_dict[CDR3_seq[pos/3]]:
                                Pi_R[0, pos] += current_PdelDldelDr*base_prob
                            else:
                                break #no longer aligned --- exit loop
                
            #start position is last nt in a codon
            for init_pos in range(-3,-len(CDR3_seq)*3-1,-3):
                Pi_R[0, init_pos] += Pi_JinsDJ_given_D[D_in][0, init_pos]*self.zeroD_given_D[D_in] 
                for delDr in range(num_delr_pos):
                    if self.min_delDl_given_DdelDr[D_in][delDr] == -1: # P(delDr | D) = 0 for this delDr --> move to next
                        continue
                    base_prob = Pi_JinsDJ_given_D[D_in][0, init_pos]
                    for pos in range(init_pos - 1, max(init_pos - l_D_seg + delDr, min_pos)-1, -1):
                        #note delDl = D_pos
                        D_pos = l_D_seg - delDr - 1 - ((init_pos - 1) - pos)
                        
                        #The dell pos may be out of range of the PdelDldelDr_given_D -- check!
                        if D_pos > self.max_delDl_given_DdelDr[D_in][delDr]:
                            current_PdelDldelDr = 0
                        else:
                            current_PdelDldelDr = self.PdelDldelDr_given_D[D_pos, delDr, D_in]
                        #Position is the first nt in codon
                        if pos%3 == 2:
                            #check alignment
                            if cutD_gen_seg[D_pos] in self.sub_codons_right[CDR3_seq[pos/3]]:
                                Pi_R[:, pos] += current_PdelDldelDr*base_prob*self.PD_nt_pos_vec[D_in][:, D_pos]
                            else:
                                break #no longer aligned -- exit loop
                        #Position is the second nt in codon
                        elif pos%3 == 1:
                            #check alignment
                            if cutD_gen_seg[D_pos:D_pos + 2] in self.sub_codons_right[CDR3_seq[pos/3]]:
                                Pi_R[:, pos] += current_PdelDldelDr*base_prob*self.PD_2nd_nt_pos_per_aa_vec[D_in][CDR3_seq[pos/3]][ :, D_pos]
                            else:
                                break #no longer aligned --- exit loop
                        #Position is the last nt in codon
                        else:
                            #check alignment
                            if cutD_gen_seg[D_pos:D_pos + 3] in self.codons_dict[CDR3_seq[pos/3]]:
                                Pi_R[0, pos] += current_PdelDldelDr*base_prob
                            else:
                                break #no longer aligned --- exit loop
                                
        return Pi_R