/*
************************************************************************
*
*   peak.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.peak.C
*   SCCS identification       : 1.4
*
************************************************************************
*/
/**************************************************************************/
/*                             peak.cc                                    */
/*                                                                        */
/* measured peaks                                                         */
/**************************************************************************/

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

class ErPeak;

/* selection distribution for expected chemical shifts */
float expectDistr[] = {
   1.0,
   0.98,
   0.95,
   0.9,
   0.8,
   0.5,
   0.1
};

/* selection distribution for assigned chemical shifts */
float assDistr[] = {
   1.0,
   0.99,
   0.98,
   0.97,
   0.96,
   0.5,
   0.001
};

/* number of entries in the distributions */
#define RESOLUTION 6

/* data structures used for n dimensional sorted  peak list ***************/

struct Limit {		/* splitting point: split list into peaks whose shift*/
   float lim;		/* in the dimension dim is smaller than lim and those*/ 
   int dim;		/* whose shift is larger than lim.                   */
};

class Sublist {		/* defines sublists [von,bis)                        */
public:
   Sublist(int v, int b) {von=v;bis=b;}
   Sublist() {von=0;bis=0;}
   int von;
   int bis;
};

/* implementation **********************************************************/

int LiPeakName::compare(const Element *p1, const Element *p2)
{
   Atom *a1,*a2;
   int res = 0;
   int i;

   i = 0;
   res = ((Peak *)p2)->p->spec->nr - ((Peak *)p1)->p->spec->nr;
   while(res == 0 && i < ((Peak *)p2)->p->spec->type->dim) {
      if(((Peak *)p1)->p->userAss[i] != 0 &&
                                      ((Peak *)p2)->p->userAss[i] != 0) {
         a1 = ((Peak *)p1)->p->userAss[i]->firstAt();
         a2 = ((Peak *)p2)->p->userAss[i]->firstAt();
         if(!res) res = a2->frag->nrExt - a1->frag->nrExt;
         if(!res) res = strcmp(a2->name, a1->name);
         if(!res) res = a2->frag->type - a1->frag->type;
      } 
      else {
         res = ((Peak *)p1)->p->userAss[i] - ((Peak *)p2)->p->userAss[i];
      }
      i++;
   }
   return(res);
}

int LiPeak::compare(const Element *p1, const Element *p2)
{
   int res;
   int i;
   res = ((Peak *)p2)->p->nr - ((Peak *)p1)->p->nr;
   if(res==0) res = ((Peak *)p2)->p->spec->nr - ((Peak *)p1)->p->spec->nr;
   i = 0;
   while(res == 0 && i < ((Peak *)p2)->p->spec->type->dim) {
      res = ((Peak *)p2)->fold[i] - ((Peak *)p1)->fold[i];
      i++;
   }
   return(res);
}

int LiPeak::load(Spectrum *spec,char *fn,int flagRefPeaks)
{
   PeakRequest pR;
   Peak *p;
   PeakInfo *pi;
   int res = 1;
   int nrPeaks;
   Peak *peak;
   Log log;
   char str[MAXLINE];
   int nrAdded = 0;
   int i;

   pR.start(fn,spec->type->dim,spec->type->noise,spec->type->permut);
   if(pR.eof) res = 0;
   while(!pR.eof) {
      pi= new PeakInfo(pR.nr, pR.vol, spec, pR.prob, pR.color, pR.comment,
                       pR.intflag);
      p = new Peak(pi, pR.w);
      if(!insert(p)) {
         sprintf(str,"\n... peak nr. %d of spectrum %s %d allready defined - not added to list",pR.nr,spec->type->name,spec->nr);
         log.warn(10,str);
      } else {
         nrAdded++;
      }
      pR.read();
   }


   nrPeaks=0;
   peak = (Peak *)start();
   while(peak != 0) {
      if(peak->p->spec == spec) nrPeaks++;
      peak = (Peak *)next();
   }

   if(!flagRefPeaks) {
       spec->peaks.init(nrPeaks,spec->type->dim);
    
       peak = (Peak *)start();
       while(peak != 0) {
	  if(peak->p->spec == spec) {
	     spec->peaks.add(peak);
	     for(i=0;i<spec->type->dim;i++) spec->liPeakDim[i].add(peak);
	  }
	  peak = (Peak *)next();
       }
    
       spec->peaks.sort();
       for(i=0;i<spec->type->dim;i++) spec->liPeakDim[i].sort();
   }

   sprintf(str,"\n... %d peaks added, %d peaks defined",nrAdded,nrPeaks);
   log.w(str);

   return(res);
}

