TPCCLIB
Loading...
Searching...
No Matches
Functions
tpcbfm.h File Reference

Header file for libtpcbfm. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include "tpcextensions.h"
#include "tpccm.h"
#include "tpctac.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)
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 bfm1TCM (TAC *input, TAC *tissue, const int bfNr, const double k2min, const double k2max, const int distr, TAC *bf, TPCSTATUS *status)
int bfmSRTM (double *t, double *cri, const int n, const int bfNr, const double t3min, const double t3max, TAC *bf, TPCSTATUS *status)

Detailed Description

Header file for libtpcbfm.

Header file for template library libtpcbfm.

Author
Vesa Oikonen
Copyright
(c) Turku PET Centre

Definition in file tpcbfm.h.

Function Documentation

◆ bfm1TCM()

int bfm1TCM ( TAC * input,
TAC * tissue,
int bfNr,
const double k2min,
const double k2max,
const int distr,
TAC * bf,
TPCSTATUS * status )

Calculate set of basis functions for generic radiowater model.

Todo
Add tests.
Returns
enum tpcerror (TPCERROR_OK when successful).
See also
bfmSRTM
Parameters
inputPointer to blood input TAC (not modified). If frame start and end times are available (->isframes!=0), then BTAC integral is calculated and used as input.
tissuePointer to PET TTAC data, which must contain the sample (frame) times.
bfNrNr of basis functions to calculate.
k2minMinimum of k2 (sec-1 or min-1, corresponding to TAC time units). If set to a negative value, |k2min| is used, but basis function with k2=0 is included.
k2maxMaximum of k2 (sec-1 or min-1, corresponding to TAC time units).
distrDistribution of k2 values: 0=log-based; 1=even.
bfPointer to output TAC, to be allocated and filled with basis functions in here.
statusPointer to status data; enter NULL if not needed.

Definition at line 27 of file bf_1tcm.c.

