/*
************************************************************************
*
*   nmrdia.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.nmrdia.C
*   SCCS identification       : 1.6
*
************************************************************************
*/
/**************************************************************************/
/*                           nmrdia.cc                                    */
/*                                                                        */
/* diagram representing nmr experiments and assignments                   */
/**************************************************************************/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "element.h"
#include "peak.h"
#include "wert.h"
#include "atom.h"
#include "co.h"
#include "global.h"
#include "log.h"
#include "bib.h"
#include "envir.h"
#include "specTypes.h"
#include "map.h"
#include "parser.h"
#include "nmrdia.h"
#include "functions.h"

void NMRdia::peakMake(Envir &)
{
   NMRTopol T(&liCo, &liErPeak);
   Spectrum *spec;
   int res;
   int number=1;
   Log wlog;
   char str[MAXLINE];

   spec = usedSpec.start();
   while(spec != 0) {
      sprintf(str,"\n... %s peaks generated: ",spec->type->name);
      wlog.w(str);
      res=T.start(spec);
      while(res) {
         number += T.match(&ss);
         T.cleanup();
         res=T.next();
      }
      sprintf(str,"%d peaks, %d coherences",liErPeak.getn()
         ,liCo.getn());
      wlog.w(str);
      spec = usedSpec.next();
   }
}

SQCo *newRs(LiCo &, Atom *, float);

void NMRdia::coMake()
{
   Atom *a,*pseudat;
   Coherence *pseudco,*co;
   Coherence *dummy;
   Iterator i;
   Log wlog;
   char str[MAXLINE];
   PseudRequest pR;
   int cnt;

   wlog.w("\n... coherences generated: ");
   a = (Atom *)ss.Atoms.start(i);
   while(a != 0) {
      if(strcmp(a->atype,"PSEUD") == 0 ||
         strcmp(a->atype,"H_ALI") == 0 ||
         strcmp(a->atype,"H_AMI") == 0 ||
         strcmp(a->atype,"H_ARO") == 0 ||
         strcmp(a->atype,"H_SUL") == 0 ||
         strcmp(a->atype,"H_OXY") == 0 ||
         strcmp(a->atype,"C_ALI") == 0 ||
         strcmp(a->atype,"C_BYL") == 0 ||
         strcmp(a->atype,"C_ARO") == 0 ||
         strcmp(a->atype,"C_VIN") == 0 ||
         strcmp(a->atype,"N_AMI") == 0 ||
         strcmp(a->atype,"N_AMO") == 0 ||
         strcmp(a->atype,"P_ALI") == 0) {
         dummy = newRs(liCo,a,a->p);
      }
      a = (Atom *)ss.Atoms.next(i);
   }
   sprintf(str,"%d coherences",liCo.getn());
   wlog.w(str);

   cnt=0;
   pR.start();
   while(!pR.eof) {
      pseudat=(Atom *)ss.Atoms.find(Atom(pR.pseudnr));
      a      =(Atom *)ss.Atoms.find(Atom(pR.atnr));
      if(pseudat != 0 && a !=0) {
         pseudco=pseudat->getCo();
         co     =a->getCo();
         if(co!=0 && pseudco != 0 && pseudco!=co) {
            pseudco->addCoToPseudo(co);
            cnt++;
         }
      }
      pR.read();
   }
   sprintf(str,"\n... %d pseudo coherences defined",cnt);
   wlog.w(str);
}

void NMRdia::listPeaks(char *t)
{
   ErPeak *peak;
   Log wlog;
   Iterator i;

   wlog.w("\n   List of expected peaks:\n");
   peak = (ErPeak *)liErPeak.Index::startPars(i,t);
   while(peak != 0) {
      peak->print(5);
      peak = (ErPeak *)liErPeak.Index::nextPars(i,t);
   }
}

void NMRdia::setPeakImp(char *t,float imp)
{
   ErPeak *peak;
   Log wlog;
   char str[MAXLINE];
   Iterator i;
   int cnt=0;

   peak = (ErPeak *)liErPeak.Index::startPars(i,t);
   while(peak != 0) {
      peak->importance = imp;
      cnt++;
      peak = (ErPeak *)liErPeak.Index::nextPars(i,t);
   }
   sprintf(str,"\n   importance set for %d peaks",cnt);
   wlog.w(str);
}

void NMRdia::initOptPos(int number)
{
   nrActive = 0;
   optPos.init(number,liErPeak,liCo,peaks);

   calcScoreTables();
}

int NMRdia::writeCo(char *name, int type)
{
   Coherence *co;
   char nameCo[MAXNAME];
   int nrAa;
   CoWriter writer;
   int res=1;
   Log wlog;
   char str[MAXLINE];
   int count;

   writer.start(name);
   co = (Coherence *)liCo.start();
   count = 0;
   while(co!=0 && writer.eof==0) {
      if(co->firstAt()->frag->type == type) {
         count++;
         co->nameGet(nameCo,nrAa);
         writer.name = nameCo;
         writer.nrAa = nrAa;
         writer.shift = co->meanF;
         if(!co->meanF.valid()) writer.shift.x = BAD_PPM;
         writer.write();
      }
      co = (Coherence *)liCo.next();
   }
   if(writer.eof == -1) {
      res = 0;
   }
   writer.close();
   if(res) {
      sprintf(str,"\n... %d coherences written to file '%s'",
            count,name);
      wlog.w(str);
   } else {
      sprintf(str,"\n... ERROR: could not write coherence file '%s'",
             name);
      wlog.w(str);
   }
   return(res);
}

void NMRdia::loadFreqCo(char *name, int typeFrag)
{
   BibAccess bib;
   FreqRequest fr;
   Coherence *co,*pseud;
   int nrSet;
   Log wlog;
   char str[MAXLINE];

   nrSet = 0;
   if(bib.coLoad(name,typeFrag)) {
      fr.start();
      while(fr.eof == 0) {
         if(fr.shift.x != BAD_PPM) {
            co = liCo.find(fr.nrExt,fr.name,typeFrag);
            if(co != 0) {
               nrSet++;
               co->assF = fr.shift;
               if(co->statPseudo == PSEUDOCO) {
       		  co->turnOnPseud();
                  pseud=(Coherence *)co->pseudo->start();
                  while(pseud != 0) {
                     pseud->assF = fr.shift;
                     pseud=(Coherence *)co->pseudo->next();
                  }
               } else if(co->statPseudo == PSEUDON) {
		  co->pseudCo->turnOffPseud();
               }
            } else {
               sprintf(str,
                 "\n... WARNING: coherence '%s(%d)' not present",
                  fr.name, fr.nrExt);
               wlog.warn(10,str);
            }
         }
         fr.read();
      }
   } else {
      sprintf(str,"\n... atom list '%s' could not be read",name);
      wlog.warn(-1,str);
   }

   sprintf(str,"\n... %d coherence frequencies set",nrSet);
   wlog.w(str);
}

