sics/hmdata.c

/*-----------------------------------------------------------------------
  This is a data handling class for histogram memory data.
  For more information see hmdata.tex.

  copyright: see file COPYRIGHT

  Mark Koennecke, January 2003
  
 Added loading HM data from file, Mark Koennecke, November 2006
-------------------------------------------------------------------------*/
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <time.h>
#include "splitter.h"
#include "fortify.h"
#include "hmdata.h"
#include <nxdataset.h>
#include "HistMem.h"
#include "HistMem.i"
#include "HistDriv.i"
#include "countdriv.h"
#include "stptok.h"
/*----------------------------------------------------------------------*/
pHMdata makeHMData(void)
{
  pHMdata self = NULL;

  self = (pHMdata) malloc(sizeof(HMdata));
  if (self == NULL) {
    return NULL;
  }
  memset(self, 0, sizeof(HMdata));
  self->nTimeChan = 1;

  return self;
}

/*---------------------------------------------------------------------*/
void killHMData(pHMdata self)
{
  if (self->localBuffer != NULL) {
    free(self->localBuffer);
  }
  free(self);
}

/*---------------------------------------------------------------------*/
void clearHMData(pHMdata self)
{
  long size;
  int i;
  size = 1;
  for (i = 0; i < self->rank; i++) {
    size *= self->iDim[i];
  }
  if (self->tofMode) {
    size *= getNoOfTimebins(self);
  }
  if(self->localBuffer != NULL){
    memset(self->localBuffer, 0, size * sizeof(HistInt));
  }
  self->nextUpdate = time(NULL) + self->updateIntervall;
}

/*----------------------------------------------------------------------*/
int resizeBuffer(pHMdata self)
{
  long size;
  int i;

  size = 1;
  for (i = 0; i < self->rank; i++) {
    size *= self->iDim[i];
  }
  if (self->tofMode) {
    size *= getNoOfTimebins(self);
  }
  if (self->localBuffer != NULL) {
    free(self->localBuffer);
    self->localBuffer = NULL;
  }
  self->localBuffer = (HistInt *) malloc(size * sizeof(HistInt));
  if (!self->localBuffer) {
    return 0;
  }
  memset(self->localBuffer, 0, size * sizeof(HistInt));
  self->nextUpdate = 0;
  return 1;
}

/*-----------------------------------------------------------------------*/
int configureHMdata(pHMdata self, pStringDict pOpt, SConnection * pCon)
{
  int status, i;
  float fVal;
  char pValue[80];
  pHistMem master = NULL;

  if (self->nTimeChan > 2) {
    self->tofMode = 1;
  } else {
    self->tofMode = 0;
  }
  status = StringDictGetAsNumber(pOpt, "rank", &fVal);
  if (!status) {
    SCWrite(pCon, "ERROR: critical configuration problem: no rank found",
            eError);
    return 0;
  }
  self->rank = (int) rint(fVal);

  for (i = 0; i < self->rank; i++) {
    snprintf(pValue,sizeof(pValue)-1, "dim%1.1d", i);
    status = StringDictGetAsNumber(pOpt, pValue, &fVal);
    if (!status) {
      snprintf(pValue,sizeof(pValue)-1, "ERROR dimension %d not found!!", i);
      return 0;
    }
    self->iDim[i] = (int) rint(fVal);
  }

  status = StringDictGetAsNumber(pOpt, "update", &fVal);
  if (!status) {
    self->updateIntervall = 0;  /* no buffering */
  } else {
    self->updateIntervall = (int) rint(fVal);
  }

  status = StringDictGet(pOpt, "timeslave", pValue, 79);
  if (status == 1) {
    master = (pHistMem) FindCommandData(pServ->pSics, pValue, "HistMem");
    if (master == NULL) {
      SCWrite(pCon, "ERROR: timeslave requested, but master HM not found",
              eError);
    } else {
      self->timeslave = master->pDriv->data;
      self->tofMode = 1;
    }
  }

  /*
     invalidate buffer
   */
  if (self->localBuffer != NULL) {
    free(self->localBuffer);
    self->localBuffer = NULL;
  }

  status = resizeBuffer(self);
  if (!status) {
    SCWrite(pCon, "ERROR: failed to resize buffer", eError);
    return 0;
  }
  return 1;
}

