molecule.java

A program to find frequent molecular substructures and discriminative fragments in a database of mol

      a = this.atoms[i];        /* traverse the molecule's atoms */
      if (a.type != t) continue;/* that have the right type */
      avec[0] =  a;             /* note the root atom and */
      a.mark  =  0;             /* mark it with its index */
      list = this.embed(seed, 0, 1, avec, list);
      a.mark  = -1;             /* match seed bonds recursively */
      if (list == found) break; /* if only to check containment, */
    }                           /* check the recursion result */
    seed.atoms[0].mark = -1;    /* unmark the first seed atom */
    return list;                /* return the recursion result */
  }  /* embed() */

  /*--------------------------------------------------------------------
  For the above function to work properly, the seed to embed must either
  have been processed with the function 'sort' with the parameter 'seed'
  set to true, or it must have been created from a fragment that itself
  was generated in the search process, so that the edges are sorted in
  a specific way (so that source nodes are always already numbered).
  --------------------------------------------------------------------*/

  protected Embedding embed    (Molecule seed)
  { return this.embed(seed, null); }

  protected boolean   contains (Molecule seed)
  { return this.embed(seed, Molecule.check) == found; }

  /*------------------------------------------------------------------*/

  protected void map (Extension ext)
  {                             /* --- build map for canonical form */
    int  i, k, n;               /* loop variables */
    Bond b;                     /* to traverse the bonds */

    k = this.bondcnt;           /* number bonds and build bond map */
    for (i = k; --i >= 0; ) ext.bonds[i].mark = i;
    for (i = k; --i >= 0; ) ext.word[i] = this.bonds[i].mark;
    System.arraycopy(ext.bonds, 0, this.bonds, 0, this.bondcnt);
    n = ext.src;                /* copy the new bonds vector */
    for (i = 0; i < k; i++) {   /* traverse the bonds and */
      b = ext.bonds[i];         /* number the atoms accordingly */
      if      (b.src.mark < 0) ext.atoms[b.src.mark = ++n] = b.src;
      else if (b.dst.mark < 0) ext.atoms[b.dst.mark = ++n] = b.dst;
    }                           /* also sort atoms into new order */
    for (k += i = this.atomcnt; --i >= 0; )  /* build atom map */
      ext.word[--k] = this.atoms[i].mark;    /* and copy atoms */
    System.arraycopy(ext.atoms, 0, this.atoms, 0, this.atomcnt);
  }  /* map() */

  /*------------------------------------------------------------------*/

  protected boolean makeCanonic (Extension ext)
  { return this.makeCanonic(ext, -1); }

  protected boolean makeCanonic (Extension ext, int keep)
  {                             /* --- turn into canonical form */
    int     i, n;               /* loop variables, buffers */
    Atom    a, root;            /* to traverse the atoms, new root */
    Bond    b;                  /* to traverse the bonds */
    boolean changed;            /* whether fragment was changed */

    n = (keep >= 0) ? keep : 0; /* get number of bonds to keep */
    if (this.bondcnt <= n)      /* if no bonds are moveable, */
      return false;             /* the molecule's word is minimal */
    this.sort(false);           /* prepare molecule for processing */
    ext.init(this, 0);          /* and init. the extension object */
    root = this.atoms[0];       /* note the old root node */
    if (keep < 0) {             /* if full freedom of reordering */
      for (i = 0; i < this.atomcnt; i++) {
        a = this.atoms[i];      /* traverse the molecule's atoms */
        if (a.type < root.type) root = a;
      }                         /* find atom with smallest type */
      changed = (root != this.atoms[0]);
      if (changed) {            /* if the root has to be replaced, */
        ext.word[0] = root.type;/* set only the initial character */
        ext.word[1] = Integer.MAX_VALUE; }
      else {                    /* if the root atom may be kept */
        System.arraycopy(this.bonds, 0, ext.bonds, 0, this.bondcnt);
        ext.makeWord(this, this.bondcnt);
      }                         /* construct an initial code word */
      this.unmark();            /* unmark all atoms and bonds */
      for (i = this.atomcnt; --i >= 0; ) {
        a = this.atoms[i];      /* traverse the molecule's atoms */
        if (a.type != root.type)/* check one letter prefix words */
          continue;             /* (must coincide with root type) */
        a.mark = 0;             /* mark the potential root and */
        ext.atoms[0] = a;       /* store it as the first atom */
        if (ext.makeCanonic(0, 0, 1)) { root = a; changed = true; }
        a.mark = -1;            /* construct code word recursively, */
      }                         /* then unmark the pot. root again */
      ext.atoms[0] = root;      /* set the (new) root atom and */
      root.mark = n = 0; }      /* mark it as the first atom */
    else {                      /* if to keep some bonds */
      System.arraycopy(this.bonds, 0, ext.bonds, 0, this.bondcnt);
      ext.makeWord(this, this.bondcnt);
      n = 0;                    /* construct an initial code word */
      for (i = keep; --i >= 0; ) {  /* traverse the bonds to keep */
        b = this.bonds[i]; b.mark = 0;
        if (b.src.mark > n) n = b.src.mark;
        if (b.dst.mark > n) n = b.dst.mark;
      }                         /* find atom with highest index */
      ext.src = n;              /* and note this atom's index */
      for (i = n+1; --i >= 0; ) /* copy the marked/visited atoms */
        ext.atoms[i] = this.atoms[i];
      for (i = this.atomcnt; --i > n; )
        this.atoms[i].mark = -1;/* unmark the unvisited atoms */
      for (i = this.bondcnt; --i >= keep; )
        this.bonds[i].mark = -1;/* unmark the bonds to reorder */
      changed = ext.makeCanonic(keep, 0, n+1);
    }                           /* construct code word recursively */
    if (changed)                /* if the molecule was changed, */
      this.map(ext);            /* copy changes and build maps */
    return changed;             /* return whether molecule is changed */
  }  /* makeCanonic() */

  /*------------------------------------------------------------------*/

  protected int isCanonic (Extension ext)
  { return this.isCanonic(ext, this.bondcnt); }

  protected int isCanonic (Extension ext, int fixed)
  {                             /* --- check for canonical form */
    int  i, k, r, t;            /* loop variable, buffers */
    Atom a;                     /* to traverse the atoms */

    if (this.bondcnt <= 0)      /* if there are no bonds, */
      return 1;                 /* the molecule's word is minimal */
    ext.init(this, fixed);      /* construct the initial code word */
    ext.makeWord(this, this.bondcnt);
    this.sort(false);           /* prepare molecule for the test */
    this.unmark();              /* (sort bonds and clear markers) */
    r = this.bondcnt;           /* set the default result */
    k = this.atoms[0].type;     /* get the root atom's type */
    for (i = this.atomcnt; --i >= 0; ) {
      a = this.atoms[i];        /* traverse the atoms of the fragment */
      if (a.type > k) continue; /* check one letter prefix words */
      if (a.type < k) break;    /* (words are e.g. A (I_s B A I_d)* ) */
      ext.atoms[a.mark = 0] = a;/* mark and note the root atom, */
      t = ext.isCanonic(0,0,1); /* check prefix words recursively, */
      a.mark = -1;              /* then unmark the root atom again */
      if (t < r) { if (t < fixed) return -1; r = t; }
    }                           /* evaluate the recursion result */
    return (r < this.bondcnt) ? 0 : 1;
  }  /* isCanonic() */          /* return the overall result */

  /*------------------------------------------------------------------*/

  protected void mark ()
  {                             /* --- set all markers */
    for (int i = this.atomcnt; --i >= 0; )
      this.atoms[i].mark = 0;   /* mark the atoms of the molecule */
    for (int i = this.bondcnt; --i >= 0; )
      this.bonds[i].mark = 0;   /* mark the bonds of the molecule */
  }  /* mark() */

  /*------------------------------------------------------------------*/

  protected void markId ()
  {                             /* --- set all markers */
    for (int i = this.atomcnt; --i >= 0; )
      this.atoms[i].mark = i;   /* number the atoms of the molecule */
    for (int i = this.bondcnt; --i >= 0; )
      this.bonds[i].mark = i;   /* number the bonds of the molecule */
  }  /* markId() */

  /*------------------------------------------------------------------*/

  protected void unmark ()
  {                             /* --- clear all markers */
    for (int i = this.atomcnt; --i >= 0; )
      this.atoms[i].mark = -1;  /* unmark the atoms of the molecule */
    for (int i = this.bondcnt; --i >= 0; )
      this.bonds[i].mark = -1;  /* unmark the bonds of the molecule */
  }  /* unmark() */

  /*------------------------------------------------------------------*/

  public String toString ()     /* --- produce a string description */
  { return (new SMILES()).format(this); }

  /*------------------------------------------------------------------*/

  public String toString (Extension ext)
  { return ext.makeString(this, true); }

  /*------------------------------------------------------------------*/

  public void toLogic (PrintStream out)
  {                             /* --- produce a Prolog description */
    int  i, elem;               /* loop variable, element buffer */
    Atom a;                     /* to traverse the atoms */
    Bond b;                     /* to traverse the bonds */

    for (i = 0; i < this.atomcnt; i++) {
      a = this.atoms[i];        /* traverse the atoms */
      a.mark = i;               /* and mark them */
      elem = a.getElement();    /* get the chemical element */
      if (this.map != null)     /* and decode it if necessary */
        elem = this.map.decode(elem);
      out.println("atom(\"" +this.id +"\", " +i +", \""
                  +Atom.names[elem] +"\").");
    }                           /* print the atom predicates */
    for (i = 0; i < this.bondcnt; i++) {
      b = this.bonds[i];        /* traverse the bonds */
      out.println("bond(\"" +this.id +"\", "
                  +b.src.mark +", " +b.dst.mark +", \""
                  +Bond.names[b.type] +"\").");
    }                           /* print the bond predicates */
  }  /* toLogic() */

  /*------------------------------------------------------------------*/

  public static final void main (String args[])
  {                             /* --- main function for testing */
    Notation  ntn;              /* molecule notation */
    Molecule  mol;              /* created molecule */
    Extension ext;              /* extension for code words */
    int[]     masks;            /* masks for atom and bond types */

    if (args.length != 1) {     /* if wrong number of arguments */
      mol = new Molecule();     /* create a dummy for the class name */
      System.err.print("usage: java " +mol.getClass().getName()
                      +" <string>");
      return;                   /* print a usage message */
    }                           /* and abort the program */

    masks = new int[4];         /* create atom and bond masks */
    masks[0] = masks[2] = Atom.TYPEMASK;
    masks[1] = masks[3] = Bond.TYPEMASK;

    try {                       /* try SMILES */
      System.out.print("SMILES: ");
      ntn = new SMILES();       /* create a SMILES object */
      mol = ntn.parse(args[0]); /* and parse the argument */
      mol.maskTypes(masks);
      mol.makeCanonic(ext = new MaxSrcExt(0, 0, null));
      System.out.println(ntn.format(mol));
      System.out.println("breadth-first code word:");
      System.out.println(ext.makeString(mol));
      mol.makeCanonic(ext = new RightExt(0, 0, null));
      System.out.println("depth-first   code word:");
      System.out.println(ext.makeString(mol)); return; }
    catch (java.io.IOException ioe) {
      System.out.println(ioe.getMessage()); }

    try {                       /* try SYBYL line notation (SLN) */
      System.out.print("SLN   : ");
      ntn = new SLN();          /* create an SLN object */
      mol = ntn.parse(args[0]); /* and parse the argument */
      mol.maskTypes(masks);
      mol.makeCanonic(ext = new MaxSrcExt(0, 0, null));
      System.out.println(ntn.format(mol));
      System.out.println("breadth-first code word:");
      System.out.println(ext.makeString(mol));
      mol.makeCanonic(ext = new RightExt(0, 0, null));
      System.out.println("depth-first   code word:");
      System.out.println(ext.makeString(mol)); }
    catch (java.io.IOException ioe) {
      System.err.println(ioe.getMessage()); }
  }  /* main() */

}  /* class Molecule */