/*
************************************************************************
*
*   erPeak.C -
*
*   Copyright (c) 1995
*
*   ETH Zuerich
*   Institut fuer Molekularbiologie und Biophysik
*   ETH-Hoenggerberg
*   CH-8093 Zuerich
*
*   All Rights Reserved
*
*   Date of last modification : 95/09/15
*   Pathname of SCCS file     : /export/home3/cb/garant-1.0/src/SCCS/s.erPeak.C
*   SCCS identification       : 1.3
*
************************************************************************
*/
/**************************************************************************/
/*                           erPeak.cc                                    */
/*                                                                        */
/* expected peaks                                                         */
/**************************************************************************/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "global.h"
#include "log.h"
#include "atom.h"
#include "co.h"
#include "frag.h"
#include "bib.h"
#include "erPeak.h"
#include "specTypes.h"
#include "parser.h"
#include "functions.h"
#include "specs.h"

ErPeak *LiErPeak::start(Iterator &i, Spectrum *S)
{
   init(i,S);
   return(next(i,S));
}

void LiErPeak::init(Iterator &it, Spectrum *S)
{
   int i,j,k,res;
   if(n==0) it.cur=n;
   else {
      i=0; j=n-1;  res=1;
      while(i<j) {
         k=(i+j)/2;
         res= ((ErPeak *)l[k])->spec->nr - S->nr;
         if(res<0) i=k+1;
         else j=k;
      }
      k=(i+j)/2; res= ((ErPeak *)l[k])->spec->nr - S->nr;
      if(res==0) it.cur=k-1;
      else it.cur=n;
   }
}

ErPeak *LiErPeak::next(Iterator &i, Spectrum *S)
{
   int res;

   i.cur++;
   if(i.cur>=n) return(0);
   res= ((ErPeak *)l[i.cur])->spec->nr - S->nr;
   if(res==0) return((ErPeak *)l[i.cur]);
   else { i.cur=n; return(0); }
}

ErPeak *LiErPeak::startPars(Iterator &i, Spectrum *S, char *t)
{
   init(i,S);
   return(nextPars(i,S,t));
}

ErPeak *LiErPeak::nextPars(Iterator &i, Spectrum *S, char *t)
{
   ErPeak *peak;

   peak = next(i,S);
   while(peak != 0 && !peak->match(t)) peak = next(i,S);
   return peak;
}


void LiErPeak::print(int off)
{
   Log log;
   char str[MAXLINE];

   sprintf(str,"%*sList of expected peaks:\n",off,"");
   log.w(str);
   Index::print(off+3);
}

int LiErPeak::compare(const Element *p1,const Element *p2)
{
   int res;
   int i=0;
   res = ((ErPeak *)p2)->spec->nr - ((ErPeak *)p1)->spec->nr;
   while(res==0 && i<((ErPeak *)p1)->spec->type->dim) {
      res = ((ErPeak *)p2)->liCo[i]->t - ((ErPeak *)p1)->liCo[i]->t;
      if(!res) {
         res = ((ErPeak *)p2)->liCo[i]->nuclei[0]->nr
               - ((ErPeak *)p1)->liCo[i]->nuclei[0]->nr;
         if(!res) res = ((ErPeak *)p2)->type - ((ErPeak *)p1)->type;
      }
      i++;
   }
   return(res);
}


ErPeak::ErPeak(Coherence **Co, float p, Spectrum *spectrum,int ty)
{
   int i;
   float pE;

   t = EXPPEAK;
   currAss=0;
   assPop=0;
   pE   = p;

   if(pE>=1.0)      group = 0;
   else if(pE>=0.8) group = 1;
   else if(pE>=0.3) group = 2;
   else if(pE>=0.1) group = 3;
   else             group = 4;

   pExist = p;
   pSelect = p;
   importance   = 1.0;
   spec = spectrum;
   type = ty;
   liEquiv = 0;
   liPseud = 0;
   status = P_UNASS;
   active = TRUE;

   for(i=0;i<spec->type->dim;i++) {liCo[i] = Co[i]; pos[i]=-1;}
}

