# assemble RNA duplexes by relative orientations from RDC 
# 
# with new rotation rules, common molecular frame, terminal loops 

import os
from glob import glob

from string import lstrip,strip, rstrip
from re import split

import numpy 
from numpy import dot, transpose, sqrt, shape, array
from numpy.linalg import det, svd  

from copy import copy 
from math import sin, cos, fabs, sqrt, radians, degrees, acos, atan2

from librna import nt_atoms_woh, base_anames, nt_atoms_woh, nt_bb_atoms, RNAPDB, Duplex
from libio import INPUT
from libaction import HomoMove, Rotation
from librdcwave import RDCDEG
from libblock import BlockOther


class LINKORDER:
    """ Generate link order to connect all structure motifs """
    
    def RNA2D2motif(self, vb=0):
        bkeys = self.sec_inf.keys()
        #print self.sec_inf 
        if 'dup_resids' in bkeys: bkeys.remove('dup_resids')
        self.rna_motifs = []
        for _bkey in bkeys:
            bks = self.sec_inf[_bkey]
            #((31, 39, 'UU'), [32, 33, 34, 35, 36, 37, 38])       
            if type(bks) == type([]):
                for i in range(len(bks)):
                    bk = bks[i]
                    bk_resids = bk[-1]
                    #print "# bks: ", bks 
                    self.rna_motifs.append(((bk_resids[0], bk_resids[-1]),(_bkey, i, len(bk_resids))))
            elif type(bks) == type({}):        
                dup_ids = bks.keys()
                for dup_id in dup_ids:
                    dup = bks[dup_id]
                    first_bp = dup[0]
                    last_bp = dup[-1]
                    self.rna_motifs.append(((first_bp[0], first_bp[1], last_bp[0], last_bp[1]), ('helix', dup_id)))
            else:
                print "# unformal blocks"
        if vb: 
            print "# self.rna_motifs: \n", self.rna_motifs          
        self.tmp_motifs = copy(self.rna_motifs)
        self.add2log([ "# num_motifs= %3d " % len(self.tmp_motifs) ] )
        

    def addBBlockID(self, bb):
        """ add building block by the direction: 5'--> 3'"""
        if self.link_order is None:
            if  bb[0] == 'helix': #
                self.link_order = [bb] # ('helix', dup_id)]
                key = "h%s" % bb[1]
                self.motifs = {key:bb}
            else:
                raise "# Not a dulex, only duplex could be set as the anchor building block:",  bb 
        elif bb in self.link_order:
            pass 
        else:
            self.link_order.append(bb)
            key = "%s%s" % (bb[0][0], bb[1])
            self.motifs[key] = bb 
                
                                            
    def getDuplexNextRes3(self, dup_id, resid=None):
        """ return the 3' resid of the duplex  3'-5'  terminal basepair"""
        dup_bps = self.duplexs[dup_id]
        first_bp = dup_bps[0]
        last_bp = dup_bps[-1]
        if resid is None or resid==first_bp[0]:
            res5, res3, bp_seq = last_bp
            next_res = res5 + 1 
        elif resid == last_bp[1]:
            next_res = first_bp[1] + 1
        else:
            self.err_exit("Duplex next res error ")
        return next_res
    
    def getDup3PreResid(self, dup_id, resid=None):
        """ return the previous resid of 3' side chain of the duplex """
        dup_bps = self.duplexs[dup_id]
        first_bp = dup_bps[0]
        last_bp = dup_bps[-1]
        if resid is None or resid == first_bp[1]:
            pre_res = last_bp[1] - 1
        elif resid == last_bp[0]:
            pre_res = first_bp[0] - 1
        else:
            self.err_exit("Duplex previous resid error ")
        return pre_res
    
    
    def getDuplexPreRes(self, dup_id):
        """ return the 3' resid of the duplex  3'-5'  terminal basepair"""
        dup_bps = self.duplexs[dup_id]
        res5, res3, bp_seq  = dup_bps[0]
        pre_res = res5 - 1
        return pre_res
    
        
            
    def setAnchorDuplexID(self):
        """ set longest duplex as the achor duplex for packing process """
        dup_lengths = []
        for _dkey in self.duplex_ids:
            dup_lengths.append(len(self.duplexs[_dkey]))
        tmp_id = dup_lengths.index(max(dup_lengths))
        dup_id = self.duplex_ids[tmp_id]
        self.addBBlockID(('helix',dup_id))
        self.anchor_dup_id = dup_id 
        self.anchor_end = 0
        return self.getDuplexNextRes3(dup_id)
    
    
    def res2bb(self, resid):
        """ resid to building block"""
        next_res = None 
        #print self.tmp_motifs 
        for _ttif  in self.tmp_motifs: 
            _mresid, _motif = _ttif
            if resid in _mresid:
                if len(_mresid) == 2:
                    #print "# mresids: ", _mresid 
                    #print "# resid: ", resid 
                    if _motif[0] == 'linker' and _motif[-1] > 4: 
                        return 'end' 
                    self.addBBlockID(_motif)
                    assert _mresid[0] < _mresid[1]
                    if resid == _mresid[0]:
                        next_res = _mresid[1]+1
                    else:
                        next_res = _mresid[0] -1 
                elif len(_mresid) == 4:
                    self.addBBlockID(_motif)
                    if resid == _mresid[0]:
                        next_res = _mresid[2] + 1
                    elif resid == _mresid[3]:
                        next_res = _mresid[1] + 1
                    elif resid == _mresid[2]:
                        next_res = _mresid[0] - 1
                    else:
                        next_res = _mresid[3] - 1
                else:
                    print "# non block"
                break 
        #print "# resid: ", resid, "next_res: ", next_res ,  "motif: ", _ttif
        return next_res
    
    def res2Prebb(self, resid):
        """ resid to building block"""
        #print "# pre: resid: ", resid 
        next_res = None 
        #print self.tmp_motifs 
        for _ttif  in self.tmp_motifs: 
            _mresid, _motif = _ttif
            if resid in _mresid:
                if len(_mresid) == 2:
                    #print "# mresids: ", _mresid 
                    #print "# resid: ", resid 
                    #print "# motif: ", _motif 
                    self.addBBlockID(_motif)
                    assert _mresid[0] < _mresid[1]
                    if resid == _mresid[0]:
                        next_res = _mresid[1]+1
                    else:
                        next_res = _mresid[0] -1 
                elif len(_mresid) == 4:
                    self.addBBlockID(_motif)
                    if resid == _mresid[0]:
                        next_res = _mresid[2] + 1
                    elif resid == _mresid[3]:
                        next_res = _mresid[1] + 1
                    elif resid == _mresid[2]:
                        next_res = _mresid[0] - 1
                    else:
                        next_res = _mresid[3] - 1
                else:
                    print "# non block"
                break 
        #print "# resid: ", resid, "next_res: ", next_res ,  "motif: ", _ttif, "next_res: ", next_res
        return next_res
    

    def setLinkOrder(self, vb=0):
        """ choose longest duplex with orientation from RDC if availabe as the achor duplex.
        Then add nucleotides from 5' to 3'.
        """
        next_res = self.setAnchorDuplexID()
        cnt = 0 
        while next_res and next_res  != "end":
            next_res = self.res2bb(next_res)
            if next_res == 'end': break
            cnt += 1
            if cnt > 12: break 
        cnt = 0 
        next_res = self.getDup3PreResid(self.anchor_dup_id)
        while True:
            next_res = self.res2Prebb(next_res)
            #print "### pre resid: ", next_res
            if next_res == 'end': break
            cnt += 1
            if cnt > 5: break 
        self.add2log( ["# link_oder: %s " % self.link_order], vb=vb )    
    