46 {
47 int verbose=0; if(status!=NULL) verbose=status->verbose;
48 if(verbose>0) printf("%s(itac, ttac, %d, %g, %g, %d, bf)\n", __func__, bfNr, k2min, k2max, distr);
49 if(input==NULL || tissue==NULL || bf==NULL) {
50 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
51 return TPCERROR_FAIL;
52 }
53 if(input->sampleNr<1 || input->tacNr<1 || tissue->sampleNr<1) {
54 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
55 return TPCERROR_NO_DATA;
56 }
57 if(input->sampleNr<3 || bfNr<2) {
58 if(verbose>1) printf("invalid sample or BF number\n");
59 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
60 return TPCERROR_TOO_FEW;
61 }
62 if(fabs(k2min)<1.0E-10 || fabs(k2min)>=k2max) { // range cannot be computed otherwise
63 if(verbose>1) printf("invalid k2 range\n");
64 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
65 return TPCERROR_INVALID_VALUE;
66 }
67
68 int ret;
69
70 /* Check units */
71 if(status==NULL || !status->forgiving) {
72 if(!unitIsTime(input->tunit) || input->tunit!=tissue->tunit) {
73 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNKNOWN_UNIT);
74 return TPCERROR_UNKNOWN_UNIT;
75 }
76 }
77
78 /* Allocate memory for basis functions */
79 if(verbose>1) printf("allocating memory for basis functions\n");
80 tacFree(bf);
81 ret=tacAllocate(bf, tissue->sampleNr, bfNr);
82 if(ret!=TPCERROR_OK) {
83 statusSet(status, __func__, __FILE__, __LINE__, ret);
84 return ret;
85 }
86 /* Copy and set information fields */
87 bf->tacNr=bfNr; bf->sampleNr=tissue->sampleNr;
88 bf->format=TAC_FORMAT_TSV_UK;
89 bf->isframe=tissue->isframe; tacXCopy(tissue, bf, 0, tissue->sampleNr-1);
90 bf->tunit=tissue->tunit; bf->cunit=tissue->cunit;
91 for(int bi=0; bi<bf->tacNr; bi++) sprintf(bf->c[bi].name, "B%5.5d", bi+1);
92 /* Compute the range of k2 values to size fields */
93 if(verbose>1) printf("computing k2 values\n");
94 int zero=0; if(k2min<0.0) {zero=1; bf->c[0].size=0.0;}
95 if(distr==0) { // printf(" log-based distribution\n");
96 double a, b, c;
97 a=log10(fabs(k2min)); b=log10(k2max); c=(b-a)/(double)(bfNr-1-zero);
98 for(int bi=zero; bi<bf->tacNr; bi++) {
99 bf->c[bi].size=pow(10.0, (double)(bi-zero)*c+a);
100 }
101 } else { // printf(" even distribution\n");
102 double c=(k2max-fabs(k2min))/(double)(bfNr-1-zero);
103 for(int bi=zero; bi<bf->tacNr; bi++) {
104 bf->c[bi].size=fabs(k2min)+(double)(bi-zero)*c;
105 }
106 }
107
108
109#if(0) // The replacement works better with image-derived input with long frames
110
111 /* Allocate memory for simulated TAC */
112 double *sim;
113 sim=(double*)malloc(input->sampleNr*sizeof(double));
114 if(sim==NULL) {
115 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
116 tacFree(bf); return TPCERROR_OUT_OF_MEMORY;
117 }
118
119 /* Simulate the basis functions at input time points, and interpolate to tissue sample times */
120 double *cbi=NULL;
121 if(input->isframe!=0) {
122 // Time frames are available:
123 // Integrate BTAC to frame middle times and use it as input */
124 if(verbose>1) printf("using BTAC integral as simulation input\n");
125 cbi=(double*)malloc(input->sampleNr*sizeof(double));
126 if(cbi==NULL) {
127 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
128 tacFree(bf); free(sim); return TPCERROR_OUT_OF_MEMORY;
129 }
130 ret=liIntegratePET(input->x1, input->x2, input->c[0].y, input->sampleNr, cbi, NULL, verbose-10);
131 if(ret) {
132 if(verbose>1) printf("liIntegratePET() = %d\n", ret);
133 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
134 tacFree(bf); free(sim); free(cbi); return TPCERROR_INVALID_VALUE;
135 }
136 }
137 if(verbose>1) printf("computing basis functions\n");
138 ret=0;
139 for(int bi=0; bi<bf->tacNr && !ret; bi++) {
140 if(input->isframe==0) { // time frames not available
141 ret=simC1(input->x, input->c[0].y, input->sampleNr, 1.0, bf->c[bi].size, sim);
142 } else { // time frames are available
143 ret=simC1_i(input->x, cbi, input->sampleNr, 1.0, bf->c[bi].size, sim);
144 }
145 if(ret) break;
146 if(verbose>100) {
147 printf("\nk2 := %g\n", bf->c[bi].size);
148 printf("simulated TAC:\n");
149 for(int i=0; i<input->sampleNr; i++) printf(" %12.6f %12.3f\n", input->x[i], sim[i]);
150 }
151 /* interpolate to PET time frames */
152 if(tissue->isframe)
153 ret=liInterpolateForPET(input->x, sim, input->sampleNr,
154 tissue->x1, tissue->x2, bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
155 else
156 ret=liInterpolate(input->x, sim, input->sampleNr, tissue->x,
157 bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
158 if(ret) break;
159 } // next basis function
160 free(sim); free(cbi);
161 if(ret) {
162 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
163 return TPCERROR_INVALID_VALUE;
164 }
165
166#else
167
168 /* Interpolate the input function to times based on both input and tissue data */
169 /* When input TAC has long frames, the last simulated tissue frame will be biased without this. */
170 TAC inp; tacInit(&inp); TPCSTATUS status2; statusInit(&status2);
171 tacInput2sim(input, tissue, &inp, &status2);
172 statusFree(&status2);
173
174 /* Allocate memory for simulated TAC */
175 double *sim;
176 sim=(double*)malloc(inp.sampleNr*sizeof(double));
177 if(sim==NULL) {
178 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
179 tacFree(bf); return TPCERROR_OUT_OF_MEMORY;
180 }
181
182 if(verbose>1) printf("computing basis functions\n");
183 ret=0;
184 for(int bi=0; bi<bf->tacNr && !ret; bi++) {
185 ret=simC1(inp.x, inp.c[0].y, inp.sampleNr, 1.0, bf->c[bi].size, sim);
186 if(ret) break;
187 if(verbose>100) {
188 printf("\nk2 := %g\n", bf->c[bi].size);
189 printf("simulated TAC:\n");
190 for(int i=0; i<inp.sampleNr; i++) printf(" %12.6f %12.3f\n", inp.x[i], sim[i]);
191 }
192 /* interpolate to PET time frames */
193 if(tissue->isframe)
194 ret=liInterpolateForPET(inp.x, sim, inp.sampleNr,
195 tissue->x1, tissue->x2, bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
196 else
197 ret=liInterpolate(inp.x, sim, inp.sampleNr, tissue->x,
198 bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
199 if(ret) break;
200 } // next basis function
201 free(sim);
202 if(ret) {
203 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
204 return TPCERROR_INVALID_VALUE;
205 }
206
207 tacFree(&inp);
208#endif
209
210
211 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
212 return(TPCERROR_OK);
213}
liIntegratePET
int liIntegratePET(double *x1, double *x2, double *y, int nr, double *ie, double *iie, const int verbose)
Calculate PET TAC AUC from start to each time frame, as averages during each frame.
Definition integrate.c:120
liInterpolate
int liInterpolate(double *x, double *y, const int nr, double *newx, double *newy, double *newyi, double *newyii, const int newnr, const int se, const int ee, const int verbose)
Linear interpolation and/or integration with trapezoidal method.
Definition interpolate.c:139
liInterpolateForPET
int liInterpolateForPET(double *x, double *y, const int nr, double *newx1, double *newx2, double *newy, double *newyi, double *newyii, const int newnr, const int se, const int ee, const int verbose)
Linear TAC interpolation and/or integration to PET frames.
Definition interpolate.c:317
tacInput2sim
int tacInput2sim(TAC *itac, TAC *ttac, TAC *stac, TPCSTATUS *status)
Modify input TAC based on tissue TAC, for use in simulations.
Definition lisim.c:29
simC1
int simC1(double *t, double *ca, const int nr, const double k1, const double k2, double *ct)
Definition sim1cm.c:93
simC1_i
int simC1_i(double *t, double *cai, const int nr, const double k1, const double k2, double *ct)
Definition sim1cm.c:164
statusInit
void statusInit(TPCSTATUS *s)
Definition statusmsg.c:104
statusSet
void statusSet(TPCSTATUS *s, const char *func, const char *srcfile, int srcline, tpcerror error)
Definition statusmsg.c:142
statusFree
void statusFree(TPCSTATUS *s)
Definition statusmsg.c:126
TACC::name
char name[MAX_TACNAME_LEN+1]
Definition tpctac.h:81
TACC::y
double * y
Definition tpctac.h:75
TACC::size
double size
Definition tpctac.h:71
TAC
Definition tpctac.h:87
TAC::x
double * x
Definition tpctac.h:97
TAC::cunit
unit cunit
Definition tpctac.h:105
TAC::format
tacformat format
Definition tpctac.h:93
TAC::sampleNr
int sampleNr
Definition tpctac.h:89
TAC::isframe
int isframe
Definition tpctac.h:95
TAC::c
TACC * c
Definition tpctac.h:117
TAC::x2
double * x2
Definition tpctac.h:101
TAC::tunit
unit tunit
Definition tpctac.h:109
TAC::x1
double * x1
Definition tpctac.h:99
TAC::tacNr
int tacNr
Definition tpctac.h:91
TPCSTATUS::verbose
int verbose
Verbose level, used by statusPrint() etc.
Definition tpcextensions.h:242
TPCSTATUS::forgiving
int forgiving
Force level, 0 for strict tests for data units etc.
Definition tpcextensions.h:244
tacFree
void tacFree(TAC *tac)
Definition tac.c:106
tacAllocate
int tacAllocate(TAC *tac, int sampleNr, int tacNr)
Definition tac.c:130
tacInit
void tacInit(TAC *tac)
Definition tac.c:24
tacXCopy
int tacXCopy(TAC *tac1, TAC *tac2, int i1, int i2)
Definition tacx.c:24
TPCERROR_INVALID_VALUE
@ TPCERROR_INVALID_VALUE
Invalid value.
Definition tpcextensions.h:220
TPCERROR_UNKNOWN_UNIT
@ TPCERROR_UNKNOWN_UNIT
Unknown data unit.
Definition tpcextensions.h:226
TPCERROR_FAIL
@ TPCERROR_FAIL
General error.
Definition tpcextensions.h:192
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
unitIsTime
int unitIsTime(int u)
Definition units.c:359
TAC_FORMAT_TSV_UK
@ TAC_FORMAT_TSV_UK
UK TSV (point as decimal separator).
Definition tpctac.h:37

◆ bfmSRTM()

int bfmSRTM ( double * t,
double * cri,
const int n,
const int bfNr,
const double t3min,
const double t3max,
TAC * bf,
TPCSTATUS * status )

Calculate set of basis functions for SRTM.

Todo
Add tests.
Returns
enum tpcerror (TPCERROR_OK when successful).
See also
bfm1TCM
Parameters
tPointer to array containing PET sample (frame middle) times.
criPointer to array containing integral of reference tissue input concentration values at each sample time, AUC 0-t. AUC must have been calculated using values that are NOT corrected for decay.
nNr of samples (array lengths).
bfNrNr of basis functions to calculate.
t3minMinimum of theta3.
t3maxMaximum of theta3.
bfPointer to output TAC, to be allocated and filled with basis functions in here.
statusPointer to status data; enter NULL if not needed

Definition at line 26 of file bf_srtm.c.

45 {
46 int verbose=0; if(status!=NULL) verbose=status->verbose;
47 if(verbose>0) printf("%s(t, cri, %d, %d, %g, %g, bf)\n", __func__, n, bfNr, t3min, t3max);
48 if(t==NULL || cri==NULL || bf==NULL) {
49 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
50 return TPCERROR_FAIL;
51 }
52 if(n<2 || bfNr<1) {
53 if(verbose>1) printf("invalid sample or BF number\n");
54 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
55 return TPCERROR_TOO_FEW;
56 }
57 if(t3min<1.0E-10 || t3min>=t3max) { // range cannot be computed otherwise
58 if(verbose>1) printf("invalid theta3 range\n");
59 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
60 return TPCERROR_INVALID_VALUE;
61 }
62
63 int ret;
64
65 /* Allocate memory for basis functions */
66 if(verbose>1) printf("allocating memory for basis functions\n");
67 tacFree(bf);
68 ret=tacAllocate(bf, n, bfNr);
69 if(ret!=TPCERROR_OK) {
70 statusSet(status, __func__, __FILE__, __LINE__, ret);
71 return ret;
72 }
73 /* Copy and set information fields */
74 bf->tacNr=bfNr; bf->sampleNr=n;
75 bf->format=TAC_FORMAT_TSV_UK;
76 bf->isframe=0;
77 for(int bi=0; bi<bf->tacNr; bi++) sprintf(bf->c[bi].name, "B%5.5d", bi+1);
78 for(int fi=0; fi<bf->sampleNr; fi++) bf->x[fi]=t[fi];
79 /* Compute theta3 values to size fields */
80 {
81 if(verbose>1) printf("computing theta3 values\n");
82 double a, b, c;
83 a=log10(t3min); b=log10(t3max); c=(b-a)/(double)(bfNr-1);
84 for(int bi=0; bi<bf->tacNr; bi++) {
85 bf->c[bi].size=pow(10.0, (double)bi*c+a);
86 }
87 }
88 if(verbose>2) {
89 printf("bf_t3_range := %g - %g\n", bf->c[0].size, bf->c[bf->tacNr-1].size);
90 printf("bf_t3_gaps := %g - %g\n", bf->c[1].size-bf->c[0].size,
91 bf->c[bf->tacNr-1].size-bf->c[bf->tacNr-2].size);
92 }
93
94 /* Calculate the basis functions */
95 if(verbose>1) printf("computing basis functions\n");
96 ret=0;
97 for(int bi=0; bi<bf->tacNr && !ret; bi++) {
98 if(verbose>99) printf(" theta3=%g\n", bf->c[bi].size);
99 ret=simC1_i(t, cri, n, 1.0, bf->c[bi].size, bf->c[bi].y);
100 }
101 if(ret) {
102 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
103 return TPCERROR_INVALID_VALUE;
104 }
105
106 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
107 return(TPCERROR_OK);
108}

◆ 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
kPointer to k array of length n; must be sorted to either order, as those usually are. Contents are not modified.
aPointer to a array of length n; not modified.
nLength of k and a arrays.
kePointer to an allocated array of extracted k values of length ne.
aePointer to an allocated array of extracted a values of length ne; not modified.
neLength of k and a arrays, and the maximum number of parameters to extract.

Definition at line 227 of file bf_dexp.c.

241 {
242 if(k==NULL || a==NULL || n<1 || ke==NULL || ae==NULL || ne<1) return(0);
243
244 /* First, extract all zero-separated basis function clusters */
245 double *buf=NULL, *ok, *oa;
246 int bfn=spectralBFNr(k, a, n); if(bfn==0) return(0);
247 if(bfn>ne) {
248 /* we need to allocate memory locally */
249 buf=(double*)malloc(2*bfn*sizeof(double)); if(buf==NULL) return(0);
250 ok=buf; oa=buf+bfn;
251 } else {
252 /* we just set pointers to user-provided arrays */
253 ok=ke; oa=ae;
254 }
255 {
256 int nn=0, i=0, ii=0;
257 do {
258 i=doubleCSpanPositives(a+ii, n-ii);
259 ii+=i; if(n-ii<1) break;
260 i=doubleSpanPositives(a+ii, n-ii);
261 if(i==1) { // cluster contains just one positive a
262 oa[nn]=a[ii]; ok[nn]=k[ii];
263 nn++;
264 } else if(i>1) { // cluster contains more than one positive a
265 /* get sum of a values in the cluster */
266 oa[nn]=doubleSum(a+ii, i);
267 /* get weighted mean of k values in the cluster */
268 ok[nn]=doubleWMean(k+ii, a+ii, i);
269 nn++;
270 }
271 ii+=i;
272 } while(ii<n && nn<bfn);
273 if(0) {
274 printf("\n\tCluster a and k values\n");
275 for(int i=0; i<bfn; i++) printf("\t%g\t%g\n", oa[i], ok[i]);
276 printf("\n"); fflush(stdout);
277 }
278 }
279 /* If the number needs not to be reduced, then we are ready */
280 if(ne>=bfn) {
281 for(int i=bfn; i<ne; i++) ke[i]=ae[i]=0.0;
282 if(buf!=NULL) free(buf);
283 return(bfn);
284 }
285
286 while(ne<bfn) {
287 /* Search the smallest a parameter */
288 int ami=doubleAbsMinIndex(oa, bfn);
289 if(0) printf("smallest |a[%d]|=%g\n", ami, oa[ami]);
290 /* Combine it with previous or next cluster? */
291 int ci;
292 if(ami==0) ci=1;
293 else if(ami==bfn-1) ci=bfn-2;
294 else {if(fabs(oa[ami-1])<fabs(oa[ami+1])) ci=ami-1; else ci=ami+1;}
295 if(0) printf("combining clusters %d and %d\n", 1+ami, 1+ci);
296 /* Replace the cluster with combination of two clusters */
297 ok[ci]=(oa[ami]*ok[ami]+oa[ci]*ok[ci])/(oa[ami]+oa[ci]);
298 oa[ci]+=oa[ami];
299 /* Delete the cluster */
300 for(int i=ami+1; i<bfn; i++) oa[i-1]=oa[i];
301 for(int i=ami+1; i<bfn; i++) ok[i-1]=ok[i];
302 bfn--;
303 if(0) {
304 printf("\n\tAfter reduction, cluster a and k values\n");
305 for(int i=0; i<bfn; i++) printf("\t%g\t%g\n", oa[i], ok[i]);
306 printf("\n"); fflush(stdout);
307 }
308 }
309 for(int i=0; i<ne; i++) ae[i]=oa[i];
310 for(int i=0; i<ne; i++) ke[i]=ok[i];
311 if(buf!=NULL) free(buf);
312
313 return(bfn);
314}
spectralBFNr
int spectralBFNr(double *k, double *a, const int n)
Definition bf_dexp.c:193
doubleSpanPositives
int doubleSpanPositives(double *a, const int n)
Definition doubleutil.c:302
doubleAbsMinIndex
unsigned int doubleAbsMinIndex(double *a, const unsigned int n)
Definition doubleutil.c:416
doubleCSpanPositives
int doubleCSpanPositives(double *a, const int n)
Definition doubleutil.c:320
doubleWMean
double doubleWMean(double *a, double *w, const unsigned int n)
Definition doubleutil.c:244
doubleSum
double doubleSum(double *a, const unsigned int n)
Definition doubleutil.c:206

◆ 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
kPointer to k array of length n; does not need to be in specific order, although those usually are. Contents are not modified.
aPointer to a array of length n; not modified.
nLength of k and a arrays.

Definition at line 193 of file bf_dexp.c.

201 {
202 if(k==NULL || a==NULL || n<1) return(0);
203 int bfn=0, i=0, ii=0;
204 do {
205 i=doubleCSpanPositives(a+ii, n-ii);
206 ii+=i; if(n-ii<1) break;
207 i=doubleSpanPositives(a+ii, n-ii);
208 ii+=i; if(i>0) bfn++;
209 } while(ii<n);
210 return(bfn);
211}

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
xPointer to TAC x data (not modified). Data must be sorted by increasing x. Negative or missing x values are not allowed.
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[] and y[]; at least 3.
kminMinimum of k (set based on sample range and unit); must be >0.
kmaxMaximum of k (set based on sample range and unit); must be >kmin.
bnrNumber of basis functions to calculate; also the length of arrays k[] and a[]; at least 4.
kPointer to array of length bnr for k values, filled in by this function.
aPointer to array of length bnr for a values, filled in by this function; notice that most values may be set to zero.
yfitPointer to array for fitted y values. Enter NULL if not needed.
statusPointer to status data; enter NULL if not needed

Definition at line 28 of file bf_dexp.c.

55 {
56 int verbose=0; if(status!=NULL) verbose=status->verbose;
57 if(verbose>0) printf("%s()\n", __func__);
58 if(x==NULL || y==NULL || k==NULL || a==NULL || snr<1 || bnr<1) {
59 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
60 return TPCERROR_NO_DATA;
61 }
62 if(snr<3 || bnr<4) {
63 if(verbose>1) fprintf(stderr, "invalid sample or BF number\n");
64 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
65 return TPCERROR_TOO_FEW;
66 }
67 if(!(kmin>0.0) || kmin>=kmax) { // range cannot be computed otherwise
68 if(verbose>1) fprintf(stderr, "invalid k range\n");
69 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
70 return TPCERROR_INVALID_VALUE;
71 }
72 int ret;
73
74 if(verbose>1) printf("computing k values\n");
75 {
76 double a, b, c;
77 a=log10(kmin); b=log10(kmax); c=(b-a)/(double)(bnr-1);
78 if(verbose>3) printf(" a := %g\n b := %g\n c := %g\n", a, b, c);
79 for(int bi=0; bi<bnr; bi++) k[bi]=pow(10.0, (double)bi*c+a);
80 }
81
82 /* Allocate memory required by NNLS */
83 int nnls_n, nnls_m, n, m, nnls_index[bnr];
84 double *nnls_a[bnr], *nnls_b, *nnls_zz, nnls_x[bnr], *nnls_mat,
85 nnls_wp[bnr], *dptr, nnls_rnorm;
86 if(verbose>1) printf("allocating memory for NNLS\n");
87 nnls_n=bnr; nnls_m=snr;
88 nnls_mat=(double*)malloc(((nnls_n+2)*nnls_m)*sizeof(double));
89 if(nnls_mat==NULL) {
90 if(verbose>1) fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
91 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
92 return TPCERROR_OUT_OF_MEMORY;
93 }
94 for(n=0, dptr=nnls_mat; n<nnls_n; n++) {nnls_a[n]=dptr; dptr+=nnls_m;}
95 nnls_b=dptr; dptr+=nnls_m; nnls_zz=dptr;
96
97 if(verbose>1) printf("filling NNLS data matrix\n");
98 /* Fill NNLS B array with measured y values */
99 for(m=0; m<nnls_m; m++) nnls_b[m]=y[m];
100
101 /* Fill NNLS A matrix with basis functions, calculated in place */
102 for(n=0; n<nnls_n; n++)
103 for(m=0; m<nnls_m; m++)
104 nnls_a[n][m]=exp(-k[n]*x[m]);
105
106 /* Apply weights, if given */
107 if(w!=NULL) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, w);
108
109 /* NNLS */
110 if(verbose>1) printf("applying NNLS\n");
111 ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm, nnls_wp, nnls_zz, nnls_index);
112 if(ret>1) {
113 if(verbose>1) fprintf(stderr, "NNLS solution not possible.\n");
114 free(nnls_mat);
115 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_SOLUTION);
116 return TPCERROR_NO_SOLUTION;
117 }
118 /* Copy a[] values */
119 if(a!=NULL) for(n=0; n<nnls_n; n++) a[n]=nnls_x[n];
120
121 /* Compute fitted y[] */
122 if(yfit!=NULL) {
123 if(verbose>1) printf("computing yfit[]\n");
124 for(m=0; m<nnls_m; m++) {
125 yfit[m]=0.0;
126 for(n=0; n<nnls_n; n++)
127 if(nnls_x[n]>0.0)
128 yfit[m]+=nnls_x[n]*exp(-k[n]*x[m]);
129 }
130 }
131 free(nnls_mat);
132
133 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
134 return TPCERROR_OK;
135}
nnls
int nnls(double **a, int m, int n, double *b, double *x, double *rnorm, double *wp, double *zzp, int *indexp)
Definition nnls.c:48
nnlsWght
int nnlsWght(int N, int M, double **A, double *b, double *weight)
Definition nnls.c:264
TPCERROR_NO_SOLUTION
@ TPCERROR_NO_SOLUTION
No solution.
Definition tpcextensions.h:227

