bf_dms.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpclinopt.h"
#include "tpcfunc.h"
#include "tpcbfm.h"
/*****************************************************************************/
 
/*****************************************************************************/
int spectralDMSurge(
  const double *x,
  const double *x2,
  const double *y,
  double *w,
  const int sNr,
  const double kMin,
  const double kMax,
  const int fNr,
  const double dtMin,
  const double dtMax,
  const double dtStep,
  double *k,
  double *a,
  double *dtEst,
  double *yfit,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  else if(verbose>1) {
    printf("%s(x, ", __func__);
    if(x2==NULL) printf("null, y, "); else printf("x2, y, ");
    if(w==NULL) printf("null"); else printf("w");
    printf(" , %d, %.1e, %.1e, %d, %g, %g, %g, *k, *a, *dtEst, *yfit\n", 
           sNr, kMin, kMax, fNr, dtMin, dtMax, dtStep);
  }
  if(x==NULL || y==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(TPCERROR_NO_DATA);
  }
  if(sNr<4) {
    if(verbose>1) fprintf(stderr, "invalid number of samples\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    return(TPCERROR_TOO_FEW);
  }
  if(fNr<4) {
    if(verbose>1) fprintf(stderr, "invalid number of functions\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    return(TPCERROR_TOO_FEW);
  }
  if(!(kMin>0.0) || !(kMax>kMin)) {
    if(verbose>1) fprintf(stderr, "invalid k range\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(TPCERROR_INVALID_VALUE);
  }
  double dtRange=dtMax-dtMin;
  if(!(dtRange>=0.0) || (dtRange>0.0 && !(dtStep<0.2*dtRange))) {
    if(verbose>1) fprintf(stderr, "invalid delay time settings\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(TPCERROR_INVALID_VALUE);
  }
 
  /* Allocate memory or set pointers for local a[] and k[] */
  double *lk, *la, *localp;
  localp=(double*)malloc(sizeof(double)*fNr*2);
  if(localp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return(TPCERROR_OUT_OF_MEMORY);
  }
  lk=localp; la=localp+fNr;
 
  if(verbose>2) printf("computing k values\n");
#if(1)
  {
    double r1, r2, s;
    r1=log(kMin); r2=log(kMax); s=(r2-r1)/(double)(fNr-1);
    if(verbose>3) printf(" r1 := %g\n r2 := %g\n s := %g\n", r1, r2, s);
    for(int bi=0; bi<fNr; bi++) lk[bi]=exp((double)bi*s+r1);
  }
#else
  {
    double r1, r2, s;
    r1=log10(kMin); r2=log10(kMax); s=(r2-r1)/(double)(fNr-1);
    if(verbose>3) printf(" r1 := %g\n r2 := %g\n s := %g\n", r1, r2, s);
    for(int bi=0; bi<fNr; bi++) lk[bi]=pow(10.0, (double)bi*s+r1);
  }
#endif
 
  if(verbose>2) printf("computing the range of delay times to test\n");
 
 
  /* Allocate memory required by NNLS */
  if(verbose>2) printf("allocating memory for LLSQ\n");
  NNLSQDATA lsq; 
  nnlsqDataInit(&lsq);
  int ret=0;
  if((ret=nnlsqDataAllocate(&lsq, fNr, sNr))!=TPCERROR_OK) {
    free(localp);
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
 
  /* Set initial delay time to the middle of the given range */
  double initDeltaT=dtMin;
  if(dtMax>dtMin) initDeltaT=0.5*(dtMin+dtMax);
  if(verbose>3) printf("initDeltaT := %g\n", initDeltaT);
 
  /* LLSQ fit with initial delay time */
  if(verbose>2) printf("LLSQ fitting\n");
  double bestDeltaT=initDeltaT;
  double bestR2=nan("");
  double deltaT=initDeltaT;
  if(verbose>2) printf("  deltaT=%g\n", deltaT);
  if(verbose>6) printf("filling data matrix\n");
  /* Fill NNLS B array with measured y values */
  for(int m=0; m<lsq.m; m++) lsq.b[m]=y[m];
  /* Fill NNLS A matrix with basis functions, calculated in place */
  double p[3]={deltaT, 1.0, 0.0};
  for(int n=0, ret=0; n<lsq.n && !ret; n++) {
    p[2]=lk[n];
    if(x2!=NULL) // frame start and end times available
      ret=mfEvalFrameY("dmsurge", 3, p, lsq.m, x, x2, lsq.a[n], 0);
    else
      ret=mfEvalY("dmsurge", 3, p, lsq.m, x, lsq.a[n], 0);
  }
  if(ret) {
    nnlsqDataFree(&lsq); free(localp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(TPCERROR_INVALID_VALUE);
  }
  /* Apply weights, if given */
  if(w!=NULL) nnlsqDataWght(&lsq, w);
  /* NNLS */
  if(verbose>6) printf("applying NNLS\n");
  lsq.iternr=3*lsq.n;
  ret=nnlsq(&lsq, verbose-1);
  if(ret>1) {
    nnlsqDataFree(&lsq); free(localp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_SOLUTION);
    return(TPCERROR_NO_SOLUTION);
  }
  if(verbose>2) printf("  -> r2=%g  iterNr=%d\n", lsq.rnorm, lsq.iternr);
  bestR2=lsq.rnorm;
  /* Copy a[] values */
  for(int n=0; n<lsq.n; n++) la[n]=lsq.x[n];
 
  /* Try moving function to the left with delay steps */
  deltaT=initDeltaT-dtStep;
  while(dtStep>0.0 && deltaT>=dtMin) {
    if(verbose>2) printf("  deltaT=%g\n", deltaT);
    for(int m=0; m<lsq.m; m++) lsq.b[m]=y[m];
    p[0]=deltaT;
    for(int n=0, ret=0; n<lsq.n && !ret; n++) {
      p[2]=lk[n];
      if(x2!=NULL) // frame start and end times available
        ret=mfEvalFrameY("dmsurge", 3, p, lsq.m, x, x2, lsq.a[n], 0);
      else
        ret=mfEvalY("dmsurge", 3, p, lsq.m, x, lsq.a[n], 0);
    }
    if(ret) {
      nnlsqDataFree(&lsq); free(localp);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return(TPCERROR_INVALID_VALUE);
    }
    if(w!=NULL) nnlsqDataWght(&lsq, w);
    lsq.iternr=3*lsq.n;
    int ret=nnlsq(&lsq, verbose-1);
    if(ret>1) break;
    if(verbose>2) printf("  -> r2=%g  iterNr=%d\n", lsq.rnorm, lsq.iternr);
    if(lsq.rnorm<bestR2) {
      bestR2=lsq.rnorm;
      bestDeltaT=deltaT;
      for(int n=0; n<lsq.n; n++) la[n]=lsq.x[n];
    } //else if(nnls_r2>bestR2) break;
    deltaT-=dtStep;
  }
  /* Try moving function to the right with delay steps */
  deltaT=initDeltaT+dtStep;
  while(dtStep>0.0 && deltaT<=dtMax) {
    if(verbose>2) printf("  deltaT=%g\n", deltaT);
    for(int m=0; m<lsq.m; m++) lsq.b[m]=y[m];
    p[0]=deltaT;
    //struct timespec ts[3];
    //timespec_get(&ts[0], TIME_UTC);
    for(int n=0, ret=0; n<lsq.n && !ret; n++) {
      p[2]=lk[n];
      if(x2!=NULL) // frame start and end times available
        ret=mfEvalFrameY("dmsurge", 3, p, lsq.m, x, x2, lsq.a[n], 0);
      else
        ret=mfEvalY("dmsurge", 3, p, lsq.m, x, lsq.a[n], 0);
    }
    if(ret) {
      nnlsqDataFree(&lsq); free(localp);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return(TPCERROR_INVALID_VALUE);
    }
    //timespec_get(&ts[1], TIME_UTC);
    if(w!=NULL) nnlsqDataWght(&lsq, w);
    lsq.iternr=3*lsq.n;
    int ret=nnlsq(&lsq, verbose-1);
    //timespec_get(&ts[2], TIME_UTC);
    if(ret>1) break;
    if(verbose>2) printf("  -> r2=%g  iterNr=%d\n", lsq.rnorm, lsq.iternr);
    if(lsq.rnorm<bestR2) {
      bestR2=lsq.rnorm;
      bestDeltaT=deltaT;
      for(int n=0; n<lsq.n; n++) la[n]=lsq.x[n];
    } //else if(nnls_r2>bestR2) break;
    //printf("  time_to-fill=%g    time_to_nnls=%g\n", timespecDifference(&ts[1], &ts[0]), timespecDifference(&ts[2], &ts[1]));
    deltaT+=dtStep;
  }
 
  /* Compute fitted y[] */
  if(yfit!=NULL) {
    if(verbose>1) printf("computing yfit[]\n");
    for(int m=0; m<lsq.m; m++) {
      yfit[m]=0.0;
      for(int n=0; n<lsq.n; n++)
        if(la[n]>0.0)
          yfit[m]+=la[n]*(x[m]+bestDeltaT)*exp(-lk[n]*(x[m]+bestDeltaT));
    }
  }
 
  nnlsqDataFree(&lsq);
 
  /* Copy a[] and b[] if requested */
  if(a!=NULL) for(int bi=0; bi<fNr; bi++) a[bi]=la[bi];
  if(k!=NULL) for(int bi=0; bi<fNr; bi++) k[bi]=lk[bi];
  free(localp);
 
  if(dtEst!=NULL) *dtEst=bestDeltaT;
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/