/*----------------------------------------------------------------------*/
int genTimeBinning(pHMdata self, float start, float step, int noSteps)
{
  int i;

  if (noSteps >= MAXCHAN || self->timeslave != NULL) {
    return 0;
  }
  for (i = 0; i < noSteps; i++) {
    self->timeBinning[i] = start + i * step;
  }
  self->tofMode = 1;
  self->nTimeChan = noSteps;
  return resizeBuffer(self);
}

/*----------------------------------------------------------------------*/
int setTimeBin(pHMdata self, int index, float value)
{
  if (self->timeslave != NULL) {
    return 0;
  }

  if (index >= 0 && index < MAXCHAN) {
    self->timeBinning[index] = value;
  } else {
    return 0;
  }
  self->tofMode = 1;
  if (index >= self->nTimeChan) {
    self->nTimeChan = index + 1;
    return resizeBuffer(self);
  }
  return 1;
}

/*-------------------------------------------------------------------*/
int isInTOFMode(pHMdata self)
{
  return self->tofMode;
}

/*---------------------------------------------------------------------*/
int getNoOfTimebins(pHMdata self)
{
  if (self->timeslave != NULL) {
    return getNoOfTimebins(self->timeslave);
  } else {
    return self->nTimeChan;
  }
}

/*---------------------------------------------------------------------*/
float *getTimeBinning(pHMdata self)
{
  if (self->timeslave != NULL) {
    return getTimeBinning(self->timeslave);
  } else {
    return self->timeBinning;
  }
}

/*-------------------------------------------------------------------*/
void clearTimeBinning(pHMdata self)
{
  if (self->timeslave == NULL) {
    self->nTimeChan = 1;
    self->tofMode = 0;
    resizeBuffer(self);
  }
}

/*--------------------------------------------------------------------*/
void getHMDataDim(pHMdata self, int iDim[MAXDIM], int *rank)
{
  memcpy(iDim, self->iDim, self->rank * sizeof(int));
  *rank = self->rank;
}

/*---------------------------------------------------------------------*/
long getHMDataLength(pHMdata self)
{
  long length = 1;
  int i;
  for (i = 0; i < self->rank; i++) {
    length *= self->iDim[i];
  }
  if (self->tofMode) {
    length *= getNoOfTimebins(self);
  }
  return length;
}

/*---------------------------------------------------------------------*/
void updateHMData(pHMdata self)
{
  self->nextUpdate = 0;
}

/*--------------------------------------------------------------------
  The idea here is that upper level code sets the updateFlag through 
  updateHMData (above) whenever the HM changes (counts). If this flag is set
  the next call to get getHMDataHistogram will read a new copy from the HM.
  After reading nextUpdate is set to time +  updateIntervall.  In order to
  prevent clients hammering the HM nextUpdate is checked as well. 
  updateIntervall can be set to a reasonable time intervall between updates       in seconds.   If updateIntervall is 0, a direct read is always perfomed.
  If this system needs to be bypassed altogether (because there is no memory
  to buffer all HM) use GetHistogramDirect (histogram.c) instead which acts 
  on the driver level.
  --------------------------------------------------------------------------*/
static int mustUpdate(pHMdata self)
{
  if (time(NULL) >= self->nextUpdate) {
    return 1;
  } else {
    return 0;
  }
}

