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 const double *p,
34 const int sampleNr,
36 const 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 /* Sum of Surge functions with delay, in form
205 f(x)=p1*(x-p0)*exp(-p2*(x-p0)) + p3*(x-p0)*exp(-p4*(x-p0)) + ... */
206 if(strcasecmp(fid, "dmsurge")==0) {
207 if(parNr<3) return(2);
208 double delayt=p[0];
209 for(int i=0; i<sampleNr; i++) {
210 double xt=x[i]-delayt;
211 y[i]=0.0;
212 if(xt>0.0)
213 for(int pi=2; pi<parNr; pi+=2) {
214 double a, b;
215 a=p[pi-1]; b=p[pi];
216 y[i] += a*xt*exp(-b*xt);
217 }
218 }
219 return(0);
220 }
221
222 /* Surge function in form f(x)=a*x*exp(-b*x)*b^2 to give AUC=a from 0 to infinity.
223 Maximum value is at 1/b. */
224 if(strcasecmp(fid, "surge")==0) {
225 if(parNr<2) return(2);
226 double f=p[0]*(p[1]*p[1]);
227 for(int i=0; i<sampleNr; i++) {
228 double xt=x[i]; y[i]=0.0;
229 if(!(xt>0.0)) continue;
230 y[i] = f*xt*exp(-p[1]*xt);
231 }
232 return(0);
233 }
234
235 /* Surge function in form f(x)=a*x*exp(-b*x). Maximum value is at 1/b. */
236 if(strcasecmp(fid, "tradsurge")==0) {
237 if(parNr<2) return(2);
238 for(int i=0; i<sampleNr; i++) {
239 double xt=x[i]; y[i]=0.0;
240 if(!(xt>0.0)) continue;
241 y[i] = p[0]*xt*exp(-p[1]*xt);
242 }
243 return(0);
244 }
245
246 /* Surge function with recirculation, in form
247 f(x)=a*(x*exp(-b*x) + (c/(b*b))*(1-(b*x+1)*exp(-b*x))). */
248 if(strcasecmp(fid, "surgerecirc")==0) {
249 if(parNr<2) return(2);
250 double c=0.0;
251 if(parNr>=3) c=p[2];
252 for(int i=0; i<sampleNr; i++) {
253 double xt=x[i]; y[i]=0.0;
254 if(!(xt>0.0)) continue;
255 double e=exp(-p[1]*xt);
256 y[i] = xt*e + (c/(p[1]*p[1]))*(1.0 - (p[1]*xt+1.0)*e);
257 y[i] *= p[0];
258 }
259 return(0);
260 }
261
262 /* Surge function with recirculation and delay, in form
263 f(x)=a*((x-d)*exp(-b*(x-d)) + (c/(b*b))*(1-(b*(x-d)+1)*exp(-b*(x-d)))). */
264 if(strcasecmp(fid, "surgerecircd")==0) {
265 if(parNr<2) return(2);
266 double c=0.0; if(parNr>=3) c=p[2];
267 double delayt=0.0; if(parNr>=4) delayt=p[3];
268 for(int i=0; i<sampleNr; i++) {
269 double xt=x[i]-delayt; y[i]=0.0;
270 if(!(xt>0.0)) continue;
271 double e=exp(-p[1]*xt);
272 y[i] = xt*e + (c/(p[1]*p[1]))*(1.0 - (p[1]*xt+1.0)*e);
273 y[i] *= p[0];
274 }
275 return(0);
276 }
277
278 /* Surge function with recirculation, for plasma-to-blood ratio.
279 RBC-to-plasma r(x)=a*(x*exp(-b*x) + (c/(b*b))*(1-(b*x+1)*exp(-b*x)))
280 Plasma-to-blood f(x)=1/(1-HCT*(1-r(x))). */
281 if(strcasecmp(fid, "p2bsrc")==0) {
282 if(parNr<3) return(2);
283 double c=0.0; if(parNr>=4) c=p[3];
284 double HCT=p[0];
285 double a=p[1];
286 double b=p[2];
287 for(int i=0; i<sampleNr; i++) {
288 double xt=x[i]; if(!(xt>0.0)) {y[i]=1.0/(1.0-HCT); continue;}
289 double e=exp(-b*xt);
290 double r=a*(xt*e + (c/(b*b))*(1.0 - (b*xt+1.0)*e));
291 y[i] = 1.0/(1.0-HCT*(1.0-r));
292 }
293 return(0);
294 }
295
296 /* Feng M2 function for plasma-to-blood ratio.
297 RBC-to-plasma ratio r(x) using Feng M2, and then Plasma-to-blood f(x)=1/(1-HCT*(1-r(x))). */
298 if(strcasecmp(fid, "p2bfm2")==0) {
299 if(parNr<3) return(2);
300 double HCT=p[0];
301 double A1=p[1];
302 double L1=p[2];
303 double A2=0.0; if(parNr>3) A2=p[3];
304 double L2=0.0; if(parNr>4) L2=p[4];
305 double A3=0.0; if(parNr>5) A3=p[5];
306 double L3=0.0; if(parNr>6) L3=p[6];
307 for(int i=0; i<sampleNr; i++) {
308 double xt=x[i]; if(!(xt>0.0)) {y[i]=1.0/(1.0-HCT); continue;}
309 double r=(A1*xt - A2 - A3)*exp(-L1*xt) + A2*exp(-L2*xt) + A3*exp(-L3*xt);
310 y[i] = 1.0/(1.0-HCT*(1.0-r));
311 }
312 return(0);
313 }
314
315 /* Surge function for late FDG AIF with time delay */
316 if(strcasecmp(fid, "surgefdgaif")==0) {
317 if(parNr<5) return(2);
318 double delayt=p[0];
319 for(int i=0; i<sampleNr; i++) {
320 double xt=x[i]-delayt; y[i]=0.0;
321 if(!(xt>0.0)) continue;
322 y[i] = p[2]*(p[3]*xt*exp(-p[4]*p[1]*xt) + exp(-p[1]*xt));
323 }
324 return(0);
325 }
326
327 /* Probability density function of Erlang distribution */
328 if(strcasecmp(fid, "erlangpdf")==0) {
329 if(parNr<3) return(2);
330 double A=p[0], k=p[1]; if(!(k>=0.0)) return(2);
331 int N=(int)round(p[2]); if(!(N>0)) return(2);
332 unsigned long long int f=lfactorial(N-1);
333 for(int i=0; i<sampleNr; i++) {
334 y[i]=0.0; if(!(x[i]>=0.0)) continue;
335 y[i] = A*pow(k, N)*pow(x[i], N-1)*exp(-k*x[i])/f;
336 }
337 return(0);
338 }
339
340 /* Rational functions */
341 if(strcasecmp(fid, "ratf11")==0) {
342 if(parNr<4) return(2);
343 for(int i=0; i<sampleNr; i++) {
344 double a=p[0]+p[2]*x[i];
345 double b=p[1]+p[3]*x[i];
346 y[i]=a/b;
347 }
348 return(0);
349 }
350 if(strcasecmp(fid, "ratf21")==0) {
351 if(parNr<5) return(2);
352 for(int i=0; i<sampleNr; i++) {
353 double a=p[0]+p[2]*x[i]+p[4]*x[i]*x[i];
354 double b=p[1]+p[3]*x[i];
355 y[i]=a/b;
356 }
357 return(0);
358 }
359 if(strcasecmp(fid, "ratf22")==0) {
360 if(parNr<6) return(2);
361 for(int i=0; i<sampleNr; i++) {
362 double sqrx=x[i]*x[i];
363 double a=p[0]+p[2]*x[i]+p[4]*sqrx;
364 double b=p[1]+p[3]*x[i]+p[5]*sqrx;
365 y[i]=a/b;
366 }
367 return(0);
368 }
369 if(strcasecmp(fid, "ratf32")==0) {
370 if(parNr<7) return(2);
371 for(int i=0; i<sampleNr; i++) {
372 double sqrx=x[i]*x[i];
373 double a=p[0]+p[2]*x[i]+p[4]*sqrx+p[6]*sqrx*x[i];
374 double b=p[1]+p[3]*x[i]+p[5]*sqrx;
375 y[i]=a/b;
376 }
377 return(0);
378 }
379 if(strcasecmp(fid, "ratf33")==0) {
380 if(parNr<8) return(2);
381 for(int i=0; i<sampleNr; i++) {
382 double sqrx=x[i]*x[i];
383 double cubx=sqrx*x[i];
384 double a=p[0]+p[2]*x[i]+p[4]*sqrx+p[6]*cubx;
385 double b=p[1]+p[3]*x[i]+p[5]*sqrx+p[7]*cubx;
386 y[i]=a/b;
387 }
388 return(0);
389 }
390 if(strcasecmp(fid, "p2brf")==0) {
391 if(parNr<7) return(2);
392 for(int i=0; i<sampleNr; i++) {
393 if(x[i]<=0.0) {
394 y[i]=1.0/(1.0-p[0]);
395 } else {
396 double sqrx=x[i]*x[i];
397 double cubx=sqrx*x[i];
398 double a=0.0+p[1]*x[i]+p[3]*sqrx+p[5]*cubx;
399 double b=1.0+p[2]*x[i]+p[4]*sqrx+p[6]*cubx;
400 y[i]=1.0/(1.0-p[0]*(1.0-(a/b)));
401 }
402 }
403 return(0);
404 }
405
406 /* Hill functions */
407 if(strcasecmp(fid, "mpfhill")==0) {
408 if(parNr<3) return(2);
409 for(int i=0; i<sampleNr; i++) {
410 y[i]=p[0]*pow(x[i], p[1]) / (pow(x[i], p[1]) + p[2]);
411 }
412 return(0);
413 }
414 if(strcasecmp(fid, "p2bhill")==0) {
415 if(parNr<4) return(2);
416 for(int i=0; i<sampleNr; i++) {
417 double cpr=p[1]*pow(x[i], p[2]) / (pow(x[i], p[2]) + p[3]);
418 y[i]=1.0/(1.0-p[0]*(1.0-cpr));
419 }
420 return(0);
421 }
422
423
424 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
425 return(10);
426}
427/*****************************************************************************/
428
429/*****************************************************************************/
436int mfEvalInt(
438 const char *fid,
440 const int parNr,
442 const double *p,
444 const int sampleNr,
446 const double *x,
448 double *v,
450 const int verbose
451) {
452 if(verbose>0) {
453 printf("%s('%s', %d, p[], %d, x[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
454 fflush(stdout);
455 }
456 if(parNr<1 || p==NULL || sampleNr<1 || x==NULL || v==NULL) return(1);
457
458 if(strcasecmp(fid, "fengm2")==0) {
459 if(parNr<6) return(2);
460 double delayt=0.0; if(parNr>6) delayt=p[6];
461 double a;
462 double A1, A2, A3, L1, L2, L3; // lambdas are negative
463 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5];
464 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20 || fabs(L3)<1.0E-20) {
465 if(verbose>1) printf("too small L parameter(s)\n");
466 return(3);
467 }
468 for(int i=0; i<sampleNr; i++) {
469 double xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
470 if(L1!=0.0) {
471 a=exp(L1*xt);
472 v[i] += (A1+L1*(A2+A3))*(1.0-a)/(L1*L1);
473 v[i] += (A1/L1)*xt*a;
474 }
475 if(L2!=0.0) v[i]-=(A2/L2)*(1.0-exp(L2*xt)); else v[i]+=A2*xt;
476 if(L3!=0.0) v[i]-=(A3/L3)*(1.0-exp(L3*xt)); else v[i]+=A3*xt;
477 }
478 return(0);
479 }
480 if(strcasecmp(fid, "fengm2e")==0) {
481 if(parNr<8) return(2);
482 double delayt=0.0; if(parNr>8) delayt=p[8];
483 double a;
484 double A1, A2, A3, A4, L1, L2, L3, L4; // lambdas are negative
485 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5]; A4=p[6]; L4=p[7];
486 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20 || fabs(L3)<1.0E-20 || fabs(L4)<1.0E-20) {
487 if(verbose>1) printf("too small L parameter(s)\n");
488 return(3);
489 }
490 for(int i=0; i<sampleNr; i++) {
491 double xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
492 if(L1!=0.0) {
493 a=exp(L1*xt);
494 v[i] += (A1+L1*(A2+A3+A4))*(1.0-a)/(L1*L1);
495 v[i] += (A1/L1)*xt*a;
496 }
497 if(L2!=0.0) v[i]-=(A2/L2)*(1.0-exp(L2*xt)); else v[i]+=A2*xt;
498 if(L3!=0.0) v[i]-=(A3/L3)*(1.0-exp(L3*xt)); else v[i]+=A3*xt;
499 if(L4!=0.0) v[i]-=(A4/L4)*(1.0-exp(L4*xt)); else v[i]+=A4*xt;
500 }
501 return(0);
502 }
503 if(strcasecmp(fid, "fengm2s")==0) {
504 if(parNr<4) return(2);
505 double delayt=0.0; if(parNr>4) delayt=p[4];
506 double a;
507 double A1, A2, L1, L2; // lambdas are negative
508 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3];
509 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20) {
510 if(verbose>1) printf("too small L parameter(s)\n");
511 return(3);
512 }
513 for(int i=0; i<sampleNr; i++) {
514 double xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
515 if(L1!=0.0) {
516 a=exp(L1*xt);
517 v[i] += (A1+L1*A2)*(1.0-a)/(L1*L1);
518 v[i] += (A1/L1)*xt*a;
519 }
520 if(L2!=0.0) v[i]-=(A2/L2)*(1.0-exp(L2*xt)); else v[i]+=A2*xt;
521 }
522 return(0);
523 }
524 /* Gamma variate function when p[1]==1 */
525 if(strcasecmp(fid, "gammav")==0 && parNr>=3 && fabs(p[1]-1.0)<1.0E-10) {
526 if(parNr<3) return(2);
527 double delayt=0.0; if(parNr>3) delayt=p[3];
528 for(int i=0; i<sampleNr; i++) {
529 double xt=x[i]-delayt; v[i]=0.0;
530 if(!(xt>0.0) || fabs(p[2])<1.0E-10) continue;
531 v[i] = p[0]*p[2]*(p[2] - (p[2]+xt)*exp(-xt/p[2]));
532 }
533 return(0);
534 }
535
536 /* Integral of sum of Surge functions with delay, in form
537 f(x)=p1*(x-p0)*exp(-p2*(x-p0)) + p3*(x-p0)*exp(-p4*(x-p0)) + ... */
538 if(strcasecmp(fid, "dmsurge")==0) {
539 if(parNr<3) return(2);
540 for(int i=0; i<sampleNr; i++) {
541 double xt=x[i]-p[0];
542 v[i]=0.0;
543 if(xt>0.0)
544 for(int pi=2; pi<parNr; pi+=2) {
545 double a, b;
546 a=p[pi-1]; b=p[pi];
547 if(a>0.0) {
548 if(fabs(b)>1.0E-12)
549 v[i] += a*(1.0-(b*xt+1.0)*exp(-b*xt))/(b*b);
550 else
551 v[i] += a*xt*xt;
552 }
553 }
554 }
555 return(0);
556 }
557
558 /* Surge function in form f(x)=a*x*exp(-b*x)*b^2, to give AUC=a from 0 to infinity. */
559 if(strcasecmp(fid, "surge")==0) {
560 if(parNr<2) return(2);
561 for(int i=0; i<sampleNr; i++) {
562 double xt=x[i]; v[i]=0.0;
563 if(!(xt>0.0)) continue;
564 v[i] = p[0]*(1.0 - (p[1]*xt + 1.0)*exp(-p[1]*xt));
565 }
566 return(0);
567 }
568
569 /* Surge function in form f(x)=a*x*exp(-b*x). Maximum value is at 1/b. */
570 if(strcasecmp(fid, "tradsurge")==0) {
571 if(parNr<2) return(2);
572 double f=p[0]/(p[1]*p[1]);
573 for(int i=0; i<sampleNr; i++) {
574 double xt=x[i]; v[i]=0.0;
575 if(!(xt>0.0)) continue;
576 v[i] = f*(1.0 - (p[1]*xt + 1.0)*exp(-p[1]*xt));
577 }
578 return(0);
579 }
580
581 /* Surge function for late FDG AIF with time delay */
582 if(strcasecmp(fid, "surgefdgaif")==0) {
583 if(parNr<5) return(2);
584 double delayt=p[0];
585 for(int i=0; i<sampleNr; i++) {
586 double xt=x[i]-delayt; v[i]=0.0;
587 if(!(xt>0.0)) continue;
588 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]);
589 v[i] += (1.0-exp(-p[1]*xt))/p[1];
590 v[i] *= p[2];
591 }
592 return(0);
593 }
594
595 /* Probability density function of Erlang distribution */
596 if(strcasecmp(fid, "erlangpdf")==0 && p[2]<20.5) { // Supports only N=0,1,2,...,20
597 if(parNr<3) return(2);
598 double A=p[0], k=p[1]; if(!(k>=0.0)) return(2);
599 int N=(int)round(p[2]); if(!(N>=0) || N>20) return(2);
600 if(N==0)
601 for(int i=0; i<sampleNr; i++) {
602 if(x[i]>=0.0) v[i]=A; else v[i]=0.0;
603 }
604 else if(N==1)
605 for(int i=0; i<sampleNr; i++) {
606 if(x[i]>=0.0) v[i]=A*(1.0-exp(-k*x[i])); else v[i]=0.0;
607 }
608 else if(N==2)
609 for(int i=0; i<sampleNr; i++) {
610 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;
611 }
612 else {
613 for(int i=0; i<sampleNr; i++) {
614 if(!(x[i]>=0.0)) {v[i]=0.0; continue;}
615 double s=1.0+k*x[i];
616 for(int j=2; j<N; j++) s+=pow(k*x[i], j)/lfactorial(j);
617 v[i]=A*(1.0-s*exp(-k*x[i]));
618 }
619 }
620
621 return(0);
622 }
623
624 /* Exponential functions */
625 if(strcasecmp(fid, "1exp")==0) {
626 if(parNr<2) return(2);
627 double A=p[0], k=p[1];
628 if(fabs(k)<1.0E-20) {
629 for(int i=0; i<sampleNr; i++) v[i]=A*x[i];
630 } else {
631 for(int i=0; i<sampleNr; i++) v[i]=(A/k)*(exp(k*x[i]) - 1.0);
632 }
633 return(0);
634 }
635 if(strcasecmp(fid, "2exp")==0) {
636 if(parNr<4) return(2);
637 for(int i=0; i<sampleNr; i++) v[i]=0.0;
638 for(int n=0; n<=1; n++) {
639 double A=p[n*2], k=p[n*2+1];
640 if(fabs(k)<1.0E-20) {
641 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
642 } else {
643 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
644 }
645 }
646 return(0);
647 }
648 if(strcasecmp(fid, "3exp")==0) {
649 if(parNr<6) return(2);
650 for(int i=0; i<sampleNr; i++) v[i]=0.0;
651 for(int n=0; n<=2; n++) {
652 double A=p[n*2], k=p[n*2+1];
653 if(fabs(k)<1.0E-20) {
654 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
655 } else {
656 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
657 }
658 }
659 return(0);
660 }
661 if(strcasecmp(fid, "4exp")==0) {
662 if(parNr<8) return(2);
663 for(int i=0; i<sampleNr; i++) v[i]=0.0;
664 for(int n=0; n<=3; n++) {
665 double A=p[n*2], k=p[n*2+1];
666 if(fabs(k)<1.0E-20) {
667 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
668 } else {
669 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
670 }
671 }
672 return(0);
673 }
674 if(strcasecmp(fid, "5exp")==0) {
675 if(parNr<10) return(2);
676 for(int i=0; i<sampleNr; i++) v[i]=0.0;
677 for(int n=0; n<=4; n++) {
678 double A=p[n*2], k=p[n*2+1];
679 if(fabs(k)<1.0E-20) {
680 for(int i=0; i<sampleNr; i++) v[i]+=A*x[i];
681 } else {
682 for(int i=0; i<sampleNr; i++) v[i]+=(A/k)*(exp(k*x[i]) - 1.0);
683 }
684 }
685 return(0);
686 }
687
688 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
689 return(10);
690}
691/*****************************************************************************/
692
693/*****************************************************************************/
700int mfEval2ndInt(
702 const char *fid,
704 const int parNr,
706 const double *p,
708 const int sampleNr,
710 const double *x,
712 double *v,
714 const int verbose
715) {
716 if(verbose>0) {
717 printf("%s('%s', %d, p[], %d, x[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
718 fflush(stdout);
719 }
720 if(parNr<1 || p==NULL || sampleNr<1 || x==NULL || v==NULL) return(1);
721
722 int i;
723 double xt, delayt=0.0;
724
725 if(strcasecmp(fid, "fengm2")==0) {
726 if(parNr<6) return(2);
727 if(parNr>6) delayt=p[6];
728 double A1, A2, A3, L1, L2, L3;
729 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5];
730 if(fabs(L1)<1.0E-06 || fabs(L2)<1.0E-06 || fabs(L3)<1.0E-06) {
731 if(verbose>1) printf("too small L parameter(s)\n");
732 return(3);
733 }
734 for(i=0; i<sampleNr; i++) {
735 xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
736 if(L1!=0.0) {
737 v[i] += (((A1*xt-A2-A3)*L1 - 2.0*A1) * exp(L1*xt) + 2.0*A1)/(L1*L1*L1);
738 v[i] += (A1*xt+A2+A3)/(L1*L1);
739 v[i] += (A2*xt+A3*xt)/L1;
740 }
741 if(L2!=0.0) v[i] -= (A2/(L2*L2))*(1.0 - exp(L2*xt)) + A2*xt/L2;
742 if(L3!=0.0) v[i] -= (A3/(L3*L3))*(1.0 - exp(L3*xt)) + A3*xt/L3;
743 }
744 return(0);
745 }
746 if(strcasecmp(fid, "fengm2e")==0) {
747 if(parNr<8) return(2);
748 if(parNr>8) delayt=p[8];
749 double A1, A2, A3, A4, L1, L2, L3, L4;
750 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5]; A4=p[6]; L4=p[7];
751 if(fabs(L1)<1.0E-15 || fabs(L2)<1.0E-15 || fabs(L3)<1.0E-15 || fabs(L4)<1.0E-15) {
752 if(verbose>1) printf("too small L parameter(s)\n");
753 return(3);
754 }
755 for(i=0; i<sampleNr; i++) {
756 xt=x[i]-delayt; v[i]=0.0; if(!(xt>0.0)) continue;
757 v[i] = A4*L1*L1*L1*L2*L2*L3*L3*L4*xt +
758 A4*L1*L1*L1*L2*L2*L3*L3*(1-exp(L4*xt)) +
759 L4*L4*(A3*L1*L1*L1*L2*L2*L3*xt + A3*L1*L1*L1*L2*L2*(1.0-exp(L3*xt)) +
760 L3*L3*(A2*L1*L1*L1*L2*xt + A2*L1*L1*L1*(1.0-exp(L2*xt)) -
761 L2*L2*((A2+A3+A4)*xt*L1*L1 + (A1*xt+A2+A3+A4)*L1 +
762 ((A1*xt-A2-A3-A4)*L1-2.0*A1)*exp(L1*xt) + 2.0*A1)));
763 v[i] /= -(L1*L1*L1*L2*L2*L3*L3*L4*L4);
764 }
765 return(0);
766 }
767 /* Gamma variate function when p[1]==1 */
768 if(strcasecmp(fid, "gammav")==0 && parNr>=3 && fabs(p[1]-1.0)<1.0E-10) {
769 if(parNr<3) return(2);
770 if(parNr>3) delayt=p[3];
771 for(i=0; i<sampleNr; i++) {
772 xt=x[i]-delayt; v[i]=0.0;
773 if(!(xt>0.0) || fabs(p[2])<1.0E-10) continue;
774 v[i] = p[0]*p[2]*p[2]*((p[2]+p[2]+xt)*exp(-xt/p[2]) - p[2] - p[2] + xt);
775 }
776 return(0);
777 }
778
779 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
780 return(10);
781}
782/*****************************************************************************/
783
784/*****************************************************************************/
790int mfEvalFrameY(
792 const char *fid,
794 const int parNr,
796 const double *p,
798 const int sampleNr,
800 const double *x1,
802 const double *x2,
804 double *y,
806 const int verbose
807) {
808 if(verbose>0) {
809 printf("%s('%s', %d, p[], %d, x1[], x2[], y[], %d)\n", __func__, fid, parNr, sampleNr, verbose);
810 fflush(stdout);
811 }
812 if(parNr<1 || p==NULL || sampleNr<1 || x1==NULL || x2==NULL || y==NULL) return(1);
813 if(verbose>10) {
814 printf("p[] := %g", p[0]);
815 for(int i=1; i<parNr; i++) printf(", %g", p[i]);
816 printf("\n");
817 }
818
819 int i, ret;
820 double xt1, xt2, v, fd, delayt=0.0;
821
822 if(strcasecmp(fid, "fengm2")==0) {
823 if(parNr<6) return(2);
824 if(parNr>6) delayt=p[6];
825 double A1, A2, A3, L1, L2, L3;
826 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5];
827 if(fabs(L1)<1.0E-12 || fabs(L2)<1.0E-12 || fabs(L3)<1.0E-12) {
828 if(verbose>1) printf("too small L parameter(s)\n");
829 return(3);
830 }
831 for(i=0; i<sampleNr; i++) {
832 xt1=x1[i]-delayt; xt2=x2[i]-delayt; y[i]=0.0;
833 if(xt1>xt2) {v=xt1; xt1=xt2; xt2=v;}
834 if(!(xt2>0.0)) continue;
835 fd=xt2-xt1; // requested frame duration
836 if(fd<1.0E-025) {
837 // if frame duration is about zero, then use the method for that
838 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
839 if(ret!=0) return(ret);
840 continue;
841 }
842 if(!(xt1>0.0)) xt1=0.0; // integral start time must not be <0
843 /* Calculate integral between xt1 and xt2 */
844 y[i] = A3*L1*L1*L2*(exp(L3*xt1) - exp(L3*xt2)) +
845 L3*( A2*L1*L1*(exp(L2*xt1) - exp(L2*xt2)) +
846 L2*(((A1*xt1-A2-A3)*L1-A1)*exp(L1*xt1) -
847 ((A1*xt2-A2-A3)*L1-A1)*exp(L1*xt2)
848 )
849 );
850 y[i] /= -(L1*L1*L2*L3);
851 /* Divide integral with the original frame duration */
852 y[i]/=fd;
853 }
854 return(0);
855 }
856
857 if(strcasecmp(fid, "fengm2e")==0) {
858 if(parNr<8) return(2);
859 if(parNr>8) delayt=p[8];
860 double A1, A2, A3, A4, L1, L2, L3, L4;
861 A1=p[0]; L1=p[1]; A2=p[2]; L2=p[3]; A3=p[4]; L3=p[5]; A4=p[6]; L4=p[7];
862 if(fabs(L1)<1.0E-20 || fabs(L2)<1.0E-20 || fabs(L3)<1.0E-20 || fabs(L4)<1.0E-20) {
863 if(verbose>1) printf("too small L parameter(s)\n");
864 return(3);
865 }
866 for(i=0; i<sampleNr; i++) {
867 xt1=x1[i]-delayt; xt2=x2[i]-delayt; y[i]=0.0;
868 if(xt1>xt2) {v=xt1; xt1=xt2; xt2=v;}
869 if(!(xt2>0.0)) continue;
870 fd=xt2-xt1; // requested frame duration
871 if(fd<1.0E-025) {
872 // if frame duration is about zero, then use the method for that
873 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
874 if(ret!=0) return(ret);
875 continue;
876 }
877 if(!(xt1>0.0)) xt1=0.0; // integral start time must not be <0
878 /* Calculate integral between xt1 and xt2 */
879 y[i] = A4*L1*L1*L2*L3*(exp(L4*xt1) - exp(L4*xt2)) +
880 L4*( A3*L1*L1*L2*(exp(L3*xt1) - exp(L3*xt2)) +
881 L3*( A2*L1*L1*(exp(L2*xt1) - exp(L2*xt2)) +
882 L2*(((A1*xt1-A2-A3-A4)*L1-A1)*exp(L1*xt1) -
883 ((A1*xt2-A2-A3-A4)*L1-A1)*exp(L1*xt2)
884 )));
885 y[i] /= -(L1*L1*L2*L3*L4);
886 /* Divide integral with the original frame duration */
887 y[i]/=fd;
888 }
889 return(0);
890 }
891 /* Gamma variate function when p[1]==1 */
892 if(strcasecmp(fid, "gammav")==0 && parNr>=3 && fabs(p[1]-1.0)<1.0E-10) {
893 if(parNr<3) return(2);
894 if(parNr>3) delayt=p[3];
895 for(i=0; i<sampleNr; i++) {
896 xt1=x1[i]-delayt; xt2=x2[i]-delayt; y[i]=0.0;
897 if(fabs(p[2])<1.0E-10) continue;
898 if(xt1>xt2) {v=xt1; xt1=xt2; xt2=v;}
899 if(!(xt2>0.0)) continue;
900 fd=xt2-xt1; // requested frame duration
901 if(fd<1.0E-025) {
902 // if frame duration is about zero, then use the method for that
903 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
904 if(ret!=0) return(ret);
905 continue;
906 }
907 if(!(xt1>=0.0)) xt1=0.0; // integral start time must not be <0
908 /* Calculate integral between xt1 and xt2 */
909 y[i] = p[0]*p[2]*( (p[2]+xt1)*exp(-xt1/p[2]) - (p[2]+xt2)*exp(-xt2/p[2]) );
910 /* Divide integral with the original frame duration */
911 y[i]/=fd;
912 }
913 return(0);
914 }
915
916
917 /* Sum of Surge functions with delay, in form
918 f(x)=p1*(x-p0)*exp(-p2*(x-p0)) + p3*(x-p0)*exp(-p4*(x-p0)) + ... */
919 if(strcasecmp(fid, "dmsurge")==0) {
920 if(parNr<3) return(2);
921 delayt=p[0];
922 for(int i=0; i<sampleNr; i++) {
923 y[i]=0.0;
924 fd=x2[i]-x1[i]; // requested frame duration
925 if(verbose>11) printf(" fd[%d]=%g\n", i, fd);
926 if(!(fd>=0.0)) continue;
927 xt1=x1[i]-delayt; xt2=x2[i]-delayt;
928 if(!(xt2>0.0)) continue;
929 if(xt1<0.0) xt1=0.0; // integral start time must not be <0
930 if(fd<1.0E-025) {
931 // if frame duration is about zero, then simply calculate f(x) for that
932 ret=mfEvalY(fid, parNr, p, 1, x2+i, y+i, verbose);
933 if(ret!=0) return(ret);
934 if(verbose>11) printf(" y[%d]=%g\n", i, y[i]);
935 continue;
936 }
937 /* Calculate integral between xt1 and xt2 */
938 for(int pi=2; pi<parNr; pi+=2) {
939 double a, b;
940 a=p[pi-1]; b=p[pi];
941 if(!(a>0.0)) {
942 } else if(fabs(b)<1.0E-10) { // exp(0)=1
943 y[i] += a*(xt2*xt2 - xt1*xt1);
944 } else { // divider (b*b) must not be zero
945 y[i] += a*((b*xt1+1.0)*exp(-b*xt1) - (b*xt2+1.0)*exp(-b*xt2))/(b*b);
946 }
947 if(verbose>11) printf(" a=%g b=%g xt1=%g xt2=%g y[%d]=%g\n", a, b, xt1, xt2, i, y[i]);
948 }
949 /* Divide integral with the original frame duration */
950 y[i]/=fd;
951 if(verbose>11) printf(" y[%d]=%g\n", i, y[i]);
952 }
953 return(0);
954 }
955
956
957 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
958 return(10);
959}
960/*****************************************************************************/
961
962/*****************************************************************************/
969int mfEvalIntToInf(
971 const char *fid,
973 const int parNr,
975 const double *p,
977 const double x,
979 double *v,
981 const int verbose
982) {
983 if(verbose>0) {
984 printf("%s('%s', %d, p[], %g, y, %d)\n", __func__, fid, parNr, x, verbose);
985 fflush(stdout);
986 }
987 if(parNr<1 || p==NULL || v==NULL) return(1);
988 if(!(x>=0.0)) return(2);
989 *v=0.0;
990
991 /* Exponential functions */
992 if(strcasecmp(fid, "1exp")==0) {
993 if(parNr<2) return(2);
994 double A=p[0], k=p[1];
995 *v=(-A/k)*exp(k*x);
996 return(0);
997 }
998 if(strcasecmp(fid, "2exp")==0) {
999 if(parNr<4) return(2);
1000 for(int n=0; n<=1; n++) {
1001 double A=p[n*2], k=p[n*2+1];
1002 *v+=(-A/k)*exp(k*x);
1003 }
1004 return(0);
1005 }
1006 if(strcasecmp(fid, "3exp")==0) {
1007 if(parNr<6) return(2);
1008 for(int n=0; n<=2; n++) {
1009 double A=p[n*2], k=p[n*2+1];
1010 *v+=(-A/k)*exp(k*x);
1011 }
1012 return(0);
1013 }
1014 if(strcasecmp(fid, "4exp")==0) {
1015 if(parNr<8) return(2);
1016 for(int n=0; n<=3; n++) {
1017 double A=p[n*2], k=p[n*2+1];
1018 *v+=(-A/k)*exp(k*x);
1019 }
1020 return(0);
1021 }
1022 if(strcasecmp(fid, "5exp")==0) {
1023 if(parNr<10) return(2);
1024 for(int n=0; n<=4; n++) {
1025 double A=p[n*2], k=p[n*2+1];
1026 *v+=(-A/k)*exp(k*x);
1027 }
1028 return(0);
1029 }
1030
1031
1032
1033 if(verbose>1) printf("function '%s' not supported by %s()\n", fid, __func__);
1034 return(10);
1035}
1036/*****************************************************************************/
1037
1038/*****************************************************************************/
mfEval2ndInt
int mfEval2ndInt(const char *fid, const int parNr, const double *p, const int sampleNr, const double *x, double *v, const int verbose)
Definition func.c:700
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
mfEvalInt
int mfEvalInt(const char *fid, const int parNr, const double *p, const int sampleNr, const double *x, double *v, const int verbose)
Definition func.c:436
mfEvalIntToInf
int mfEvalIntToInf(const char *fid, const int parNr, const double *p, const double x, double *v, const int verbose)
Definition func.c:969
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
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 for TPCCLIB by doxygen 1.15.0