◆ 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

◆ 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
kPointer to k array of length n; does not need to be in specific order, although those usually are. Contents are not modified.
aPointer to a array of length n; not modified.
nLength of k and a arrays.
kminPointer min k value. Enter NULL if not needed.
kmaxPointer max k value. Enter NULL if not needed.
statusPointer to status data; enter NULL if not needed

Definition at line 145 of file bf_dexp.c.

159 {
160 int verbose=0; if(status!=NULL) verbose=status->verbose;
161 if(verbose>0) printf("%s(k[], a[], %d, *kmin, *kmax)\n", __func__, n);
162 if(k==NULL || a==NULL || n<1) {
163 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
164 return TPCERROR_NO_DATA;
165 }
166 double max=nan(""), min=nan("");
167 for(int i=0; i<n; i++) if(a[i]>0.0 && isfinite(k[i])) {
168 if(isnan(max) || k[i]>max) max=k[i];
169 if(isnan(min) || k[i]<min) min=k[i];
170 }
171 if(verbose>2) {
172 printf(" kmin := %g\n", min);
173 printf(" kmax := %g\n", max);
174 }
175 if(kmin!=NULL) *kmin=min;
176 if(kmax!=NULL) *kmax=max;
177 if(isnan(min) || isnan(max)) {
178 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_VALUE);
179 return TPCERROR_NO_VALUE;
180 }
181 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
182 return TPCERROR_OK;
183}
TPCERROR_NO_VALUE
@ TPCERROR_NO_VALUE
Value not found.
Definition tpcextensions.h:219
Generated on for TPCCLIB by doxygen 1.15.0