/*---------------------------------------------------------------------*/
static int updateHMbuffer(pHistMem hist, int bank, SConnection * pCon)
{
  int status, iErr, i;
  char pError[80], pBueffel[256];
  pHMdata self = hist->pDriv->data;

  assert(self);

  if (self->timeslave != NULL) {
    resizeBuffer(self);
  }

  for (i = 0; i < 3; i++) {
    status = hist->pDriv->GetHistogram(hist->pDriv, pCon,
                                       bank, 0, getHMDataLength(self),
                                       self->localBuffer);
    if (status == OKOK) {
      self->nextUpdate = time(NULL) + self->updateIntervall;
      tracePar(hist->pDes->name,"!!datachange!!");
      break;
    } else {
      status = hist->pDriv->GetError(hist->pDriv, &iErr, pError, 79);
      snprintf(pBueffel,sizeof(pBueffel)-1, "ERROR: %s ", pError);
      SCWrite(pCon, pBueffel, eError);
      status = hist->pDriv->TryAndFixIt(hist->pDriv, iErr);
      if (status == COTERM) {
        return 0;
      }
    }
  }
  if (status == OKOK) {
    return 1;
  } else { 
    return HWFault;
  }
}

/*----------------------------------------------------------------------*/
int getHMDataHistogram(pHistMem hist, SConnection * pCon,
                       int bank, int start, int length, HistInt * lData)
{
  int status;
  pHMdata self = hist->pDriv->data;
  HistInt *lStart;

  assert(self);

  if (self->localBuffer == NULL) {
    resizeBuffer(self);
  }

  /*
     update buffer if necessary
   */
  if (mustUpdate(self)) {
    status = updateHMbuffer(hist, bank, pCon);
    if (status != OKOK) {
      return status;
    }
  }

  /*
     copy buffered data to lData
   */
  lStart = self->localBuffer + start;
  if (start + length > getHMDataLength(self)) {
    length = getHMDataLength(self) - start - 1;
  }
  memcpy(lData, lStart, length * sizeof(HistInt));
  return 1;
}

/*-----------------------------------------------------------------------*/
HistInt *getHMDataBufferPointer(pHistMem hist, SConnection * pCon)
{
  int status;
  pHMdata self = hist->pDriv->data;

  assert(self);

  if (self->localBuffer == NULL || self->timeslave != NULL) {
    resizeBuffer(self);
  }
  /*
     update buffer if necessary
   */
  if (mustUpdate(self)) {
    status = updateHMbuffer(hist, 0, pCon);
    if (status != OKOK) {
      return NULL;
    }
  }
  return self->localBuffer;
}

/*------------------------------------------------------------------------*/
static long SumRow(HistInt * iData, int iDataLength, int iStart, int iEnd)
{
  int i;
  long lSum;

  if (iEnd > iDataLength) {
    return -1;
  }

  lSum = 0;
  for (i = iStart; i < iEnd; i++) {
    lSum += iData[i];
  }
  return lSum;
}

