TPCCLIB
Loading...
Searching...
No Matches
func.c
Go to the documentation of this file.
1
5/*****************************************************************************/
6#include "tpcclibConfig.h"
7/*****************************************************************************/
8#include <stdio.h>
9#include <stdlib.h>
10#include <math.h>
11#include <time.h>
12#include <string.h>
13/*****************************************************************************/
14#include "tpcextensions.h"
15/*****************************************************************************/
16#include "tpcfunc.h"
17/*****************************************************************************/
18
19/*****************************************************************************/
26int mfEvalY(
28 const char *fid,
30 const int parNr,
32 double *p,
34 const int sampleNr,
36 double *x,
38 double *y,
40 const int verbose
41) {
42 if(verbose>0) {
43 printf("%s('%s', %d, p[], %d, x[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
44 fflush(stdout);
45 }
46 if(parNr<1 || p==NULL || sampleNr<1 || x==NULL || y==NULL) return(1);
47
48 if(strcasecmp(fid, "step")==0) {
49 /* Parameter number must be an even number and >=2 */
50 if(parNr<2 || parNr&1) return(2);
51 for(int i=0; i<sampleNr; i++) {
52 y[i]=0.0;
53 for(int j=0; j<parNr; j+=2) if(x[i]>=p[j]) y[i]=p[j+1];
54 }
55 return(0);
56 }
57
58 if(strcasecmp(fid, "level")==0) {
59 for(int i=0; i<sampleNr; i++) y[i]=p[0];
60 return(0);
61 }
62
63 if(strcasecmp(fid, "fengm2")==0) {
64 if(parNr<6) return(2);
65 double delayt=0.0; if(parNr>6) delayt=p[6];
66 for(int i=0; i<sampleNr; i++) {
67 double xt=x[i]-delayt; y[i]=0.0; if(!(xt>0.0)) continue;
68 y[i] += (p[0]*xt-p[2]-p[4])*exp(p[1]*xt);
69 y[i] += p[2]*exp(p[3]*xt);
70 y[i] += p[4]*exp(p[5]*xt);
71 }
72 return(0);
73 }
74 if(strcasecmp(fid, "fengm2s")==0) {
75 if(parNr<4) return(2);
76 double delayt=0.0; if(parNr>4) delayt=p[4];
77 for(int i=0; i<sampleNr; i++) {
78 double xt=x[i]-delayt; y[i]=0.0; if(!(xt>0.0)) continue;
79 y[i] += (p[0]*xt-p[2])*exp(p[1]*xt);
80 y[i] += p[2]*exp(p[3]*xt);
81 }
82 return(0);
83 }
84 if(strcasecmp(fid, "fengm2e")==0) {
85 if(parNr<8) return(2);
86 double delayt=0.0; if(parNr>8) delayt=p[8];
87 for(int i=0; i<sampleNr; i++) {
88 double xt=x[i]-delayt; y[i]=0.0; if(!(xt>0.0)) continue;
89 y[i] += (p[0]*xt-p[2]-p[4]-p[6])*exp(p[1]*xt);
90 y[i] += p[2]*exp(p[3]*xt);
91 y[i] += p[4]*exp(p[5]*xt);
92 y[i] += p[6]*exp(p[7]*xt);
93 }
94 return(0);
95 }
96 /* Gamma variate function */
97 if(strcasecmp(fid, "gammav")==0) {
98 if(parNr<3) return(2);
99 double delayt=0.0; if(parNr>3) delayt=p[3];
100 for(int i=0; i<sampleNr; i++) {
101 double xt=x[i]-delayt; y[i]=0.0;
102 if(!(xt>0.0) || fabs(p[2])<1.0E-10) continue;
103 y[i] = p[0]*pow(xt, p[1])*exp(-xt/p[2]);
104 }
105 return(0);
106 }
107 /* Exponential function */
108 if(strcasecmp(fid, "1exp")==0) {
109 if(parNr<2) return(2);
110 for(int i=0; i<sampleNr; i++) {
111 y[i]=0.0;
112 y[i] += p[0]*exp(p[1]*x[i]);
113 }
114 return(0);
115 }
116 /* Sum of two exponentials */
117 if(strcasecmp(fid, "2exp")==0) {
118 if(parNr<4) return(2);
119 for(int i=0; i<sampleNr; i++) {
120 y[i]=0.0;
121 y[i] += p[0]*exp(p[1]*x[i]);
122 y[i] += p[2]*exp(p[3]*x[i]);
123 }
124 return(0);
125 }
126 /* Sum of three exponentials */
127 if(strcasecmp(fid, "3exp")==0) {
128 if(parNr<6) return(2);
129 for(int i=0; i<sampleNr; i++) {
130 y[i]=0.0;
131 y[i] += p[0]*exp(p[1]*x[i]);
132 y[i] += p[2]*exp(p[3]*x[i]);
133 y[i] += p[4]*exp(p[5]*x[i]);
134 }
135 return(0);
136 }
137 /* Sum of four exponentials */
138 if(strcasecmp(fid, "4exp")==0) {
139 if(parNr<8) return(2);
140 for(int i=0; i<sampleNr; i++) {
141 y[i]=0.0;
142 y[i] += p[0]*exp(p[1]*x[i]);
143 y[i] += p[2]*exp(p[3]*x[i]);
144 y[i] += p[4]*exp(p[5]*x[i]);
145 y[i] += p[6]*exp(p[7]*x[i]);
146 }
147 return(0);
148 }
149 /* Sum of five exponentials */
150 if(strcasecmp(fid, "5exp")==0) {
151 if(parNr<10) return(2);
152 for(int i=0; i<sampleNr; i++) {
153 y[i]=0.0;
154 y[i] += p[0]*exp(p[1]*x[i]);
155 y[i] += p[2]*exp(p[3]*x[i]);
156 y[i] += p[4]*exp(p[5]*x[i]);
157 y[i] += p[6]*exp(p[7]*x[i]);
158 y[i] += p[8]*exp(p[9]*x[i]);
159 }
160 return(0);
161 }
162
163 /* Inverted gamma cdf for plasma parent fraction */
164 if(strcasecmp(fid, "ppfigamc")==0) {
165 if(parNr<4) return(2);
166 double delayt=0.0; if(parNr>4) delayt=p[4];
167 double a=p[0];
168 double b=p[1];
169 double c=p[2]; if(c<0.0) return(2);
170 double d=p[3]; if(d<=0.0) return(2);
171 for(int i=0; i<sampleNr; i++) {
172 double t=x[i]-delayt; if(t<0.0) t=0.0;
173 y[i]=a*(1.0-b*igam(d, c*t));
174 }
175 return(0);
176 }
177
178 /* Weibull cdf with time delay */
179 if(strcasecmp(fid, "weibullcdfd")==0) {
180 if(parNr<4) return(2);
181 double delayt=p[0];
182 for(int i=0; i<sampleNr; i++) {
183 double xt=x[i]-delayt; y[i]=0.0;
184 if(!(xt>0.0)) continue;
185 y[i] = p[1]*(1.0-exp(-pow(xt/p[2], p[3])));
186 }
187 return(0);
188 }
189
190 /* Weibull cdf with time delay and cdf derivative */
191 if(strcasecmp(fid, "weibullcdfdd")==0) {
192 if(parNr<5) return(2);
193 double delayt=p[0];
194 for(int i=0; i<sampleNr; i++) {
195 double xt=x[i]-delayt; y[i]=0.0;
196 if(!(xt>0.0)) continue;
197 double a=exp(-pow(xt/p[2], p[3]));
198 double b=(p[3]/p[2]) * pow(xt/p[2], p[3]-1.0) * a;
199 y[i] = p[1]*(b + p[4]*(1.0-a));
200 }
201 return(0);
202 }
203
204 /* Surge function in form f(x)=a*x*exp(-b*x)*b^2 to give AUC=a from 0 to infinity.
205 Maximum value is at 1/b. */
206 if(strcasecmp(fid, "surge")==0) {
207 if(parNr<2) return(2);
208 double f=p[0]*(p[1]*p[1]);
209 for(int i=0; i<sampleNr; i++) {
210 double xt=x[i]; y[i]=0.0;
211 if(!(xt>0.0)) continue;
212 y[i] = f*xt*exp(-p[1]*xt);
213 }
214 return(0);
215 }
216
217 /* Surge function in form f(x)=a*x*exp(-b*x). Maximum value is at 1/b. */
218 if(strcasecmp(fid, "tradsurge")==0) {
219 if(parNr<2) return(2);
220 for(int i=0; i<sampleNr; i++) {
221 double xt=x[i]; y[i]=0.0;
222 if(!(xt>0.0)) continue;
223 y[i] = p[0]*xt*exp(-p[1]*xt);
224 }
225 return(0);
226 }
227
228 /* Surge function with recirculation, in form
229 f(x)=a*(x*exp(-b*x) + (c/(b*b))*(1-(b*x+1)*exp(-b*x))). */
230 if(strcasecmp(fid, "surgerecirc")==0) {
231 if(parNr<2) return(2);
232 double c=0.0;
233 if(parNr>=3) c=p[2];
234 for(int i=0; i<sampleNr; i++) {
235 double xt=x[i]; y[i]=0.0;
236 if(!(xt>0.0)) continue;
237 double e=exp(-p[1]*xt);
238 y[i] = xt*e + (c/(p[1]*p[1]))*(1.0 - (p[1]*xt+1.0)*e);
239 y[i] *= p[0];
240 }
241 return(0);
242 }
243
244 /* Surge function with recirculation and delay, in form
245 f(x)=a*((x-d)*exp(-b*(x-d)) + (c/(b*b))*(1-(b*(x-d)+1)*exp(-b*(x-d)))). */
246 if(strcasecmp(fid, "surgerecircd")==0) {
247 if(parNr<2) return(2);
248 double c=0.0; if(parNr>=3) c=p[2];
249 double delayt=0.0; if(parNr>=4) delayt=p[3];
250 for(int i=0; i<sampleNr; i++) {
251 double xt=x[i]-delayt; y[i]=0.0;
252 if(!(xt>0.0)) continue;
253 double e=exp(-p[1]*xt);
254 y[i] = xt*e + (c/(p[1]*p[1]))*(1.0 - (p[1]*xt+1.0)*e);
255 y[i] *= p[0];
256 }
257 return(0);
258 }
259
260 /* Surge function with recirculation, for plasma-to-blood ratio.
261 RBC-to-plasma r(x)=a*(x*exp(-b*x) + (c/(b*b))*(1-(b*x+1)*exp(-b*x)))
262 Plasma-to-blood f(x)=1/(1-HCT*(1-r(x))). */
263 if(strcasecmp(fid, "p2bsrc")==0) {
264 if(parNr<3) return(2);
265 double c=0.0; if(parNr>=4) c=p[3];
266 double HCT=p[0];
267 double a=p[1];
268 double b=p[2];
269 for(int i=0; i<sampleNr; i++) {
270 double xt=x[i]; if(!(xt>0.0)) {y[i]=1.0/(1.0-HCT); continue;}
271 double e=exp(-b*xt);
272 double r=a*(xt*e + (c/(b*b))*(1.0 - (b*xt+1.0)*e));
273 y[i] = 1.0/(1.0-HCT*(1.0-r));
274 }
275 return(0);
276 }
277
278 /* Feng M2 function for plasma-to-blood ratio.
279 RBC-to-plasma ratio r(x) using Feng M2, and then Plasma-to-blood f(x)=1/(1-HCT*(1-r(x))). */
280 if(strcasecmp(fid, "p2bfm2")==0) {
281 if(parNr<3) return(2);
282 double HCT=p[0];
283 double A1=p[1];
284 double L1=p[2];
285 double A2=0.0; if(parNr>3) A2=p[3];
286 double L2=0.0; if(parNr>4) L2=p[4];
287 double A3=0.0; if(parNr>5) A3=p[5];
288 double L3=0.0; if(parNr>6) L3=p[6];
289 for(int i=0; i<sampleNr; i++) {
290 double xt=x[i]; if(!(xt>0.0)) {y[i]=1.0/(1.0-HCT); continue;}
291 double r=(A1*xt - A2 - A3)*exp(-L1*xt) + A2*exp(-L2*xt) + A3*exp(-L3*xt);
292 y[i] = 1.0/(1.0-HCT*(1.0-r));
293 }
294 return(0);
295 }
296
297 /* Surge function for late FDG AIF with time delay */
298 if(strcasecmp(fid, "surgefdgaif")==0) {
299 if(parNr<5) return(2);
300 double delayt=p[0];
301 for(int i=0; i<sampleNr; i++) {
302 double xt=x[i]-delayt; y[i]=0.0;
303 if(!(xt>0.0)) continue;
304 y[i] = p[2]*(p[3]*xt*exp(-p[4]*p[1]*xt) + exp(-p[1]*xt));
305 }
306 return(0);
307 }
308
309 /* Probability density function of Erlang distribution */
310 if(strcasecmp(fid, "erlangpdf")==0) {
311 if(parNr<3) return(2);
312 double A=p[0], k=p[1]; if(!(k>=0.0)) return(2);
313 int N=(int)round(p[2]); if(!(N>0)) return(2);
314 unsigned long long int f=lfactorial(N-1);
315 for(int i=0; i<sampleNr; i++) {
316 y[i]=0.0; if(!(x[i]>=0.0)) continue;
317 y[i] = A*pow(k, N)*pow(x[i], N-1)*exp(-k*x[i])/f;
318 }
319 return(0);
320 }
321
322 /* Rational functions */
323 if(strcasecmp(fid, "ratf11")==0) {
324 if(parNr<4) return(2);
325 for(int i=0; i<sampleNr; i++) {
326 double a=p[0]+p[2]*x[i];
327 double b=p[1]+p[3]*x[i];
328 y[i]=a/b;
329 }
330 return(0);
331 }
332 if(strcasecmp(fid, "ratf21")==0) {
333 if(parNr<5) return(2);
334 for(int i=0; i<sampleNr; i++) {
335 double a=p[0]+p[2]*x[i]+p[4]*x[i]*x[i];
336 double b=p[1]+p[3]*x[i];
337 y[i]=a/b;
338 }
339 return(0);
340 }
341 if(strcasecmp(fid, "ratf22")==0) {
342 if(parNr<6) return(2);
343 for(int i=0; i<sampleNr; i++) {
344 double sqrx=x[i]*x[i];
345 double a=p[0]+p[2]*x[i]+p[4]*sqrx;
346 double b=p[1]+p[3]*x[i]+p[5]*sqrx;
347 y[i]=a/b;
348 }
349 return(0);
350 }
351 if(strcasecmp(fid, "ratf32")==0) {
352 if(parNr<7) return(2);
353 for(int i=0; i<sampleNr; i++) {
354 double sqrx=x[i]*x[i];
355 double a=p[0]+p[2]*x[i]+p[4]*sqrx+p[6]*sqrx*x[i];
356 double b=p[1]+p[3]*x[i]+p[5]*sqrx;
357 y[i]=a/b;
358 }
359 return(0);
360 }
361 if(strcasecmp(fid, "ratf33")==0) {
362 if(parNr<8) return(2);
363 for(int i=0; i<sampleNr; i++) {
364 double sqrx=x[i]*x[i];
365 double cubx=sqrx*x[i];
366 double a=p[0]+p[2]*x[i]+p[4]*sqrx+p[6]*cubx;
367 double b=p[1]+p[3]*x[i]+p[5]*sqrx+p[7]*cubx;
368 y[i]=a/b;
369 }
370 return(0);
371 }
372 if(strcasecmp(fid, "p2brf")==0) {
373 if(parNr<7) return(2);
374 for(int i=0; i<sampleNr; i++) {
375 if(x[i]<=0.0) {
376 y[i]=1.0/(1.0-p[0]);
377 } else {
378 double sqrx=x[i]*x[i];
379 double cubx=sqrx*x[i];
380 double a=0.0+p[1]*x[i]+p[3]*sqrx+p[5]*cubx;
381 double b=1.0+p[2]*x[i]+p[4]*sqrx+p[6]*cubx;
382 y[i]=1.0/(1.0-p[0]*(1.0-(a/b)));
383 }
384 }
385 return(0);
386 }
387
388 /* Hill functions */
389 if(strcasecmp(fid, "mpfhill")==0) {
390 if(parNr<3) return(2);
391 for(int i=0; i<sampleNr; i++) {
392 y[i]=p[0]*pow(x[i], p[1]) / (pow(x[i], p[1]) + p[2]);
393 }
394 return(0);
395 }
396 if(strcasecmp(fid, "p2bhill")==0) {
397 if(parNr<4) return(2);
398 for(int i=0; i<sampleNr; i++) {
399 double cpr=p[1]*pow(x[i], p[2]) / (pow(x[i], p[2]) + p[3]);
400 y[i]=1.0/(1.0-p[0]*(1.0-cpr));
401 }
402 return(0);
403 }
404
405
406 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
407 return(10);
408}
409/*****************************************************************************/
410
411/*****************************************************************************/
418int mfEvalInt(
420 const char *fid,
422 const int parNr,
424 double *p,
426 const int sampleNr,
428 double *x,
430 double *v,
432 const int verbose
433) {
434 if(verbose>0) {
435 printf("%s('%s', %d, p[], %d, x[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
436 fflush(stdout);
437 }
438 if(parNr<1 || p==NULL || sampleNr<1 || x==NULL || v==NULL) return(1);
439
440 if(strcasecmp(fid, "fengm2")==0) {
441 if(parNr<6) return(2);
442 double delayt=0.0; if(parNr>6) delayt=p[6];
443 double a;
444 double A1, A2, A3, L1, L2, L3; // lambdas are negative
445 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5];
446 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20 || fabs(L3)<1.0E-20) {
447 if(verbose>1) printf("too small L parameter(s)\n");
448 return(3);
449 }
450 for(int i=0; i<sampleNr; i++) {
451 double xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
452 if(L1!=0.0) {
453 a=exp(L1*xt);
454 v[i] += (A1+L1*(A2+A3))*(1.0-a)/(L1*L1);
455 v[i] += (A1/L1)*xt*a;
456 }
457 if(L2!=0.0) v[i]-=(A2/L2)*(1.0-exp(L2*xt)); else v[i]+=A2*xt;
458 if(L3!=0.0) v[i]-=(A3/L3)*(1.0-exp(L3*xt)); else v[i]+=A3*xt;
459 }
460 return(0);
461 }
462 if(strcasecmp(fid, "fengm2e")==0) {
463 if(parNr<8) return(2);
464 double delayt=0.0; if(parNr>8) delayt=p[8];
465 double a;
466 double A1, A2, A3, A4, L1, L2, L3, L4; // lambdas are negative
467 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5]; A4=p[6]; L4=p[7];
468 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20 || fabs(L3)<1.0E-20 || fabs(L4)<1.0E-20) {
469 if(verbose>1) printf("too small L parameter(s)\n");
470 return(3);
471 }
472 for(int i=0; i<sampleNr; i++) {
473 double xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
474 if(L1!=0.0) {
475 a=exp(L1*xt);
476 v[i] += (A1+L1*(A2+A3+A4))*(1.0-a)/(L1*L1);
477 v[i] += (A1/L1)*xt*a;
478 }
479 if(L2!=0.0) v[i]-=(A2/L2)*(1.0-exp(L2*xt)); else v[i]+=A2*xt;
480 if(L3!=0.0) v[i]-=(A3/L3)*(1.0-exp(L3*xt)); else v[i]+=A3*xt;
481 if(L4!=0.0) v[i]-=(A4/L4)*(1.0-exp(L4*xt)); else v[i]+=A4*xt;
482 }
483 return(0);
484 }
485 if(strcasecmp(fid, "fengm2s")==0) {
486 if(parNr<4) return(2);
487 double delayt=0.0; if(parNr>4) delayt=p[4];
488 double a;
489 double A1, A2, L1, L2; // lambdas are negative
490 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3];
491 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20) {
492 if(verbose>1) printf("too small L parameter(s)\n");
493 return(3);
494 }
495 for(int i=0; i<sampleNr; i++) {
496 double xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
497 if(L1!=0.0) {
498 a=exp(L1*xt);
499 v[i] += (A1+L1*A2)*(1.0-a)/(L1*L1);
500 v[i] += (A1/L1)*xt*a;
501 }
502 if(L2!=0.0) v[i]-=(A2/L2)*(1.0-exp(L2*xt)); else v[i]+=A2*xt;
503 }
504 return(0);
505 }
506 /* Gamma variate function when p[1]==1 */
507 if(strcasecmp(fid, "gammav")==0 && parNr>=3 && fabs(p[1]-1.0)<1.0E-10) {
508 if(parNr<3) return(2);
509 double delayt=0.0; if(parNr>3) delayt=p[3];
510 for(int i=0; i<sampleNr; i++) {
511 double xt=x[i]-delayt; v[i]=0.0;
512 if(!(xt>0.0) || fabs(p[2])<1.0E-10) continue;
513 v[i] = p[0]*p[2]*(p[2] - (p[2]+xt)*exp(-xt/p[2]));
514 }
515 return(0);
516 }
517
518 /* Surge function in form f(x)=a*x*exp(-b*x)*b^2, to give AUC=a from 0 to infinity. */
519 if(strcasecmp(fid, "surge")==0) {
520 if(parNr<2) return(2);
521 for(int i=0; i<sampleNr; i++) {
522 double xt=x[i]; v[i]=0.0;
523 if(!(xt>0.0)) continue;
524 v[i] = p[0]*(1.0 - (p[1]*xt + 1.0)*exp(-p[1]*xt));
525 }
526 return(0);
527 }
528
529 /* Surge function in form f(x)=a*x*exp(-b*x). Maximum value is at 1/b. */
530 if(strcasecmp(fid, "tradsurge")==0) {
531 if(parNr<2) return(2);
532 double f=p[0]/(p[1]*p[1]);
533 for(int i=0; i<sampleNr; i++) {
534 double xt=x[i]; v[i]=0.0;
535 if(!(xt>0.0)) continue;
536 v[i] = f*(1.0 - (p[1]*xt + 1.0)*exp(-p[1]*xt));
537 }
538 return(0);
539 }
540
541 /* Surge function for late FDG AIF with time delay */
542 if(strcasecmp(fid, "surgefdgaif")==0) {
543 if(parNr<5) return(2);
544 double delayt=p[0];
545 for(int i=0; i<sampleNr; i++) {
546 double xt=x[i]-delayt; v[i]=0.0;
547 if(!(xt>0.0)) continue;
548 v[i] += p[3]*(1.0-(p[4]*p[1]*xt + 1.0)*exp(-p[4]*p[1]*xt))/(p[1]*p[1]*p[4]*p[4]);
549 v[i] += (1.0-exp(-p[1]*xt))/p[1];
550 v[i] *= p[2];
551 }
552 return(0);
553 }
554
555 /* Probability density function of Erlang distribution */
556 if(strcasecmp(fid, "erlangpdf")==0 && p[2]<20.5) { // Supports only N=0,1,2,...,20
557 if(parNr<3) return(2);
558 double A=p[0], k=p[1]; if(!(k>=0.0)) return(2);
559 int N=(int)round(p[2]); if(!(N>=0) || N>20) return(2);
560 if(N==0)
561 for(int i=0; i<sampleNr; i++) {
562 if(x[i]>=0.0) v[i]=A; else v[i]=0.0;
563 }
564 else if(N==1)
565 for(int i=0; i<sampleNr; i++) {
566 if(x[i]>=0.0) v[i]=A*(1.0-exp(-k*x[i])); else v[i]=0.0;
567 }
568 else if(N==2)
569 for(int i=0; i<sampleNr; i++) {
570 if(x[i]>=0.0) v[i]=A*(1.0 - exp(-k*x[i]) - k*x[i]*exp(-k*x[i])); else v[i]=0.0;
571 }
572 else {
573 for(int i=0; i<sampleNr; i++) {
574 if(!(x[i]>=0.0)) {v[i]=0.0; continue;}
575 double s=1.0+k*x[i];
576 for(int j=2; j<N; j++) s+=pow(k*x[i], j)/lfactorial(j);
577 v[i]=A*(1.0-s*exp(-k*x[i]));
578 }
579 }
580
581 return(0);
582 }
583
584 /* Exponential functions */
585 if(strcasecmp(fid, "1exp")==0) {
586 if(parNr<2) return(2);
587 double A=p[0], k=p[1];
588 if(fabs(k)<1.0E-20) {
589 for(int i=0; i<sampleNr; i++) v[i]=A*x[i];
590 } else {
591 for(int i=0; i<sampleNr; i++) v[i]=(A/k)*(exp(k*x[i]) - 1.0);
592 }
593 return(0);
594 }
595 if(strcasecmp(fid, "2exp")==0) {
596 if(parNr<4) return(2);
597 for(int i=0; i<sampleNr; i++) v[i]=0.0;
598 for(int n=0; n<=1; n++) {
599 double A=p[n*2], k=p[n*2+1];
600 if(fabs(k)<1.0E-20) {
601 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
602 } else {
603 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
604 }
605 }
606 return(0);
607 }
608 if(strcasecmp(fid, "3exp")==0) {
609 if(parNr<6) return(2);
610 for(int i=0; i<sampleNr; i++) v[i]=0.0;
611 for(int n=0; n<=2; n++) {
612 double A=p[n*2], k=p[n*2+1];
613 if(fabs(k)<1.0E-20) {
614 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
615 } else {
616 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
617 }
618 }
619 return(0);
620 }
621 if(strcasecmp(fid, "4exp")==0) {
622 if(parNr<8) return(2);
623 for(int i=0; i<sampleNr; i++) v[i]=0.0;
624 for(int n=0; n<=3; n++) {
625 double A=p[n*2], k=p[n*2+1];
626 if(fabs(k)<1.0E-20) {
627 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
628 } else {
629 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
630 }
631 }
632 return(0);
633 }
634 if(strcasecmp(fid, "5exp")==0) {
635 if(parNr<10) return(2);
636 for(int i=0; i<sampleNr; i++) v[i]=0.0;
637 for(int n=0; n<=4; n++) {
638 double A=p[n*2], k=p[n*2+1];
639 if(fabs(k)<1.0E-20) {
640 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
641 } else {
642 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
643 }
644 }
645 return(0);
646 }
647
648 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
649 return(10);
650}
651/*****************************************************************************/
652
653/*****************************************************************************/
660int mfEval2ndInt(
662 const char *fid,
664 const int parNr,
666 double *p,
668 const int sampleNr,
670 double *x,
672 double *v,
674 const int verbose
675) {
676 if(verbose>0) {
677 printf("%s('%s', %d, p[], %d, x[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
678 fflush(stdout);
679 }
680 if(parNr<1 || p==NULL || sampleNr<1 || x==NULL || v==NULL) return(1);
681
682 int i;
683 double xt, delayt=0.0;
684
685 if(strcasecmp(fid, "fengm2")==0) {
686 if(parNr<6) return(2);
687 if(parNr>6) delayt=p[6];
688 double A1, A2, A3, L1, L2, L3;
689 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5];
690 if(fabs(L1)<1.0E-06 || fabs(L2)<1.0E-06 || fabs(L3)<1.0E-06) {
691 if(verbose>1) printf("too small L parameter(s)\n");
692 return(3);
693 }
694 for(i=0; i<sampleNr; i++) {
695 xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
696 if(L1!=0.0) {
697 v[i] += (((A1*xt-A2-A3)*L1 - 2.0*A1) * exp(L1*xt) + 2.0*A1)/(L1*L1*L1);
698 v[i] += (A1*xt+A2+A3)/(L1*L1);
699 v[i] += (A2*xt+A3*xt)/L1;
700 }
701 if(L2!=0.0) v[i] -= (A2/(L2*L2))*(1.0 - exp(L2*xt)) + A2*xt/L2;
702 if(L3!=0.0) v[i] -= (A3/(L3*L3))*(1.0 - exp(L3*xt)) + A3*xt/L3;
703 }
704 return(0);
705 }
706 if(strcasecmp(fid, "fengm2e")==0) {
707 if(parNr<8) return(2);
708 if(parNr>8) delayt=p[8];
709 double A1, A2, A3, A4, L1, L2, L3, L4;
710 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5]; A4=p[6]; L4=p[7];
711 if(fabs(L1)<1.0E-15 || fabs(L2)<1.0E-15 || fabs(L3)<1.0E-15 || fabs(L4)<1.0E-15) {
712 if(verbose>1) printf("too small L parameter(s)\n");
713 return(3);
714 }
715 for(i=0; i<sampleNr; i++) {
716 xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
717 v[i] = A4*L1*L1*L1*L2*L2*L3*L3*L4*xt +
718 A4*L1*L1*L1*L2*L2*L3*L3*(1-exp(L4*xt)) +
719 L4*L4*(A3*L1*L1*L1*L2*L2*L3*xt + A3*L1*L1*L1*L2*L2*(1.0-exp(L3*xt)) +
720 L3*L3*(A2*L1*L1*L1*L2*xt + A2*L1*L1*L1*(1.0-exp(L2*xt)) -
721 L2*L2*((A2+A3+A4)*xt*L1*L1 + (A1*xt+A2+A3+A4)*L1 +
722 ((A1*xt-A2-A3-A4)*L1-2.0*A1)*exp(L1*xt) + 2.0*A1)));
723 v[i] /= -(L1*L1*L1*L2*L2*L3*L3*L4*L4);
724 }
725 return(0);
726 }
727 /* Gamma variate function when p[1]==1 */
728 if(strcasecmp(fid, "gammav")==0 && parNr>=3 && fabs(p[1]-1.0)<1.0E-10) {
729 if(parNr<3) return(2);
730 if(parNr>3) delayt=p[3];
731 for(i=0; i<sampleNr; i++) {
732 xt=x[i]-delayt; v[i]=0.0;
733 if(!(xt>0.0) || fabs(p[2])<1.0E-10) continue;
734 v[i] = p[0]*p[2]*p[2]*((p[2]+p[2]+xt)*exp(-xt/p[2]) - p[2] - p[2] + xt);
735 }
736 return(0);
737 }
738
739 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
740 return(10);
741}
742/*****************************************************************************/
743
744/*****************************************************************************/
750int mfEvalFrameY(
752 const char *fid,
754 const int parNr,
756 double *p,
758 const int sampleNr,
760 double *x1,
762 double *x2,
764 double *y,
766 const int verbose
767) {
768 if(verbose>0) {
769 printf("%s('%s', %d, p[], %d, x1[], x2[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
770 fflush(stdout);
771 }
772 if(parNr<1 || p==NULL || sampleNr<1 || x1==NULL || x2==NULL || y==NULL) return(1);
773
774 int i, ret;
775 double xt1, xt2, v, fd, delayt=0.0;
776
777 if(strcasecmp(fid, "fengm2")==0) {
778 if(parNr<6) return(2);
779 if(parNr>6) delayt=p[6];
780 double A1, A2, A3, L1, L2, L3;
781 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5];
782 if(fabs(L1)<1.0E-12 || fabs(L2)<1.0E-12 || fabs(L3)<1.0E-12) {
783 if(verbose>1) printf("too small L parameter(s)\n");
784 return(3);
785 }
786 for(i=0; i<sampleNr; i++) {
787 xt1=x1[i]-delayt; xt2=x2[i]-delayt; y[i]=0.0;
788 if(xt1>xt2) {v=xt1; xt1=xt2; xt2=v;}
789 if(!(xt2>0.0)) continue;
790 fd=xt2-xt1; // requested frame duration
791 if(fd<1.0E-025) {
792 // if frame duration is about zero, then use the method for that
793 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
794 if(ret!=0) return(ret);
795 }
796 if(!(xt1>0.0)) xt1=0.0; // integral start time must not be <0
797 /* Calculate integral between xt1 and xt2 */
798 y[i] = A3*L1*L1*L2*(exp(L3*xt1) - exp(L3*xt2)) +
799 L3*( A2*L1*L1*(exp(L2*xt1) - exp(L2*xt2)) +
800 L2*(((A1*xt1-A2-A3)*L1-A1)*exp(L1*xt1) -
801 ((A1*xt2-A2-A3)*L1-A1)*exp(L1*xt2)
802 )
803 );
804 y[i] /= -(L1*L1*L2*L3);
805 /* Divide integral with the original frame duration */
806 y[i]/=fd;
807 }
808 return(0);
809 }
810 if(strcasecmp(fid, "fengm2e")==0) {
811 if(parNr<8) return(2);
812 if(parNr>8) delayt=p[8];
813 double A1, A2, A3, A4, L1, L2, L3, L4;
814 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5]; A4=p[6]; L4=p[7];
815 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20 || fabs(L3)<1.0E-20 || fabs(L4)<1.0E-20) {
816 if(verbose>1) printf("too small L parameter(s)\n");
817 return(3);
818 }
819 for(i=0; i<sampleNr; i++) {
820 xt1=x1[i]-delayt; xt2=x2[i]-delayt; y[i]=0.0;
821 if(xt1>xt2) {v=xt1; xt1=xt2; xt2=v;}
822 if(!(xt2>0.0)) continue;
823 fd=xt2-xt1; // requested frame duration
824 if(fd<1.0E-025) {
825 // if frame duration is about zero, then use the method for that
826 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
827 if(ret!=0) return(ret);
828 }
829 if(!(xt1>0.0)) xt1=0.0; // integral start time must not be <0
830 /* Calculate integral between xt1 and xt2 */
831 y[i] = A4*L1*L1*L2*L3*(exp(L4*xt1) - exp(L4*xt2)) +
832 L4*( A3*L1*L1*L2*(exp(L3*xt1) - exp(L3*xt2)) +
833 L3*( A2*L1*L1*(exp(L2*xt1) - exp(L2*xt2)) +
834 L2*(((A1*xt1-A2-A3-A4)*L1-A1)*exp(L1*xt1) -
835 ((A1*xt2-A2-A3-A4)*L1-A1)*exp(L1*xt2)
836 )));
837 y[i] /= -(L1*L1*L2*L3*L4);
838 /* Divide integral with the original frame duration */
839 y[i]/=fd;
840 }
841 return(0);
842 }
843 /* Gamma variate function when p[1]==1 */
844 if(strcasecmp(fid, "gammav")==0 && parNr>=3 && fabs(p[1]-1.0)<1.0E-10) {
845 if(parNr<3) return(2);
846 if(parNr>3) delayt=p[3];
847 for(i=0; i<sampleNr; i++) {
848 xt1=x1[i]-delayt; xt2=x2[i]-delayt; y[i]=0.0;
849 if(fabs(p[2])<1.0E-10) continue;
850 if(xt1>xt2) {v=xt1; xt1=xt2; xt2=v;}
851 if(!(xt2>0.0)) continue;
852 fd=xt2-xt1; // requested frame duration
853 if(fd<1.0E-025) {
854 // if frame duration is about zero, then use the method for that
855 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
856 if(ret!=0) return(ret);
857 }
858 if(!(xt1>0.0)) xt1=0.0; // integral start time must not be <0
859 /* Calculate integral between xt1 and xt2 */
860 y[i] = p[0]*p[2]*( (p[2]+xt1)*exp(-xt1/p[2]) - (p[2]+xt2)*exp(-xt2/p[2]) );
861 /* Divide integral with the original frame duration */
862 y[i]/=fd;
863 }
864 return(0);
865 }
866
867 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
868 return(10);
869}
870/*****************************************************************************/
871
872/*****************************************************************************/
879int mfEvalIntToInf(
881 const char *fid,
883 const int parNr,
885 double *p,
887 double x,
889 double *v,
891 const int verbose
892) {
893 if(verbose>0) {
894 printf("%s('%s', %d, p[], %g, y, %d)\n", __func__, fid, parNr, x, verbose);
895 fflush(stdout);
896 }
897 if(parNr<1 || p==NULL || v==NULL) return(1);
898 if(!(x>=0.0)) return(2);
899 *v=0.0;
900
901 /* Exponential functions */
902 if(strcasecmp(fid, "1exp")==0) {
903 if(parNr<2) return(2);
904 double A=p[0], k=p[1];
905 *v=(-A/k)*exp(k*x);
906 return(0);
907 }
908 if(strcasecmp(fid, "2exp")==0) {
909 if(parNr<4) return(2);
910 for(int n=0; n<=1; n++) {
911 double A=p[n*2], k=p[n*2+1];
912 *v+=(-A/k)*exp(k*x);
913 }
914 return(0);
915 }
916 if(strcasecmp(fid, "3exp")==0) {
917 if(parNr<6) return(2);
918 for(int n=0; n<=2; n++) {
919 double A=p[n*2], k=p[n*2+1];
920 *v+=(-A/k)*exp(k*x);
921 }
922 return(0);
923 }
924 if(strcasecmp(fid, "4exp")==0) {
925 if(parNr<8) return(2);
926 for(int n=0; n<=3; n++) {
927 double A=p[n*2], k=p[n*2+1];
928 *v+=(-A/k)*exp(k*x);
929 }
930 return(0);
931 }
932 if(strcasecmp(fid, "5exp")==0) {
933 if(parNr<10) return(2);
934 for(int n=0; n<=4; n++) {
935 double A=p[n*2], k=p[n*2+1];
936 *v+=(-A/k)*exp(k*x);
937 }
938 return(0);
939 }
940
941
942
943 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
944 return(10);
945}
946/*****************************************************************************/
947
948/*****************************************************************************/
mfEvalY
int mfEvalY(const char *fid, const int parNr, double *p, const int sampleNr, double *x, double *y, const int verbose)
Definition func.c:26
mfEvalInt
int mfEvalInt(const char *fid, const int parNr, double *p, const int sampleNr, double *x, double *v, const int verbose)
Definition func.c:418
mfEvalFrameY
int mfEvalFrameY(const char *fid, const int parNr, double *p, const int sampleNr, double *x1, double *x2, double *y, const int verbose)
Definition func.c:750
mfEvalIntToInf
int mfEvalIntToInf(const char *fid, const int parNr, double *p, double x, double *v, const int verbose)
Definition func.c:879
mfEval2ndInt
int mfEval2ndInt(const char *fid, const int parNr, double *p, const int sampleNr, double *x, double *v, const int verbose)
Definition func.c:660
lfactorial
unsigned long long int lfactorial(unsigned long long int n)
Definition intutil.c:63
igam
double igam(double a, double x)
Definition rgamma.c:29
tpcextensions.h
Header file for library libtpcextensions.
tpcfunc.h
Header file for libtpcfunc.
Generated on Fri Sep 19 2025 13:32:18 for TPCCLIB by doxygen 1.12.0