void NMRdia::setFreqTol(char *name, float value)
{
   Coherence *co;
   int nrSet;
   Log wlog;
   char str[MAXLINE];
   Iterator i;

   nrSet = 0;
   co = (Coherence *)liCo.startPars(i,name);
   while(co != 0) {
      if(co->assF.x != BAD_PPM) {
         co->assF.v = value;
         nrSet++;
      }
      co = (Coherence *)liCo.nextPars(i,name);
   }
   sprintf(str,"\n... %d tolerances of frequencies set",nrSet);
   wlog.w(str);
}

void NMRdia::defSecShifts(char *fn)
{
   FILE *fp;
   Coherence *co;
   int nrSet,resNum,warned;
   Log log;
   char str[MAXLINE],resNam[MAXLINE],secNam[MAXLINE],s1[MAXLINE];
   char coName[MAXLINE];
   SecShiftRequest shiftRequest;

   nrSet = 0;
   fp=fopen(fn,"r");
   if(fp == NULL) {
      sprintf(str, "\n ... can not open file %s\n",fn);
      log.w(str);
   } else {
      while(!feof(fp)) {
         fgets(s1,MAXLINE,fp);
         if(sscanf(s1,"%d%s%s",&resNum,resNam,secNam)!=3) {
           log.w(" ... could not parse line\n");
           sprintf(str, "     %s",s1);
           log.w(str);
         } else {
           warned=0;
	   co = (Coherence *)liCo.start();
	   while(co != 0) {
	      if(co->firstAt()->frag->nrExt==resNum) {
		 if(strcmp(co->firstAt()->frag->name, resNam)!=0 && !warned) {
		    sprintf(str, "... residue %d has different name (%s and %s)\n", 
		            resNum, co->firstAt()->frag->name, resNam);
                    log.w(str);
                    warned=1;
		 } 
	         if(co->firstAt()->frag->type == DESTFRAG &&
                    shiftRequest.get(co->firstAt()->frag->nrInt,
                                     co->firstAt()->nr,
                                     secNam))
                 {
                     co->assF=shiftRequest.shift;
                     co->fullNameGet(coName);
                     sprintf(str,"%12s expected shift %7.3f +/- %5.3f\n",
                             coName,co->assF.x,co->assF.v);
                     log.w(str);
                     nrSet++;
		 }
	      }
	      co = (Coherence *)liCo.next();
	   }
         }
      }
   }
   sprintf(str,"\n... %d frequencies set\n",nrSet);
   log.w(str);
}

void NMRdia::writeExpPeaks(Spectrum *spec, char *pname,
                         char *coname, int typeFrag, int maxUndef)
{
   BibAccess bib;
   ErPeak *peak;
   PeakWriter writer;
   int i;
   Iterator it;
   int nrPeak = 1;
   Log wlog;
   char str[MAXLINE];

   if(bib.coLoad(coname,typeFrag)) {
      writer.dim = spec->type->dim;
      strcpy(writer.comment,"");
      for(i=0;i<writer.dim;i++) writer.perm[i] = spec->type->permut[i];
      writer.start(pname);
      peak = (ErPeak *)liErPeak.start(it,spec);
      while(peak != 0 && writer.eof == 0) {
         if(peak->active) {
            writer.nrPeak = nrPeak;
            if(peak->write(writer,typeFrag,maxUndef)) nrPeak++;
         }
         peak = (ErPeak *)liErPeak.next(it,spec);
      }
      writer.close();
   } else {
      sprintf(str,"\n... atom list '%s' could not be read",coname);
      wlog.warn(-1,str);
   }

   sprintf(str,"\n... %d peaks written to file '%s'",
       nrPeak-1,pname);
   wlog.w(str);
}

void NMRdia::reportMissPeaks(Spectrum *spec, int typeFrag, int undef)
{
   Coherence *co;
   ErPeak *erPeak;
   Peak *peak;
   PeakInfo newInf=PeakInfo(spec);
   Peak     newP  =Peak(&newInf);
   Iterator it;
   Wert freq[DIMENSION];
   Log wlog;
   int i;
   LiPeakName liSort = LiPeakName(1024);

   peak = (Peak *)peaks.start();
   while(peak != 0) {
      if(peak->p->spec == spec) {
         liSort.add(peak);
      }
      peak = (Peak *)peaks.next();
   }
   liSort.sort(-1);

   wlog.w("\n... expected peaks not present in list of measured peaks:\n");
   erPeak = (ErPeak *)liErPeak.start(it,spec);
   while(erPeak != 0) {
      if(erPeak->getWriteFreq(freq) == undef && erPeak->type == typeFrag) {
         for(i=0;i<erPeak->spec->type->dim;i++) {
            co=erPeak->getCo(i);
            if(co->statPseudo==PSEUDON) co=co->pseudCo;
            newP.p->userAss[i] = co;
         }
         if(liSort.find(newP) == 0) {
            erPeak->print(3);
         }
      }
      erPeak = (ErPeak *)liErPeak.next(it,spec);
   }
}

void NMRdia::appendExpPeaks(Spectrum *spec, int typeFrag, int undef,
                            int color)
{
   Coherence *co;
   ErPeak *erPeak;
   Peak *peak;
   Peak *newP;
   PeakInfo *newInf;
   Iterator it;
   float w[DIMENSION];
   Wert freq[DIMENSION];
   int added = 0;
   int nrPeaks = 0;
   int i;
   Log wlog;
   char str[MAXLINE];
   LiPeakName liSort = LiPeakName(1024);
   LiPeakName liNew  = LiPeakName(1024);

   peak = (Peak *)peaks.start();
   while(peak != 0) {
      if(peak->p->spec == spec) {
         liSort.add(peak);
         nrPeaks++;
      }
      peak = (Peak *)peaks.next();
   }
   liSort.sort(-1);

   erPeak = (ErPeak *)liErPeak.start(it,spec);
   while(erPeak != 0) {
      if(erPeak->getWriteFreq(freq) == undef && erPeak->type == typeFrag) {
         for(i=0;i<DIMENSION;i++) w[i] = freq[i].x;
         newInf = new PeakInfo(spec);
         newP   = new Peak(newInf,w);
         for(i=0;i<erPeak->spec->type->dim;i++) {
            co=erPeak->getCo(i);
            if(co->statPseudo==PSEUDON) co=co->pseudCo;
            newP->p->userAss[i] = co;
         }
         newP->p->color = color;
         if(liSort.find(*newP) != 0 || !liNew.insert(newP)) {
            delete newP;
            delete newInf;
         }
      }
      erPeak = (ErPeak *)liErPeak.next(it,spec);
   }

   peak = (Peak *)liNew.start();
   while(peak != 0) {
      peaks.add(peak);
      added++;
      peak = (Peak *)liNew.next();
   }
   peaks.sort(-1);

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

   peak = (Peak *)peaks.start();
   while(peak != 0) {
      if(peak->p->spec == spec) spec->peaks.add(peak);
      peak = (Peak *)peaks.next();
   }
   spec->peaks.sort();

   sprintf(str,"\n... %d peaks appended to %d measured peaks",
       added,nrPeaks);
   wlog.w(str);
}

