path: root/cad/src/outtakes/BuildDuplex_old/DnaGenerator.py

# Copyright 2006-2008 Nanorex, Inc.  See LICENSE file for details. 
"""
DnaGenerator.py

WARNING: this file has been mostly superseded by DnaDuplexGenerator.py.

@author: Will Ware
@version: $Id$
@copyright: 2006-2008 Nanorex, Inc.  See LICENSE file for details.

History:

Jeff 2007-06-13:
- Moved Dna class (and subclasses) to file Dna.py
  [subsequently renamed DnaGenHelper.py].

Mark 2007-08-23:
- Heavily restructured and cleaned up.

Mark 2007-10-18: 
- Did a major rewrite of this module, superseding it --
  DnaDuplexGenerator.py.
"""

# To do:
# 1) Use two endpoints to create an arbitrarily positioned duplex.
# 2) Remove support for Atomistic DNA models.

import foundation.env as env
import random

from utilities.constants import darkred, blue, lightgray

from foundation.Group import Group
from utilities.Log  import redmsg, greenmsg ##, orangemsg
from geometry.VQT import Q, V, angleBetween, cross, vlen, Veq
from math           import pi
from model.chunk import Chunk
from utilities.constants import gensym    

from dna.model.Dna_Constants import basesDict, getReverseSequence
from dna.commands.BuildDuplex_old.DnaGenHelper import B_Dna_PAM3, B_Dna_PAM5
from dna.commands.BuildDuplex_old.DnaGenHelper import basepath_ok

from command_support.GeneratorBaseClass import GeneratorBaseClass
from utilities.exception_classes import CadBug, PluginBug, UserError
from dna.commands.BuildDuplex_old.DnaGeneratorPropertyManager import DnaGeneratorPropertyManager

############################################################################