void LiPeak::unfold(Spectrum *spec,int d, int von, int bis)
{
   Peak *p;
   Peak *peak;
   Log log;
   char str[MAXLINE];
   int nrPeaks;
   int nrAdded = 0;
   short int fold[DIMENSION];
   float freq[DIMENSION];
   int nrFold,i;
   Spec_Param_Data *s;
   LiPeak templist;
   float minFreq=0.0,maxFreq=0.0;
   int minmaxSet = 0;

   d--;
   s = spec->specdata;
   if(s == 0) {
      sprintf(str,
          "\n... WARNING: no parameter file (*.3D.param) specified");
      log.warn(-1,str);
   } else if (d >= spec->type->dim || d < 0)  {
      log.warn(-1,"\n... WARNING: invalid dimension of spectrum");
   } else {
      for(nrFold = von; nrFold<=bis; nrFold++) {
         if(nrFold != 0) {
            for(i=0;i<DIMENSION;i++) fold[i] = 0;
            fold[d] = nrFold;
            peak = (Peak *)start();
            while(peak != 0) {
               if(peak->p->spec == spec && peak->fold[d] == 0) {
                  for(i=0;i<DIMENSION;i++) freq[i] = peak->w[i];
                  spec->type->unpermute(freq);
                  if(!minmaxSet) {
                     minFreq = freq[d]; maxFreq = freq[d];
                     minmaxSet = 1;
                  } else {
                     if(minFreq > freq[d]) minFreq = freq[d];
                     if(maxFreq < freq[d]) maxFreq = freq[d];
                  }
                  UnFold(s,freq,fold);
                  if(minFreq > freq[d]) minFreq = freq[d];
                  if(maxFreq < freq[d]) maxFreq = freq[d];
                  spec->type->permute(freq);
                  p = new Peak(peak->p, freq);
                  for(i=0;i<spec->type->dim;i++)
                     p->fold[i] = peak->fold[i] + fold[i];
                  if(!templist.insert(p)) {
                     delete p;
                     sprintf(str,"\n... WARNING: peak %d allready defined - not added to list",peak->p->nr);
                     log.warn(10,str);
                  } else {
                     nrAdded++;
                  }
               }
               peak = (Peak *)next();
            }
         }
      }
   }

   peak = (Peak *)templist.start();
   while(peak != 0) {
      if(!insert(peak)) {
         sprintf(str,"\n... WARNING: peak %d allready defined - not added to list",peak->p->nr);
         log.warn(10,str);
         delete peak;
      }
      peak = (Peak *)templist.next();
   }

   nrPeaks=0;
   peak = (Peak *)start();
   while(peak != 0) {
      if(peak->p->spec == spec) nrPeaks++;
      peak = (Peak *)next();
   }

   spec->peaks.init(nrPeaks,spec->type->dim);

   peak = (Peak *)start();
   while(peak != 0) {
      if(peak->p->spec == spec) {
         spec->peaks.add(peak);
         for(i=0;i<spec->type->dim;i++) spec->liPeakDim[i].add(peak);
      }
      peak = (Peak *)next();
   }

   spec->peaks.sort();
   for(i=0;i<spec->type->dim;i++) spec->liPeakDim[i].sort();

   sprintf(str,"\n... %d peaks added, %d peaks defined, frequency range %7.3g...%7.3g",nrAdded,nrPeaks,minFreq,maxFreq);
   log.w(str);

}

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

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

void LiPeak::listPeaks()
{
   Log log;

   log.w("\n   List of all measured peaks:\n");
   print(3);
}

