/*
************************************************************************
*
*   match.C -
*
*   Copyright (c) 1995
*
*   ETH Zuerich
*   Institut fuer Molekularbiologie und Biophysik
*   ETH-Hoenggerberg
*   CH-8093 Zuerich
*
*   All Rights Reserved
*
*   Date of last modification : 95/03/03
*   Pathname of SCCS file     : /export/home3/cb/garant-1.0/src/SCCS/s.atom.C
*   SCCS identification       : 1.1
*
************************************************************************
*/
/**************************************************************************/
/*                             match.C                                    */
/*                                                                        */
/* matching of two peak lists                                             */
/**************************************************************************/

#include <string.h>
#include <stdio.h>
#include <math.h>
#include "match.h"
#include "global.h"			/* global constants               */
#include "log.h"			/* general log file               */
#include "specTypes.h"
#include "peak.h"
#include "co.h"


static float delt[DIMENSION];		/* offset between the peak lists  */
static float power;			/* norm to calculate distances    */


class Node : public Element {
public:
   Peak *p;
   float s;				/* offset for score               */
   int m;				/* matching                       */
   int pred;				/* pointer to                     */
   float score;
   float calcScore(Node &);
};

float Node::calcScore(Node &n)
{
    float d, t;
    int i;
    
    d=0.0;
    for(i=0;i<p->p->spec->type->dim;i++) {
        t = (p->w[i] - n.p->w[i] - delt[i]) * p->p->spec->type->matchFac[i];
	d += t*t;
    }
    d=exp(log(d)*power/2.0);
    return -d;
}


void Matching::match(Spectrum *spec, LiPeak &peaks1, LiPeak &peaks2, 
                     float p, float f)
				    /* peaks1 ... normal peak list    */
				    /* peaks2 ... reference peak list */
				    /* p ... power of norm for dist.  */
				    /* f ... factor for max. shift    */
{
   int np1,np2, ip1, ip2, i;
   Peak *peak;
   Node *liNode1,*liNode2;

   dim=spec->type->dim;
   power=p;

   np1=0;
   peak=(Peak *)peaks1.start();
   while(peak!=0) {
      if(peak->p->spec == spec) np1++;
      peak=(Peak *)peaks1.next();
   }
   liNode1 = new Node[np1];

   np2=0;
   peak=(Peak *)peaks2.start();
   while(peak!=0) {
      if(peak->p->spec == spec) np2++;
      peak=(Peak *)peaks2.next();
   }
   liNode2 = new Node[np2];

   ip1=0;
   peak=(Peak *)peaks1.start();
   while(peak!=0) {
      if(peak->p->spec == spec) {
         liNode1[ip1].p=peak;
	 liNode1[ip1].s=-f;
	 ip1++;
      }
      peak=(Peak *)peaks1.next();
   }

   ip2=0;
   peak=(Peak *)peaks2.start();
   while(peak!=0) {
      if(peak->p->spec == spec) {
         liNode2[ip2].p=peak;
	 liNode2[ip2].s=-f;
	 ip2++;
      }
      peak=(Peak *)peaks2.next();
   }
   
   for(i=0;i<dim;i++) delt[i]=0.0;
   
   domatch(np1, liNode1, np2, liNode2);
   sort(np1, liNode1, np2, liNode2);
   listmatch(np1, liNode1, np2, liNode2);
   
   calcdelt(np1, liNode1, np2, liNode2);
   
   domatch(np1, liNode1, np2, liNode2);
   sort(np1, liNode1, np2, liNode2);
   listmatch(np1, liNode1, np2, liNode2);
   copytoref(np1, liNode1, np2, liNode2);
   copytoco(np1, liNode1, np2, liNode2);

//   delete liNode1;
//   delete liNode2;
}

void Matching::sort(int n1, Node *m1, int n2, Node *m2)
{
   int i, j;
   Node tmp;
   
   for(i=0;i<n1;i++) {
      if(m1[i].m>=0 && m1[i].m<n2) {
         m1[i].score=m2[m1[i].m].calcScore(m1[i]);
      } else {
	 m1[i].score=0.0;
      }
   }

   for(i=1;i<n1;i++) for(j=0;j<i;j++) {
      if(m1[i].score>m1[j].score) {
	  tmp=m1[i]; m1[i]=m1[j]; m1[j]=tmp;
      }
   }  
}

void Matching::copytoref(int n1, Node *m1, int n2, Node *m2)
{
   int i, j;
   
   for(i=0;i<n1;i++) {
      if(m1[i].m>=0 && m1[i].m<n2) {
         for(j=0;j<dim;j++) {
	    m2[m1[i].m].p->p->userAss[j] = m1[i].p->p->userAss[j];
            sprintf(m2[m1[i].m].p->p->comment,"#%d",m1[i].p->p->nr);
         }
      } 
   }
}

void Matching::copytoco(int n1, Node *m1, int n2, Node *m2)
{
   int i, j;
   
   for(i=0;i<n1;i++) {
      if(m1[i].m>=0 && m1[i].m<n2) {
         for(j=0;j<dim;j++) {
	    if(m1[i].p->p->userAss[j]) {
	       m1[i].p->p->userAss[j]->notifyNewFreq(
	          m2[m1[i].m].p->w[j], 0, 0, 
	          m2[m1[i].m].p->p->spec,  ADAPT);
	    }
         }
      } 
   }
}