# DnaGeneratorPropertyManager must come BEFORE GeneratorBaseClass in this list
class DnaGenerator(DnaGeneratorPropertyManager, GeneratorBaseClass):

    cmd              =  greenmsg("Build DNA: ")
    sponsor_keyword  =  'DNA'
    prefix           =  'Dna'   # used for gensym

    # Generators for DNA, nanotubes and graphene have their MT name 
    # generated (in GeneratorBaseClass) from the prefix.
    create_name_from_prefix  =  True 
    
    # pass window arg to constructor rather than use a global, wware 051103
    def __init__(self, win):
        DnaGeneratorPropertyManager.__init__(self)
        GeneratorBaseClass.__init__(self, win)
        self._random_data = []
    
    # ##################################################
    # How to build this kind of structure, along with
    # any necessary helper functions.

    def change_random_seed(self):
        if 0:
            print "change_random_seed() called."
        self._random_data  =  []

    def _random_data_for_index(self, inIndex):
        while len( self._random_data ) < (inIndex + 1):
            self._random_data.append( random.randrange(12) )
        return self._random_data[inIndex]

    def gather_parameters(self):
        """
        Return the parameters from the property manager UI.
        
        @return: All the parameters:
                 - dnaSequence
                 - dnaType
                 - basesPerTurn
                 - chunkOption
        @rtype:  tuple
        """
        if not basepath_ok:
            raise PluginBug("The cad/plugins/DNA directory is missing.")
        
        dnaModel = str(self.modelComboBox.currentText())
        dnaType  = str(self.conformationComboBox.currentText())

        assert dnaType in ('B-DNA')

        # Get bases per turn.
        basesPerTurnString = str(self.basesPerTurnComboBox.currentText())
        basesPerTurn = float(basesPerTurnString)

        chunkOption = str(self.createComboBox.currentText())
        resolve_random = False
            # Later this flag may depend on a new checkbox in that case;
            # for now it doesn't matter, since sequence info is 
            # discarded for reduced bases anyway.
                

        (dnaSequence, allKnown) = \
            self._getSequence( resolve_random = resolve_random)
        
        x1 = self.x1SpinBox.value()
        y1 = self.y1SpinBox.value()
        z1 = self.z1SpinBox.value()
        
        x2 = self.x2SpinBox.value()
        y2 = self.y2SpinBox.value()
        z2 = self.z2SpinBox.value()
        
        endpoint1 = V(x1, y1, z1)
        endpoint2 = V(x2, y2, z2)

        return (dnaSequence, 
                dnaModel,
                dnaType,
                basesPerTurn,
                chunkOption,
                endpoint1, 
                endpoint2)
    
    def checkParameters( self, inParams ):
        """
        Verify that the strand sequence contains no unknown/invalid bases.
        """
        theSequence, isValid  =  self._getSequence()
        
        return isValid

    # This is never called, which is good, because neither result, nor
    # choice is defined anywhere
    #def correctParameters( self, inParams):
    #    """
    #    Alert the user that the entered sequence is invalid. Give them
    #    some options for how to correct the sequence.
    #    """
    #    #theDialog  =  Ui_InvalidSequenceDialog()
    #    
    #    #optionsButtonGroup  =  theDialog.findChild( 'buttonbox_options' )
    #    #result  =  theDialog.exec()
    #    #choice  =  optionsButtonGroup.checkedid()
    #    
    #    if result == QDialog.Accepted:
    #        print 'choice: ', choice
    #
    #    return inParams
    
    def build_struct(self, name, params, position):
        """
        Build the DNA helix based on parameters in the UI.
        
        @param name: The name to assign the node in the model tree.
        @type  name: str
        
        @param params: The list of parameters gathered from the PM.
        @type  params: tuple
        
        @param position: The position in 3d model space at which to
                         create the DNA strand. This is always 0, 0, 0.
        @type position:  position
        """
        # No error checking in build_struct, do all your error
        # checking in gather_parameters
        theSequence, \
        dnaModel, \
        dnaType, \
        basesPerTurn, \
        chunkOption, \
        endpoint1, \
        endpoint2 = params
        
        if Veq(endpoint1, endpoint2):
            raise CadBug("Dna endpoints cannot be the same point.")
            return
        
        if len(theSequence) < 1:
            msg = redmsg("Enter a strand sequence to preview/insert DNA")
            self.MessageGroupBox.insertHtmlMessage(msg, setAsDefault=False)
            self.dna = None # Fixes bug 2530. Mark 2007-09-02
            return None

        if dnaModel == 'PAM3':
            dna = B_Dna_PAM3()
        else:
            dna = B_Dna_PAM5()
        
        self.dna = dna  # needed for done msg
        
        # Create the model tree group node. 
        rawDnaGroup = Group(self.name, 
                            self.win.assy,
                            self.win.assy.part.topnode)
        try:
            # Make the DNA duplex. <rawDnaGroup> returns a different 
            # grouping arrangement for atomistic vs. PAM5. This 'issue'
            # is resolved when we regroup the atoms into strand chunks
            # below.
            dna.make(rawDnaGroup, theSequence, basesPerTurn)
             
            self._orientRawDnaGroup(rawDnaGroup, endpoint1, endpoint2)
        
            # Now group the DNA atoms based on the grouping option selected
            # (i.e. "Strand chunks" or "Single Chunk").
            dnaGroup = self._makePAMStrandAndAxisChunks(rawDnaGroup)
            
            if chunkOption == 'Single chunk':
                return self._makeDuplexChunk(dnaGroup)
            
            return dnaGroup
        
        except (PluginBug, UserError):
            # Why do we need UserError here? Mark 2007-08-28
            rawDnaGroup.kill()
            raise PluginBug("Internal error while trying to create DNA duplex.")
            return None
        
    def _orientRawDnaGroup(self, rawDnaGroup, pt1, pt2):
        """
        Orients the raw DNA group based on two endpoints.
        
        @param rawDnaGroup: The raw DNA group created by make().
        @type  rawDnaGroup: L{Group}
        
        @param pt1: The first endpoint of the DNA strand.
        @type  pt1: L{V}
        
        @param pt2: The second endpoint of the DNA strand.
        @type  pt2: L{V}
        
        @attention: Only works for PAM5 models.
        """
        
        a = V(0.0, 0.0, -1.0)
        # <a> is the unit vector pointing down the center axis of the default
        # rawDnaGroup structure which is aligned along the Z axis.
        bLine = pt2 - pt1
        bLength = vlen(bLine)
        b = bLine/bLength
        # <b> is the unit vector parallel to the line (i.e. pt1, pt2).
        axis = cross(a, b)
        # <axis> is the axis of rotation.
        theta = angleBetween(a, b)
        # <theta> is the angle (in degress) to rotate about <axis>.
        scalar = self.dna.getBaseRise() * self.getSequenceLength() * 0.5
        rawOffset = b * scalar
        
        if 0: # Debugging code.
            print ""
            print "uVector  a = ", a
            print "uVector  b = ", b
            print "cross(a,b) =", axis
            print "theta      =", theta
            print "baserise   =", self.dna.getBaseRise()
            print "seqLength  =", self.getSequenceLength()
            print "scalar     =", scalar
            print "rawOffset  =", rawOffset
        
        if theta == 0.0 or theta == 180.0:
            axis = V(0, 1, 0)
            # print "Now cross(a,b) =", axis
            
        rot =  (pi / 180.0) * theta  # Convert to radians
        qrot = Q(axis, rot) # Quat for rotation delta.
        
        # Move and rotate the base chunks into final orientation.
        for m in rawDnaGroup.members:
            m.move(qrot.rot(m.center) - m.center + rawOffset + pt1)
            m.rot(qrot)

    def _getSequenceLength_OBS(self):
        """
        Returns the number of bases of the current sequence
        (from the Property Manager).
        
        @note: This is duplicated from DnaGeneratorPropert
        """
        (sequence, allKnown) = self._getSequence()
        return len(sequence)
    
    def _getSequence( self, 
                      reverse = False, 
                      complement = False, 
                      resolve_random = False, 
                      cdict = {} ):
        """
        Get the current DNA sequence from the Property Manager.
        
        This method is not fully private. It's used repeatedly to get the 
        same sequence when making the DNA (which means its return value 
        should be deterministic, even when making sequences with randomly 
        chosen bases [nim]), and it's also called from class 
        DnaGeneratorPropertyManager to return data to be stored back into the 
        Property Manager, for implementing the reverse and complement actions.
        (Ideally it would preserve whitespace and capitalization when used 
        that way, but it doesn't.)
        
        @param reverse: If true, returns the reverse sequence.
        @type  reverse: bool
        
        @param complement: If true, returns the complement sequence.
        @type  complement: bool
        
        @param resolve_random:
        @type  resolve_random: True
        
        @param cdict:
        @type  cdict: dictionary
        
        @return: (sequence, allKnown) where I{sequence} is a string in which 
                 each letter describes one base of the sequence currently
                 described by the UI, as modified by the passed reverse, 
                 complement, and resolve_random flags, and I{allKnown} is a 
                 boolean which says whether every base in the return value 
                 has a known identity.
        @rtype:  tuple
        
        @note: All punctuation/symbols are purged from the sequence and
               any bogus/unknown bases are substituted as 'N' (unknown).
        """
        
        sequence = ''
        allKnown = True
        
        cdict    =  basesDict
        
        # (Note: I think this code implies that it can no longer be a 
        #        number of bases. [bruce 070518 comment])
        currentSequence  =  str(self.getPlainSequence(inOmitSymbols = True))

        for ch in currentSequence:
            if ch in cdict.keys():  #'CGATN':
                properties = cdict[ch]
                if ch == 'N': ###e soon: or any other letter indicating a random base
                    if resolve_random: #bruce 070518 new feature
                        i    = len(sequence)
                        data = self._random_data_for_index(i) 
                        # a random int in range(12), in a lazily extended cache
                        ch   = list(cdict)[data%4]  
                        # modulus must agree with number of valid entries in 
                        # cdict.
                    else:
                        allKnown = False
                if complement:
                    try:
                        ch = properties['Complement']
                    except (KeyError):
                        raise KeyError("DNA dictionary doesn't have a \
                        'Complement' key for '%r'." % ch)
                        ch = 'N'
            elif ch in self.validSymbols: #'\ \t\r\n':
                ch  =  ''
            else:              
                allKnown  =  False

            sequence += ch

        if reverse: 
            sequence = getReverseSequence(sequence)
        
        return (sequence, allKnown)
    
    def _makeDuplexChunk(self, dnaGroup):
        """
        Returns a single DNA chunk given a dnaGroup containing multiple strand
        chunks.
        
        @param dnaGroup: The group object containing the DNA strand chunks.
        @type  dnaGroup: L{Group}
        
        @return: The DNA chunk.
        @rtype:  L{Chunk}
        """
        
        if not isinstance(dnaGroup.members[0], Chunk):
            env.history.message(redmsg(
                "Internal error in creating a single chunk DNA"))
            return
        
        for m in dnaGroup.members[1:]:
            if isinstance(m, Chunk):
                dnaGroup.members[0].merge(m)
                
        # Rename the merged chunk 
        dnaGroup.members[0].name = dnaGroup.name
        
        dnaChunk = dnaGroup.members[0]
        dnaChunk.setcolor(None)
        dnaGroup.ungroup()
        
        return dnaChunk
    
    def _getStrandName(self, strandNumber, numberOfBasesToDisplay = 0):
        """
        Returns a strand name given a strand number and the number of base
        letters to display in the name.
        
        @param strandNumber: The strand number, where:
                             - 0 = Strand1
                             - 1 = Strand2
                             - 2 = Axis (PAM5 only)
        @type  strandNumber: int
        
        @param numberOfBasesToDisplay: The number of base letters to display
                                       in the name. The default is 0.
        @type  numberOfBasesToDisplay: int
        
        @return: The strand name. (i.e. "StrandA:ATCG...")
        @rtype:  str
        """
        assert (strandNumber >= 0) and (strandNumber <= 2), \
               "strandNumber is %d. It can only be 0, 1, or 2." % strandNumber
        
        if strandNumber == 0:
            (sequence, allKnown) = self._getSequence()
            strandName = 'Strand1'
        elif strandNumber == 1:
            (sequence, allKnown) = self._getSequence(complement=True)
            strandName = 'Strand2'
        else:
            strandName = "Axis"
            
        if numberOfBasesToDisplay: 
            # Add strand letters to MT node name.
            numberOfLetters = min(len(sequence), numberOfBasesToDisplay)
            strandName += ":" + sequence[0:numberOfLetters]
        
            if len(sequence) > numberOfBasesToDisplay:
                # Add "..." if the sequence is longer than
                # <numberOfBasesToDisplay>.
                strandName += '...'
                          
        return strandName
    
    def _makePAMStrandAndAxisChunks(self, rawDnaGroup):
        """
        Returns a group containing the three strand chunks I{StrandA},
        I{StrandB} and I{Axis} of the current DNA sequence.
        
        @param rawDnaGroup: The raw Dna group which contains the 
                            base-pair chunks representing the sequence.
        @type  grawDnaGrouprp: L{Group}
        
        @return: The new DNA group that contains the three strand chunks
                 I{StrandA}, I{StrandB} and I{Axis}.
        @rtype:  L{Group}
        """
        
        startBasePair = rawDnaGroup.members[0]
        
        if not isinstance(startBasePair, Chunk):
            env.history.message(redmsg(
                "Internal error in creating a chunks for strands and axis"
                ))
            return rawDnaGroup
        
        # <startAtoms> are the PAM atoms that start StrandA, StrandB and Axis.
        # If the sequence is a single base, then we have 2 Pe atoms (one for 
        # strandA and one for StrandB.
        if self.getSequenceLength() == 1:
            startAtoms = ('Se3', 'Ae3', 'Pe5', 'Ae5')
        else:
            startAtoms = ('Se3', 'Sh3', 'Ae3', 'Pe5', 'Sh5', 'Ae5')
        
        Pe_count = 0
        tempList      =  []
        
        for atm in startBasePair.atoms.values():  
            if atm.element.symbol in startAtoms:                        
                tempList.append(atm)
                atomList = self.win.assy.getConnectedAtoms(tempList)
                tempList = []
                if atm.element.symbol in ('Se3', 'Pe5'):
                    Pe_count += 1
                    if Pe_count == 1:
                        strandAChunk = self._makeChunkFromAtomList(atomList)
                        strandAChunk.name = self._getStrandName(0)
                        first_Pe_found = True
                    else: # Pe_count == 2
                        # Only happens if the user entered a single letter
                        # for the sequence.
                        strandBChunk = self._makeChunkFromAtomList(atomList)
                        strandBChunk.name = self._getStrandName(1)
                elif atm.element.symbol in ('Sh3', 'Sh5'):
                    strandBChunk = self._makeChunkFromAtomList(atomList)
                    strandBChunk.name = self._getStrandName(1)
                elif atm.element.symbol in ('Ae3', 'Ae5'):
                    axisChunk = self._makeChunkFromAtomList(atomList)
                    axisChunk.name = self._getStrandName(2)
        
        # Assign default colors to strand and axis chunks.
        strandAChunk.setcolor(darkred)
        strandBChunk.setcolor(blue)
        axisChunk.setcolor(lightgray)
        
        # Place strand and axis chunks in this order: StrandA, StrandB, Axis.
        rawDnaGroup.addmember(strandAChunk)
        rawDnaGroup.addmember(strandBChunk)
        rawDnaGroup.addmember(axisChunk)
                
        self.win.win_update() # Needed?
                
        return rawDnaGroup
             
    def _makeChunkFromAtomList(self, atomList):
        """
        Creates a new chunk from the given atom list.
        
        @param atomList: List of atoms from which to create the chunk.
        @type  atomList: list
        
        @return: The new chunk.
        @rtype:  L{Chunk}
        
        @deprecated: use ops_rechunk.makeChunkFromAtomsList() instead.
        
        """
        if not atomList:
            print "bug in creating chunks from the given atom list"
            return
        
        newChunk = Chunk(self.win.assy, gensym("Chunk", self.win.assy))
        for a in atomList:            
            # leave the moved atoms picked, so still visible
            a.hopmol(newChunk)
        return newChunk   
    
    ###################################################
    # The done message

    def done_msg(self):
        
        if not self.dna: # Mark 2007-06-01
            return "No DNA added."

        return "Done creating a strand of %s." % (self.dna.form)