/*
************************************************************************
*
*   prep.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
*
************************************************************************
*/
/**************************************************************************/
/*                             prep.C                                     */
/*                                                                        */
/* class to define prprocessing of peak lists                             */
/**************************************************************************/

#include <string.h>
#include <stdio.h>
#include "global.h"			/* global constants               */
#include "log.h"			/* general log file               */
#include "bib.h"
#include "prep.h"
#include "specTypes.h"

Compon::Compon(SpecType *s, int *d)
{
   int i;

   for(i=0;i<DIMENSION;i++) {dims[i]=d[i];}
   type=s;
}

EquivDim::EquivDim(char *n, Compon *r)
{
   strcpy(name,n);
//   ref=r;
//   printf("initialized with %s\n",r->type->name);
   compo.add(r);
}

void EquivDim::addDim(int is, int d1, int d2, int d3)
{
   Compon *compon;
   LiSpecType liSpecType;

   int d[DIMENSION];
   int i;

   for(i=0;i<DIMENSION;i++) {d[i]=-1;}
   d[0]=d1-1;
   d[1]=d2-1;
   d[2]=d3-1;

   compon=new Compon(&(liSpecType.l[is]),d);
//   printf("added %s\n",compon->type->name);
   compo.add(compon);
}

void EquivDim::filter(Spectrum *rspec, Spectrum *spec, char *pname, char *coname)
{
   BibAccess bib;
   PeakWriter writer;
   Log log;
   Compon *com, *c, *ref;
   Peak *peak,*refp,*wpeak;
   char str[MAXLINE];
   int ok,i,id,cnt,peaknr,t,lim,title;
   int nrUnass=0, nrDegen=0, nrPeak=0, nrUserAss=0;
   float low[DIMENSION],up[DIMENSION],w[DIMENSION];

   com=0;
   for(i=0;i<usedSpec.used;i++) {
      c=(Compon *)compo.start();
      while(c!=0) {
         if( spec->type == c->type) com=c;
         if(rspec->type == c->type) ref=c;
         c=(Compon *)compo.next();
      }
   }
   ok=1;
   if(com==0) {
      sprintf(str,"ERROR: spectrum %s does not belong to %s\n",spec->type->name,name);
      ok=0;
      log.warn(-1,str);
   }
   if(ref==0) {
      sprintf(str,"ERROR: spectrum %s does not belong to %s\n",rspec->type->name,name);
      ok=0;
      log.warn(-1,str);
   }
   if(ok) {
//      /* report for each peak of spec */
//      peak=(Peak *)spec->liPeakDim[1].start();
//      while(peak!=0) {
//         /* set limits for search */
//         for(id=0;id<rspec->type->dim;id++) {
//            low[id]=rspec->peaks.minshift[id]-0.001;
//            up[id] =rspec->peaks.maxshift[id]+0.001;
//         }
//         for(id=0;id<rspec->type->dim;id++) {
//            if(com->dims[id]>=0 && ref->dims[id]>=0) {
//             low[ref->dims[id]]=peak->w[com->dims[id]] 
//                                  - spec->type->accSpec[com->dims[id]];
//               up[ref->dims[id]]=peak->w[com->dims[id]] 
//                                  + spec->type->accSpec[com->dims[id]];
//            }
//         }
//         /* for all matching peaks in reference spectrum */
//         cnt=0;
//	 refp=rspec->peaks.start(low,up);
//         while(refp!=0) {
//            cnt++;
//            refp=rspec->peaks.next();
//         }
//	 if(cnt==0) {
//           log.w("No matching peak for\n");
//           peak->print(3);
//	 } else if (cnt>1) {  
//           log.w("Multiple possiblities for\n");
//           peak->print(3);
//           refp=rspec->peaks.start(low,up);
//           while(refp!=0) {
//              refp->print(3);
//              refp=rspec->peaks.next();
//           }
//	 }
//         peak=(Peak *)spec->liPeakDim[1].next();
//      }
      for(lim=0;lim<=10;lim++) {
          title=0;
	  /* report for each peak of reference */
	  peak=(Peak *)rspec->liPeakDim[1].start();
	  while(peak!=0) {
	     /* set limits for search */
	     for(id=0;id<spec->type->dim;id++) {
		low[id]=spec->peaks.minshift[id]-0.001;
		up[id] =spec->peaks.maxshift[id]+0.001;
	     }
	     for(id=0;id<spec->type->dim;id++) {
		if(com->dims[id]>=0 && ref->dims[id]>=0) {
		   low[com->dims[id]]=peak->w[ref->dims[id]] 
				      - rspec->type->accSpec[ref->dims[id]];
		   up[com->dims[id]]=peak->w[ref->dims[id]] 
				      + rspec->type->accSpec[ref->dims[id]];
		}
	     }
	     /* for all matching peaks in spectrum */
	     cnt=0;
	     refp=spec->peaks.start(low,up);
	     while(refp!=0) {
		cnt++;
		refp=spec->peaks.next();
	     }
	     if(cnt==lim || (lim==10 && cnt>=10)) {
               if(!title) {
	          title=1;
		  if(lim==0) {
		      log.w("\n\nNo matching peaks for\n");
		  } else if(lim<10) {
		      sprintf(str, "\n\n%d possibilities for\n", lim);
		      log.w(str);
		  } else {
		      log.w("\n\nTen or more possibilities for\n");
		  }
               }
	       peak->print(3);
	       refp=spec->peaks.start(low,up);
	       while(refp!=0) {
		  refp->print(3);
		  refp=spec->peaks.next();
	       }
	       log.w("\n");
	     }
	     peak=(Peak *)rspec->liPeakDim[1].next();
	  }
	  /* write peak list */
      }
      peaknr=1;
      if(strlen(pname)>0) {
       if(bib.coLoad(coname,DESTFRAG)) {
         writer.dim = spec->type->dim;
         for(i=0;i<writer.dim;i++) writer.perm[i] = spec->type->permut[i];
         writer.start(pname);
         /* for all peaks of spec */
         peak=(Peak *)spec->liPeakDim[1].start();
         while(peak!=0) {
            /* set limits for search */
            for(id=0;id<rspec->type->dim;id++) {
               low[id]=rspec->peaks.minshift[id]-0.001;
               up[id] =rspec->peaks.maxshift[id]+0.001;
	       w[id]=peak->w[id];
            }
            for(id=0;id<rspec->type->dim;id++) {
               if(com->dims[id]>=0 && ref->dims[id]>=0) {
                  low[ref->dims[id]]=peak->w[com->dims[id]] 
                                  - spec->type->accSpec[com->dims[id]];
                  up[ref->dims[id]]=peak->w[com->dims[id]] 
                                  + spec->type->accSpec[com->dims[id]];
               }
            }
            /* for all matching peaks in reference spectrum */
            cnt=0;
	    refp=rspec->peaks.start(low,up);
            while(refp!=0) {
               cnt++;
	       for(id=0;id<rspec->type->dim;id++) {
                 if(com->dims[id]>=0 && ref->dims[id]>=0) {
		    w[com->dims[id]]=refp->w[ref->dims[id]];
                 }
               }
               wpeak = new Peak(peak->p, w);
               t=peak->p->nr;
               peak->p->nr=peaknr;
               peaknr++;
               wpeak->write(writer,DESTFRAG,nrUnass,nrDegen,nrUserAss);
	       delete wpeak;
               peak->p->nr=t;
               nrPeak++;
               refp=rspec->peaks.next();
            }
            peak=(Peak *)spec->liPeakDim[1].next();
         }
         writer.close();
       } else {
          sprintf(str,"\n... atom list '%s' could not be read",coname);
          log.warn(-1,str);
       }
       sprintf(str,"\n... %d peaks written to file '%s'",nrPeak,pname);
       log.w(str);
       sprintf(str,"\n... %d peaks without assignments",nrUnass);
       log.w(str);
       sprintf(str,"\n... %d peaks with multiple assignments",nrDegen);
       log.w(str);
       sprintf(str,"\n... %d peaks with unchanged user assignments",nrUserAss);
       log.w(str);
      }
   }
}


int LiEquivDim::compare(const Element *a1,const Element *a2)
{
   return(
     strcmp(((EquivDim *)a2)->name,((EquivDim *)a1)->name)
   );
}

void LiEquivDim::addEq(char *n, int is, int d1, int d2, int d3)
{
   Compon *compon;
   LiSpecType liSpecType;
   EquivDim *equivDim;
   int d[DIMENSION];
   int i;

   for(i=0;i<DIMENSION;i++) {d[i]=-1;}
   d[0]=d1-1;
   d[1]=d2-1;
   d[2]=d3-1;

   compon=new Compon(&(liSpecType.l[is]),d);
   equivDim=new EquivDim(n,compon);
   add(equivDim);
}

EquivDim *LiEquivDim::start(char *n, int &i)
{
   i=0;
   return next(n,i);
}

EquivDim *LiEquivDim::next(char *nam, int &i)
{
   int len;
   EquivDim *res = 0;

   len =strlen(nam);
   while(i<n && !res) {
      if(strncmp(((EquivDim *)l[i])->name,nam,len)==0) res=(EquivDim *)l[i];
      i++;
   } 
   return res;
}