1 // $Id: ssm_graph.cpp $
2 // =================================================================
3 //
4 // CCP4 SSM Library: Protein structure superposition based on
5 // SSM algorithm: E. Krissinel & K. Henrick (2004) Acta Cryst. D60,
6 // 2256-2268.
7 //
8 // Copyright (C) Eugene Krissinel 2002-2013.
9 //
10 // This library is free software: you can redistribute it and/or
11 // modify it under the terms of the GNU Lesser General Public
12 // License version 3, modified in accordance with the provisions
13 // of the license to address the requirements of UK law.
14 //
15 // You should have received a copy of the modified GNU Lesser
16 // General Public License along with this library. If not, copies
17 // may be downloaded from http://www.ccp4.ac.uk/ccp4license.php
18 //
19 // This program is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU Lesser General Public License for more details.
23 //
24 // =================================================================
25 //
26 // 03.02.14 <-- Date of Last Modification.
27 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
28 // -----------------------------------------------------------------
29 //
30 // **** Module : ssm_graph <implementation>
31 // ~~~~~~~~~
32 // **** Classes : ssm::Graph ( secondary structure graph )
33 // ~~~~~~~~~
34 //
35 // E. Krissinel 2002-2014
36 //
37 // =================================================================
38 //
39
40
41 #include <string.h>
42
43 #include "ssm_graph.h"
44 #include "mmdb2/mmdb_math_linalg.h"
45
46
47 // ========================== ssm::Graph ===========================
48
Graph()49 ssm::Graph::Graph() : mmdb::io::Stream() {
50 InitGraph();
51 }
52
Graph(mmdb::io::RPStream Object)53 ssm::Graph::Graph ( mmdb::io::RPStream Object )
54 : mmdb::io::Stream(Object) {
55 InitGraph();
56 }
57
~Graph()58 ssm::Graph::~Graph() {
59 FreeMemory();
60 }
61
InitGraph()62 void ssm::Graph::InitGraph() {
63
64 name = NULL;
65 mmdb::CreateCopy ( name,"" );
66
67 V = NULL;
68 E = NULL;
69 graph = NULL;
70 nVertices = 0;
71 nEdges = 0;
72 nVAlloc = 0;
73 nEAlloc = 0;
74 nGAlloc = 0;
75
76 nHelices = 0;
77 nStrands = 0;
78
79 strcpy ( devChain," " );
80 devChain[0] = char(1);
81
82 }
83
FreeMemory()84 void ssm::Graph::FreeMemory() {
85 int i;
86
87 if (name) {
88 delete[] name;
89 name = NULL;
90 }
91
92 if (V) {
93 for (i=0;i<nVAlloc;i++)
94 if (V[i]) delete V[i];
95 delete[] V;
96 V = NULL;
97 }
98 nVertices = 0;
99 nVAlloc = 0;
100
101 ReleaseEdges();
102
103 }
104
SetGraphName(mmdb::cpstr gname)105 void ssm::Graph::SetGraphName ( mmdb::cpstr gname ) {
106 mmdb::CreateCopy ( name,gname );
107 }
108
109
110 namespace ssm {
111
SelectDomain(mmdb::PManager MMDB,int & selHnd,mmdb::cpstr select,mmdb::SELECTION_TYPE selType)112 int SelectDomain ( mmdb::PManager MMDB, int & selHnd,
113 mmdb::cpstr select,
114 mmdb::SELECTION_TYPE selType ) {
115 // select is of the following format:
116 // "*", "(all)" - take all file
117 // "-" - take chain without chain ID
118 // "a:Ni-Mj,b:Kp-Lq,..." - take chain a residue number N
119 // insertion code i to residue number M
120 // insertion code j plus chain b
121 // residue number K insertion code p to
122 // residue number L insertion code q and
123 // so on.
124 // "a:,b:..." - take whole chains a and b and so on
125 // "a:,b:Kp-Lq,..." - any combination of the above.
126 int rc;
127
128 selHnd = MMDB->NewSelection();
129 rc = MMDB->SelectDomain ( selHnd,select,selType,mmdb::SKEY_NEW,1 );
130 if ((!rc) && (selType==mmdb::STYPE_ATOM)) {
131 // only C-alphas are needed
132 MMDB->Select ( selHnd,selType,MMDB->GetFirstModelNum(),
133 "*",mmdb::ANY_RES,"*",mmdb::ANY_RES,"*",
134 "*","[ CA ]","*","*",mmdb::SKEY_AND );
135 }
136
137 return rc;
138
139 }
140
CutOutDomain(mmdb::PManager MMDB,mmdb::cpstr select)141 int CutOutDomain ( mmdb::PManager MMDB, mmdb::cpstr select ) {
142 // select is of the following format:
143 // "*", "(all)" - take all file
144 // "-" - chain without chain ID
145 // "a:Ni-Mj,b:Kp-Lq,..." - take chain a residue number N
146 // insertion code i to residue number M
147 // insertion code j plus chain b
148 // residue number K insertion code p to
149 // residue number L insertion code q and
150 // so on.
151 // "a:,b:..." - take whole chains a and b and so on
152 // "a:,b:Kp-Lq,..." - any combination of the above.
153 int selHnd,rc;
154
155 if (!select) return 0;
156 if ((!select[0]) || (select[0]=='*')) return 0;
157 else if (!strcasecmp(select,"(all)")) return 0;
158
159 rc = SelectDomain ( MMDB,selHnd,select,mmdb::STYPE_RESIDUE );
160
161 if (!rc) {
162 MMDB->Select ( selHnd,mmdb::STYPE_RESIDUE,0,"*",
163 mmdb::ANY_RES,"*",mmdb::ANY_RES,"*",
164 "*","*","*","*",mmdb::SKEY_XOR );
165 MMDB->DeleteSelObjects ( selHnd );
166 MMDB->FinishStructEdit ();
167 MMDB->DeleteSelection ( selHnd );
168 }
169
170 return rc;
171
172 }
173
174 }
175
SelectCalphas(mmdb::PManager MMDB,int & selHnd,mmdb::cpstr selstring)176 void ssm::Graph::SelectCalphas ( mmdb::PManager MMDB, int & selHnd,
177 mmdb::cpstr selstring ) {
178 mmdb::PChainID chain;
179 mmdb::PPAtom A;
180 mmdb::ChainID chID;
181 mmdb::ResName rName;
182 mmdb::InsCode iCode;
183 mmdb::AltLoc aLoc;
184 mmdb::pstr S;
185 int i,j,k,nchains,natoms,seqNum;
186 bool B;
187
188 if (selstring) {
189
190 S = NULL;
191 mmdb::CreateCopy ( S,selstring );
192
193 } else {
194
195 chain = NULL;
196 GetAllChains ( chain,nchains );
197
198 i = (sizeof(mmdb::ChainID)+4)*nchains + 5;
199 S = new char[i];
200 S[0] = char(0);
201 for (i=0;i<nchains;i++) {
202 if (i>0) strcat ( S,"," );
203 if (chain[i][0]!=' ') strcat ( S,chain[i] );
204 }
205 if (chain) delete[] chain;
206
207 }
208
209 selHnd = MMDB->NewSelection();
210 MMDB->Select ( selHnd,mmdb::STYPE_ATOM,1,S,
211 mmdb::ANY_RES,"*",mmdb::ANY_RES,"*",
212 "*","[ CA ]","*","*",mmdb::SKEY_NEW );
213
214 // deselect extra altlocs, if any found
215 MMDB->GetSelIndex ( selHnd,A,natoms );
216
217 i = 0;
218 k = 0;
219 while (i<natoms) {
220 if (A[i]) {
221 seqNum = A[i]->GetSeqNum();
222 strcpy ( chID ,A[i]->GetChainID() );
223 strcpy ( rName,A[i]->GetResName() );
224 strcpy ( iCode,A[i]->GetInsCode() );
225 strcpy ( aLoc ,A[i]->altLoc );
226 B = true;
227 j = i;
228 i++;
229 while ((i<natoms) && B)
230 if (A[i]) {
231 if ((A[i]->GetSeqNum()==seqNum) &&
232 (!strcmp(A[i]->GetInsCode(),iCode)) &&
233 (!strcmp(A[i]->GetResName(),rName)) &&
234 (!strcmp(A[i]->GetChainID(),chID ))) {
235 if (A[i]->occupancy<A[j]->occupancy) {
236 MMDB->SelectAtom ( selHnd,A[i],mmdb::SKEY_CLR,false );
237 k++;
238 } else if (A[i]->occupancy>A[j]->occupancy) {
239 strcpy ( aLoc,A[i]->altLoc );
240 MMDB->SelectAtom ( selHnd,A[j],mmdb::SKEY_CLR,false );
241 k++;
242 j = i;
243 } else if (!aLoc[0]) {
244 MMDB->SelectAtom ( selHnd,A[i],mmdb::SKEY_CLR,false );
245 k++;
246 } else if (!A[i]->altLoc[0]) {
247 strcpy ( aLoc,A[i]->altLoc );
248 MMDB->SelectAtom ( selHnd,A[j],mmdb::SKEY_CLR,false );
249 k++;
250 j = i;
251 } else if (strcmp(aLoc,A[i]->altLoc)<=0) {
252 MMDB->SelectAtom ( selHnd,A[i],mmdb::SKEY_CLR,false );
253 k++;
254 } else {
255 strcpy ( aLoc,A[i]->altLoc );
256 MMDB->SelectAtom ( selHnd,A[j],mmdb::SKEY_CLR,false );
257 k++;
258 j = i;
259 }
260 i++;
261 } else
262 B = false;
263 } else
264 i++;
265 } else
266 i++;
267 }
268
269 if (k) MMDB->MakeSelIndex ( selHnd );
270
271 if (S) delete[] S;
272
273 }
274
inRange(mmdb::cpstr chainID,int initPos,int endPos)275 bool ssm::Graph::inRange ( mmdb::cpstr chainID,
276 int initPos, int endPos ) {
277 int i;
278 for (i=0;i<nVertices;i++)
279 if (V[i]->inRange(chainID,initPos,endPos)) return true;
280 return false;
281 }
282
GetChainList(mmdb::pstr S)283 mmdb::pstr ssm::Graph::GetChainList ( mmdb::pstr S ) {
284 // returns comma-separated list of chains in graph's vertices
285 char N[100];
286 int i;
287 if (nVertices>0) {
288 if (!V[0]->chainID[0]) strcpy ( S,"''" );
289 else strcpy ( S,V[0]->chainID );
290 strcat ( S,"," );
291 for (i=1;i<nVertices;i++) {
292 if (!V[i]->chainID[0]) strcpy ( N,"''" );
293 else strcpy ( N,V[i]->chainID );
294 strcat ( N,"," );
295 if (!strstr(S,N)) strcat ( S,N );
296 }
297 if (!strcmp(S,"'',")) S[0] = char(0);
298 else S[strlen(S)-1] = char(0);
299 } else
300 S[0] = char(0);
301 return S;
302 }
303
Reset()304 void ssm::Graph::Reset() {
305 FreeMemory();
306 }
307
AddVertex(PVertex Vx)308 void ssm::Graph::AddVertex ( PVertex Vx ) {
309 int i,nV1;
310 PPVertex V1;
311 if (nVertices>=nVAlloc) {
312 nV1 = nVertices + 20;
313 V1 = new PVertex[nV1];
314 for (i=0;i<nVAlloc;i++)
315 V1[i] = V[i];
316 for (i=nVAlloc;i<nV1;i++)
317 V1[i] = NULL;
318 if (V) delete[] V;
319 V = V1;
320 nVAlloc = nV1;
321 }
322 V[nVertices++] = Vx;
323 }
324
GetSSEType(mmdb::pstr chainID,int atomPos)325 ssm::VERTEX_TYPE ssm::Graph::GetSSEType ( mmdb::pstr chainID,
326 int atomPos ) {
327 // Returns SSE type (V_HELIX or V_STRAND) for atom at atomPos
328 // sequence position (0.. on ).
329 VERTEX_TYPE sse;
330 int i;
331 sse = V_UNKNOWN;
332 for (i=0;i<nVertices;i++)
333 if ((!strcmp(V[i]->chainID,chainID)) &&
334 (V[i]->initPos<=atomPos) &&
335 (atomPos<=V[i]->endPos)) {
336 sse = V[i]->type;
337 break;
338 }
339 return sse;
340 }
341
GetSSEType(mmdb::PAtom A)342 ssm::VERTEX_TYPE ssm::Graph::GetSSEType ( mmdb::PAtom A ) {
343 // Returns SSE type (V_HELIX or V_STRAND) for the atom
344 mmdb::pstr chID;
345 int i,apos;
346 VERTEX_TYPE sse;
347 sse = V_UNKNOWN;
348 if (A) {
349 chID = A->GetChainID();
350 if (chID) {
351 apos = A->GetResidueNo();
352 for (i=0;i<nVertices;i++)
353 if ((!strcmp(V[i]->chainID,chID)) &&
354 (V[i]->initPos<=apos) &&
355 (apos<=V[i]->endPos)) {
356 sse = V[i]->type;
357 break;
358 }
359 }
360 }
361 return sse;
362 }
363
MakeGraph(mmdb::PManager MMDB)364 int ssm::Graph::MakeGraph ( mmdb::PManager MMDB ) {
365 mmdb::PModel model;
366 mmdb::PPResidue res;
367 PVertex vertex;
368 int rc,nresidues,i,j,k;
369 VERTEX_TYPE vtype;
370
371 FreeMemory();
372
373 model = MMDB->GetModel(MMDB->GetFirstModelNum());
374 if (!model) return mmdb::SSERC_noResidues;
375
376 rc = model->CalcSecStructure ( true );
377 if (rc!=mmdb::SSERC_Ok) return rc;
378
379 res = NULL;
380 // #### was MMDB->GetResidueTable ( res,nresidues );
381 model->GetResidueTable ( res,nresidues );
382
383 i = 0;
384 k = 0;
385 while (i<nresidues) {
386 while ((i<nresidues) && (res[i]->SSE!=mmdb::SSE_Strand) &&
387 (res[i]->SSE!=mmdb::SSE_Helix)) i++;
388 if (i<nresidues) {
389 j = i;
390 while ((i<nresidues) && (res[i]->SSE==res[j]->SSE) &&
391 (!strcmp(res[i]->GetChainID(),res[j]->GetChainID())))
392 i++;
393 if (res[j]->SSE==mmdb::SSE_Strand) vtype = V_STRAND;
394 else vtype = V_HELIX;
395 k++;
396 vertex = new ssm::Vertex();
397 if (vertex->SetVertex(MMDB,vtype,k,1,res[j]->GetChainID(),
398 res[j ]->GetSeqNum(),res[j ]->GetInsCode(),
399 res[i-1]->GetSeqNum(),res[i-1]->GetInsCode())) {
400 k--;
401 delete vertex;
402 } else
403 AddVertex ( vertex );
404 }
405 }
406
407 if (res) delete[] res;
408
409 if (nVertices<=0) return SSGE_NoVertices;
410
411 RepairSS ( MMDB );
412
413 return SSGE_Ok;
414
415 }
416
417
CalcVertexOrder()418 void ssm::Graph::CalcVertexOrder() {
419 int i,i0,VNo;
420 bool Done;
421
422 for (i=0;i<nVertices;i++)
423 V[i]->VNo = 0;
424
425 do {
426 i0 = 0;
427 while (i0<nVertices)
428 if (V[i0]->VNo==0) break;
429 else i0++;
430 Done = (i0>=nVertices);
431 if (!Done) {
432 VNo = 0;
433 for (i=0;i<nVertices;i++)
434 if (!strcmp(V[i]->chainID,V[i0]->chainID)) {
435 if (V[i]->VNo>VNo) VNo = V[i]->VNo;
436 if ((V[i]->VNo==0) &&
437 (V[i]->endPos<=V[i0]->initPos))
438 i0 = i;
439 }
440 V[i0]->VNo = VNo + 1;
441 }
442 } while (!Done);
443
444 }
445
446 #define VXREP_SD 3
447 #define VXREP_HX 3
448 #define VXREP_DEV1 (mmdb::Pi/9.0)
449 #define VXREP_DEV2 (mmdb::Pi/3.0)
450
RepairSS(mmdb::PManager MMDB)451 void ssm::Graph::RepairSS ( mmdb::PManager MMDB ) {
452 // Reunites strands and helices separated by just a few
453 // residues if there are indications that that could
454 // be a single strand or helix.
455 PPVertex Vx;
456 mmdb::ChainID chID;
457 mmdb::realtype adev1,adev2, vx,vy,vz;
458 int i,j,k,j0,k0,k1,iPos;
459 bool B;
460
461 if (nVertices<=1) return;
462
463 adev1 = VXREP_DEV1;
464 adev2 = VXREP_DEV2;
465
466 Vx = new PVertex[nVAlloc];
467 k = 0;
468
469 for (i=0;i<nVertices;i++)
470 if (V[i]) {
471 k0 = k; // starting vertex for the chain
472 strcpy ( chID,V[i]->chainID );
473 // Put all vertices of this chain into order
474 // of increasing their position in the chain
475 do {
476 j0 = -1;
477 iPos = mmdb::MaxInt;
478 for (j=i;j<nVertices;j++)
479 if (V[j]) {
480 if ((!strcmp(chID,V[j]->chainID)) &&
481 (V[j]->initPos<=iPos)) {
482 iPos = V[j]->initPos;
483 j0 = j;
484 }
485 }
486 if (j0>=0) {
487 Vx[k++] = V[j0];
488 V[j0] = NULL;
489 }
490 } while (j0>=0);
491 // Check pairs of neighbouring vertices for gaps
492 // between them, and if the gaps are small while
493 // the general direction is kept, merge the
494 // vertices
495 while (k0<k-1) {
496 k1 = k0 + 1;
497 B = (Vx[k0]->type==Vx[k1]->type);
498 if (B) {
499 if (Vx[k0]->type==V_STRAND)
500 B = ((Vx[k1]->initPos-Vx[k0]->endPos)<=VXREP_SD);
501 else if (Vx[k0]->classID==Vx[k1]->classID)
502 B = ((Vx[k1]->initPos-Vx[k0]->endPos)<=VXREP_HX);
503 if (B) {
504 vx = Vx[k1]->GetX2() - Vx[k0]->GetX1();
505 vy = Vx[k1]->GetY2() - Vx[k0]->GetY1();
506 vz = Vx[k1]->GetZ2() - Vx[k0]->GetZ1();
507 if (((Vx[k0]->GetAngle(vx,vy,vz)<=adev1) ||
508 (Vx[k1]->GetAngle(vx,vy,vz)<=adev1)) &&
509 (Vx[k0]->GetAngle(Vx[k1])<=adev2)) {
510 // Evidently k0 and k1 represent a broken strand/helix.
511 // Restore it.
512 Vx[k0]->SetVertex ( MMDB,Vx[k0]->type,Vx[k0]->serNum,
513 Vx[k0]->classID,chID,Vx[k0]->initSeqNum,
514 Vx[k0]->initICode,Vx[k1]->endSeqNum,
515 Vx[k1]->endICode );
516 delete Vx[k1];
517 k--; // one vertex less in the chain
518 while (k1<k) {
519 Vx[k1] = Vx[k1+1];
520 k1++;
521 }
522 } else
523 B = false;
524 }
525 }
526 if (!B) k0++;
527 }
528 }
529
530 delete[] V;
531
532 nVertices = k;
533 while (k<nVAlloc)
534 Vx[k++] = NULL;
535
536 V = Vx;
537
538 }
539
BuildGraph()540 void ssm::Graph::BuildGraph() {
541 int i,j;
542
543 ReleaseEdges ();
544 CalcVertexOrder();
545
546 nHelices = 0;
547 nStrands = 0;
548
549 if (nVertices>1) {
550
551 nGAlloc = nVertices;
552 mmdb::GetMatrixMemory ( graph,nGAlloc,nGAlloc,1,1 );
553
554 // Connect all vertices with edges. The following
555 // is unsophisticated, but simple and reliable:
556 for (i=1;i<=nVertices;i++) {
557 V[i-1]->id = i;
558 if (V[i-1]->type==V_HELIX) nHelices++;
559 else nStrands++;
560 graph[i][i] = -1;
561 for (j=i+1;j<=nVertices;j++) {
562 graph[i][j] = nEdges++;
563 graph[j][i] = graph[i][j];
564 }
565 }
566
567 if (nEdges>0) {
568 nEAlloc = nEdges;
569 E = new PEdge[nEAlloc];
570 nEdges = 0;
571 for (i=1;i<=nVertices;i++)
572 for (j=i+1;j<=nVertices;j++) {
573 E[nEdges] = new ssm::Edge();
574 E[nEdges]->SetEdge ( V[i-1],V[j-1] );
575 nEdges++;
576 }
577 if (nEdges!=nEAlloc)
578 printf ( "\n #### PROGRAM ERROR IN ssm::Graph::BuildGraph()\n" );
579 }
580
581 }
582
583 }
584
isBuild()585 bool ssm::Graph::isBuild() {
586 if (!graph) return false;
587 if (!E) return false;
588 return true;
589 }
590
591
calcVTypes()592 void ssm::Graph::calcVTypes() {
593 int i;
594 nHelices = 0;
595 nStrands = 0;
596 for (i=0;i<nVertices;i++)
597 if (V[i]->type==V_HELIX) nHelices++;
598 else nStrands++;
599 }
600
ReleaseEdges()601 void ssm::Graph::ReleaseEdges() {
602 int i;
603
604 mmdb::FreeMatrixMemory ( graph,nGAlloc,1,1 );
605 nGAlloc = 0;
606
607 for (i=0;i<nEAlloc;i++)
608 if (E[i]) delete E[i];
609 if (E) delete[] E;
610 E = NULL;
611 nEdges = 0;
612 nEAlloc = 0;
613
614 }
615
GetEdgeDirection(int v1,int v2,mmdb::vect3 & v)616 bool ssm::Graph::GetEdgeDirection ( int v1, int v2, mmdb::vect3 & v ) {
617 if (graph) {
618 if ((1<=v1) && (v1<=nVertices) &&
619 (1<=v2) && (v2<=nVertices) && (v1!=v2)) {
620 E[graph[v1][v2]]->GetDirection ( v );
621 if (v1>v2) {
622 v[0] = -v[0];
623 v[1] = -v[1];
624 v[2] = -v[2];
625 }
626 return true;
627 }
628 }
629 return false;
630 }
631
CompareEdges(int i,int j,PGraph G,int k,int l)632 int ssm::Graph::CompareEdges ( int i, int j, PGraph G,
633 int k, int l ) {
634 // CompareEdges(..) compares edge (ij) of the graph with
635 // edge (kl) of graph G. i may be either less or greater
636 // than j, same about k and l. If edges compare, the function
637 // returns 0. Edges with equal indices (i.e. (ii) and (kk))
638 // are considered as comparable (returns 0).
639 // The function may be used only after both graphs have
640 // been built.
641 if (i==j) {
642 if (k==l) return 0;
643 else return 7;
644 } else if (k==l)
645 return 7;
646 return E[graph[i][j]]->Compare (
647 (i>j),G->E[G->graph[k][l]],(k>l) );
648 }
649
CheckEdgeConnectivity(int i,int j,PGraph G,int k,int l)650 int ssm::Graph::CheckEdgeConnectivity ( int i, int j, PGraph G,
651 int k, int l ) {
652 if (i==j) return -1;
653 if (k==l) return -1;
654 return E[graph[i][j]]->CheckConnectivity (
655 (i>j),G->E[G->graph[k][l]],(k>l) );
656 }
657
RemoveShortVertices(int nmin_hx,int nmin_sd)658 void ssm::Graph::RemoveShortVertices ( int nmin_hx, int nmin_sd ) {
659 PPVertex Vx;
660 int i,n;
661 n = 0;
662 for (i=0;i<nVertices;i++)
663 if (V[i]) {
664 if (((V[i]->type==V_HELIX) && (V[i]->nres>nmin_hx)) ||
665 ((V[i]->type==V_STRAND) && (V[i]->nres>nmin_sd)))
666 n++;
667 }
668 if (n<nVertices) {
669 if (n>0) {
670 Vx = new PVertex[n];
671 n = 0;
672 for (i=0;i<nVertices;i++)
673 if (V[i]) {
674 if (((V[i]->type==V_HELIX) && (V[i]->nres>nmin_hx)) ||
675 ((V[i]->type==V_STRAND) && (V[i]->nres>nmin_sd)))
676 Vx[n++] = V[i];
677 else delete V[i];
678 V[i] = NULL;
679 }
680 for (i=nVertices;i<nVAlloc;i++)
681 if (V[i]) delete V[i];
682 delete[] V;
683 V = Vx;
684 nVertices = n;
685 nVAlloc = 0;
686 } else if (V) {
687 for (i=0;i<nVAlloc;i++)
688 if (V[i]) delete V[i];
689 delete[] V;
690 V = NULL;
691 nVertices = 0;
692 nVAlloc = 0;
693 }
694 }
695 BuildGraph();
696 }
697
698
LeaveVertices(mmdb::ivector vlist,int vllen)699 void ssm::Graph::LeaveVertices ( mmdb::ivector vlist, int vllen ) {
700 int i,j,n;
701 bool B;
702 n = 0;
703 for (i=0;i<nVertices;i++) {
704 B = false;
705 for (j=1;(j<=vllen) && (!B);j++)
706 B = (vlist[j]==i+1);
707 if (!B) {
708 if (V[i])
709 delete V[i];
710 V[i] = NULL;
711 } else if (n<i) {
712 V[n++] = V[i];
713 V[i] = NULL;
714 } else
715 n++;
716 }
717 nVertices = n;
718 }
719
720
LeaveVertices(mmdb::cpstr select,mmdb::PManager M)721 void ssm::Graph::LeaveVertices ( mmdb::cpstr select,
722 mmdb::PManager M ) {
723 // select is of the following format:
724 // "*", "(all)" - take all file
725 // "-" - take chain without chain ID
726 // "a:Ni-Mj,b:Kp-Lq,..." - take chain a residue number N
727 // insertion code i to residue number M
728 // insertion code j plus chain b
729 // residue number K insertion code p to
730 // residue number L insertion code q and
731 // so on.
732 // "a:,b:..." - take whole chains a and b and so on
733 // "a:,b:Kp-Lq,..." - any combination of the above.
734 int rc,selHnd1;
735 rc = ssm::SelectDomain ( M,selHnd1,select,mmdb::STYPE_RESIDUE );
736 if (!rc) _leaveVertices ( M,selHnd1 );
737 M->DeleteSelection ( selHnd1 );
738 }
739
740
LeaveVertices(int selHnd,mmdb::PManager M)741 void ssm::Graph::LeaveVertices ( int selHnd, mmdb::PManager M ) {
742 // Leaves only vertices that are covered by the given selection.
743 // selHnd may refer to the selection of atoms, residues or chains.
744 int stype,selHnd1;
745
746 stype = M->GetSelType ( selHnd );
747 if ((stype==mmdb::STYPE_INVALID) || (stype==mmdb::STYPE_UNDEFINED))
748 return;
749
750 if (stype==mmdb::STYPE_RESIDUE) selHnd1 = selHnd;
751 else {
752 selHnd1 = M->NewSelection();
753 M->Select ( selHnd1,mmdb::STYPE_RESIDUE,selHnd,mmdb::SKEY_NEW );
754 }
755
756 _leaveVertices ( M,selHnd1 );
757
758 if (stype!=mmdb::STYPE_RESIDUE) M->DeleteSelection ( selHnd1 );
759
760 }
761
762
_leaveVertices(mmdb::PManager M,int selHnd1)763 void ssm::Graph::_leaveVertices ( mmdb::PManager M, int selHnd1 ) {
764 int selHnd2, i,mdl,n;
765
766 mdl = M->GetFirstModelNum();
767 selHnd2 = M->NewSelection();
768
769 n = 0;
770 for (i=0;i<nVertices;i++)
771 if (V[i]) {
772 M->Select ( selHnd2,mmdb::STYPE_RESIDUE,mdl,V[i]->chainID,
773 V[i]->initSeqNum,V[i]->initICode,
774 V[i]->endSeqNum ,V[i]->endICode ,
775 "*","*","*","*",mmdb::SKEY_NEW );
776 M->Select ( selHnd2,mmdb::STYPE_RESIDUE,selHnd1,mmdb::SKEY_AND );
777 if (M->GetSelLength(selHnd2)<=0) {
778 delete V[i];
779 V[i] = NULL;
780 } else if (n<i) {
781 V[n++] = V[i];
782 V[i] = NULL;
783 } else
784 n++;
785 }
786
787 nVertices = n;
788
789 M->DeleteSelection ( selHnd2 );
790
791 }
792
793
RemoveVertex(int vertex_no)794 void ssm::Graph::RemoveVertex ( int vertex_no ) {
795 int i;
796 if ((0<vertex_no) && (vertex_no<=nVertices)) {
797 if (V[vertex_no-1]) delete V[vertex_no-1];
798 for (i=vertex_no;i<nVertices;i++)
799 V[i-1] = V[i];
800 V[nVertices-1] = NULL;
801 nVertices--;
802 }
803 }
804
805
GetVertexType(int vertex_no)806 ssm::VERTEX_TYPE ssm::Graph::GetVertexType ( int vertex_no ) {
807 if ((vertex_no>0) && (vertex_no<=nVertices)) {
808 if (V[vertex_no-1]) return V[vertex_no-1]->type;
809 }
810 return V_UNKNOWN;
811 }
812
GetVertexClass(int vertex_no)813 int ssm::Graph::GetVertexClass ( int vertex_no ) {
814 if ((vertex_no>0) && (vertex_no<=nVertices)) {
815 if (V[vertex_no-1]) return V[vertex_no-1]->classID;
816 }
817 return 0;
818 }
819
GetVertexDirection(int vertex_no,mmdb::vect3 & v)820 bool ssm::Graph::GetVertexDirection ( int vertex_no, mmdb::vect3 & v ) {
821 if ((vertex_no>0) && (vertex_no<=nVertices)) {
822 if (V[vertex_no-1]) {
823 V[vertex_no-1]->GetDirection ( v );
824 return true;
825 }
826 }
827 return false;
828 }
829
GetSeqLength(int vertex_no)830 int ssm::Graph::GetSeqLength ( int vertex_no ) {
831 if ((vertex_no>0) && (vertex_no<=nVertices)) {
832 if (V[vertex_no-1]) return V[vertex_no-1]->nres;
833 }
834 return V_UNKNOWN;
835 }
836
GetMass(int vertex_no)837 mmdb::realtype ssm::Graph::GetMass ( int vertex_no ) {
838 if ((vertex_no>0) && (vertex_no<=nVertices)) {
839 if (V[vertex_no-1]) return V[vertex_no-1]->mass;
840 }
841 return 0.0;
842 }
843
GetGraphVertex(int vertex_no)844 ssm::PVertex ssm::Graph::GetGraphVertex ( int vertex_no ) {
845 if ((vertex_no>0) && (vertex_no<=nVertices))
846 return V[vertex_no-1];
847 else return NULL;
848 }
849
GetVertexChainID(int vertex_no)850 mmdb::pstr ssm::Graph::GetVertexChainID ( int vertex_no ) {
851 if ((vertex_no>0) && (vertex_no<=nVertices)) {
852 if (V[vertex_no-1]) return V[vertex_no-1]->chainID;
853 }
854 return NULL;
855 }
856
GetVertexInitRes(int vertex_no)857 mmdb::pstr ssm::Graph::GetVertexInitRes ( int vertex_no ) {
858 if ((vertex_no>0) && (vertex_no<=nVertices)) {
859 if (V[vertex_no-1]) return V[vertex_no-1]->initResName;
860 }
861 return NULL;
862 }
863
GetVertexEndRes(int vertex_no)864 mmdb::pstr ssm::Graph::GetVertexEndRes ( int vertex_no ) {
865 if ((vertex_no>0) && (vertex_no<=nVertices)) {
866 if (V[vertex_no-1]) return V[vertex_no-1]->endResName;
867 }
868 return NULL;
869 }
870
GetVertexRange(int vertex_no,mmdb::ChainID chID,int & initSeqNum,mmdb::InsCode initICode,int & endSeqNum,mmdb::InsCode endICode)871 void ssm::Graph::GetVertexRange ( int vertex_no, mmdb::ChainID chID,
872 int & initSeqNum, mmdb::InsCode initICode,
873 int & endSeqNum, mmdb::InsCode endICode ) {
874 int vn;
875 if ((vertex_no>0) && (vertex_no<=nVertices)) {
876 vn = vertex_no - 1;
877 if (V[vn]) {
878 strcpy ( chID,V[vn]->chainID );
879 initSeqNum = V[vn]->initSeqNum;
880 endSeqNum = V[vn]->endSeqNum;
881 strcpy ( initICode,V[vn]->initICode );
882 strcpy ( endICode ,V[vn]->endICode );
883 return;
884 }
885 }
886 chID[0] = char(0);
887 initSeqNum = mmdb::ANY_RES;
888 endSeqNum = mmdb::ANY_RES;
889 initICode[0] = char(0);
890 endICode[0] = char(0);
891 }
892
GetVertexRange(int vertex_no,mmdb::ChainID chID,int & initPos,int & endPos)893 void ssm::Graph::GetVertexRange ( int vertex_no, mmdb::ChainID chID,
894 int & initPos, int & endPos ) {
895 int vn;
896 if ((vertex_no>0) && (vertex_no<=nVertices)) {
897 vn = vertex_no - 1;
898 if (V[vn]) {
899 strcpy ( chID,V[vn]->chainID );
900 initPos = V[vn]->initPos;
901 endPos = V[vn]->endPos;
902 return;
903 }
904 }
905 chID[0] = char(0);
906 initPos = mmdb::ANY_RES;
907 endPos = mmdb::ANY_RES;
908 }
909
GetGraphEdge(int edge_no)910 ssm::PEdge ssm::Graph::GetGraphEdge ( int edge_no ) {
911 if ((edge_no>0) && (edge_no<=nEdges))
912 return E[edge_no-1];
913 else return NULL;
914 }
915
GetGraphEdge(int v1,int v2)916 ssm::PEdge ssm::Graph::GetGraphEdge ( int v1, int v2 ) {
917 if (graph && (1<=v1) && (v1<=nVertices) &&
918 (1<=v2) && (v2<=nVertices) &&
919 (v1!=v2))
920 return E[graph[v1][v2]];
921 else return NULL;
922 }
923
924
CalcCombinations(mmdb::ivector F,int nm)925 mmdb::realtype ssm::Graph::CalcCombinations ( mmdb::ivector F, int nm ) {
926 //
927 // F contains list of nm graph vertices (e.g. those matched),
928 // the function returns the number of combinations these
929 // vertices may be picked from the graph, taking into account
930 // their type, length and sequence order.
931 //
932 mmdb::rmatrix C;
933 mmdb::realtype nCombs;
934 int i,j,k,k0;
935
936 if ((nm<=0) || (nVertices<nm)) return 1.0;
937
938 mmdb::GetMatrixMemory ( C,nm,nVertices,1,1 );
939 for (i=1;i<=nm;i++)
940 for (j=1;j<=nVertices;j++)
941 C[i][j] = 0.0;
942
943 k0 = 0;
944 for (i=1;i<=nm;i++) {
945 k = mmdb::MaxInt4;
946 for (j=1;j<=nm;j++)
947 if ((F[j]>k0) && (F[j]<k)) k = F[j];
948 if (k<mmdb::MaxInt4) {
949 k0 = k;
950 k--;
951 for (j=i;j<=nVertices-nm+i;j++)
952 if (V[k]->Compare(V[j-1])) C[i][j] = 1.0;
953 }
954 }
955
956 for (j=nVertices-1;j>=nm;j--)
957 C[nm][j] += C[nm][j+1];
958
959 for (i=nm-1;i>=1;i--)
960 for (j=nVertices-nm+i;j>=i;j--)
961 if (C[i+1][j+1]<=0.01) C[i][j] = 0.0;
962 else if (C[i][j]<=0.01) C[i][j] = C[i][j+1];
963 else C[i][j] = C[i][j+1] + C[i+1][j+1];
964
965 nCombs = C[1][1];
966 mmdb::FreeMatrixMemory ( C,nm,1,1 );
967
968 return nCombs;
969
970 }
971
972
GetAllChains(mmdb::PChainID & chain,int & nchains)973 void ssm::Graph::GetAllChains ( mmdb::PChainID & chain,
974 int & nchains ) {
975 // returns all chain IDs found in the graph's vertices
976 int i,j,k;
977 nchains = 0;
978 if (chain) {
979 delete[] chain;
980 chain = NULL;
981 }
982 if (nVertices>0) {
983 chain = new mmdb::ChainID[nVertices];
984 for (i=0;i<nVertices;i++) {
985 k = 0;
986 // is chain of this vertex already counted?
987 for (j=0;(j<nchains) && (k==0);j++)
988 if (!strcmp(chain[j],V[i]->chainID)) k = 1;
989 if (k==0) {
990 // register the chain
991 strcpy ( chain[nchains],V[i]->chainID );
992 nchains++;
993 }
994 }
995 }
996 }
997
998
GetNofChains()999 int ssm::Graph::GetNofChains() {
1000 // counts number of chains == number of single-chain graphs
1001 mmdb::PChainID chain;
1002 int nchains;
1003 chain = NULL;
1004 GetAllChains ( chain,nchains );
1005 return nchains;
1006 }
1007
DevelopChainGraphs(PPGraph & G,int & nGraphs)1008 void ssm::Graph::DevelopChainGraphs ( PPGraph & G, int & nGraphs ) {
1009 mmdb::PChainID S;
1010 int i,j,k;
1011 PVertex Vx;
1012
1013 DisposeGraphs ( G,nGraphs );
1014
1015 if (nVertices>0) {
1016
1017 // count number of chains == number of single-chain graphs
1018 S = new mmdb::ChainID[nVertices];
1019 for (i=0;i<nVertices;i++) {
1020 k = 0;
1021 // is chain of this vertex already counted?
1022 for (j=0;(j<nGraphs) && (k==0);j++)
1023 if (!strcmp(S[j],V[i]->chainID)) k = 1;
1024 if (k==0) {
1025 // register the chain
1026 strcpy ( S[nGraphs],V[i]->chainID );
1027 nGraphs++;
1028 }
1029 }
1030
1031 if (nGraphs>0) {
1032 G = new PGraph[nGraphs];
1033 for (i=0;i<nGraphs;i++) {
1034 G[i] = new ssm::Graph();
1035 mmdb::CreateCopy ( G[i]->name ,name );
1036 mmdb::CreateConcat ( G[i]->name ,":",S[i] );
1037 strcpy ( G[i]->devChain,S[i] );
1038 for (j=0;j<nVertices;j++)
1039 if (!strcmp(S[i],V[j]->chainID)) {
1040 Vx = new Vertex();
1041 Vx->Copy ( V[j] );
1042 G[i]->AddVertex ( Vx );
1043 }
1044 G[i]->BuildGraph();
1045 }
1046 }
1047
1048 delete[] S;
1049
1050 }
1051
1052 }
1053
1054
Superpose(PGraph G,mmdb::ivector F1,mmdb::ivector F2,int nMatch,mmdb::mat44 & TMatrix)1055 void ssm::Graph::Superpose ( PGraph G, mmdb::ivector F1,
1056 mmdb::ivector F2,
1057 int nMatch, mmdb::mat44 & TMatrix ) {
1058 // Returns TMatrix - a transformation matrix for G's coordinates,
1059 // such that TMatrix*{G} ~= {this}
1060 // F1 is for this graph, F2 = for G.
1061 SuperposeGraphs ( G,F2,this,F1,nMatch,TMatrix );
1062 }
1063
1064
Copy(PGraph G)1065 void ssm::Graph::Copy ( PGraph G ) {
1066 int i;
1067
1068 FreeMemory();
1069
1070 mmdb::CreateCopy ( name ,G->name );
1071 strcpy ( devChain,G->devChain );
1072
1073 nVertices = G->nVertices;
1074 if (nVertices>0) {
1075 nVAlloc = nVertices;
1076 V = new PVertex[nVertices];
1077 for (i=0;i<nVertices;i++) {
1078 V[i] = new ssm::Vertex();
1079 V[i]->Copy ( G->V[i] );
1080 }
1081 }
1082
1083 }
1084
1085
write(mmdb::io::RFile f)1086 void ssm::Graph::write ( mmdb::io::RFile f ) {
1087 int i;
1088 int Version=1;
1089
1090 f.WriteInt ( &Version );
1091 f.CreateWrite ( name );
1092 f.WriteTerLine ( devChain,false );
1093 f.WriteInt ( &nVertices );
1094 for (i=0;i<nVertices;i++)
1095 StreamWrite ( f,V[i] );
1096
1097 }
1098
read(mmdb::io::RFile f)1099 void ssm::Graph::read ( mmdb::io::RFile f ) {
1100 int i,Version;
1101
1102 FreeMemory();
1103
1104 f.ReadInt ( &Version );
1105 f.CreateRead ( name );
1106 f.ReadTerLine ( devChain,false );
1107 f.ReadInt ( &nVertices );
1108 if (nVertices>0) {
1109 nVAlloc = nVertices;
1110 V = new PVertex[nVertices];
1111 for (i=0;i<nVertices;i++) {
1112 V[i] = NULL;
1113 StreamRead ( f,V[i] );
1114 }
1115 }
1116
1117 }
1118
1119 namespace ssm {
1120 MakeStreamFunctions(Graph)
1121 }
1122
1123
1124 // ==================================================================
1125
1126
1127 namespace ssm {
1128
GetSSGraph(mmdb::PManager M,int selHnd,int & rc)1129 PGraph GetSSGraph ( mmdb::PManager M, int selHnd, int & rc ) {
1130 PGraph G;
1131
1132 G = new Graph();
1133 rc = G->MakeGraph ( M );
1134 if (!rc) {
1135 if (selHnd>0) {
1136 G->LeaveVertices ( selHnd,M );
1137 if (G->GetNofVertices()<=0) {
1138 delete G;
1139 rc = RC_NoVertices;
1140 return NULL;
1141 }
1142 }
1143 G->BuildGraph();
1144 return G;
1145 } else {
1146 rc = RC_NoGraph;
1147 if (G) delete G;
1148 return NULL;
1149 }
1150
1151 }
1152
DisposeGraphs(PPGraph & G,int & nGraphs)1153 void DisposeGraphs ( PPGraph & G, int & nGraphs ) {
1154 int i;
1155 if (G) {
1156 for (i=0;i<nGraphs;i++)
1157 if (G[i]) delete G[i];
1158 delete[] G;
1159 }
1160 G = NULL;
1161 nGraphs = 0;
1162 }
1163
1164
SuperposeGraphs(PGraph G1,mmdb::ivector F1,PGraph G2,mmdb::ivector F2,int matchlen,mmdb::mat44 & TMatrix)1165 int SuperposeGraphs ( PGraph G1, mmdb::ivector F1,
1166 PGraph G2, mmdb::ivector F2,
1167 int matchlen,
1168 mmdb::mat44 & TMatrix ) {
1169 PVertex Vx;
1170 mmdb::rmatrix A,U,V;
1171 mmdb::rvector W,RV1;
1172 mmdb::vect3 v1,v2;
1173 mmdb::realtype det,B, x01,y01,z01, x02,y02,z02, mass,mass1,mass2;
1174 int i,j,k,l, nE1,nE2;
1175
1176 nE1 = G1->GetNofEdges();
1177 if (!nE1) G1->BuildGraph();
1178
1179 nE2 = G2->GetNofEdges();
1180 if (!nE2) G2->BuildGraph();
1181
1182 mmdb::GetMatrixMemory ( A ,3,3,1,1 );
1183 mmdb::GetMatrixMemory ( U ,3,3,1,1 );
1184 mmdb::GetMatrixMemory ( V ,3,3,1,1 );
1185 mmdb::GetVectorMemory ( W ,3,1 );
1186 mmdb::GetVectorMemory ( RV1,3,1 );
1187
1188 for (j=1;j<=3;j++)
1189 for (k=1;k<=3;k++)
1190 A[j][k] = 0.0;
1191
1192 for (i=1;i<=matchlen;i++) {
1193 Vx = G1->GetGraphVertex ( F1[i] );
1194 Vx->GetDirection ( v1 );
1195 Vx = G2->GetGraphVertex ( F2[i] );
1196 Vx->GetDirection ( v2 );
1197 for (j=1;j<=3;j++)
1198 for (k=1;k<=3;k++)
1199 A[j][k] += v1[k-1]*v2[j-1];
1200 }
1201
1202 for (i=1;i<matchlen;i++)
1203 for (l=i+1;l<=matchlen;l++)
1204 if (G1->GetEdgeDirection(F1[i],F1[l],v1) &&
1205 G2->GetEdgeDirection(F2[i],F2[l],v2))
1206 for (j=1;j<=3;j++)
1207 for (k=1;k<=3;k++)
1208 A[j][k] += v1[k-1]*v2[j-1];
1209
1210 det = A[1][1]*A[2][2]*A[3][3] +
1211 A[1][2]*A[2][3]*A[3][1] +
1212 A[2][1]*A[3][2]*A[1][3] -
1213 A[1][3]*A[2][2]*A[3][1] -
1214 A[1][1]*A[2][3]*A[3][2] -
1215 A[3][3]*A[1][2]*A[2][1];
1216
1217 mmdb::math::SVD ( 3,3,3,A,U,V,W,RV1,true,true,i );
1218
1219 if (i!=0) {
1220 for (j=0;j<4;j++) {
1221 for (k=0;k<4;k++)
1222 TMatrix[j][k] = 0.0;
1223 TMatrix[j][j] = 1.0;
1224 }
1225 return 1;
1226 }
1227
1228 if (det<0.0) {
1229 k = 0;
1230 B = mmdb::MaxReal;
1231 for (j=1;j<=3;j++)
1232 if (W[j]<B) {
1233 B = W[j];
1234 k = j;
1235 }
1236 for (j=1;j<=3;j++)
1237 V[k][j] = -V[k][j];
1238 }
1239
1240 for (j=1;j<=3;j++)
1241 for (k=1;k<=3;k++) {
1242 B = 0.0;
1243 for (i=1;i<=3;i++)
1244 B += U[j][i]*V[k][i];
1245 TMatrix[j-1][k-1] = B;
1246 }
1247
1248
1249 // 9. Add translation
1250 x01 = 0.0; y01 = 0.0; z01 = 0.0; mass1 = 0.0;
1251 x02 = 0.0; y02 = 0.0; z02 = 0.0; mass2 = 0.0;
1252 for (i=1;i<=matchlen;i++) {
1253 Vx = G1->GetGraphVertex ( F1[i] );
1254 mass = Vx->GetMass();
1255 Vx->GetPosition ( v1 );
1256 x01 += v1[0]*mass;
1257 y01 += v1[1]*mass;
1258 z01 += v1[2]*mass;
1259 mass1 += mass;
1260 Vx = G2->GetGraphVertex ( F2[i] );
1261 mass = Vx->GetMass();
1262 Vx->GetPosition ( v2 );
1263 x02 += v2[0]*mass;
1264 y02 += v2[1]*mass;
1265 z02 += v2[2]*mass;
1266 mass2 += mass;
1267 }
1268 x01 /= mass1; y01 /= mass1; z01 /= mass1;
1269 x02 /= mass2; y02 /= mass2; z02 /= mass2;
1270 TMatrix[0][3] = x02 - TMatrix[0][0]*x01 - TMatrix[0][1]*y01 -
1271 TMatrix[0][2]*z01;
1272 TMatrix[1][3] = y02 - TMatrix[1][0]*x01 - TMatrix[1][1]*y01 -
1273 TMatrix[1][2]*z01;
1274 TMatrix[2][3] = z02 - TMatrix[2][0]*x01 - TMatrix[2][1]*y01 -
1275 TMatrix[2][2]*z01;
1276
1277 mmdb::FreeMatrixMemory ( A ,1,1 );
1278 mmdb::FreeMatrixMemory ( U ,1,1 );
1279 mmdb::FreeMatrixMemory ( V ,1,1 );
1280 mmdb::FreeVectorMemory ( W ,1 );
1281 mmdb::FreeVectorMemory ( RV1,1 );
1282
1283 if (!nE1) G1->ReleaseEdges();
1284 if (!nE2) G2->ReleaseEdges();
1285
1286 return 0;
1287
1288 }
1289
1290
CalcCombinations(mmdb::rvector & combs,int & vlen,PGraph G1,PGraph G2)1291 void CalcCombinations ( mmdb::rvector & combs, int & vlen,
1292 PGraph G1, PGraph G2 ) {
1293 // combs[i], i=1..vlen, returns the number of common
1294 // substructures of size i of graphs G1 and G2. The
1295 // sequential order of graph vertices is taken into
1296 // account, however the actual edges are completely
1297 // neglected.
1298 PPVertex V1,V2;
1299 mmdb::rmatrix3 P;
1300 mmdb::imatrix C;
1301 mmdb::realtype q;
1302 int n,m, i,j,k;
1303
1304 n = G1->GetNofVertices();
1305 m = G2->GetNofVertices();
1306 if (n<=m) {
1307 V1 = G1->GetVertices();
1308 V2 = G2->GetVertices();
1309 } else {
1310 m = G1->GetNofVertices();
1311 n = G2->GetNofVertices();
1312 V2 = G1->GetVertices();
1313 V1 = G2->GetVertices();
1314 }
1315
1316 vlen = 0;
1317 mmdb::FreeVectorMemory ( combs,1 );
1318 if (n<=0) return;
1319
1320 mmdb::GetMatrix3Memory ( P,n,m,n,1,1,1 );
1321 mmdb::GetMatrixMemory ( C,n,m,1,1 );
1322 for (i=1;i<=n;i++)
1323 for (j=1;j<=m;j++) {
1324 if (V1[i-1]->Compare(V2[j-1])) C[i][j] = 1;
1325 else C[i][j] = 0;
1326 for (k=1;k<=n;k++)
1327 P[i][j][k] = 0.0;
1328 }
1329
1330 q = 0.0;
1331 for (j=1;j<=m;j++) {
1332 q += C[1][j];
1333 P[1][j][1] = q;
1334 }
1335
1336 for (i=2;i<=n;i++) {
1337
1338 q = 0.0;
1339 for (j=1;j<=m;j++) {
1340 q += C[i][j];
1341 P[i][j][1] = P[i-1][j][1] + q;
1342 }
1343
1344 for (k=2;k<=i;k++) {
1345 for (j=k;j<=m;j++)
1346 if (C[i][j]==0) P[i][j][k] = P[i][j-1][k];
1347 else P[i][j][k] = P[i][j-1][k] + P[i-1][j-1][k-1];
1348 for (j=k;j<=m;j++)
1349 P[i][j][k] += P[i-1][j][k];
1350 }
1351
1352 }
1353
1354 vlen = n;
1355 mmdb::GetVectorMemory ( combs,n,1 );
1356 for (k=1;k<=n;k++)
1357 combs[k] = P[n][m][k];
1358
1359 mmdb::FreeMatrix3Memory ( P,n,m,1,1,1 );
1360 mmdb::FreeMatrixMemory ( C,n,1,1 );
1361
1362 }
1363
1364 }
1365