TPCCLIB
Loading...
Searching...
No Matches
Functions
bf_dms.c File Reference

BFM for the sum of surge functions with delay. 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 "tpcfunc.h"
#include "tpcbfm.h"

Go to the source code of this file.

Functions

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)

Detailed Description

BFM for the sum of surge functions with delay.

Definition in file bf_dms.c.

Function Documentation

◆ spectralDMSurge()

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 )

Spectral x,y-data fitting to the sum of surge functions and delay time.

f(x) = a1*x*exp(-k1*x) + a2*x*exp(-k2*x) + a3*x*exp(-k3*x) + ... where a and k parameters are larger than zero.

Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
See also
spectralDExp, spectralKRange, spectralBFNr
Parameters
xPointer to TAC x (sample times or frame mid times), or x1 (frame start times) when x2 are given, too. Data is not modified in this function. Data must be sorted by increasing x. Negative or missing x values are not allowed.
x2Pointer to TAC x2 data (frame end times), or NULL if frame start and end times are not available. Data is not modified.
yPointer to TAC y data (not modified). Data must be sorted by increasing x. Missing y values are not allowed.
wPointer to TAC sample weights (not modified). Enter NULL if not needed.
sNrNumber of samples in x[], y[], and possibly x2[] and w[]; at least 4.
kMinMinimum of k; 0<kMin<kMax.
kMaxMaximum of k; 0<kMin<kMax.
fNrNumber of basis functions to calculate; also the length of arrays k[] and a[]; at least 4.
dtMinMin delay time. To fix delay time enter dtMin=dtMax, and dtStep=0.
dtMaxMax delay time. To fix delay time enter dtMin=dtMax, and dtStep=0.
dtStepDelay time step size. Enter zero to use delay time fixed to dtMin&dtMax.
kPointer to array of length fNr for k values, filled in by this function. Enter NULL if not needed.
aPointer to array of length fNr for a values, filled in by this function; notice that most values may be set to zero. Enter NULL if not needed.
dtEstPointer to delay time estimate; NULL if not needed.
yfitPointer to array for fitted y values. Enter NULL if not needed.
statusPointer to status data; enter NULL if not needed

Definition at line 27 of file bf_dms.c.

