/*---------------------------------------------------------------------------
                          P E A K I N T E G R A T E
                          
   Scan Peak integration along the scheme described by Grant & Gabe in
   J. Appl. Cryst (1978), 11, 114-120
 
 
   Mark Koennecke, March 1998
   
--------------------------------------------------------------------------*/
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "fortify.h"
#include "integrate.h"

/* define PRINT if you want more info about the peak limits and what is
in between

#define PRINT 1
*/

#define PROBABILITY 0.01
#define LEFT -1
#define RIGHT 1

#define ABS(x) (x < 0 ? -(x) : (x))

/*----------------- a private data  structure to keep track of things ----*/

typedef struct {
  long *lCounts;
  int iCounts;
  int m;
  int iLeft;
  int iMaxCount;
  int iRight;
} IntegData, *pIntegData;

/*-------------------------------------------------------------------------*/
static int FindMaxCount(long lCounts[], int iCounts)
{
  int i, iPtr;
  long l = -100000;

  for (i = 0, iPtr = 0; i < iCounts; i++) {
    if (lCounts[i] > l) {
      l = lCounts[i];
      iPtr = i;
    }
  }
  return iPtr;
}

/*------------------------------------------------------------------------
  Finds the index where the counts are less then half peak height. Assumes
  that iMaxCount has been set to a proper value beforehand. This helps to 
  solve a problem: Given very broad peaks, the window condition would be
  fulfilled when starting for background counts from the top. The solution
  is to search for the half width point and start searchnig for background 
  from there. Saves calculation time as well. FindHalf returns an index on
  success, else - 1 if it goes beyond limits!
  -----------------------------------------------------------------------*/
static int FindHalf(pIntegData self, int iSign)
{
  int i;
  long lHalf;

  assert(self);

  lHalf = self->lCounts[self->iMaxCount] / 2;
  i = self->iMaxCount;
  while (self->lCounts[i] > lHalf) {
    i += iSign;
    if ((i < 0) || (i > self->iCounts)) {
      return -1;
    }
  }
  return i;
}

/*--------------------------------------------------------------------------
  FindWindowHits find the number of hits in the window defined by
  0.67 sigma of the value at iPos and m. Returns -1 if insufficient data
  is available to do the trick. iSign can  be -1 for  search to left or
  1 for search to right. 
---------------------------------------------------------------------------*/
static int FindWindowHits(pIntegData self, int iPos, int iSign)
{
  float fWinMin, fWinMax, fVar;
  int i, iTest, iCount;

  assert(self);
  assert((iSign == -1) || (iSign == 1));

  /* check limits of the show */
  iTest = iPos + iSign * self->m;
  if ((iTest < 0) || (iTest >= self->iCounts)) {
    return -1;
  }

  /* calculate limits */
  if (self->lCounts[iPos] > 0) {
    fVar = sqrt((double) self->lCounts[iPos]);
    fWinMin = (float) self->lCounts[iPos] - 0.67 * fVar;
    fWinMax = (float) self->lCounts[iPos] + 0.67 * fVar;
  } else {
    fWinMin = 0.;
    fWinMax = 1.;
  }

  /* count the window hits */
  for (i = 0, iCount = 0; i < self->m; i++) {
    iTest = iPos + iSign * i;
    if ((self->lCounts[iTest] >= fWinMin) &&
        (self->lCounts[iTest] <= fWinMax)) {
      iCount++;
    }
  }
  return iCount;
}

/*----------------------------------------------------------------------------
  AdvanceBackground advances in the background until the predefined
  Probability has been reached. The probability is defined at the top
  of the file. iSign and errors are the same as above.
----------------------------------------------------------------------------*/
static int AdvanceBackground(pIntegData self, int iPos, int iSign)
{
  float fProb;
  int iRet, iTest;

  assert(self);
  assert((iSign == 1) || (iSign == -1));

  /* initialize, calculate initial Probability */
  iRet = FindWindowHits(self, iPos, iSign);
  if (iRet < 0) {
    return iRet;
  }
  fProb = (self->m - iRet) / self->m;
  iTest = iPos;

  /* loop further on */
  while (fProb > PROBABILITY) {
    iTest += iSign;
    iRet = FindWindowHits(self, iTest, iSign);
    if (iRet < 0) {
      return iRet;
    }
    fProb *= (self->m - iRet) / self->m;
  }

  if ((iTest < 0) || (iTest >= self->iCounts)) {
    return -1;
  } else {
    return iTest;
  }
}