ErPeak::ErPeak(Coherence **Co, Spectrum *spectrum,int ty)
{
   int i;

   t = EXPPEAK;
   spec = spectrum;
   type = ty;
   liEquiv = 0;
   liPseud = 0;
   for(i=0;i<spec->type->dim;i++) liCo[i] = Co[i];
}

ErPeak::ErPeak()
{
   t=EXPPEAK;
   currAss=0;
   assPop=0;
   liEquiv=0;
   liPseud = 0;
}

void ErPeak::getIntensity()
{
   float ppm[DIMENSION];	/* assigned frequency in ppm */
   float points[DIMENSION];	/* position in spectrum      */
   short int fold[DIMENSION];	/* necessary folding         */
   int res;
   int i;

   if(spec->specdata==0) return;
   for(i=0;i<spec->type->dim;i++) {
      ppm[i] = (liCo[i]->meanFreq(spec,i)).x;
      if(ppm[i]==BAD_PPM) {
        pos[i]=-1;
        intensity=0.0;
        return;
      }
   }
//   printf("unpermuted: ");
//   for(i=0;i<spec->type->dim;i++) {
//      printf("%6.2f ", ppm[i]);
//   }
   spec->type->unpermute(ppm);
//   printf("\npermuted:   ");
//   for(i=0;i<spec->type->dim;i++) {
//      printf("%6.2f ", ppm[i]);
//   }
   Fold(spec->specdata,ppm,fold);
//   printf("\n  folded:   ");
//   for(i=0;i<spec->type->dim;i++) {
//      printf("%6.2f ", ppm[i]);
//   }
   P_Sp(spec->specdata,ppm,points);
//   printf("\n  pointd:   ");
//   for(i=0;i<spec->type->dim;i++) {
//      printf("%6.2f ", points[i]);
//   }
   res=1;
   for(i=0;i<spec->type->dim;i++) {
      if((int)points[i] != pos[i]) {
         res = 0;
         pos[i]=(int)points[i];
      }
   }
   if(res) return;
   intensity = Read_Spec_Point(spec->specdata,pos);
//   printf("  --> %7d\n",intensity);
//   print(5);
}

int ErPeak::isPresent()
{
//   int i;

   if(spec->type->noise == 0) return(0);

   getIntensity();
   if(spec->specdata==0) return(0);
   if(nabs(intensity) > spec->type->noise) {
//     printf("Peak at: ");
//     for(i=0;i<spec->type->dim;i++) {
//	printf("%6.2f ", (liCo[i]->meanFreq(spec,i)).x);
//     }
//     printf("\n");
//     print(5);
     return(1);
   } else {
     return(0);
   }
}

ErPeak *newRe(LiErPeak &liP, Coherence **Co, float prob, Spectrum *spec, int ty)
{
   ErPeak *ret;

   ret = (ErPeak *)liP.find(ErPeak(Co,spec,ty));
   if(!ret) {
      ret = new ErPeak(Co,prob,spec,ty);
      (void)liP.insert(ret);
   }
   return(ret);
}

void ErPeak::print(int off)
{
   Log log;
   char str[MAXLINE];
   char n[MAXLINE];
   int nrEq=0;

   name(n);
   if(liEquiv != 0) nrEq=liEquiv->getn();
   if(type == DESTFRAG) {
      sprintf(str,
   "%*s%-32s p(exist)=%7.3g imp=%7.3f nrEq=%d\n",
         off,"",n,pExist,imp.sig,nrEq);
      log.w(str);
      if(getAss() != 0) getAss()->print(off+5);
      //printEquiv();
   } else {
      sprintf(str,"%*s%-32s subfragment peak\n",
         off,"",n);
      log.w(str);
   }
}