void NMRdia::loadPeakAss(LiPeak &peaks, Spectrum *sp, char *an, char *cn, int t,
			 float tol)
{
   BibAccess bib;
   AssRequest aR;
   ErPeak *erPeak,*pseudPeak,*foundPeak;
   Peak *peak;
   Log wlog;
   char str[MAXLINE];
   char coName[MAXLINE];
   int count =0;
   int i,ipeak;
   int nrUnass;
   int ok,res,nrmatch;
   int missCo;
   int nopeak;
   int presentCo;
   Coherence *co1,*co2;
   Coherence *userAss[DIMENSION];
   float low[DIMENSION],up[DIMENSION];
   int flag[DIMENSION];
   int cntmatches,match[3];

   for(i=0;i<DIMENSION;i++) flag[i]=ADAPT;
   for(i=0;i<3;i++) match[i]=0;
   if(bib.coLoad(cn,t)) {
      aR.start(an,sp->type->dim,sp->type->noise,sp->type->permut);
      while(!aR.eof) {
         nrUnass = 0;
         missCo = -1;
         ok = 1;
         nopeak=0;
         /* check coherences */
         for(i=0;i<sp->type->dim;i++) {
            if(aR.nrFrag[i] == -1) {
               nrUnass++;
               userAss[i] = 0;
            } else {
               userAss[i] = liCo.find(aR.nrFrag[i],aR.getName(i),t);
               if(userAss[i] == 0) missCo = i;
               else {
                  presentCo=i;
                  co1=userAss[i];
       		  if(co1->statPseudo == PSEUDOCO) {
                     co1->turnOnPseud();
                  } else if(co1->statPseudo == PSEUDON) {
                     co1->pseudCo->turnOffPseud();
                  }
               }
            }
         }
         /* set user assignment */
         if(tol<=0.0) {
            peak = peaks.find(aR.nr,sp);
            if(peak == 0) nopeak=1;
         } else {
            for(i=0;i<sp->type->dim;i++) {
               low[i]= aR.w[i] - sp->type->accPick[i]*tol;
               up[i] = aR.w[i] + sp->type->accPick[i]*tol;
            }
            peak = sp->peaks.start(low,up);
            if(peak == 0) nopeak=2;
         }
         cntmatches=0;
         while(peak!=0) {
            cntmatches++;
            if(t == DESTFRAG) {
               for(i=0;i<sp->type->dim;i++)
                  peak->p->userAss[i] = userAss[i];
            }
            if(tol<=0.0) peak=0;
            else peak = sp->peaks.next();
         }
         if(cntmatches>2) cntmatches=2;
         match[cntmatches]++;

         /* search for matching expected peaks */
         if(nrUnass <= sp->type->dim-2 && missCo == -1) {
            ok=0;
            co1 = userAss[presentCo];
            if(co1->statPseudo==PSEUDOCO) {
               co2=(Coherence *)co1->pseudo->start();
            } else co2=co1;
            while(co2 != 0) {
               erPeak = co2->startPeak(ipeak); res=0; nrmatch=0;
               while(erPeak != 0 && !res) {
                  if(erPeak->spec==sp && erPeak->type==t) {res=1;}
                  else {res=0;}
                  for(i=0;i<sp->type->dim && res;i++) {
                     if(userAss[i] != 0 &&
                        !userAss[i]->
                           pseudIdentical(erPeak->getCo(i))) res=0;
                  }
                  if(res) {
                     nrmatch++;
                     foundPeak=erPeak;;
                  }
                  erPeak = co2->nextPeak(ipeak);
               }
               if(nrmatch==1) {   /* matching peak found, set assignments */
                  ok = 1;
                  erPeak=foundPeak;
                  if(tol<=0.0) {
                     peak = peaks.find(aR.nr,sp);
                     if(peak == 0) nopeak=1;
                  } else {
                     for(i=0;i<sp->type->dim;i++) {
                        low[i]= aR.w[i] - sp->type->accPick[i]*tol;
                        up[i] = aR.w[i] + sp->type->accPick[i]*tol;
                     }
                     peak = sp->peaks.start(low,up);
                     if(peak == 0) nopeak=2;
                  }
                  while(peak!=0) {
                     peak->p->ass.addAss(erPeak);
                     if(t == DESTFRAG) erPeak->makeAss(peak,flag,0,0);
                     count++;
                     if(erPeak->liPseud!=0 && erPeak->liPseud->getn()>0) {
                        pseudPeak=(ErPeak *)erPeak->liPseud->start();
                        while(pseudPeak != 0) {
                           peak->p->ass.addAss(pseudPeak);
                           if(t == DESTFRAG)
				pseudPeak->makeAss(peak,flag,0,0);
                           pseudPeak=(ErPeak *)erPeak->liPseud->next();
                        }
                     }
                     if(tol<=0.0) peak=0;
                     else peak = sp->peaks.next();
                  }
               }
               if(co1->statPseudo==PSEUDOCO) {
                  co2=(Coherence *)co1->pseudo->next();
               } else co2=0;
            }
         } /* end search for matching expected peaks */
         /* treat errors */
         if(missCo != -1) {
            sprintf(str,
               "\nWARNING: unknown coherence for peak %d: %s %d",
               aR.nr,aR.getName(missCo),aR.nrFrag[missCo]);
            wlog.warn(-1,str);
         } else if(!ok) {
            sprintf(str,
               "\nWARNING: unknown assignment for peak %d: ",aR.nr);
            wlog.warn(-1,str);
            for(i=0;i<sp->type->dim;i++) {
               if(userAss[i] == 0) {
                  wlog.warn(-1,"-");
               } else {
                  userAss[i]->fullNameGet(coName);
                  wlog.warn(-1,coName);
               }
               if(i<sp->type->dim - 1) wlog.warn(-1,"/");
            }
         } else if(nopeak==1) {
            sprintf(str,"\nWARNING: peak #%d not defined",aR.nr);
                    wlog.warn(10,str);
         } else if(nopeak==2) {
            sprintf(str,
               "\nWARNING: no peak present at position of peak #%d",
               aR.nr);
            wlog.warn(10,str);
         }
         aR.read();
      }
   } else {
      sprintf(str,"\n... atom list '%s' could not be read",cn);
      wlog.warn(-1,str);
   }

   sprintf(str,"\n... %d assignment possibilities read",count);
   wlog.w(str);
   sprintf(str,"\n       %d assignments with no corresponding peak",match[0]);
   wlog.w(str);
   sprintf(str,"\n       %d assignments with a unique peak",match[1]);
   wlog.w(str);
   sprintf(str,"\n       %d assignments with multiple peaks",match[2]);
   wlog.w(str);
}