class DUPLEX3D(Duplex, RDCDEG):
    def __init__(self):
        Duplex.__init__(self)
        
    def rotBackPub(self, tpr, hex_xyz):
        the, phi, rho = tpr 
        z_axis = array([0., 0., 1.])
        y_axis = array([0.,1., 0.])
        print "# back to command moleuclar frame \n phi= %6.1f the= %6.1f rho= %6.1f " % (phi, the, rho)
        new_pdb_xyz = {}
        resids = hex_xyz.keys()
        for resid in resids:
            xyz = hex_xyz[resid]
            new_xyz = []
            for _xyz in xyz:
                _xyz = _xyz.copy()
                _xyz = self.new_arb_rot(_xyz, z_axis, -phi)
                _xyz = self.new_arb_rot(_xyz, y_axis, -the) #, hex_axis)
                _xyz = self.new_arb_rot(_xyz, z_axis, -rho)#, hex_axis)
                new_xyz.append(_xyz)
            new_pdb_xyz[resid] = array(new_xyz)
        #self.new_pdb_xyz = new_pdb_xyz 
        #pdb_lns = self.resDict2pdb()
        #return pdb_lns  
        return new_pdb_xyz
    
    def thephi2rho(self, the, phi,  phase_resid, rev, phase_nt_name, duplex_xyz):
        """ calculate phase angle by rotate duplex to public frame """
        #print "# duplex  the= %6.1f phi= %6.1f" % (phi, the)
        #self.prnDict(self.seq_dict)
        #nt_name =  self.seq_dict[self.phase_iresid]
        xyz = duplex_xyz[phase_resid]
        rNH = self.resXYZ2Imino(numpy.array(xyz), phase_nt_name)
        if rev: rNH = -rNH
        rNH = self.new_arb_rot(rNH, self.z_axis, -phi)
        rNH = self.new_arb_rot(rNH, self.y_axis, -the) #, hex_axis)
        the, phi = self.vector2thePhi(rNH, vb=0)
        rho = phi # first imino proton' projector on xy plane. 
        #print "# phase rho= %6.1f " % rho 
        return rho 
    
    def rotDup2Pub(self, the, phi, rho, dup_xyz, trans=None):
        """ rotate duplex to public frame by duplex orientation(the, phi, rho)"""
        lns = [  "# rotDup2Pub(): duplex DRO --> Pub :  the= %6.1f phi= %6.1f  rho= %6.1f " % ( the, phi, rho) ]
        self.add2log(lns)
        resids = dup_xyz.keys()
        new_pdb_xyz = {}
        for resid in resids:
            xyz = dup_xyz[resid]
            new_xyz = []
            for _xyz in xyz:
                _xyz = _xyz.copy()
                _xyz = self.new_arb_rot(_xyz, self.z_axis, -phi)
                _xyz = self.new_arb_rot(_xyz, self.y_axis, -the) 
                _xyz = self.new_arb_rot(_xyz, self.z_axis, -rho)
                new_xyz.append(_xyz)
            new_pdb_xyz[resid] = numpy.array(new_xyz)
        return new_pdb_xyz
    
    def rotPubDup2TPR(self, the, phi, rho, dup_xyz, trans=None):
        """ rotate duplex from public from to its dros(the, phi, rho)"""
        lns = [  "# rotPubDup2TPR():  pub -> DROs:  the= %6.1f phi= %6.1f  rho= %6.1f " % ( the, phi, rho) ]
        resids = dup_xyz.keys()
        new_pdb_xyz = {}
        for resid in resids:
            xyz = dup_xyz[resid]
            new_xyz = []
            for _xyz in xyz:
                _xyz = _xyz.copy()
                _xyz = self.new_arb_rot(_xyz, self.z_axis, rho)
                _xyz = self.new_arb_rot(_xyz, self.y_axis, the) 
                _xyz = self.new_arb_rot(_xyz, self.z_axis, phi)
                new_xyz.append(_xyz)
            new_pdb_xyz[resid] = numpy.array(new_xyz)
        self.add2log(lns)
        return new_pdb_xyz
    
    
    def rotPub2refTPR(self, ref_tpr, hex_xyz, vb=0):
        """ rotate duplex from  helix axis tpr(the, phi, rho) to ref_tpr(calculated by experimental RDCs)"""
        the, phi, rho = ref_tpr 
        if vb: 
            print "# rot duplex with public orientation to target orientation"
            print "# the= %6.1f  phi= %6.1f rho= %6.1f " % ref_tpr
        return self.rotPubDup2TPR(the, phi, rho, hex_xyz)
        
    
    def resid2bp(self, resids):
        """ resids to basepair """
        resids.sort()
        bps = []
        num_bp = len(resids)/2 
        for i in range(num_bp):
            bps.append((resids[i],resids[-(i+1)]))
        return bps 
                        
    
    def rotXYZ_by_refTPR(self, rot, xyz):
        hex_axis, rot_axis, z_axis, d_the, d_rho, d_phi = rot 
        new_xyz = []
        for _xyz in xyz:
            #print _xyz
            _xyz = _xyz.copy()
            _xyz = self.new_arb_rot(_xyz, hex_axis, d_rho) #, hex_axis)
            _xyz = self.new_arb_rot(_xyz, rot_axis, d_the)
            _xyz = self.new_arb_rot(_xyz, z_axis, d_phi) #, hex_axis)
            new_xyz.append(_xyz)
        #print len(new_xyz)     
        return array(new_xyz)
            
    
    
    def chkDupXYZDRO(self, dup_id, dup_xyz, vb=0):
        """ check 3D duplex' DRO """
        resids = self.helices[dup_id]
        bps = self.resid2BPresid(resids)
        phase_resid = self.dup_phresids[dup_id]
        rev = self.isReverseImino(phase_resid, bps)
        
        dup_axis = self.pp2axis(bps, dup_xyz, vb=vb)
        the, phi = self.axis2ThePhi(dup_axis)
        rho = self.thephi2rho(the, phi, phase_resid, rev, self.seq_dict[phase_resid], dup_xyz)
        if vb: 
            print "#bps: ", bps 
            print "# phase: %3s rev: %s " % (phase_resid, rev)
        ln =  "# dup %s: the= %6.1f phi= %6.1f rho= %6.1f " % (dup_id, the, phi, rho)
        self.add2log([ln], vb=vb)
        
    def dupTPT2pub(self, dup_id):
        """ rotate duplex to the public frame (the=0, phi=0, rho=0) """
        resids = self.helices[dup_id]
        bps = self.resid2BPresid(resids)
        phase_resid = self.dup_phresids[dup_id]
        rev = self.isReverseImino(phase_resid, bps)
        self.residues = self.all_helix_residues
        #print "# duplex %3s  phase resid= %3d--%s " % (dup_id, phase_resid, self.seq_dict[phase_resid])
        #self.prnDpSeq(bps)
        duplex_xyz = self.getSegResXYZ_WoH_atom(resids)
        dup_axis = self.pp2axis(bps, vb=0)
        #self.analyzSlantAngle(dup_axis, resids)
        the, phi = self.axis2ThePhi(dup_axis)
        rho = self.thephi2rho(the, phi, phase_resid, rev,  self.seq_dict[phase_resid],  duplex_xyz)
        dup_pub_xyz = self.rotDup2Pub(the, phi, rho, duplex_xyz)
        dup_pub_xyz = self.transMove(dup_pub_xyz, dup_pub_xyz[resids[0]][0,:])
        
        #dup_axis = self.pp2axis(bps, dup_pub_xyz, vb=0)
        #the, phi = self.axis2ThePhi(dup_axis)
        #rho = self.thephi2rho(the, phi, phase_resid, rev, self.seq_dict[phase_resid], dup_pub_xyz)
        #print "# pub dup: the= %6.1f phi= %6.1f rho= %6.1f " % (the, phi, rho)
        ##dup_pub_xyz = self.anchor_by_firstP(dup_pub_xyz)
        #lns = self.resPdbDt2pdb(dup_pub_xyz, pdb_resids = self.pdb_resids)
        #self.saveLines(lns, "%s_dup_rpub_%s.pdb" % (fn_pre, dup_id), vb=1)
        return dup_pub_xyz
    
    
    def transMove(self, seg_xyz, trans):
        resids = seg_xyz.keys()
        for resid in resids:
            seg_xyz[resid] = seg_xyz[resid] - trans
        return seg_xyz 
    
    
    def isReverseImino(self, resid, bps):
        for _bp in bps:
            if resid in _bp:
                if resid == _bp[0]:
                    return 0
                else:
                    return 1
                
   
        
    def genPubDupPDB(self, save=0):
        self.all_dup_pub_xyz = {}
        for dup_id in self.duplex_ids:
            self.all_dup_pub_xyz[dup_id] = self.dupTPT2pub(dup_id)
        if save:
            fn_pub_dup = "%s_pub_dup.bin" % self.fn_pre
            self.dumpObj(self.all_dup_pub_xyz, fn_pub_dup)
    
    
    def savePubDupPDB(self):
        for dup_id in self.duplex_ids:
            dup_xyz = self.all_dup_pub_xyz[dup_id]
            lns =  self.resXYZ2pdb(pdb_xyz=dup_xyz)
            self.saveLines(lns, "%s_pub_dup%s.pdb" % (self.fn_pre, dup_id), vb=0)
    
            
    def genDup4Orients(self, vb=0):
        self.dup4dros = {}
        print self.dup_orients 
        for dup_id in self.duplex_ids:
            the, phi, rho, delta = self.dup_orients[dup_id]
            self.dup4dros[dup_id] = self.gen4ThePhiRho(the, phi, rho, vb=vb)
        #print self.dup_orients 
        
        
    
            
    def loadDulexPDB(self, vb=0): # pass
        """ read duplex's pdb files """
        self.all_helix_residues = {}
        self.helices = {}
        for dup_id in self.duplex_ids:
            fn_dup = "%s_duplex%s.pdb" % (self.fn_pre, dup_id) 
            self.pdb_lns = self.readTxtFile(fn_dup, vb=vb)
            self.simpleParsePDB()
            self.all_helix_residues.update(self.residues)
            self.helices[dup_id] = self.resids
        if vb: 
            print "# self.all_helix_residues.key(): \n", self.all_helix_residues.keys()
            print "# self.helices ", self.helices
        
    
    def genDupDRO_ids(self, vb=0):
        num_dup  = len(self.duplex_ids)
        num_total_dros = pow(4, num_dup -1 )
        total_dros = [[0] * num_dup ]
        tmp = [0] * num_dup 
        for i in range(1, num_total_dros):
            cnt = -1
            while True:
                tmp[cnt] += 1
                if tmp[cnt] == 4:
                    tmp[cnt] = 0
                    cnt -= 1 
                else:
                    break 
            total_dros.append([j for j in tmp])
        if vb: 
            print "# num_duplex: %3d " % num_dup 
            print "# num_totla_dros: ", num_total_dros     
            print "# total_dros: ", total_dros
        self.tot_dros = total_dros 
        self.num_tot_dros = num_total_dros 
        