int ErPeak::match(char *t)
{
   int nrRes,tb,ib;
   Parser pars;
   char str[MAXLINE];
   int res;

   nrRes = UNDEFINED_RES_NR;
   name(str);
   res = pars.peak(t,0,tb,str,0,ib,nrRes);
   res &= (strlen(str) == ib);
   return res;
}


void ErPeak::name(char *name)
{
   char str[MAXLINE];
   int i;

   sprintf(name,"%s: ",spec->type->name);
   for(i=0;i<spec->type->dim;i++) {
      liCo[i]->fullNameGet(str);
      strcat(name,str);
      if(i != spec->type->dim - 1) strcat(name," ");
   }
}

int ErPeak::getWriteFreq(Wert *w)
{
   int res = 0;
   int i;

   for(i=0;i<spec->type->dim;i++) {
      w[i] = liCo[i]->meanFreq(spec,i);
      if(w[i].x == BAD_PPM) w[i] = liCo[i]->assF;
      if(w[i].x == BAD_PPM) {
         res++;
         w[i] = liCo[i]->expFreq();
      }
   }
   return res;
}

int ErPeak::write(PeakWriter &writer, int typeFrag, int maxUndef)
{
   int i;
   int res = 0;

   if(type == typeFrag) {
      writer.v = Wert(pExist,0);
      writer.intflag = '-';
      res = getWriteFreq(writer.w);
      for(i=0; i<writer.dim; i++) {
         if(liCo[i]->statPseudo==PSEUDON)
            liCo[i]->pseudCo->nameGet(writer.coAtNam[i], writer.coResNum[i]);
         else
            liCo[i]->nameGet(writer.coAtNam[i], writer.coResNum[i]);
      }
      if(res <= maxUndef) {
         writer.color=res+1;
         strcpy(writer.comment,"");
         if(minDist >= 0) {
             writer.color = xnint(minDist);
             if(writer.color < 1) writer.color = 1;
             if(writer.color > 6) writer.color = 6;
             sprintf(writer.comment," %3.1f",minDist);
         } 
         writer.write();
      }
   }

   return (res <= maxUndef);
}

void ErPeak::resetAss()
{
   int i;

   // printf("resetAss: "); print(0);
   if(status == P_ASS) {
      for(i=0;i<spec->type->dim;i++) {
         liCo[i]->notifyFreqReset(currAss->w[i],spec);
      }
      currAss->p->nrAss--;
      currAss->p->retractAsslink(asslinkindex);
      currAss = 0;
   }
   status = P_UNASS;
   intensity = 0;
}

void ErPeak::rejectAss()
{
   // printf("rejectAss: "); print(0);
   if(currAss != 0) {
      currAss->p->nrAss--;
      currAss->p->retractAsslink(asslinkindex);
      currAss = 0;
   }
   status = P_REJECTED;
   // getIntensity();
}

void ErPeak::makeAss(Peak *peak, const int *decision,ErPeak *parPeak,int parNr)
{
   int i;

   for(i=0;i<spec->type->dim;i++) {
      if(parPeak==0) {
         liCo[i]->notifyNewFreq(peak->w[i],peak->w[i],0,spec,decision[i]);
      } else {
         liCo[i]->notifyNewFreq(
            peak->w[i],
            parPeak->getCo(i)->assPop[parNr].meanF,
            &(parPeak->getCo(i)->assPop[parNr].measF[0]),
            spec,decision[i]);
      }
   }
   currAss = peak;
   peak->p->nrAss++;
   peak->p->establishAsslink(this);
   status = P_ASS;
}

void ErPeak::statCo(int *s)
{
   int i;
   for(i=0;i<spec->type->dim;i++)  s[i] = liCo[i]->status;
}

int ErPeak::nrUnAssCo()
{
   int i,count=0;
   for(i=0;i<spec->type->dim;i++) if(liCo[i]->status==UNASSIGNED)count++;
   return count;
}

