#
import os,pickle, sys 
import numpy as npy
from numpy import dot, transpose, sqrt, shape, matrix
from numpy.linalg import det, svd, eig
from math import  log, sin, cos, sqrt, pi, fabs, radians,degrees, acos, asin, atan2, sqrt, fmod, modf
from re import split
from string import lstrip,strip, rstrip
from libaction import Rotation
LNSEP = "#\n#" + "-" * 60
STARSEP = "#" + "*" * 60
from libmath  import getNormVect

        
class RDC(Rotation):
    """ functions for calculate RDC """
    # library for reading and outputting RDCs in xplor format 
    def normVect(self, vect):
        """ return a normalized vector)"""
        return vect/sqrt(vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2])    

    def bond2rdc(self, rAB):
        rAB = self.normVect(rAB)
        _theta  = acos(rAB[2])
        _phi = atan2(rAB[1],rAB[0])
        cos_theta = cos(_theta)
        cos2theta = cos_theta * cos_theta
        sin2theta = 1 - cos2theta
        rdc = self.Da * ((3*cos2theta -1) + 1.5*self.Rh *(sin2theta*cos(2*_phi)))
        return rdc 
    
    
    def thePhi2Vector(self, the, phi, vb=0):
        """ the, phi in units of degree """
        if vb: print "# thi,phi-->vec: the= %6.1f phi= %6.1f " % (the, phi)
        the = radians(the)
        phi = radians(phi)
        x = cos(phi)*sin(the)
        y = sin(phi) * sin(the)
        z = cos(the)
        return npy.array([x,y,z])
    
    def vec2ThePhe(self, rAB, **rest):
        """ return the, phi angles in units of Deg"""
        return self.bond2ThePhi(rAB, **rest)
    
    def bond2ThePhi(self, rAB, deg=0,  vb=0):
        """ return the, phi angles in units of Deg"""
        rAB = self.normVect(rAB)
        _theta  = acos(rAB[2])
        if fabs(rAB[1]) < 1.0e-6 and fabs(rAB[0]) < 1.0e-6:
            _phi = 0
        else:
            _phi = atan2(rAB[1],rAB[0])
        if vb:
            print "\n# bond(vector) to The, phi: "
            print "# vec:  %8.4f  %8.4f %8.4f " % (rAB[0], rAB[1], rAB[2])
            print "# bond-->the,phi: theta= %6.1f  phi= %6.1f " % (degrees(_theta), degrees(_phi))
            print "# -----------------------------\n"
        if deg:
            return degrees(_theta), degrees(_phi)
        else:
            return _theta, _phi
        
    
           
    def euler2rdc(self, rdcs, euler, vb=0):
        _euler = euler[0]
        lns = [ "# back calculating RDC from Euler angles, "]
        lns += [ "# euler: phi= %6.1f  theta= %6.1f  psi= %6.1f " % (_euler[0], _euler[1], _euler[2]) ]
        rot_m = self.getRotM(radians(_euler[0]), radians(_euler[1]), radians(_euler[2]))
        sum_sq = 0
        cnt = 0
        #plot_lns = ["# resid  exp_rdc  cal_rdc"]
        if len(rdcs[0]) == 5:
            RDC_ERR = 1
        else:
            RDC_ERR = 0
        rdc_dif_lns = []
        for i in range(len(self.rdc_vecs)):
            _bond_rdc = self.rdc_vecs[i]
            bond_vec, _rdc =  _bond_rdc
            if RDC_ERR:
                resid, atom_name1, atom_name2, xplor_rdc, _err = rdcs[i] 
            else:
                resid, atom_name1, atom_name2, xplor_rdc  = rdcs[i] 
            assert xplor_rdc == _rdc
            rot_vec =  self.rotVec(rot_m, bond_vec)
            back_rdc = self.bond2rdc(rot_vec)
            rdc_dif = _rdc - back_rdc
            sum_sq += rdc_dif * rdc_dif
            cnt += 1
            rdc_dif_lns += [ "# resid= %3d  %s - %s rdc= %8.3f   back_rdc= %8.3f  rdc_dif= %8.3f" % (resid, atom_name1, atom_name2, _rdc,  back_rdc, rdc_dif) ]
            #plot_lns += [ "%4d  %8.4f  %8.4f " % (resid, _rdc, back_rdc)]
        rmsd = sqrt(sum_sq/float(cnt))
        ln = "# svd: num_rdc= %3d  rdc rmsd= %8.3f " % (cnt, rmsd) 
        if vb: 
            lns += rdc_dif_lns
            lns += [ln]
            self.add2log(lns)
        else:
            lns += [ln]
            self.add2log(lns + rdc_dif_lns)
        self.prnLns(lns)
        
        
    def getVecAngle(self, vec1, vec2):
        """ return angle between vec1 and vec2 in radians """
        vec1 = self.normVect(vec1)
        vec2 = self.normVect(vec2)
        return acos(vec1[0] * vec2[0] + vec1[1]*vec2[1] + vec1[2]*vec2[2])
    
