bf_dexp.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 "tpcbfm.h"
/*****************************************************************************/
 
/*****************************************************************************/
int spectralDExp(
  const double *x,
  const double *y,
  double *w,
  const int snr,
  const double kmin,
  const double kmax,
  const int bnr,
  double *k,
  double *a,
  double *yfit,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(x==NULL || y==NULL || k==NULL || a==NULL || snr<1 || bnr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(snr<3 || bnr<4) {
    if(verbose>1) fprintf(stderr, "invalid sample or BF number\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    return TPCERROR_TOO_FEW;
  }
  if(!(kmin>0.0) || kmin>=kmax) { // range cannot be computed otherwise
    if(verbose>1) fprintf(stderr, "invalid k range\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  int ret;
 
  if(verbose>1) printf("computing k values\n");
  {
    double a, b, c;
    a=log10(kmin); b=log10(kmax); c=(b-a)/(double)(bnr-1);
    if(verbose>3) printf(" a := %g\n b := %g\n c := %g\n", a, b, c);
    for(int bi=0; bi<bnr; bi++) k[bi]=pow(10.0, (double)bi*c+a);
  }
 
  /* Allocate memory required by NNLS */
  int     nnls_n, nnls_m, n, m, nnls_index[bnr];
  double *nnls_a[bnr], *nnls_b, *nnls_zz, nnls_x[bnr], *nnls_mat,
          nnls_wp[bnr], *dptr, nnls_rnorm;
  if(verbose>1) printf("allocating memory for NNLS\n");
  nnls_n=bnr; nnls_m=snr;
  nnls_mat=(double*)malloc(((nnls_n+2)*nnls_m)*sizeof(double));
  if(nnls_mat==NULL) {
    if(verbose>1) fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return TPCERROR_OUT_OF_MEMORY;
  }
  for(n=0, dptr=nnls_mat; n<nnls_n; n++) {nnls_a[n]=dptr; dptr+=nnls_m;}
  nnls_b=dptr; dptr+=nnls_m; nnls_zz=dptr;
 
  if(verbose>1) printf("filling NNLS data matrix\n");
  /* Fill NNLS B array with measured y values */
  for(m=0; m<nnls_m; m++) nnls_b[m]=y[m];
 
  /* Fill NNLS A matrix with basis functions, calculated in place */
  for(n=0; n<nnls_n; n++)
    for(m=0; m<nnls_m; m++)
      nnls_a[n][m]=exp(-k[n]*x[m]);
 
  /* Apply weights, if given */
  if(w!=NULL) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, w);
 
  /* NNLS */
  if(verbose>1) printf("applying NNLS\n");
  ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm, nnls_wp, nnls_zz, nnls_index);
  if(ret>1) {
    if(verbose>1) fprintf(stderr, "NNLS solution not possible.\n");
    free(nnls_mat);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_SOLUTION);
    return TPCERROR_NO_SOLUTION;
  }
  /* Copy a[] values */
  if(a!=NULL) for(n=0; n<nnls_n; n++) a[n]=nnls_x[n];
 
  /* Compute fitted y[] */
  if(yfit!=NULL) {
    if(verbose>1) printf("computing yfit[]\n");
    for(m=0; m<nnls_m; m++) {
      yfit[m]=0.0;
      for(n=0; n<nnls_n; n++)
        if(nnls_x[n]>0.0)
          yfit[m]+=nnls_x[n]*exp(-k[n]*x[m]);
    }
  }
  free(nnls_mat);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int spectralKRange(
  double *k,
  double *a,
  const int n,
  double *kmin,
  double *kmax,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(k[], a[], %d, *kmin, *kmax)\n", __func__, n);
  if(k==NULL || a==NULL || n<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  double max=nan(""), min=nan("");
  for(int i=0; i<n; i++) if(a[i]>0.0 && isfinite(k[i])) {
    if(isnan(max) || k[i]>max) max=k[i];
    if(isnan(min) || k[i]<min) min=k[i];
  }
  if(verbose>2) {
    printf("  kmin := %g\n", min);
    printf("  kmax := %g\n", max);
  }
  if(kmin!=NULL) *kmin=min;
  if(kmax!=NULL) *kmax=max;
  if(isnan(min) || isnan(max)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_VALUE);
    return TPCERROR_NO_VALUE;
  }
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int spectralBFNr(
  double *k,
  double *a,
  const int n
) {
  if(k==NULL || a==NULL || n<1) return(0);
  int bfn=0, i=0, ii=0;
  do {
    i=doubleCSpanPositives(a+ii, n-ii);
    ii+=i; if(n-ii<1) break;
    i=doubleSpanPositives(a+ii, n-ii);
    ii+=i; if(i>0) bfn++;
  } while(ii<n);
  return(bfn);
}
/*****************************************************************************/
 
/*****************************************************************************/
int spectralBFExtract(
  double *k,
  double *a,
  const int n,
  double *ke,
  double *ae,
  const int ne
) {
  if(k==NULL || a==NULL || n<1 || ke==NULL || ae==NULL || ne<1) return(0);
 
  /* First, extract all zero-separated basis function clusters */
  double *buf=NULL, *ok, *oa;
  int bfn=spectralBFNr(k, a, n); if(bfn==0) return(0);
  if(bfn>ne) {
    /* we need to allocate memory locally */
    buf=(double*)malloc(2*bfn*sizeof(double)); if(buf==NULL) return(0);
    ok=buf; oa=buf+bfn;
  } else {
    /* we just set pointers to user-provided arrays */
    ok=ke; oa=ae;
  }
  {
    int nn=0, i=0, ii=0;
    do {
      i=doubleCSpanPositives(a+ii, n-ii);
      ii+=i; if(n-ii<1) break;
      i=doubleSpanPositives(a+ii, n-ii);
      if(i==1) { // cluster contains just one positive a
        oa[nn]=a[ii]; ok[nn]=k[ii];
        nn++;
      } else if(i>1) { // cluster contains more than one positive a
        /* get sum of a values in the cluster */
        oa[nn]=doubleSum(a+ii, i);
        /* get weighted mean of k values in the cluster */
        ok[nn]=doubleWMean(k+ii, a+ii, i);
        nn++;
      }
      ii+=i;
    } while(ii<n && nn<bfn);
    if(0) {
      printf("\n\tCluster a and k values\n");
      for(int i=0; i<bfn; i++) printf("\t%g\t%g\n", oa[i], ok[i]);
      printf("\n"); fflush(stdout);
    }
  }
  /* If the number needs not to be reduced, then we are ready */
  if(ne>=bfn) {
    for(int i=bfn; i<ne; i++) ke[i]=ae[i]=0.0;
    if(buf!=NULL) free(buf); 
    return(bfn);
  }
 
  while(ne<bfn) {
    /* Search the smallest a parameter */
    int ami=doubleAbsMinIndex(oa, bfn);
    if(0) printf("smallest |a[%d]|=%g\n", ami, oa[ami]);
    /* Combine it with previous or next cluster? */
    int ci;
    if(ami==0) ci=1;
    else if(ami==bfn-1) ci=bfn-2;
    else {if(fabs(oa[ami-1])<fabs(oa[ami+1])) ci=ami-1; else ci=ami+1;}
    if(0) printf("combining clusters %d and %d\n", 1+ami, 1+ci);
    /* Replace the cluster with combination of two clusters */
    ok[ci]=(oa[ami]*ok[ami]+oa[ci]*ok[ci])/(oa[ami]+oa[ci]);
    oa[ci]+=oa[ami];
    /* Delete the cluster */
    for(int i=ami+1; i<bfn; i++) oa[i-1]=oa[i];
    for(int i=ami+1; i<bfn; i++) ok[i-1]=ok[i];
    bfn--;
    if(0) {
      printf("\n\tAfter reduction, cluster a and k values\n");
      for(int i=0; i<bfn; i++) printf("\t%g\t%g\n", oa[i], ok[i]);
      printf("\n"); fflush(stdout);
    }
  }
  for(int i=0; i<ne; i++) ae[i]=oa[i];
  for(int i=0; i<ne; i++) ke[i]=ok[i];
  if(buf!=NULL) free(buf); 
 
  return(bfn);
}
/*****************************************************************************/
 
/*****************************************************************************/