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");
  {
    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);
  }
 
 
  /* Allocate memory required by NNLS */
  if(verbose>2) printf("allocating memory for LLSQ\n");
  int     nnls_n, nnls_m, nnls_index[fNr];
  double *nnls_a[fNr], nnls_x[fNr], nnls_wp[fNr], *nnls_b, *nnls_zz, *nnls_mat, *dptr, nnls_r2;
  nnls_n=fNr; nnls_m=sNr;
  nnls_mat=(double*)malloc(((nnls_n+2)*nnls_m)*sizeof(double));
  if(nnls_mat==NULL) {
    free(localp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return(TPCERROR_OUT_OF_MEMORY);
  }
  dptr=nnls_mat;
  for(int n=0; n<nnls_n; n++) {nnls_a[n]=dptr; dptr+=nnls_m;}
  nnls_b=dptr; dptr+=nnls_m; nnls_zz=dptr;
 
  /* 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<nnls_m; m++) nnls_b[m]=y[m];
  /* Fill NNLS A matrix with basis functions, calculated in place */
  int ret=0;
  double p[3]={deltaT, 1.0, 0.0};
  for(int n=0; n<nnls_n && !ret; n++) {
    p[2]=lk[n];
    if(x2!=NULL) // frame start and end times available
      ret=mfEvalFrameY("dmsurge", 3, p, nnls_m, x, x2, nnls_a[n], 0);
    else
      ret=mfEvalY("dmsurge", 3, p, nnls_m, x, nnls_a[n], 0);
  }
  if(ret) {
    free(nnls_mat); free(localp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(TPCERROR_INVALID_VALUE);
  }
  /* Apply weights, if given */
  if(w!=NULL) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, w);
  /* NNLS */
  if(verbose>6) printf("applying NNLS\n");
  ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_r2, nnls_wp, nnls_zz, nnls_index);
  if(ret>1) {
    free(nnls_mat); free(localp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_SOLUTION);
    return(TPCERROR_NO_SOLUTION);
  }
  if(verbose>2) printf("  -> r2=%g\n", nnls_r2);
  bestR2=nnls_r2;
  /* Copy a[] values */
  for(int n=0; n<nnls_n; n++) la[n]=nnls_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<nnls_m; m++) nnls_b[m]=y[m];
    p[0]=deltaT;
    for(int n=0; n<nnls_n && !ret; n++) {
      p[2]=lk[n];
      if(x2!=NULL) // frame start and end times available
        ret=mfEvalFrameY("dmsurge", 3, p, nnls_m, x, x2, nnls_a[n], 0);
      else
        ret=mfEvalY("dmsurge", 3, p, nnls_m, x, nnls_a[n], 0);
    }
    if(ret) {
      free(nnls_mat); free(localp);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return(TPCERROR_INVALID_VALUE);
    }
    if(w!=NULL) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, w);
    int ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_r2, nnls_wp, nnls_zz, nnls_index);
    if(ret>1) break;
    if(verbose>2) printf("  -> r2=%g\n", nnls_r2);
    if(nnls_r2<bestR2) {
      bestR2=nnls_r2;
      bestDeltaT=deltaT;
      for(int n=0; n<nnls_n; n++) la[n]=nnls_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<nnls_m; m++) nnls_b[m]=y[m];
    p[0]=deltaT;
    //struct timespec ts[3];
    //timespec_get(&ts[0], TIME_UTC);
    for(int n=0; n<nnls_n && !ret; n++) {
      p[2]=lk[n];
      if(x2!=NULL) // frame start and end times available
        ret=mfEvalFrameY("dmsurge", 3, p, nnls_m, x, x2, nnls_a[n], 0);
      else
        ret=mfEvalY("dmsurge", 3, p, nnls_m, x, nnls_a[n], 0);
    }
    if(ret) {
      free(nnls_mat); free(localp);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return(TPCERROR_INVALID_VALUE);
    }
    //timespec_get(&ts[1], TIME_UTC);
    if(w!=NULL) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, w);
    int ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_r2, nnls_wp, nnls_zz, nnls_index);
    //timespec_get(&ts[2], TIME_UTC);
    if(ret>1) break;
    if(verbose>2) printf("  -> r2=%g\n", nnls_r2);
    if(nnls_r2<bestR2) {
      bestR2=nnls_r2;
      bestDeltaT=deltaT;
      for(int n=0; n<nnls_n; n++) la[n]=nnls_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;
  }
 
  free(nnls_mat);
 
  /* Compute fitted y[] */
  if(yfit!=NULL) {
    if(verbose>1) printf("computing yfit[]\n");
    for(int m=0; m<nnls_m; m++) {
      yfit[m]=0.0;
      for(int n=0; n<nnls_n; n++)
        if(la[n]>0.0)
          yfit[m]+=la[n]*(x[m]+bestDeltaT)*exp(-lk[n]*(x[m]+bestDeltaT));
    }
  }
 
  /* 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;
}
/*****************************************************************************/
 
/*****************************************************************************/