TPCCLIB
Loading...
Searching...
No Matches
mathfunc.c
Go to the documentation of this file.
1
5int MATHFUNC_TEST;
6char fiterrmsg[64];
7/*****************************************************************************/
8#include "libtpccurveio.h"
9#include <unistd.h>
10/*****************************************************************************/
11#include "libtpcmodel.h"
12/*****************************************************************************/
13
14/*****************************************************************************/
18void fitEmpty(
20 FIT *fit
21) {
22 if(fit==NULL) return;
23 if(fit->_voidataNr>0) {
24 free((char*)(fit->voi));
25 fit->_voidataNr=0;
26 }
27 fit->voiNr=0;
28 fit->datafile[0]=fit->studynr[0]=fit->unit[0]=fit->program[0]=(char)0;
29 fit->timeunit=0;
30 fit->time=0;
31}
32/*****************************************************************************/
33
34/*****************************************************************************/
38void fitInit(
39 FIT *fit
40) {
41 if(fit==NULL) return;
42 memset(fit, 0, sizeof(FIT));
43 fit->_voidataNr=0; fit->voiNr=0;
44}
45/*****************************************************************************/
46
47/*****************************************************************************/
54int fitWrite(
56 FIT *fit,
58 char *filename
59) {
60 int i, j, n, savedNr=0;
61 char tmp[1024], is_stdout=0;
62 FILE *fp;
63
64
65 if(MATHFUNC_TEST>0) printf("fitWrite(FIT, %s)\n", filename);
66 /* Check that there is some data to write */
67 if(fit==NULL || fit->voiNr<1) {strcpy(fiterrmsg, "no data"); return 1;}
68 for(i=0; i<fit->voiNr; i++)
69 if(!isnan(fit->voi[i].wss) && fit->voi[i].type>0) savedNr++;
70 if(savedNr<1) {strcpy(fiterrmsg, "no fitted data"); return 1;}
71
72 /* Check if writing to stdout */
73 if(!strcmp(filename, "stdout")) is_stdout=1;
74 if(MATHFUNC_TEST>1) {
75 if(is_stdout) printf(" output is stdout()\n");
76 else printf(" output in file\n");
77 }
78
79 /* Check if file exists; backup, if necessary */
80 if(!is_stdout) (void)backupExistingFile(filename, NULL, NULL);
81
82 /* Open output file */
83 if(is_stdout) fp=(FILE*)stdout;
84 else {
85 if(MATHFUNC_TEST>2) printf("opening file %s for write\n", filename);
86 if((fp = fopen(filename, "w")) == NULL) {
87 strcpy(fiterrmsg, "cannot open file"); return 2;
88 }
89 }
90
91 /* Fit file format */
92 n=fprintf(fp, "%-11.11s %s\n", FIT_VER, fit->program);
93 if(n==0) {
94 strcpy(fiterrmsg, "disk full");
95 if(!is_stdout) fclose(fp);
96 return 3;
97 }
98 /* Write fit date and time */
99 if(!ctime_r_int(&fit->time, tmp)) strcpy(tmp, "");
100 fprintf(fp, "Date:\t%s\n", tmp);
101 /* Write the name of the original datafile */
102 fprintf(fp, "Data file:\t%s\n", fit->datafile);
103 /* Write the studynr */
104 if(fit->studynr[0]) fprintf(fp, "Study number:\t%s\n", fit->studynr);
105 /* Write the 'activity' unit */
106 fprintf(fp, "Data unit:\t%s\n", fit->unit);
107 /* Write the time unit */
108 if(fit->timeunit==TUNIT_UM || fit->timeunit==TUNIT_MM)
109 fprintf(fp, "Distance unit:\t%s\n", petTunit(fit->timeunit));
110 else
111 fprintf(fp, "Time unit:\t%s\n", petTunit(fit->timeunit));
112
113 /* Write the voiNr to be saved */
114 fprintf(fp, "Nr of VOIs:\t%d\n", savedNr);
115 /* Write the Fit title */
116 fprintf(fp, "%s %s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n",
117 "Region", "Plane", "Start", "End", "dataNr", "WSS", "parNr", "Type", "Parameters");
118 /* Write regional fits */
119 for(i=0; i<fit->voiNr; i++) {
120 if(isnan(fit->voi[i].wss) || fit->voi[i].type<=0) continue;
121 if(fit->voi[i].voiname[0]) strcpy(tmp, fit->voi[i].voiname);
122 else strcpy(tmp, ".");
123 fprintf(fp, "%.*s ", MAX_REGIONSUBNAME_LEN, tmp);
124 if(fit->voi[i].hemisphere[0]) strcpy(tmp, fit->voi[i].hemisphere);
125 else strcpy(tmp, ".");
126 fprintf(fp, "%.*s ", MAX_REGIONSUBNAME_LEN, tmp);
127 if(fit->voi[i].place[0]) strcpy(tmp, fit->voi[i].place);
128 else strcpy(tmp, ".");
129 fprintf(fp, "%.*s\t", MAX_REGIONSUBNAME_LEN, tmp);
130 fprintf(fp, "%.3f\t%.3f\t%d\t%.2E\t%d\t%04d",
131 fit->voi[i].start, fit->voi[i].end, fit->voi[i].dataNr,
132 fit->voi[i].wss, fit->voi[i].parNr, fit->voi[i].type );
133 for(j=0; j<fit->voi[i].parNr; j++) fprintf(fp, "\t%.6E", fit->voi[i].p[j]);
134 fprintf(fp, "\n");
135 }
136
137 /* Close file */
138 if(!is_stdout) {
139 if(MATHFUNC_TEST>2) printf("closing file %s\n", filename);
140 fflush(fp); fclose(fp);
141 }
142 strcpy(fiterrmsg, "");
143
144 if(MATHFUNC_TEST>0) printf("done with fitWrite()\n");
145 return(0);
146}
147/*****************************************************************************/
148
149/*****************************************************************************/
154int fitSetmem(
156 FIT *fit,
158 int voiNr
159) {
160 /* Check that there is something to do */
161 if(fit==NULL || voiNr<1) return 1;
162
163 /* Clear previous data, but only if necessary */
164 if(fit->_voidataNr>0 || fit->voiNr>0) fitEmpty(fit);
165
166 /* Allocate memory for regional curves */
167 fit->voi=(FitVOI*)calloc(voiNr, sizeof(FitVOI));
168 if(fit->voi==NULL) return 2;
169 fit->_voidataNr=voiNr;
170
171 return 0;
172}
173/*****************************************************************************/
174
175/*****************************************************************************/
180void fitPrint(
182 FIT *fit
183) {
184 if(fit==NULL) {printf("FIT = Null\n"); return;}
185 if(MATHFUNC_TEST>0) printf("Number of curves: %d\n", fit->voiNr);
186 if(MATHFUNC_TEST>0) printf("_voidataNr = %d\n", fit->_voidataNr);
187 fitWrite(fit, "stdout");
188}
189/*****************************************************************************/
190
191/*****************************************************************************/
196int fitRead(
198 char *filename,
200 FIT *fit,
202 int verbose
203) {
204 char *cptr, line[1024], *lptr;
205 int i, ri, n, type=0;
206 struct tm st;
207
208
209 if(verbose>0) printf("fitRead(%s)\n", filename);
210 if(fit==NULL) return 1;
211 /* Empty data */
212 fitEmpty(fit);
213
214 /* Open file */
215 FILE *fp=fopen(filename, "r");
216 if(fp==NULL) {strcpy(fiterrmsg, "cannot open file"); return 1;}
217
218 /* Read data */
219 strcpy(fiterrmsg, "wrong format");
220 /* Read file type and program name */
221 while(fgets(line, 1024, fp)!=NULL) {
222 /* Ignore empty and comment lines */
223 if(strlen(line)<4 || line[0]=='#') continue;
224 /* Check for id string */
225 if(!strncmp(line, FIT_VER, strlen(FIT_VER))) type=1; else type=0;
226 break;
227 }
228 if(type!=1) {fclose(fp); return 3;}
229 lptr=line; cptr=strtok(lptr, " \t\n\r");
230 if(cptr!=NULL) cptr=strtok(NULL, " \t\n\r");
231 if(cptr!=NULL && strlen(cptr)<1024) strcpy(fit->program, cptr);
232#if(1)
233 /* Read optional fields until "Nr of VOIs:" */
234 while(1) {
235 /* Read title field */
236 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
237 /* Stop when "Nr of VOIs" is reached */
238 if(strncasecmp(line, "Nr of VOIs:", 11)==0) break;
239 /* Read fit date and time */
240 if(strncasecmp(line, "Date:", 5)==0) {
241 cptr=line+5; while(*cptr && !isdigit(*cptr)) cptr++;
242 if(get_datetime(cptr, &st, verbose-3)==0) fit->time=timegm(&st);
243 continue;
244 }
245 /* Read the name of the original datafile */
246 if(strncasecmp(line, "Data file:", 10)==0) {
247 lptr=&line[10]; cptr=strtok(lptr, " \t\n\r");
248 if(cptr!=NULL && strlen(cptr)<FILENAME_MAX) strcpy(fit->datafile, cptr);
249 continue;
250 }
251 /* Read study number */
252 if(strncasecmp(line, "Study number:", 13)==0) {
253 lptr=&line[13]; cptr=strtok(lptr, " \t\n\r");
254 if(cptr!=NULL && strlen(cptr)<MAX_STUDYNR_LEN+1) strcpy(fit->studynr, cptr);
255 continue;
256 }
257 /* Read the activity unit */
258 if(strncasecmp(line, "Data unit:", 10)==0) {
259 lptr=&line[10]; cptr=strtok(lptr, " \t\n\r");
260 if(cptr!=NULL && strlen(cptr)<1024) strcpy(fit->unit, cptr);
261 continue;
262 }
263 /* Read the time unit */
264 if(strncasecmp(line, "Time unit:", 10)==0) {
265 lptr=&line[10]; cptr=strtok(lptr, " \t\n\r");
266 fit->timeunit=petTunitId(cptr);
267 continue;
268 }
269 if(strncasecmp(line, "Distance unit:", 14)==0) {
270 lptr=&line[14]; cptr=strtok(lptr, " \t\n\r");
271 fit->timeunit=petTunitId(cptr);
272 continue;
273 }
274 fclose(fp); return 8;
275 }
276 /* Read the nr of regions */
277 if(strncasecmp(line, "Nr of VOIs:", 11)) {fclose(fp); return 8;}
278 lptr=&line[11]; cptr=strtok(lptr, " \t\n\r");
279 n=atoi(cptr); if(n<1 || n>32000) {fclose(fp); return 8;}
280#else
281 /* Read fit date and time */
282 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
283 if(strncasecmp(line, "Date:", 5)) {fclose(fp); return 4;}
284 cptr=line+5; while(*cptr && !isdigit(*cptr)) cptr++;
285 if(get_datetime(cptr, &st, verbose-3)==0) fit->time=timegm(&st);
286 /* Read the name of the original datafile */
287 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
288 if(strncasecmp(line, "Data file:", 10)) {fclose(fp); return 5;}
289 lptr=&line[10]; cptr=strtok(lptr, " \t\n\r");
290 if(cptr!=NULL && strlen(cptr)<FILENAME_MAX) strcpy(fit->datafile, cptr);
291 /* Read the activity unit */
292 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
293 if(strncasecmp(line, "Data unit:", 10)) {fclose(fp); return 6;}
294 lptr=&line[10]; cptr=strtok(lptr, " \t\n\r");
295 if(cptr!=NULL && strlen(cptr)<1024) strcpy(fit->unit, cptr);
296 /* Read the time unit */
297 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
298 if(strncasecmp(line, "Time unit:", 10)==0) lptr=&line[10];
299 else if(strncasecmp(line, "Distance unit:", 14)==0) lptr=&line[14];
300 else {fclose(fp); return 7;}
301 cptr=strtok(lptr, " \t\n\r");
302 fit->timeunit=petTunitId(cptr);
303 /* Read the nr of regions */
304 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
305 if(strncasecmp(line, "Nr of VOIs:", 11)) {fclose(fp); return 8;}
306 lptr=&line[11]; cptr=strtok(lptr, " \t\n\r");
307 n=atoi(cptr); if(n<1 || n>32000) {fclose(fp); return 8;}
308#endif
309 /* Allocate memory for regions */
310 if(fitSetmem(fit, n)) {strcpy(fiterrmsg, "out of memory"); fclose(fp); return 9;}
311 fit->voiNr=n;
312 /* Read (and ignore) title line */
313 strcpy(fiterrmsg, "wrong format");
314 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
315 if(strncasecmp(line, "Region", 6)) {fclose(fp); return 10;}
316 /* Read regional data */
317 for(ri=0; ri<fit->voiNr; ri++) {
318 /* Read line of data */
319 line[0]=(char)0;
320 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
321 if(!line[0]) break;
322 int pindx=0; char seps[8];
323 int sn=strTokenNr(line, " \t\n\r"); if(sn<8) {fclose(fp); fitEmpty(fit); return 11;}
324 int tn=strTokenNr(line, "\t\n\r");
325 if(tn==sn || tn==sn-1 || tn==sn-2) { // tab as separator
326 strcpy(seps, "\t\n\r"); n=tn;
327 char *s=strTokenDup(line, seps, NULL);
328 char *cptr=s;
329 int i=0, n=strlen(s);
330 strReplaceChar(cptr, ' ', (char)0);
331 strlcpy(fit->voi[ri].voiname, cptr, MAX_REGIONSUBNAME_LEN+1);
332 i+=strlen(fit->voi[ri].voiname);
333 cptr=s+i; if(i<n) {cptr++; i++;}
334 strReplaceChar(cptr, ' ', (char)0);
335 strlcpy(fit->voi[ri].hemisphere, cptr, MAX_REGIONSUBNAME_LEN+1);
336 i+=strlen(fit->voi[ri].hemisphere);
337 cptr=s+i; if(i<n) {cptr++; i++;}
338 strReplaceChar(cptr, ' ', (char)0);
339 strlcpy(fit->voi[ri].place, cptr, MAX_REGIONSUBNAME_LEN+1);
340 free(s);
341 pindx=2;
342 } else { // spaces as separators
343 strcpy(seps, " \t\n\r"); n=sn;
344 strTokenNCpy(line, seps, 1, fit->voi[ri].voiname, MAX_REGIONSUBNAME_LEN+1);
345 strTokenNCpy(line, seps, 2, fit->voi[ri].hemisphere, MAX_REGIONSUBNAME_LEN+1);
346 strTokenNCpy(line, seps, 3, fit->voi[ri].place, MAX_REGIONSUBNAME_LEN+1);
347 pindx=4;
348 }
349 rnameCatenate(fit->voi[ri].name, MAX_REGIONNAME_LEN+1,
350 fit->voi[ri].voiname, fit->voi[ri].hemisphere, fit->voi[ri].place, ' ');
351 /* Fit start and end times, and original data nr */
352 char s[128];
353 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].start=atof_dpi(s);
354 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].end=atof_dpi(s);
355 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].dataNr=atoi(s);
356 /* Fit error, parameter nr and function number (type) */
357 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].wss=atof_dpi(s);
358 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].parNr=atoi(s);
359 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].type=atoi(s);
360 /* Parameters */
361 for(i=0; i<fit->voi[ri].parNr; i++) {
362 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].p[i]=atof_dpi(s);
363 }
364 }
365 if(ri==0) {fclose(fp); fitEmpty(fit); return(12);}
366 if(ri<fit->voiNr) fit->voiNr=ri;
367
368 /* Close file */
369 fclose(fp);
370 strcpy(fiterrmsg, "");
371
372 if(verbose>1) printf("done fitRead()\n");
373 return(0);
374}
375/*****************************************************************************/
376
377/*****************************************************************************/
381int fitFunctionformat(
383 int type,
385 char *str
386) {
387 strcpy(str, "");
388 switch(type) {
389 /* Polynomials, including line */
390 case 100: strcpy(str, "f(x)=A"); break;
391 case 101: strcpy(str, "f(x)=A+B*x"); break;
392 case 102: strcpy(str, "f(x)=A+B*x+C*x^2"); break;
393 case 103: strcpy(str, "f(x)=A+B*x+C*x^2+D*x^3"); break;
394 case 104: strcpy(str, "f(x)=A+B*x+C*x^2+D*x^3+E*x^4"); break;
395 case 105: strcpy(str, "f(x)=A+B*x+C*x^2+D*x^3+E*x^4+F*x^5"); break;
396 case 106: strcpy(str, "f(x)=A+B*x+C*x^2+D*x^3+E*x^4+F*x^5+G*x^6"); break;
397 case 107: strcpy(str, "f(x)=A+B*x+C*x^2+D*x^3+E*x^4+F*x^5+G*x^6+H*x^7"); break;
398 case 108: strcpy(str, "f(x)=A+B*x+C*x^2+D*x^3+E*x^4+F*x^5+G*x^6+H*x^7+I*x^8"); break;
399 case 109: strcpy(str, "f(x)=A+B*x+C*x^2+D*x^3+E*x^4+F*x^5+G*x^6+H*x^7+I*x^8+J*x^9"); break;
400 /* Rational functions */
401 case 211: strcpy(str, "f(x)=(A+C*x)/(B+D*x)"); break;
402 case 221: strcpy(str, "f(x)=(A+C*x+E*x^2)/(B+D*x)"); break;
403 case 222: strcpy(str, "f(x)=(A+C*x+E*x^2)/(B+D*x+F*x^2)"); break;
404 case 232: strcpy(str, "f(x)=(A+C*x+E*x^2+G*x^3)/(B+D*x+F*x^2)"); break;
405 case 233: strcpy(str, "f(x)=(A+C*x+E*x^2+G*x^3)/(B+D*x+F*x^2+H*x^3)"); break;
406 case 1232: strcpy(str, "f(x)=(A+C*(x-t)+E*(x-t)^2+G*(x-t)^3)/(B+D*(x-t)+F*(x-t)^2)"); break;
407 /* Rational function for plasma-to-blood ratio */
408 case 2233: strcpy(str, "f(x)=1/(1-H*(1-r(x))), where r(x) is 3/3 order rational function for RBC-to-plasma"); break;
409 /* Exponential functions */
410 case 301: strcpy(str, "f(x)=A*exp(B*x)"); break;
411 case 302: strcpy(str, "f(x)=A*exp(B*x)+C*exp(D*x)"); break;
412 case 303: strcpy(str, "f(x)=A*exp(B*x)+C*exp(D*x)+E*exp(F*x)"); break;
413 case 304: strcpy(str, "f(x)=A*exp(B*x)+C*exp(D*x)+E*exp(F*x)+G*exp(H*x)"); break;
414 case 305: strcpy(str, "f(x)=A*exp(B*x)+C*exp(D*x)+E*exp(F*x)+G*exp(H*x)+I*exp(J*x)"); break;
415 /* Feng function */
416 case 1312: strcpy(str, "f(x)=(A*(x-t)-C)*exp(B*(x-t))+C*exp(D*(x-t))"); break;
417 case 1313: strcpy(str, "f(x)=(A*(x-t)-C-E)*exp(B*(x-t))+C*exp(D*(x-t))+E*exp(F*(x-t))"); break;
418 case 1314: strcpy(str, "f(x)=(A*(x-t)-C-E-G)*exp(B*(x-t))+C*exp(D*(x-t))+E*exp(F*(x-t))+G*exp(H*(x-t))"); break;
419 case 2313: strcpy(str, "f(x)=1/(1-A*(1-((B*x-D-F)*exp(C*x)+D*exp(E*x)+F*exp(G*x))))"); break;
420 /* Lundqvist function */
421 case 321: strcpy(str, "f(x)=A*exp(B*x)*(1-exp(C*x))"); break;
422 case 322: strcpy(str, "f(x)=A*exp(B*x)*(1-exp(C*x))+D*exp(E*x)*(1-exp(F*x))"); break;
423 case 323: strcpy(str, "f(x)=A*exp(B*x)*(1-exp(C*x))+D*exp(E*x)*(1-exp(F*x))+G*exp(H*x)*(1-exp(I*x))"); break;
424 case 1321: strcpy(str, "f(x)=A*exp(-B*(x-t))*(1-exp(-C*(x-t))) + D*(A/(B*(B+C)))*(C-((B+C)*exp(C*(x-t))-B)*exp(-(B+C)*(x-t))) "); break;
425 /* Exponential bolus infusion functions */
426 case 331: strcpy(str, "f(x)=(1/Ti)*Sum[i=1..n, (Ai/Li)*(exp(-Li*(x-tA-Ti)) - exp(-Li*(x-Ta)))], when x>=Ta+Ti\nf(x)=(1/Ti)*Sum[i=1..n, (Ai/Li)*(1-exp(-Li*(t-Ta)))], when x>Ta and x<Ta+Ti\nf(x)=0, when t<=Ta");
427 break;
428 /* Kudomi's function for radiowater */
429 case 332: strcpy(str, "f(x)=0, when x<=Ta \nf(x)=(A/L^2)*(1-exp(-L(x-Ta))), when Ta<x<=Ta+Ti \nf(x)=(A/L^2)*(exp(-L*Ti)+exp(-L(x-Ta-Ti))-2exp(-L(x-t1))), when x>Ta+Ti");
430 break;
431 /* Bolus infusion approaching zero */
432 case 334: strcpy(str, "f(x)=0, when x<=t1 \nf(x)=A*(1-exp(L(t1-x)))/(1-exp(L*(t1-t2))), when t1<x<=t2 \nf(x)=A*(exp(L(t2-x))-exp(L(t1-x)))/(1-exp(L*(t1-t2))), when x>t2");
433 break;
434 /* Exponential functions for plasma fractions */
435 case 351: strcpy(str, "f(x)=1-a*(2-exp(-b*x)-exp(-c*x))"); break;
436 /* Inverted gamma cdf for plasma parent fraction */
437 case 403: strcpy(str, "f(x)=A*(1-B*gammacdf(D,C*(x-t)))"); break;
438 /* Gamma variate function */
439 case 1401: strcpy(str, "f(x)=A*((x-D)^B)*exp(-(x-D)/C) , when x>=D, else f(x)=0"); break;
440 /* Gamma variate function with background */
441 case 1402: strcpy(str, "f(x)=A*((x-D)^B)*exp(-(x-D)/C) + E , when x>=D, else f(x)=E"); break;
442 /* Gamma variate bolus plus recirculation function */
443 case 1403: strcpy(str, "f(x)=B*((x-A)^C)*exp(-(x-A)/D) + E*(1-exp(-(x-A)/D))*exp(-(x-A)/F) , when x>A, else f(x)=0"); break;
444 /* Weibull cdf */
445 case 1421: strcpy(str, "f(x)=A*(1-exp(-((x-t)/B)^C) , when x>t, else f(x)=0"); break;
446 /* Weibull cdf plus pdf (derivative of cdf) */
447 case 1423: strcpy(str, "f(x)=A*[C*((x-t)/B)^(C-1)*exp(-((x-t)/B)^C))/B + k*(1-exp(-((x-t)/B)^C))] , when x>t, else f(x)=0"); break;
448 /* Surge function with AUC=A */
449 case 1431: strcpy(str, "f(x)=A*x*exp(-B*x)*B^2 , when x>0, else f(x)=0"); break;
450 /* Traditional Surge function */
451 case 1432: strcpy(str, "f(x)=A*x*exp(-B*x) , when x>0, else f(x)=0"); break;
452 /* Surge function with recirculation */
453 case 1433: strcpy(str, "f(x)=A*[x*exp(-B*x) + (C/B^2)*(1-(B*x+1)*exp(-B*x))], when x>0, else f(x)=0"); break;
454 /* Surge function with recirculation for plasma-to-blood ratio */
455 case 1434: strcpy(str, "f(x)=1/(1-H*(1-r(x))), where r(x) is function for RBC-to-plasma"); break;
456 /* Surge function for late FDG input */
457 case 1435: strcpy(str, "f(x)=b*(m*(x-t)*exp(-n*a*(x-t)) + exp(-a*(-t)x)) , when x>t, else f(x)=0"); break;
458 /* Probability density function of Erlang distribution */
459 case 1441: strcpy(str, "f(x)=A*(k^n)*(x^(n-1))*exp(-k*x)/(n-1)! , when x>0, else f(x)=0"); break;
460 /* Hill functions for TACs */
461 case 1801: strcpy(str, "f(x)=[A*(x-t)^B]/[(x-t)^B+C^B]"); break;
462 case 1811: strcpy(str, "f(x)=A*{[B*(x-t)^(B-1)]/[C^B+(x-t)^B] - [B*(x-t)^(2*B-1)]/[C^B+(x-t)^B]^2}"); break;
463 case 1821: strcpy(str, "f(x)=A*{[B*(x-t)^(B-1)]/[C^B+(x-t)^B] - [B*(x-t)^(2*B-1)]/[C^B+(x-t)^B]^2 + K*(x-t)^B]/[(x-t)^B+C^B}"); break;
464 /* Surge function with recirculation and delay */
465 case 1833: strcpy(str, "f(x)=A*[(x-D)*exp(-B*(x-D)) + (A*C/B^2)*(1-(B*(x-D)+1)*exp(-B*(x-D)))], when x>D, else f(x)=0"); break;
466 /* Hill functions for dose-response curves */
467 case 2801: strcpy(str, "f(x)=B+[A-B]/[1+(C/x)^D]"); break;
468 case 2802: strcpy(str, "f(x)=B+[A-B]/[1+10^{(C-x)*D}]"); break;
469 /* Exponential/power type functions for parent fractions */
470 case 831: strcpy(str, "f(x)=B*(C*exp(-D*(x-A)^3/((x-A)^2+E))+(1-C)*exp(-F*(x-A))), when x>A, else f(x)=B"); break;
471 /* Hill type functions for fractions */
472 case 841: strcpy(str, "f(x)=(A*x^B)/(x^B+C)"); break;
473 case 842: strcpy(str, "f(x)=1-((A*x^B)/(x^B+C))"); break;
474 case 843: strcpy(str, "f(x)=1-((A*(1+D*x)*x^B)/(x^B+C))"); break;
475 case 844: strcpy(str, "f(x)=(A*(x-t)^B)/((x-t)^B+C)+D, when x>t, else f(x)=D"); break;
476 case 845: strcpy(str, "f(x)=A-(A*x^B)/(x^B+C))"); break;
477 case 846: strcpy(str, "f(x)=D+((A-D)*(x-t)^B)/((x-t)^B+C), when x>t, else f(x)=D"); break;
478 case 847: strcpy(str, "f(x)=1-D-((A-D)*(x-t)^B)/((x-t)^B+C), when x>t, else f(x)=1-A"); break;
479 case 848: strcpy(str, "f(x)=D*((1-A)*(x-t)^B)/((x-t)^B+C), when x>t, else f(x)=D"); break;
480 case 849: strcpy(str, "f(x)=1-D*((1-A)*(x-t)^B)/((x-t)^B+C), when x>t, else f(x)=1-A"); break;
481 /* Hill function for plasma-to-blood ratio */
482 case 2841: strcpy(str, "f(x)=1/(1-H*(1-r(x))), where r(x) is Hill function for RBC-to-plasma"); break;
483 /* Mamede/Watabe function for fractions */
484 case 851: strcpy(str, "f(x)=1/(1+(A*x)^2)^B"); break;
485 case 852: strcpy(str, "f(x)=1-1/(1+(A*x)^2)^B"); break;
486 /* Mamede/Watabe function for fractions, as extended by Meyer */
487 case 861: strcpy(str, "f(x)=(1+(A*(x-t))^B)^(-C), when x>t, else f(x)=1"); break;
488 case 862: strcpy(str, "f(x)=1-(1+(A*(x-t))^B)^(-C), when x>t, else f(x)=0"); break;
489 /* ... and further extended by letting fraction start somewhere between 0 and 1 */
490 case 863: strcpy(str, "f(x)=(D^(-1/C)+(A*(x-t))^B)^(-C), when x>t, else f(x)=D"); break;
491 case 864: strcpy(str, "f(x)=1-(D^(-1/C)+(A*(x-t))^B)^(-C), when x>t, else f(x)=1-D"); break;
492 /* Functions for fitting plasma fractions via separate metabolite fractions */
493 case 871:
494 case 881:
495 strcpy(str, "f(x)=1-f1(x)-f2(x)-f3(x)"); break;
496 case 872:
497 case 882:
498 strcpy(str, "f1(x)=a(x)(1-b(x)-c(x)+b(x)c(x))/(1-a(x)b(x)-a(x)c(x)-b(x)c(x)+2a(x)b(x)c(x))"); break;
499 case 873:
500 case 883:
501 strcpy(str, "f2(x)=b(x)(1-a(x)-c(x)+a(x)c(x))/(1-a(x)b(x)-a(x)c(x)-b(x)c(x)+2a(x)b(x)c(x))"); break;
502 case 874:
503 case 884:
504 strcpy(str, "f3(x)=c(x)(1-a(x)-b(x)+a(x)b(x))/(1-a(x)b(x)-a(x)c(x)-b(x)c(x)+2a(x)b(x)c(x))"); break;
505 case 1010:
506 strcpy(str, "f(x)=p2 when x>=p1, f(x)=p4 when x>=p2, ..."); break;
507 break;
508 /* PET profile functions */
509 case 2111: strcpy(str, "P(x)=(C/2)*(erf((x-d+R)/(sqrt(2)*FWHM/2355))-erf((x-d-R)/(sqrt(2)*FWHM/2355)))+bkg");
510 break;
511 /* Combined functions and models */
512 case 3331: strcpy(str, "f(x)=(1/Ti)*Sum[i=1..n, (Ai/Li)*(exp(-Li*(x-tA-Ti)) - exp(-Li*(x-Ta)))], when x>=Ta+Ti\nf(x)=(1/Ti)*Sum[i=1..n, (Ai/Li)*(1-exp(-Li*(t-Ta)))], when x>Ta and x<Ta+Ti\nf(x)=0, when t<=Ta, with additional delay and dispersion");
513 break;
514 /* Compartmental model functions */
515 /* Graham's plasma curve smoothing function */
516 case 9501: strcpy(str, "Cp(t)<=>Ci(t)<=>Ct(t)"); break;
517 /* Extended Graham's plasma curve smoothing function */
518 case 9502: strcpy(str, "Ce(t)<=>Cp(t)<=>Ci(t)<=>Ct(t)"); break;
519 /* Extended Graham's plasma curve smoothing function with metabolite */
520 case 9503: strcpy(str, "Cpa(t)<=>Cia(t)<=>Cta(t)->Ctm(t)<=>Cim(t)<=>Cpm(t)"); break;
521 /* Huang's plasma metabolite model */
522 case 9601: strcpy(str, "C4(t)<=>C3(t)<-C0(t)->C1(t)<=>C2(t)"); break;
523 /* Extended Carson's plasma metabolite model */
524 case 9602: strcpy(str, "Cpa(t)<=>Cta(t)->Ctm(t)<=>Cpm(t)"); break;
525 /* New plasma metabolite model */
526 case 9603: strcpy(str, "Cpa(t)->Ct1(t)<=>Cpm(t)<=>Ct2(t)"); break;
527 /* Multi-linear multi-compartmental TAC fitting model */
528 case 9701: strcpy(str, "Ideal bolus -> n compartments"); break;
529 default: return(1);
530 }
531 return(0);
532}
533/*****************************************************************************/
534
535/*****************************************************************************/
539int fitFunctionname(
541 int type,
543 char *str
544) {
545 strcpy(str, "");
546 switch(type) {
547 case 100: strcpy(str, "f(x)=A"); break;
548 case 101: strcpy(str, "line"); break;
549 case 102: strcpy(str, "2nd order polynomial"); break;
550 case 103: strcpy(str, "3rd order polynomial"); break;
551 case 104: strcpy(str, "4th order polynomial"); break;
552 case 105: strcpy(str, "5th order polynomial"); break;
553 case 106: strcpy(str, "6th order polynomial"); break;
554 case 107: strcpy(str, "7th order polynomial"); break;
555 case 108: strcpy(str, "8th order polynomial"); break;
556 case 109: strcpy(str, "9th order polynomial"); break;
557 case 211: strcpy(str, "1/1 order rational function"); break;
558 case 221: strcpy(str, "2/1 order rational function"); break;
559 case 222: strcpy(str, "2/2 order rational function"); break;
560 case 232: strcpy(str, "3/2 order rational function"); break;
561 case 233: strcpy(str, "3/3 order rational function"); break;
562 case 1232: strcpy(str, "3/2 order rational function with delay"); break;
563 case 2233: strcpy(str, "3/3 order rational function for plasma-to-blood ratio"); break;
564 case 301: strcpy(str, "exponential function"); break;
565 case 302: strcpy(str, "sum of 2 exponential functions"); break;
566 case 303: strcpy(str, "sum of 3 exponential functions"); break;
567 case 304: strcpy(str, "sum of 4 exponential functions"); break;
568 case 305: strcpy(str, "sum of 5 exponential functions"); break;
569 case 1312: strcpy(str, "Feng model 2 function with 2 exponentials"); break;
570 case 1313: strcpy(str, "Feng model 2 function"); break;
571 case 1314: strcpy(str, "Feng model 2 function with 4 exponentials"); break;
572 case 2313: strcpy(str, "Feng M2 function for plasma-to-blood ratio"); break;
573 case 321: strcpy(str, "Lundqvist function"); break;
574 case 322: strcpy(str, "sum of 2 Lundqvist functions"); break;
575 case 323: strcpy(str, "sum of 3 Lundqvist functions"); break;
576 case 1321: strcpy(str, "Lundqvist function with integral and delay"); break;
577 case 331: strcpy(str, "Exponential bolus infusion function"); break;
578 case 332: strcpy(str, "Kudomi's exponential bolus infusion function for radiowater"); break;
579 case 334: strcpy(str, "Exponential bolus function approaching zero"); break;
580 case 351: strcpy(str, "Exponential function for [C-11]PK11195 plasma fractions"); break;
581 case 403: strcpy(str, "Inverted gamma cdf for plasma parent fraction"); break;
582 case 1401: strcpy(str, "Gamma variate function"); break;
583 case 1402: strcpy(str, "Gamma variate with background"); break;
584 case 1403: strcpy(str, "Gamma variate bolus plus recirculation"); break;
585 case 1421: strcpy(str, "Weibull cdf with delay"); break;
586 case 1423: strcpy(str, "Weibull cdf and derivative with delay"); break;
587 case 1431: strcpy(str, "Surge function"); break;
588 case 1432: strcpy(str, "Surge function (trad)"); break;
589 case 1433: strcpy(str, "Surge function with recirculation"); break;
590 case 1434: strcpy(str, "Surge function with recirculation for plasma-to-blood ratio"); break;
591 case 1435: strcpy(str, "Surge function for late FDG AIF"); break;
592 case 1441: strcpy(str, "Probability density function of Erlang distribution"); break;
593 case 1801: strcpy(str, "Hill function with delay"); break;
594 case 1811: strcpy(str, "Derivative of Hill function with delay"); break;
595 case 1821: strcpy(str, "Sum of Hill function and derivative with delay"); break;
596 case 1833: strcpy(str, "Surge function with recirculation and delay"); break;
597 case 2801: strcpy(str, "Hill function for dose-response curve on linear scale"); break;
598 case 2802: strcpy(str, "Hill function for dose-response curve on log scale"); break;
599 case 831: strcpy(str, "Exponential/Power function for parent fractions"); break;
600 case 841: strcpy(str, "Hill function"); break;
601 case 842: strcpy(str, "Hill function (1-f(x))"); break;
602 case 843: strcpy(str, "Hill function (1-f(x)) with ascending or descending end"); break;
603 case 844: strcpy(str, "Hill function with background"); break;
604 case 845: strcpy(str, "Hill function (A-f(x))"); break;
605 case 846: strcpy(str, "Extended Hill function for plasma parent fraction"); break;
606 case 847: strcpy(str, "Extended Hill function for plasma metabolite fraction"); break;
607 case 848: strcpy(str, "Extended Hill function #2 for plasma parent fraction"); break;
608 case 849: strcpy(str, "Extended Hill function #2 for plasma metabolite fraction"); break;
609 case 851: strcpy(str, "Mamede function"); break;
610 case 852: strcpy(str, "Mamede function (1-f(x)"); break;
611 case 861: strcpy(str, "Meyer parent fraction function"); break;
612 case 862: strcpy(str, "Meyer metabolite fraction function"); break;
613 case 863: strcpy(str, "Extended Meyer parent fraction function"); break;
614 case 864: strcpy(str, "Extended Meyer metabolite fraction function"); break;
615 case 871: strcpy(str, "1-3 metabolite Hill function for parent"); break;
616 case 872: strcpy(str, "1-3 metabolite Hill function for metab1"); break;
617 case 873: strcpy(str, "1-3 metabolite Hill function for metab2"); break;
618 case 874: strcpy(str, "1-3 metabolite Hill function for metab3"); break;
619 case 881: strcpy(str, "1-3 metabolite power function for parent"); break;
620 case 882: strcpy(str, "1-3 metabolite power function for metab1"); break;
621 case 883: strcpy(str, "1-3 metabolite power function for metab2"); break;
622 case 884: strcpy(str, "1-3 metabolite power function for metab3"); break;
623 case 1010: strcpy(str, "Step function"); break;
624 case 2111: strcpy(str, "Image profile function"); break;
625 case 2841: strcpy(str, "Hill function for plasma-to-blood ratio"); break;
626 case 3331: strcpy(str, "Exponential bolus infusion function with delay and dispersion"); break;
627 case 9501: strcpy(str, "Graham's input function"); break;
628 case 9502: strcpy(str, "Extended Graham's input function"); break;
629 case 9503: strcpy(str, "Graham's input function with metabolite"); break;
630 case 9601: strcpy(str, "Huang's plasma metabolite model"); break;
631 case 9602: strcpy(str, "Extended Carson's plasma metabolite model"); break;
632 case 9603: strcpy(str, "New plasma metabolite model"); break;
633 case 9701: strcpy(str, "Multilinear multicompartmental TAC fitting model"); break;
634 default: return(1);
635 }
636 return(0);
637}
638/*****************************************************************************/
639
640/*****************************************************************************/
645int fitEval(
647 FitVOI *r,
649 double x,
651 double *y
652) {
653 double a, b, f, sqrx, cubx, xt;
654 int i, j, m, n;
655 double mp[3][5];
656 double mf[]={0.0,0.0,0.0};
657
658 if(r==NULL) return(1);
659 if(MATHFUNC_TEST>0) {
660 printf("fitEval(r, %g, %g): type=%d ;", x, *y, r->type);
661 for(i=0; i<r->parNr; i++) printf(" %g", r->p[i]);
662 printf("\n");
663 }
664 if(y==NULL) return(1);
665 *y=nan("");
666 switch(r->type) {
667 case 100:
668 case 101:
669 case 102:
670 case 103:
671 case 104:
672 case 105:
673 case 106:
674 case 107:
675 case 108:
676 case 109:
677 n=r->type-99; f=0.0; for(i=n-1; i>0; i--) {f+=r->p[i]; f*=x;} f+=r->p[0];
678 *y=f;
679 break;
680 case 211:
681 a=r->p[0]+r->p[2]*x; b=r->p[1]+r->p[3]*x; if(b==0.0) break;
682 f=a/b; *y=f;
683 break;
684 case 221:
685 sqrx=x*x;
686 a=r->p[0]+r->p[2]*x+r->p[4]*sqrx; b=r->p[1]+r->p[3]*x; if(b==0.0) break;
687 f=a/b; *y=f;
688 break;
689 case 222:
690 sqrx=x*x;
691 a=r->p[0]+r->p[2]*x+r->p[4]*sqrx;
692 b=r->p[1]+r->p[3]*x+r->p[5]*sqrx; if(b==0.0) break;
693 f=a/b; *y=f;
694 break;
695 case 232:
696 sqrx=x*x; cubx=sqrx*x;
697 a=r->p[0]+r->p[2]*x+r->p[4]*sqrx+r->p[6]*cubx;
698 b=r->p[1]+r->p[3]*x+r->p[5]*sqrx; if(b==0.0) break;
699 f=a/b; *y=f;
700 break;
701 case 233:
702 sqrx=x*x; cubx=sqrx*x;
703 a=r->p[0]+r->p[2]*x+r->p[4]*sqrx+r->p[6]*cubx;
704 b=r->p[1]+r->p[3]*x+r->p[5]*sqrx+r->p[7]*cubx; if(b==0.0) break;
705 f=a/b; *y=f;
706 break;
707 case 301:
708 f=r->p[0]*exp(r->p[1]*x); *y=f;
709 break;
710 case 302:
711 f=0;
712 a=r->p[0]*exp(r->p[1]*x); f+=a;
713 a=r->p[2]*exp(r->p[3]*x); f+=a;
714 *y=f;
715 break;
716 case 303:
717 f=0;
718 a=r->p[0]*exp(r->p[1]*x); f+=a;
719 a=r->p[2]*exp(r->p[3]*x); f+=a;
720 a=r->p[4]*exp(r->p[5]*x); f+=a;
721 *y=f;
722 break;
723 case 304:
724 f=0;
725 a=r->p[0]*exp(r->p[1]*x); f+=a;
726 a=r->p[2]*exp(r->p[3]*x); f+=a;
727 a=r->p[4]*exp(r->p[5]*x); f+=a;
728 a=r->p[6]*exp(r->p[7]*x); f+=a;
729 *y=f;
730 break;
731 case 305:
732 f=0;
733 a=r->p[0]*exp(r->p[1]*x); f+=a;
734 a=r->p[2]*exp(r->p[3]*x); f+=a;
735 a=r->p[4]*exp(r->p[5]*x); f+=a;
736 a=r->p[6]*exp(r->p[7]*x); f+=a;
737 a=r->p[8]*exp(r->p[9]*x); f+=a;
738 *y=f;
739 break;
740 case 321:
741 f=r->p[0]*exp(r->p[1]*x)*(1.0-exp(r->p[2]*x));
742 *y=f;
743 break;
744 case 322:
745 f=0.0;
746 a=r->p[0]*exp(r->p[1]*x)*(1.0-exp(r->p[2]*x)); f+=a;
747 a=r->p[3]*exp(r->p[4]*x)*(1.0-exp(r->p[5]*x)); f+=a;
748 *y=f;
749 break;
750 case 323:
751 f=0.0;
752 a=r->p[0]*exp(r->p[1]*x)*(1.0-exp(r->p[2]*x)); f+=a;
753 a=r->p[3]*exp(r->p[4]*x)*(1.0-exp(r->p[5]*x)); f+=a;
754 a=r->p[6]*exp(r->p[7]*x)*(1.0-exp(r->p[8]*x)); f+=a;
755 *y=f;
756 break;
757 case 1321:
758 /* A=p0 B=p1 C=p2 D=k=p3 dT=p4 */
759 xt=x-r->p[4]; if(xt<=0.0) {
760 f=0.0;
761 } else {
762 a=exp(-r->p[2]*xt); f=r->p[0]*exp(-r->p[1]*xt)*(1.0-a);
763 if(r->p[3]>0.0)
764 f+=r->p[3]*(r->p[0]/(r->p[1]*(r->p[1]+r->p[2])))
765 *(r->p[2]-((r->p[1]+r->p[2])/a-r->p[1])*exp(-(r->p[1]+r->p[2])*xt));
766 }
767 *y=f;
768 break;
769 case 331:
770 n=(r->parNr-2)/2;
771 if(x<=r->p[0])
772 *y=0.0;
773 else if(x<r->p[0]+r->p[1]) {
774 for(i=0, f=0.0; i<n; i++) {
775 if(r->p[2*i+3]>1.0E-12) b=r->p[2*i+2]/r->p[2*i+3]; else b=r->p[2*i+2];
776 f+=b*(1.0-exp(-r->p[2*i+3]*(x-r->p[0])));
777 }
778 if(r->p[1]>0.0) f*=(1.0/r->p[1]);
779 *y=f;
780 } else {
781 for(i=0, f=0.0; i<n; i++) {
782 if(r->p[2*i+3]>1.0E-12) b=r->p[2*i+2]/r->p[2*i+3]; else b=r->p[2*i+2];
783 f+=b*(exp(-r->p[2*i+3]*(x-r->p[0]-r->p[1])) - exp(-r->p[2*i+3]*(x-r->p[0])));
784 }
785 if(r->p[1]>0.0) f*=(1.0/r->p[1]);
786 *y=f;
787 }
788 break;
789 case 332:
790 if(x<=r->p[0])
791 *y=0.0;
792 else if(x<=r->p[0]+r->p[1]) {
793 f=exp(-r->p[3]*(x-r->p[0]));
794 *y=(r->p[2]/(r->p[3]*r->p[3]))*(1.0-f);
795 } else {
796 f=exp(-r->p[3]*r->p[1]);
797 f+=exp(-r->p[3]*(x-r->p[0]-r->p[1]));
798 f-=2.0*exp(-r->p[3]*(x-r->p[0]));
799 *y=(r->p[2]/(r->p[3]*r->p[3]))*f;
800 }
801 break;
802 case 334:
803 if(x<=r->p[0])
804 *y=0.0;
805 else if(x>r->p[0] && x<=r->p[1]) {
806 *y=r->p[2]*(1.0-exp(r->p[3]*(r->p[0]-x)))/(1.0-exp(r->p[3]*(r->p[0]-r->p[1])));
807 } else {
808 *y=r->p[2]*(exp(r->p[3]*(r->p[1]-x))-exp(r->p[3]*(r->p[0]-x)))/(1.0-exp(r->p[3]*(r->p[0]-r->p[1])));
809 }
810 break;
811 case 351:
812 f=1.0-r->p[0]*(2.0-exp(-r->p[1]*x)-exp(-r->p[2]*x)); *y=f;
813 break;
814 case 403:
815 xt=x-r->p[4]; if(xt<0.0) xt=0.0;
816 f=r->p[0]*(1.0 - r->p[1]*igam(r->p[3], r->p[2]*xt));
817 *y=f;
818 break;
819 case 831:
820 xt=x-r->p[0]; if(xt<=0.0) {
821 f=r->p[1];
822 } else {
823 f=r->p[2]*exp(-r->p[3]*xt*xt*xt/(xt*xt+r->p[4])) + (1.0-r->p[2])*exp(-r->p[5]*xt);
824 f*=r->p[1];
825 }
826 *y=f;
827 break;
828 case 841:
829 f=r->p[0]*pow(x, r->p[1]) / (pow(x, r->p[1]) + r->p[2]);
830 *y=f;
831 break;
832 case 842:
833 f=1.0-r->p[0]*pow(x, r->p[1]) / (pow(x, r->p[1]) + r->p[2]);
834 *y=f;
835 break;
836 case 843:
837 f=1.0 - r->p[0]*(1.0+r->p[3]*x)*pow(x, r->p[1]) / (pow(x, r->p[1]) + r->p[2]);
838 *y=f;
839 break;
840 case 844:
841 if(r->parNr>4) xt=x-r->p[4]; else xt=x; // delay time accounted, if avail
842 if(xt<=0.0) f=r->p[3];
843 else {a=pow(xt, r->p[1]); f=r->p[0]*a/(a + r->p[2]) + r->p[3];}
844 *y=f;
845 break;
846 case 845:
847 f=r->p[0] - r->p[0]*pow(x, r->p[1]) / (pow(x, r->p[1]) + r->p[2]);
848 *y=f;
849 break;
850 case 846:
851 xt=x-r->p[4]; if(xt<=0.0) {
852 f=r->p[3];
853 } else {
854 a=pow(xt, r->p[1]);
855 f=r->p[3]+(r->p[0]-r->p[3])*a/(a+r->p[2]);
856 }
857 *y=f;
858 break;
859 case 847:
860 xt=x-r->p[4]; if(xt<=0.0) {
861 f=1.0-r->p[3];
862 } else {
863 a=pow(xt, r->p[1]);
864 f=1.0-r->p[3]-(r->p[0]-r->p[3])*a/(a+r->p[2]);
865 }
866 *y=f;
867 break;
868 case 848:
869 xt=x-r->p[4]; if(xt<=0.0) {
870 f=r->p[3];
871 } else {
872 a=pow(xt, r->p[1]);
873 f=r->p[3] * (1.0 - (1.0-r->p[0])*a/(r->p[2]+a));
874 }
875 *y=f;
876 break;
877 case 849:
878 xt=x-r->p[4]; if(xt<=0.0) {
879 f=1.0-r->p[3];
880 } else {
881 a=pow(xt, r->p[1]);
882 f=1.0 - r->p[3] * (1.0 - (1.0-r->p[0])*a/(r->p[2]+a));
883 }
884 *y=f;
885 break;
886 case 851:
887 a=r->p[0]*x; a*=a; f=1.0/pow(1.0+a, r->p[1]);
888 *y=f;
889 break;
890 case 852:
891 a=r->p[0]*x; a*=a; f=1.0 - 1.0/pow(1.0+a, r->p[1]);
892 *y=f;
893 break;
894 case 861:
895 xt=x-r->p[3]; if(xt<=0.0) {
896 f=1.0;
897 } else {
898 f=pow(1.0+pow(r->p[0]*xt, r->p[1]), -r->p[2]);
899 }
900 *y=f;
901 break;
902 case 862:
903 xt=x-r->p[3]; if(xt<=0.0) {
904 f=1.0;
905 } else {
906 f=pow(1.0+pow(r->p[0]*xt, r->p[1]), -r->p[2]);
907 }
908 *y=1.0-f;
909 break;
910 case 863:
911 xt=x-r->p[4]; if(xt<=0.0) {
912 f=r->p[3];
913 } else {
914 f=pow(pow(r->p[3],-1.0/r->p[2])+pow(r->p[0]*xt, r->p[1]), -r->p[2]);
915 }
916 *y=f;
917 break;
918 case 864:
919 xt=x-r->p[4]; if(xt<=0.0) {
920 f=1.0-r->p[3];
921 } else {
922 f=1.0-pow(pow(r->p[3],-1.0/r->p[2])+pow(r->p[0]*xt, r->p[1]), -r->p[2]);
923 }
924 *y=f;
925 break;
926 case 871:
927 case 872:
928 case 873:
929 case 874:
930 for(m=0; m<3; m++) for(j=0; j<5; j++) mp[m][j]=0.0;
931 for(i=m=j=0; i<r->parNr; i++) {mp[m][j]=r->p[i]; j++; if(j>4) {j=0; m++;}}
932 n=r->parNr/5;
933 for(m=0; m<n; m++) {
934 xt=x-mp[m][4];
935 if(xt<=0.0) mf[m]=1.0-mp[m][3];
936 else {
937 f=pow(xt, mp[m][1]);
938 mf[m]=1.0-(mp[m][3]+(mp[m][0]-mp[m][3])*f/(mp[m][2]+f));
939 }
940 }
941 a=1.0-mf[0]*mf[1]-mf[0]*mf[2]-mf[1]*mf[2]+2.0*mf[0]*mf[1]*mf[2];
942 if(r->type==871) {
943 f=1.0;
944 f-=mf[0]*(1.0-mf[1]-mf[2]+mf[1]*mf[2])/a;
945 f-=mf[1]*(1.0-mf[0]-mf[2]+mf[0]*mf[2])/a;
946 f-=mf[2]*(1.0-mf[0]-mf[1]+mf[0]*mf[1])/a;
947 *y=f;
948 } else if(r->type==872) {
949 *y=mf[0]*(1.0-mf[1]-mf[2]+mf[1]*mf[2])/a;
950 } else if(r->type==873) {
951 *y=mf[1]*(1.0-mf[0]-mf[2]+mf[0]*mf[2])/a;
952 } else if(r->type==874) {
953 *y=mf[2]*(1.0-mf[0]-mf[1]+mf[0]*mf[1])/a;
954 }
955 break;
956 case 881:
957 case 882:
958 case 883:
959 case 884:
960 for(m=0; m<3; m++) for(j=0; j<5; j++) mp[m][j]=0.0;
961 for(i=m=j=0; i<r->parNr; i++) {mp[m][j]=r->p[i]; j++; if(j>4) {j=0; m++;}}
962 n=r->parNr/5;
963 for(m=0; m<n; m++) {
964 xt=x-mp[m][4];
965 if(xt<=0.0) mf[m]=1.0-mp[m][3];
966 else {
967 mf[m]=1.0 - pow(pow(mp[m][3], -1.0/mp[m][2]) + pow(mp[m][0]*(xt), mp[m][1]), -mp[m][2]);
968 }
969 }
970 a=1.0-mf[0]*mf[1]-mf[0]*mf[2]-mf[1]*mf[2]+2.0*mf[0]*mf[1]*mf[2];
971 if(r->type==881) {
972 f=1.0;
973 f-=mf[0]*(1.0-mf[1]-mf[2]+mf[1]*mf[2])/a;
974 f-=mf[1]*(1.0-mf[0]-mf[2]+mf[0]*mf[2])/a;
975 f-=mf[2]*(1.0-mf[0]-mf[1]+mf[0]*mf[1])/a;
976 *y=f;
977 } else if(r->type==882) {
978 *y=mf[0]*(1.0-mf[1]-mf[2]+mf[1]*mf[2])/a;
979 } else if(r->type==883) {
980 *y=mf[1]*(1.0-mf[0]-mf[2]+mf[0]*mf[2])/a;
981 } else if(r->type==884) {
982 *y=mf[2]*(1.0-mf[0]-mf[1]+mf[0]*mf[1])/a;
983 }
984 break;
985 case 1010: // step functions
986 /* Parameter number must be an even number and >=2 */
987 if(r->parNr<2 || r->parNr&1) return(2);
988 *y=0.0;
989 for(int i=0; i<r->parNr; i+=2) if(x>=r->p[i]) *y=r->p[i+1];
990 break;
991 case 2111:
992 xt=x-r->p[3]; a=sqrt(2.0)*(r->p[2]/2.355);
993 f=r->p[4]+(r->p[0]/2.0)*(erf((xt+r->p[1])/a)-erf((xt-r->p[1])/a));
994 *y=f;
995 break;
996 case 1232:
997 xt=x; if(r->parNr>7) xt-=r->p[7];
998 if(xt<=0.0) {
999 f=0.0;
1000 } else {
1001 sqrx=xt*xt; cubx=sqrx*xt;
1002 a=r->p[0]+r->p[2]*xt+r->p[4]*sqrx+r->p[6]*cubx;
1003 b=r->p[1]+r->p[3]*xt+r->p[5]*sqrx; if(b==0.0) break;
1004 f=a/b;
1005 }
1006 *y=f;
1007 break;
1008 case 1312:
1009 xt=x; if(r->parNr>4) xt-=r->p[4];
1010 if(xt<=0.0) {
1011 f=0.0;
1012 } else {
1013 f=0.0;
1014 a=(r->p[0]*(xt)-r->p[2])*exp(r->p[1]*(xt)); f+=a;
1015 a=r->p[2]*exp(r->p[3]*(xt)); f+=a;
1016 }
1017 *y=f;
1018 break;
1019 case 1313:
1020 xt=x; if(r->parNr>6) xt-=r->p[6];
1021 if(xt<=0.0) {
1022 f=0.0;
1023 } else {
1024 f=0.0;
1025 a=(r->p[0]*(xt)-r->p[2]-r->p[4])*exp(r->p[1]*(xt)); f+=a;
1026 a=r->p[2]*exp(r->p[3]*(xt)); f+=a;
1027 a=r->p[4]*exp(r->p[5]*(xt)); f+=a;
1028 }
1029 *y=f;
1030 break;
1031 case 1314:
1032 xt=x; if(r->parNr>8) xt-=r->p[8];
1033 if(xt<=0.0) {
1034 f=0.0;
1035 } else {
1036 f=0.0;
1037 a=(r->p[0]*(xt)-r->p[2]-r->p[4]-r->p[6])*exp(r->p[1]*(xt)); f+=a;
1038 a=r->p[2]*exp(r->p[3]*(xt)); f+=a;
1039 a=r->p[4]*exp(r->p[5]*(xt)); f+=a;
1040 a=r->p[6]*exp(r->p[7]*(xt)); f+=a;
1041 }
1042 *y=f;
1043 break;
1044 case 1401:
1045 xt=x-r->p[3]; if(xt<=0.0 || r->p[2]==0.0) {
1046 f=0.0;
1047 } else {
1048 f=r->p[0]*pow(xt, r->p[1])*exp(-xt/r->p[2]);
1049 }
1050 *y=f;
1051 break;
1052 case 1402:
1053 xt=x-r->p[3]; if(xt<=0.0 || r->p[2]==0.0) {
1054 f=r->p[4];
1055 } else {
1056 f=r->p[0]*pow(xt, r->p[1])*exp(-xt/r->p[2]) + r->p[4];
1057 }
1058 *y=f;
1059 break;
1060 case 1403:
1061 xt=x-r->p[0];
1062 f=0.0;
1063 if(xt>0.0) {
1064 a=exp(-xt/r->p[3]);
1065 if(r->p[1]>0.0) f += r->p[1]*pow(xt, r->p[2])*a;
1066 if(r->parNr==6 && r->p[4]>0.0) f += r->p[4]*(1.0-a)*exp(-xt/r->p[5]);
1067 }
1068 *y=f;
1069 break;
1070 case 1421:
1071 xt=x-r->p[0]; if(xt<=0.0) {
1072 f=0.0;
1073 } else {
1074 f=r->p[1]*(1.0-exp(-pow(xt/r->p[2], r->p[3])));
1075 }
1076 *y=f;
1077 break;
1078 case 1423:
1079 xt=x-r->p[0]; if(xt<=0.0) {
1080 f=0.0;
1081 } else {
1082 a=xt/r->p[2]; b=pow(a, r->p[3]-1.0); f=exp(-b*a);
1083 a=r->p[3]*b*f/r->p[2]; b=1.0-f;
1084 f=r->p[1]*(a+r->p[4]*b);
1085 }
1086 *y=f;
1087 break;
1088 case 1431:
1089 if(x<=0.0) {
1090 f=0.0;
1091 } else {
1092 f=r->p[0]*x*exp(-r->p[1]*x)*r->p[1]*r->p[1];
1093 }
1094 *y=f;
1095 break;
1096 case 1432:
1097 if(x<=0.0) {
1098 f=0.0;
1099 } else {
1100 f=r->p[0]*x*exp(-r->p[1]*x);
1101 }
1102 *y=f;
1103 break;
1104 case 1433:
1105 if(x<=0.0) {
1106 f=0.0;
1107 } else {
1108 double e=exp(-r->p[1]*x);
1109 f=x*e + (r->p[2]/(r->p[1]*r->p[1]))*(1.0-(r->p[1]*x+1.0)*e);
1110 f*=r->p[0];
1111 }
1112 *y=f;
1113 break;
1114 case 1434:
1115 if(x<=0.0) {
1116 f=1.0/(1.0-r->p[0]);
1117 } else {
1118 double e=exp(-r->p[2]*x);
1119 double rcp=x*e + (r->p[3]/(r->p[2]*r->p[2]))*(1.0-(r->p[2]*x+1.0)*e);
1120 rcp*=r->p[1];
1121 f=1.0/(1.0-r->p[0]*(1.0-rcp));
1122 }
1123 *y=f;
1124 break;
1125 case 1435:
1126 xt=x-r->p[0]; if(xt<=0.0) {
1127 f=0.0;
1128 } else {
1129 f=r->p[2]*(r->p[3]*xt*exp(-r->p[4]*r->p[1]*xt) + exp(-r->p[1]*xt) );
1130 }
1131 *y=f;
1132 break;
1133 case 1441: /* Probability density function of Erlang distribution */
1134 f=0.0;
1135 if(x>=0.0) {
1136 double A=r->p[0], k=r->p[1]; if(!(k>=0.0)) return(2);
1137 int N=(int)round(r->p[2]); if(!(N>0)) return(2);
1138 unsigned long long int f=lfactorial(N-1);
1139 f = A*pow(k, N)*pow(x, N-1)*exp(-k*x)/f;
1140 }
1141 *y=f;
1142 break;
1143 case 1801:
1144 xt=x-r->p[0]; if(xt<=0.0) {
1145 f=0.0;
1146 } else {
1147 f=r->p[1]*pow(xt, r->p[3])/(pow(r->p[2], r->p[3])+pow(xt, r->p[3]));
1148 }
1149 *y=f;
1150 break;
1151 case 1811:
1152 xt=x-r->p[0]; if(xt<=0.0) {
1153 f=0.0;
1154 } else {
1155 a=pow(r->p[2], r->p[3])+pow(xt, r->p[3]);
1156 f=r->p[1]*r->p[3]*
1157 (pow(xt, r->p[3]-1.0)/a - pow(xt, 2.0*r->p[3]-1.0)/(a*a));
1158 }
1159 *y=f;
1160 break;
1161 case 1821:
1162 xt=x-r->p[0]; if(xt<=0.0) {
1163 f=0.0;
1164 } else {
1165 double a, c, b, k, xb, cb, cbxb, hill, hilld;
1166 a=r->p[1]; c=r->p[2]; b=r->p[3]; k=r->p[4];
1167 cb=pow(c, b); xb=pow(xt, b); cbxb=cb+xb;
1168 hilld= b*pow(xt, b-1.0)/cbxb - b*pow(xt, 2.0*b-1.0)/(cbxb*cbxb);
1169 hill=k*xb/cbxb;
1170 f=a*(hilld+hill);
1171 }
1172 *y=f;
1173 break;
1174 case 1833:
1175 xt=x; if(r->parNr>3) xt-=r->p[3];
1176 if(xt<=0.0) {
1177 f=0.0;
1178 } else {
1179 double c=0.0; if(r->parNr>2) c=r->p[2];
1180 double e=exp(-r->p[1]*xt);
1181 f=xt*e + (c/(r->p[1]*r->p[1]))*(1.0-(r->p[1]*xt+1.0)*e);
1182 f*=r->p[0];
1183 }
1184 *y=f;
1185 break;
1186 case 2233:
1187 if(x<=0.0) {
1188 f=1.0/(1.0-r->p[0]);
1189 } else {
1190 sqrx=x*x; cubx=sqrx*x;
1191 a=0.0+r->p[1]*x+r->p[3]*sqrx+r->p[5]*cubx;
1192 b=1.0+r->p[2]*x+r->p[4]*sqrx+r->p[6]*cubx;
1193 f=1.0/(1.0-r->p[0]*(1.0-(a/b)));
1194 }
1195 *y=f;
1196 break;
1197 case 2313:
1198 if(r->parNr<3) return(2);
1199 if(x<=0.0) {
1200 f=0.0;
1201 } else {
1202 double A1=r->p[1];
1203 double L1=r->p[2];
1204 double A2=0.0; if(r->parNr>3) A2=r->p[3];
1205 double L2=0.0; if(r->parNr>4) L2=r->p[4];
1206 double A3=0.0; if(r->parNr>5) A3=r->p[5];
1207 double L3=0.0; if(r->parNr>6) L3=r->p[6];
1208 f=(A1*x - A2 - A3)*exp(-L1*x) + A2*exp(-L2*x) + A3*exp(-L3*x);
1209 }
1210 *y=1.0/(1.0-r->p[0]*(1.0-f));
1211 break;
1212 case 2801:
1213 {
1214 double top, bottom, ec50, hillslope;
1215 top=r->p[0]; bottom=r->p[1]; ec50=r->p[2]; hillslope=r->p[3];
1216 if(x<=0.0) f=bottom;
1217 else f=bottom+(top-bottom)/(1.0+pow(ec50/x,hillslope));
1218 }
1219 *y=f;
1220 break;
1221 case 2802:
1222 {
1223 double top, bottom, logec50, hillslope;
1224 top=r->p[0]; bottom=r->p[1]; logec50=r->p[2]; hillslope=r->p[3];
1225 f=bottom+(top-bottom)/(1.0+pow(10.0, (logec50-x)*hillslope));
1226 }
1227 *y=f;
1228 break;
1229 case 2841:
1230 a=r->p[1]*pow(x, r->p[2]) / (pow(x, r->p[2]) + r->p[3]);
1231 f=1.0/(1.0-r->p[0]*(1.0-a));
1232 *y=f;
1233 break;
1234 case 3331: *y=nan(""); break; /* cannot be applied to one point only */
1235 case 9501: *y=nan(""); break; /* cannot be applied to one point only */
1236 case 9502: *y=nan(""); break; /* cannot be applied to one point only */
1237 case 9503: *y=nan(""); break; /* cannot be applied to one point only */
1238 case 9701: *y=nan(""); break; /* cannot be applied to one point only */
1239 default:
1240 *y=nan("");
1241 }
1242 if(isnan(*y)) return 1;
1243 return 0;
1244}
1245/*****************************************************************************/
1246
1247/*****************************************************************************/
1252int fitEvaltac(
1254 FitVOI *r,
1256 double *x,
1258 double *y,
1260 int dataNr
1261) {
1262 if(r==NULL || x==NULL || y==NULL || dataNr<1) return(1);
1263
1264 /* Special cases that can only be computed as TACs */
1265 if(r->type==3331) {
1266 double Ta=r->p[0];
1267 double Ti=r->p[1];
1268 double dT=r->p[r->parNr-2]; Ta+=dT;
1269 double tau=r->p[r->parNr-1];
1270 int n=(r->parNr-4)/2; if(n<1) return(2);
1271 for(int i=0; i<dataNr; i++) {
1272 if(x[i]<=Ta) y[i]=0.0;
1273 else if(x[i]<Ta+Ti) {
1274 double f=0.0;
1275 for(int j=0; j<n; j++) {
1276 double b;
1277 if(r->p[2*j+3]>1.0E-12) b=r->p[2*j+2]/r->p[2*j+3]; else b=r->p[2*j+2];
1278 f+=b*(1.0-exp(-r->p[2*j+3]*(x[i]-Ta)));
1279 }
1280 if(Ti>0.0) f*=(1.0/Ti);
1281 y[i]=f;
1282 } else {
1283 double f=0.0;
1284 for(int j=0; j<n; j++) {
1285 double b;
1286 if(r->p[2*j+3]>1.0E-12) b=r->p[2*j+2]/r->p[2*j+3]; else b=r->p[2*j+2];
1287 f+=b*(exp(-r->p[2*j+3]*(x[i]-Ta-Ti)) - exp(-r->p[2*j+3]*(x[i]-Ta)));
1288 }
1289 if(Ti>0.0) f*=(1.0/Ti);
1290 y[i]=f;
1291 }
1292 }
1293 if(simDispersion(x, y, dataNr, tau, 0.0, NULL)) return(2);
1294 return(0);
1295 }
1296
1297 /* Usual functions */
1298 for(int i=0; i<dataNr; i++) if(fitEval(r, x[i], y+i)) return(2);
1299 return 0;
1300}
1301/*****************************************************************************/
1302
1303/*****************************************************************************/
1308int fitIntegralEval(
1310 FitVOI *r,
1312 double x,
1314 double *yi
1315) {
1316 double a, f, xt, t;
1317
1318 if(r==NULL || yi==NULL) return 1;
1319 *yi=nan("");
1320 switch(r->type) {
1321 case 301:
1322 if(fabs(r->p[1])>1.0e-12) f=(r->p[0]/r->p[1])*(exp(r->p[1]*x)-1.0);
1323 else f=r->p[0]*x;
1324 *yi=f;
1325 break;
1326 case 302:
1327 f=0;
1328 if(fabs(r->p[1])>1.0e-12) a=(r->p[0]/r->p[1])*(exp(r->p[1]*x)-1.0);
1329 else a=r->p[0]*x;
1330 f+=a;
1331 if(fabs(r->p[3])>1.0e-12) a=(r->p[2]/r->p[3])*(exp(r->p[3]*x)-1.0);
1332 else a=r->p[2]*x;
1333 f+=a;
1334 *yi=f;
1335 break;
1336 case 303:
1337 f=0;
1338 if(fabs(r->p[1])>1.0e-12) a=(r->p[0]/r->p[1])*(exp(r->p[1]*x)-1.0);
1339 else a=r->p[0]*x;
1340 f+=a;
1341 if(fabs(r->p[3])>1.0e-12) a=(r->p[2]/r->p[3])*(exp(r->p[3]*x)-1.0);
1342 else a=r->p[2]*x;
1343 f+=a;
1344 if(fabs(r->p[5])>1.0e-12) a=(r->p[4]/r->p[5])*(exp(r->p[5]*x)-1.0);
1345 else a=r->p[4]*x;
1346 f+=a;
1347 *yi=f;
1348 break;
1349 case 304:
1350 f=0;
1351 if(fabs(r->p[1])>1.0e-12) a=(r->p[0]/r->p[1])*(exp(r->p[1]*x)-1.0);
1352 else a=r->p[0]*x;
1353 f+=a;
1354 if(fabs(r->p[3])>1.0e-12) a=(r->p[2]/r->p[3])*(exp(r->p[3]*x)-1.0);
1355 else a=r->p[2]*x;
1356 f+=a;
1357 if(fabs(r->p[5])>1.0e-12) a=(r->p[4]/r->p[5])*(exp(r->p[5]*x)-1.0);
1358 else a=r->p[4]*x;
1359 f+=a;
1360 if(fabs(r->p[7])>1.0e-12) a=(r->p[6]/r->p[7])*(exp(r->p[7]*x)-1.0);
1361 else a=r->p[6]*x;
1362 f+=a;
1363 *yi=f;
1364 break;
1365 case 305:
1366 f=0;
1367 if(fabs(r->p[1])>1.0e-12) a=(r->p[0]/r->p[1])*(exp(r->p[1]*x)-1.0);
1368 else a=r->p[0]*x;
1369 f+=a;
1370 if(fabs(r->p[3])>1.0e-12) a=(r->p[2]/r->p[3])*(exp(r->p[3]*x)-1.0);
1371 else a=r->p[2]*x;
1372 f+=a;
1373 if(fabs(r->p[5])>1.0e-12) a=(r->p[4]/r->p[5])*(exp(r->p[5]*x)-1.0);
1374 else a=r->p[4]*x;
1375 f+=a;
1376 if(fabs(r->p[7])>1.0e-12) a=(r->p[6]/r->p[7])*(exp(r->p[7]*x)-1.0);
1377 else a=r->p[6]*x;
1378 f+=a;
1379 if(fabs(r->p[9])>1.0e-12) a=(r->p[8]/r->p[9])*(exp(r->p[8]*x)-1.0);
1380 else a=r->p[8]*x;
1381 f+=a;
1382 *yi=f;
1383 break;
1384 case 321:
1385 f=(r->p[0]/r->p[1])*exp(r->p[1]*x)
1386 - r->p[0]*exp((r->p[1]+r->p[2])*x)/(r->p[1]+r->p[2]);
1387 *yi=f;
1388 break;
1389 case 322:
1390 f=0.0;
1391 a=(r->p[0]/r->p[1])*exp(r->p[1]*x)
1392 - r->p[0]*exp((r->p[1]+r->p[2])*x)/(r->p[1]+r->p[2]); f+=a;
1393 a=(r->p[3]/r->p[4])*exp(r->p[4]*x)
1394 - r->p[3]*exp((r->p[4]+r->p[5])*x)/(r->p[4]+r->p[5]); f+=a;
1395 *yi=f;
1396 break;
1397 case 323:
1398 f=0.0;
1399 a=(r->p[0]/r->p[1])*exp(r->p[1]*x)
1400 - r->p[0]*exp((r->p[1]+r->p[2])*x)/(r->p[1]+r->p[2]); f+=a;
1401 a=(r->p[3]/r->p[4])*exp(r->p[4]*x)
1402 - r->p[3]*exp((r->p[4]+r->p[5])*x)/(r->p[4]+r->p[5]); f+=a;
1403 a=(r->p[6]/r->p[7])*exp(r->p[7]*x)
1404 - r->p[6]*exp((r->p[7]+r->p[8])*x)/(r->p[7]+r->p[8]); f+=a;
1405 *yi=f;
1406 break;
1407 case 331: *yi=nan(""); break; /* not yet applied */
1408 case 351: *yi=nan(""); break; /* not yet applied */
1409 case 1312:
1410 t=r->p[4]; xt=x-t; f=0.0;
1411 if(xt>0.0) {
1412 double A1, A2, L1, L2;
1413 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3];
1414 if(L1!=0.0) {
1415 a=exp(L1*xt);
1416 f+=(A2/L1)*(1.0 - a);
1417 f+=(A1/(L1*L1))*(1.0 + a*(L1*xt - 1.0));
1418 }
1419 if(L2!=0.0) f+=(A2/L2)*(exp(L2*xt) - 1.0); else f+=A2*xt;
1420 }
1421 *yi=f;
1422 break;
1423 case 1313:
1424 t=r->p[6]; xt=x-t; f=0.0;
1425 if(xt>0.0) {
1426 double A1, A2, A3, L1, L2, L3;
1427 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3]; A3=r->p[4]; L3=r->p[5];
1428 if(L1!=0.0) {
1429 a=exp(L1*xt);
1430 f+=(A2/L1)*(1.0 - a);
1431 f+=(A3/L1)*(1.0 - a);
1432 f+=(A1/(L1*L1))*(1.0 + a*(L1*xt - 1.0));
1433 }
1434 if(L2!=0.0) f+=(A2/L2)*(exp(L2*xt) - 1.0); else f+=A2*xt;
1435 if(L3!=0.0) f+=(A3/L3)*(exp(L3*xt) - 1.0); else f+=A3*xt;
1436 }
1437 *yi=f;
1438 break;
1439 case 1314:
1440 t=r->p[8]; xt=x-t; f=0.0;
1441 if(xt>0.0) {
1442 double A1, A2, A3, A4, L1, L2, L3, L4;
1443 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3]; A3=r->p[4]; L3=r->p[5];
1444 A4=r->p[6]; L4=r->p[7];
1445 if(L1!=0.0) {
1446 a=exp(L1*xt);
1447 f+=(A2/L1)*(1.0 - a);
1448 f+=(A3/L1)*(1.0 - a);
1449 f+=(A4/L1)*(1.0 - a);
1450 f+=(A1/(L1*L1))*(1.0 + a*(L1*xt - 1.0));
1451 }
1452 if(L2!=0.0) f+=(A2/L2)*(exp(L2*xt) - 1.0); else f+=A2*xt;
1453 if(L3!=0.0) f+=(A3/L3)*(exp(L3*xt) - 1.0); else f+=A3*xt;
1454 if(L4!=0.0) f+=(A4/L4)*(exp(L4*xt) - 1.0); else f+=A4*xt;
1455 }
1456 *yi=f;
1457 break;
1458 case 1401: *yi=nan(""); break; /* not yet applied */
1459 case 1402: *yi=nan(""); break; /* not yet applied */
1460 case 1431:
1461 f=0.0;
1462 if(x>0.0) f=r->p[0]*(1.0 - (r->p[1]*x + 1.0)*exp(-r->p[1]*x));
1463 *yi=f;
1464 break;
1465 case 1432:
1466 f=0.0;
1467 if(x>0.0) {
1468 a=r->p[0]/(r->p[1]*r->p[1]);
1469 f=a*(1.0 - (r->p[1]*x + 1.0)*exp(-r->p[1]*x));
1470 }
1471 *yi=f;
1472 break;
1473 case 1435:
1474 xt=x-r->p[0]; f=0.0;
1475 if(xt>0.0) {
1476 double a, b, m, n; a=r->p[1]; b=r->p[2]; m=r->p[3]; n=r->p[4];
1477 f+=m*((1.0-(n*a*xt+1.0)*exp(-n*a*xt)))/(n*n*a*a) + (1.0-exp(-a*xt))/a;
1478 f*=b;
1479 }
1480 *yi=f;
1481 break;
1482 case 1441: /* Probability density function of Erlang distribution, or actually,
1483 its integral is the cumulative distribution function of Erlang distribution */
1484 f=0.0;
1485 if(x>=0.0) {
1486 double A, k; A=r->p[0]; k=r->p[1]; if(!(k>=0.0)) return(2);
1487 int N=(int)round(r->p[2]); if(!(N>=0) || N>20) return(2); // Supports only N=0,1,2,..,20
1488 if(N==0) f=A;
1489 else if(N==1) f=A*(1.0-exp(-k*x));
1490 else if(N==2) f=A*(1.0 - exp(-k*x) - k*x*exp(-k*x));
1491 else {
1492 double s=1.0+k*x;
1493 for(int j=2; j<N; j++) s+=pow(k*x, j)/lfactorial(j);
1494 f=A*(1.0-s*exp(-k*x));
1495 }
1496 }
1497 *yi=f;
1498 break;
1499 case 2111: *yi=nan(""); break; /* no application */
1500 case 9501: *yi=nan(""); break; /* cannot be applied to one point only */
1501 case 9502: *yi=nan(""); break; /* cannot be applied to one point only */
1502 case 9503: *yi=nan(""); break; /* cannot be applied to one point only */
1503 case 9701: *yi=nan(""); break; /* cannot be applied to one point only */
1504 default:
1505 *yi=nan("");
1506 }
1507 if(isnan(*yi)) return 3;
1508 return 0;
1509}
1510/*****************************************************************************/
1511
1512/*****************************************************************************/
1517int fitIntegralEvaltac(
1519 FitVOI *r,
1521 double *x,
1523 double *yi,
1525 int dataNr
1526) {
1527 int i;
1528
1529 if(r==NULL || x==NULL || yi==NULL || dataNr<1) return(1);
1530 for(i=0; i<dataNr; i++) if(fitIntegralEval(r, x[i], yi+i)) return(2);
1531 return(0);
1532}
1533/*****************************************************************************/
1534
1535/*****************************************************************************/
1540int fitDerivEval(
1542 FitVOI *r,
1544 double x,
1546 double *yd
1547) {
1548 double a, f, xt, t;
1549
1550 if(r==NULL || yd==NULL) return 1;
1551 *yd=nan("");
1552 switch(r->type) {
1553 case 301:
1554 break;
1555 case 302:
1556 break;
1557 case 303:
1558 break;
1559 case 304:
1560 break;
1561 case 305:
1562 break;
1563 case 321:
1564 f=r->p[0]*r->p[1]*exp(r->p[1]*x)*(1.0-exp(r->p[2]*x)) - r->p[0]*r->p[2]*exp((r->p[1]+r->p[2])*x);
1565 *yd=f;
1566 break;
1567 case 322:
1568 f=0.0;
1569 a=r->p[0]*r->p[1]*exp(r->p[1]*x)*(1.0-exp(r->p[2]*x)) - r->p[0]*r->p[2]*exp((r->p[1]+r->p[2])*x); f+=a;
1570 a=r->p[3]*r->p[4]*exp(r->p[4]*x)*(1.0-exp(r->p[5]*x)) - r->p[3]*r->p[5]*exp((r->p[4]+r->p[5])*x); f+=a;
1571 *yd=f;
1572 break;
1573 case 323:
1574 f=0.0;
1575 a=r->p[0]*r->p[1]*exp(r->p[1]*x)*(1.0-exp(r->p[2]*x)) - r->p[0]*r->p[2]*exp((r->p[1]+r->p[2])*x); f+=a;
1576 a=r->p[3]*r->p[4]*exp(r->p[4]*x)*(1.0-exp(r->p[5]*x)) - r->p[3]*r->p[5]*exp((r->p[4]+r->p[5])*x); f+=a;
1577 a=r->p[6]*r->p[7]*exp(r->p[7]*x)*(1.0-exp(r->p[8]*x)) - r->p[6]*r->p[8]*exp((r->p[7]+r->p[8])*x); f+=a;
1578 *yd=f;
1579 break;
1580 case 331: break; /* not yet applied */
1581 case 351: break; /* not yet applied */
1582 case 1312:
1583 t=r->p[4]; xt=x-t; f=0.0;
1584 if(xt>0.0) {
1585 double A1, A2, L1, L2;
1586 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3];
1587 f= A1*exp(L1*xt) + A2*L2*exp(L2*xt) + (A1*xt -A2)*L1*exp(L1*xt);
1588 }
1589 *yd=f;
1590 break;
1591 case 1313:
1592 t=r->p[6]; xt=x-t; f=0.0;
1593 if(xt>0.0) {
1594 double A1, A2, A3, L1, L2, L3;
1595 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3]; A3=r->p[4]; L3=r->p[5];
1596 f= A1*exp(L1*xt) + A2*L2*exp(L2*xt) + A3*L3*exp(L3*xt) +
1597 (A1*xt -A2 -A3)*L1*exp(L1*xt);
1598 }
1599 *yd=f;
1600 break;
1601 case 1314: break;
1602 t=r->p[8]; xt=x-t; f=0.0;
1603 if(xt>0.0) {
1604 double A1, A2, A3, A4, L1, L2, L3, L4;
1605 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3]; A3=r->p[4]; L3=r->p[5];
1606 A4=r->p[6]; L4=r->p[7];
1607 f= A1*exp(L1*xt) + A2*L2*exp(L2*xt) + A3*L3*exp(L3*xt) +
1608 A4*L4*exp(L4*xt) + (A1*xt -A2 -A3 -A4)*L1*exp(L1*xt);
1609 }
1610 *yd=f;
1611 break;
1612 case 1401:
1613 case 1402:
1614 xt=x-r->p[3]; if(xt<=0.0 || r->p[2]==0.0) {
1615 f=0.0;
1616 } else {
1617 f=r->p[0]*pow(xt, r->p[1]-1.0)*exp(-xt/r->p[2])*(r->p[1]-(xt/r->p[2]));
1618 }
1619 *yd=f;
1620 break;
1621 case 1421:
1622 xt=x-r->p[0]; if(xt<=0.0) {
1623 *yd=0.0;
1624 } else {
1625 a=pow(xt/r->p[2], r->p[3]-1.0);
1626 *yd=r->p[1]*r->p[3]*a*exp(-a*xt/r->p[2])/r->p[2];
1627 }
1628 break;
1629 case 1801:
1630 xt=x-r->p[0]; if(xt<=0.0) {
1631 *yd=0.0;
1632 } else {
1633 a=pow(r->p[2], r->p[3])+pow(xt, r->p[3]);
1634 *yd=r->p[1]*r->p[3]*(pow(xt, r->p[3]-1.0)/a - pow(xt, 2.0*r->p[3]-1.0)/(a*a));
1635 }
1636 break;
1637 case 2111: *yd=nan(""); break; /* no application */
1638 case 9501: *yd=nan(""); break; /* cannot be applied to one point only */
1639 case 9502: *yd=nan(""); break; /* cannot be applied to one point only */
1640 case 9503: *yd=nan(""); break; /* cannot be applied to one point only */
1641 case 9701: *yd=nan(""); break; /* cannot be applied to one point only */
1642 default:
1643 *yd=nan("");
1644 }
1645 if(isnan(*yd)) return(3);
1646 return(0);
1647}
1648/*****************************************************************************/
1649
1650/*****************************************************************************/
1655int fitDerivEvaltac(
1657 FitVOI *r,
1659 double *x,
1661 double *yd,
1663 int dataNr
1664) {
1665 int i;
1666
1667 if(r==NULL || x==NULL || yd==NULL || dataNr<1) return(1);
1668 for(i=0; i<dataNr; i++) if(fitDerivEval(r, x[i], yd+i)) return(2);
1669 return(0);
1670}
1671/*****************************************************************************/
1672
1673/*****************************************************************************/
1678unsigned int factorial(
1682 unsigned int n
1683) {
1684 if(n<1) return(1);
1685 if(n>12) return(0);
1686 return(n*factorial(n-1));
1687}
1688/*****************************************************************************/
1689
1690/*****************************************************************************/
1695unsigned long long int lfactorial(
1699 unsigned long long int n
1700) {
1701 if(n<1) return(1);
1702 if(n>20) return(0);
1703 return(n*lfactorial(n-1));
1704}
1705/*****************************************************************************/
1706
1707/*****************************************************************************/
1715double igam(
1717 double a,
1719 double x
1720) {
1721 /* Check parameters */
1722 if(x==0.0) return(0.0);
1723 if((x<0.0) || (a<=0.0)) return(nan(""));
1724
1725 if((x>1.0) && (x>a)) return(1.0 - igamc(a, x));
1726
1727 /* Left tail of incomplete Gamma function:
1728 x^a * e^-x * Sum(k=0,Inf)(x^k / Gamma(a+k+1)) */
1729
1730 double ans, ax, c, r;
1731
1732 /* Compute x**a * exp(-x) / Gamma(a) */
1733 ax=a*log(x) - x - lgamma(a);
1734 if(ax<-DBL_MAX_10_EXP) return(0.0); // underflow
1735 ax=exp(ax);
1736
1737 /* power series */
1738 r=a; c=1.0; ans=1.0;
1739 do {
1740 r+=1.0;
1741 c*=x/r;
1742 ans+=c;
1743 } while(c/ans > DBL_EPSILON);
1744 return(ans*ax/a);
1745}
1746/*****************************************************************************/
1747
1748/*****************************************************************************/
1756double igamc(
1758 double a,
1760 double x
1761) {
1762 /* Check parameters */
1763 if((x<0.0) || (a<=0.0)) return(nan(""));
1764
1765 if((x<1.0) || (x<a)) return(1.0-igam(a, x));
1766
1767 double ans, ax, c, yc, r, t, y, z;
1768 double pk, pkm1, pkm2, qk, qkm1, qkm2;
1769 double big=4.503599627370496E+015;
1770 double biginv=2.22044604925031308085E-016;
1771
1772 ax = a*log(x) - x - lgamma(a);
1773 if(ax < -DBL_MAX_10_EXP) return(0.0); // underflow
1774 ax=exp(ax);
1775
1776 /* continued fraction */
1777 y=1.0-a; z=x+y+1.0; c=0.0;
1778 pkm2=1.0; qkm2=x; pkm1=x+1.0; qkm1=z*x;
1779 ans=pkm1/qkm1;
1780 do {
1781 c+=1.0; y+=1.0; z+=2.0;
1782 yc=y*c; pk=pkm1*z - pkm2*yc; qk=qkm1*z - qkm2*yc;
1783 if(qk!=0.0) {r=pk/qk; t=fabs((ans-r)/r); ans=r;}
1784 else t=1.0;
1785 pkm2=pkm1; pkm1=pk; qkm2=qkm1; qkm1=qk;
1786 if(fabs(pk)>big) {pkm2*=biginv; pkm1*=biginv; qkm2*=biginv; qkm1*=biginv;}
1787 } while(t>DBL_EPSILON);
1788 return(ans*ax);
1789}
1790/*****************************************************************************/
1791
1792/*****************************************************************************/
backupExistingFile
int backupExistingFile(char *filename, char *backup_ext, char *status)
Definition backup.c:14
get_datetime
int get_datetime(char *str, struct tm *date, int verbose)
Definition datetime.c:322
ctime_r_int
char * ctime_r_int(const time_t *t, char *buf)
Convert calendard time t into a null-terminated string of the form YYYY-MM-DD hh:mm:ss,...
Definition datetime.c:110
timegm
time_t timegm(struct tm *tm)
Inverse of gmtime, converting struct tm to time_t.
Definition datetime.c:69
atof_dpi
double atof_dpi(char *str)
Definition decpoint.c:59
libtpccurveio.h
Header file for libtpccurveio.
rnameCatenate
int rnameCatenate(char *rname, int max_rname_len, char *name1, char *name2, char *name3, char space)
Definition rname.c:189
MAX_REGIONNAME_LEN
#define MAX_REGIONNAME_LEN
Definition libtpcmisc.h:154
strReplaceChar
void strReplaceChar(char *str, char c1, char c2)
Definition strext.c:159
strTokenNCpy
int strTokenNCpy(const char *str1, const char *str2, int i, char *str3, int count)
Definition strext.c:45
strTokenNr
int strTokenNr(const char *str1, const char *str2)
Definition strext.c:17
strlcpy
size_t strlcpy(char *dst, const char *src, size_t dstsize)
Definition strext.c:245
strTokenDup
char * strTokenDup(const char *s1, const char *s2, int *next)
Definition strext.c:89
petTunit
char * petTunit(int tunit)
Definition petunits.c:226
MAX_STUDYNR_LEN
#define MAX_STUDYNR_LEN
Definition libtpcmisc.h:163
petTunitId
int petTunitId(const char *timeunit)
Definition petunits.c:187
MAX_REGIONSUBNAME_LEN
#define MAX_REGIONSUBNAME_LEN
Definition libtpcmisc.h:158
libtpcmodel.h
Header file for libtpcmodel.
simDispersion
int simDispersion(double *x, double *y, int n, double tau1, double tau2, double *tmp)
Definition simulate.c:1911
igam
double igam(double a, double x)
Definition mathfunc.c:1715
lfactorial
unsigned long long int lfactorial(unsigned long long int n)
Definition mathfunc.c:1695
fitDerivEval
int fitDerivEval(FitVOI *r, double x, double *yd)
Definition mathfunc.c:1540
fitSetmem
int fitSetmem(FIT *fit, int voiNr)
Definition mathfunc.c:154
fitEvaltac
int fitEvaltac(FitVOI *r, double *x, double *y, int dataNr)
Definition mathfunc.c:1252
fitEval
int fitEval(FitVOI *r, double x, double *y)
Definition mathfunc.c:645
igamc
double igamc(double a, double x)
Definition mathfunc.c:1756
fitEmpty
void fitEmpty(FIT *fit)
Definition mathfunc.c:18
fitRead
int fitRead(char *filename, FIT *fit, int verbose)
Definition mathfunc.c:196
fitWrite
int fitWrite(FIT *fit, char *filename)
Definition mathfunc.c:54
fitIntegralEvaltac
int fitIntegralEvaltac(FitVOI *r, double *x, double *yi, int dataNr)
Definition mathfunc.c:1517
fitIntegralEval
int fitIntegralEval(FitVOI *r, double x, double *yi)
Definition mathfunc.c:1308
MATHFUNC_TEST
int MATHFUNC_TEST
Definition mathfunc.c:5
fitFunctionformat
int fitFunctionformat(int type, char *str)
Definition mathfunc.c:381
fiterrmsg
char fiterrmsg[64]
Definition mathfunc.c:6
fitInit
void fitInit(FIT *fit)
Definition mathfunc.c:38
fitDerivEvaltac
int fitDerivEvaltac(FitVOI *r, double *x, double *yd, int dataNr)
Definition mathfunc.c:1655
factorial
unsigned int factorial(unsigned int n)
Definition mathfunc.c:1678
fitFunctionname
int fitFunctionname(int type, char *str)
Definition mathfunc.c:539
fitPrint
void fitPrint(FIT *fit)
Definition mathfunc.c:180
FIT::timeunit
int timeunit
Definition libtpccurveio.h:319
FIT::voiNr
int voiNr
Definition libtpccurveio.h:311
FIT::_voidataNr
int _voidataNr
Definition libtpccurveio.h:327
FIT::datafile
char datafile[FILENAME_MAX]
Definition libtpccurveio.h:313
FIT::program
char program[1024]
Definition libtpccurveio.h:325
FIT::unit
char unit[MAX_UNITS_LEN+1]
Definition libtpccurveio.h:317
FIT::studynr
char studynr[MAX_STUDYNR_LEN+1]
Definition libtpccurveio.h:315
FIT::voi
FitVOI * voi
Definition libtpccurveio.h:321
FIT::time
time_t time
Definition libtpccurveio.h:323
FitVOI::wss
double wss
Definition libtpccurveio.h:299
FitVOI::end
double end
Definition libtpccurveio.h:293
FitVOI::place
char place[MAX_REGIONSUBNAME_LEN+1]
Definition libtpccurveio.h:285
FitVOI::hemisphere
char hemisphere[MAX_REGIONSUBNAME_LEN+1]
Definition libtpccurveio.h:283
FitVOI::voiname
char voiname[MAX_REGIONSUBNAME_LEN+1]
Definition libtpccurveio.h:281
FitVOI::name
char name[MAX_REGIONNAME_LEN+1]
Definition libtpccurveio.h:279
FitVOI::start
double start
Definition libtpccurveio.h:291
FitVOI::p
double p[MAX_FITPARAMS]
Definition libtpccurveio.h:297
Generated on Fri Sep 19 2025 13:32:04 for TPCCLIB by doxygen 1.12.0