void NMRdia::addPeakAss()
{
   Log wlog;
   char str[MAXLINE];
   int count =0;
   Iterator it;
   ErPeak *erPeak;
   Peak *ass;

   erPeak = (ErPeak *)liErPeak.start(it);
   while(erPeak != 0) {
      ass = erPeak->getAss();
      if(ass != 0) {
         ass->p->ass.addAss(erPeak);
         count++;
      }
      erPeak = (ErPeak *)liErPeak.next(it);
   }

   sprintf(str,"\n... %d assignment possibilities added",count);
   wlog.w(str);
}

void NMRdia::addPossAss()
{
   Log wlog;
   char str[MAXLINE];
   int count =0;
   Iterator it;
   ErPeak *erPeak;
   Peak *peak;
   int nrAss,i;
   float low[DIMENSION],up[DIMENSION];
   Spectrum *sp;
   Coherence *co;

   calcMeanFreq();
   erPeak = (ErPeak *)liErPeak.start(it);
   while(erPeak != 0) {
      sp=erPeak->spec;
      nrAss=0;
      for(i=0;i<sp->type->dim;i++) {
         co=erPeak->getCo(i);
         if(co->measF[sp->nr].valid()) {
            nrAss++;
            low[i]=co->measF[sp->nr].x - sp->type->accPick[i];
            up[i] =co->measF[sp->nr].x + sp->type->accPick[i];
         } else if(co->meanF.valid()) {
            nrAss++;
            low[i]=co->meanF.x - sp->type->accSpec[i];
            up[i] =co->meanF.x + sp->type->accSpec[i];
         } else if(co->assF.x != BAD_PPM) {
            nrAss++;
            low[i]=co->assF.x - sp->type->accSpec[i];
            up[i] =co->assF.x + sp->type->accSpec[i];
         }
      }
      if(nrAss==sp->type->dim) {
         peak = sp->peaks.start(low,up);
         while(peak!=0) {
            peak->p->ass.addAss(erPeak);
            count++;
            peak = sp->peaks.next();
         }
      }
      erPeak = (ErPeak *)liErPeak.next(it);
   }

   sprintf(str,"\n... %d assignment possibilities added",count);
   wlog.w(str);
}


void NMRdia::markPosPeaks(Spectrum *spec, int d1, int d2)
{
   ErPeak *erPeak;
   Peak *pos;
   int status[DIMENSION];
   float low[DIMENSION], up[DIMENSION];
   Log wlog;
   char str[MAXLINE];
   int count = 0;
   int i;
   float diff;

   d1--; d2--;
   if(d1>spec->type->dim || d2>spec->type->dim) {
      wlog.warn(-1,"\n... WARNING: invalid dimension specified");
      return;
   }

   pos = (Peak *)peaks.start();
   while(pos != 0) {
      pos->setP(0.0);
      pos = (Peak *)peaks.next();
   }

   erPeak = liErPeak.start();
   while(erPeak != 0) {
      if(erPeak->type == DESTFRAG && erPeak->spec == spec) {
         erPeak->statCo(status);
         if((d1<0 || status[d1] == ASSIGNED)&&(d2<0 || status[d2] == ASSIGNED)){
            erPeak->range(status,low,up);
            for(i=0;i<spec->type->dim;i++) {
               if(status[i] == UNASSIGNED) {
                  diff = (up[i] - low[i]) * 5;
                  low[i] -= diff;
                  up[i]  += diff;
               }
            }
            pos = spec->peaks.start(low,up);
            while(pos != 0) {
               if(pos->getP() == 0.0) count++;
               pos->setP(1.0);
               pos = spec->peaks.next();
            }
         }
      }
      erPeak = liErPeak.next();
   }

   sprintf(str,"\n... %d measured peaks marked",count);
   wlog.w(str);
}

void NMRdia::resetAss()
{
   ErPeak *erPeak;
   int count = 0;
   Log wlog;
   char str[MAXLINE];
   Coherence *co;

   erPeak = liErPeak.start();
   while(erPeak != 0) {
      if(erPeak->status != P_FIXED) {
         erPeak->resetAss();
         count++;
      }
      erPeak = liErPeak.next();
   }
   sprintf(str,"\n... %d peak assignments reset",count);
   wlog.w(str);

//   co = (Coherence *)liCo.start();
//   while(co != 0) {
//      co->maxScore=0.0;
//      co = (Coherence *)liCo.next();
//   }


   calcMeanFreq();
}

void NMRdia::resetAss(int size,int res1,int res2)
{
   ErPeak *erPeak;
   Frag *frag;
   int count = 0;
   Log wlog;
   char str[MAXLINE];
   int i,res,nrRes,equivNr;
   Iterator it;

   /* starting residue should have same spin system type */
   frag = (Frag *)ss.Frags.find(Frag(res1));
   if(frag != 0)  equivNr=frag->equivNr;
   else equivNr=0;

   if(equivNr != 0) {
      frag = (Frag *)ss.Frags.find(Frag(res2));
      frag = (Frag *)ss.Frags.start(it);
      while(frag != 0 && frag->nrInt != res2) frag=(Frag *)ss.Frags.next(it);
      while(frag != 0 && frag->equivNr != equivNr)
         frag=(Frag *)ss.Frags.nextCyc(it);
      if(frag) res2=frag->nrInt;
   }

   /* reset assignments */
   erPeak = liErPeak.start();
   while(erPeak != 0) {
      if(erPeak->type == DESTFRAG) {
         res=0;
         for(i=0;i<erPeak->spec->type->dim;i++) {
            nrRes=erPeak->liCo[i]->firstAt()->frag->nrInt;
            if((nrRes >= res1 && nrRes < res1+size) ||
               (nrRes >= res2 && nrRes < res2+size)) res=1;
         }
         if(res) {
            if(erPeak->status != P_FIXED) erPeak->resetAss();
            count++;
         }
      }
      erPeak = liErPeak.next();
   }
   sprintf(str,
       "\n... %d peak assignments reset for residues %d - %d and %d - %d",
       count,res1,res1+size,res2,res2+size);
   wlog.w(str);

   calcMeanFreq();
}

void NMRdia::calcMeanFreq()
{
   Coherence *co;
   Peak *peak;
   ErPeak *erPeak;
   Log wlog;
   char str[MAXLINE];
   int count = 0;

   co = (Coherence *)liCo.start();
   while(co != 0) {
      if(co->updateFreq()) count++;
      co = (Coherence *)liCo.next();
   }

   /* set number of assignments */
   peak = peaks.start();
   while(peak != 0) {
      peak->resetAss();
      peak = peaks.next();
   }
   erPeak = (ErPeak *)liErPeak.start();
   while(erPeak != 0) {
      peak = erPeak->getAss();
      if(peak != 0) {
         peak->p->nrAss++;
         peak->p->establishAsslink(erPeak);
      }
      erPeak = (ErPeak *)liErPeak.next();
   }

   sprintf(str,"\n... %d coherence frequencies calculated",count);
   wlog.w(str);
}

