Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

MeasureCluster.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2019 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: MeasureCluster.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.20 $       $Date: 2020/07/23 03:27:52 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   Code to find clusters in MD trajectories.
00019  * Current implementation is based on the quality threshold (QT) algorithm:
00020  *   http://dx.doi.org/10.1101/gr.9.11.1106
00021  *   http://en.wikipedia.org/wiki/Cluster_analysis#QT_clustering_algorithm
00022  *
00023  ***************************************************************************/
00024 
00025 #include <stdio.h>
00026 #include <stdlib.h>
00027 #include <math.h>
00028 #include "Measure.h"
00029 #include "AtomSel.h"
00030 #include "utilities.h"
00031 #include "ResizeArray.h"
00032 #include "MoleculeList.h"
00033 #include "Inform.h"
00034 #include "Timestep.h"
00035 #include "VMDApp.h"
00036 #include "WKFThreads.h"
00037 #include "WKFUtils.h"
00038 #include "SpatialSearch.h"
00039 
00040 class AtomSelThr : public AtomSel
00041 {
00042 public:
00043   AtomSelThr(VMDApp *vmdapp, AtomSel *osel, wkf_mutex_t *olock)
00044       : AtomSel(vmdapp, NULL, osel->molid()),
00045         sel(osel), lock(olock) {
00046     if (sel) {
00047       selected=sel->selected;
00048       num_atoms=sel->num_atoms;
00049       which_frame=sel->which_frame;
00050       if (sel->on) {
00051         on = new int[num_atoms];
00052         memcpy(on,sel->on,num_atoms*sizeof(int));
00053       }
00054     } else {
00055       selected=-1;
00056       num_atoms=-1;
00057       which_frame=-1;
00058     }
00059   }
00060 
00061   ~AtomSelThr() {
00062     sel=NULL;
00063   }
00064 
00065   // disable these methods
00066 private:
00067 //  AtomSelThr() : AtomSel(NULL, NULL, -1) {};
00068   AtomSelThr& operator=(const AtomSelThr &) { return *this; };
00069 //  AtomSelThr(AtomSelThr &) : AtomSel(NULL, NULL, -1) {};
00070   int change(const  char *newcmd, DrawMolecule *mol) { return NO_PARSE; }
00071 
00072 public:
00073 
00074   /* thread safe selection update */
00075   void update(/* const */ DrawMolecule *mol, const int frame) {
00076     if (!sel) return;
00077     
00078     wkf_mutex_lock(lock);
00079 
00080     sel->which_frame=frame;
00081     which_frame=frame;
00082 
00083     if (sel->change(NULL, mol) != AtomSel::PARSE_SUCCESS)
00084       msgErr << "AtomSelThr::update(): failed to evaluate atom selection update";
00085 
00086     num_atoms=sel->num_atoms;
00087     selected=sel->selected;
00088     if (!on) on = new int[num_atoms];
00089     memcpy(on,sel->on,num_atoms*sizeof(int));
00090     
00091     wkf_mutex_unlock(lock);
00092   }
00093   
00094 protected:
00095   AtomSel *sel;
00096   wkf_mutex_t *lock;
00097 };
00098 
00099 /* 
00100    XXX: below is a custom version of MatrixFitRMS. unlike the
00101         original in fitrms.c, this one computes and provides
00102         the RMS and does not output the transformation matrix
00103         (not needed below).
00104 
00105         this needs to go away as soon as an improved general
00106         version of MatrixFitRMS is available, where this feature
00107         would be made an option. 
00108 */
00109 
00110 /*
00111 
00112 Code in this file was taken from PyMol v0.90 and used by permissing under
00113 the following license agreement contained in the PyMol distribution.  
00114 Trivial modifications have been made to permit incorporation into VMD.
00115 
00116 
00117 
00118 PyMOL Copyright Notice
00119 ======================
00120 
00121 The PyMOL source code is copyrighted, but you can freely use and copy
00122 it as long as you don't change or remove any of the copyright notices.
00123 
00124 ----------------------------------------------------------------------
00125 PyMOL is Copyright 1998-2003 by Warren L. DeLano of 
00126 DeLano Scientific LLC, San Carlos, CA, USA (www.delanoscientific.com).
00127 
00128                         All Rights Reserved
00129 
00130 Permission to use, copy, modify, distribute, and distribute modified 
00131 versions of this software and its documentation for any purpose and 
00132 without fee is hereby granted, provided that the above copyright 
00133 notice appear in all copies and that both the copyright notice and 
00134 this permission notice appear in supporting documentation, and that 
00135 the names of Warren L. DeLano or DeLano Scientific LLC not be used in 
00136 advertising or publicity pertaining to distribution of the software 
00137 without specific, written prior permission.
00138 
00139 WARREN LYFORD DELANO AND DELANO SCIENTIFIC LLC DISCLAIM ALL WARRANTIES 
00140 WITH REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
00141 MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL WARREN LYFORD DELANO
00142 OR DELANO SCIENTIFIC LLC BE LIABLE FOR ANY SPECIAL, INDIRECT OR 
00143 CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
00144 OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE
00145 OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE 
00146 USE OR PERFORMANCE OF THIS SOFTWARE.
00147 ----------------------------------------------------------------------
00148 
00149 Where indicated, portions of the PyMOL system are instead protected
00150 under the copyrights of the respective authors.  However, all code in
00151 the PyMOL system is released as non-restrictive open-source software
00152 under the above license or an equivalent license.  
00153 
00154 PyMOL Trademark Notice
00155 ======================
00156 
00157 PyMOL(TM) is a trademark of DeLano Scientific LLC.  Derivate software
00158 which contains PyMOL source code must be plainly distinguished from
00159 the PyMOL package distributed by DeLano Scientific LLC in all publicity,
00160 advertising, and documentation.
00161 
00162 The slogans, "Includes PyMOL(TM).", "Based on PyMOL(TM) technology.",
00163 "Contains PyMOL(TM) source code.", and "Built using PyMOL(TM).", may
00164 be used in advertising, publicity, and documentation of derivate
00165 software provided that the notice, "PyMOL is a trademark of DeLano
00166 Scientific LLC.", is included in a footnote or at the end of the document.
00167 
00168 All other endorsements employing the PyMOL trademark require specific,
00169 written prior permission.
00170 
00171 --Warren L. DeLano (warren@delanoscientific.com)
00172 
00173 */
00174 
00175 #ifdef __cplusplus
00176 extern "C" {
00177 #endif
00178 
00179 #ifdef R_SMALL
00180 #undef R_SMALL
00181 #endif
00182 #define R_SMALL 0.000000001
00183 
00184 static void normalize3d(double *v) {
00185   double vlen;
00186   vlen = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
00187   if (vlen > R_SMALL) {
00188     v[0] /= vlen;
00189     v[1] /= vlen;
00190     v[2] /= vlen;
00191   } else {
00192     v[0] = 0;
00193     v[1] = 0;
00194     v[2] = 0;
00195   }
00196 }
00197 
00198 /*========================================================================*/
00199   static float MyMatrixFitRMS(int n, float *v1, float *v2, const float *wt, const double tol)
00200 {
00201   /*
00202         Subroutine to do the actual RMS fitting of two sets of vector coordinates
00203         This routine does not rotate the actual coordinates, but instead returns 
00204         the RMS fitting value, along with the center-of-mass translation vectors 
00205         T1 and T2 and the rotation vector M, which rotates the translated 
00206         coordinates of molecule 2 onto the translated coordinates of molecule 1.
00207   */
00208 
00209   float *vv1,*vv2;
00210   double m[3][3],aa[3][3];
00211   double sumwt, sig, gam;
00212   double sg, bb, cc, tmp, err, etmp;
00213   int a, b, c, maxiter, iters, iy, iz;
00214   double t1[3],t2[3];
00215   double aatmp[9];
00216 
00217   /* Initialize arrays. */
00218 
00219   for(a=0;a<3;a++) {
00220     for(b=0;b<3;b++) {
00221       m[a][b] = 0.0;
00222       aa[a][b] = 0.0;
00223       aatmp[3*a+b] = 0;
00224     }
00225     m[a][a] = 1.0;
00226     t1[a]=0.0;
00227     t2[a]=0.0;
00228   }
00229 
00230   /* maximum number of fitting iterations */
00231   maxiter = 1000;
00232 
00233   /* Calculate center-of-mass vectors */
00234 
00235   vv1=v1;
00236   vv2=v2;
00237   sumwt = 0.0;
00238 
00239   for(c=0;c<n;c++) {
00240     double w = wt ? wt[c] : 1;
00241     t1[0] += w * vv1[0];
00242     t1[1] += w * vv1[1];
00243     t1[2] += w * vv1[2];
00244     t2[0] += w * vv2[0];
00245     t2[1] += w * vv2[1];
00246     t2[2] += w * vv2[2];
00247     sumwt += w;
00248     vv1+=3;
00249     vv2+=3;
00250   }
00251   for(a=0;a<3;a++) {
00252     t1[a] /= sumwt;
00253     t2[a] /= sumwt;
00254   }
00255 
00256   /* Calculate correlation matrix */
00257   vv1=v1;
00258   vv2=v2;
00259   for(c=0;c<n;c++) {
00260     double w = wt ? wt[c] : 1;
00261     double x1 = w * (vv1[0] - t1[0]);
00262     double y1 = w * (vv1[1] - t1[1]);
00263     double z1 = w * (vv1[2] - t1[2]);
00264 
00265     /* don't multply x2/y2/z2 by w, otherwise weights get squared */
00266     double x2 =     (vv2[0] - t2[0]); 
00267     double y2 =     (vv2[1] - t2[1]);
00268     double z2 =     (vv2[2] - t2[2]);
00269     aatmp[0] += x2 * x1;
00270     aatmp[1] += x2 * y1;
00271     aatmp[2] += x2 * z1;
00272     aatmp[3] += y2 * x1;
00273     aatmp[4] += y2 * y1;
00274     aatmp[5] += y2 * z1;
00275     aatmp[6] += z2 * x1;
00276     aatmp[7] += z2 * y1;
00277     aatmp[8] += z2 * z1;
00278     vv1+=3;
00279     vv2+=3;
00280   }
00281 
00282   for (a=0; a<3; a++) 
00283     for (b=0; b<3; b++) 
00284       aa[a][b] = aatmp[3*a+b];
00285 
00286   if(n>1) {
00287     /* Primary iteration scheme to determine rotation matrix for molecule 2 */
00288     iters = 0;
00289     while(1) {
00290       /* IX, IY, and IZ rotate 1-2-3, 2-3-1, 3-1-2, etc.*/
00291       iz = (iters+1) % 3;
00292       iy = (iz+1) % 3;
00293       // unused...
00294       // ix = (iy+1) % 3;
00295       sig = aa[iz][iy] - aa[iy][iz];
00296       gam = aa[iy][iy] + aa[iz][iz];
00297 
00298       if(iters>=maxiter) {
00299         fprintf(stderr,
00300                 " Matrix: Warning: no convergence (%.15f>%.15f after %d iterations).\n",
00301                 fabs(sig),tol*fabs(gam),iters);
00302         break;
00303       }
00304 
00305       /* Determine size of off-diagonal element.  If off-diagonals exceed the
00306          diagonal elements * tolerance, perform Jacobi rotation. */
00307       tmp = sig*sig + gam*gam;
00308       sg = sqrt(tmp);
00309       if( (sg > 0.0) && (fabs(sig)>(tol*fabs(gam))) ) {
00310         sg = 1.0 / sg;
00311         for(a=0;a<3;a++) {
00312           bb = gam*aa[iy][a] + sig*aa[iz][a];
00313           cc = gam*aa[iz][a] - sig*aa[iy][a];
00314           aa[iy][a] = bb*sg;
00315           aa[iz][a] = cc*sg;
00316 
00317           bb = gam*m[iy][a] + sig*m[iz][a];
00318           cc = gam*m[iz][a] - sig*m[iy][a];
00319           m[iy][a] = bb*sg;
00320           m[iz][a] = cc*sg;
00321         }
00322       } else {
00323         break;
00324       }
00325       iters++;
00326     }
00327   }
00328   /* At this point, we should have a converged rotation matrix (M).  Calculate
00329          the weighted RMS error. */
00330   err=0.0;
00331   vv1=v1;
00332   vv2=v2;
00333 
00334   normalize3d(m[0]);
00335   normalize3d(m[1]);
00336   normalize3d(m[2]);
00337 
00338   for(c=0;c<n;c++) {
00339         etmp = 0.0;
00340         for(a=0;a<3;a++) {
00341           tmp = m[a][0]*(vv2[0]-t2[0])
00342                 + m[a][1]*(vv2[1]-t2[1])
00343                 + m[a][2]*(vv2[2]-t2[2]);
00344           tmp = (vv1[a]-t1[a])-tmp;
00345           etmp += tmp*tmp;
00346         }
00347 
00348         if(wt)
00349           err += wt[c] * etmp;
00350         else 
00351           err += etmp;
00352 
00353         vv1+=3;
00354         vv2+=3;
00355   }
00356 
00357   err=err/sumwt;
00358   err=sqrt(err);
00359   return (float)err;
00360 }
00361 
00362 #ifdef __cplusplus
00363 }
00364 #endif
00365 
00366 /* XXX: end of customized MatrixFitRMS */
00367 
00368 
00369 // compute weighted RMSD between selected atoms in two frames
00370 // this is a simple wrapper around measure_rmsd.
00371 static float cluster_get_rmsd(const float *Frame1Pos, const float *Frame2Pos, 
00372                               AtomSel *sel, float *weights) {
00373   float distance = 0.0f;
00374   measure_rmsd(sel, sel, sel->num_atoms, Frame1Pos, Frame2Pos, weights, &distance);
00375   return distance;
00376 }
00377 
00378 
00379 // compute weighted difference between radius of gyration
00380 // of the selected atoms in the two frames.
00381 static float cluster_get_rgyrd(const float *Frame1Pos, const float *Frame2Pos, 
00382                                AtomSel *sel, float *weights) {
00383 
00384   float distance = 10000000.0f;
00385 
00386   // compute the center of mass with the current weights
00387   float com1[3], com2[3];
00388   int ret_val;
00389 
00390   ret_val = measure_center(sel, Frame1Pos, weights, com1);
00391   if (ret_val < 0) 
00392     return distance;
00393 
00394   ret_val = measure_center(sel, Frame2Pos, weights, com2);
00395   if (ret_val < 0) 
00396     return distance;
00397   
00398   // measure center of gyration
00399   int i, j;
00400   float total_w, w, sum1, sum2;
00401   total_w=sum1=sum2=0.0f;
00402   for (j=0,i=sel->firstsel; i<=sel->lastsel; i++) {
00403     if (sel->on[i]) {
00404       w = weights[j];
00405       total_w += w;
00406       sum1 += w * distance2(Frame1Pos + 3*i, com1);
00407       sum2 += w * distance2(Frame2Pos + 3*i, com2);
00408       j++;
00409     }
00410   }
00411 
00412   if (total_w == 0.0f)
00413     return distance;
00414 
00415   // and finalize the computation
00416   distance = sqrtf(sum1/total_w) - sqrtf(sum2/total_w);
00417   return fabsf(distance);
00418 }
00419 
00420 
00421 // This is a stripped down version of measure fit supporting only
00422 // selections with 4 atoms or larger. not much value in clustering
00423 // smaller systems with fit and rmsd.
00424 // This algorithm comes from Kabsch, Acta Cryst. (1978) A34, 827-828.
00425 static float cluster_get_fitrmsd(const float *Frame1Pos, const float *Frame2Pos, 
00426                                  AtomSel *sel, float *weights, const double tol) {
00427 
00428   int num = sel->selected;
00429   
00430   // failure of the fit+rmsd is indicated by a very large distance.
00431   float distance = 10000000.0f;
00432 
00433   // use the new RMS fit implementation only
00434   // fit+clustering with 3 or less atoms doesn't make much sense.
00435   // the Kabsch method won't work of the number of atoms is less than 4
00436   // (and won't work in some cases of n > 4; I think it works so long as
00437   // three or more planes are needed to intersect all the data points
00438 
00439   if (sel->selected < 4)
00440     return distance;
00441     
00442   int i, j, k;
00443   float *v1, *v2, *wt;
00444   v1 = new float[3*num];
00445   v2 = new float[3*num];
00446   wt = new float[num];
00447   for (j=0,k=0,i=sel->firstsel; i<=sel->lastsel; i++) {
00448     if (sel->on[i]) {
00449       int ind = 3 * i;
00450       wt[j] = weights[i];
00451       ++j;
00452       v1[k] = Frame1Pos[ind];
00453       v2[k] = Frame2Pos[ind];
00454       v1[k+1] = Frame1Pos[ind+1];
00455       v2[k+1] = Frame2Pos[ind+1];
00456       v1[k+2] = Frame1Pos[ind+2];
00457       v2[k+2] = Frame2Pos[ind+2];
00458       k+=3;
00459     }
00460   }
00461   distance = MyMatrixFitRMS(num, v1, v2, wt, tol);
00462 
00463   delete [] v1;
00464   delete [] v2;
00465   delete [] wt;
00466 
00467   return distance;
00468 }
00469 
00470 
00471 
00472 
00473 typedef struct {
00474   int threadid;
00475   int threadcount;
00476 
00477   int max_cluster_size;
00478   const int *skip_list;
00479   int *new_skip_list;
00480   int *max_cluster;
00481 
00482   int istart;
00483   int iend;
00484   int *frames_list;
00485   int numframes;
00486 
00487   AtomSelThr *sel;
00488   Molecule *mol;
00489   int selupdate;
00490   float cutoff;
00491   int likeness;
00492   float *weights;
00493 } clusterparms_t;
00494 
00495 
00496 // cluster search thread worker function
00497 extern "C" void * find_cluster_thr(void *voidparms)
00498 {
00499 
00500   clusterparms_t *parms = (clusterparms_t *)voidparms;
00501   const int istart = parms->istart;
00502   const int iend   = parms->iend;
00503   int *framesList = parms->frames_list;
00504   const int numframes = parms->numframes;
00505 
00506   const int selupdate = parms->selupdate;
00507   const int likeness = parms->likeness;
00508   float cutoff = parms->cutoff;
00509   float *weights = parms->weights;
00510 
00511   AtomSelThr *sel = parms->sel;
00512   Molecule *mymol = parms->mol;
00513   const int *skipList = parms->skip_list;
00514 
00515   int *maxCluster = parms->max_cluster;
00516   memset(maxCluster, 0, numframes*sizeof(int));
00517   int *newSkipList = parms->new_skip_list;
00518   memset(newSkipList, 0, numframes*sizeof(int));
00519 
00520   int maxClusterSize = 0, tempClusterSize = 0;
00521   int *tempCluster = new int[numframes];
00522   int *tempSkipList = new int[numframes];
00523 
00524 
00525   // MatrixFitRMS returns RMS distance of fitted molecule. 
00526   /* RMS fit tolerance */
00527   double tol = 1e-15;
00528   const char *TOL = getenv( "VMDFITRMSTOLERANCE" );
00529   if (TOL)
00530     tol = atof(TOL);
00531 
00532   // Loops through assigned frames find the one with the max cluster size
00533   int i,j;
00534   for (i = istart; i < iend; i++) {
00535     memset(tempSkipList, 0, numframes*sizeof(int));
00536     memset(tempCluster, 0, numframes*sizeof(int));
00537 
00538     if (skipList[i]==0) {
00539       if (selupdate)
00540         sel->update(mymol,framesList[i]);
00541       
00542       const Timestep *tsMain = mymol->get_frame(framesList[i]);
00543       const float *framePos = tsMain->pos;
00544 
00545       tempCluster[0] = i;
00546       tempSkipList[i] = 1;
00547       tempClusterSize = 1;
00548 
00549       // Loops through all frames other then frame i and computes frame i's distance to them
00550       for (j = 0; j < numframes; j++) {
00551         if (skipList[j]==0 && j != i) {
00552           const Timestep *ts2;
00553           ts2 = mymol->get_frame(framesList[j]);
00554           float distance;
00555 
00556           // branch to the implemented likeness algorithms
00557           switch(likeness) {
00558             case MEASURE_DIST_RMSD:
00559               distance = cluster_get_rmsd(framePos, ts2->pos, sel, weights);
00560               break;
00561 
00562             case MEASURE_DIST_FITRMSD:
00563               distance = cluster_get_fitrmsd(framePos, ts2->pos, sel, weights, tol);
00564               break;
00565               
00566             case MEASURE_DIST_RGYRD:
00567               distance = cluster_get_rgyrd(framePos, ts2->pos, sel, weights);
00568               break;
00569               
00570             default:
00571               distance = 10000000.0f;
00572           }
00573 
00574           if (distance <= cutoff) {
00575             tempCluster[tempClusterSize] = j;
00576             ++tempClusterSize;
00577             tempSkipList[j] = 1;
00578           }
00579         }
00580       }
00581 
00582       // If size of temp cluster > max cluster, temp cluster becomes max cluster
00583       if (tempClusterSize > maxClusterSize) {
00584         int *temp;
00585         maxClusterSize = tempClusterSize;
00586 
00587         temp = maxCluster;
00588         maxCluster = tempCluster;
00589         tempCluster = temp;
00590 
00591         temp = newSkipList;
00592         newSkipList = tempSkipList;
00593         tempSkipList = temp;
00594       }
00595     }
00596   }
00597 
00598   // update parameter struct with results
00599   parms->max_cluster_size = maxClusterSize;
00600   parms->max_cluster = maxCluster;
00601   parms->new_skip_list = newSkipList;
00602   
00603   // cleanup
00604   delete[] tempCluster;
00605   delete[] tempSkipList;
00606 
00607   return MEASURE_NOERR;
00608 }
00609     
00610 
00612 static int *find_next_cluster(Molecule *mymol, int *framesList, const int numframes, 
00613                               const int remframes, const int *skipList, int **newSkipList,
00614                               const int likeness, AtomSel *sel, const int selupdate, 
00615                               const double cutoff, float *weights)
00616 {
00617   int i,j;
00618   
00619   // threading setup.
00620   wkf_thread_t   *threads;
00621   clusterparms_t *parms;
00622 
00623 #if defined(VMDTHREADS)
00624   int numprocs = wkf_thread_numprocessors();
00625 #else
00626   int numprocs = 1;
00627 #endif
00628 
00629   int delta = remframes / numprocs;
00630   int istart = 0;
00631   int iend = 0;
00632   
00633   // not enough work to do, force serial execution
00634   if (delta < 1) {
00635     numprocs=1;
00636     delta=numframes;
00637   }
00638 
00639   threads = new wkf_thread_t[numprocs];
00640   memset(threads, 0, numprocs * sizeof(wkf_thread_t));
00641   wkf_mutex_t *atomsel_lock = new wkf_mutex_t;
00642   wkf_mutex_init(atomsel_lock);
00643 
00644   // allocate and (partially) initialize array of per-thread parameters
00645   parms = new clusterparms_t[numprocs];
00646   for (i=0; i<numprocs; ++i) {
00647     parms[i].threadid = i;
00648     parms[i].threadcount = numprocs;
00649 
00650     parms[i].max_cluster_size = 1;
00651     parms[i].skip_list = skipList;
00652     parms[i].new_skip_list = new int[numframes];
00653     parms[i].max_cluster = new int[numframes];
00654 
00655     // use a thread-safe wrapper to access "the one" 
00656     // AtomSel class. The wrapper uses mutexes to 
00657     // prevent from updating the global selection from
00658     // multiple threads at the same time. The whole data 
00659     // access infrastructure in VMD is currently not thread-safe.
00660     parms[i].sel = new AtomSelThr(mymol->app, sel, atomsel_lock);
00661     parms[i].mol = mymol;
00662 
00663     // load balancing. scatter the remaining frames evenly
00664     // by skipping over eliminated frames.
00665     parms[i].istart = istart;
00666     int nframe=0;
00667     for (j=istart; (j < numframes) && (nframe < delta); ++j) {
00668       if (skipList[framesList[j]]==0) 
00669         ++nframe;
00670       iend=j;
00671     }
00672     parms[i].iend = iend;
00673     istart=iend;
00674                    
00675     parms[i].frames_list = framesList;
00676     parms[i].numframes = numframes;
00677     parms[i].selupdate = selupdate;
00678     parms[i].likeness = likeness;
00679     parms[i].cutoff = (float) cutoff;
00680     parms[i].weights = weights;
00681   }
00682   parms[numprocs-1].iend=numframes;
00683 
00684 
00685 #if defined(VMDTHREADS)
00686   if (numprocs > 1) {
00687     for (i=0; i<numprocs; ++i) {
00688       wkf_thread_create(&threads[i], find_cluster_thr, &parms[i]);
00689     }
00690     for (i=0; i<numprocs; ++i) {
00691       wkf_thread_join(threads[i], NULL);
00692     }
00693   } else
00694 #endif
00695   find_cluster_thr(&parms[0]);
00696   
00697   int maxClusterSize = parms[0].max_cluster_size;
00698   int *maxCluster = parms[0].max_cluster;
00699   delete[] *newSkipList;
00700   *newSkipList= parms[0].new_skip_list;
00701 
00702   // retrieve results from additional threads,
00703   // override, if needed, and free temporary storage.
00704   if (numprocs > 1) {
00705     for (i = 1; i < numprocs; i++) {
00706       if (parms[i].max_cluster_size > maxClusterSize) {
00707         maxClusterSize = parms[i].max_cluster_size;
00708         delete[] maxCluster;
00709         maxCluster = parms[i].max_cluster;
00710         delete[] *newSkipList;
00711         *newSkipList = parms[i].new_skip_list;
00712       } else {
00713         delete[] parms[i].max_cluster;
00714         delete[] parms[i].new_skip_list;
00715       }
00716     }
00717   }
00718     
00719   // Transform cluster list back to real frame numbers
00720   for (i = 0; i < numframes; i++) {
00721     maxCluster[i] = framesList[maxCluster[i]];
00722   }
00723 
00724   // cleanup.
00725   wkf_mutex_destroy(atomsel_lock);
00726   delete atomsel_lock;
00727 
00728   if (selupdate) {
00729     for (i=0; i<numprocs; ++i)
00730       delete parms[i].sel;
00731   }
00732   delete[] threads;
00733   delete[] parms;
00734 
00735   return maxCluster;
00736 }
00737 
00738 int measure_cluster(AtomSel *sel, MoleculeList *mlist,
00739                     const int numcluster, const int algorithm,
00740                     const int likeness, const double cutoff,
00741                     int *clustersize, int **clusterlist,
00742                     int first, int last, int step, int selupdate, 
00743                     float *weights)
00744 {
00745   Molecule *mymol = mlist->mol_from_id(sel->molid());
00746   int maxframe = mymol->numframes()-1;
00747 
00748   if (last == -1) last = maxframe;
00749 
00750   if ((last < first) || (last < 0) || (step <=0) || (first < 0)
00751       || (last > maxframe)) {
00752     msgErr << "measure cluster: bad frame range given."
00753            << " max. allowed frame#: " << maxframe << sendmsg;
00754     return MEASURE_ERR_BADFRAMERANGE;
00755   }
00756 
00757   int numframes = (last-first+1)/step;
00758   int remframes = numframes;
00759 
00760   // create list with frames numbers selected to process
00761   int *framesList = new int[numframes];
00762   int frame_count = 0;
00763   int n;
00764 
00765   for(n = first; n <= last; n += step)
00766     framesList[frame_count++] = n;
00767 
00768   // accumulated list of frames to skip because they belong to a cluster
00769   int *skipList = new int[numframes];
00770   // new list of frames to skip to be added to the existing skip list.
00771   int *newSkipList = new int[numframes];
00772   // initially we want all frames.
00773   memset(skipList, 0, numframes*sizeof(int));
00774   // timer for progress messages
00775   wkfmsgtimer *msgtp = wkf_msg_timer_create(5);
00776 
00777   // compute the numcluster largest clusters
00778   for(n = 0; n < numcluster; ++n){
00779 
00780     // wipe out list of frames to be added to the global skiplist
00781     memset(newSkipList, 0, numframes*sizeof(int));
00782     clusterlist[n] = find_next_cluster(mymol, framesList, numframes, remframes,
00783                                        skipList, &newSkipList, likeness,
00784                                        sel, selupdate, cutoff, weights);
00785     int n_cluster=0;
00786     for(int i = 0; i < numframes; ++i){
00787       if (newSkipList[i] == 1) {
00788         skipList[i] = 1;
00789         n_cluster++;
00790       }
00791     }
00792     clustersize[n]=n_cluster;
00793     remframes -= n_cluster;
00794 
00795     // print progress messages for long running tasks
00796     if (msgtp && wkf_msg_timer_timeout(msgtp)) {
00797       char tmpbuf[1024];
00798       sprintf(tmpbuf, "cluster %d of %d: (%6.2f%% complete). %d frames of %d left.", 
00799               n+1, numcluster, 100.0f*(n+1)/((float) numcluster), remframes, numframes);
00800       msgInfo << "measure cluster: " << tmpbuf << sendmsg;
00801     }
00802   }
00803 
00804   // combine unclustered frames to form the last cluster
00805   int *unclustered = new int[numframes];
00806   int numunclustered = 0;
00807   for (n = 0; n < numframes; ++n) {
00808     if (skipList[n] == 0) {
00809       unclustered[numunclustered] = framesList[n];
00810       ++numunclustered;
00811     }
00812   }
00813   // NOTE: both lists have been allocated
00814   // to the length of numcluster+1
00815   clusterlist[numcluster] = unclustered;
00816   clustersize[numcluster] = numunclustered;
00817 
00818   // Cleanup
00819   delete[] newSkipList;
00820   delete[] skipList;
00821   wkf_msg_timer_destroy(msgtp);
00822 
00823   return MEASURE_NOERR;
00824 }
00825 
00826 
00827 /**************************************************************************/
00828 
00829 typedef ResizeArray<int> intlist;
00830 
00831 // helper function for cluster size analysis
00832 // build index list for the next cluster.
00833 static void assemble_cluster(intlist &cluster_list, intlist &candidate_list,
00834                intlist **neighbor_grid, int atom, int numshared, int *idxmap) {
00835   int idx, nn, i,j;
00836   
00837   // clear lists and add initial atom pairs to candidates list.
00838   candidate_list.clear();
00839   cluster_list.clear();
00840 
00841   idx = idxmap[atom];
00842   nn  = neighbor_grid[idx]->num();
00843   for (i = 0; i < nn; i++) {
00844     int bn = (*neighbor_grid[idx])[i];
00845     if (neighbor_grid[idxmap[bn]]) {
00846       candidate_list.append(atom);
00847       candidate_list.append(bn);
00848     }
00849   }
00850 
00851   // pointer to the currently processed cluster candidate list entry.
00852   int curidx=0;
00853   
00854   while (curidx < candidate_list.num()) {
00855 
00856     // a pair of neighbors has to share at least numshared
00857     // neighbors to be added to the cluster.
00858     // at least numshared neighbors.
00859     int count = 0;
00860 
00861     if (numshared > 0) {
00862       int ida = idxmap[candidate_list[curidx]];
00863       int idb = idxmap[candidate_list[curidx+1]];
00864       int nna  = neighbor_grid[ida]->num();
00865       int nnb  = neighbor_grid[idb]->num();
00866 
00867       for (i = 0; i < nna; i++) {
00868         if (neighbor_grid[ida]) {
00869           for (j = 0; j < nnb; j++) {
00870             if (neighbor_grid[idb]) {
00871               if ( (*neighbor_grid[ida])[i] == (*neighbor_grid[idb])[j] ) {
00872                 ++count;
00873                 if (count == numshared) 
00874                   goto exit;
00875               }
00876             }
00877           }
00878         }
00879       }
00880     }
00881     exit:
00882     
00883     if (count == numshared) {
00884 
00885       // add central atom of group of neighbors.
00886       // its neighbors had already been added to
00887       // the candidate list.
00888       int atma = candidate_list[curidx];
00889       if (cluster_list.find(atma) < 0) 
00890         cluster_list.append(atma);
00891 
00892       // add neighbor of central atom to cluster and
00893       // add neighbors of this atom to candidate list, 
00894       // if they are not already in it.
00895       int atmb = candidate_list[curidx+1];
00896       idx = idxmap[atmb];
00897 
00898       if (cluster_list.find(atmb) < 0) {
00899         cluster_list.append(atmb);
00900         
00901         int nnb = neighbor_grid[idx]->num();
00902         for (i = 0; i < nnb; i++) {
00903           int bn = (*neighbor_grid[idx])[i];
00904           if ((neighbor_grid[idxmap[bn]]) && (cluster_list.find(bn) < 0)) {
00905             candidate_list.append(atmb);
00906             candidate_list.append(bn);
00907           }
00908         }
00909       }
00910     }
00911 
00912     ++curidx;++curidx; // next candidate pair
00913   }
00914   
00915   return;
00916 }
00917 
00918 // perform cluster size analysis
00919 int measure_clustsize(const AtomSel *sel, MoleculeList *mlist,
00920                       const double cutoff, int *clustersize,
00921                       int *clusternum, int *clusteridx, 
00922                       int minsize, int numshared, int usepbc) {
00923   int i,j;
00924 
00925   const float *framepos  = sel->coordinates(mlist);
00926   const int num_selected = sel->selected;
00927   const int num_atoms    = sel->num_atoms;
00928   const int *selected    = sel->on;
00929 
00930   // forward and reverse index maps for relative position
00931   // of atoms in the selection and in the global arrays.
00932   int *idxmap      = new int[num_atoms];
00933   int *idxrev      = new int[num_selected];
00934 
00935   for (i=0; i<sel->firstsel; i++) 
00936     idxmap[i]=-1;
00937 
00938   for(j=0,i=sel->firstsel; i<=sel->lastsel; i++) {
00939     if (sel->on[i]) {
00940       idxrev[j]=i;
00941       idxmap[i]=j++;
00942     } else {
00943       idxmap[i]=-1;
00944     }
00945   }
00946 
00947   for (i=sel->lastsel+1; i<sel->num_atoms; i++) 
00948     idxmap[i]=-1;
00949 
00950   // allocate list of neighbor lists.
00951   intlist **neighbor_grid;
00952   neighbor_grid = new intlist *[num_selected];
00953   for (i = 0; i < num_selected; i++)
00954     neighbor_grid[i] = new intlist;
00955 
00956   // compile list of pairs for selection.
00957   GridSearchPair *pairlist, *currentPair, *nextPair;
00958   pairlist = vmd_gridsearch1(framepos, num_atoms, selected, (float)cutoff, 0, -1);
00959 
00960   // populate the neighborlist grid.
00961   for (currentPair = pairlist; currentPair != NULL; currentPair = nextPair) {
00962     neighbor_grid[idxmap[currentPair->ind1]]->append(currentPair->ind2);
00963     neighbor_grid[idxmap[currentPair->ind2]]->append(currentPair->ind1);
00964     nextPair = currentPair->next;
00965     free(currentPair);
00966   }
00967 
00968   // collect the cluster size information.
00969   int currentClusterNum = 0;
00970   int currentPosition = 0;
00971   intlist cluster_list(64);
00972   intlist candidate_list(128);
00973 
00974   while (currentPosition < num_selected) {
00975     if (neighbor_grid[currentPosition]) {
00976       // pick next atom that has not been processed yet and build list of 
00977       // all indices for this cluster. by looping over the neighbors of 
00978       // the first atom and adding all pairs of neighbors that share at
00979       // least numshared neighbors. continue with neighbors of neighbors 
00980       // accordingly until no more new unique neighbors are found.
00981       // entries of atoms added to a cluster are removed from neighbor_grid
00982       if (neighbor_grid[currentPosition]->num() > numshared) {
00983         assemble_cluster(cluster_list, candidate_list, neighbor_grid, 
00984                          idxrev[currentPosition], numshared, idxmap);
00985 
00986         if (minsize <= cluster_list.num()) {
00987           // these atoms have been processed. remove from global list
00988           for (i = 0; i < cluster_list.num(); i++) {
00989             int idx = idxmap[cluster_list[i]];
00990             delete neighbor_grid[idx];
00991             neighbor_grid[idx] = 0;
00992           }
00993       
00994           // store the cluster size, cluster index and atom index information in 
00995           // the designated arrays.
00996           for (i = 0; i < cluster_list.num(); i++) {
00997             int atom = idxmap[cluster_list[i]];
00998             clusteridx[atom] = cluster_list[i];
00999             clusternum[atom] = currentClusterNum;
01000             clustersize[atom] = cluster_list.num();
01001           }
01002           currentClusterNum++;
01003         }
01004       }
01005     }
01006     ++currentPosition;
01007   }
01008 
01009   for(i=0; i < num_selected; ++i) {
01010     if (neighbor_grid[i])
01011       delete neighbor_grid[i];
01012   }
01013   delete[] neighbor_grid;
01014 
01015   return MEASURE_NOERR;
01016 }

Generated on Thu Mar 28 02:43:23 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002