Logo Search packages:      
Sourcecode: adun.app version File versions  Download package

Selection.m

/* Copyright 2003  Alexander V. Diemand

    This file is part of MolTalk.

    MolTalk is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    MolTalk is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with MolTalk; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 */

/* vim: set filetype=objc: */


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>

#include "Selection.oh"
#include "privateSelection.oh"
#include "Matrix53.oh"
#include "Residue.oh"
#include "Chain.oh"


/*  s e t  t o  g e t  memory allocation debugging */
#undef MEMDEBUG


@implementation Selection


-(id)init   //@nodoc
{
      [super init];
      chain = nil;
      selection = RETAIN([NSMutableArray new]);
      return self;
}


-(void)dealloc    //@nodoc
{
      //printf("Selection_dealloc\n");
      [selection removeAllObjects];
      RELEASE(selection);
      if (chain) 
      {
            RELEASE(chain);
      }
      [super dealloc];
}


/*
 *   returns the number of residues in selection
 */
-(unsigned long)count
{
      return [selection count];
}


/*
 *   return string, which describes this selection
 */
-(NSString*)description
{
      NSString *res = [NSString stringWithFormat: @"Selection of '%@' with %d residues",chain, [self count]];
      return res;
}


/*
 *   give access to selection
 */
-(NSMutableArray*)getSelection
{
      return selection;
}


/*
 *   return an enumerator over the selected residues
 */
-(NSEnumerator*)selectedResidues
{
      return [selection objectEnumerator];
}


/*
 *   test a residue in selection
 */
-(BOOL)containsResidue: (Residue*)r
{
      return [selection containsObject: r];
}


/*
 *   include residue in selection
 */
-(id)addResidue: (Residue*)r
{
      [selection addObject: r];
      //NSLog(@"Selection-addResidue: %@",r);
      return self;
}


/*
 *   remove a residue from this selection
 */
-(id)removeResidue: (Residue*)r
{
      [selection removeObjectIdenticalTo: r];
      return self;
}


/*
 *   exclude other selection from this one (same chain!)
 */
-(id)difference:(Selection*)sel2
{
      if (sel2->chain != chain)
      {
            return self;
      }
      int ct = [sel2->selection count];
      int i;
      Residue *res;
      for (i=0; i<ct; i++)
      {
            res = [sel2->selection objectAtIndex:i];
            if ([selection containsObject:res])
            {
                  [self removeResidue:res];
            }
      }
      return self;
}


/*
 *   include other selection into this one (same chain!)
 */
-(id)union:(Selection*)sel2
{
      if (sel2->chain != chain)
      {
            return self;
      }
      int ct = [sel2->selection count];
      int i;
      Residue *res;
      for (i=0; i<ct; i++)
      {
            res = [sel2->selection objectAtIndex:i];
            if (![selection containsObject:res])
            {
                  [self addResidue:res];
            }
      }
      return self;
}


/*
 *   structurally align two selections (must be 2 different chains)
 *   size of selections must match
 */
