bf_dexp.c File Reference

BFM for the sum of decaying exponential functions. More...

#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"

Go to the source code of this file.

Functions
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	spectralKRange (double *k, double *a, const int n, double *kmin, double *kmax, TPCSTATUS *status)
int	spectralBFNr (double *k, double *a, const int n)
int	spectralBFExtract (double *k, double *a, const int n, double *ke, double *ae, const int ne)

Detailed Description

BFM for the sum of decaying exponential functions.

Definition in file bf_dexp.c.

Function Documentation

spectralBFExtract()

int spectralBFExtract	(	double *	k,
		double *	a,
		const int	n,
		double *	ke,
		double *	ae,
		const int	ne )

After spectral x,y-data fitting to the sum of decaying exponential functions, extract the zero-separated (positive a) basis functions we actually got.

The maximum number of basis function parameters to be extracted must be specified; if more are found, then the ones with smallest a are averaged. If less are found, then the rest of output arrays are filled with zeroes, and the returned actual number is smaller than the requested number.

Returns

The number of extracted basis functions, or 0 in case of an error.

Author

Vesa Oikonen

See also

spectralDExp, spectralBFNr, spectralKRange

Parameters

k	Pointer to k array of length n; must be sorted to either order, as those usually are. Contents are not modified.
a	Pointer to a array of length n; not modified.
n	Length of k and a arrays.
ke	Pointer to an allocated array of extracted k values of length ne.
ae	Pointer to an allocated array of extracted a values of length ne; not modified.
ne	Length of k and a arrays, and the maximum number of parameters to extract.

Definition at line 227 of file bf_dexp.c.

  {
  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);
}

spectralBFNr()

int spectralBFNr	(	double *	k,
		double *	a,
		const int	n )

After spectral x,y-data fitting to the sum of decaying exponential functions, determine the number of zero-separated (positive a) basis functions we actually got.

Returns

The number of basis functions, or 0 in case of an error.

Author

Vesa Oikonen

See also

spectralDExp, spectralBFExtract, spectralKRange

Parameters

k	Pointer to k array of length n; does not need to be in specific order, although those usually are. Contents are not modified.
a	Pointer to a array of length n; not modified.
n	Length of k and a arrays.

Definition at line 193 of file bf_dexp.c.

  {
  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);
}

Referenced by spectralBFExtract().

spectralDExp()

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 )

Spectral x,y-data fitting to the sum of decaying exponential functions.

f(x) = a1*exp(-k1*x) + a2*exp(-k2*x) + a3*exp(-k3*x) + ... where aN>=0 and kN>0.

Todo

Add tests.

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

See also

spectralKRange, spectralBFNr, nnls, nnlsWght, tacExtractRange

Parameters

x	Pointer to TAC x data (not modified). Data must be sorted by increasing x. Negative or missing x values are not allowed.
y	Pointer to TAC y data (not modified). Data must be sorted by increasing x. Missing y values are not allowed.
w	Pointer to TAC sample weights (not modified). Enter NULL if not needed.
snr	Number of samples in x[] and y[]; at least 3.
kmin	Minimum of k (set based on sample range and unit); must be >0.
kmax	Maximum of k (set based on sample range and unit); must be >kmin.
bnr	Number of basis functions to calculate; also the length of arrays k[] and a[]; at least 4.
k	Pointer to array of length bnr for k values, filled in by this function.
a	Pointer to array of length bnr for a values, filled in by this function; notice that most values may be set to zero.
yfit	Pointer to array for fitted y values. Enter NULL if not needed.
status	Pointer to status data; enter NULL if not needed

Definition at line 28 of file bf_dexp.c.

  {
  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;
}

spectralKRange()

int spectralKRange	(	double *	k,
		double *	a,
		const int	n,
		double *	kmin,
		double *	kmax,
		TPCSTATUS *	status )

After spectral x,y-data fitting to the sum of decaying exponential functions, determine the min and max of k values with positive a value.

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

See also

spectralDExp, spectralBFNr

Parameters

k	Pointer to k array of length n; does not need to be in specific order, although those usually are. Contents are not modified.
a	Pointer to a array of length n; not modified.
n	Length of k and a arrays.
kmin	Pointer min k value. Enter NULL if not needed.
kmax	Pointer max k value. Enter NULL if not needed.
status	Pointer to status data; enter NULL if not needed

Definition at line 145 of file bf_dexp.c.

  {
  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;
}