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 {
113 double r1, r2, s;
114 r1=log10(kMin); r2=log10(kMax); s=(r2-r1)/(double)(fNr-1);
115 if(verbose>3) printf(" r1 := %g\n r2 := %g\n s := %g\n", r1, r2, s);
116 for(int bi=0; bi<fNr; bi++) lk[bi]=pow(10.0, (double)bi*s+r1);
117 }
118
119
120 /* Allocate memory required by NNLS */
121 if(verbose>2) printf("allocating memory for LLSQ\n");
122 int nnls_n, nnls_m, nnls_index[fNr];
123 double *nnls_a[fNr], nnls_x[fNr], nnls_wp[fNr], *nnls_b, *nnls_zz, *nnls_mat, *dptr, nnls_r2;
124 nnls_n=fNr; nnls_m=sNr;
125 nnls_mat=(double*)malloc(((nnls_n+2)*nnls_m)*sizeof(double));
126 if(nnls_mat==NULL) {
127 free(localp);
128 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
129 return(TPCERROR_OUT_OF_MEMORY);
130 }
131 dptr=nnls_mat;
132 for(int n=0; n<nnls_n; n++) {nnls_a[n]=dptr; dptr+=nnls_m;}
133 nnls_b=dptr; dptr+=nnls_m; nnls_zz=dptr;
134
135 /* Set initial delay time to the middle of the given range */
136 double initDeltaT=dtMin;
137 if(dtMax>dtMin) initDeltaT=0.5*(dtMin+dtMax);
138 if(verbose>3) printf("initDeltaT := %g\n", initDeltaT);
139
140 /* LLSQ fit with initial delay time */
141 if(verbose>2) printf("LLSQ fitting\n");
142 double bestDeltaT=initDeltaT;
143 double bestR2=nan("");
144 double deltaT=initDeltaT;
145 if(verbose>2) printf(" deltaT=%g\n", deltaT);
146 if(verbose>6) printf("filling data matrix\n");
147 /* Fill NNLS B array with measured y values */
148 for(int m=0; m<nnls_m; m++) nnls_b[m]=y[m];
149 /* Fill NNLS A matrix with basis functions, calculated in place */
150 int ret=0;
151 double p[3]={deltaT, 1.0, 0.0};
152 for(int n=0; n<nnls_n && !ret; n++) {
153 p[2]=lk[n];
154 if(x2!=NULL) // frame start and end times available
155 ret=mfEvalFrameY("dmsurge", 3, p, nnls_m, x, x2, nnls_a[n], 0);
156 else
157 ret=mfEvalY("dmsurge", 3, p, nnls_m, x, nnls_a[n], 0);
158 }
159 if(ret) {
160 free(nnls_mat); free(localp);
161 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
162 return(TPCERROR_INVALID_VALUE);
163 }
164 /* Apply weights, if given */
165 if(w!=NULL) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, w);
166 /* NNLS */
167 if(verbose>6) printf("applying NNLS\n");
168 ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_r2, nnls_wp, nnls_zz, nnls_index);
169 if(ret>1) {
170 free(nnls_mat); free(localp);
171 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_SOLUTION);
172 return(TPCERROR_NO_SOLUTION);
173 }
174 if(verbose>2) printf(" -> r2=%g\n", nnls_r2);
175 bestR2=nnls_r2;
176 /* Copy a[] values */
177 for(int n=0; n<nnls_n; n++) la[n]=nnls_x[n];
178
179 /* Try moving function to the left with delay steps */
180 deltaT=initDeltaT-dtStep;
181 while(dtStep>0.0 && deltaT>=dtMin) {
182 if(verbose>2) printf(" deltaT=%g\n", deltaT);
183 for(int m=0; m<nnls_m; m++) nnls_b[m]=y[m];
184 p[0]=deltaT;
185 for(int n=0; n<nnls_n && !ret; n++) {
186 p[2]=lk[n];
187 if(x2!=NULL) // frame start and end times available
188 ret=mfEvalFrameY("dmsurge", 3, p, nnls_m, x, x2, nnls_a[n], 0);
189 else
190 ret=mfEvalY("dmsurge", 3, p, nnls_m, x, nnls_a[n], 0);
191 }
192 if(ret) {
193 free(nnls_mat); free(localp);
194 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
195 return(TPCERROR_INVALID_VALUE);
196 }
197 if(w!=NULL) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, w);
198 int ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_r2, nnls_wp, nnls_zz, nnls_index);
199 if(ret>1) break;
200 if(verbose>2) printf(" -> r2=%g\n", nnls_r2);
201 if(nnls_r2<bestR2) {
202 bestR2=nnls_r2;
203 bestDeltaT=deltaT;
204 for(int n=0; n<nnls_n; n++) la[n]=nnls_x[n];
205 } //else if(nnls_r2>bestR2) break;
206 deltaT-=dtStep;
207 }
208 /* Try moving function to the right with delay steps */
209 deltaT=initDeltaT+dtStep;
210 while(dtStep>0.0 && deltaT<=dtMax) {
211 if(verbose>2) printf(" deltaT=%g\n", deltaT);
212 for(int m=0; m<nnls_m; m++) nnls_b[m]=y[m];
213 p[0]=deltaT;
214 //struct timespec ts[3];
215 //timespec_get(&ts[0], TIME_UTC);
216 for(int n=0; n<nnls_n && !ret; n++) {
217 p[2]=lk[n];
218 if(x2!=NULL) // frame start and end times available
219 ret=mfEvalFrameY("dmsurge", 3, p, nnls_m, x, x2, nnls_a[n], 0);
220 else
221 ret=mfEvalY("dmsurge", 3, p, nnls_m, x, nnls_a[n], 0);
222 }
223 if(ret) {
224 free(nnls_mat); free(localp);
225 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
226 return(TPCERROR_INVALID_VALUE);
227 }
228 //timespec_get(&ts[1], TIME_UTC);
229 if(w!=NULL) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, w);
230 int ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_r2, nnls_wp, nnls_zz, nnls_index);
231 //timespec_get(&ts[2], TIME_UTC);
232 if(ret>1) break;
233 if(verbose>2) printf(" -> r2=%g\n", nnls_r2);
234 if(nnls_r2<bestR2) {
235 bestR2=nnls_r2;
236 bestDeltaT=deltaT;
237 for(int n=0; n<nnls_n; n++) la[n]=nnls_x[n];
238 } //else if(nnls_r2>bestR2) break;
239 //printf(" time_to-fill=%g time_to_nnls=%g\n", timespecDifference(&ts[1], &ts[0]), timespecDifference(&ts[2], &ts[1]));
240 deltaT+=dtStep;
241 }
242
243 free(nnls_mat);
244
245 /* Compute fitted y[] */
246 if(yfit!=NULL) {
247 if(verbose>1) printf("computing yfit[]\n");
248 for(int m=0; m<nnls_m; m++) {
249 yfit[m]=0.0;
250 for(int n=0; n<nnls_n; n++)
251 if(la[n]>0.0)
252 yfit[m]+=la[n]*(x[m]+bestDeltaT)*exp(-lk[n]*(x[m]+bestDeltaT));
253 }
254 }
255
256 /* Copy a[] and b[] if requested */
257 if(a!=NULL) for(int bi=0; bi<fNr; bi++) a[bi]=la[bi];
258 if(k!=NULL) for(int bi=0; bi<fNr; bi++) k[bi]=lk[bi];
259 free(localp);
260
261 if(dtEst!=NULL) *dtEst=bestDeltaT;
262
263 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
264 return TPCERROR_OK;
265}
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
nnls
int nnls(double **a, int m, int n, double *b, double *x, double *rnorm, double *wp, double *zzp, int *indexp)
Definition nnls.c:43
nnlsWght
int nnlsWght(int N, int M, double **A, double *b, double *weight)
Definition nnls.c:259
statusSet
void statusSet(TPCSTATUS *s, const char *func, const char *srcfile, int srcline, tpcerror error)
Definition statusmsg.c:142
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.15.0