void NMRdia::checkIncon()
{
   Coherence *co;
   int ip,idim;
   int title,incSpec,incAll;
   ErPeak *erPeak;
   Peak *peak;
   Log wlog;
   char str[MAXLINE];
   char name[MAXLINE];

   co = (Coherence *)liCo.start();
   while(co != 0) {
      title=0;
      if(co->updateFreq()) {
         erPeak=co->startPeak(ip);
         while(erPeak != 0) {
            incSpec=0;incAll=0;
            peak=erPeak->getAss();
            if(peak) {
               for(idim=0;idim<erPeak->spec->type->dim &&
                   erPeak->getCo(idim)!=co;idim++);
               if(idim<erPeak->spec->type->dim) {
                   if(co->measF[erPeak->spec->nr].valid() &&
                      fabs(peak->w[idim] - co->measF[erPeak->spec->nr].x)>
                      erPeak->spec->type->accPick[idim]) incSpec=1;
                   if(fabs(peak->w[idim] - co->meanF.x)>
                      erPeak->spec->type->accSpec[idim]) incAll=1;
               }
            }
            if(incSpec || incAll) {
               if(!title) {
                  title=1;
                  co->fullNameGet(name);
                  sprintf(str,"\n%s at %7.3f has inconsistent assignments involving",
                     name,co->meanF.x);
                  wlog.w(str);
               }
               erPeak->name(name);
               if(incSpec)
                  sprintf(str,"\n   %4d %7.3f ppm,  %-45s",
                     peak->p->nr,peak->w[idim],name);
               else
                  sprintf(str,"\n   %4d %7.3f ppm,  %-45s (spec. ok)",peak->p->nr,peak->w[idim],name);
               wlog.w(str);
            }
            erPeak=co->nextPeak(ip);
         }
      }
      co = (Coherence *)liCo.next();
   }

}

void NMRdia::resetAssF()
{
   Coherence *co;
   Log wlog;
   char str[MAXLINE];
   int count = 0;

   co = (Coherence *)liCo.start();
   while(co != 0) {
      co->resetAssF();
      count++;
      co = (Coherence *)liCo.next();
   }

   sprintf(str,"\n... %d coherence frequencies reset",count);
   wlog.w(str);
}

void NMRdia::genSsEquiv()
{
   ErPeak *p1; Iterator i1;
   ErPeak *p2; Iterator i2;

   p1 = (ErPeak *)liErPeak.start(i1);
   p2 = (ErPeak *)liErPeak.start(i2);

   while(p1 != 0 && p2 != 0) {
      if(p1->type == DESTFRAG && p2->type == DESTFRAG && p1->spec == p2->spec &&
         (p2->equiv(p1) || p1->equiv(p2)) && p1!=p2) {
         p2->addEquiv(p1);
      }
      p1 = (ErPeak *)liErPeak.next(i1);
      if(p1 == 0) {
         p2 = (ErPeak *)liErPeak.next(i2);
         while(p2 != 0 && p2->type != DESTFRAG)
            p2 = (ErPeak *)liErPeak.next(i2);
         p1=p2; i1=i2;
      }
   }
}
   
void NMRdia::reportFrags(LiCrit &criteria)
{
   Log wlog;
   char str[MAXLINE];
   Crit *crit;
   Frag *frag;
   Score *score;
   int i;
   int maxlen;

   wlog.r("\n");
   maxlen = criteria.lenNam();
   for(i=0;i<maxlen;i++) {
      wlog.r("\n           ");
      crit = (Crit *)criteria.start();
      while(crit != 0) {
         sprintf(str," %c",crit->nam(i,maxlen));
         wlog.r(str);
         crit = (Crit *)criteria.next();
      }
   }

   crit = (Crit *)criteria.start();
   while(crit != 0) {
      crit->initSums();
      crit = (Crit *)criteria.next();
   }

   frag = (Frag *)ss.Frags.start();
   while(frag != 0) {
      if(frag->type == DESTFRAG) {
         sprintf(str,"\n%5s %-3d: ",frag->name,frag->nrExt);
         wlog.r(str);
         crit = (Crit *)criteria.start();
         score= (Score *)frag->scores.start();
         while(crit != 0) {
            sprintf(str," %c",crit->decision(score));
            crit->addSum(score);
            wlog.r(str);
            crit = (Crit *)criteria.next();
            score= (Score *)frag->scores.next();
         }
         wlog.r("   |");
         crit = (Crit *)criteria.start();
         score= (Score *)frag->scores.start();
         while(crit != 0) {
            if(crit->comment(score) && crit->writeCom) {
               sprintf(str," %s",crit->comment(score));
               wlog.r(str);
            }
            crit = (Crit *)criteria.next();
            score= (Score *)frag->scores.next();
         }
      }
      frag = (Frag *)ss.Frags.next();
   }

   wlog.r("\n\n");
   crit = (Crit *)criteria.start();
   while(crit != 0) {
      if(crit->printSum()) {
         sprintf(str,"%s\n",crit->printSum());
         wlog.r(str);
      }
      crit = (Crit *)criteria.next();
   }

}

void NMRdia::reportCos(LiCrit &criteria, char *templ)
{
   Log wlog;
   char str[MAXLINE];
   char name[MAXLINE];
   Crit *crit;
   Score *score;
   Coherence *co;
   int i;
   int maxlen;
   Iterator it;

   wlog.r("\n");
   maxlen = criteria.lenNam();
   for(i=0;i<maxlen;i++) {
      if(i==maxlen-1)
         sprintf(str,"\n%-16s%-8s%-8s%-8s%-8s%","name","freq",
         "score","mutInf","Signif");
//         sprintf(str,"\n%-16s%-8s%","name","bonus");
      else sprintf(str,"\n%64s","");
//      else sprintf(str,"\n%24s","");
      wlog.r(str);
      crit = (Crit *)criteria.start();
      while(crit != 0) {
         sprintf(str," %c",crit->nam(i,maxlen));
         wlog.r(str);
         crit = (Crit *)criteria.next();
      }
   }

   co = (Coherence *)liCo.startPars(it,templ);
   while(co != 0) {
      if(co->firstAt()->frag->type == DESTFRAG) {
         co->fullNameGet(name);
         sprintf(str,
            "\n%-15s %-7.3f %-7.3f %-7.1f %-7.1f:",
            name,co->meanF.x,co->localScore,
            co->mutInf,co->imp.sig);
//         sprintf(str,
//            "\n%-15s %-7.3f:",name,co->meanF.x);
         wlog.r(str);
         crit = (Crit *)criteria.start();
         score= (Score *)co->scores.start();
         while(crit != 0) {
            sprintf(str," %c",crit->decision(score));
            wlog.r(str);
            crit = (Crit *)criteria.next();
            score= (Score *)co->scores.next();
         }
         wlog.r("   |");
         crit = (Crit *)criteria.start();
         score= (Score *)co->scores.start();
         while(crit != 0) {
            if(crit->comment(score) && crit->writeCom) {
               sprintf(str," %s",crit->comment(score));
               wlog.r(str);
            }
            crit = (Crit *)criteria.next();
            score= (Score *)co->scores.next();
         }
      }
      co = (Coherence *)liCo.nextPars(it,templ);
   }
}