-(Matrix53*)alignTo: (Selection*)sel2
{
      /* sel2 is the stationary structure, this structure (chain) is being transformed */
      
      if ([sel2 count] != [self count])
      {
            NSLog(@"Selection-alignTo: both selections must be of the same size (%d != %d).",[self count],[sel2 count]);
            return nil;
      }

      CREATE_AUTORELEASE_POOL(pool);
#ifdef MEMDEBUG
      GSDebugAllocationActive(YES);
      NSLog(@"allocated objects on entering method\n%s",GSDebugAllocationList(NO));
#endif

      /* clean up selections: residues must have coordinates for CA, otherwise remove pair */
      Residue *dummy_residue = [Residue new];
      NSMutableArray *res_enum1 = [self getSelection];
      NSMutableArray *res_enum2 = [sel2 getSelection];
      int idx;
      Residue *res1, *res2;
      for (idx=0;idx<[res_enum1 count]; idx++)
      {
            res1 = [res_enum1 objectAtIndex:idx];
            res2 = [res_enum2 objectAtIndex:idx];
            //printf("%@ <-> %@\n",res1,res2);
            if (!([res1 getCA] && [res2 getCA]))
            {
                  NSLog(@"due to missing CA coordinates the pair %@ <-> %@ has been removed from the superposition.",res1, res2);
                  [res_enum1 replaceObjectAtIndex:idx withObject:dummy_residue];
                  [res_enum2 replaceObjectAtIndex:idx withObject:dummy_residue];
            }
      }
      [res_enum1 removeObjectIdenticalTo:dummy_residue];
      [res_enum2 removeObjectIdenticalTo:dummy_residue];
      //printf("selection 1: %@\n", self);
      //printf("selection 2: %@\n", sel2);

      Matrix *m1 = [self matrixWithCACoords];
      //NSLog(@"m1: %@",m1);
      Matrix *m2 = [sel2 matrixWithCACoords];
      //NSLog(@"m2: %@",m2);
      Matrix *com1 = [m1 centerOfMass];
      //NSLog(@"com1: %@",com1);
      Matrix *com2 = [m2 centerOfMass];
      //NSLog(@"com2: %@",com2);
      
      int irow;
      double o11,o12,o13;
      double o21,o22,o23;
      o11 = [com1 atRow: 0 col: 0];
      o12 = [com1 atRow: 0 col: 1];
      o13 = [com1 atRow: 0 col: 2];
      o21 = [com2 atRow: 0 col: 0];
      o22 = [com2 atRow: 0 col: 1];
      o23 = [com2 atRow: 0 col: 2];
      for (irow=0; irow<[m1 rows]; irow++)
      { /* move to origin */
            [m1 atRow: irow col: 0 subtract: o11];
            [m1 atRow: irow col: 1 subtract: o12];
            [m1 atRow: irow col: 2 subtract: o13];
            [m2 atRow: irow col: 0 subtract: o21];
            [m2 atRow: irow col: 1 subtract: o22];
            [m2 atRow: irow col: 2 subtract: o23];
      }
      Matrix *subm = [m2 msubtract: m1];
      //NSLog(@"subm: %@",subm);
      Matrix *subp = [m1 madd: m2];
      //NSLog(@"subp: %@",subp);
      Matrix *xm = [subm matrixOfColumn: 0];
      Matrix *ym = [subm matrixOfColumn: 1];
      Matrix *zm = [subm matrixOfColumn: 2];
      Matrix *xp = [subp matrixOfColumn: 0];
      Matrix *yp = [subp matrixOfColumn: 1];
      Matrix *zp = [subp matrixOfColumn: 2];

      Matrix *xmym = [xm mmultiply: ym];
      Matrix *xmyp = [xm mmultiply: yp];
      Matrix *xpym = [xp mmultiply: ym];
      Matrix *xpyp = [xp mmultiply: yp];
      Matrix *xmzm = [xm mmultiply: zm];
      Matrix *xmzp = [xm mmultiply: zp];
      Matrix *xpzm = [xp mmultiply: zm];
      Matrix *xpzp = [xp mmultiply: zp];
      Matrix *ymzm = [ym mmultiply: zm];
      Matrix *ymzp = [ym mmultiply: zp];
      Matrix *ypzm = [yp mmultiply: zm];
      Matrix *ypzp = [yp mmultiply: zp];
      Matrix *mdiag = [Matrix matrixWithRows: 4 cols: 4];
      
      double sumall;
      sumall = [[ypzm msubtract: ymzp] sum];
      [mdiag atRow: 0 col: 1 value: sumall];
      sumall = [[xmzp msubtract: xpzm] sum];
      [mdiag atRow: 0 col: 2 value: sumall];
      sumall = [[xpym msubtract: xmyp] sum];
      [mdiag atRow: 0 col: 3 value: sumall];
      sumall = [[ypzm msubtract: ymzp] sum];
      [mdiag atRow: 1 col: 0 value: sumall];
      sumall = [[xmym msubtract: xpyp] sum];
      [mdiag atRow: 1 col: 2 value: sumall];
      sumall = [[xmzm msubtract: xpzp] sum];
      [mdiag atRow: 1 col: 3 value: sumall];
      sumall = [[xmzp msubtract: xpzm] sum];
      [mdiag atRow: 2 col: 0 value: sumall];
      sumall = [[xmym msubtract: xpyp] sum];
      [mdiag atRow: 2 col: 1 value: sumall];
      sumall = [[ymzm msubtract: ypzp] sum];
      [mdiag atRow: 2 col: 3 value: sumall];
      sumall = [[xpym msubtract: xmyp] sum];
      [mdiag atRow: 3 col: 0 value: sumall];
      sumall = [[xmzm msubtract: xpzp] sum];
      [mdiag atRow: 3 col: 1 value: sumall];
      sumall = [[ymzm msubtract: ypzp] sum];
      [mdiag atRow: 3 col: 2 value: sumall];

      [xm square];
      [xp square];
      [ym square];
      [yp square];
      [zm square];
      [zp square];
      
      sumall = [[[xm madd: ym] madd: zm] sum];
      [mdiag atRow: 0 col: 0 value: sumall];
      sumall = [[[yp madd: zp] madd: xm] sum];
      [mdiag atRow: 1 col: 1 value: sumall];
      sumall = [[[xp madd: zp] madd: ym] sum];
      [mdiag atRow: 2 col: 2 value: sumall];
      sumall = [[[xp madd: yp] madd: zm] sum];
      [mdiag atRow: 3 col: 3 value: sumall];
      
      Matrix *eigen = [mdiag jacobianDiagonalizeWithMaxError: 0.0001];
      double q1,q2,q3,q4,val;
      double vmax = FLT_MAX;
      for (irow=0; irow<[eigen cols]; irow++)
      { /* find smallest eigenvalue and corresponding eigenvector */
            val = [eigen atRow: 0 col: irow];
            if (val < vmax)
            {
                  q1 = [eigen atRow: 1 col: irow];
                  q2 = [eigen atRow: 2 col: irow];
                  q3 = [eigen atRow: 3 col: irow];
                  q4 = [eigen atRow: 4 col: irow];
                  vmax = val;
            }
      }
      //NSLog(@"eigenvalues/eigenvectors: %@",eigen);
      
      /* return RT operator */
      Matrix53 *res = [Matrix53 new];
      /* enter rotation */
      [res atRow: 0 col: 0 value: (q1*q1 + q2*q2 - q3*q3 - q4*q4)];
      [res atRow: 0 col: 1 value: (2.0*(q2*q3 + q1*q4))];
      [res atRow: 0 col: 2 value: (2.0*(q2*q4 - q1*q3))];
      [res atRow: 1 col: 0 value: (2.0*(q2*q3 - q1*q4))];
      [res atRow: 1 col: 1 value: (q1*q1 - q2*q2 + q3*q3 - q4*q4)];
      [res atRow: 1 col: 2 value: (2.0*(q3*q4 + q1*q2))];
      [res atRow: 2 col: 0 value: (2.0*(q2*q4 + q1*q3))];
      [res atRow: 2 col: 1 value: (2.0*(q3*q4 - q1*q2))];
      [res atRow: 2 col: 2 value: (q1*q1 - q2*q2 - q3*q3 + q4*q4)];
      /* enter origin */
      [res atRow: 3 col: 0 value: [com1 atRow: 0 col: 0]];
      [res atRow: 3 col: 1 value: [com1 atRow: 0 col: 1]];
      [res atRow: 3 col: 2 value: [com1 atRow: 0 col: 2]];
      /* enter translation */
      [res atRow: 4 col: 0 value: [com2 atRow: 0 col: 0]];
      [res atRow: 4 col: 1 value: [com2 atRow: 0 col: 1]];
      [res atRow: 4 col: 2 value: [com2 atRow: 0 col: 2]];

      RELEASE(pool);

#ifdef MEMDEBUG
      NSLog(@"change of allocated objects since last list\n%s",GSDebugAllocationList(YES));
      NSLog(@"allocated objects\n%s",GSDebugAllocationList(NO));
      GSDebugAllocationActive(NO);
#endif
            
      return AUTORELEASE(res);
}


/*
 *   create a selection for a chain
 */
+(Selection*)selectionWithChain: (Chain*)ch
{
      Selection *res = [Selection new];
      res->chain = RETAIN(ch);
      return AUTORELEASE(res);
}

@end


Generated by  Doxygen 1.6.0   Back to index