/*--------------------------------------------------------------------------*/
long sumHMDataRectangle(pHistMem hist, SConnection * pCon,
                        int iStart[MAXDIM], int iEnd[MAXDIM])
{
  HistInt *iData;
  pHMdata self = hist->pDriv->data;
  int i, iHistLength, status, iIndex, myrank;
  char pBueffel[256];
  unsigned long lSum, lRowSum;

  assert(self);

  /*
     error checking
   */
  myrank = self->rank;
  if (isInTOFMode(self)) {
    self->iDim[self->rank] = getNoOfTimebins(self);
    myrank++;
  }
  for (i = 0; i < myrank; i++) {
    if ((iStart[i] < 0) || (iStart[i] > self->iDim[i])) {
      snprintf(pBueffel,sizeof(pBueffel)-1, "ERROR: %d is out of data dimension range",
              iStart[i]);
      SCWrite(pCon, pBueffel, eError);
      return -1;
    }
    if ((iEnd[i] < 0) || (iEnd[i] > self->iDim[i])) {
      snprintf(pBueffel,sizeof(pBueffel)-1, "ERROR: %d is out of data dimension range",
              iEnd[i]);
      SCWrite(pCon, pBueffel, eError);
      return -1;
    }
  }

  if (self->localBuffer == NULL) {
    resizeBuffer(self);
  }

  /*
     get an update of the HM if necessary
   */
  if (mustUpdate(self)) {
    status = updateHMbuffer(hist, 0, pCon);
    if (status != OKOK) {
      return -1;
    }
  }

  iHistLength = getHMDataLength(self);
  /* actually sum */
  switch (myrank) {
  case 1:
    lSum = SumRow(self->localBuffer, iHistLength, iStart[0], iEnd[0]);
    break;
  case 2:
    if (isInTOFMode(self)) {
      lSum = 0;
      for (i = iStart[0]; i < iEnd[0]; i++) {
        iIndex = i * self->iDim[1];
        lRowSum = SumRow(self->localBuffer, iHistLength,
                         iIndex + iStart[1], iIndex + iEnd[1]);
        lSum += lRowSum;
      }
    } else {
      lSum = 0;
      for (i = iStart[1]; i < iEnd[1]; i++) {
        iIndex = i * self->iDim[0];
        lRowSum = SumRow(self->localBuffer, iHistLength,
                         iIndex + iStart[0], iIndex + iEnd[0]);
        lSum += lRowSum;
      }
    }
    break;
  default:
    snprintf(pBueffel,sizeof(pBueffel)-1,
            "ERROR: summing in %d dimensions not yet implemented", myrank);
    SCWrite(pCon, pBueffel, eError);
    return -1;
    break;
  }
  return lSum;
}

/*--------------------------------------------------------------------------*/
int loadHMData(pHMdata self, SConnection * pCon, char *filename)
{
  FILE *fd = NULL;
  char buffer[1024], pNumber[80], *pPtr;
  long i = 0, length;
  HistInt *data = NULL;

  fd = fopen(filename, "r");
  if (fd == NULL) {
    snprintf(buffer, 1023, "ERROR: failed to open file %s", filename);
    SCWrite(pCon, buffer, eError);
    return 0;
  }
  length = getHMDataLength(self);
  if (self->localBuffer == NULL || self->timeslave != NULL) {
    resizeBuffer(self);
  }
  data = self->localBuffer;
  if (data == NULL) {
    SCWrite(pCon, "ERROR: failed to allocate HM", eError);
    fclose(fd);
    return 0;
  }
  while (i < length && fgets(buffer, 1024, fd) != NULL) {
    pPtr = buffer;
    while (pPtr != NULL) {
      pPtr = sicsNextNumber(pPtr, pNumber);
      if (pPtr != NULL) {
        data[i] = atoi(pNumber);
        i++;
      }
    }
  }
  if (i < length - 1) {
    SCWrite(pCon, "WARNING: not enough data in file to fill HM", eWarning);
  }
  fclose(fd);
  return 1;
}

/*==========================================================================
 * subsampling was stolen from the SinqHTTP histogram memory code and 
 * thus contains some additional indirections.
 * =========================================================================*/
static pNXDS hmDataToNXDataset(pHMdata self)
{
  pNXDS result = NULL;
  int i;

  result = malloc(sizeof(NXDS));
  if (result == NULL) {
    return NULL;
  }
  memset(result, 0, sizeof(NXDS));
  result->magic = MAGIC;
  result->type = NX_INT32;
  result->rank = self->rank;
  if (isInTOFMode(self)) {
    result->rank++;
  }
  result->dim = malloc(self->rank * sizeof(int64_t));
  if (result->dim == NULL) {
    free(result);
    return NULL;
  }
  for (i = 0; i < self->rank; i++) {
    result->dim[i] = self->iDim[i];
  }
  if (isInTOFMode(self)) {
    result->dim[result->rank - 1] = getNoOfTimebins(self);
  }
  if (self->localBuffer == NULL) {
    resizeBuffer(self);
  }
  result->u.iPtr = self->localBuffer;
  return result;
}