class LINK(Rotation, RDCDEG, HomoMove):
    """ all functions for linking building blocks """
    z_axis = numpy.array([0., 0., 1.])
    y_axis = numpy.array([0.,1., 0.])
    def resid2BPresid(self, bp_resids,vb=0):
        """ input a duplex resids  and return basepair resids """
        num_bp = len(bp_resids)/2 
        bps = []
        for i in range(num_bp):
            bps.append((bp_resids[i], bp_resids[-(i+1)]))
        if vb: print "# bps: ", bps 
        return bps 
    
    def prnDpSeq(self, bps):
        """ print duplex sequence """
        for _bp in bps:
            print "%3d--%3d  %s%s " % (_bp[0], _bp[1], self.seq_dict[_bp[0]], self.seq_dict[_bp[1]])
            
            
    def getResRefFrame(self, res_xyz): 
        """ return residue reference frame: A,G: P, C4',C1',N9,  C,U:, P, C4',C1',N1 """
        ref_xyz = []
        for atom_id in [0,5,7,8]:
            ref_xyz.append(res_xyz[atom_id,:].copy())
        return array(ref_xyz)
   
    def getSegRefFrame(self, seg_xyz, ref_resids):
        """ return reference frame with more than 1 nucleotides """
        ref_xyz = []
        for resid in ref_resids:
            res_xyz = seg_xyz[resid]
            #print res_xyz
            if type(res_xyz) == type([]):
                res_xyz = array(res_xyz)
            for atom_id in [0,5,7,8]:    
            #for atom_id in range(len(res_xyz)): #range(9):#[0,5,7,8]:
                #print "# atom_id: ", atom_id
                ref_xyz.append(res_xyz[atom_id,:].copy())
        return array(ref_xyz)        
            

    def getResC5O5P(self, res_xyz):
        """ return residue C5', O5', P' from residue's XYZ 
        # resname: A
        # id=   4  atom= C5'
        # id=   3  atom= O5'
        # id=   0  atom= P
        # id=  21  atom= O3'
        # resname: G
        # id=   4  atom= C5'
        # id=   3  atom= O5'
        # id=   0  atom= P
        # id=  22  atom= O3'
        # resname: C
        # id=   4  atom= C5'
        # id=   3  atom= O5'
        # id=   0  atom= P
        # id=  19  atom= O3'
        # resname: U
        # id=   4  atom= C5'
        # id=   3  atom= O5'
        # id=   0  atom= P
        # id=  19  atom= O3'
        """
        return res_xyz[4,:].copy(), res_xyz[3,:].copy(), res_xyz[0,:].copy()
    
    
    def getResP_OP1_OP2(self, res_xyz):
        """ return residue P, OP1, OP2 from residue's XYZ 
        # resname: A# resname: A
        # id=   0  atom= P
        # id=   1  atom= OP1
        # id=   2  atom= OP2
        """
        return res_xyz[0,:].copy(),  res_xyz[1,:].copy(),  res_xyz[2,:].copy(), res_xyz[3,:].copy()

    
    def getResP(self, res_xyz):
        """ return residue P, OP1, OP2 from residue's XYZ 
        # id=   0  atom= P
        """
        return res_xyz[0,:].copy()
    
    
    def predictO3(self, res_xyz):
        """ predict O3'(n+1) from P(n) tetralhedral geometry """
        P, OP1, OP2, O5P = self.getResP_OP1_OP2(res_xyz)
        OPP = (OP1 - P) + (OP2 - P)
        OPP = self.normVect(OPP)
        P_O5 = self.normVect(O5P - P)
        rot_axis = self.cross(P_O5, OPP)
        P_O3  = 1.608 * self.new_arb_rot(P_O5, rot_axis, -105.0)
        #self.prnVec(P, "P_pos")
        #self.prnVec(OP1, "OP1")
        #self.prnVec(OP2, "OP2")
        #self.prnVec(O5P, "O5")
        #self.prnVec(P_O3, "P_O3")
        O3 = P_O3 + P 
        #self.prnVec(P, "P")
        #self.prnVec(O3, "O3")
        return O3 
        #P_O3 = 1.608 * self.new_arb_rot(nP_O3, nO5_P, 180.0) # alpha: -60, 180, 60 
        
    def getResO3(self, res_xyz, res_nt):
        """ return residue C5', O5', P' from residue's XYZ 
        # resname: A
        # id=   4  atom= C5'
        # id=   3  atom= O5'
        # id=   0  atom= P
        # id=  21  atom= O3'
        # resname: G
        # id=   4  atom= C5'
        # id=   3  atom= O5'
        # id=   0  atom= P
        # id=  22  atom= O3'
        # resname: C
        # id=   4  atom= C5'
        # id=   3  atom= O5'
        # id=   0  atom= P
        # id=  19  atom= O3'
        # resname: U
        # id=   4  atom= C5'
        # id=   3  atom= O5'
        # id=   0  atom= P
        # id=  19  atom= O3'
        """
        if res_nt == 'A':
            aid = 21
        elif res_nt == 'G':
            aid = 22
        elif res_nt in ['U',"C"]:
            aid = 19 
        else:
            raise "# %s not in A,U,G,C" % res_nt 
        return res_xyz[aid,:].copy()
    