class XplorRDC(RDC):
    """ Read and Write RDC in Xplor format """
    num_media = 0
    def getBondXplorRDC(self,  bond_rdc,  atom_name1, resid1, atom_name2, resid2=None, segid=None):
        """ generate one bond RDC in Xplor format """
        if resid2 is None: resid2 = resid1
        rdc_lns = ["", "assign ( resid 500  and name OO   )" ,
                        "       ( resid 500  and name Z    )",
                        "       ( resid 500  and name X    )",
                        "       ( resid 500  and name Y    )"]
        if segid is not None:
            rdc_lns += ["       ( segid %s and resid %3d  and name %-2s   )" % (segid, resid1, atom_name1),
                        "       ( segid %s and resid %3d  and name %-2s   ) %8.4f   1.2000 1.2000 " % (segid, resid2, atom_name2, bond_rdc)]
        else:
            rdc_lns += ["       ( resid %3d  and name %-2s   )" % (resid1, atom_name1),
                        "       ( resid %3d  and name %-2s   ) %8.4f   1.2000 1.2000 " % (resid2, atom_name2, bond_rdc)]
        return rdc_lns 


    def readXplorRDC(self, fn_rdc_xp=None,rdc_cid=7, rdc_err_cid=8, vb=0):
        """ rdc in Xplor format:
        assign ( resid 500  and name OO  )
       ( resid 500  and name Z   )
       ( resid 500  and name X   )
       ( resid 500  and name Y   )
       ( resid 44   and name N   )
       ( resid 44   and name HN  )   -7.30  0.5000
       return rdc in dictionary format 
       """
        if fn_rdc_xp is None: fn_rdc_xp = self.fn_rdc_xp 
        if vb: print "# fn_rdc_xp: ", fn_rdc_xp 
        self.num_media += 1
        txt = self.readTxtFile(fn_rdc_xp)
        num_lns = len(txt)
        i = 0
        bond_lengths = []
        rdcs_dt = {}
        rdc_lns = []
        while i < num_lns:
            ln = txt[i]
            ln = ln.strip()
            if ln and ln[0:6] == "assign":
                i += 4 
                tmp = split("\s+", txt[i].strip())
                resid1 = int(tmp[2])
                aname = tmp[5]
                i+=1 
                tmp = split("\s+", txt[i].strip())
                resid2 = int(tmp[2])
                aname2 = tmp[5]
                #print "# rdc_cid= %4d " % rdc_cid, tmp[rdc_cid]
                rdc = float(tmp[rdc_cid])
                rdc_err = float(tmp[rdc_err_cid])
                assert resid1 == resid2 
                if vb: print "%4d  %4s - %4s   %8.2f  %8.2f" % (resid1, aname, aname2, rdc, rdc_err)
                rdcs_dt[resid1] = (resid1, aname, aname2, rdc, rdc_err)
            i+=1
        return rdcs_dt
    
    
    
    def rdc_dt2li(self, rdc_dt):
        """ rdcs in list format from dict """
        resids = rdc_dt.keys()
        resids.sort()
        rdcs = []
        for resid in resids:
            rdcs.append(rdc_dt[resid])
        return rdcs
    
    
    def backRDCs(self, rdcs):
        rdc_lns = ["!RDC_restraints"]
        for _rdc_ob in rdcs:
            resid = _rdc_ob[0]
            atom1 = _rdc_ob[1]
            atom2 = _rdc_ob[2]
            rBond = self.residues[resid].getBondVect(atom1, atom2)
            rdc = self.bond2rdc(rBond)
            base_name = self.residues[resid].name
            print "base: %s%-3d  bond: %s-%s rdc= %8.4f"  % (base_name, resid, atom1, atom2, rdc)
            rdc_lns += self.getBondXplorRDC(rdc, atom1, resid, atom2)
        #prnLns(rdc_lns)
        #self.saveLines(rdc_lns, "sim_rdc.tbl",vb=1, cmt="back calculated RDC")
        return rdc_lns 
    
    def rdcs2XplorTable(self, rdcs):
        """ save rdcs in Xplor format: rcds(resid1, aname, aname2, rdc)"""
        rdc_lns = ["!RDC_restraints"]
        for _rdc_ob in rdcs:
            resid, atom1, atom2, rdc  = _rdc_ob
            rdc_lns += self.getBondXplorRDC(rdc, atom1, resid, atom2)
        return rdc_lns 
    
    
    def rdcs2ColData(self, rdcs_dt, vb=0):
        """ output rdc in colum format: (resid, rdc) for Powell fitting """
        resids = rdcs_dt.keys()
        resids.sort()
        lns = []
        for resid in resids:
            lns.append("%4d %8.2f " % (resid, rdcs_dt[resid][-2]))
        
        if vb:
            print "# num_rdc= %3d " % len(lns)
            self.prnLns(lns)
        return lns 
    
    def getRDCID4Fit(self, rdcs_dt, resids=None, helix_init_id=1, vb=0):
        """ return rdc with index list for RDC fitting:
        rdcs: [(2, 'N', 'HN', 2.1899999999999999, 0.5), 
               (3, 'N', 'HN', 0.90000000000000002, 0.5),
               ...]
            """
        rdc_ids = []
        rdc_vals = []
        rdc_resids = []
        resids_dt = rdcs_dt.keys()
        if resids is None: 
            resids = resids_dt
        resids.sort()
        resids_order = range(resids[0], resids[-1]+1)
        for resid in resids_order:
            if resid in resids and resid in resids_dt:
                _rdc = rdcs_dt[resid]
                rdc_ids.append(helix_init_id)
                rdc_vals.append(_rdc[3])
                rdc_resids.append(_rdc[0])
                if vb: print "#hid= %3d  resid= %3d  rdc= %8.2f " % ( helix_init_id,  _rdc[0], _rdc[3])
            helix_init_id += 1
        #self.bpids[i], self.dp_rdc_resids[i]
        self.bpids = rdc_ids
        self.dp_rdc_resids = rdc_resids
        return rdc_ids, rdc_vals

    
    def prnXpRDCs(self, rdcs_dt):
        """ print rdcs from xplor file: (40, 'N3', 'H3', 8.37) """
        print "# rdc from xplor table "
        print "# No.     resid  Atom--Atom   RDC  RDC_error"
        resids = rdcs_dt.keys()
        resids.sort()
        num =0 
        for resid in resids:
            num += 1
            print "# %3d " % num, 
            print "  %3d  %s -- %s  %7.2f %7.2f " % rdcs_dt[resid]
    
    
    def getRDCValues(self, rdcs, vb=0):
        """ return rdc values from rdc sets """
        rdc_vals = []
        for _rdc in rdcs:
            rdc_vals.append(_rdc[3])
        if vb: print "# rdc_vals: ", rdc_vals
        return rdc_vals
    
    
    def rmRDCs(self, rdcs, resids):
        """ remove RDCs not included in resids: rcds(resid1, aname, aname2, rdc) """
        left_rdcs = []
        rdc_resids = []
        for _rdc in rdcs:
            if _rdc[0] in resids:
                left_rdcs.append(_rdc)
                rdc_resids.append(_rdc[0])
        return left_rdcs, rdc_resids
  
    
    
    def readXpRDCs4Tensor(self, fn_xp_rdc):
        """ read rdc in xplor format and remove rdcs not in residues in memory """
        rdc_dt = self.readXplorRDC(fn_xp_rdc)
        rdcs = self.rdc_dt2li(rdc_dt)
        self.rdcs, self.rdc_resids = self.rmRDCs(rdcs, self.resids)
        self.rdc_dt = rdc_dt 
        
    
    def rdc2Euler_svd(self, vb=0):
        """ read RDC from file fn_xp_rdc in XPlor format and return Euler angles """
        euler = self.rdc2euler(self.rdcs)
        self.Rh = self.R 
        lns = self.prnEuler_Xplor_Form(euler, DaR=(self.Da, self.R))
        self.add2log(lns)
        self.euler2rdc(self.rdcs, euler, vb=vb) 
        self.svd_tensor = lns 
    
    def NH2Cos(self, resids=None, vb=0, li=0):
        """ return NH vector's projectors in xy, yz, xz planes:
        cosine(cos(alpha), cos(beta), cos(gamma))         
        """
        lns = ["# NH vector unit projections in xy,yz,xz planes"]
        if li:
            xyz_cos = []
        else:
            xyz_cos = {}
        if resids is None: resids = self.resids
        
        for resid in resids:
            ressid = self.resid2res_sid(resid)
            NH_vec = self.residues[ressid].getImino() #getNH()
            #print "# resid= %3d ressid= %3d  NH: %7.2f  %7.2f %7.2f " % (resid, ressid, NH_vec[0], NH_vec[1], NH_vec[2])
            if NH_vec is not None:
                #print NH_vec
                norm_NH = getNormVect(NH_vec) #(NH_vec)
                #print "# norm_NH: ", norm_NH
                x = norm_NH[0] #NH_vec[0]/norm_NH
                y = norm_NH[1]#NH_vec[1]/norm_NH
                z = norm_NH[2] #NH_vec[2]/norm_NH
                lns.append( "%3d\t%3.3f\t%3.3f\t%3.3f" % ( resid, x, y, z))
                if li:
                    xyz_cos.append( (x,y,z) )
                else:
                    xyz_cos[resid] = (x,y,z)
        lns = ["# **********    End of projections *********"]
        if vb: 
            self.prnLns(lns)        
            self.add2log(lns)
        return npy.array(xyz_cos)
    
    
    def chkNorm(self, xyz_cos):
        xyz_resids = xyz_cos.keys()
        xyz_resids.sort()
        for resid in  xyz_resids:
            x,y,z = xyz_cos[resid]
            nomr = x*x + y*y + z*z 
            print "# resid= %4d nomr=  %8.4f "  % (resid, nomr)
            
    
    
    def prnPowFitRes(self, fit_res, cmt= ""):
        ln =  "# Powell %s : phi= %6.1f the= %6.1f psi= %6.1f  Da= %6.2f  R= %6.2f Dc= %6.2f " %   \
              (cmt, degrees(fit_res[0]), degrees(fit_res[1]), degrees(fit_res[2]), fit_res[3], fit_res[4], fit_res[5])
        print ln 
        return ln 
    
    def getPowellFitInit(self):
        ptp =  self.euler_angles[0]
        init_pow = [radians(ptp[0]), radians(ptp[1]),  radians(ptp[2]),  self.Da,  self.R, 0] 
        return init_pow
    
    
        
    def rdc2Euler_powell(self, init_pm):
        
        rdc_vals = self.getRDCValues(self.rdcs, vb=0)
        
        run_pow = npy.array(init_pm)
        xyz_cos_arr = self.NH2Cos(resids=self.rdc_resids, li=1)
        c_pow_align(6, xyz_cos_arr, npy.array(rdc_vals), run_pow)
        lns = ["# *********Begin of powell fiting "]
        lns += [ self.prnPowFitRes(init_pm, cmt=" in")]
        lns += [self.prnPowFitRes(run_pow, cmt="out") ]
        lns += ["#  ******** End of powell fiting "]
        self.add2log(lns)
        self.xyz_cos_arr = xyz_cos_arr 
        self.rdc_vals  = rdc_vals
        self.powell_pm = run_pow
        self.pow_euler = self.gen4ptp(degrees(run_pow[0]), degrees(run_pow[1]), degrees(run_pow[2]))
        self.pow_Da = run_pow[3]
        self.pow_R = run_pow[4]
        self.pow_DC = run_pow[5]
        self.pow_tensor = self.prnEuler_Xplor_Form(self.pow_euler, DaR=(self.pow_Da, self.pow_R), cmt=self.prnPowFitRes(run_pow))
        self.add2log(self.pow_tensor)
        return run_pow
        
        
    def euler2rdc_powell(self, pm, vb=0):
        back_rdc_vals = c_solve(6, self.xyz_cos_arr, npy.array(pm))
        lns = [ "\n# back calculating RDC from Euler angles from Powell fitting "]
        lns += [ self.prnPowFitRes(pm, cmt="result")] 
        sum_sq = 0
        num_rdc = len(back_rdc_vals)
        rdc_dif_lns = []
        for i in range(num_rdc):
            b_rdc = back_rdc_vals[i]
            e_rdc = self.rdc_vals[i]
            _rdc = self.rdcs[i] # 
            rdc_dif = e_rdc - b_rdc
            sum_sq += rdc_dif * rdc_dif
            rdc_dif_lns += [ "# resid= %3d  %s - %s rdc= %8.3f   back_rdc= %8.3f  rdc_dif= %8.3f" % (_rdc[0], _rdc[1], _rdc[2], e_rdc,  b_rdc, rdc_dif) ]
            #plot_lns += [ "%4d  %8.4f  %8.4f " % (resid, _rdc, back_rdc)]
        rmsd = sqrt(sum_sq/float(num_rdc))
        ln = "# svd: num_rdc= %3d  rdc rmsd= %8.3f " % (num_rdc, rmsd) 
        if vb: 
            lns += rdc_dif_lns
            lns += [ln]
            self.add2log(lns)
        else:
            lns += [ln]
            self.add2log(lns + rdc_dif_lns)
        self.prnLns(lns)
        
        
        
         
    