void ErPeak::updateScoreCo(float *scores, float *mean, float *deviation)
{
   int i;

   for(i=0;i<spec->type->dim;i++) {
      liCo[i]->updateScores();
      scores[i]=liCo[i]->localScore;
      mean[i] = liCo[i]->getScoreMean();
      deviation[i] = liCo[i]->getScoreDev();
   }
}

void ErPeak::setSumPeakP()
{
   int i;

   // sumPeakP = spec->bonusPeak * pExist;
   sumPeakP = 0;
   for(i=0;i<spec->type->dim;i++) {
      sumPeakP += spec->bonusAss[i][0] * pExist;
      sumPeakP += liCo[i]->bonusFreq;
   }
}

void ErPeak::freq(int *s, Wert *w)
{
   int i;
   for(i=0;i<spec->type->dim;i++) {
      if(s[i] == ASSIGNED) w[i] = liCo[i]->meanFreq(spec,i);
      else w[i] = liCo[i]->expFreq();
   }
}

void ErPeak::expFreq(Wert *w)
{
   int i;
   for(i=0;i<spec->type->dim;i++) w[i] = liCo[i]->expFreq();
}

void ErPeak::assFreq(Wert *w)
{
   int i;
   for(i=0;i<spec->type->dim;i++) w[i] = liCo[i]->meanFreq(spec,i);
}

int ErPeak::range(int *s, float *low, float *up)
{
   int i;
   int nrAssigned=0;

   for(i=0;i<spec->type->dim;i++) {
      if(s[i] == ASSIGNED) {
         liCo[i]->meanRange(spec,i,low[i],up[i]);
         nrAssigned++;
      } else {
         liCo[i]->expRange(low[i],up[i]);
      }
   }
   return nrAssigned;
}

int ErPeak::testDiagUnDef(int *flag)
{
   int i,j;
   int res;
   res = 0;
   for(i=0;i<spec->type->dim && !res;i++)
      if(flag[i] == UNASSIGNED)
         for(j=0;j<spec->type->dim && !res;j++)
            if(i != j && liCo[i] == liCo[j]) res = 1;

   return res;
}

int ErPeak::testSeqRange()
{
   int i,j;
   int res;
   res = 0;
   for(i=0;i<spec->type->dim;i++) {
      for(j=i+1;j<spec->type->dim;j++) {
         res=max(res,nabs(liCo[i]->firstAt()->frag->nrExt - 
                          liCo[j]->firstAt()->frag->nrExt));
      }
   }
   return res;
}

int ErPeak::equiv(ErPeak *peak)
{
   int res=1;
   int i;
   if(spec != peak->spec || testSeqRange() != peak->testSeqRange()) {
      res=0;
   } else {
      for(i=0;res && i<spec->type->dim;i++)
         res = liCo[i]->equiv(peak->liCo[i]);
   }
   return res;
}

int ErPeak::pseudIdentical(ErPeak *peak)
{
   int res=1;
   int i;
   if(spec != peak->spec) {
      res=0;
   } else {
      for(i=0;res && i<spec->type->dim;i++)
         res = liCo[i]->pseudIdentical(peak->liCo[i]);
   }
   return res;
}

int ErPeak::floatIdentical(ErPeak *peak)
{
   int res=1;
   int i;
   if(spec != peak->spec) {
      res=0;
   } else {
      for(i=0;res && i<spec->type->dim;i++)
         res = liCo[i]->floatIdentical(peak->liCo[i]);
   }
   return res;
}

void ErPeak::addEquiv(ErPeak *peak)
{
   if(liEquiv == 0) {
      if(peak->liEquiv == 0) liEquiv = new LiEquiv();
      else liEquiv=peak->liEquiv;
   }
   if(peak->liEquiv == 0) peak->liEquiv = liEquiv;
   liEquiv->insert(this);
   liEquiv->insert(peak);
}