66 {
67 int verbose=0; if(status!=NULL) verbose=status->verbose;
68 if(verbose>0) printf("%s()\n", __func__);
69 else if(verbose>1) {
70 printf("%s(x, ", __func__);
71 if(x2==NULL) printf("null, y, "); else printf("x2, y, ");
72 if(w==NULL) printf("null"); else printf("w");
73 printf(" , %d, %.1e, %.1e, %d, %g, %g, %g, *k, *a, *dtEst, *yfit\n",
74 sNr, kMin, kMax, fNr, dtMin, dtMax, dtStep);
75 }
76 if(x==NULL || y==NULL) {
77 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
78 return(TPCERROR_NO_DATA);
79 }
80 if(sNr<4) {
81 if(verbose>1) fprintf(stderr, "invalid number of samples\n");
82 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
83 return(TPCERROR_TOO_FEW);
84 }
85 if(fNr<4) {
86 if(verbose>1) fprintf(stderr, "invalid number of functions\n");
87 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
88 return(TPCERROR_TOO_FEW);
89 }
90 if(!(kMin>0.0) || !(kMax>kMin)) {
91 if(verbose>1) fprintf(stderr, "invalid k range\n");
92 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
93 return(TPCERROR_INVALID_VALUE);
94 }
95 double dtRange=dtMax-dtMin;
96 if(!(dtRange>=0.0) || (dtRange>0.0 && !(dtStep<0.2*dtRange))) {
97 if(verbose>1) fprintf(stderr, "invalid delay time settings\n");
98 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
99 return(TPCERROR_INVALID_VALUE);
100 }
101
102 /* Allocate memory or set pointers for local a[] and k[] */
103 double *lk, *la, *localp;
104 localp=(double*)malloc(sizeof(double)*fNr*2);
105 if(localp==NULL) {
106 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
107 return(TPCERROR_OUT_OF_MEMORY);
108 }
109 lk=localp; la=localp+fNr;
110
111 if(verbose>2) printf("computing k values\n");
112#if(1)
113 {
114 double r1, r2, s;
115 r1=log(kMin); r2=log(kMax); s=(r2-r1)/(double)(fNr-1);
116 if(verbose>3) printf(" r1 := %g\n r2 := %g\n s := %g\n", r1, r2, s);
117 for(int bi=0; bi<fNr; bi++) lk[bi]=exp((double)bi*s+r1);
118 }
119#else
120 {
121 double r1, r2, s;
122 r1=log10(kMin); r2=log10(kMax); s=(r2-r1)/(double)(fNr-1);
123 if(verbose>3) printf(" r1 := %g\n r2 := %g\n s := %g\n", r1, r2, s);
124 for(int bi=0; bi<fNr; bi++) lk[bi]=pow(10.0, (double)bi*s+r1);
125 }
126#endif
127
128 if(verbose>2) printf("computing the range of delay times to test\n");
129
130
131 /* Allocate memory required by NNLS */
132 if(verbose>2) printf("allocating memory for LLSQ\n");
133 NNLSQDATA lsq;
134 nnlsqDataInit(&lsq);
135 int ret=0;
136 if((ret=nnlsqDataAllocate(&lsq, fNr, sNr))!=TPCERROR_OK) {
137 free(localp);
138 statusSet(status, __func__, __FILE__, __LINE__, ret);
139 return(ret);
140 }
141
142 /* Set initial delay time to the middle of the given range */
143 double initDeltaT=dtMin;
144 if(dtMax>dtMin) initDeltaT=0.5*(dtMin+dtMax);
145 if(verbose>3) printf("initDeltaT := %g\n", initDeltaT);
146
147 /* LLSQ fit with initial delay time */
148 if(verbose>2) printf("LLSQ fitting\n");
149 double bestDeltaT=initDeltaT;
150 double bestR2=nan("");
151 double deltaT=initDeltaT;
152 if(verbose>2) printf(" deltaT=%g\n", deltaT);
153 if(verbose>6) printf("filling data matrix\n");
154 /* Fill NNLS B array with measured y values */
155 for(int m=0; m<lsq.m; m++) lsq.b[m]=y[m];
156 /* Fill NNLS A matrix with basis functions, calculated in place */
157 double p[3]={deltaT, 1.0, 0.0};
158 for(int n=0, ret=0; n<lsq.n && !ret; n++) {
159 p[2]=lk[n];
160 if(x2!=NULL) // frame start and end times available
161 ret=mfEvalFrameY("dmsurge", 3, p, lsq.m, x, x2, lsq.a[n], 0);
162 else
163 ret=mfEvalY("dmsurge", 3, p, lsq.m, x, lsq.a[n], 0);
164 }
165 if(ret) {
166 nnlsqDataFree(&lsq); free(localp);
167 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
168 return(TPCERROR_INVALID_VALUE);
169 }
170 /* Apply weights, if given */
171 if(w!=NULL) nnlsqDataWght(&lsq, w);
172 /* NNLS */
173 if(verbose>6) printf("applying NNLS\n");
174 lsq.iternr=3*lsq.n;
175 ret=nnlsq(&lsq, verbose-1);
176 if(ret>1) {
177 nnlsqDataFree(&lsq); free(localp);
178 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_SOLUTION);
179 return(TPCERROR_NO_SOLUTION);
180 }
181 if(verbose>2) printf(" -> r2=%g iterNr=%d\n", lsq.rnorm, lsq.iternr);
182 bestR2=lsq.rnorm;
183 /* Copy a[] values */
184 for(int n=0; n<lsq.n; n++) la[n]=lsq.x[n];
185
186 /* Try moving function to the left with delay steps */
187 deltaT=initDeltaT-dtStep;
188 while(dtStep>0.0 && deltaT>=dtMin) {
189 if(verbose>2) printf(" deltaT=%g\n", deltaT);
190 for(int m=0; m<lsq.m; m++) lsq.b[m]=y[m];
191 p[0]=deltaT;
192 for(int n=0, ret=0; n<lsq.n && !ret; n++) {
193 p[2]=lk[n];
194 if(x2!=NULL) // frame start and end times available
195 ret=mfEvalFrameY("dmsurge", 3, p, lsq.m, x, x2, lsq.a[n], 0);
196 else
197 ret=mfEvalY("dmsurge", 3, p, lsq.m, x, lsq.a[n], 0);
198 }
199 if(ret) {
200 nnlsqDataFree(&lsq); free(localp);
201 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
202 return(TPCERROR_INVALID_VALUE);
203 }
204 if(w!=NULL) nnlsqDataWght(&lsq, w);
205 lsq.iternr=3*lsq.n;
206 int ret=nnlsq(&lsq, verbose-1);
207 if(ret>1) break;
208 if(verbose>2) printf(" -> r2=%g iterNr=%d\n", lsq.rnorm, lsq.iternr);
209 if(lsq.rnorm<bestR2) {
210 bestR2=lsq.rnorm;
211 bestDeltaT=deltaT;
212 for(int n=0; n<lsq.n; n++) la[n]=lsq.x[n];
213 } //else if(nnls_r2>bestR2) break;
214 deltaT-=dtStep;
215 }
216 /* Try moving function to the right with delay steps */
217 deltaT=initDeltaT+dtStep;
218 while(dtStep>0.0 && deltaT<=dtMax) {
219 if(verbose>2) printf(" deltaT=%g\n", deltaT);
220 for(int m=0; m<lsq.m; m++) lsq.b[m]=y[m];
221 p[0]=deltaT;
222 //struct timespec ts[3];
223 //timespec_get(&ts[0], TIME_UTC);
224 for(int n=0, ret=0; n<lsq.n && !ret; n++) {
225 p[2]=lk[n];
226 if(x2!=NULL) // frame start and end times available
227 ret=mfEvalFrameY("dmsurge", 3, p, lsq.m, x, x2, lsq.a[n], 0);
228 else
229 ret=mfEvalY("dmsurge", 3, p, lsq.m, x, lsq.a[n], 0);
230 }
231 if(ret) {
232 nnlsqDataFree(&lsq); free(localp);
233 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
234 return(TPCERROR_INVALID_VALUE);
235 }
236 //timespec_get(&ts[1], TIME_UTC);
237 if(w!=NULL) nnlsqDataWght(&lsq, w);
238 lsq.iternr=3*lsq.n;
239 int ret=nnlsq(&lsq, verbose-1);
240 //timespec_get(&ts[2], TIME_UTC);
241 if(ret>1) break;
242 if(verbose>2) printf(" -> r2=%g iterNr=%d\n", lsq.rnorm, lsq.iternr);
243 if(lsq.rnorm<bestR2) {
244 bestR2=lsq.rnorm;
245 bestDeltaT=deltaT;
246 for(int n=0; n<lsq.n; n++) la[n]=lsq.x[n];
247 } //else if(nnls_r2>bestR2) break;
248 //printf(" time_to-fill=%g time_to_nnls=%g\n", timespecDifference(&ts[1], &ts[0]), timespecDifference(&ts[2], &ts[1]));
249 deltaT+=dtStep;
250 }
251
252 /* Compute fitted y[] */
253 if(yfit!=NULL) {
254 if(verbose>1) printf("computing yfit[]\n");
255 for(int m=0; m<lsq.m; m++) {
256 yfit[m]=0.0;
257 for(int n=0; n<lsq.n; n++)
258 if(la[n]>0.0)
259 yfit[m]+=la[n]*(x[m]+bestDeltaT)*exp(-lk[n]*(x[m]+bestDeltaT));
260 }
261 }
262
263 nnlsqDataFree(&lsq);
264
265 /* Copy a[] and b[] if requested */
266 if(a!=NULL) for(int bi=0; bi<fNr; bi++) a[bi]=la[bi];
267 if(k!=NULL) for(int bi=0; bi<fNr; bi++) k[bi]=lk[bi];
268 free(localp);
269
270 if(dtEst!=NULL) *dtEst=bestDeltaT;
271
272 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
273 return TPCERROR_OK;
274}
mfEvalY
int mfEvalY(const char *fid, const int parNr, const double *p, const int sampleNr, const double *x, double *y, const int verbose)
Definition func.c:26
mfEvalFrameY
int mfEvalFrameY(const char *fid, const int parNr, const double *p, const int sampleNr, const double *x1, const double *x2, double *y, const int verbose)
Definition func.c:790
nnlsqDataWght
int nnlsqDataWght(NNLSQDATA *d, double *weight)
Definition nnlsq.c:513
nnlsqDataFree
void nnlsqDataFree(NNLSQDATA *d)
Definition nnlsq.c:379
nnlsqDataAllocate
int nnlsqDataAllocate(NNLSQDATA *d, const int n, const int m)
Definition nnlsq.c:412
nnlsqDataInit
void nnlsqDataInit(NNLSQDATA *d)
Definition nnlsq.c:358
nnlsq
int nnlsq(NNLSQDATA *d, int verbose)
Definition nnlsq.c:55
statusSet
void statusSet(TPCSTATUS *s, const char *func, const char *srcfile, int srcline, tpcerror error)
Definition statusmsg.c:142
NNLSQDATA::rnorm
double rnorm
Definition tpclinopt.h:136
NNLSQDATA::x
double * x
Definition tpclinopt.h:123
NNLSQDATA::b
double * b
Definition tpclinopt.h:121
NNLSQDATA::a
double ** a
Definition tpclinopt.h:118
TPCSTATUS::verbose
int verbose
Verbose level, used by statusPrint() etc.
Definition tpcextensions.h:242
TPCERROR_INVALID_VALUE
@ TPCERROR_INVALID_VALUE
Invalid value.
Definition tpcextensions.h:220
TPCERROR_NO_SOLUTION
@ TPCERROR_NO_SOLUTION
No solution.
Definition tpcextensions.h:227
TPCERROR_OUT_OF_MEMORY
@ TPCERROR_OUT_OF_MEMORY
Cannot allocate memory.
Definition tpcextensions.h:193
TPCERROR_OK
@ TPCERROR_OK
No error.
Definition tpcextensions.h:191
TPCERROR_NO_DATA
@ TPCERROR_NO_DATA
File contains no data.
Definition tpcextensions.h:203
TPCERROR_TOO_FEW
@ TPCERROR_TOO_FEW
File contains too few samples.
Definition tpcextensions.h:208
Generated on for TPCCLIB by doxygen 1.17.0