fpocket: cluster.c Source File

00001 
00002 #include "../headers/cluster.h"
00003 
00004 /*
00005 
00006 ## GENERAL INFORMATION
00007 ##
00008 ## FILE                                         cluster.h
00009 ## AUTHORS                                      P. Schmidtke and V. Le Guilloux
00010 ## LAST MODIFIED                        28-11-08
00011 ##
00012 ## SPECIFICATIONS
00013 ##
00014 ##      This file contains currently only one function, providing
00015 ##      a mutliple linkage clustering algorithm performed on a list
00016 ##      of pockets.
00017 ##
00018 ## MODIFICATIONS HISTORY
00019 ##
00020 ##      19-11-08        (p)  Extension of comments, change in multiple linkage clustering
00021 ##      28-11-08        (v)  Comments UTD + minor relooking
00022 ##      11-05-08        (v)  singleLinkageClustering -> pck_sl_clust
00023 ##      01-04-08        (v)  Added comments and creation of history
00024 ##      01-01-08        (vp) Created (random date...)
00025 ##
00026 ## TODO or SUGGESTIONS
00027 ##
00028 ##      (v) Possible improvement:
00029 ##              Use the sorted structure to find neighbors in a more
00030 ##              efficient way.
00031 ##      (v) Rename the file ! (mlcluster.c for example...)
00032 ##              Or maybe move this function into pocket.c, as the
00033 ##              algorithm deals with pockets only...
00034 ##      (v) Check and update if necessary comments of each function!!
00035 
00036 */
00037 
00038 /*
00039     COPYRIGHT DISCLAIMER
00040 
00041     Vincent Le Guilloux, Peter Schmidtke and Pierre Tuffery, hereby
00042         disclaim all copyright interest in the program “fpocket” (which
00043         performs protein cavity detection) written by Vincent Le Guilloux and Peter
00044         Schmidtke.
00045 
00046     Vincent Le Guilloux  28 November 2008
00047     Peter Schmidtke      28 November 2008
00048     Pierre Tuffery       28 November 2008
00049 
00050     GNU GPL
00051 
00052     This file is part of the fpocket package.
00053 
00054     fpocket is free software: you can redistribute it and/or modify
00055     it under the terms of the GNU General Public License as published by
00056     the Free Software Foundation, either version 3 of the License, or
00057     (at your option) any later version.
00058 
00059     fpocket is distributed in the hope that it will be useful,
00060     but WITHOUT ANY WARRANTY; without even the implied warranty of
00061     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00062     GNU General Public License for more details.
00063 
00064     You should have received a copy of the GNU General Public License
00065     along with fpocket.  If not, see <http://www.gnu.org/licenses/>.
00066 
00067 */
00068 
00069 
00070 void pck_final_clust(c_lst_pockets *pockets, s_fparams *params,s_pdb *pdb,s_pdb *pdb_w_lig)
00071 {
00072         node_pocket *pcur = NULL,
00073                                 *pnext = NULL ,
00074                                 *curMobilePocket = NULL ;
00075 
00076         node_vertice *vcur = NULL ;
00077         node_vertice *curMobileVertice = NULL ;
00078 
00079         s_vvertice *vvcur = NULL,
00080                            *mvvcur = NULL ;
00081         float vcurx,
00082                   vcury,
00083                   vcurz ;
00084         float curdist;
00085         float curasphdens,dens1,dens2;
00086         float **dmat;   /*distance matrix*/
00087         size_t i,j;
00088         dmat=(float **)malloc(sizeof(float *)*pockets->n_pockets);
00089         for(i=0;i<pockets->n_pockets;i++) {
00090             dmat[i]=(float *) malloc(sizeof(float)*pockets->n_pockets);
00091             for(j=0;j<pockets->n_pockets;j++) dmat[i][j]=0.0;
00092         }
00093 
00094         /* Flag to know if two clusters are next to each other by single linkage
00095          * clustering...or not */
00096         int nflag ;
00097 
00098         if(!pockets) {
00099                 fprintf(stderr, "! Incorrect argument during Multiple Linkage Clustering.\n") ;
00100                 return ;
00101         }
00102 
00103         /* Set the first pocket */
00104         set_pockets_descriptors(pockets,pdb,params,pdb_w_lig);
00105         pcur = pockets->first ;
00106         //fprintf(stdout,"\nHaving %d comparisons\n",pockets->n_pockets*pockets->n_pockets);
00107         i=0;
00108         size_t n_slist=0;
00109         while(pcur) {
00110             j=i+1;
00111                 /* Set the second pocket */
00112                 curMobilePocket = pcur->next ;
00113                 while(curMobilePocket) {
00114                         curdist=0.0;
00115 
00116                         nflag = 0 ;
00117                         /* Set the first vertice/alpha sphere center of the first pocket */
00118                         vcur = pcur->pocket->v_lst->first ;
00119                         while(vcur){
00120                                 /* Set the first vertice/alpha sphere center of the second pocket */
00121                                 curMobileVertice = curMobilePocket->pocket->v_lst->first ;
00122                                 vvcur = vcur->vertice ;
00123                                 vcurx = vvcur->x ;
00124                                 vcury = vvcur->y ;
00125                                 vcurz = vvcur->z ;
00126 
00127                                 /* Double loop for vertices -> if not near */
00128                                 while(curMobileVertice){
00129                                         mvvcur = curMobileVertice->vertice ;
00130                                         if(dist(vcurx, vcury, vcurz, mvvcur->x, mvvcur->y, mvvcur->z)<params->sl_clust_max_dist) curdist-=1.0;
00131                                         curMobileVertice = curMobileVertice->next;
00132                                 }
00133                                 vcur = vcur->next ;
00134                         }
00135 
00136                         pnext =  curMobilePocket->next ;
00137                         /* If the distance flag has counted enough occurences of near neighbours, merge pockets*/
00138                                 /* If they are next to each other, merge them */
00139                                 //mergePockets(pcur,curMobilePocket,pockets);
00140                         //fprintf(stdout,"\ni %d j %d\n",i,j),
00141 
00142                         dens1=pcur->pocket->pdesc->as_density;
00143                         dens1=((isnan(dens1)) ? 0.0 : dens1);
00144 
00145                         dens2=curMobilePocket->pocket->pdesc->as_density;
00146                         dens2=((isnan(dens2)) ? 0.0 : dens2);
00147                         curasphdens=0.01*(dens1-dens2)*(dens1-dens2);
00148                         curasphdens+=(exp(0.1*dens1)+exp(0.1*dens2))/2.0-1.0;
00149                         //curasphdens=0.0;
00150                         dmat[i][j]=curdist;
00151                         dmat[j][i]=curdist;
00152                         curMobilePocket = pnext ;
00153                         n_slist++;
00154                         j++;
00155                 }
00156                 pcur = pcur->next ;
00157                 i++;
00158         }
00159 
00160         /* Now we have to merge nearby pockets without loosing track*/
00161 
00162         /*create a chained list with track on */
00163         n_slist=((pockets->n_pockets*pockets->n_pockets)-pockets->n_pockets)/2;
00164         s_sorted_pocket_list *slist=NULL;
00165         slist=(s_sorted_pocket_list *)my_malloc(sizeof(s_sorted_pocket_list)*n_slist);
00166         //s_sorted_pocket_list slist[n_slist];
00167         for(i=0;i<n_slist;i++) slist[i].dist=0.0;
00168         //for(i=0;i<n_slist;i++)slist[i]=(s_sorted_pocket_list *)my_malloc(sizeof(s_sorted_pocket_list));
00169         s_sorted_pocket_list *el=my_malloc(sizeof(s_sorted_pocket_list));
00170 //        pcur=pockets->first;
00171 
00172 
00173         int c=0;
00174         for(i=0;i<pockets->n_pockets-1;i++){
00175             //curMobilePocket=pcur->next;
00176             for(j=i+1;j<pockets->n_pockets;j++){
00177 
00178                 el->dist=dmat[i][j];
00179                 memcpy(&(slist[c].dist),&(dmat[i][j]),sizeof(float));
00180                 //slist[i+j-1].dist=dmat[i][j];
00181                 slist[c].pid1=i;
00182                 slist[c].pid2=j;
00183                 //if(i==173 && j==299)fprintf(stdout,"\ncurdist %d %d %f %f\n",j+i-1, n_slist,dmat[i][j], slist[j+i-1].dist);
00184                 //fprintf(stdout,"%f %d %d\n",slist[c].dist,slist[c].pid1,slist[c].pid2);
00185                 c++;
00186             }
00187   //          pcur=pcur->next;
00188         }
00189         //for(i=0;i<n_slist;i++) if(slist[i].dist<-480.0)fprintf(stdout,"%f %d %d\n",slist[i].dist,slist[i].pid1,slist[i].pid2);
00190         //fflush(stdout);
00191         qsort((void *)slist,n_slist,sizeof(s_sorted_pocket_list),comp_pocket);
00192 
00193         //for(i=0;i<n_slist;i++) fprintf(stdout,"%f %d %d\n",slist[i].dist,slist[i].pid1,slist[i].pid2);
00194         /*TODO : debug here there are still some neighbours with nan pid's after the qsort*/
00195         i=0;
00196 
00197         /*create a tmp pocket list for updating pointers*/
00198         node_pocket **pock_list=my_malloc(sizeof(node_pocket *)*pockets->n_pockets);
00199         pcur=pockets->first;
00200         i=0;
00201         /*get a list of pointers to nodes*/
00202         while(pcur) {
00203             pock_list[i]=pcur;
00204             pcur=pcur->next;
00205             i++;
00206         }
00207 
00208         int idx1,idx2;
00209         node_pocket *p1,*p2;
00210         i=0;
00211         size_t init_n_pockets=pockets->n_pockets;
00212 
00213         while((slist[i].dist<=-params->sl_clust_min_nneigh) && (i<n_slist)){
00214             /*for all nearby pockets merge*/
00215             //fprintf(stdout,"%f\n",slist[i].dist);
00216             idx1=slist[i].pid1;
00217             idx2=slist[i].pid2;
00218 
00219             /*fprintf(stdout,"%d %d\n",idx1,idx2);
00220             fflush(stdout);*/
00221             if(pock_list[idx1]!=pock_list[idx2]){
00222                 p1=*(pock_list+idx1);
00223                 p2=*(pock_list+idx2);
00224                 dens1=((isnan(p1->pocket->pdesc->as_density)) ? 0.0 : p1->pocket->pdesc->as_density)/p1->pocket->pdesc->nb_asph;
00225                 dens2=((isnan(p2->pocket->pdesc->as_density)) ? 0.0 : p2->pocket->pdesc->as_density)/p2->pocket->pdesc->nb_asph;
00226                 //printf("%f vs %f\n",dens1/p1->pocket->pdesc->nb_asph, dens2/p1->pocket->pdesc->nb_asph);
00227                 if(dens1 < 0.1 && dens2 < 0.1){
00228                     for(j=0;j<init_n_pockets;j++){
00229                         if(pock_list[j]==p2){ //j!=idx2 &&
00230                             pock_list[j]=p1;
00231                             //fprintf(stdout,"update %d to %d\n",j, idx1);
00232                         }
00233                     }
00234 
00235                     mergePockets(p1,p2,pockets);
00236                     pock_list[idx2]=p1;
00237                 }
00238             }
00239             i++;
00240         }
00241 
00242         //for(i=0;i<pockets->n_pockets-1;i++){
00243             //fprintf(stdout,"dist %d vs %d = %f\n",slist[i].p1->pocket->rank,slist[i].p2->pocket->rank,slist[i].dist);
00244         //}
00245 
00246         /*free dmat*/
00247         for(i=0;i<pockets->n_pockets;i++) free(dmat[i]);
00248         free(dmat);
00249 
00250 }
00251 
00252 
00253 
00254 
00255 int comp_pocket(const void *el1, const void *el2)
00256 {
00257 
00258     s_sorted_pocket_list *ia = (s_sorted_pocket_list *)el1;
00259     s_sorted_pocket_list *ib = (s_sorted_pocket_list *)el2;
00260 
00261     //if (ia->dist<0.0)printf("dist %f\n",((s_sorted_pocket_list *)el1)->dist);
00262   if (ia->dist <  ib->dist) return -1;
00263   if (ia->dist == ib->dist) return  0;
00264   if (ia->dist >  ib->dist) return  1;
00265     return 0;
00266 }
00267 
00268 
00269 
00270 
00271 /**
00272    ## FONCTION:
00273         void pck_ml_clust(c_lst_pockets *pockets, s_fparams *params)
00274 
00275    ## SPECIFICATION:
00276         This function will apply a mutliple linkage clustering algorithm on the given
00277         list of pockets. Considering two pockets, if params->ml_clust_min_nneigh
00278         alpha spheres are separated by a distance lower than params->ml_clust_max_dist,
00279         then merge the two pockets.
00280 
00281    ## PARAMETRES:
00282         @ c_lst_pockets *pockets  : The list of pockets
00283         @ s_fparams *params       : Parameters of the program, including single
00284                                                                 linkage parameters
00285 
00286    ## RETURN:
00287         void
00288 
00289 */
00290 void pck_ml_clust(c_lst_pockets *pockets, s_fparams *params)
00291 {
00292         node_pocket *pcur = NULL,
00293                                 *pnext = NULL ,
00294                                 *curMobilePocket = NULL ;
00295 
00296         node_vertice *vcur = NULL ;
00297         node_vertice *curMobileVertice = NULL ;
00298 
00299         s_vvertice *vvcur = NULL,
00300                            *mvvcur = NULL ;
00301         float vcurx,
00302                   vcury,
00303                   vcurz ;
00304 
00305         /* Flag to know if two clusters are next to each other by single linkage
00306          * clustering...or not */
00307         int nflag ;
00308 
00309         if(!pockets) {
00310                 fprintf(stderr, "! Incorrect argument during Single Linkage Clustering.\n") ;
00311                 return ;
00312         }
00313 
00314         /* Set the first pocket */
00315         pcur = pockets->first ;
00316         while(pcur) {
00317                 /* Set the second pocket */
00318                 curMobilePocket = pcur->next ;
00319                 while(curMobilePocket) {
00320                         nflag = 0 ;
00321                         /* Set the first vertice/alpha sphere center of the first pocket */
00322                         vcur = pcur->pocket->v_lst->first ;
00323                         while(vcur && nflag <= params->sl_clust_min_nneigh){
00324                                 /* Set the first vertice/alpha sphere center of the second pocket */
00325                                 curMobileVertice = curMobilePocket->pocket->v_lst->first ;
00326                                 vvcur = vcur->vertice ;
00327                                 vcurx = vvcur->x ;
00328                                 vcury = vvcur->y ;
00329                                 vcurz = vvcur->z ;
00330 
00331                                 /* Double loop for vertices -> if not near */
00332                                 while(curMobileVertice && nflag <= params->sl_clust_min_nneigh){
00333                                         mvvcur = curMobileVertice->vertice ;
00334                                         if(dist(vcurx, vcury, vcurz, mvvcur->x, mvvcur->y, mvvcur->z)
00335                                                 < params->sl_clust_max_dist) {
00336                                                                                                         /*if beneath the clustering max distance, increment the distance flag*/
00337                                                 nflag++;
00338                                         }
00339                                         curMobileVertice = curMobileVertice->next;
00340                                 }
00341                                 vcur = vcur->next ;
00342                         }
00343 
00344                         pnext =  curMobilePocket->next ;
00345                         /* If the distance flag has counted enough occurences of near neighbours, merge pockets*/
00346                         if(nflag >= params->sl_clust_min_nneigh) {
00347                                 /* If they are next to each other, merge them */
00348                                 mergePockets(pcur,curMobilePocket,pockets);
00349                         }
00350                         curMobilePocket = pnext ;
00351                 }
00352 
00353                 pcur = pcur->next ;
00354         }
00355 }
00356 
00357 /**
00358    ## FONCTION:
00359         void pck_ml_clust(c_lst_pockets *pockets, s_fparams *params)
00360 
00361    ## SPECIFICATION:
00362         This function will apply a mutliple linkage clustering algorithm on the given
00363         list of pockets. Considering two pockets, if params->ml_clust_min_nneigh
00364         alpha spheres are separated by a distance lower than params->ml_clust_max_dist,
00365         then merge the two pockets.
00366 
00367    ## PARAMETRES:
00368         @ c_lst_pockets *pockets  : The list of pockets
00369         @ s_fparams *params       : Parameters of the program, including single
00370                                                                 linkage parameters
00371 
00372    ## RETURN:
00373         void
00374 
00375 */
00376 void pck_ml_clust_test(c_lst_pockets *pockets, s_fparams *params)
00377 {
00378         node_pocket *pcur = NULL,
00379                                 *pnext = NULL ,
00380                                 *curMobilePocket = NULL ;
00381 
00382         node_vertice *vcur = NULL ;
00383         node_vertice *curMobileVertice = NULL ;
00384 
00385         s_vvertice *vvcur = NULL,
00386                            *mvvcur = NULL ;
00387         float vcurx,
00388                   vcury,
00389                   vcurz ;
00390 
00391         /* Flag to know if two clusters are next to each other by single linkage
00392          * clustering...or not */
00393         int nflag,
00394                 restart = 0;
00395 
00396         if(!pockets) {
00397                 fprintf(stderr, "! Incorrect argument during Single Linkage Clustering.\n") ;
00398                 return ;
00399         }
00400         printf("ML starting\n") ;
00401         /* Set the first pocket */
00402         pcur = pockets->first ;
00403         while(pcur) {
00404                 /* Set the second pocket */
00405                 curMobilePocket = pcur->next ;
00406                 while(curMobilePocket) {
00407                         nflag = 0 ;
00408                         /* Set the first vertice/alpha sphere center of the first pocket */
00409                         vcur = pcur->pocket->v_lst->first ;
00410                         while(vcur && nflag <= params->sl_clust_min_nneigh){
00411                                 /* Set the first vertice/alpha sphere center of the second pocket */
00412                                 curMobileVertice = curMobilePocket->pocket->v_lst->first ;
00413                                 vvcur = vcur->vertice ;
00414                                 vcurx = vvcur->x ;
00415                                 vcury = vvcur->y ;
00416                                 vcurz = vvcur->z ;
00417 
00418                                 /* Double loop for vertices -> if not near */
00419                                 while(curMobileVertice && nflag <= params->sl_clust_min_nneigh){
00420                                         mvvcur = curMobileVertice->vertice ;
00421                                         if(dist(vcurx, vcury, vcurz, mvvcur->x, mvvcur->y, mvvcur->z)
00422                                                 < params->sl_clust_max_dist) {
00423                                                                                                         /*if beneath the clustering max distance, increment the distance flag*/
00424                                                 nflag++;
00425                                                 break ; /* Ensure that at least N vertice in one of the
00426                                                                    two pockets have N vertice at a distance
00427                                                                    <= SL_MAX_DST */
00428                                         }
00429                                         curMobileVertice = curMobileVertice->next;
00430                                 }
00431                                 vcur = vcur->next ;
00432                         }
00433 
00434                         pnext =  curMobilePocket->next ;
00435                         /* If the distance flag has counted enough occurences of near neighbours, merge pockets*/
00436                         if(nflag >= params->sl_clust_min_nneigh) {
00437                                 /* If they are next to each other, merge them */
00438                                 mergePockets(pcur,curMobilePocket,pockets);
00439                                 restart = 1 ; printf("Merging\n") ;
00440                                 break ;
00441                         }
00442                         curMobilePocket = pnext ;
00443                 }
00444 
00445                 /* Restart the algorithm if two pockets have been merged */
00446                 if(restart == 1) {
00447                         restart = 0 ;
00448                         pcur = pockets->first ;
00449                 }
00450                 else pcur = pcur->next ;
00451         }
00452         printf("ML ending\n") ;
00453 }
00454 
00455