void Matching::calcdelt(int n1, Node *m1, int n2, Node *m2)
{
   int i, j, n, cnt;
   Node tmp;
   char str[MAXLINE], dist[MAXLINE],  t[MAXLINE];
   Log log;
   
   n=n1/3;
   cnt=0;
   for(j=0;j<dim;j++) delt[j]=0.0;
   for(i=0;i<n;i++) {
      if(m1[i].m>=0 && m1[i].m<n2) {
         cnt++;
         for(j=0;j<dim;j++) {
            delt[j]+=m2[m1[i].m].p->w[j] - m1[i].p->w[j];
         }
      } 
   }
   
   strcpy(dist, "");
   for(j=0;j<dim;j++) {
      if(cnt>0) delt[j]/=cnt;
      if(j==0) sprintf(t, "%6.3f", delt[j]);
      else     sprintf(t, " / %6.3f", delt[j]);
      strcat(dist, t);
   }
   sprintf(str, "\nReference list shifted by ( %s ) ppm\n", dist);
   log.w(str);
}

void Matching::listmatch(int n1, Node *m1, int n2, Node *m2)
{
   int i, cnt;
   float score;
   char str[MAXLINE];
   Log log;

   cnt=0;
   for(i=0;i<n1;i++) if(m1[i].m>=0 && m1[i].m<n2) cnt++;
   
   log.w("\n\n");
   score=calcscore(n1, m1, n2, m2);
   sprintf(str, "Number of matched peaks: %3d\n", cnt);
   log.w(str); 
   sprintf(str, "Total score            : %7.3f\n", score);
   log.w(str); 
   if(cnt>0){
     sprintf(str, "Score per matched peak : %7.3f\n", score/cnt);
     log.w(str);      
   }
   
   log.w("\n\nMatched peaks:\n");
   for(i=0;i<n1;i++) {
      if(m1[i].m>=0 && m1[i].m<n2) {
         sprintf(str, " %7.3f : ", m1[i].score);
	 log.w(str);
         m1[i].p->print(0);
         m2[m1[i].m].p->print(11);
      }
   }  
   
   log.w("\n\nUnmatched peaks:\n");
   for(i=0;i<n1;i++) {
      if(m1[i].m<0 || m1[i].m>n2) {
         m1[i].p->print(3);
      }
   }  
   
   log.w("\n\nUnmatched reference peaks:\n");
   for(i=0;i<n2;i++) {
      if(m2[i].m<0 || m2[i].m>n1) {
         m2[i].p->print(3);
      }
   }  
   
   
}

float Matching::calcscore(int n1, Node *m1, int n2, Node *m2)
{
   float score = 0.0;
   int i;

   for(i=0;i<n1;i++) {
      if(m1[i].m>=0 && m1[i].m<n2) {
         score += m2[m1[i].m].calcScore(m1[i]);
      }
   }
   return score;
}

float Matching::alternatingPath(int n1, Node *m1, int n2, Node *m2)
{
   int i,i1,i2;
   int finished = 0;
   float s;
   float max=0.0;

   /* init labels */
   for(i=0;i<n1;i++) {m1[i].pred=-1; m1[i].score=0.0; }
   for(i=0;i<n2;i++) {m2[i].pred=-1; m2[i].score=0.0; }

   /* for all edge */
   while(!finished) {
      finished=1;
      /* edges from m1 to m2 */
      for(i1=0;i1<n1;i1++) for(i2=0;i2<n2;i2++) {
         if(m1[i1].m!=-1 && m1[i1].pred==-1) continue;
         s = m2[i2].calcScore(m1[i1]) + m2[i2].s;
         if(m1[i1].score + s > m2[i2].score || 
            m2[i2].pred==-1) {
            m2[i2].pred = i1;
            m2[i2].score=m1[i1].score+s;
            if(m2[i2].score>max && m2[i2].m==-1) max=m2[i2].score;
            finished = 0;
         }
      }
      /* edges back to m1 */
      for(i1=0;i1<n1;i1++) {
         i2 = m1[i1].m;
         if(i2==-1) continue;
         s = m2[i2].calcScore(m1[i1]) + m2[i2].s;
         if(m2[i2].score - s > m1[i1].score || m1[i1].pred==-1) {
            m1[i1].pred = i2;
            m1[i1].score= m2[i2].score - s;
            if(m1[i1].score>max) max=m1[i1].score;
            finished = 0;
         }
      }
//      printf("Max: %7.3f\n", max);
   }
   return max;
}

void Matching::updateCycle(int n1, Node *m1, int n2, Node *m2)
{
   float max=0.0;
   int k1=-1;
   int k2=-1;
   int i;
   int toggle;

   for(i=0;i<n2;i++)
      if(m2[i].score>max && m2[i].m==-1) {max=m2[i].score;k1=-1;k2=i;}
   for(i=0;i<n1;i++)
      if(m1[i].score>max) {max=m1[i].score;k1=i; /*printf("\nfound\n"); */}

   if(k1!=-1) {toggle=1; m1[k1].m = -1;}
   else       {toggle=2; }
   do{
      if(toggle==1) {
         k2=m1[k1].pred;
         toggle=2;
      } else  {
         k1=m2[k2].pred;
         m1[k1].m = k2;
         m2[k2].m = k1;
         toggle=1;
      }
   } while (k1!=-1 && k2!=-1);
}

void Matching::domatch(int n1, Node *m1, int n2, Node *m2)
{
   int i;
   float score;
//   char str[MAXLINE];
//   Log log;

   for(i=0;i<n1;i++) m1[i].m=-1;
   for(i=0;i<n2;i++) m2[i].m=-1;

   while(alternatingPath(n1,m1,n2,m2)>0.0) {
       updateCycle(n1,m1,n2,m2);
       score=calcscore(n1,m1,n2,m2);
//       sprintf(str,"\n    current score: %7.3f",score);
//       log.w(str);
   }
}