/*--------------------------------------------------------------------------
  FindLimit assumes, that iMaxCount has been set to something sensible
  in advance!
--------------------------------------------------------------------------*/
static int FindLimit(pIntegData self, int iStart, int iSign)
{
  int i, iHalf, iHit, iLeft, iTest;

  assert(self);
  assert((iSign == 1) || (iSign == -1));

  iHit = 0;
  iHalf = self->m / 2;
  i = iStart;

  /* loop until iHalf is reached */
  while ((iHit = FindWindowHits(self, i, iSign)) < iHalf) {
    i += iSign;
    if (iHit < 0) {
      if (iSign == LEFT) {
        return INTEGLEFT;
      } else {
        return INTEGRIGHT;
      }
    }
  }

  /* step further into the backgound */
  iTest = AdvanceBackground(self, i, iSign);
  if (iTest < 0) {
    if (iSign == LEFT) {
      return INTEGLEFT;
    } else {
      return INTEGRIGHT;
    }
  }

  /* store it away and return success */
  if (iSign == LEFT) {
    self->iLeft = iTest;
  } else {
    self->iRight = iTest;
  }
  return 1;
}

/*--------------------------------------------------------------------------
  This routine does the actual integration once the limits have been 
  determined.
 --------------------------------------------------------------------------*/
static int DoIntegrate(pIntegData self, float *fSum, float *fVariance)
{
  long lLeftBack = 0, lRightBack, lBackMean;
  int i, iRes = 1;
  float fScan, fBackVar, fLeft, fRight, fN, fNP, fNB, fSumm;
  float fMLeft, fMRight;

  assert(self);
  assert((self->iLeft != 0) && (self->iRight != 0));


  /* sum the left background */
  for (i = 0, lLeftBack = 0; i < self->iLeft; i++) {
    lLeftBack += self->lCounts[i];
  }

  /* sum the right background */
  for (i = self->iRight, lRightBack = 0; i < self->iCounts; i++) {
    lRightBack += self->lCounts[i];
  }


  /* sum all of it  */
  fScan = 0.;
  for (i = 0; i < self->iCounts; i++) {
    fScan += self->lCounts[i];
  }

  /* sum the peak */
  fSumm = 0.;
  for (i = self->iLeft; i < self->iRight + 1; i++) {
    fSumm += (float) self->lCounts[i];
  }

  /* calculate the values */
  fLeft = (float) self->iLeft;
  fRight = (float) self->iRight;
  fN = (float) self->iCounts;

  fNP = fRight - fLeft + 1.;
  fNB = fLeft + (fN - fRight);

  *fSum = fSumm - (fNP / fNB) * ((float) lLeftBack + (float) lRightBack);
  *fVariance = sqrt(fSumm + (fNP / fNB) * (fNP / fNB) * ((float) lLeftBack
                                                         +
                                                         (float)
                                                         lRightBack));

  /* check the background */
  fMLeft = lLeftBack / fLeft;
  fMRight = lRightBack / fRight;
  if (ABS(fMLeft - fMRight) > ((fMLeft + fMRight) / 2)) {
    iRes = INTEGFUNNYBACK;
  }

  return iRes;
}

/*---------------------------------------------------------------------------*/
int GabePeakIntegrate(int m, int iCounts, long lCounts[],
                      float *fIntens, float *fVariance)
{
  IntegData self;
  int iTest, i, iStart;

  /* check input */
  assert(iCounts > 0);
  assert(lCounts);


  /* initialize Data Structure */
  self.iCounts = iCounts;
  self.lCounts = lCounts;       /* no need to duplicate data */
  self.iLeft = 0;
  self.iRight = 0;
  self.iMaxCount = 0;
  self.m = m;

  /* find maximum */
  self.iMaxCount = FindMaxCount(self.lCounts, self.iCounts);

  /* find Left Limit */
  iStart = FindHalf(&self, LEFT);
  if (iStart < 0) {
    return INTEGLEFT;
  }
  iTest = FindLimit(&self, iStart, LEFT);
  if (iTest < 0) {
    return iTest;
  }

  /* find right limit */
  iStart = FindHalf(&self, RIGHT);
  if (iStart < 0) {
    return INTEGRIGHT;
  }
  iTest = FindLimit(&self, iStart, RIGHT);
  if (iTest < 0) {
    return iTest;
  }
#ifdef PRINT
  printf("Limits: %d %d\n", self.iLeft, self.iRight);
  for (i = self.iLeft; i < self.iRight; i++) {
    printf(" %d\n", self.lCounts[i]);
  }
#endif
  /* integrate */
  return DoIntegrate(&self, fIntens, fVariance);

}