class PackRNA(INPUT, RNAPDB, LINKORDER, DUPLEX3D, LINK, BlockOther):
    def __init__(self):
        DUPLEX3D.__init__(self)
        self.new_pdb_xyz = None     
        self.link_order = None 
        self.link_resids = None
    
    def prn2ndRNA_in_mem(self):
        self.prnDict(self.sec_inf)
            
    
    def addTloopRDC(self, _seg):
        """ link both ends of Tloop by rotating on the P(lowest resid of terminal loop) """
        rna_seg, seq_dt  = self.rna_segs[_seg[0]][_seg[1]]
        #print rna_seg 
        self.seq_dict.update(seq_dt)
        seg_xyz = rna_seg[2]
        resids = seg_xyz.keys()
        resids.sort()
        link_ref = self.getResRefFrame(rna_seg[1][0]) 
        if self.anchor_ref is None:
            anchor_ref = self.getResRefFrame(self.new_pdb_xyz[resids[0]-1])
        else:
            anchor_ref = self.anchor_ref 
        HM = self.getHM(anchor_ref, link_ref, is_center=True, vb=1)
        for resid in resids:
            self.new_pdb_xyz[resid] = self.homoMove(seg_xyz[resid], HM)  # center2 is the center of ref_frame, 
        Pivot = self.getResP(self.new_pdb_xyz[resids[0]])
        
        resid = resids[-1] # last resid in tloop 
        O3Pd = self.getResO3(self.new_pdb_xyz[resid], self.seq_dict[resid])
        O3Pt = self.predictO3(self.new_pdb_xyz[resid+1])
        
        P_O3Pd = O3Pd - Pivot
        P_O3Pt = O3Pt - Pivot
        ang = self.getVecAngle(P_O3Pd, P_O3Pt)
        #print "# ang= %6.1f " % degrees(ang)
        rot_axis = self.cross(P_O3Pt, P_O3Pd)
        ang_deg = degrees(ang)
        #new_pivot =  self.new_arb_rot(Pivot, rot_axis, ang_deg)
        po3d = self.getVecLen(P_O3Pd)
        po3t = self.getVecLen(P_O3Pt)
        ratio = po3t/po3d
        #print "# po3d= %8.3f po3t= %8.3f ratio(po3t/po3d)= %8.3f  " % (po3d, po3t, ratio)
        #self.prnVec(O3Pd, cmt="O3Pd")
        #self.prnVec(O3Pt, cmt='O3pt')
        new_pv =  self.new_arb_rot(P_O3Pd, rot_axis, -ang_deg)
        new_Po3 = new_pv*ratio + Pivot
        self.prnVec(new_Po3, cmt="new_po3")
        for resid in resids:
            res_xyz = self.new_pdb_xyz[resid].copy()
            num_atom, dim = res_xyz.shape 
            new_res_xyz = []
            for i in range(num_atom):
                atom_xyz = res_xyz[i,:]
                pa = atom_xyz - Pivot
                new_pv =  self.new_arb_rot(pa, rot_axis, -ang_deg)
                #if fabs(pa[0]) < 0.01: 
                    #self.prnVec(pa, cmt='pa') #print pa 
                    #self.prnVec(new_pv, cmt='new_pv')
                new_atom_xyz  = new_pv*ratio + Pivot
                new_res_xyz.append(new_atom_xyz)
            #print new_res_xyz
            self.new_pdb_xyz[resid] =  array(new_res_xyz) 
            #print "# resid: ", resid 
        self.anchor_ref =  None #self.getResRefFrame(self.homoMove(rna_seg[1][1], ROTM, center1, center2))
    
    
    
    
            
    def prnResArrAids(self):
        """ print residue XYZ coordinates' array atom ids """
        #seek_atoms = ["C5'","O5'", "P", "O3'"]
        seek_atoms = ['P','OP1','OP2']
        for res_nt in ['A','G',"C","U"]:
            _nt_atoms = nt_atoms_woh[res_nt]
            print "# resname: %s " % res_nt 
            #for i in range(len(_nt_atoms)):
            #    print "# id= %3d  atom= %s  " % (i, _nt_atoms[i])
            aids = []
            for _satom in seek_atoms:
                aids.append(_nt_atoms.index(_satom))
            for i in aids:
                print "# id= %3d  atom= %s  " % (i, _nt_atoms[i])    


    
    def orient2pdb(self, orient):
        """ pack duplexes to pdb with different orientations """
        self.new_pdb_xyz = None 
        for seg_id  in range(self.num_lnks): #link_order:
            _seg = self.link_order[seg_id]
            print "# linking seg: ", _seg 
            if self.new_pdb_xyz is None:
                self.setAnchorSegRDC(_seg)
            elif _seg[0] == 'linker':
                self.addOss(_seg)
            elif _seg[0] == 'helix':
                if _seg[1] == 1: orient_id = orient[0]   #h2_rot
                elif _seg[1] == 2: orient_id = orient[1] # h3_rot
                self.addHelixRDC(_seg, orient_id=orient_id)
            elif _seg[0] == 'tloop':
                self.addTloopRDC(_seg)
        lns =  self.resDict2pdb()
        return lns 

    def correctFinalLnk(self, _seg):
        """ correct last segment orientation  """
        rna_seg = self.rna_segs[_seg[0]][_seg[1]]
        seg_xyz = rna_seg[2]
        resids = seg_xyz.keys()
        resids.sort()
        seg_first_resid, seg_last_resid = resids[0], resids[-1] 
        Pivot = self.getResP(self.new_pdb_xyz[resids[0]])
        resid = resids[-1] # last resid in tloop 
        O3Pd = self.getResO3(self.new_pdb_xyz[resid], self.seq_dict[resid])
        O3Pt = self.predictO3(self.new_pdb_xyz[resid+1])
        
        P_O3Pd = O3Pd - Pivot
        P_O3Pt = O3Pt - Pivot
        ang = self.getVecAngle(P_O3Pd, P_O3Pt)
        #print "# ang= %6.1f " % degrees(ang)
        rot_axis = self.cross(P_O3Pt, P_O3Pd)
        ang_deg = degrees(ang)
        #new_pivot =  self.new_arb_rot(Pivot, rot_axis, ang_deg)
        po3d = self.getVecLen(P_O3Pd)
        po3t = self.getVecLen(P_O3Pt)
        ratio = po3t/po3d
        #print "# po3d= %8.3f po3t= %8.3f ratio(po3t/po3d)= %8.3f  " % (po3d, po3t, ratio)
        #self.prnVec(O3Pd, cmt="O3Pd")
        #self.prnVec(O3Pt, cmt='O3pt')
        new_pv =  self.new_arb_rot(P_O3Pd, rot_axis, -ang_deg)
        new_Po3 = new_pv*ratio + Pivot
        self.prnVec(new_Po3, cmt="new_po3")
        for resid in resids:
            res_xyz = self.new_pdb_xyz[resid].copy()
            num_atom, dim = res_xyz.shape 
            new_res_xyz = []
            for i in range(num_atom):
                atom_xyz = res_xyz[i,:]
                pa = atom_xyz - Pivot
                new_pv =  self.new_arb_rot(pa, rot_axis, -ang_deg)
                #if fabs(pa[0]) < 0.01: 
                    #self.prnVec(pa, cmt='pa') #print pa 
                    #self.prnVec(new_pv, cmt='new_pv')
                new_atom_xyz  = new_pv*ratio + Pivot
                new_res_xyz.append(new_atom_xyz)
            #print new_res_xyz
            self.new_pdb_xyz[resid] =  array(new_res_xyz) 
            #print "# resid: ", resid 
        self.anchor_ref =  None #self.getResRefFrame(self.homoMove(rna_seg[1][1], ROTM, center1, center2))



  
    
            
    def n_setAnchorSeg(self, _seg, vb=0):
        """ set anchor with only duplex rotated back to public axis """
        self.new_pdb_xyz = {}
        self.new_2nd_str = []
        if _seg[0] == 'helix':
            hex_id = _seg[1]
            if vb: print "# achor helix:  %s " % hex_id 
            pub_xyz = self.all_dup_pub_xyz[hex_id]
            #self.chkDupXYZDRO(hex_id, pub_xyz)
            _tpt = self.dup_orients[hex_id]
            if vb: print "# the, phi, rho: ", _tpt 
            _tpt_ref = (_tpt[0], _tpt[1], _tpt[2])
            tmp_xyz = self.rotPub2refTPR(_tpt_ref, pub_xyz)  # helix 1:  orientation 3 is correct orientation to crystall structure 
            #self.chkDupXYZDRO(hex_id, tmp_xyz)
            self.addSeg3D2PDB(tmp_xyz)
            if vb: print "# set duplex %s as anchor duplex" % hex_id 
        else:
            self.err_exit("# anchor segment must be a duplex ")
   
    
    def getRefResids(self, link_resids, pdb_resids, vb=0): #pass
        """ return resids of reference frame( one or two residues ) """
        ref_resids = []
        for resid in link_resids:
            if resid in pdb_resids:
                ref_resids.append(resid)
        if vb: 
            print "# get refresids: "
            print "# link resids:", link_resids
            print "# pdb resids: ", pdb_resids
            print "# refence resids: ", ref_resids
        return ref_resids
    
        
    def n_addHelix(self, _seg, orient_id=0):
        """ add O3' 3-> 5: using P tetralhedral geometry  """
        hex_id = _seg[1]
        #print "# add helix to pdb: %s " % hex_id
        dup_xyz = self.all_dup_pub_xyz[hex_id]
        if orient_id > 3:
            #print " # add by reference frame "
            self.addSeg3D2PDB(dup_xyz)
        else:
            #print "# add by bond link"
            tmp_xyz = self.rotPub2refTPR(self.dup4dros[hex_id][orient_id], dup_xyz)
            #self.chkDupXYZDRO(hex_id, tmp_xyz)
            #self.addSeg3D2PDB_BondLink(tmp_xyz)
            self.addSeg3D2PDB_BLP(tmp_xyz)
    
    
    def addSeg3D2PDB_BondLink(self, seg_xyz, vb=0):
        pdb_resids = self.new_pdb_xyz.keys()
        resids = seg_xyz.keys()
        resids.sort()
        seg_first_resid = resids[0]
        seg_last_resid = resids[-1]
        anchor_resid = seg_first_resid - 1
        lnk_resid = seg_first_resid
        if anchor_resid not in pdb_resids:
            anchor_resid = seg_last_resid + 1
            lnk_resid = seg_last_resid
            if anchor_resid not in pdb_resids: 
                raise "# could not find anchor resid for seg: %s %d  "   % _seg
        print "# add by bond linkage: "
        if anchor_resid < lnk_resid: # link 5'--> 3'
            print "# 5'--3': achor resid: %3d  lnk_resid: %3d  " % (anchor_resid, lnk_resid)
            ref_P = self.new_pdb_xyz[lnk_resid][0]
            lnk_P = seg_xyz[lnk_resid][0]
            trans   = ref_P - lnk_P
            ##print seg_xyz.keys()
            ##print self.new_pdb_xyz.keys()
            #self.prnVec(ref_P, cmt='ref_P')
            #self.prnVec(lnk_P, cmt='lnk_P')
            #self.prnVec(trans, cmt='tran')
            
        else: # 3' -- > 5' 
            ref_O3 = self.getResO3(self.new_pdb_xyz[lnk_resid], self.seq_dict[lnk_resid]) #self.predictO3(self.new_pdb_xyz[anchor_resid])
            lnk_O3 = self.getResO3(seg_xyz[lnk_resid], self.seq_dict[lnk_resid])
            trans = ref_O3 - lnk_O3
        for resid in resids:
            self.new_pdb_xyz[resid] = seg_xyz[resid] + trans 
        
            
    def addSeg3D2PDB_BLP(self, seg_xyz, vb=0):    # pass 
        """ add segment 3D structure to pdb by bond linkage, only translation movement, no rotation
        add O3' 3-> 5: using P tetralhedral geometry  
        """
        pdb_resids = self.new_pdb_xyz.keys()
        resids = seg_xyz.keys()
        resids.sort()
        seg_first_resid = resids[0]
        seg_last_resid = resids[-1]
        anchor_resid = seg_first_resid - 1
        lnk_resid = seg_first_resid
        if anchor_resid not in pdb_resids:
            anchor_resid = seg_last_resid + 1
            lnk_resid = seg_last_resid
            if anchor_resid not in pdb_resids: 
                raise "# could not find anchor resid for seg: %s %d  "   % _seg
        #print "# add by bond linkage: "
        if anchor_resid < lnk_resid: # link 5'--> 3'
            #print "# 5'--3': achor resid: %3d  lnk_resid: %3d  " % (anchor_resid, lnk_resid)
            ref_O3 = self.getResO3(self.new_pdb_xyz[anchor_resid], self.seq_dict[anchor_resid])
            lnk_O3 = self.predictO3(seg_xyz[lnk_resid])
            trans_O3 = ref_O3 - lnk_O3
            #self.prnVec(ref_O3, cmt='ref_o3')
            #self.prnVec(lnk_O3, cmt='lnk_o3')
            #self.prnVec(trans_O3, cmt='trans_o3')
            
        else: # 3' -- > 5' 
            ref_O3 = self.predictO3(self.new_pdb_xyz[anchor_resid])
            lnk_O3 = self.getResO3(seg_xyz[lnk_resid], self.seq_dict[lnk_resid])
            trans_O3 = ref_O3 - lnk_O3
        for resid in resids:
            self.new_pdb_xyz[resid] = seg_xyz[resid] + trans_O3 
        
            
    def addSeg3D2PDB(self, seg_xyz, LINK=99):
        """ rotation and traslation """
        resids = seg_xyz.keys()
        pdb_resids = self.new_pdb_xyz.keys()
        ref_resids = self.getRefResids(resids, pdb_resids, vb=0)
        if len(ref_resids) > 0:  #fixed link point
            if LINK < 2: ref_resids = [ref_resids[LINK]]
            link_ref = self.getSegRefFrame(seg_xyz, ref_resids) 
            anchor_ref = self.getSegRefFrame( self.new_pdb_xyz, ref_resids)
            hm = self.tgetHM(anchor_ref, link_ref, is_center=True, vb=0)
            for resid in resids:
                if resid in pdb_resids: continue
                #print "# HM add %s to pdb" % resid
                self.new_pdb_xyz[resid] = self.homoMove(seg_xyz[resid], hm)  # center2 is the center of ref_frame, 
        else:
            for resid in resids:
                if resid in pdb_resids: continue
                #print "# directly add %s to pdb" % resid
                self.new_pdb_xyz[resid] =  seg_xyz[resid].copy() #tmp_xyz[resid] #self.rotXYZ_by_refTPR(rot, hex_xyz[resid]) #hex_xyz[resid]
        

    
    def n_addBulge(self, _seg):
        rna_seg = self.rna_segs[_seg[0]]
        #print rna_seg 
        _bg = rna_seg[_seg[1]] #  ['A', 'AA', 'GU'], [114, 206, 115, 116, 205], {114: array([[ 30.574,  40.001,  64.725],
        #print "bugle: ", _bg
        _bg_xyz = _bg[2]
        _bg_resids = _bg[1]
        pdb_resids = self.new_pdb_xyz.keys()
        if _bg_resids[0] in pdb_resids and _bg_resids[1] in pdb_resids:
            ref_resids = _bg_resids[0], _bg_resids[1]
        elif _bg_resids[3] in pdb_resids and _bg_resids[4] in pdb_resids:    
            ref_resids = _bg_resids[3], _bg_resids[4]
        else:
            print "# bg: ", _bg[0], _bg[1]
            print "# current pdb_resids: ", pdb_resids
            raise " # could not find link ref for the bulge: "
        
        link_ref = self.getSegRefFrame( _bg_xyz, ref_resids) #RefXYZ() ResRefFrame(rna_seg[1][0]) 
        anchor_ref = self.getSegRefFrame( self.new_pdb_xyz, ref_resids)
        
        hm = self.getHM(anchor_ref, link_ref, is_center=True, vb=1)
        for resid in _bg_resids:
            #print "# add %s to pdb" % resid
            self.new_pdb_xyz[resid] = self.homoMove(_bg_xyz[resid], hm)  # center2 is the center of ref_frame, 
        #self.anchor_ref =  None #self.getResRefFrame(self.homoMove(rna_seg[1][1], ROTM, center1, center2))
    
    def n_addTloop(self, _seg):
        obk_tpt = self.tpt_tloops[_seg[1]]
        seg_xyz = obk_tpt[1]
        self.addSeg3D2PDB(seg_xyz)
        
    def n_addLinker(self, _seg):
        #obk_2str  = self.sec_inf[_seg[0]][_seg[1]]
        obk_tpt = self.tpt_linkers[_seg[1]]
        seg_xyz = obk_tpt[1]
        self.addSeg3D2PDB(seg_xyz, LINK=0)
        
    
        
    def link_Duplex(self, link_order=None, inc_dro_ids=None):
        """ link two duplexes, one duplex orientation is fixed, the other one has four orientations """
        if link_order is None: link_order = self.link_order
        num_segs = len(link_order) 
        dup_pre = ''
        for _seg in link_order:
            if _seg[0] == 'helix':
                dup_pre +="H%s" % _seg[1]
        for i in range(4):
            self.add2log(["# Linking two duplexes: %s - %s "  % (link_order[0], link_order[-1]), "# duplex: %s dro_id: %s " % (link_order[-1], i)])
            self.new_pdb_xyz = None 
            for seg_id  in range(num_segs): #link_order:
                _seg = link_order[seg_id]
                if self.new_pdb_xyz is None:
                    self.n_setAnchorSeg(_seg, vb=0)
                elif _seg[0] == 'linker':
                    self.n_addLinker(_seg)
                elif _seg[0] == 'helix':
                    self.n_addHelix(_seg, orient_id=i)
                elif _seg[0] == 'tloop':
                    self.n_addTloop(_seg)
                elif _seg[0] == 'bulge':
                    self.n_addBulge(_seg)
            #for dup_id in self.duplex_ids:
            #    self.chkDupXYZDRO(dup_id, self.new_pdb_xyz)
            fn_mid_key = "h%s_r%s" % (link_order[-1][1],i)
            if hasattr(self, 'dup_names_in_paper'):
                fn_mid = self.dup_names_in_paper[fn_mid_key]
            else:
                fn_mid = fn_mid_key
            lns =  self.resXYZ2pdb()
            self.saveLines(lns, "%s_%s_%s.pdb" % (self.out_pdb_pre, dup_pre, fn_mid), vb=1)
        
    def isExDROIDs(self, dro_ids, link_order):
        if not hasattr(self, 'allowed_dup_dro_ids'): 
            #print "# no self.allowed_dup_dro_ids"
            return False
        adp_keys = self.allowed_dup_dro_ids.keys()
        for _seg in link_order:
            if _seg[0] == 'helix':
                dup_key = 'h%s' % _seg[1]
                orient_id= dro_ids[int(_seg[1])-1]
                if dup_key in adp_keys and orient_id not in self.allowed_dup_dro_ids[dup_key]:
                    return True
        return False
                
    
    
    def pdbFnmid(self, i, dro_ids, test=0):
        if self.fn_pre == 'riboAs' and not test:
            h2_key = 'h2_r%s' % dro_ids[1]
            h3_key = 'h3_r%s' % dro_ids[2]
            fn_mid= "%s%s" % (self.dup_names_in_paper[h2_key], self.dup_names_in_paper[h3_key])
        else:
            fn_mid = "dro%s" % i 
        return fn_mid
    
    
    def build_RNA(self, link_order=None, inc_dro_ids=None):
        """ pack duplexes to pdb with different orientations """
        if link_order is None: link_order = self.link_order
        num_segs = len(link_order)    
        if not hasattr(self, 'tot_dros'):
            self.tot_dros = [[0] * len(self.duplex_ids)]
            self.num_tot_dros = 1
        for i in range(self.num_tot_dros):
            if inc_dro_ids is not None and i not in inc_dro_ids: continue
            dro_ids = self.tot_dros[i]
            if self.isExDROIDs(dro_ids, link_order): continue 
            self.add2log(["# packing duplexes:", "# duplex: %s dro_id: %s " % (i, dro_ids)])
            self.new_pdb_xyz = None 
            for seg_id  in range(num_segs): #link_order:
                _seg = link_order[seg_id]
                if self.new_pdb_xyz is None:
                    self.n_setAnchorSeg(_seg, vb=0)
                elif _seg[0] == 'linker':
                    self.n_addLinker(_seg)
                elif _seg[0] == 'helix':
                    self.n_addHelix(_seg, orient_id= dro_ids[int(_seg[1])-1])
                elif _seg[0] == 'tloop':
                    self.n_addTloop(_seg)
                elif _seg[0] == 'bulge':
                    self.n_addBulge(_seg)
            for dup_id in self.duplex_ids:
                self.chkDupXYZDRO(dup_id, self.new_pdb_xyz)        
            lns =  self.resXYZ2pdb()
            fn_mid = self.pdbFnmid(i, dro_ids)
            self.saveLines(lns, "%s_%s.pdb" % (self.out_pdb_pre, fn_mid), vb=1)
        #self.correctFinalLnk(link_order[-1]) 
    
    
    def rotAllDup2Pub(self):
        self.genPubDupPDB()
        self.all_hex_resids = self.helices
         
    
    def loadLnkPDB(self, vb=0): # pass
        """ read duplex's pdb files """
        mem_residues = self.residues
        for lnk_id_str, fn_pdb in self.linker_pdbs:#lnk_id in self.link_ids:
            #fn_pdb = '%s_lnk%s.pdb' %  (self.fn_pre, lnk_id) 
            self.pdb_lns = self.readTxtFile(fn_pdb, vb=vb)
            self.simpleParsePDB()
            seg_xyz = self.getSegResXYZ_WoH_dict(self.resids)
            lnk_resids = (self.resids[0], self.resids[-1])
            self.tpt_linkers[int(lnk_id_str[1:])] = (lnk_resids, seg_xyz)
        self.residues = mem_residues    
    
    
    def loadTloopPDB(self, vb=0): # pass
        """ read duplex's pdb files """
        mem_residues = self.residues
        for seg_id_str, fn_pdb in self.tloop_pdbs: 
            #fn_pdb = '%s_tloop%s.pdb' %  (self.fn_pre, seg_id) 
            self.pdb_lns = self.readTxtFile(fn_pdb, vb=vb)
            self.simpleParsePDB()
            seg_xyz = self.getSegResXYZ_WoH_dict(self.resids)
            lnk_resids = (self.resids[0], self.resids[-1])
            #print "# lnk resids: ",  lnk_resids
            self.tpt_tloops[int(seg_id_str[1:])] = (lnk_resids, seg_xyz)
        
        self.residues = mem_residues    
    
    def getLinkOrder(self, lnk_str):
        lnks = lnk_str.split('-')
        link_oder = []
        for _lnk in lnks:
            link_oder.append(self.motifs[_lnk])
        return link_oder 
    
    def setAllowedDupDRO(self, allowed_droids):
        self.allowed_dup_dro_ids = {}
        dro_ids = allowed_droids.split(';')
        for _droids in dro_ids:
            _dro_ids = _droids.strip()
            if _dro_ids and ':' in _dro_ids:
                dup_id, dup_droids = _dro_ids.split(':')
                dup_droids = dup_droids.split(',') 
                dup_droids = [ int(_droid) for _droid in dup_droids]
                self.allowed_dup_dro_ids[dup_id] = dup_droids 
                    
        
    def stepPack(self):
        self.STEP = 0
        if hasattr(self, 'all_link_orders'):
            for lnk_oder_str, allowed_droids in self.all_link_orders:
                if allowed_droids:
                    self.setAllowedDupDRO(allowed_droids)
                if lnk_oder_str:
                    link_order = self.getLinkOrder(lnk_oder_str)
                    self.link_Duplex(link_order=link_order)    
                    if not self.STEP: self.STEP = 1
                
    
    
    def pack(self):
        self.genLinker(vb=0)
        self.genTloop()
        if hasattr(self, 'LOADLINKERPDB') and  self.LOADLINKERPDB:
            self.loadLnkPDB()
        if hasattr(self, 'LOADTLPPDB') and self.LOADTLPPDB:
            self.loadTloopPDB()
        self.rotAllDup2Pub()
        self.stepPack()
        if not self.STEP:
            self.genDupDRO_ids(vb=0)    
            self.build_RNA()

            
    def loadPackInp(self, fn_inp):
        lns = self.readTxtFile(fn_inp)
        seg_lns = self.getSegLns(lns)
        self.parsePackInput(seg_lns, vb=0)
        if hasattr(self, 'fn_seq_inp') and self.fn_seq_inp is not None:
            self.load2str(fn_inp = self.fn_seq_inp) 
            
        
    def parsePackInput(self, all_seg_lns, vb=0):
        for seg_lns in all_seg_lns:
            if seg_lns:
                first_ln = seg_lns[0]
                if 'input' in first_ln and 'sequence' in first_ln:
                    #self.fn_seq_inp =
                    name, prop = self.getInpNameProp(seg_lns[1], vb=vb)
                    if name == 'fn_seq': self.fn_seq_inp = prop 
                elif 'load' in first_ln:
                    if 'linker' in first_ln:
                        self.linker_pdbs = self.parseLoadPDB(seg_lns, vb=vb)
                    elif 'hairpin' in first_ln:
                        self.tloop_pdbs = self.parseLoadPDB(seg_lns, vb=vb)
                    else:
                        print "# non parsed load files"
                        self.prnLines(seg_lns)
                elif 'Orientations' in first_ln:
                    name, prop = self.getInpNameProp(seg_lns[1], vb=vb)
                    setattr(self, name, prop)
                    seg_lns = self.readTxtFile(prop)
                    self.prnLines(seg_lns)
                    self.parseDupOrients(seg_lns, vb=vb)
                elif 'control' in first_ln:
                    self.parseControl(seg_lns)
                else:
                    self.prnLines(["#******       None Formated text in the input file ****"]+ seg_lns)
            else:
                self.err_exit("# input block text format error, please check input file format ")
     
        
    def parseLoadPDB(self, seg_lns, vb=0):
        name, prop = self.getInpNameProp(seg_lns[1])
        setattr(self, name, prop)
        _pdbs = []
        for _ln in seg_lns[2:]:
            _pdbs.append(self.getInpNameProp(_ln))
        if vb: 
            print "# self.%s: %s " % (name, prop)
            print _pdbs 
        return _pdbs 
    
                
    def parseControl(self, seg_lns):
        self.prnLines(seg_lns)
        self.all_link_orders = []
        for ln in seg_lns:
            ln = ln.strip()
            if ln[0] != "#":
                tmp = split('\s+=\s+', ln)
                if tmp[0] == 'LOADLINKERPDB':
                    if tmp[1] in ['1','ON']:
                        self.LOADLINKERPDB = True
                        self.stack = 0
                    else:
                        self.LOADLINKERPDB = False
                elif tmp[0] == 'LOADTLPPDB':
                    if tmp[1] in ['1','ON']:
                        self.LOADTLPPDB = True
                        self.stack = 0
                    else:
                        self.LOADTLPPDB = False
                elif tmp[0] == 'link_ids':
                    self.link_ids = split('\s*,\s*', tmp[1])
                elif tmp[0] == 'tlp_ids':
                    self.tloop_ids = split('\s*,\s*', tmp[1])
                elif tmp[0] == 'link_order':
                    self.all_link_orders.append([tmp[1], ""])
                elif tmp[0] == 'dup_names_in_paper':
                    self.__parseDupNames(tmp[1])
                elif tmp[0] == 'allowed_dup_dro_ids':
                    if self.all_link_orders:
                        self.all_link_orders[-1][1] = tmp[1]
                    else:
                        self.all_link_orders.append(["",tmp[1]])
                else:
                    print "# nonparsed control: ", ln 
        
    
    def __parseDupNames(self, dup_names):
        #print dup_names
        """   self.dup_names_in_paper = {'h2_r0':'a','h2_r1':'d','h2_r2':'c','h2_r3':'b', 
                                  'h3_r0':'e','h3_r1':'f','h3_r2':'g','h3_r3':'h'} """
        nms = split('\s*,\s+', dup_names)
        self.dup_names_in_paper = {}
        for _nm in nms:
            _nm = _nm.strip()
            tmp = _nm.split(':')
            key = "%s_%s" % (tmp[0][:2], tmp[0][2:])
            self.dup_names_in_paper[key] = tmp[1]
        #self.prnDict(self.dup_names_in_paper)
        #print self.dup_names_in_paper
        
     
    def parseDupOrients(self, lns, vb=0):
        """ parse duplex's orientation
        {1: '[20.9 165.2 122.0 98.4]', 2: '[21.1 50.3 253.2 107.1]', 3: '[158.6 55.1 61.7 97.0 ]'}
        """
        self.dup_orients = {}
        self.dup_phresids = {}
        for ln in lns:
            ln = ln.strip()
            if not ln or ln[0] == "#": continue 
            tmp = split("\s+", ln)
            if len(tmp) == 6:
                dup_orients = [ float(_val) for _val in tmp[1:-1]]
                dup_id = tmp[0]
                self.dup_orients[dup_id] = dup_orients
                self.dup_phresids[dup_id] = int(tmp[-1])
        if vb:
            print "# duplex orientations:(self.dup_orients) "
            self.prnDict(self.dup_orients)
            print "# self.dup_phresids:(phase_resids): ", self.dup_phresids
        self.add2log(lns)
    
        