void NMRdia::reportCalc(LiCrit &critCo, LiCrit &critFrag, LiCrit &)
{
   Log wlog;
   Frag *frag;
   Coherence *co;

   frag = (Frag *)ss.Frags.start();
   while(frag != 0) {
      frag->scores.reset();
      frag = (Frag *)ss.Frags.next();
   }

   co = (Coherence *)liCo.start();
   while(co != 0) {
      co->scores.reset();
      co = (Coherence *)liCo.next();
   }

   scoreCos(critCo);
   wlog.w("\n... scoring for coherences calculated");
   scoreFrags(critFrag);
   wlog.w("\n... scoring for fragments calculated");
}

void NMRdia::scoreFrags(LiCrit &criteria)
{
   ErPeak *peak;
   Atom *atom;
   Iterator i,j,k;

   peak = (ErPeak *)liErPeak.start(i);
   while(peak != 0) {
      atom = peak->startAt(j,k);
      while(atom != 0 && peak->type == DESTFRAG) {
         criteria.upErPeak(atom->frag->scores,peak,atom->frag,0);
         atom = peak->nextAt(j,k);
      }
      peak = liErPeak.next(i);
   }
}

void NMRdia::scoreCos(LiCrit &criteria)
{
   ErPeak *peak;
   Coherence *co;
   Iterator i,j;

   peak = (ErPeak *)liErPeak.start(i);
   while(peak != 0) {
      co = peak->startCo(j);
      while(co != 0 && peak->type == DESTFRAG) {
         criteria.upErPeak(co->scores,peak,co->firstAt()->frag,co);
         co = peak->nextCo(j);
      }
      peak = liErPeak.next(i);
   }

   criteria.upCos(liCo);
}

void NMRdia::sumScore(int nr,float &s, float &v)
{
   Coherence *co;
   ErPeak *erPeak;
   float inf;

   v=0.0;
   if(nr<0) {
      s=0.0;
      co = (Coherence *)liCo.start();
      while(co != 0) {
         if(co->firstAt()->frag->type == DESTFRAG) {
            inf= co->mutInf;
            v += inf * co->imp.sig;
            s += inf;
         }
         co = (Coherence *)liCo.next();
      }

      erPeak = (ErPeak *)liErPeak.start();
      while(erPeak != 0) {
         if(erPeak->type == DESTFRAG) {
            inf = erPeak->mutInf;
            v += inf * erPeak->imp.sig;
            s += inf;
         }
         erPeak = (ErPeak *)liErPeak.next();
      }
      v /= (float)liCo.getn();
      s /= (float)liCo.getn();
   } else {
      co = (Coherence *)liCo.start();
      while(co != 0) {
         if(co->firstAt()->frag->type == DESTFRAG) {
            v += co->getMutInf(nr) * co->imp.sig;;
         }
         co = (Coherence *)liCo.next();
      }

      erPeak = (ErPeak *)liErPeak.start();
      while(erPeak != 0) {
         if(erPeak->type == DESTFRAG) {
            v += erPeak->getMutInf(nr) * erPeak->imp.sig;
         }
         erPeak = (ErPeak *)liErPeak.next();
      }
      v /= (float)liCo.getn();
   }

}


void NMRdia::scoreGlobal()
{
   Coherence *co;
   Log wrlog;
   char str[MAXLINE];
   ErPeak *erPeak;

   /* reset mutInf of peaks */
   erPeak = (ErPeak *)liErPeak.start();
   while(erPeak != 0) {
      if(erPeak->type == DESTFRAG) erPeak->mutInf = 0.0;
      erPeak = (ErPeak *)liErPeak.next();
   }

   /* calc. prob. of coherences */
   co = (Coherence *)liCo.start();
   while(co != 0) {
      co->updateScores();
      co->calcMeanScore();
      co = (Coherence *)liCo.next();
   }

   sumScore(-1,score,value);

   sprintf(str,
      "\n    TOTAL VALUE FOR SELECTION ... %8.3f",value);
   wrlog.w(str);
   sprintf(str,
      "\n    MUTUAL INFORMATION .......... %8.3f",score);
   wrlog.w(str);

   distHamilton();
}



#define CL_EXP  6
#define CL_MEAS 6

void NMRdia::assStatistic(Spectrum *spec)
{
   Log wrlog;
   char str[MAXLINE];
   ErPeak *erPeak;
   Peak *peak;
   int dim,i,j;
   int cnts[CL_EXP][CL_MEAS];
   int totExp[CL_EXP];
   int totMeas[CL_MEAS];
   int nrLocalAss,nrGlobalAss;

   /* reset cnts */
   for(i=0;i<CL_EXP;i++) for(j=0;j<CL_MEAS;j++) cnts[i][j]=0;

   /* count from expected peaks */
   erPeak = (ErPeak *)liErPeak.start();
   while(erPeak != 0) {
      if(erPeak->type == DESTFRAG && erPeak->spec==spec) {
         dim = erPeak->spec->type->dim;
         peak = erPeak->getAss();
         if(peak != 0) {
            nrLocalAss = peak->p->nrIntraDegen(erPeak);
            nrGlobalAss= peak->p->nrAss - nrLocalAss - 1;
            if(nrLocalAss==0 && nrGlobalAss==0) {
               cnts[erPeak->group][0]++;
            } else if(nrLocalAss!=0 && nrGlobalAss==0) {
               cnts[erPeak->group][1]++;
            } else if(nrLocalAss==0 && nrGlobalAss!=0) {
               cnts[erPeak->group][2]++;
            } else {
               cnts[erPeak->group][3]++;
            }
         } else if(erPeak->isPresent()) {
            cnts[erPeak->group][4]++;
         } else {
            cnts[erPeak->group][5]++;
         }
      }
      erPeak = (ErPeak *)liErPeak.next();
   }

   /* count from measured peaks */
   peak = (Peak *)peaks.start();
   while(peak != 0) {
      if(peak->p->spec==spec && peak->p->nrAss==0) cnts[5][0]++;
      peak = (Peak *)peaks.next();
   }

   /* sum up */
   for(i=0;i<CL_EXP;i++) totExp[i]=0;
   for(j=0;j<CL_MEAS;j++) totMeas[j]=0;
   for(i=0;i<CL_EXP;i++) for(j=0;j<CL_MEAS;j++) {
      totExp[i] += cnts[i][j];
      totMeas[j] += cnts[i][j];
   }

   /* report */
   wrlog.w("\n");
   wrlog.w("\n");
   for(i=-1;i<=CL_EXP;i++) { for(j=-1;j<=CL_MEAS;j++) {
      /* headings */
      if(j==-1 && i != -1) {sprintf(str,"%4d   ",i);}
      else if(j!=-1 && i == -1) {sprintf(str,"%4d   ",j); }
      else if(j==-1 && i == -1) {sprintf(str,"       "); }

      /* sums */
      else if(j==CL_MEAS && i != CL_EXP){sprintf(str,"%4d   ",totExp[i]); }
      else if(j!=CL_MEAS && i == CL_EXP){sprintf(str,"%4d   ",totMeas[j]); }
      else if(j==CL_MEAS && i == CL_EXP) {sprintf(str,"       "); }

      /* entries */
      else {sprintf(str,"%4d   ",cnts[i][j]); }

      wrlog.w(str);
      if(j==CL_MEAS) wrlog.w("\n");
   }}
}