/*---------------------------------------------------------------------------*/
static pNXDS subSampleCommand(pNXDS source, char *command,
                              char *error, int errLen)
{
  int64_t startDim[NX_MAXRANK], endDim[NX_MAXRANK];
  int dim = 0, start = 0, i;
  char *pPtr = NULL, token[80];


  pPtr = stptok(command, token, 79, ":\0");
  while ((pPtr = stptok(pPtr, token, 79, ":\0")) != NULL) {
    if (start == 0) {
      startDim[dim] = atoi(token);
      start = 1;
    } else {
      endDim[dim] = atoi(token);
      start = 0;
      dim++;
    }
  }

  if (dim < source->rank - 1) {
    strlcpy(error,
            "ERROR: Not enough border values specified for subsampling",
            errLen);
    return NULL;
  }
  for (i = 0; i < source->rank; i++) {
    if (startDim[i] < 0 || startDim[i] >= source->dim[i]) {
      snprintf(error, errLen,
               "ERROR: invalid start value %d for dimension %d",
               (int)startDim[1], i);
      return NULL;
    }
    if (endDim[i] < startDim[i] || endDim[i] >= source->dim[i]) {
      snprintf(error, errLen,
               "ERROR: invalid end value %d for dimension %d", (int)endDim[1],
               i);
      return NULL;
    }
  }

  return cutNXDataset(source, startDim, endDim);
}

/*-----------------------------------------------------------------------------*/
static pNXDS sumCommand(pNXDS source, char *command, char *error,
                        int errlen)
{
  int dimNo = -1, start = -1, end = -1;
  char *pPtr = NULL;
  char token[80];

  pPtr = stptok(command, token, 79, ":\0");
  pPtr = stptok(pPtr, token, 79, ":\0");
  if (pPtr != NULL) {
    dimNo = atoi(token);
  }
  pPtr = stptok(pPtr, token, 79, ":\0");
  if (pPtr != NULL) {
    start = atoi(token);
  }
  pPtr = stptok(pPtr, token, 79, ":\0");
  if (pPtr != NULL) {
    end = atoi(token);
  }
  if (dimNo < 0 || dimNo > source->rank - 1) {
    snprintf(error, errlen, "ERROR: invalid dimension %d requestd to sum",
             dimNo);
    return NULL;
  }
  if (end < 0 || end > source->dim[dimNo] || start < 0 || start > end) {
    snprintf(error, errlen,
             "ERROR: invalid summing limits %d to %d requested", start,
             end);
    return NULL;
  }

  return sumNXDataset(source, dimNo, start, end);
}

/*--------------------------------------------------------------------------*/
HistInt *subSample(pHMdata self, char *command, char *error, int errLen)
{
  pNXDS source = NULL, start = NULL;
  pNXDS result = NULL;
  char *pPtr = NULL;
  char subCommand[132];
  HistInt *data = NULL;
  int length;


  start = hmDataToNXDataset(self);
  if (start == NULL) {
    strlcpy(error, "Out-Of-Memory or no data while subsampling ", errLen);
    return NULL;
  }

  source = start;
  pPtr = command;
  while ((pPtr = stptok(pPtr, subCommand, 131, ";\0\r\n")) != NULL) {
    if (strstr(subCommand, "sample") != NULL) {
      result = subSampleCommand(source, subCommand, error, errLen);
    } else if (strstr(subCommand, "sum") != NULL) {
      result = sumCommand(source, subCommand, error, errLen);
    } else {
      strlcpy(error, "ERROR: invalid subcommand to process requested",
              errLen);
      return NULL;
    }
    if (result == NULL) {
      return NULL;
    }
    if (source != start) {
      dropNXDataset(source);
    }
    source = result;
  }
  length = getNXDatasetLength(result);
  data = malloc((length + 1) * sizeof(int));
  if (data == NULL) {
    strlcpy(error, "Out-Of-Mmeory in Subsample", errLen);
    dropNXDataset(result);
    return NULL;
  }
  data[0] = length;
  memcpy(data + 1, result->u.iPtr, length * sizeof(int));
  dropNXDataset(result);
  return data;
}