void LiPeak::writeMeasPeaks(Spectrum *spec, char *pname,
                         char *coname, int typeFrag, int flag)
{
   BibAccess bib;
   Peak *peak;
   PeakWriter writer;
   int i;
   int res;
   Log log;
   char str[MAXLINE];
   int nrUnass=0, nrDegen=0, nrPeak=0, nrUserAss=0;

   if(bib.coLoad(coname,typeFrag)) {
      writer.dim = spec->type->dim;
      for(i=0;i<writer.dim;i++) writer.perm[i] = spec->type->permut[i];
      writer.start(pname);
      peak = (Peak *)start();
      while(peak != 0 && writer.eof == 0) {
         res = 1; for(i=0;i<spec->type->dim;i++) if(peak->fold[i] != 0) res = 0;
         if(res && peak->p->spec == spec && (!flag || peak->getP() > 0.0)) {
            peak->write(writer,typeFrag,nrUnass,nrDegen,nrUserAss);
            nrPeak++;
         }
         peak = (Peak *)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);
}

void LiPeak::resetUserAss(Spectrum *spec)
{
   Peak *peak;
   int i;
   Log log;
   char str[MAXLINE];
   int count=0;
   int dim;

   dim = spec->type->dim;
   peak = (Peak *)start();
   while(peak != 0) {
      if(peak->p->spec == spec) {
         for(i=0;i<dim;i++) {
            if(peak->p->userAss[i] != 0) {
                peak->p->userAss[i]  = 0;
                count++;
            }
         }
      }
      peak = (Peak *)next();
   }
   sprintf(str,"\n... %d user assignments reset ",count);
   log.w(str);
}

void LiPeak::writeAssPos(Spectrum *spec, char *pname,
                         char *coname, int typeFrag)
{
   BibAccess bib;
   Peak *peak;
   AssWriter writer;
   int i;
   Log log;
   char str[MAXLINE];
   int nrUnass=0, nrDegen=0, nrPeak=0;

   if(bib.coLoad(coname,typeFrag)) {
      writer.dim = spec->type->dim;
      for(i=0;i<writer.dim;i++) {
         writer.perm[i] = spec->type->permut[i];
         writer.tolerance[i] = spec->type->accPick[i];
      }
      strcpy(writer.namCoLi,coname);
      writer.start(pname);
      peak = (Peak *)start();
      while(peak != 0 && writer.eof == 0) {
         if(peak->p->spec == spec) {
            peak->writeAssPos(writer,typeFrag,nrUnass,nrDegen);
            nrPeak++;
         }
         peak = (Peak *)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);
}

void LiPeak::resetAssPos(int type)
{
   Peak *peak;
   Log log;

   peak = (Peak *)start();
   while(peak != 0) {
      peak->p->ass.removeAss(type);
      peak = (Peak *)next();
   }
   log.w("\n... all assignments to measured peaks reset");
}

void LiPeak::mapAssPos(LiErPeak &liErPeak)
{
   Log log;
   char str[MAXLINE];
   Peak *peak;
   ErPeak *poss, *compatibel;
   Iterator iComp;
   int iPoss;
   int count=0;

   peak = (Peak *)start();
   while(peak != 0) {
      poss = peak->p->ass.start(iPoss);
      while(poss != 0) {
         if(poss->type == SUBFRAG) {
            compatibel = liErPeak.start(iComp,poss->spec);
            while(compatibel != 0) {
               if(compatibel->type == DESTFRAG && poss->equiv(compatibel)) {
                  peak->p->ass.addAss(compatibel);
                  count++;
               }
               compatibel = liErPeak.next(iComp,poss->spec);
            }
         }
         poss = peak->p->ass.next(iPoss);
      }
      peak = (Peak *)next();
   }
   sprintf(str,"\n... %d assignment possibilities added",count);
   log.w(str);
}

void LiPeak::cmpAssPos()
{
   Log log;
   char str[MAXLINE];
   char name[MAXLINE];
   Peak *peak;
   ErPeak *poss1, *poss2;
   int iPoss1, iPoss2;
   int countCompatibel;
   int count;
   int error;
   int okPeaks=0,errPeaks=0;
   int titlewritten=0;

   peak = (Peak *)start();
   while(peak != 0) {
      error = 0;
      poss1 = peak->p->ass.start(iPoss1);
      while(poss1 != 0) {
         poss2 = peak->p->ass.start(iPoss2);
         count = 0; countCompatibel=0;
         while(poss2 != 0) {
            if(poss1->type == DESTFRAG && poss2->type == SUBFRAG) {
               count++;
               if(poss2->floatIdentical(poss1)) countCompatibel++;
            }
            poss2 = peak->p->ass.next(iPoss2);
         }
         if(countCompatibel == 0 && count > 0) {
            error = 1;
            if(!titlewritten) {
               log.w("\n... only incompatibel subfragment ass. for following dest. peak ass.:");
               titlewritten = 1;
            }
            poss1->name(name);
            sprintf(str,"\n    peak nr. %4d assigned to %s allowed are",
                    peak->p->nr,name);
            log.w(str);
         }
         poss1 = peak->p->ass.next(iPoss1);
      }
      if(error) {
         errPeaks++;
         poss1 = peak->p->ass.start(iPoss1);
         while(poss1 != 0) {
            if(poss1->type == SUBFRAG) {
               poss1->name(name);
               sprintf(str,"\n                              %s",name);
               log.w(str);
            }
            poss1 = peak->p->ass.next(iPoss1);
         }
      } else {
         okPeaks++;
      }
      peak = (Peak *)next();
   }
   sprintf(str,"\n... %d correct peaks, %d with no corresponding assignment",
           okPeaks,errPeaks);
   log.w(str);
}

Peak *LiPeak::find(int nr,Spectrum *s)
{
   static PeakInfo pi= PeakInfo(nr,s);
   pi  = PeakInfo(nr,s);
   return (Peak *)Index::find(Peak(&pi));
}

PeakInfo::PeakInfo(int peaknr, Wert v, Spectrum *s, float pickProb,
                   int c, char *com, char iflag)
{
   int i;
   nr = peaknr;
   t = PEAK;
   nrAss = 0;
   spec = s;
   vol = v;
   prob= pickProb;
   for(i=0;i<DIMENSION;i++) userAss[i] = 0;
   if(s->maxPeakNr < peaknr) s->maxPeakNr = peaknr;
   color = c;
   intflag = iflag;
   strcpy(comment,com);
   maxAsslinks = 0;
   nrAsslinks = 0;
   fixed=0;
}

PeakInfo::PeakInfo(Spectrum *s)
{
   int i;

   s->maxPeakNr++;
   nr = s->maxPeakNr;
   t = PEAK;
   nrAss = 0;
   spec = s;
   vol = Wert(0.0);
   prob= 1.0;
   for(i=0;i<DIMENSION;i++) userAss[i] = 0;
   color = 1;
   comment[0] = '\0';
   intflag = '-';
   maxAsslinks = 0;
   nrAsslinks = 0;
   fixed=0;
}

void PeakInfo::establishAsslink(ErPeak *erPeak)
{
   ErPeak **newlinks;
   int i;

   if(nrAsslinks>=maxAsslinks) {
      newlinks = new ErPeak *[maxAsslinks+5];
      for(i=0;i<nrAsslinks;i++) newlinks[i] = asslink[i];
      if(maxAsslinks>0) delete asslink;
      maxAsslinks += 5;
      asslink = newlinks;
   }
   /* put newest assignment into asslink[0] */
   if(nrAsslinks>0) {
      asslink[nrAsslinks] = asslink[0];
      asslink[nrAsslinks]->asslinkindex=nrAsslinks;
   }
   asslink[0] = erPeak;
   erPeak->asslinkindex = 0;
   nrAsslinks++;

   /* debugging */
   //for(i=0;i<nrAsslinks-1;i++) {
   //   if(asslink[i]==erPeak) {
   //      printf("found  "); erPeak->print(0);
   //   }
   //}
}

void PeakInfo::retractAsslink(int index)
{
   nrAsslinks--;
   //printf("retracting: %d ",nr); asslink[index]->print(0);
   if(index != nrAsslinks) {
      asslink[nrAsslinks]->asslinkindex = index;
      asslink[index] = asslink[nrAsslinks];
   }
}

int PeakInfo::nrIntraDegen(ErPeak *erPeak)
{
   int i,j;
   int cnt=0;
   int res;

   /* count how many assignments which are equal in dim-1 dimensions */
   for(i=0;i<nrAsslinks;i++) {
      if(asslink[i] == erPeak) continue;
      res=0;
      for(j=0;j<spec->type->dim;j++) {
         if(asslink[i]->liCo[j] == erPeak->liCo[j] ) res++;
      }
      if(res==spec->type->dim-1)cnt++;
   }
   return(cnt);
}

int PeakInfo::getAmbigiuety(ErPeak *erPeak)
{
   int i,j;
   int ambig;
   int ident,equiv;

   ambig = INDISTINGUISHABLE;
   for(i=0;i<nrAsslinks;i++) {
      ident=0; equiv=0;
      for(j=0;j<spec->type->dim;j++) {
         if(asslink[i]->liCo[j] == erPeak->liCo[j] ) {
            ident++;
         } else if(asslink[i]->liCo[j]->floatIdentical(erPeak->liCo[j])) {
            equiv++;
         }
      }
      if(ident+equiv<spec->type->dim) {
         if(ident == spec->type->dim - 1 && ambig != DISTINGUISHABLE) {
            ambig = AMBIGOUS;
         } else {
            ambig = DISTINGUISHABLE;
         }
      }
   }
   return(ambig);
}

Peak::Peak(PeakInfo *inf, float *w_) : p(inf)
{
   int i;
   for(i=0;i<p->spec->type->dim;i++) {
      w[i] = w_[i];
      fold[i] = 0;
   }
}

Peak::Peak(PeakInfo *inf) : p(inf)
{
   int i;
   for(i=0;i<p->spec->type->dim;i++) {
      fold[i] = 0;
   }
}

void Peak::print(int off)
{
   Log log;
   char str[MAXLINE],s2[MAXLINE];
   int i;
   int ass=0;

   sprintf(str,"%*s%-10s peak %4d at",
      off,"",p->spec->type->name,p->nr);
   log.w(str);
   for(i=0;i<p->spec->type->dim;i++) {
      sprintf(str," %7.3f",w[i]);
      log.w(str);
   }
   for(i=0;i<p->spec->type->dim;i++) if(p->userAss[i]) ass=1;
   if(ass) {
      log.w(" assigned to ");
      for(i=0;i<p->spec->type->dim;i++) {
         if(p->userAss[i])
            p->userAss[i]->shortNameGet(str);
         else
           strcpy(str,"unass.");
         sprintf(s2,"%-8s",str);
         log.w(s2);
         if(i!=p->spec->type->dim-1) log.w("/ ");
      }
   }
   log.w("\n");
}

void Peak::write(PeakWriter &writer, int type, int &nrUnass, int &nrDegen,
		 int &nrUserAss)
{
   int i,j;
   ErPeak *erPeak,*cur;
   int cntAss = 0;
   int isUserAss = 0;

   writer.v = p->vol;
   writer.nrPeak = p->nr;
   writer.color = p->color;
   writer.intflag = p->intflag;
   strcpy(writer.comment,p->comment);
   for(i=0; i<p->spec->type->dim; i++) writer.w[i] = w[i];
   for(i=0; i<p->spec->type->dim; i++) if(p->userAss[i] != 0) isUserAss=1;
   if(!isUserAss) {
      erPeak = p->ass.start(j);
      while(erPeak != 0) {
         if(erPeak->type == type) {
            cntAss++;
            cur = erPeak;
         }
         erPeak = p->ass.next(j);
      }
      if(cntAss == 0) {
         for(i=0; i<p->spec->type->dim; i++) {
            writer.coAtNam[i][0] = '\0';
            writer.coResNum[i] = 0;
         }
         nrUnass++;
      } else {

         for(i=0; i<p->spec->type->dim; i++) {
            if(cur->getCo(i)->statPseudo==PSEUDON)
               cur->getCo(i)->pseudCo->
                     nameGet(writer.coAtNam[i], writer.coResNum[i]);
            else
               cur->getCo(i)->
                     nameGet(writer.coAtNam[i], writer.coResNum[i]);
         }
         if(cntAss>1) { nrDegen++; }
      }
   } else {
      nrUserAss++;
      for(i=0; i<p->spec->type->dim; i++) {
          if(p->userAss[i] == 0) {
             writer.coAtNam[i][0] = '\0';
             writer.coResNum[i] = 0;
          } else {
             p->userAss[i]->nameGet(writer.coAtNam[i], writer.coResNum[i]);
          }
      }
   }
   writer.write();
}

void Peak::writeAssPos(AssWriter &writer, int type, int &nrUnass, int &nrDegen)
{
   int i,j;
   ErPeak *erPeak;
   int cntAss = 0;

   writer.nrPeak = p->nr;
   erPeak = p->ass.start(j);
   while(erPeak != 0) {
      if(erPeak->type == type && erPeak->active) {
         cntAss++;
      }
      erPeak = p->ass.next(j);
   }
   if(cntAss != 0) {
      if(cntAss>1) nrDegen++;
      writer.startPeak();
      erPeak = p->ass.start(j);
      while(erPeak != 0) {
         if(erPeak->type == type && erPeak->active) {
            for(i=0; i<p->spec->type->dim; i++) {
               if(erPeak->getCo(i)->statPseudo==PSEUDON)
                  erPeak->getCo(i)->pseudCo->
                     nameGet(writer.coAtNam[i], writer.coResNum[i]);
               else
                  erPeak->getCo(i)->
                     nameGet(writer.coAtNam[i], writer.coResNum[i]);
            }
            writer.writeAssPos();
         }
         erPeak = p->ass.next(j);
      }
   } else {
      nrUnass++;
   }
}

void Peak::setP(float prob)
{
   p->p = prob;
}

void Peak::addP(float prob)
{
   if(prob>p->p) p->p = prob;
}

float Peak::getP()
{
   return(p->p);
}

float Peak::calcPSelect(Wert *ass,Wert *exp,int *)
{
   float sumP,prob;
   float x,v,value;
   int i;
   int index;

   sumP = 1.0;
   for(i=0;i<p->spec->type->dim;i++) {
      value = w[i];
      if(ass[i].x != BAD_PPM) {
         x=ass[i].x; v=ass[i].v; 
         index = (int)floor(fabs((x-value)/v) * RESOLUTION);
         if(index > RESOLUTION) index = RESOLUTION;
         prob = assDistr[index];
      } else {
         prob = expDistr(exp[i],value);
      }
      sumP *= prob;
   }
   sumP *= p->prob;
   return(sumP);
}


void Assignment::addAss(ErPeak *newAss)
{
   int i;
   int res = 0;
   ErPeak **newList;

   for(i=0;i<nrAss && !res;i++) if(newAss == ass[i]) res = 1;
   if(!res) {
      if(nrAss>=maxAss) {
         newList = new ErPeak*[maxAss+10];
         for(i=0;i<maxAss;i++) newList[i] = ass[i];
         if(maxAss>0) delete ass;
         ass = newList;
         maxAss += 10;
      }
      ass[nrAss] = newAss;
      nrAss++;
   }
}

void Assignment::removeAss(int type)
{
   int to=0, from=0;
   while(from<nrAss) {
      if(ass[from]->type == type) {
         from++;
      } else {
         if(from != to) ass[to]=ass[from];
         from++;
         to++;
      }
   }
   nrAss = nrAss + to - from;
}

void Assignment::resetAss()
{
   nrAss = 0;
}

//void Assignment::writeAss(AssWriter &writer, Spectrum *spec)
//{
//   int i;
//   for(i=0;i<nrAss;i++) {
//      ass[i]->writeAssPos(writer,spec);
//   }
//}

Assignment::~Assignment()
{
   if(ass!=0) delete ass;
}

ErPeak *Assignment::start(int &cur)
{
   cur = 0;
   return next(cur);
}

ErPeak *Assignment::next(int &cur)
{
   ErPeak *res;

   if(cur<nrAss) {
      res = ass[cur];
      cur++;
   } else {
      res = 0;
   }
   return res;
}


/***************************************************************************/
/* Data structure to search in n dimensional peaklists                     */
/* Partititioning of peaks into groubs of peaks that lie on the same side  */
/* of a limit                                                              */
/***************************************************************************/
 
void NdLiPeak::init(int n,int d)
{
   if(max != 0) {
      delete li;
      delete liLim;
      delete stack;
   }
   max=n;
   dim=d;
   nrStack=0;
   nr=0;
   if(max != 0) {
      li = new Peak*[max];
      liLim = new Limit[max];
      stack = new Sublist[(int)log(max) * 2 + 1];
   }
}

Peak *NdLiPeak::start(float *l, float *u)
{
   low = l; up = u;
   stack[0] = Sublist(0,nr);
   nrStack++;
   return next();
}

Peak *NdLiPeak::next()
{
   Peak *res;
   res = 0;
   while(res==0 && nrStack>0) {
      nrStack--;
      res=search(stack[nrStack]);
   }
   return(res);
}

Peak *NdLiPeak::search(Sublist sub)
{
   int i;
   int g;			// splitting point
   int d;			// splitting dimension
   float lim;			// splitting limit
   int res=-1;
   Peak *peak;

   while(res==-1) {
      if(sub.bis-sub.von <= 1) {
         res=1;
         for(i=0;i<dim;i++)
            if(li[sub.von]->w[i] < low[i] ||
               li[sub.von]->w[i] >  up[i]) res=0;
         if(res) peak=li[sub.von];
         else peak=0;
      } else {
         g = (sub.von+sub.bis)/2;
         d = liLim[g].dim;
         lim=liLim[g].lim;
         if(lim>=low[d] && lim<=up[d]) {
            stack[nrStack] = Sublist(sub.von,g);
            nrStack++;
            sub.von = g;
         } else if(lim<low[d]) {
            sub.von = g;
         } else {
            sub.bis = g;
         }
      }
   }
   return peak;
}

void NdLiPeak::add(Peak *peak)
{
   int i;

   for(i=0;i<dim;i++) {
      if(nr==0) {
         maxshift[i]=peak->w[i];
         minshift[i]=peak->w[i];
      } else {
         if(peak->w[i] > maxshift[i]) maxshift[i]=peak->w[i];
         if(peak->w[i] < minshift[i]) minshift[i]=peak->w[i];
      }
   }
   li[nr]=peak;
   nr++;
}

void NdLiPeak::sort()
{
   recSort(Sublist(0,nr),0);
}

void NdLiPeak::recSort(const Sublist sub,int d)
{
   int g;
   if(sub.bis - sub.von >1) {
      quicksort(Sublist(sub.von,sub.bis-1),d);
      g=(sub.von+sub.bis)/2;
      liLim[g].lim = (li[g]->w[d] + li[g-1]->w[d]) / 2;
      liLim[g].dim = d;
      d++; if(d>=dim) d=0;
      recSort(Sublist(sub.von,g),d);
      recSort(Sublist(g,sub.bis),d);
   }
}

void NdLiPeak::quicksort(const Sublist sub,int d)
{
   int i,j;
   Peak *x;
   Peak *w;

   i=sub.von; j=sub.bis;
   x=li[(i+j)/2];
   do {
      while(li[i]->w[d] < x->w[d]) i++;
      while(x->w[d] < li[j]->w[d]) j--;
      if(i<=j) {
         w=li[i]; li[i]=li[j]; li[j]=w;
         i++; j--;
      }
   } while(i<=j);
   if(sub.von < j) quicksort(Sublist(sub.von,j),d);
   if(i < sub.bis) quicksort(Sublist(i,sub.bis),d);
}

Peak *NdLiPeak::selectUnusedClose(Wert *w) 
{
   float shift[DIMENSION];
   int i;

   for(i=0;i<dim;i++) shift[i] = w[i].random(5.0);
   return(searchUnusedClose(shift,Sublist(0,nr)));
}

Peak *NdLiPeak::selectUnusedClose(Wert *w, float *l, float *u) 
{
   float shift[DIMENSION];
   int i;

   low = l; up = u;
   for(i=0;i<dim;i++) shift[i] = w[i].random(1.0);

//   printf("\nsearching: ");
//   for(i=0;i<dim;i++) printf("%5.2f %5.2f %5.2f /",shift[i],l[i],u[i]);  
//   printf("\n");

   return(searchRangeClose(shift,Sublist(0,nr)));
}

Peak *NdLiPeak::searchRangeClose(float *w,Sublist sub)
{
   int g;			// splitting point
   int d;			// splitting dimension
   float lim;			// splitting limit
   Peak *peak=0;
   int i;

   if(sub.bis-sub.von == 1) {
      peak=li[sub.von];
      if(peak->p->nrAss > 0 || peak->p->fixed) peak = 0;
      for(i=0;i<dim && peak;i++)
         if(peak->w[i]<low[i] || peak->w[i]>up[i]) peak=0;
//      if(peak) {
//         printf("found: ");
//         for(i=0;i<dim;i++) printf("%5.2f /",peak->w[i]);  
//         printf("\n");
//      }
   } else {
      g = (sub.von+sub.bis)/2;
      d = liLim[g].dim;
      lim=liLim[g].lim;
//      printf("%1d %5.2f\n",d,lim);
      if(lim<w[d]) {
         if(lim<up[d])
            peak = searchRangeClose(w,Sublist(g,sub.bis));
         if(peak==0 && lim>low[d])
            peak = searchRangeClose(w,Sublist(sub.von,g));
      } else {
         if(lim>low[d])
            peak = searchRangeClose(w,Sublist(sub.von,g));
         if(peak == 0 && lim<up[d])
            peak = searchRangeClose(w,Sublist(g,sub.bis));
      }
   }
   return peak;
}

Peak *NdLiPeak::searchUnusedClose(float *w,Sublist sub)
{
   int g;			// splitting point
   int d;			// splitting dimension
   float lim;			// splitting limit
   Peak *peak;

   if(sub.bis-sub.von == 1) {
      peak=li[sub.von];
      if(peak->p->nrAss > 0) peak = 0;
   } else {
      g = (sub.von+sub.bis)/2;
      d = liLim[g].dim;
      lim=liLim[g].lim;
      if(lim<w[d]) {
         peak = searchUnusedClose(w,Sublist(g,sub.bis));
         if(peak == 0) peak = searchUnusedClose(w,Sublist(sub.von,g));
      } else {
         peak = searchUnusedClose(w,Sublist(sub.von,g));
         if(peak == 0) peak = searchUnusedClose(w,Sublist(g,sub.bis));
      }
   }
   return peak;
}


int LiPeakDim::compare(const Element *p1, const Element *p2)
{
   int res;
   float d;
   int i;

   d =  ((Peak *)p1)->w[dim] - ((Peak *)p2)->w[dim];
   res = ((Peak *)p2)->p->nr - ((Peak *)p1)->p->nr;
   if(d>0.0) res=1;
   else if(d<0.0) res=-1;
   if(res==0) res = ((Peak *)p2)->p->spec->nr - ((Peak *)p1)->p->spec->nr;
   i = 0;
   while(res == 0 && i < ((Peak *)p2)->p->spec->type->dim) {
      res = ((Peak *)p2)->fold[i] - ((Peak *)p1)->fold[i];
      i++;
   }
   return(res);
}

Peak *LiPeakDim::startRange(float low, float up)
					    /* find first peak bigger than */
					    /* low, call next              */
{
   int i,j;

   limit=up;
   if(n==0) return(0);

   i=0;
   j=n-1;
   cur=0;
   while(i<=j) {
      cur = (i+j)/2;
      // printf("%d",cur); ((Peak *)l[cur])->print(0);
      if(((Peak *)l[cur])->w[dim] < low) i=cur+1;
      else j=cur-1;
   }
   if(((Peak *)l[cur])->w[dim] < low) cur++;
   return nextRange();
}

Peak *LiPeakDim::nextRange()
{
   Peak *p=0;

   if(cur<n && ((Peak *)l[cur])->w[dim] < limit) {
      // printf("%d",cur); ((Peak *)l[cur])->print(0);
      p=(Peak *)l[cur];
      cur++;
   }
   return p;
}