void NMRdia::distHamilton()
{
   Log wrlog;
   static int *cnt=0;
   static int *distr=0;
   static int maxCnt = 0;
   int nrPar;
   char str[MAXLINE * 10];
   char str2[MAXLINE];
   int unass,nrpeak,totcnt,mutations,ismut,totPeaksAss;
   ErPeak *erPeak;
   Peak *peak;
   int i;

   if(maxCnt < nrActive) {
      if(cnt) delete cnt;
      cnt=new int [nrActive];
   }
   if(maxCnt < nrActive+1) {
      if(distr) delete distr;
      distr= new int [nrActive+1];
   }

   for(i=0;i<nrActive;i++) { cnt[i]=0; distr[i]=0; }
   distr[nrActive]=0;
   unass=0;
   nrpeak=0;
   totcnt=0;
   mutations=0;
   totPeaksAss=0;

   erPeak = (ErPeak *)liErPeak.start();
   while(erPeak != 0) {		// bonus for assigned peaks
      if(erPeak->type == DESTFRAG) {
         nrPar=0;
         nrpeak++;
         peak=erPeak->getAss();
         if(peak==0) { unass++; }
         else {
            totPeaksAss++;
            ismut=1;
            for(i=0;i<nrActive;i++) {
               if(peak == erPeak->getAss(activePos[i])) {ismut=0; nrPar++;}
               else {
                  cnt[i]++;
               }
            }
            if(ismut) {
               //printf("mutation: "); erPeak->print(0);
               //printf("instead of:");
               //if(erPeak->getAss(activePos[0])==0) printf(" nothing\n");
               //else erPeak->getAss(activePos[0])->print(0);
               mutations++;
            }
            distr[nrPar]++;
         }
      }
      erPeak = (ErPeak *)liErPeak.next();
   }
   for(i=0;i<nrActive;i++) totcnt+=cnt[i];
   ratioMut = (float)mutations/(float)totPeaksAss*100.0;

   sprintf(str,"\n    Percent of unass. peaks .... %8.3f%%",
          (float)unass/(float)totPeaksAss*100.0);
   wrlog.w(str);
   sprintf(str,"\n    Percent of mutations ....... %8.3f%%",
          ratioMut);
   wrlog.w(str);
   if(nrActive>0) {
      sprintf(str,"\n    Average 'Hamming' distance . %8.3f",
         (float)totcnt/(float)totPeaksAss/(float)(nrActive)*1000.0);
      wrlog.w(str);
   }
   sprintf(str,"\n    'Hamming' distances:");
   for(i=0;i<nrActive;i++) {
      sprintf(str2," %4d",xnint((float)cnt[i]/(float)totPeaksAss*1000.0));
      strcat(str,str2);
      if(i%10==9) strcat(str,"\n                       ");
   }
   wrlog.w(str);
   sprintf(str,"\n    Distr. of Mutations:");
   for(i=0;i<nrActive+1;i++) {
      sprintf(str2," %4d",distr[i]);
      strcat(str,str2);
      if(i%10==9) strcat(str,"\n                       ");
   }
   wrlog.w(str);

}

void NMRdia::fixAss(char *t, int fix)
{
   Coherence *co;
   ErPeak *erPeak;
   Log wlog;
   char str[MAXLINE];
   Iterator i;
   int cntCo=0,cntPeaks=0;
   int idim;
   int res;

   calcMeanFreq();
   co = (Coherence *)liCo.startPars(i,t);
   while(co != 0) {
      if(co->firstAt()->frag->type == DESTFRAG && co->meanF.valid()) {
         co->fixed =fix;
         cntCo++;
      }
      co = (Coherence *)liCo.nextPars(i,t);
   }

   erPeak = (ErPeak *)liErPeak.start(i);
   while(erPeak != 0) {
      res=1;
      for(idim=0;idim<erPeak->spec->type->dim;idim++) {
         if(!erPeak->getCo(idim)->fixed) res=0;
      }
      if(res && erPeak->status==P_ASS && erPeak->type == DESTFRAG) {
         erPeak->status=P_FIXED;
         (erPeak->getAss())->p->fixed=1;
         cntPeaks++;
      }
      erPeak = (ErPeak *)liErPeak.next(i);
   }

   sprintf(str,"\n... %d coherences and %d peaks fixed",
	   cntCo,cntPeaks);
   wlog.w(str);
}

void NMRdia::sortFloatAss()
{
   Atom *at1,*at2;
   Coherence *co1,*co2;
   Iterator i1,i2;
   Log wlog;
   char str[MAXLINE];
   int cnt=0;

   /* treats only pairs of covalent equivalent atoms */
   at1=(Atom *)ss.Atoms.start(i1);
   while(at1!=0) {
      if(at1->pseud.getn()==2) {
         at2=(Atom *)at1->pseud.start(i2);
         if(at2 != at1) {
            co1=at1->getCo();
            co2=at2->getCo();
            if(co1!=0 && co2 !=0) cnt+=co1->sortFloatAss(co1,co2);
         }
      }
      at1=(Atom *)ss.Atoms.next(i1);
   }
   sprintf(str,"\n... %d floating assignments exchanged",cnt);
   wlog.w(str);
}

/* mut. inf. are averaged over +/- AVERGING standard deviations */
#define AVERAGING 2.0

