in contrib/champ/champ.c [2134:2769]
int ChampMatch2(CChamp *I,int template,int target,
int start_tmpl,int start_targ,int n_wanted,
int *match_start,int tag_mode)
{ /* does the template covalent tree match the target? */
/* basic algorithm is a multilayer mark-and-sweep traversal
through all atoms and bonds starting from the
two nucleating atoms, which are assumed to be equivalent */
int n_match = 0;
int stereo_template = false;
#ifdef MATCHDEBUG
printf("\n\n ChampMatch2: %d %d\n",start_tmpl,start_targ);
#endif
{
/* first, initialize the marks */
ListAtom *at;
ListBond *bd;
int cur_atom,cur_bond;
/* template uses mark_tmpl */
cur_atom = I->Pat[template].atom;
while(cur_atom) {
at = I->Atom+cur_atom;
at->mark_tmpl=0;
if(at->stereo) stereo_template = true;
cur_atom = at->link;
}
cur_bond = I->Pat[template].bond;
while(cur_bond) {
bd = I->Bond+cur_bond;
bd->mark_tmpl=0;
cur_bond = bd->link;
}
/* target uses mark_targ */
cur_atom = I->Pat[target].atom;
while(cur_atom) {
at = I->Atom+cur_atom;
at->mark_targ=0;
cur_atom = at->link;
}
cur_bond = I->Pat[target].bond;
while(cur_bond) {
bd = I->Bond+cur_bond;
bd->mark_targ=0;
cur_bond = bd->link;
}
}
/* okay, now start our sweep... */
{
int tmpl_stack=0;
int targ_stack=0;
int tmpl_par;
int tmpl_idx,targ_idx;
int bond_off;
int bond_idx,atom_idx;
int mode=0; /* 0: iterate template atom, 1: need target atom */
int done_flag=false;
int match_flag;
int bond_start;
int atom_start;
int stereo_match;
ListTmpl *tmpl_ent,*parent_ent;
ListTarg *targ_ent;
ListAtom *tmpl_at,*targ_at,*base_at;
ListMatch *match_ent;
ListInt2 *int2;
/* stack entries contain:
0: the respective bond & atom indices
1: current bond index in atom
2: backward entry index (branch point)
*/
/* initialize target stack */
targ_stack = ListElemPush(&I->Targ,targ_stack);
targ_ent = I->Targ + targ_stack;
targ_ent->atom = start_targ;
targ_ent->bond = 0; /* to get here */
/* initalize template stack */
tmpl_stack = ListElemPush(&I->Tmpl,tmpl_stack);
tmpl_ent = I->Tmpl + tmpl_stack;
tmpl_ent->atom = start_tmpl;
tmpl_ent->parent = 0;
tmpl_ent->bond = 0; /* bond index to get here */
tmpl_ent->match = targ_stack; /* matching target index */
tmpl_ent->targ_start = 0; /* for matches */
/* initialize marks */
I->Atom[start_tmpl].mark_tmpl=1;
I->Atom[start_tmpl].first_tmpl = tmpl_stack;
I->Atom[start_targ].mark_targ=1;
I->Atom[start_targ].first_targ = targ_stack;
while(!done_flag) {
if(!targ_stack) break;
#ifdef MATCHDEBUG
printf("============================\n");
printf("tmpl: ");
tmpl_ent = I->Tmpl + tmpl_stack;
while(1) {
ChampAtomDump(I,tmpl_ent->atom);
printf("(%2d) ",tmpl_ent->atom);
if(tmpl_ent->link)
tmpl_ent=I->Tmpl + tmpl_ent->link;
else
break;
}
printf("\ntarg: ");
targ_ent = I->Targ + targ_stack;
while(1) {
ChampAtomDump(I,targ_ent->atom);
printf("(%2d) ",targ_ent->atom);
if(targ_ent->link)
targ_ent=I->Targ + targ_ent->link;
else
break;
}
printf("\n");
#endif
switch(mode) {
case 0: /* iterate template atom */
/* first, see if there is an open bond somewhere in the tree... */
tmpl_par = tmpl_stack; /* start from last entry */
while(tmpl_par) {
tmpl_ent = I->Tmpl + tmpl_par;
tmpl_at = I->Atom + tmpl_ent->atom; /* get atom record */
bond_off = 0;
bond_idx = tmpl_at->bond[bond_off];
while(bond_idx) { /* iterate over all bonds */
if(I->Bond[bond_idx].mark_tmpl) { /* skip over marked bonds...*/
bond_off++;
bond_idx = tmpl_at->bond[bond_off];
} else
break; /* found an open bond */
}
if(bond_idx) { /* there is an open bond */
#ifdef MATCHDEBUG
printf(" tmpl: bond_idx %d is open (%2d,%2d)\n",bond_idx,
I->Bond[bond_idx].atom[0],I->Bond[bond_idx].atom[1]);
#endif
atom_idx = I->Bond[bond_idx].atom[1];
if(atom_idx==tmpl_ent->atom)
atom_idx = I->Bond[bond_idx].atom[0];
I->Bond[bond_idx].mark_tmpl++; /* mark bond "in use" */
tmpl_stack = ListElemPush(&I->Tmpl,tmpl_stack); /* allocate new record for atom */
tmpl_ent = I->Tmpl + tmpl_stack;
tmpl_ent->atom = atom_idx; /* record which atom */
tmpl_ent->parent = tmpl_par; /* record parent identity */
tmpl_ent->bond = bond_idx; /* record bond index */
tmpl_ent->targ_start = 0; /* launch searches from this bond in parent */
tmpl_ent->match = 0; /* no match yet */
tmpl_at = I->Atom + atom_idx;
tmpl_at->mark_tmpl++; /* mark atom "in use" */
if(tmpl_at->mark_tmpl==1) { /* record where it was used */
tmpl_at->first_tmpl=tmpl_stack;
}
#ifdef MATCHDEBUG
printf(" tmpl: created tmpl_idx %d tmpl_par %d\n",
tmpl_stack,tmpl_par);
printf(" targ: tmpl atom %d:",tmpl_ent->atom);
ChampAtomDump(I,tmpl_ent->atom);
printf("\n");
#endif
mode = 1; /* now we need to find matching template bond/atom */
break; /* ...leave loop to hunt */
} else {
#ifdef MATCHDEBUG
printf(" tmpl: nothing found in level %d...\n",tmpl_par);
#endif
tmpl_par = tmpl_ent->parent; /* otherwise, proceed up tree looking for
nearest open bond */
}
}
if(!tmpl_par) { /* no open bonds-> complete match */
#ifdef MATCHDEBUG2
printf(" tmpl: EXACT MATCH DETECTED\n");
printf(" %d %d %d %d\n",template,target,start_tmpl,start_targ);
printf("tmpl: ");
tmpl_ent = I->Tmpl + tmpl_stack;
while(1) {
ChampAtomDump(I,tmpl_ent->atom);
printf("(%2d) ",tmpl_ent->atom);
if(tmpl_ent->link)
tmpl_ent=I->Tmpl + tmpl_ent->link;
else
break;
}
printf("\ntarg: ");
targ_ent = I->Targ + targ_stack;
while(1) {
ChampAtomDump(I,targ_ent->atom);
printf("(%2d) ",targ_ent->atom);
if(targ_ent->link)
targ_ent=I->Targ + targ_ent->link;
else
break;
}
printf("\n");
#endif
stereo_match = true;
if(stereo_template) { /* must check stereochemistry */
int tmpl_idx2;
int targ_idx2;
int bond_idx2;
int bond_off2;
int n_pri;
int tmpl_pri[MAX_BOND];
int targ_pri[MAX_BOND];
ListTmpl *tmpl_ent2;
ListTarg *targ_ent2;
ListBond *tmpl_bd,*targ_bd;
ListAtom *tmpl_at1,*targ_at1;
tmpl_idx = tmpl_stack; /* prepare to traverse... */
while(tmpl_idx) {
tmpl_ent = I->Tmpl + tmpl_idx;
I->Atom[tmpl_ent->atom].mark_read = false;
tmpl_idx = tmpl_ent->link;
}
tmpl_idx = tmpl_stack;
targ_idx = targ_stack;
while(tmpl_idx&&targ_idx&&stereo_match) {
tmpl_ent = I->Tmpl + tmpl_idx;
targ_ent = I->Targ + targ_idx;
tmpl_at1 = I->Atom + tmpl_ent->atom;
targ_at1 = I->Atom + targ_ent->atom;
if(!tmpl_at1->mark_read) { /* only compare non-virtual atoms */
tmpl_at1->mark_read = true;
if(tmpl_at1->stereo) {
if(!targ_at1->stereo) {
/*
printf("failure 1 %d %d %d %d \n",tmpl_at1->stereo,targ_at1->stereo,
tmpl_ent->atom,targ_ent->atom);*/
stereo_match = false;
}
/* this is a stereo-center -- so locate the corresponding bond priorities in
template and the target */
n_pri = 0;
if(tmpl_at1->imp_hydro==1) { /* deal with implicit hydrogen case */
tmpl_pri[n_pri]=0; /* considered first, since atom priorities always > 0 */
if(targ_at1->imp_hydro==1)
targ_pri[n_pri]=0;
else if(targ_at1->tot_hydro==1) {
/* track down the wayward hydrogen priority */
bond_off2 = 0;
bond_idx2 = targ_at1->bond[bond_off2];
while(bond_idx2) { /* iterate over all bonds */
targ_bd = I->Bond + bond_idx2;
if(I->Atom[targ_bd->atom[0]].atom&cH_H) {
targ_pri[n_pri] = targ_bd->pri[0];
break;
} else if(I->Atom[targ_bd->atom[1]].atom&cH_H) {
targ_pri[n_pri] = targ_bd->pri[1];
break;
}
bond_off2++;
bond_idx2 = targ_at1->bond[bond_off2];
if(!bond_idx2)
err_message("ChampMatch2","can't locate explicit hydrogen on target!");
}
/* failure? */
} else {
printf("failure 2\n");
stereo_match = false;
}
n_pri++;
}
if(!stereo_match)
break;
tmpl_idx2 = tmpl_stack; /* prepare to traverse bonds... */
while(tmpl_idx2) {
tmpl_ent2 = I->Tmpl + tmpl_idx2;
if(tmpl_ent2->bond)
I->Bond[tmpl_ent2->bond].mark_read = false;
tmpl_idx2 = tmpl_ent2->link;
}
tmpl_idx2 = tmpl_stack;
targ_idx2 = targ_stack;
while(tmpl_idx2&&targ_idx2) {
tmpl_ent2 = I->Tmpl + tmpl_idx2;
targ_ent2 = I->Targ + targ_idx2;
if(tmpl_ent2->bond&&targ_ent2->bond) {
tmpl_bd = I->Bond + tmpl_ent2->bond;
if(!tmpl_bd->mark_read) {
tmpl_bd->mark_read = true;
targ_bd = I->Bond + targ_ent2->bond;
/* if matched bond touches this atom, then
get bond priority for the atom */
if(tmpl_bd->atom[0]==tmpl_ent->atom) {
tmpl_pri[n_pri]=tmpl_bd->pri[0];
if(targ_bd->atom[0]==targ_ent->atom) {
targ_pri[n_pri]=targ_bd->pri[0];
} else { /* assuming targ_bd->atom[1]==targ_ent->atom */
targ_pri[n_pri]=targ_bd->pri[1];
}
if(n_pri<4)
n_pri++;
else {
err_message("ChampMatch2","too many connections on chiral atom!");
}
} else if(tmpl_bd->atom[1]==tmpl_ent->atom) {
tmpl_pri[n_pri]=tmpl_bd->pri[1];
if(targ_bd->atom[0]==targ_ent->atom) {
targ_pri[n_pri]=targ_bd->pri[0];
} else { /* assuming targ_bd->atom[1]==targ_ent->atom */
targ_pri[n_pri]=targ_bd->pri[1];
}
if(n_pri<4)
n_pri++;
else {
err_message("ChampMatch2","too many connections on chiral atom!");
}
}
}
}
tmpl_idx2 = tmpl_ent2->link;
targ_idx2 = targ_ent2->link;
}
if(n_pri<4) {
err_message("ChampMatch2","stereo comparison on achiral atom.");
stereo_match = false;
} else {
/* now test for a handedness match */
int tmpl_hand,targ_hand;
tmpl_hand = ChiralHandedness(tmpl_pri);
targ_hand = ChiralHandedness(targ_pri);
if(tmpl_at1->stereo==targ_at1->stereo) { /* same clockwizedness */
if(tmpl_hand!=targ_hand) { /* opposite handedness */
stereo_match=false;
}
} else { /* opposite clockwizedness */
if(tmpl_hand==targ_hand) { /* same handedness */
stereo_match = false;
}
}
}
}
}
tmpl_idx = tmpl_ent->link;
targ_idx = targ_ent->link;
}
}
if(stereo_match) {
n_match++;
if(n_wanted>0)
if(n_match>=n_wanted) done_flag=true;
if(match_start) {
(*match_start) = ListElemPush(&I->Match,*match_start);
match_ent = I->Match + (*match_start);
atom_start=0;
bond_start=0;
tmpl_idx = tmpl_stack; /* prepare to traverse... */
while(tmpl_idx) {
tmpl_ent = I->Tmpl + tmpl_idx;
I->Atom[tmpl_ent->atom].mark_read = false;
tmpl_idx = tmpl_ent->link;
}
tmpl_idx = tmpl_stack;
targ_idx = targ_stack;
while(tmpl_idx&&targ_idx) {
tmpl_ent = I->Tmpl + tmpl_idx;
targ_ent = I->Targ + targ_idx;
if(!I->Atom[tmpl_ent->atom].mark_read) { /* save non-virtual atom match */
I->Atom[tmpl_ent->atom].mark_read = true;
atom_start = ListElemPush(&I->Int2,atom_start);
int2 = I->Int2 + atom_start;
int2->value[0] = tmpl_ent->atom;
int2->value[1] = targ_ent->atom;
}
if(tmpl_ent->bond) { /* record bond match */
bond_start = ListElemPush(&I->Int2,bond_start);
int2 = I->Int2 + bond_start;
int2->value[0] = tmpl_ent->bond;
int2->value[1] = targ_ent->bond;
}
tmpl_idx = tmpl_ent->link;
targ_idx = targ_ent->link;
}
match_ent->atom = atom_start;
match_ent->bond = bond_start;
#ifdef MATCHDEBUG2
/* ChampPatDump(I,template);
ChampPatDump(I,target);*/
ChampMatchDump(I,*match_start);
#endif
}
if(tag_mode) { /* are we using tags to mark atoms and bonds? */
tmpl_idx = tmpl_stack; /* prepare to read... */
while(tmpl_idx) {
tmpl_ent = I->Tmpl + tmpl_idx;
I->Atom[tmpl_ent->atom].mark_read = false;
tmpl_idx = tmpl_ent->link;
}
tmpl_idx = tmpl_stack;
targ_idx = targ_stack;
while(tmpl_idx&&targ_idx) {
tmpl_ent = I->Tmpl + tmpl_idx;
targ_ent = I->Targ + targ_idx;
if(!I->Atom[tmpl_ent->atom].mark_read) { /* save non-virtual atom match */
I->Atom[tmpl_ent->atom].mark_read = true;
if(tag_mode==cTag_merge) { /* merge */
I->Atom[targ_ent->atom].tag |= I->Atom[tmpl_ent->atom].tag;
I->Atom[targ_ent->atom].tag &= (0xFFFFFFFF^I->Atom[tmpl_ent->atom].not_tag);
} else if(tag_mode==cTag_copy) { /* copy */
I->Atom[targ_ent->atom].tag = I->Atom[tmpl_ent->atom].tag;
}
}
if(tmpl_ent->bond) { /* record bond match */
if(tag_mode==cTag_merge) {
I->Bond[targ_ent->bond].tag |= I->Bond[tmpl_ent->bond].tag;
I->Bond[targ_ent->bond].tag &= (0xFFFFFFFF^I->Bond[tmpl_ent->bond].not_tag);
} else if(tag_mode==cTag_copy) {
I->Bond[targ_ent->bond].tag = I->Bond[tmpl_ent->bond].tag;
}
}
tmpl_idx = tmpl_ent->link;
targ_idx = targ_ent->link;
}
}
}
/* back-out the last target match... */
targ_ent = I->Targ + targ_stack;
atom_idx = targ_ent->atom;
I->Atom[atom_idx].mark_targ--; /* free-up atom */
bond_idx = targ_ent->bond;
I->Bond[bond_idx].mark_targ--; /* free-up bond */
targ_stack = ListElemPop(I->Targ,targ_stack);
mode = 1; /* now we need to find another matching template bond/atom */
}
break;
case 1:
/* try to locate matching target atom & bond */
tmpl_ent = I->Tmpl + tmpl_stack;
bond_off = tmpl_ent->targ_start; /* start with this bond */
parent_ent = I->Tmpl + tmpl_ent->parent; /* get parent of current template atom */
targ_ent = I->Targ + parent_ent->match; /* start from target atom which
matches the parent in template */
base_at = I->Atom + targ_ent->atom;
#ifdef MATCHDEBUG
printf(" targ: tmpl_stack %d bond_off %d\n",tmpl_stack,bond_off);
printf(" targ: match %d targ root atom %d:",parent_ent->match,targ_ent->atom);
ChampAtomDump(I,targ_ent->atom);
printf("\n");
printf(" targ: bonds: ");
bond_start = bond_off;
bond_idx = base_at->bond[bond_start];
while(bond_idx) { /* iterate over all bonds */
bond_start++;
printf("%d ",bond_idx);
bond_idx = base_at->bond[bond_start];
}
printf("\n");
#endif
bond_idx = base_at->bond[bond_off];
while(bond_idx) { /* iterate over all bonds */
#ifdef MATCHDEBUG
printf(" targ: trying %d bond_idx %d (%2d-%2d)\n",bond_off,bond_idx,
I->Bond[bond_idx].atom[0],I->Bond[bond_idx].atom[1]);
#endif
match_flag=false;
tmpl_ent->targ_start = bond_off + 1; /* insure we never duplicate attempt */
if(!I->Bond[bond_idx].mark_targ) { /* is bond unmarked...*/
atom_idx = I->Bond[bond_idx].atom[1]; /* find other atom */
if(atom_idx==targ_ent->atom)
atom_idx = I->Bond[bond_idx].atom[0];
tmpl_at = I->Atom+tmpl_ent->atom;
targ_at = I->Atom+atom_idx;
if((tmpl_at->mark_tmpl-targ_at->mark_targ)==1) {/* is atom available? */
if(ChampBondMatch(I->Bond+tmpl_ent->bond,
I->Bond+bond_idx)) /* does bond match? */ {
if(ChampAtomMatch(tmpl_at,targ_at)) /* does atom match? */
{
match_flag=true;
if(targ_at->mark_targ==1) { /* virtual */
/* verify that matching virtual atoms correspond to matching real atoms */
#ifdef MATCHDEBUG
printf(" targ: first_tmpl %d first_targ %d\n",
tmpl_at->first_tmpl,targ_at->first_targ);
#endif
if(I->Tmpl[tmpl_at->first_tmpl].match!=targ_at->first_targ)
match_flag=false;
}
if(match_flag) {
I->Bond[bond_idx].mark_targ++; /* mark bond "in use" */
targ_stack = ListElemPush(&I->Targ,targ_stack); /* allocate new record for atom */
targ_ent = I->Targ + targ_stack;
targ_ent->atom = atom_idx;
targ_ent->bond = bond_idx;
targ_at->mark_targ++; /* mark atom "in use" */
if(targ_at->mark_targ==1) { /* record where it was used */
targ_at->first_targ = targ_stack;
}
/* inform template entry about match */
tmpl_ent->match = targ_stack;
match_flag=true;
mode = 0; /* we have a candidate atom/bond,
so return to template */
#ifdef MATCHDEBUG
printf(" targ: MATCHED atoms %d & %d (bond %d):",
tmpl_ent->atom,
targ_ent->atom,
bond_idx);
ChampAtomDump(I,tmpl_ent->atom);
ChampAtomDump(I,targ_ent->atom);
printf("\n");
#endif
break;
}
} else {
#ifdef MATCHDEBUG
printf(" atom match failed %d vs %d ",tmpl_ent->atom,atom_idx);
ChampAtomDump(I,tmpl_ent->atom);
ChampAtomDump(I,atom_idx);
printf("\n");
#endif
}
} else {
#ifdef MATCHDEBUG
printf(" bond match failed\n");
#endif
}
}
}
if(!match_flag) {
bond_off++;
bond_idx = base_at->bond[bond_off];
}
}
if(mode) { /* unable to locate a match */
/* so back-off the previous template atom... */
tmpl_ent = I->Tmpl + tmpl_stack;
atom_idx = tmpl_ent->atom;
I->Atom[atom_idx].mark_tmpl--; /* free-up atom */
bond_idx = tmpl_ent->bond;
I->Bond[bond_idx].mark_tmpl--; /* free-up bond */
tmpl_stack = ListElemPop(I->Tmpl,tmpl_stack);
/* and back-off the previous target match */
targ_ent = I->Targ + targ_stack;
atom_idx = targ_ent->atom;
I->Atom[atom_idx].mark_targ--; /* free-up atom */
bond_idx = targ_ent->bond;
I->Bond[bond_idx].mark_targ--; /* free-up bond */
targ_stack = ListElemPop(I->Targ,targ_stack);
}
}
}
ListElemFreeChain(I->Tmpl,tmpl_stack);
ListElemFreeChain(I->Targ,targ_stack);
}
#ifdef MATCHDEBUG
printf(" ChampMatch2: returning n_match = %d\n",
n_match);
#endif
return n_match;
}