text: ###############################################################################
  #
  # pttableau.py - Object to represent protein tableaux and functions to
  #                parse output of TableauCreator program into tableau object.
  #
  #
  # File:    pttableau.py
  # Author:  Alex Stivala
  # Created: October 2007
  #
  #
  # $Id: pttableau.py 2703 2009-07-27 06:01:05Z astivala $
  #
  ###############################################################################
  """
  Ths module contains routines to generate protein tableaux using the
  axes fitted to SSEs by functions in the ptnode.py module, and the relative
  angle calculate in that module.
  
  This module also contains functions to parse the output of
  Arun Konargurthu's TableauCreator program, which was used before
  tableaux were re-implemented internally to this code.
  NOTE: not yet published or available (as of October 2007).
  IMPORTANT: this requires my modified versin of TableauCreator,
  and also a patched version of Bio.PDB file PDBIO.py :
  context diff for patching it is Bio.PDB.PDBIO.py.diff
  This patch was made relative to BioPython release 1.43.
  
  Tableaux are described by Kamat and Lesk 2007
  'Contact Patterns Between Helices and Strands of Sheet Define Protein
   Folding Patterns' Proteins 66:869-876
  and Lesk 2003 'From Electrons to Proteins and Back Again'
  Int. J. Quant. Chem. 95:678-682
  and Lesk 1995 'Systematic representation of folding patterns'
  J. Mol. Graph. 13:159-164.
  
  Two classes are provided, PTTableau and PTTableauPacked. The latter is
  a more compact format based on LAPACK style symmetric matrix packed
  array storage, useful for holding whole databse of tableaux in memory
  (see buildtableauxdb.py) and dumping/loading it.  They can be used
  interchangeably when it comes to getting and setting with []
  (__getitem__ and __setitem__). (They should probably both inherit from
  some tableau base class to make this explicity, but it doesn't really
  matter with Python 'duck typing' of which this is a use (abuse?).
  """
  
  import os,sys
  from math import pi
  import numpy.oldnumeric as Numeric
  from Bio.PDB import *
  
  from ptnode import *
  from ptdomain import PTDomain
  from ptsecstruct import pdb_chainid_to_stride_chainid
  from ptutils import cleanup_tmpdir
  
  #-----------------------------------------------------------------------------
  #
  # Module globals
  #
  #-----------------------------------------------------------------------------
  
  
  # constants
  
  TABLEAU_MIN_HELIX_LEN  = 4 # min length of a helix in TableauCreator is 4
  TABLEAU_MIN_STRAND_LEN = 2 # min length of a strand in TableauCreator is 2
  
  # global variables
  
  verbose = False
  
  #-----------------------------------------------------------------------------
  #
  # Class definitions 
  #
  #-----------------------------------------------------------------------------
  
  class PTTableau:
      """
      The PTTableau class is a protein tableau, as per Kamat and Lesk 2007
      'Contact Patterns Between Helices and Strands of Sheet Define Protein
      Folding Patterns' Proteins 66:869-876.
  
      The tableau is a 2 dimensional symmetric matrix indexed by SSEs
      in the protein where each entry
      is a two character code representing the angle between those SSEs.
      (See paper(s) for details).
  
      We implement it as a mapping container, i.e. using __getitem__ and
      __setitem__ so that elements can bet get/set with dictionary/array
      type syntax e.g. tableau[(helix1,strand2)]. (index is a tuple of SSEs
      represented by PTNode objects)- NB must have in parens to ensure tuple.
      This is implemented with a standard dictionary object,
      and since it is symmetric, only one copy is stored, the one where
      i < j in (i,j) index; however either can be get/set, they are swapped
      internally if necessary. Accessing (i,i) returns
      'xa', 'xi', 'xg' for respectively alpha,pi,310 helices and
      'e ' for strand.
  
  
      """
  
      def __init__(self, nodelist):
          """
          Intialize a PTTableau with no tableau entries set yet.
  
          Parameters:
             nodelist - list of nodes that will be in the tableau, ins
                        residue sequence number order.
          """
          self.tabdict = {} # { (res1, res2) : code }; see class documentation\
          self.nodelist = nodelist # ptnodes with those not in tableau removed
          
  
      def __str__(self):
          """
          Return string representation of the tableau; we will write a full matrix
          just like TableauCreator does.
          """
          s = ""
          for sse1 in self.nodelist:
              for sse2 in self.nodelist:
                  try:
                      s += self[(sse1, sse2)] + ' '
                  except KeyError:
                      s += "?? "
              s += '\n'
          return s
      
          
  
      #
      # methods defined to implement container type
      #
  
      def __len__(self):
          """
          Return number of SSEs in the tableau
  
          Parameters: None
          Return value: Number of SSEs in nodelist for building tableau
          """
          return len(self.nodelist)
  
  
      def __getitem__(self, ssepair):
          """
          Return the entry in the tableau for the pair of SSEs
          (sse1, sse2) where sse1 and sse2 are PTNode objects;
          or if ssepair is (i,j) where i,j are integers, the corresponding
          tableau entry for (nodelist[i], nodelist[j]).
  
          Parameters:
             ssepair - tuple (sse1,sse2) (PTNode objects) or
                       tuple (i,j) (integers) to look up tableau entry for
          Return value:
             two character tableau string e.g. 'RD' or 'HH', or '  ' (2 spaces).
             On the main diagonal (self-orientation) since this has no menaing
             we return a (two-char) encoding of the SSE type instead:
             'xa', 'xi', 'xg' for respectively alpha,pi,310 helices and
             'e ' for strand.
  
          Raises Exceptions:
              TypeError if ssepair is not PTNode pair or int pair.
          """
          ssespec1 = ssepair[0]
          ssespec2 = ssepair[1]
          if isinstance(ssespec1, PTNode) and isinstance(ssespec2, PTNode):
              sse1 = ssespec1
              sse2 = ssespec2
          elif isinstance(ssespec1, int) and isinstance(ssespec2, int):
              sse1 = self.nodelist[ssespec1]
              sse2 = self.nodelist[ssespec2]
          else:
              raise TypeError("bad tuple type in PTTableau getitem")
          if sse1 == sse2:
              if isinstance(sse1, PTNodeHelix):
                  if sse1.get_type() == "ALPHA":
                      return "xa"
                  elif sse1.get_type() == "PI":
                      return "xi"
                  elif sse1.get_type() == "310":
                      return "xg"
                  else:
                      return "??" # should not happen
              elif isinstance(sse1, PTNodeStrand):
                  return "e "
              else:
                  return "??" # should not happen
          elif sse1 < sse2:
              ssepair = (sse1,sse2)
          else:
              ssepair = (sse2,sse1)
          return self.tabdict[ssepair]
  
  
      def __setitem__(self, ssepair, tabcode):
          """
          Set the entry in the tableau for the pair of SSEs (sse1,sse2)
          specified as the key (ssepair) parameter to the tabcode value.
  
          Parameters:
             ssepair - tuple (sse1,sse2) to set.
             tabccode - two character tableau string e.g. 'RD' or 'HH', or '  '.
  
          Return value: None
  
          Raises exceptions:
             TypeError  if tabcode is not a valid 2 char uppercase string or '  '
          """
          if len(tabcode) != 2 or not tabcode.isupper() and not tabcode.isspace():
              raise TypeError("bad tableau code '" + tabcode + "'\n")
  
          if (tabcode[0] not in ['L','R','P','O'] or   \
              tabcode[1] not in ['E','D','S','T']) and \
             tabcode != 'HH' and tabcode != 'KK':
              raise TypeError("bad tableau code '" + tabcode + "'\n")
  
          sse1 = ssepair[0]
          sse2 = ssepair[1]
          if sse1 == sse2:
              return
          elif sse1 < sse2:
              ssepair = (sse1, sse2)
          else:
              ssepair = (sse2, sse1)
          self.tabdict[ssepair] = tabcode
          
  
      # have not implemented: __delitem__, __iter__, __contains__
  
  
      # TODO: work out how to implement things like tab[2:] to get row 2,
      # ust like Numeric.array etc.
      def getrow(self, i):
          """
          Return a row of the tableau as a list of tableau codes.
  
          Parameters:
             i - row to get 0 <= i < len(self)
  
          Return value:
             list of two-character tableau codes for row i.
          """
          return [self[(i,j)] for j in xrange(len(self))]
              
  
  class PTTableauPacked:
      """
      The PTTableauPacked class is a compact representation
      of a  protein tableau, as per Kamat and Lesk 2007
      'Contact Patterns Between Helices and Strands of Sheet Define Protein
      Folding Patterns' Proteins 66:869-876.
  
      The tableau is a 2 dimensional symmetric matrix indexed by SSEs
      in the protein where each entry
      is a two character code representing the angle between those SSEs.
      (See paper(s) for details).
  
      We implement it as a mapping container, i.e. using __getitem__ and
      __setitem__ so that elements can bet get/set with dictionary/array
      type syntax e.g. tableau[(1,2)]. (index is a pair of sequential
      SSE numbers, from 0 to n-1 where n is the order of tableau ie number
      of SSEs) - NB must have in parens to ensure tuple.
  
      This is the compact respresentation, storing tableau as simply
      a linear string of two-character tableau codes, in the same as as the
      LAPACK 'packed' format for triangular/symmetric arrays. i.e.
      each column of the matrix is stored in sequence.
      We could save even more space by using only 4 bits for each tableau
      code (since there are only 16 possible codes), but in Python
      it doesn't really make sense to try to be so efficient - but
      we are trying to save space to some degree so that the entire ASTRAL
      PDB non-redundant set or similar can be loaded as tableaux in
      memory.
      As it happens, strings in python don't even support item assignemnt,
      so we have to store it as a list anyway  i.e ['xa','OT',...]
      instead of 'xaOT...'
  
      Unlike PTTableau, this format contains just the tableau codes
      and diagonal SSE type entries, i.e. just character data. there
      are no PTNode object references or anything, so it is simple and
      quick to dump/load with Python pickle module (or similar) with
      no need to build all sorts of other objects (PTNode, Bio.PDB.Structure,
      etc.).
  
      Accessing (i,i) returns 'xa', 'xi', 'xg' for respectively
      alpha,pi,310 helices and 'e ' for strand.
      """
  
      def __init__(self, tableau):
          """
          Intialize a PTTableauPacked given an already built tableaux in
          the full PTTableau format.
  
          Parameters:
              tableau - an already built PTTableau object
          """
          #self.n = 3 
          #self.uplist = ['xa','OS','e ','OT','PE','xa']
          self.n = len(tableau)  # order of tableau (number of SSEs)
          self.uplist = []       # packed format of matrix upper triangle
                                 # NB COLUMN-MAJOR (LAPACK style)
          for j in range(self.n):
              for i in range(j+1):
                  try:
                      tabcode = tableau[(i,j)]
                  except:
                      tabcode = '??'
                  self.uplist.append(tabcode)
          assert(len(self.uplist) ==  self.n * (self.n + 1) / 2)
  
      def __str__(self):
          """
          Return string representation of the tableau; we will write a full matrix
          just like TableauCreator does.
          """
          s = ""
          for i in range(self.n):
              for j in range(self.n):
                  s += self[(i,j)] + ' '
              s += '\n'
          return s
          
  
      #
      # methods defined to implement container type
      #
  
      def __len__(self):
          """
          Return number of SSEs respresented the tableau
  
          Parameters: None
          Return value: order of tableau
          """
          return self.n
      
  
      def __getitem__(self, ssepair):
          """
          Return the entry in the tableau for the pair of SSEs
          (i,j) where i,j are integers, 0 <= i,j < n.
  
          Parameters:
             ssepair - tuple (i,j) (integers) to look up tableau entry for
          Return value:
             two character tableau string e.g. 'RD' or 'HH', or '  ' (2 spaces).
             On the main diagonal (self-orientation) since this has no menaing
             we return a (two-char) encoding of the SSE type instead:
             'xa', 'xi', 'xg' for respectively alpha,pi,310 helices and
             'e ' for strand.
  
          Raises Exceptions:
              TypeError if ssepair is not PTNode pair or int pair.
          """
          i = ssepair[0] + 1
          j = ssepair[1] + 1 # more convenient to have 1 < i,j <= n internally
          if j < i:
              tmp = i
              i = j
              j = tmp
          r = i + j*(j-1)/2
          r -= 1             # back to zero-based for list indexing
          return self.uplist[r]
  
  
      def __setitem__(self, ssepair, tabcode):
          """
          Set the entry in the tableau for the pair of SSEs (sse1,sse2)
          specified as the key (ssepair) parameter to the tabcode value.
  
          Parameters:
             ssepair - tuple (i,j) to set, 0 <= i,j < n.
             tabccode - two character tableau string e.g. 'RD' or 'HH',
                        or 'xa','e ' etc. if i==j for SSE type on diagonal.
  
          Return value: None
  
          Raises exceptions:
             TypeError  if tabcode is not a valid 2 char uppercase string
                         or lowercase type code for diagonal.
          """
          i = ssepair[0]
          j = ssepair[1]
          if i == j:
              if tabcode not in ['xa','xi','xg','e ']:
                  raise TypeError('bad tableau sse type code ' + tabcode + '\n')
          else:
              if len(tabcode) != 2 or not tabcode.isupper():
                  raise TypeError("bad tableau code '" + tabcode + "'\n")
              
              if ( (tabcode[0] not in ['L','R','P','O'] or
                    tabcode[1] not in ['E','D','S','T']) and
                   tabcode != 'HH' and tabcode != 'KK' ):
                  raise TypeError("bad tableau code '" + tabcode + "'\n")
  
              if j < i:
                  tmp = i
                  i = j
                  j = tmp
          i += 1
          j += 1  # more convenient to have 1 < i,j <= n internally
          r = i + j*(j-1)/2  # location in packed rep assuming each entry len 1
          r -= 1             # back to zero-based for list indexing
          self.uplist[r] = tabcode
          
      # have not implemented: __delitem__, __iter__, __contains__
      
      # TODO: work out how to implement things like tab[2:] to get row 2,
      # ust like Numeric.array etc.
      def getrow(self, i):
          """
          Return a row of the tableau as a list of tableau codes.
  
          Parameters:
             i - row to get 0 <= i < len(self)
  
          Return value:
             list of two-character tableau codes for row i.
          """
          # TODO: we could do this more efficiently for packed tableau
          return [self[(i,j)] for j in xrange(len(self))]
  
  
  
  #-----------------------------------------------------------------------------
  #
  # Function definitions
  #
  #-----------------------------------------------------------------------------
  
  def angle_to_tabcode(omega):
      """
      Convert an angle (radians in (-pi, pi]) to a two-character tableau
      code (double quadrant encoding) as described in the papers cited
      at top of module.
  
      Parmaters:
         omega - relative angle to encode
      Return value:
         two-character tableua code (OS, PD, etc.)
      Raises exceptions:
         ValueError if angle is out of range 
      """
      if omega > -pi/4 and omega <= pi/4:
          tabcode = 'P'                           # parallel
      elif omega > pi/4 and omega <= 3*pi/4:
          tabcode = 'R'                           # crossing-right
      elif (omega > 3*pi/4 and omega <= pi) or omega > -pi and omega <= -3*pi/4:
          tabcode = 'O'                           # antiparallel (opposite)
      elif omega > -3*pi/4 and omega <= -pi/4:
          tabcode = 'L'                           # crossing-left
      else:
          raise ValueError('bad omega value ' + str(omega) + '\n')
  
      if omega > 0 and omega <= pi/2:
          tabcode += 'D'                          # dinner
      elif omega > pi/2 and omega <= pi:
          tabcode += 'T'                          # tea
      elif omega > -pi and omega <= -pi/2:
          tabcode += 'S'                          # supper
      elif omega > -pi/2 and omega <= 0:
          tabcode += 'E'                          # elevenses
      else:
          raise ValueError('bad omega value ' + str(omega) + '\n')
  
      return tabcode
  
  
  
  def compute_tableau(ptnode_list, pdb_structure, use_hk=True):
      """
      Build a PTTableau object for the tableau by computing relative angles
      between all SSEs in the ptnode_list.
  
      Parameters:
          ptnode_list - list of PTNode objects (ie iterable of PTNode)
                           representing the SSEs (helices,strands) the
                           tabelau is for.
          pdb_structure - parsed Bio.PDB structure
          use_hk - If True, use the HH and KK codes for respectively
                   antiparallel and parallel strands of the same sheet.
                   Default True.
  
      Return value:
         PTTableau object with entry for each pair of SSEs.
      """
      tableau = PTTableau(ptnode_list)
      for i in range(len(ptnode_list)):
          for j in range(i+1, len(ptnode_list)):
              omega = ptnode_list[i].relative_angle(ptnode_list[j], pdb_structure)
              if omega != None:
                  try:
                      tabcode = angle_to_tabcode(omega)
                  except ValueError:
                      sys.stderr.write('WARNING: catch bad tableau angle, seting Parallel (%d,%d)\n' %(i,j))
                      tabcode = "PE" # NaN -> 0.0 -> parallel: should not happen but does e.g. d7pcka_ -35
                  # set tabcode to HH for antiparallel strands and
                  # KK for parallel strands
                  if (use_hk and
                      isinstance(ptnode_list[i], PTNodeStrand) and
                      isinstance(ptnode_list[j], PTNodeStrand) and
                      (ptnode_list[i].get_sheet_id() != None and
                       ptnode_list[i].get_sheet_id() ==
                       ptnode_list[j].get_sheet_id())):
                      if tabcode[0] == 'O':
                          tableau[(ptnode_list[i], ptnode_list[j])] = 'HH'
                      elif tabcode[0] == 'P':
                          tableau[(ptnode_list[i], ptnode_list[j])] = 'KK'
                      else:
                          tableau[(ptnode_list[i], ptnode_list[j])] = tabcode
                  else:
                      tableau[(ptnode_list[i], ptnode_list[j])] = tabcode
  
      if verbose:
          sys.stderr.write(str(tableau))
  
      return tableau
  
  
  def compute_omega_matrix(ptnode_list, pdb_structure):
      """
      Return the omega (relative angle, in radians) matrix as a 2D Numeric.array
      by computing relative angles between all SSEs in the ptnode_list
  
      Parameters:
          ptnode_list - list of PTNode objects (ie iterable of PTNode)
                           representing the SSEs (helices,strands) the
                           tabelau is for.
          pdb_structure - parsed Bio.PDB structure
  
      Return value:
         Numeric.array square symmetric (order length of ptnode_list) where
         each entry is relative angle between SSEs in radians.
         Main diagonal entries set to 0.
      
      """
      n = len(ptnode_list)
      omega_array = Numeric.zeros((n, n), Numeric.Float)
      for i in range(n):
          for j in range(i+1, n):
              omega = ptnode_list[i].relative_angle(ptnode_list[j], pdb_structure)
              if omega == None:
                  omega_array[i, j] = float('NaN')
              else:
                  omega_array[i, j] = omega
              omega_array[j, i] = omega_array[i, j]
  
      # set the diagonal as follows:
      # 0.00 for strand
      # 1.00 for alpha helix
      # 2.00 for pi helix
      # 3.00 for 3_10 helix
      for i in range(n):
          if isinstance(ptnode_list[i], PTNodeHelix):
              if ptnode_list[i].get_type() == "ALPHA":
                  v = 1.00
              elif ptnode_list[i].get_type() == "PI":
                  v = 2.00
              elif ptnode_list[i].get_type() == "310":
                  v = 3.00
              else:
                  pass # should not happen
          elif isinstance(ptnode_list[i], PTNodeStrand):
              v = 0.00
          omega_array[i,i] = v
  
      return omega_array
  
      
  #-----------------------------------------------------------------------------
  #
  # Classes and functions for running external TableauCreator
  #
  #-----------------------------------------------------------------------------
  
  
  # Inherit from the PDBIO.Select class for writing only parts of PDB to file
  # See the Bio.PDB documentation: biopython-1.43/Doc/biopdb_faq.pdf
  class DomainSelect(Select):
      """
      The DomainSelect class inherits from the PDBIO.Select class
      and overrides function to select only certain residues for writing
      ATOM records in the domain we are interested in to the
      simplified PDB file for TableauCreator.
  
      See the Bio.PDB documentation by Thomas Hamelryck:
        biopython-1.43/Doc/biopdb_faq.pdf
      """
      def __init__(self, domain):
          """
          Constructor for the DomainSelect class, sets the domain member
          used to accept only residues in that domain.
          Parameters:
             domain - ptdomain object of domain to select residues from
          """
          self.domain = domain
  
      def __repr__(self):
          """
          Overrides the base __repr__ to write out the domain we have
          """
          return "<DomainSelect: " + str(self.domain) + ">"
      
      def accept_residue(self, residue):
          """
          overrides the base accept_residue() function to accept only
          residues in our domain. Also reject HETATMS.
          Paramteters:
             residue - Bio.PDB Residue object of residue to test
          Return value:
             1 to accept residue, 0 to reject.
          """
          chain = residue.get_parent()
          chainid = pdb_chainid_to_stride_chainid(chain.get_id())
          # id of a residue in Bio.PDB is tuple (hetatm, resseqnum, icode)
          # so we choose those where chain and residue PDB number
          # is in the domain.
          resnum = residue.get_id()[1]
          if (self.domain.is_in_domain(chainid, resnum) and
              residue.get_id()[0] == ' '):
              return 1
          else:
              return 0
  
  
  
  def get_tableau_from_pdbstruct(pdbid, domain,
                                 pdb_structure, ptnode_list):
      """
      Build a PTTableau object for the tableau by first creating a
      simple PDB file with only the ATOM records for residues in the
      domain we are processing, and also a .SSEsInfo file containing the
      secnodary structure assignments we already have, then running
      TableauCreator on it (using our simple PDB file and SSEsInfo) and
      parsing the output.
  
      Parameters:
          pdbid - PDB identifier of the strucutre
          domain - The PTDomain object for our current domain
          pdb_structure - parsed Bio.PDB structure
          ptnode_list - list of PTNode objects (ie iterable of PTNode)
                           representing the SSEs (helices,strands) the
                           tabelau is for.
      Return value:
         PTTableau object built from TableauCreator output
  
      """
      TMPDIR = os.tempnam(None, "pttabin")
      os.mkdir(TMPDIR)
      try:
          filename = pdbid
          if domain.domainid != None:
              filename += '-' + domain.domainid
          filename += '.pdb'
          domain_pdb_filename = os.path.join(TMPDIR, filename)
          io = PDBIO()
          io.set_structure(pdb_structure)
          io.save(domain_pdb_filename, DomainSelect(domain))
  
          ssesinfo_filename = os.path.join(TMPDIR, filename + ".input-SSEsInfo")
          write_ssesinfo(ssesinfo_filename, ptnode_list)
  
          tableau =  read_tableau_from_tableaucreator(domain_pdb_filename,
                                                      ptnode_list,
                                                      ssesinfo_filename)
          os.unlink(domain_pdb_filename)
          os.unlink(ssesinfo_filename)
      finally:
          cleanup_tmpdir(TMPDIR)
      return tableau
  
      
  def read_tableau_from_tableaucreator(pdb_filename, ptnode_list,
                                       ssesinfo_filename):
      """
      Run Arun's TableauCreator program on the supplied pdb_filename
      using SSEsInfo file.
  
      Parameters:
         pdb_filename - PDB file to run TableauCreator on
         ptnode_list - list of PTNode objects (ie iterable of PTNode)
                           representing the SSEs (helices,strands) the
                           tabelau is for.
         ssesinfo_filename - filename of the .SSEsInfo file that was written
                             to define SSEs for TableauCreator.
      Return value:
         PTTableau object built from TableauCreator output
         
      NB: TableauCreator is not yet published or available (October 2007)
      and I am using a private version which Arun sent me, which I modified
      to add the -s option to use STRIDE rather than DSSP
      and to have the -i option to parse .SSEsInfo files.
      """
  
      # TableauCreator needs an output directory where it writes all its
      # intermediate/output files, only puts progress information/errors
      # to stdout/stderr.
  
      tmpdir = os.tempnam(None, "pttab")
      os.mkdir(tmpdir)
      command = "TableauCreator "
      command += "-i " + ssesinfo_filename + " "
      command += pdb_filename + " " + tmpdir
      command += " >/dev/null"
      if verbose:
          sys.stderr.write("running '" + command + "'...")
      os.system(command)
      if verbose:
          sys.stderr.write("done\n")
      # output files are:
      #   <pdbfilename>.angles
      #   <pdbfilename>.SSEsInfo
      #   <pdbfilename>.stride or <pdbfilename>.dssp
      #   <pdbfilename>.tableau
      outfile_prefix = os.path.join(tmpdir, os.path.basename(pdb_filename))
      if not os.path.isfile(os.path.join(tmpdir, "TABCREATE_OK")):
          sys.stderr.write("ERROR: TableauCreator failed\n")
          cleanup_tmpdir(tmpdir)
          return None
      # Now the tricky thing is TableauCreator indexes its matrix just with
      # purely sequential numbers from 0 (as conventional)
      # assuming all SSEs in one domain and in fact one chain
      # (so we handle this by creating our own simple PDB file with only
      # ATOM records for our current domain, and only one TER record on
      # end so chains concatenated effectively).
      # And also (as in comments above functions) we have the dodginess of
      # doing the same thing in different ways in multiple places (DSSP/STRIDE
      # parsing, PDB parsing, etc.).
      # So let's check that the TableauCreator SSE info lines up with ours
      # (otherwise we can't use the tableau data).
  
      # parse the SSEsInfo file and check lines up with ptnodes,
      # returns list of ptnodes corresponding to Tableau entries (may be shorter
      # than our input node list; some removed as no equivalent in tableua).
      nodelist = parse_tableaucreator_ssesinfo(outfile_prefix + '.SSEsInfo',
                                               ptnode_list)
      if nodelist != None:
          tableau_filename = outfile_prefix + ".tableau"
          tableau = parse_tableaucreator_output(tableau_filename, nodelist)
          if tableau != None:
              if verbose:
                  sys.stderr.write(str(tableau))
          else:
              sys.stderr.write('WARNING: problem parsing TableauCreator output;\n'
                           '         tableau information will not be used\n')
      else:
          sys.stderr.write('WARNING: problem with TableauCreator output;\n'
                           '         tableau information will not be used\n')
          tableau = None
  
      cleanup_tmpdir(tmpdir)
      return tableau
   
  
  def parse_tableaucreator_ssesinfo(filename, nodelist):
      """
      Parse the .SSEsInfo file created by TableauCreator and check that
      it lines up with our SSE info in the form of the list of helix/strand
      PTNodes.
  
      Parameters:
        filename - filename of the .SSEsInfo file
        nodelist - list of PTNodes, in order of residue sequence number
  
      Return value:
        nodelist where nodes with no tableau entry removed (too short) if
        OK (they line up) else None (different number of nodes,
        residue sequence numbers/types don't match, etc).
      """
  
      # first remove all nodes with len < TABLEAU_MIN_SSE_LEN, since Tableau
      # Creator won't have entries for them; we will have to set them
  #     ptnodelist = [ node for node in nodelist \
  #                    if ( (isinstance(node, PTNodeHelix) and
  #                          node.get_span() >= TABLEAU_MIN_HELIX_LEN) or
  #                         (isinstance(node, PTNodeStrand) and
  #                          node.get_span() >= TABLEAU_MIN_STRAND_LEN) )
  #                  ]
  
      ptnodelist = nodelist 
      # FIXME: no longer need this filtering of too short SSEs,
      # now that .SSEsInfo input is being used
      if len(ptnodelist) != len(nodelist):
          sys.stderr.write('WARNING: no tableau entry for '
                           + str(len(nodelist)-len(ptnodelist)) +
                           ' nodes due to length to small' 
                           + '\n')
      
      fh = open(filename)
      # first line is number of SSEs, subsequent lines are 
      # type start_resnum end_resnum
      # where type is E or H (DSSP code) and resnums are PDB residue numbers
      # 
      numlines = int(fh.readline())
      if numlines < 2:
          sys.stderr.write('ERROR: bad SSEsInfo data\n')
          fh.close()
          return None
      linenum = 1
      sseinfo = [] # tuple (type, start, end)
      line = fh.readline()
      while line != "":
          fields = line.split()
          sseinfo.append((fields[0], int(fields[1]), int(fields[2])))
          linenum += 1
          line = fh.readline()
      fh.close()
      if len(sseinfo) != numlines:
          sys.stderr.write('ERROR: TableauCreator SSEsInfo file specified ' \
                           + str(numlines) + ' entries but ' \
                           + str(len(sseinfo)) + ' read\n')
          return None
      if len(sseinfo) != len(ptnodelist):
          sys.stderr.write('ERROR: TableauCreator SSEsInfo file has ' \
                           + str(numlines) + ' entries but ' \
                           + 'we have ' + str(len(ptnodelist)) + ' SSEs\n')
          return None
      for i in range(len(sseinfo)):
          sse = sseinfo[i]
          ptnode = ptnodelist[i]
          if sse[0] == 'H' and not isinstance(ptnode, PTNodeHelix) or \
             sse[0] == 'E' and not isinstance(ptnode, PTNodeStrand) or \
             sse[1] != ptnode.get_start_res_seq() or \
             sse[2] != ptnode.get_end_res_seq():
              sys.stderr.write('ERROR: TableauCreator SSEInfo entry ' + \
                               str(sse) + ' does not match node ' +
                               str(ptnode) + '\n')
              return None
      return ptnodelist
  
  
  def parse_tableaucreator_output(filename, nodelist):
      """
      Parse the .tableau file created by TableauCreator
  
      Parameters:
         filename - filename of the .tableau file
         nodelist - list of PTNodes, in order of residue sequence number
  
      Return value:
        PTTableau for the tableau parsed
      """
      tableau = PTTableau(nodelist)
      # First line of file is number of SSEs (order of square matrix)
      # The whole matrix is in the file (it is symmetric), diagonal elements
      # and other unset (ie non-contact) elements are set to '--', or '  '
      # (two spaces) when cannot be calculated.
      fh = open(filename)
      numlines = int(fh.readline())
      if numlines < 2:
          sys.stderr.write('ERROR: bad tableau data\n')
          fh.close()
          return None
      linenum = 1
      line = fh.readline()
      i = 0
      while line != "":
          # we will use the fact that fields are fixed length rather than
          # splitting on space separator as fields may be set to '  ' (2 spaces)
          # when error calculating (helices too short etc.) in TableauCreator.
          # fields are two chars with two spaces between each
          node_i = nodelist[i]
          if len(line) < len(nodelist) * 4:
              sys.stderr.write('ERROR: bad line in tableau; line too short:\n')
              sys.stderr.write(line)
              fh.close()
              return None
          for j in range(i+1, len(nodelist)): # no need to store both i,j and j,i
              node_j = nodelist[j]
              col = 4*j  # each field as 2 char tabcode then 2 spaces
              tabcode = line[col:col+2]
              if tabcode !=  '  ' and tabcode != '--':
                  tableau[(node_i, node_j)] = line[col:col+2]
          i += 1
          linenum += 1
          line = fh.readline()
      fh.close()
      return tableau
          
  
  
  def write_ssesinfo(filename, nodelist):
      """
      Write a TableauCreator .SSEsInfo file describing the SSE asignments
      we have to the specified filename. This is so that we avoid having
      TableauCreator re-run DSSP or STRIDE for the assignments, which is
      inefficient and leads to inconsistencies. TableauCreator has been
      modified to be able to read this .SSEsInfo file instead, allowing
      the same assignments we have to be re-used by TableauCreator.
  
      WARNING: file is overwritten if it exists
      
      Parameters:
         filenanme - filename to write SSEsInfo to
         nodelist - list of PTNodes defining the SSEs
  
      Return value: None
      """
      # The format of the .SSEsInfo file is (see writeTableauAnglesSSEinfo())
      # that the first line has number of records and each subsequent
      # line (record)
      # is whitespace-separated:
      # dssp-code start end chainid
      # e.g.
      # H   10  21  A
      # Only H and E codes are used.
      # blank chainid is not allowed, '-' used instead.
      fh = open(filename, 'w')
      fh.write(str(len(nodelist)) + "\n")
      for node in nodelist:
          if isinstance(node, PTNodeHelix):
              typecode = 'H'
          elif isinstance(node, PTNodeStrand):
              typecode = 'E'
          else:
              assert(False)
          fh.write(typecode + " " + str(node.get_start_res_seq()) + " " +
                   str(node.get_end_res_seq()) + " " + node.get_chainid() + "\n")
      fh.close()
  
      
  def pttableau_set_verbose(verb):
      """
      set the module global verbose flag in this module to supplied value
      Parameters: verb - True (for verbose output) or False
      Return value: None
      Uses globals: verbose (in this module)
      """
      global verbose
      verbose = verb
      
  