void NMRdia::calcScoreTables()
{
   Log wlog;
   char str[MAXLINE];
   char name[MAXLINE];
   int i,j;
   ErPeak *erPeak;
   Coherence *co;
   Peak *peak;

   /* coherences */
   float minFr[NR_NUCLEI], maxFr[NR_NUCLEI];

   wlog.w("\n\n");
   for(i=0;i<NR_NUCLEI;i++) {
      minFr[i]=0.0; maxFr[i]=0.0;
      co = (Coherence *)liCo.start();
      while(co!=0 && co->typeNucleus!=i) co = (Coherence *)liCo.next();
      if(co!=0) {
         minFr[i]=co->expFreq().x;
	 maxFr[i]=co->expFreq().x;
         while(co != 0) {
            if(co->typeNucleus==i && co->expFreq().x!=BAD_PPM && 
               co->expFreq().x < minFr[i])
               minFr[i]=co->expFreq().x;
            if(co->typeNucleus==i && co->expFreq().x!=BAD_PPM &&
               co->expFreq().x > maxFr[i])
               maxFr[i]=co->expFreq().x;
            co = (Coherence *)liCo.next();
         }
      }
      sprintf(str,"Frequency limits for nuclei of type %s: [%5.2f,%5.2f]\n",
         name_nuclei[i],minFr[i],maxFr[i]);
      wlog.w(str);
   }

   co = (Coherence *)liCo.start();
   while(co != 0) {
      if(maxFr[co->typeNucleus]>minFr[co->typeNucleus] && 
	 co->expFreq().v>0.0) {
         co->bonusFreq = log((maxFr[co->typeNucleus]-minFr[co->typeNucleus])/
                             sqrt(6.28) / (co->expFreq().v));
      } else {
         co->bonusFreq = 0.0;
      }
//      co->print(0);
//      sprintf(str,"--> %7.3f %7.3f\n",co->expFreq().v,co->expFreq().x);
//      wlog.w(str);
      co = (Coherence *)liCo.next();
   }

   /* contribution of peaks */
   float cntPeaks[NREXP];
   float totPeaks;

   wlog.w("\n");
   for(i=0;i<usedSpec.used;i++) cntPeaks[i]=0.0;
   totPeaks = 0.0;
   erPeak = (ErPeak *)liErPeak.start();
   while(erPeak != 0) {
      if(erPeak->type == DESTFRAG) {
         cntPeaks[erPeak->spec->nr] += erPeak->pExist; 
         totPeaks += erPeak->pExist;
      }
      erPeak = (ErPeak *)liErPeak.next();
   }

   for(i=0;i<usedSpec.used;i++) {
      if(cntPeaks[i]>0.0 && totPeaks>0.0)
         usedSpec.li[i].bonusPeak = log(totPeaks/cntPeaks[i]);
      else
         usedSpec.li[i].bonusPeak = 0.0;
      usedSpec.li[i].name(name);
      sprintf(str,"Bonus for assigned %s peak: %4.2f\n",
         name,usedSpec.li[i].bonusPeak);
      wlog.w(str);
   }

   /* contribution of assignments */
   float flow,fup;
   float nf;
   float sd;
   int idim, res;
   float penWrong;

   wlog.w("\n");
   for(i=0;i<usedSpec.used;i++) {	/* for each spectrum */
      usedSpec.li[i].name(name);
      for(idim=0;idim<usedSpec.li[i].type->dim;idim++) { /* and each dimension*/
         wlog.w("\n");
         sprintf(str,"w%d of %s:\n",idim+1,name);
         wlog.w(str);
         /* count number of different coherences: nf */
         nf=0.0;
         co = (Coherence *)liCo.start();
         while(co != 0) {
            res=0;
            erPeak = co->startPeak(j);
            while(erPeak != 0 && !res) {
               if(erPeak->spec->nr == usedSpec.li[i].nr &&
                  erPeak->getCo(idim)==co && erPeak->type == DESTFRAG)
		  res=1;
               erPeak = co->nextPeak(j);
            }
            if(res) nf+=1.0;
            co = (Coherence *)liCo.next();
         }
         sprintf(str,"   %d different coherences assigned\n",(int)nf);
         wlog.w(str);
	 /* determine range of frequencies */
         flow=0.0; fup=0.0;
         peak=(Peak *)usedSpec.li[i].liPeakDim[idim].start();
         if(peak != 0) {flow=peak->w[idim]; fup=flow;}
         while(peak != 0) {
            if(peak->w[idim] < flow) flow=peak->w[idim];
            if(peak->w[idim] > fup)  fup=peak->w[idim];
            peak=(Peak *)usedSpec.li[i].liPeakDim[idim].next();
         }
         sprintf(str,"   range of measured peaks: [%5.2f,%5.2f]\n",flow,fup);
         wlog.w(str);
         sd=usedSpec.li[i].type->accPick[idim]/MAX_DEV;
         sprintf(str,"   standard deviation of peak picking: %6.4f\n",sd);
         wlog.w(str);
	 wlog.w("   bonus for assignments:\n      ");
         for(j=0;j<MAX_DEV*RES_SCORE;j++) {
            usedSpec.li[i].bonusAss[idim][j] =
                log(1.0 / (1.0/nf + sqrt(6.28) * sd / (fup-flow) *
                                    (1.0 - 1.0/nf) *
                                    exp(0.5*(float)j/(float)RES_SCORE*
  				            (float)j/(float)RES_SCORE)
                          )
                   );
             sprintf(str,"%6.4f ",usedSpec.li[i].bonusAss[idim][j]);
             wlog.w(str);
         }
         penWrong = log(0.000317/(0.000317/nf + (1-MAX_DEV*2.0*sd/(fup-flow)) *
					        (1.0-1.0/nf)
			         )
                       );
         sprintf(str,"| %6.4f\n",penWrong);
         wlog.w(str);
	 wlog.w("   penalties for degeneracies:\n      ");
         for(j=0;j<MAX_DEV*RES_SCORE;j++) {
            usedSpec.li[i].penDegen[idim][j] =
                log(1.0/( (fup-flow)/sd/nf/sqrt(6.28)*
                          exp(-0.5*(float)j/(float)RES_SCORE*
  				  (float)j/(float)RES_SCORE
                             ) +
                          (1.0-1.0/nf)
                        )
                   );
             sprintf(str,"%6.4f ",usedSpec.li[i].penDegen[idim][j]);
             wlog.w(str);
         }
         usedSpec.li[i].bonusDegen[idim] =
            log(1.0/(0.000317/nf*(fup-flow)/(fup-flow-MAX_DEV*2.0*sd) +
		     (1.0-1.0/nf)
		    )
                );
         sprintf(str,"| %6.4f\n",usedSpec.li[i].bonusDegen[idim]);
         wlog.w(str);
         /* Deviation of shift in spectrum from average value */
         sd=usedSpec.li[i].type->accSpec[idim]/MAX_DEV;
         usedSpec.li[i].bonusDevMean[idim]=
               log((fup-flow) / sqrt(6.28) / sd);
         sprintf(str,"   bonus spectrum frequencies: %5.2f\n",
                    usedSpec.li[i].bonusDevMean[idim]);
         wlog.w(str);
     }
   }

   /* calc. max. possible score for each coherence */
   co = (Coherence *)liCo.start();
   while(co != 0) {
      co->setSumPeakP();
      co = (Coherence *)liCo.next();
   }
   erPeak = (ErPeak *)liErPeak.start();
   while(erPeak != 0) {
      if(erPeak->type == DESTFRAG) {
          erPeak->setSumPeakP();
      }
      erPeak = (ErPeak *)liErPeak.next();
   }

}