class DuplexRDC(XplorRDC):
    """ read and parse rdcs in simple table text file """
        
    def rmRDC_byID(self, rm_id):
        rm_val = self.bp_ids[rm_id]
        self.bp_ids.remove(rm_val)
        rm_val = self.bp_rdcs[rm_id]
        self.bp_rdcs.remove(rm_val)
        
        
    def getDelID(self, rdc_difs):
        """ return the index of the rdc with largest absolute error """
        rdc_difs = [ fabs(_rdc) for _rdc in rdc_difs]
        rdc_diff_max = max(rdc_difs)
        max_dif_id =  rdc_difs.index(rdc_diff_max)
        #print "# max_dif_id= %3d " % max_dif_id
        #print "# rdc_difs: ", rdc_difs 
        return max_dif_id 
    
    
    def parseExpRDC(self, vb=0):
        """ parse input experimental RDCs """
        lns = self.readTxtFile(self.fn_exp_rdc)
        rdcs = self.txt2col(lns)
        exp_rdcs = {}
        unknow_rdcs = []
        for _rdc in rdcs:
            if "X" not in _rdc[1]:
                nt_name = _rdc[1][0]
                nt_resid = int(_rdc[1][1:])
                rdc_val = float(_rdc[3])
                rdc_quality = _rdc[2]
                print "# %3d %s  %s  %8.4f "  % (nt_resid, nt_name,  rdc_quality, rdc_val)
                exp_rdcs[nt_resid] = (nt_name, nt_resid, rdc_val, rdc_quality)
            else:
                unknow_rdcs.append(_rdc)
        if vb: 
            print "# experimental RDCs: "
            self.prnDict(exp_rdcs)
            print "# unknow rdcs: ", unknow_rdcs
        #self.exp_rdc = rdcs
        return exp_rdcs
    
    
    def parseXpRDC(self, vb=0):
        """ parse input experimental RDCs in Xplor format  """
        rdcs_dt = self.readXplorRDC(self.fn_exp_rdc)
        
        #print self.seq_dict 
        # rdcs_dt[resid1] = (resid1, aname, aname2, rdc, rdc_err)
        resids = rdcs_dt.keys()
        exp_rdcs = {}
        ex_rdcs = getattr(self, 'ex_rdcs', [])
        for resid in resids:
            if resid in ex_rdcs: continue 
            _rdc_xp = rdcs_dt[resid]
            if _rdc_xp[1] == 'N1':
                nt_name = 'G'
            else:
                nt_name = 'U'
            rdc_val = _rdc_xp[3]
            rdc_quality = _rdc_xp[4]
            if vb: print "# %3d %s  %s  %8.4f "  % (resid, nt_name,  rdc_quality, rdc_val)
            exp_rdcs[resid] = (nt_name, resid, rdc_val, rdc_quality)
        if vb: 
            print "# experimental RDCs: "
            self.prnDict(exp_rdcs)
        self.exp_rdc_resids = exp_rdcs.keys()
        self.exp_rdcs = exp_rdcs
        self.xp_rdcs = exp_rdcs
    
    
    def parseExpRDC3Col(self, vb=0): #passed
        """ parse input experimental RDCs: rdc in 3 colum format: 
        # index name_resid  RDC   
        #duplex #1
        1	     G14         11.115
        """
        lns = self.readTxtFile(self.fn_exp_rdc)
        rdcs = self.txt2col(lns)
        exp_rdcs = {}
        lns = ["# Parsing experimental RDCs: fn_rdc= %s "  % self.fn_exp_rdc]
        lns += [ "# Index  resid  name   rdc "]
        cnt = 1
        for _rdc in rdcs:
            nt_name = _rdc[1][0]
            nt_resid = int(_rdc[1][1:])
            rdc_val = float(_rdc[2])
            lns += [ "#%4d   %5d    %s   %8.4f "  % (cnt, nt_resid, nt_name,   rdc_val) ]
            cnt += 1
            exp_rdcs[nt_resid] = (nt_name, nt_resid, rdc_val)
        self.add2log(lns)
        self.exp_rdcs = exp_rdcs
        self.exp_rdc_resids = self.exp_rdcs.keys()
        #return exp_rdcs
    
    
    def loadExpRDC(self, vb=0):
        """ read input RDC and transformed into block txt """
        lns = self.readTxtFile(self.fn_exp_rdc)
        rdcs = self.txt2col(lns)
        if vb: 
            print "# read experimental RDCs into block text "
            print rdcs
        return rdcs 
    
        
    def setDup_Bps(self, vb=0):
        """ set duplex basepair resids """
        self.dup_bps = {}
        for hkey in self.duplexes.keys():
            self.dup_bps[hkey] = self.resid2BPresid(self.duplexes[hkey], vb=vb)
            
    
    def old_getBpId(self, hkey, resids, vb=0):
        """ return basepair ids for rdc fitting
            resids: residue Ids with experimental RDCs 
        """
        dbps = self.dup_bps[hkey]
        bp_ids = []
        bp_rdcs = []
        dp_rdc_resids = []
        bp_id = None 
        for i in range(len(dbps)):
            _bp = dbps[i] 
            if _bp[0] in resids:
                if bp_id is None: bp_id = 1
                bp_ids.append(bp_id)
                bp_rdcs.append(self.exp_rdcs[_bp[0]][2])
                dp_rdc_resids.append(_bp[0])
                if vb: print "# bp_id= %3d  rdc= %8.4f  resid= %3d " % (bp_ids[-1], bp_rdcs[-1], _bp[0])
            elif  _bp[1] in resids:
                if bp_id is None: bp_id = 1
                bp_ids.append(bp_id)
                bp_rdcs.append(self.exp_rdcs[_bp[1]][2])
                dp_rdc_resids.append(_bp[1])
                if vb: print "# bp_id= %3d  rdc= %8.4f  resid= %3d " % (bp_ids[-1], bp_rdcs[-1], _bp[1])
            if bp_id is not None: bp_id += 1 
        if vb: 
            print "# bp_ids: ", bp_ids
            print "# bp_rdcs: ", bp_rdcs
        return bp_ids, bp_rdcs, dp_rdc_resids
    
    def setDupPhaseResid(self, resid, bp_id):
        if bp_id is None and self.dup_phase_resid is None: 
            bp_id = 1
            self.dup_phase_resid = resid
        return bp_id     
    
        
    def setBpIds(self, hkey, vb=0): # passed
        """ return basepair ids for rdc fitting
            resids: residue Ids with experimental RDCs 
        """
        dbps = self.duplexs[hkey]
        self.bp_ids = []
        self.bp_rdcs = []
        self.dp_rdc_resids = []
        self.dup_phase_resid = None 
        lns = ["# duplex %s basepair index and rdcs " % hkey]
        bp_id = None 
        for i in range(len(dbps)):
            _bp = dbps[i] 
            if _bp[0] in self.exp_rdc_resids:
                if _bp[1] in self.exp_rdc_resids:
                    if self.xp_rdcs[_bp[1]][0] == 'G':
                        bp_id = self.setDupPhaseResid(_bp[1], bp_id)
                        self.bp_ids.append(bp_id)
                        self.bp_rdcs.append(self.exp_rdcs[_bp[1]][2])
                        self.dp_rdc_resids.append(_bp[1])
                        lns += [ "# bp_id= %3d  rdc= %8.4f  resid= %3d " % (self.bp_ids[-1], self.bp_rdcs[-1], _bp[1]) ]
                        bp_id += 1
                        continue 
                bp_id = self.setDupPhaseResid(_bp[0], bp_id)
                self.bp_ids.append(bp_id)
                self.bp_rdcs.append(self.exp_rdcs[_bp[0]][2])
                self.dp_rdc_resids.append(_bp[0])
                lns += [ "# bp_id= %3d  rdc= %8.4f  resid= %3d " % (self.bp_ids[-1], self.bp_rdcs[-1], _bp[0]) ]
            elif  _bp[1] in self.exp_rdc_resids:
                bp_id = self.setDupPhaseResid(_bp[1], bp_id)
                self.bp_ids.append(bp_id)
                self.bp_rdcs.append(self.exp_rdcs[_bp[1]][2])
                self.dp_rdc_resids.append(_bp[1])
                lns += [ "# bp_id= %3d  rdc= %8.4f  resid= %3d " % (self.bp_ids[-1], self.bp_rdcs[-1], _bp[1]) ]
            if bp_id: bp_id += 1
            
        lns += [ "# self.bp_ids: %s "  %  self.bp_ids ]
        lns += [ "# self.bp_rdcs: %s" % self.bp_rdcs ]
        self.add2log(lns)
        



if (__name__ == "__main__"):
    pass 
    
    