void ErPeak::addPseudo(ErPeak *peak)
{
   if(liPseud == 0) liPseud = new LiErPeak(4);
   liPseud->insert(this);
   liPseud->insert(peak);
   peak->liPseud = liPseud;
}

void ErPeak::printEquiv()
{
   ErPeak *eq;
   Log log;
   char str[MAXLINE];
   char n[MAXLINE];

   name(n);
   sprintf(str,"\n... %s:",n);
   log.w(str);
   if(liEquiv != 0) {
      eq = (ErPeak *)liEquiv->start();
      while(eq != 0) {
         eq->name(n);
         sprintf(str,"\n       %s:",n);
         log.w(str);
         eq = (ErPeak *)liEquiv->next();
      }
   }
}


Coherence *ErPeak::startCo(Iterator &i)
{
   i.cur = 0;
   return nextCo(i);
}

Coherence *ErPeak::nextCo(Iterator &i)
{
   Coherence *res=0;
   if(i.cur < spec->type->dim) {
      res = liCo[i.cur];
      i.cur++;
   }
   return res;
}

Coherence *ErPeak::getCo(Iterator i)
{
   Coherence *res=0;
   if(i.cur < spec->type->dim) {
      res = liCo[i.cur];
   }
   return res;
}
      
Atom *ErPeak::startAt(Iterator &iAt,  Iterator &iCo)
{
   Coherence *co;
   Atom *at;

   co = startCo(iCo);
   if(co != 0) at = co->start(iAt);
   else at = 0;

   return at;
}


Atom *ErPeak::nextAt(Iterator &iAt,  Iterator &iCo)
{
   Coherence *co;
   Atom *at = 0;

   co = getCo(iCo);
   while(co!=0 && at==0) {
      at = co->next(iAt);
      if(at==0) {
         co=nextCo(iCo);
         at=co->start(iAt);
      }
   }
   return at;
}

ErPeak *ErPeak::startNeighbour(int &iP,  int &iCo, int *decision)
{
   iCo = 0;
   while((iCo<spec->type->dim) && (decision[iCo] != RESET)) {iCo++; }
   iP = 0;
   return nextNeighbour(iP, iCo, decision);
}


ErPeak *ErPeak::nextNeighbour(int &iP,  int &iCo, int *decision)
{
   Coherence *co;
   ErPeak *p = 0;

   while(iCo < spec->type->dim && p==0) {
      co = liCo[iCo];
      while(iP < co->nrErPeak && p==0) {
         if(co->peaks[iP]->type ==  type) p = co->peaks[iP];
         iP++;
      }
      if(p==0) {
         iP = 0;
         iCo++;
         while(iCo < spec->type->dim && decision[iCo] != RESET) iCo++;
      }
   }
   return p;
}

ErPeak *ErPeak::startNeighbour(int &iP,  int &iCo)
{
   iCo = 0;
   iP = 0;
   return nextNeighbour(iP, iCo);
}


ErPeak *ErPeak::nextNeighbour(int &iP,  int &iCo)
{
   Coherence *co;
   ErPeak *p = 0;

   while(iCo < spec->type->dim && p==0) {
      co = liCo[iCo];
      while(iP < co->nrErPeak && p==0) {
         if(co->peaks[iP]->type ==  type) p = co->peaks[iP];
         iP++;
      }
      if(p==0) {
         iP = 0;
         iCo++;
      }
   }
   return p;
}

int ErPeak::cmpAss(AssRequest &aR)
{
   int i,d;
   int res = 1;

   d=spec->type->dim;
   for(i=0;i<d && res;i++) {
      if(aR.nrFrag[i] != -1 &&
         liCo[i]->firstAt()->frag->nrExt != aR.nrFrag[i]) res = 0;
   }
   for(i=0;i<d && res;i++) {
      if(!aR.equivName(i,liCo[i]->firstAt()->name)) res = 0;
   }
   return res;
}