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 /* Read optional fields until "Nr of VOIs:" */
233 while(1) {
234 /* Read title field */
235 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
236 /* Stop when "Nr of VOIs" is reached */
237 if(strncasecmp(line, "Nr of VOIs:", 11)==0) break;
238 /* Read fit date and time */
239 if(strncasecmp(line, "Date:", 5)==0) {
240 cptr=line+5; while(*cptr && !isdigit(*cptr)) cptr++;
241 if(get_datetime(cptr, &st, verbose-3)==0) fit->time=timegm(&st);
242 continue;
243 }
244 /* Read the name of the original datafile */
245 if(strncasecmp(line, "Data file:", 10)==0) {
246 lptr=&line[10]; cptr=strtok(lptr, " \t\n\r");
247 if(cptr!=NULL && strlen(cptr)<FILENAME_MAX) strcpy(fit->datafile, cptr);
248 continue;
249 }
250 /* Read study number */
251 if(strncasecmp(line, "Study number:", 13)==0) {
252 lptr=&line[13]; cptr=strtok(lptr, " \t\n\r");
253 if(cptr!=NULL && strlen(cptr)<MAX_STUDYNR_LEN+1) strcpy(fit->studynr, cptr);
254 continue;
255 }
256 /* Read the activity unit */
257 if(strncasecmp(line, "Data unit:", 10)==0) {
258 lptr=&line[10]; cptr=strtok(lptr, " \t\n\r");
259 if(cptr!=NULL && strlen(cptr)<1024) strcpy(fit->unit, cptr);
260 continue;
261 }
262 /* Read the time unit */
263 if(strncasecmp(line, "Time unit:", 10)==0) {
264 lptr=&line[10]; cptr=strtok(lptr, " \t\n\r");
265 fit->timeunit=petTunitId(cptr);
266 continue;
267 }
268 if(strncasecmp(line, "Distance unit:", 14)==0) {
269 lptr=&line[14]; cptr=strtok(lptr, " \t\n\r");
270 fit->timeunit=petTunitId(cptr);
271 continue;
272 }
273 fclose(fp); return 8;
274 }
275 /* Read the nr of regions */
276 if(strncasecmp(line, "Nr of VOIs:", 11)) {fclose(fp); return 8;}
277 lptr=&line[11]; cptr=strtok(lptr, " \t\n\r");
278 n=atoi(cptr); if(n<1 || n>32000) {fclose(fp); return 8;}
279 /* Allocate memory for regions */
280 if(fitSetmem(fit, n)) {strcpy(fiterrmsg, "out of memory"); fclose(fp); return 9;}
281 fit->voiNr=n;
282 /* Read (and ignore) title line */
283 strcpy(fiterrmsg, "wrong format");
284 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
285 if(strncasecmp(line, "Region", 6)) {fclose(fp); return 10;}
286 /* Read regional data */
287 for(ri=0; ri<fit->voiNr; ri++) {
288 /* Read line of data */
289 line[0]=(char)0;
290 while(fgets(line, 1024, fp)!=NULL) if(strlen(line)>2 && line[0]!='#') break;
291 if(!line[0]) break;
292 int pindx=0; char seps[8];
293 int sn=strTokenNr(line, " \t\n\r"); if(sn<8) {fclose(fp); fitEmpty(fit); return 11;}
294 int tn=strTokenNr(line, "\t\n\r");
295 if(tn==sn || tn==sn-1 || tn==sn-2) { // tab as separator
296 strcpy(seps, "\t\n\r"); n=tn;
297 char *s=strTokenDup(line, seps, NULL);
298 char *cptr=s;
299 int i=0, n=strlen(s);
300 strReplaceChar(cptr, ' ', (char)0);
301 strlcpy(fit->voi[ri].voiname, cptr, MAX_REGIONSUBNAME_LEN+1);
302 i+=strlen(fit->voi[ri].voiname);
303 cptr=s+i; if(i<n) {cptr++; i++;}
304 strReplaceChar(cptr, ' ', (char)0);
305 strlcpy(fit->voi[ri].hemisphere, cptr, MAX_REGIONSUBNAME_LEN+1);
306 i+=strlen(fit->voi[ri].hemisphere);
307 cptr=s+i; if(i<n) {cptr++; i++;}
308 strReplaceChar(cptr, ' ', (char)0);
309 strlcpy(fit->voi[ri].place, cptr, MAX_REGIONSUBNAME_LEN+1);
310 free(s);
311 pindx=2;
312 } else { // spaces as separators
313 strcpy(seps, " \t\n\r"); n=sn;
314 strTokenNCpy(line, seps, 1, fit->voi[ri].voiname, MAX_REGIONSUBNAME_LEN+1);
315 strTokenNCpy(line, seps, 2, fit->voi[ri].hemisphere, MAX_REGIONSUBNAME_LEN+1);
316 strTokenNCpy(line, seps, 3, fit->voi[ri].place, MAX_REGIONSUBNAME_LEN+1);
317 pindx=4;
318 }
319 rnameCatenate(fit->voi[ri].name, MAX_REGIONNAME_LEN+1,
320 fit->voi[ri].voiname, fit->voi[ri].hemisphere, fit->voi[ri].place, ' ');
321 /* Fit start and end times, and original data nr */
322 char s[128];
323 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].start=atof_dpi(s);
324 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].end=atof_dpi(s);
325 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].dataNr=atoi(s);
326 /* Fit error, parameter nr and function number (type) */
327 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].wss=atof_dpi(s);
328 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].parNr=atoi(s);
329 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].type=atoi(s);
330 /* Parameters */
331 for(i=0; i<fit->voi[ri].parNr; i++) {
332 strTokenNCpy(line, seps, pindx++, s, 128); fit->voi[ri].p[i]=atof_dpi(s);
333 }
334 }
335 if(ri==0) {fclose(fp); fitEmpty(fit); return(12);}
336 if(ri<fit->voiNr) fit->voiNr=ri;
337
338 /* Close file */
339 fclose(fp);
340 strcpy(fiterrmsg, "");
341
342 if(verbose>1) printf("done fitRead()\n");
343 return(0);
344}
345/*****************************************************************************/
346
347/*****************************************************************************/
351int fitFunctionformat(
353 int type,
355 char *str
356) {
357 strcpy(str, "");
358 switch(type) {
359 /* Polynomials, including line */
360 case 100: strcpy(str, "f(x)=A"); break;
361 case 101: strcpy(str, "f(x)=A+B*x"); break;
362 case 102: strcpy(str, "f(x)=A+B*x+C*x^2"); break;
363 case 103: strcpy(str, "f(x)=A+B*x+C*x^2+D*x^3"); break;
364 case 104: strcpy(str, "f(x)=A+B*x+C*x^2+D*x^3+E*x^4"); break;
365 case 105: strcpy(str, "f(x)=A+B*x+C*x^2+D*x^3+E*x^4+F*x^5"); break;
366 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;
367 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;
368 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;
369 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;
370 /* Rational functions */
371 case 211: strcpy(str, "f(x)=(A+C*x)/(B+D*x)"); break;
372 case 221: strcpy(str, "f(x)=(A+C*x+E*x^2)/(B+D*x)"); break;
373 case 222: strcpy(str, "f(x)=(A+C*x+E*x^2)/(B+D*x+F*x^2)"); break;
374 case 232: strcpy(str, "f(x)=(A+C*x+E*x^2+G*x^3)/(B+D*x+F*x^2)"); break;
375 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;
376 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;
377 /* Rational function for plasma-to-blood ratio */
378 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;
379 /* Exponential functions */
380 case 301: strcpy(str, "f(x)=A*exp(B*x)"); break;
381 case 302: strcpy(str, "f(x)=A*exp(B*x)+C*exp(D*x)"); break;
382 case 303: strcpy(str, "f(x)=A*exp(B*x)+C*exp(D*x)+E*exp(F*x)"); break;
383 case 304: strcpy(str, "f(x)=A*exp(B*x)+C*exp(D*x)+E*exp(F*x)+G*exp(H*x)"); break;
384 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;
385 /* Feng function */
386 case 1312: strcpy(str, "f(x)=(A*(x-t)-C)*exp(B*(x-t))+C*exp(D*(x-t))"); break;
387 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;
388 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;
389 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;
390 /* Lundqvist function */
391 case 321: strcpy(str, "f(x)=A*exp(B*x)*(1-exp(C*x))"); break;
392 case 322: strcpy(str, "f(x)=A*exp(B*x)*(1-exp(C*x))+D*exp(E*x)*(1-exp(F*x))"); break;
393 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;
394 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;
395 /* Exponential bolus infusion functions */
396 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");
397 break;
398 /* Kudomi's function for radiowater */
399 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");
400 break;
401 /* Bolus infusion approaching zero */
402 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");
403 break;
404 /* Exponential functions for plasma fractions */
405 case 351: strcpy(str, "f(x)=1-a*(2-exp(-b*x)-exp(-c*x))"); break;
406 /* Inverted gamma cdf for plasma parent fraction */
407 case 403: strcpy(str, "f(x)=A*(1-B*gammacdf(D,C*(x-t)))"); break;
408 /* Gamma variate function */
409 case 1401: strcpy(str, "f(x)=A*((x-D)^B)*exp(-(x-D)/C) , when x>=D, else f(x)=0"); break;
410 /* Gamma variate function with background */
411 case 1402: strcpy(str, "f(x)=A*((x-D)^B)*exp(-(x-D)/C) + E , when x>=D, else f(x)=E"); break;
412 /* Gamma variate bolus plus recirculation function */
413 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;
414 /* Weibull cdf */
415 case 1421: strcpy(str, "f(x)=A*(1-exp(-((x-t)/B)^C) , when x>t, else f(x)=0"); break;
416 /* Weibull cdf plus pdf (derivative of cdf) */
417 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;
418 /* Sum of Surge functions with delay */
419 case 1430: strcpy(str, "f(x)=A*(x-t)*exp(-B*(x-t)) + C*(x-t)*exp(-D*(x-t)) + ... , when x-t>0, else f(x)=0"); break;
420 /* Surge function with AUC=A */
421 case 1431: strcpy(str, "f(x)=A*x*exp(-B*x)*B^2 , when x>0, else f(x)=0"); break;
422 /* Traditional Surge function */
423 case 1432: strcpy(str, "f(x)=A*x*exp(-B*x) , when x>0, else f(x)=0"); break;
424 /* Surge function with recirculation */
425 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;
426 /* Surge function with recirculation for plasma-to-blood ratio */
427 case 1434: strcpy(str, "f(x)=1/(1-H*(1-r(x))), where r(x) is function for RBC-to-plasma"); break;
428 /* Surge function for late FDG input */
429 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;
430 /* Probability density function of Erlang distribution */
431 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;
432 /* Hill functions for TACs */
433 case 1801: strcpy(str, "f(x)=[A*(x-t)^B]/[(x-t)^B+C^B]"); break;
434 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;
435 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;
436 /* Surge function with recirculation and delay */
437 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;
438 /* Hill functions for dose-response curves */
439 case 2801: strcpy(str, "f(x)=B+[A-B]/[1+(C/x)^D]"); break;
440 case 2802: strcpy(str, "f(x)=B+[A-B]/[1+10^{(C-x)*D}]"); break;
441 /* Exponential/power type functions for parent fractions */
442 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;
443 /* Hill type functions for fractions */
444 case 841: strcpy(str, "f(x)=(A*x^B)/(x^B+C)"); break;
445 case 842: strcpy(str, "f(x)=1-((A*x^B)/(x^B+C))"); break;
446 case 843: strcpy(str, "f(x)=1-((A*(1+D*x)*x^B)/(x^B+C))"); break;
447 case 844: strcpy(str, "f(x)=(A*(x-t)^B)/((x-t)^B+C)+D, when x>t, else f(x)=D"); break;
448 case 845: strcpy(str, "f(x)=A-(A*x^B)/(x^B+C))"); break;
449 case 846: strcpy(str, "f(x)=D+((A-D)*(x-t)^B)/((x-t)^B+C), when x>t, else f(x)=D"); break;
450 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;
451 case 848: strcpy(str, "f(x)=D*((1-A)*(x-t)^B)/((x-t)^B+C), when x>t, else f(x)=D"); break;
452 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;
453 /* Hill function for plasma-to-blood ratio */
454 case 2841: strcpy(str, "f(x)=1/(1-H*(1-r(x))), where r(x) is Hill function for RBC-to-plasma"); break;
455 /* Mamede/Watabe function for fractions */
456 case 851: strcpy(str, "f(x)=1/(1+(A*x)^2)^B"); break;
457 case 852: strcpy(str, "f(x)=1-1/(1+(A*x)^2)^B"); break;
458 /* Mamede/Watabe function for fractions, as extended by Meyer */
459 case 861: strcpy(str, "f(x)=(1+(A*(x-t))^B)^(-C), when x>t, else f(x)=1"); break;
460 case 862: strcpy(str, "f(x)=1-(1+(A*(x-t))^B)^(-C), when x>t, else f(x)=0"); break;
461 /* ... and further extended by letting fraction start somewhere between 0 and 1 */
462 case 863: strcpy(str, "f(x)=(D^(-1/C)+(A*(x-t))^B)^(-C), when x>t, else f(x)=D"); break;
463 case 864: strcpy(str, "f(x)=1-(D^(-1/C)+(A*(x-t))^B)^(-C), when x>t, else f(x)=1-D"); break;
464 /* Functions for fitting plasma fractions via separate metabolite fractions */
465 case 871:
466 case 881:
467 strcpy(str, "f(x)=1-f1(x)-f2(x)-f3(x)"); break;
468 case 872:
469 case 882:
470 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;
471 case 873:
472 case 883:
473 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;
474 case 874:
475 case 884:
476 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;
477 case 1010:
478 strcpy(str, "f(x)=p2 when x>=p1, f(x)=p4 when x>=p2, ..."); break;
479 break;
480 /* PET profile functions */
481 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");
482 break;
483 /* Combined functions and models */
484 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");
485 break;
486 /* Compartmental model functions */
487 /* Graham's plasma curve smoothing function */
488 case 9501: strcpy(str, "Cp(t)<=>Ci(t)<=>Ct(t)"); break;
489 /* Extended Graham's plasma curve smoothing function */
490 case 9502: strcpy(str, "Ce(t)<=>Cp(t)<=>Ci(t)<=>Ct(t)"); break;
491 /* Extended Graham's plasma curve smoothing function with metabolite */
492 case 9503: strcpy(str, "Cpa(t)<=>Cia(t)<=>Cta(t)->Ctm(t)<=>Cim(t)<=>Cpm(t)"); break;
493 /* Huang's plasma metabolite model */
494 case 9601: strcpy(str, "C4(t)<=>C3(t)<-C0(t)->C1(t)<=>C2(t)"); break;
495 /* Extended Carson's plasma metabolite model */
496 case 9602: strcpy(str, "Cpa(t)<=>Cta(t)->Ctm(t)<=>Cpm(t)"); break;
497 /* New plasma metabolite model */
498 case 9603: strcpy(str, "Cpa(t)->Ct1(t)<=>Cpm(t)<=>Ct2(t)"); break;
499 /* Multi-linear multi-compartmental TAC fitting model */
500 case 9701: strcpy(str, "Ideal bolus -> n compartments"); break;
501 default: return(1);
502 }
503 return(0);
504}
505/*****************************************************************************/
506
507/*****************************************************************************/
511int fitFunctionname(
513 int type,
515 char *str
516) {
517 strcpy(str, "");
518 switch(type) {
519 case 100: strcpy(str, "f(x)=A"); break;
520 case 101: strcpy(str, "line"); break;
521 case 102: strcpy(str, "2nd order polynomial"); break;
522 case 103: strcpy(str, "3rd order polynomial"); break;
523 case 104: strcpy(str, "4th order polynomial"); break;
524 case 105: strcpy(str, "5th order polynomial"); break;
525 case 106: strcpy(str, "6th order polynomial"); break;
526 case 107: strcpy(str, "7th order polynomial"); break;
527 case 108: strcpy(str, "8th order polynomial"); break;
528 case 109: strcpy(str, "9th order polynomial"); break;
529 case 211: strcpy(str, "1/1 order rational function"); break;
530 case 221: strcpy(str, "2/1 order rational function"); break;
531 case 222: strcpy(str, "2/2 order rational function"); break;
532 case 232: strcpy(str, "3/2 order rational function"); break;
533 case 233: strcpy(str, "3/3 order rational function"); break;
534 case 1232: strcpy(str, "3/2 order rational function with delay"); break;
535 case 2233: strcpy(str, "3/3 order rational function for plasma-to-blood ratio"); break;
536 case 301: strcpy(str, "exponential function"); break;
537 case 302: strcpy(str, "sum of 2 exponential functions"); break;
538 case 303: strcpy(str, "sum of 3 exponential functions"); break;
539 case 304: strcpy(str, "sum of 4 exponential functions"); break;
540 case 305: strcpy(str, "sum of 5 exponential functions"); break;
541 case 1312: strcpy(str, "Feng model 2 function with 2 exponentials"); break;
542 case 1313: strcpy(str, "Feng model 2 function"); break;
543 case 1314: strcpy(str, "Feng model 2 function with 4 exponentials"); break;
544 case 2313: strcpy(str, "Feng M2 function for plasma-to-blood ratio"); break;
545 case 321: strcpy(str, "Lundqvist function"); break;
546 case 322: strcpy(str, "sum of 2 Lundqvist functions"); break;
547 case 323: strcpy(str, "sum of 3 Lundqvist functions"); break;
548 case 1321: strcpy(str, "Lundqvist function with integral and delay"); break;
549 case 331: strcpy(str, "Exponential bolus infusion function"); break;
550 case 332: strcpy(str, "Kudomi's exponential bolus infusion function for radiowater"); break;
551 case 334: strcpy(str, "Exponential bolus function approaching zero"); break;
552 case 351: strcpy(str, "Exponential function for [C-11]PK11195 plasma fractions"); break;
553 case 403: strcpy(str, "Inverted gamma cdf for plasma parent fraction"); break;
554 case 1401: strcpy(str, "Gamma variate function"); break;
555 case 1402: strcpy(str, "Gamma variate with background"); break;
556 case 1403: strcpy(str, "Gamma variate bolus plus recirculation"); break;
557 case 1421: strcpy(str, "Weibull cdf with delay"); break;
558 case 1423: strcpy(str, "Weibull cdf and derivative with delay"); break;
559 case 1430: strcpy(str, "Sum of surge functions with delay"); break;
560 case 1431: strcpy(str, "Surge function"); break;
561 case 1432: strcpy(str, "Surge function (trad)"); break;
562 case 1433: strcpy(str, "Surge function with recirculation"); break;
563 case 1434: strcpy(str, "Surge function with recirculation for plasma-to-blood ratio"); break;
564 case 1435: strcpy(str, "Surge function for late FDG AIF"); break;
565 case 1441: strcpy(str, "Probability density function of Erlang distribution"); break;
566 case 1801: strcpy(str, "Hill function with delay"); break;
567 case 1811: strcpy(str, "Derivative of Hill function with delay"); break;
568 case 1821: strcpy(str, "Sum of Hill function and derivative with delay"); break;
569 case 1833: strcpy(str, "Surge function with recirculation and delay"); break;
570 case 2801: strcpy(str, "Hill function for dose-response curve on linear scale"); break;
571 case 2802: strcpy(str, "Hill function for dose-response curve on log scale"); break;
572 case 831: strcpy(str, "Exponential/Power function for parent fractions"); break;
573 case 841: strcpy(str, "Hill function"); break;
574 case 842: strcpy(str, "Hill function (1-f(x))"); break;
575 case 843: strcpy(str, "Hill function (1-f(x)) with ascending or descending end"); break;
576 case 844: strcpy(str, "Hill function with background"); break;
577 case 845: strcpy(str, "Hill function (A-f(x))"); break;
578 case 846: strcpy(str, "Extended Hill function for plasma parent fraction"); break;
579 case 847: strcpy(str, "Extended Hill function for plasma metabolite fraction"); break;
580 case 848: strcpy(str, "Extended Hill function #2 for plasma parent fraction"); break;
581 case 849: strcpy(str, "Extended Hill function #2 for plasma metabolite fraction"); break;
582 case 851: strcpy(str, "Mamede function"); break;
583 case 852: strcpy(str, "Mamede function (1-f(x)"); break;
584 case 861: strcpy(str, "Meyer parent fraction function"); break;
585 case 862: strcpy(str, "Meyer metabolite fraction function"); break;
586 case 863: strcpy(str, "Extended Meyer parent fraction function"); break;
587 case 864: strcpy(str, "Extended Meyer metabolite fraction function"); break;
588 case 871: strcpy(str, "1-3 metabolite Hill function for parent"); break;
589 case 872: strcpy(str, "1-3 metabolite Hill function for metab1"); break;
590 case 873: strcpy(str, "1-3 metabolite Hill function for metab2"); break;
591 case 874: strcpy(str, "1-3 metabolite Hill function for metab3"); break;
592 case 881: strcpy(str, "1-3 metabolite power function for parent"); break;
593 case 882: strcpy(str, "1-3 metabolite power function for metab1"); break;
594 case 883: strcpy(str, "1-3 metabolite power function for metab2"); break;
595 case 884: strcpy(str, "1-3 metabolite power function for metab3"); break;
596 case 1010: strcpy(str, "Step function"); break;
597 case 2111: strcpy(str, "Image profile function"); break;
598 case 2841: strcpy(str, "Hill function for plasma-to-blood ratio"); break;
599 case 3331: strcpy(str, "Exponential bolus infusion function with delay and dispersion"); break;
600 case 9501: strcpy(str, "Graham's input function"); break;
601 case 9502: strcpy(str, "Extended Graham's input function"); break;
602 case 9503: strcpy(str, "Graham's input function with metabolite"); break;
603 case 9601: strcpy(str, "Huang's plasma metabolite model"); break;
604 case 9602: strcpy(str, "Extended Carson's plasma metabolite model"); break;
605 case 9603: strcpy(str, "New plasma metabolite model"); break;
606 case 9701: strcpy(str, "Multilinear multicompartmental TAC fitting model"); break;
607 default: return(1);
608 }
609 return(0);
610}
611/*****************************************************************************/
612
613/*****************************************************************************/
618int fitEval(
620 FitVOI *r,
622 double x,
624 double *y
625) {
626 double a, b, f, sqrx, cubx, xt;
627 int i, j, m, n;
628 double mp[3][5];
629 double mf[]={0.0,0.0,0.0};
630
631 if(r==NULL) return(1);
632 if(MATHFUNC_TEST>0) {
633 printf("fitEval(r, %g, %g): type=%d ;", x, *y, r->type);
634 for(i=0; i<r->parNr; i++) printf(" %g", r->p[i]);
635 printf("\n");
636 }
637 if(y==NULL) return(1);
638 *y=nan("");
639 switch(r->type) {
640 case 100:
641 case 101:
642 case 102:
643 case 103:
644 case 104:
645 case 105:
646 case 106:
647 case 107:
648 case 108:
649 case 109:
650 n=r->type-99; f=0.0; for(i=n-1; i>0; i--) {f+=r->p[i]; f*=x;} f+=r->p[0];
651 *y=f;
652 break;
653 case 211:
654 a=r->p[0]+r->p[2]*x; b=r->p[1]+r->p[3]*x; if(b==0.0) break;
655 f=a/b; *y=f;
656 break;
657 case 221:
658 sqrx=x*x;
659 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;
660 f=a/b; *y=f;
661 break;
662 case 222:
663 sqrx=x*x;
664 a=r->p[0]+r->p[2]*x+r->p[4]*sqrx;
665 b=r->p[1]+r->p[3]*x+r->p[5]*sqrx; if(b==0.0) break;
666 f=a/b; *y=f;
667 break;
668 case 232:
669 sqrx=x*x; cubx=sqrx*x;
670 a=r->p[0]+r->p[2]*x+r->p[4]*sqrx+r->p[6]*cubx;
671 b=r->p[1]+r->p[3]*x+r->p[5]*sqrx; if(b==0.0) break;
672 f=a/b; *y=f;
673 break;
674 case 233:
675 sqrx=x*x; cubx=sqrx*x;
676 a=r->p[0]+r->p[2]*x+r->p[4]*sqrx+r->p[6]*cubx;
677 b=r->p[1]+r->p[3]*x+r->p[5]*sqrx+r->p[7]*cubx; if(b==0.0) break;
678 f=a/b; *y=f;
679 break;
680 case 301:
681 f=r->p[0]*exp(r->p[1]*x); *y=f;
682 break;
683 case 302:
684 f=0;
685 a=r->p[0]*exp(r->p[1]*x); f+=a;
686 a=r->p[2]*exp(r->p[3]*x); f+=a;
687 *y=f;
688 break;
689 case 303:
690 f=0;
691 a=r->p[0]*exp(r->p[1]*x); f+=a;
692 a=r->p[2]*exp(r->p[3]*x); f+=a;
693 a=r->p[4]*exp(r->p[5]*x); f+=a;
694 *y=f;
695 break;
696 case 304:
697 f=0;
698 a=r->p[0]*exp(r->p[1]*x); f+=a;
699 a=r->p[2]*exp(r->p[3]*x); f+=a;
700 a=r->p[4]*exp(r->p[5]*x); f+=a;
701 a=r->p[6]*exp(r->p[7]*x); f+=a;
702 *y=f;
703 break;
704 case 305:
705 f=0;
706 a=r->p[0]*exp(r->p[1]*x); f+=a;
707 a=r->p[2]*exp(r->p[3]*x); f+=a;
708 a=r->p[4]*exp(r->p[5]*x); f+=a;
709 a=r->p[6]*exp(r->p[7]*x); f+=a;
710 a=r->p[8]*exp(r->p[9]*x); f+=a;
711 *y=f;
712 break;
713 case 321:
714 f=r->p[0]*exp(r->p[1]*x)*(1.0-exp(r->p[2]*x));
715 *y=f;
716 break;
717 case 322:
718 f=0.0;
719 a=r->p[0]*exp(r->p[1]*x)*(1.0-exp(r->p[2]*x)); f+=a;
720 a=r->p[3]*exp(r->p[4]*x)*(1.0-exp(r->p[5]*x)); f+=a;
721 *y=f;
722 break;
723 case 323:
724 f=0.0;
725 a=r->p[0]*exp(r->p[1]*x)*(1.0-exp(r->p[2]*x)); f+=a;
726 a=r->p[3]*exp(r->p[4]*x)*(1.0-exp(r->p[5]*x)); f+=a;
727 a=r->p[6]*exp(r->p[7]*x)*(1.0-exp(r->p[8]*x)); f+=a;
728 *y=f;
729 break;
730 case 1321:
731 /* A=p0 B=p1 C=p2 D=k=p3 dT=p4 */
732 xt=x-r->p[4]; if(xt<=0.0) {
733 f=0.0;
734 } else {
735 a=exp(-r->p[2]*xt); f=r->p[0]*exp(-r->p[1]*xt)*(1.0-a);
736 if(r->p[3]>0.0)
737 f+=r->p[3]*(r->p[0]/(r->p[1]*(r->p[1]+r->p[2])))
738 *(r->p[2]-((r->p[1]+r->p[2])/a-r->p[1])*exp(-(r->p[1]+r->p[2])*xt));
739 }
740 *y=f;
741 break;
742 case 331:
743 n=(r->parNr-2)/2;
744 if(x<=r->p[0])
745 *y=0.0;
746 else if(x<r->p[0]+r->p[1]) {
747 for(i=0, f=0.0; i<n; i++) {
748 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];
749 f+=b*(1.0-exp(-r->p[2*i+3]*(x-r->p[0])));
750 }
751 if(r->p[1]>0.0) f*=(1.0/r->p[1]);
752 *y=f;
753 } else {
754 for(i=0, f=0.0; i<n; i++) {
755 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];
756 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])));
757 }
758 if(r->p[1]>0.0) f*=(1.0/r->p[1]);
759 *y=f;
760 }
761 break;
762 case 332:
763 if(x<=r->p[0])
764 *y=0.0;
765 else if(x<=r->p[0]+r->p[1]) {
766 f=exp(-r->p[3]*(x-r->p[0]));
767 *y=(r->p[2]/(r->p[3]*r->p[3]))*(1.0-f);
768 } else {
769 f=exp(-r->p[3]*r->p[1]);
770 f+=exp(-r->p[3]*(x-r->p[0]-r->p[1]));
771 f-=2.0*exp(-r->p[3]*(x-r->p[0]));
772 *y=(r->p[2]/(r->p[3]*r->p[3]))*f;
773 }
774 break;
775 case 334:
776 if(x<=r->p[0])
777 *y=0.0;
778 else if(x>r->p[0] && x<=r->p[1]) {
779 *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])));
780 } else {
781 *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])));
782 }
783 break;
784 case 351:
785 f=1.0-r->p[0]*(2.0-exp(-r->p[1]*x)-exp(-r->p[2]*x)); *y=f;
786 break;
787 case 403:
788 xt=x-r->p[4]; if(xt<0.0) xt=0.0;
789 f=r->p[0]*(1.0 - r->p[1]*igam(r->p[3], r->p[2]*xt));
790 *y=f;
791 break;
792 case 831:
793 xt=x-r->p[0]; if(xt<=0.0) {
794 f=r->p[1];
795 } else {
796 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);
797 f*=r->p[1];
798 }
799 *y=f;
800 break;
801 case 841:
802 f=r->p[0]*pow(x, r->p[1]) / (pow(x, r->p[1]) + r->p[2]);
803 *y=f;
804 break;
805 case 842:
806 f=1.0-r->p[0]*pow(x, r->p[1]) / (pow(x, r->p[1]) + r->p[2]);
807 *y=f;
808 break;
809 case 843:
810 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]);
811 *y=f;
812 break;
813 case 844:
814 if(r->parNr>4) xt=x-r->p[4]; else xt=x; // delay time accounted, if avail
815 if(xt<=0.0) f=r->p[3];
816 else {a=pow(xt, r->p[1]); f=r->p[0]*a/(a + r->p[2]) + r->p[3];}
817 *y=f;
818 break;
819 case 845:
820 f=r->p[0] - r->p[0]*pow(x, r->p[1]) / (pow(x, r->p[1]) + r->p[2]);
821 *y=f;
822 break;
823 case 846:
824 xt=x-r->p[4]; if(xt<=0.0) {
825 f=r->p[3];
826 } else {
827 a=pow(xt, r->p[1]);
828 f=r->p[3]+(r->p[0]-r->p[3])*a/(a+r->p[2]);
829 }
830 *y=f;
831 break;
832 case 847:
833 xt=x-r->p[4]; if(xt<=0.0) {
834 f=1.0-r->p[3];
835 } else {
836 a=pow(xt, r->p[1]);
837 f=1.0-r->p[3]-(r->p[0]-r->p[3])*a/(a+r->p[2]);
838 }
839 *y=f;
840 break;
841 case 848:
842 xt=x-r->p[4]; if(xt<=0.0) {
843 f=r->p[3];
844 } else {
845 a=pow(xt, r->p[1]);
846 f=r->p[3] * (1.0 - (1.0-r->p[0])*a/(r->p[2]+a));
847 }
848 *y=f;
849 break;
850 case 849:
851 xt=x-r->p[4]; if(xt<=0.0) {
852 f=1.0-r->p[3];
853 } else {
854 a=pow(xt, r->p[1]);
855 f=1.0 - r->p[3] * (1.0 - (1.0-r->p[0])*a/(r->p[2]+a));
856 }
857 *y=f;
858 break;
859 case 851:
860 a=r->p[0]*x; a*=a; f=1.0/pow(1.0+a, r->p[1]);
861 *y=f;
862 break;
863 case 852:
864 a=r->p[0]*x; a*=a; f=1.0 - 1.0/pow(1.0+a, r->p[1]);
865 *y=f;
866 break;
867 case 861:
868 xt=x-r->p[3]; if(xt<=0.0) {
869 f=1.0;
870 } else {
871 f=pow(1.0+pow(r->p[0]*xt, r->p[1]), -r->p[2]);
872 }
873 *y=f;
874 break;
875 case 862:
876 xt=x-r->p[3]; if(xt<=0.0) {
877 f=1.0;
878 } else {
879 f=pow(1.0+pow(r->p[0]*xt, r->p[1]), -r->p[2]);
880 }
881 *y=1.0-f;
882 break;
883 case 863:
884 xt=x-r->p[4]; if(xt<=0.0) {
885 f=r->p[3];
886 } else {
887 f=pow(pow(r->p[3],-1.0/r->p[2])+pow(r->p[0]*xt, r->p[1]), -r->p[2]);
888 }
889 *y=f;
890 break;
891 case 864:
892 xt=x-r->p[4]; if(xt<=0.0) {
893 f=1.0-r->p[3];
894 } else {
895 f=1.0-pow(pow(r->p[3],-1.0/r->p[2])+pow(r->p[0]*xt, r->p[1]), -r->p[2]);
896 }
897 *y=f;
898 break;
899 case 871:
900 case 872:
901 case 873:
902 case 874:
903 for(m=0; m<3; m++) for(j=0; j<5; j++) mp[m][j]=0.0;
904 for(i=m=j=0; i<r->parNr; i++) {mp[m][j]=r->p[i]; j++; if(j>4) {j=0; m++;}}
905 n=r->parNr/5;
906 for(m=0; m<n; m++) {
907 xt=x-mp[m][4];
908 if(xt<=0.0) mf[m]=1.0-mp[m][3];
909 else {
910 f=pow(xt, mp[m][1]);
911 mf[m]=1.0-(mp[m][3]+(mp[m][0]-mp[m][3])*f/(mp[m][2]+f));
912 }
913 }
914 a=1.0-mf[0]*mf[1]-mf[0]*mf[2]-mf[1]*mf[2]+2.0*mf[0]*mf[1]*mf[2];
915 if(r->type==871) {
916 f=1.0;
917 f-=mf[0]*(1.0-mf[1]-mf[2]+mf[1]*mf[2])/a;
918 f-=mf[1]*(1.0-mf[0]-mf[2]+mf[0]*mf[2])/a;
919 f-=mf[2]*(1.0-mf[0]-mf[1]+mf[0]*mf[1])/a;
920 *y=f;
921 } else if(r->type==872) {
922 *y=mf[0]*(1.0-mf[1]-mf[2]+mf[1]*mf[2])/a;
923 } else if(r->type==873) {
924 *y=mf[1]*(1.0-mf[0]-mf[2]+mf[0]*mf[2])/a;
925 } else if(r->type==874) {
926 *y=mf[2]*(1.0-mf[0]-mf[1]+mf[0]*mf[1])/a;
927 }
928 break;
929 case 881:
930 case 882:
931 case 883:
932 case 884:
933 for(m=0; m<3; m++) for(j=0; j<5; j++) mp[m][j]=0.0;
934 for(i=m=j=0; i<r->parNr; i++) {mp[m][j]=r->p[i]; j++; if(j>4) {j=0; m++;}}
935 n=r->parNr/5;
936 for(m=0; m<n; m++) {
937 xt=x-mp[m][4];
938 if(xt<=0.0) mf[m]=1.0-mp[m][3];
939 else {
940 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]);
941 }
942 }
943 a=1.0-mf[0]*mf[1]-mf[0]*mf[2]-mf[1]*mf[2]+2.0*mf[0]*mf[1]*mf[2];
944 if(r->type==881) {
945 f=1.0;
946 f-=mf[0]*(1.0-mf[1]-mf[2]+mf[1]*mf[2])/a;
947 f-=mf[1]*(1.0-mf[0]-mf[2]+mf[0]*mf[2])/a;
948 f-=mf[2]*(1.0-mf[0]-mf[1]+mf[0]*mf[1])/a;
949 *y=f;
950 } else if(r->type==882) {
951 *y=mf[0]*(1.0-mf[1]-mf[2]+mf[1]*mf[2])/a;
952 } else if(r->type==883) {
953 *y=mf[1]*(1.0-mf[0]-mf[2]+mf[0]*mf[2])/a;
954 } else if(r->type==884) {
955 *y=mf[2]*(1.0-mf[0]-mf[1]+mf[0]*mf[1])/a;
956 }
957 break;
958 case 1010: // step functions
959 /* Parameter number must be an even number and >=2 */
960 if(r->parNr<2 || r->parNr&1) return(2);
961 *y=0.0;
962 for(int i=0; i<r->parNr; i+=2) if(x>=r->p[i]) *y=r->p[i+1];
963 break;
964 case 2111:
965 xt=x-r->p[3]; a=sqrt(2.0)*(r->p[2]/2.355);
966 f=r->p[4]+(r->p[0]/2.0)*(erf((xt+r->p[1])/a)-erf((xt-r->p[1])/a));
967 *y=f;
968 break;
969 case 1232:
970 xt=x; if(r->parNr>7) xt-=r->p[7];
971 if(xt<=0.0) {
972 f=0.0;
973 } else {
974 sqrx=xt*xt; cubx=sqrx*xt;
975 a=r->p[0]+r->p[2]*xt+r->p[4]*sqrx+r->p[6]*cubx;
976 b=r->p[1]+r->p[3]*xt+r->p[5]*sqrx; if(b==0.0) break;
977 f=a/b;
978 }
979 *y=f;
980 break;
981 case 1312:
982 xt=x; if(r->parNr>4) xt-=r->p[4];
983 if(xt<=0.0) {
984 f=0.0;
985 } else {
986 f=0.0;
987 a=(r->p[0]*(xt)-r->p[2])*exp(r->p[1]*(xt)); f+=a;
988 a=r->p[2]*exp(r->p[3]*(xt)); f+=a;
989 }
990 *y=f;
991 break;
992 case 1313:
993 xt=x; if(r->parNr>6) xt-=r->p[6];
994 if(xt<=0.0) {
995 f=0.0;
996 } else {
997 f=0.0;
998 a=(r->p[0]*(xt)-r->p[2]-r->p[4])*exp(r->p[1]*(xt)); f+=a;
999 a=r->p[2]*exp(r->p[3]*(xt)); f+=a;
1000 a=r->p[4]*exp(r->p[5]*(xt)); f+=a;
1001 }
1002 *y=f;
1003 break;
1004 case 1314:
1005 xt=x; if(r->parNr>8) xt-=r->p[8];
1006 if(xt<=0.0) {
1007 f=0.0;
1008 } else {
1009 f=0.0;
1010 a=(r->p[0]*(xt)-r->p[2]-r->p[4]-r->p[6])*exp(r->p[1]*(xt)); f+=a;
1011 a=r->p[2]*exp(r->p[3]*(xt)); f+=a;
1012 a=r->p[4]*exp(r->p[5]*(xt)); f+=a;
1013 a=r->p[6]*exp(r->p[7]*(xt)); f+=a;
1014 }
1015 *y=f;
1016 break;
1017 case 1401:
1018 xt=x-r->p[3]; if(xt<=0.0 || r->p[2]==0.0) {
1019 f=0.0;
1020 } else {
1021 f=r->p[0]*pow(xt, r->p[1])*exp(-xt/r->p[2]);
1022 }
1023 *y=f;
1024 break;
1025 case 1402:
1026 xt=x-r->p[3]; if(xt<=0.0 || r->p[2]==0.0) {
1027 f=r->p[4];
1028 } else {
1029 f=r->p[0]*pow(xt, r->p[1])*exp(-xt/r->p[2]) + r->p[4];
1030 }
1031 *y=f;
1032 break;
1033 case 1403:
1034 xt=x-r->p[0];
1035 f=0.0;
1036 if(xt>0.0) {
1037 a=exp(-xt/r->p[3]);
1038 if(r->p[1]>0.0) f += r->p[1]*pow(xt, r->p[2])*a;
1039 if(r->parNr==6 && r->p[4]>0.0) f += r->p[4]*(1.0-a)*exp(-xt/r->p[5]);
1040 }
1041 *y=f;
1042 break;
1043 case 1421:
1044 xt=x-r->p[0]; if(xt<=0.0) {
1045 f=0.0;
1046 } else {
1047 f=r->p[1]*(1.0-exp(-pow(xt/r->p[2], r->p[3])));
1048 }
1049 *y=f;
1050 break;
1051 case 1423:
1052 xt=x-r->p[0]; if(xt<=0.0) {
1053 f=0.0;
1054 } else {
1055 a=xt/r->p[2]; b=pow(a, r->p[3]-1.0); f=exp(-b*a);
1056 a=r->p[3]*b*f/r->p[2]; b=1.0-f;
1057 f=r->p[1]*(a+r->p[4]*b);
1058 }
1059 *y=f;
1060 break;
1061 case 1430:
1062 if(r->parNr<3) {*y=nan(""); break;}
1063 f=0.0;
1064 xt=x-r->p[0];
1065 if(xt>0.0)
1066 for(i=2; i<r->parNr; i+=2)
1067 f+=r->p[i-1]*xt*exp(-r->p[i]*xt);
1068 *y=f;
1069 break;
1070 case 1431:
1071 if(x<=0.0) {
1072 f=0.0;
1073 } else {
1074 f=r->p[0]*x*exp(-r->p[1]*x)*r->p[1]*r->p[1];
1075 }
1076 *y=f;
1077 break;
1078 case 1432:
1079 if(x<=0.0) {
1080 f=0.0;
1081 } else {
1082 f=r->p[0]*x*exp(-r->p[1]*x);
1083 }
1084 *y=f;
1085 break;
1086 case 1433:
1087 if(x<=0.0) {
1088 f=0.0;
1089 } else {
1090 double e=exp(-r->p[1]*x);
1091 f=x*e + (r->p[2]/(r->p[1]*r->p[1]))*(1.0-(r->p[1]*x+1.0)*e);
1092 f*=r->p[0];
1093 }
1094 *y=f;
1095 break;
1096 case 1434:
1097 if(x<=0.0) {
1098 f=1.0/(1.0-r->p[0]);
1099 } else {
1100 double e=exp(-r->p[2]*x);
1101 double rcp=x*e + (r->p[3]/(r->p[2]*r->p[2]))*(1.0-(r->p[2]*x+1.0)*e);
1102 rcp*=r->p[1];
1103 f=1.0/(1.0-r->p[0]*(1.0-rcp));
1104 }
1105 *y=f;
1106 break;
1107 case 1435:
1108 xt=x-r->p[0]; if(xt<=0.0) {
1109 f=0.0;
1110 } else {
1111 f=r->p[2]*(r->p[3]*xt*exp(-r->p[4]*r->p[1]*xt) + exp(-r->p[1]*xt) );
1112 }
1113 *y=f;
1114 break;
1115 case 1441: /* Probability density function of Erlang distribution */
1116 f=0.0;
1117 if(x>=0.0) {
1118 double A=r->p[0], k=r->p[1]; if(!(k>=0.0)) return(2);
1119 int N=(int)round(r->p[2]); if(!(N>0)) return(2);
1120 unsigned long long int f=lfactorial(N-1);
1121 f = A*pow(k, N)*pow(x, N-1)*exp(-k*x)/f;
1122 }
1123 *y=f;
1124 break;
1125 case 1801:
1126 xt=x-r->p[0]; if(xt<=0.0) {
1127 f=0.0;
1128 } else {
1129 f=r->p[1]*pow(xt, r->p[3])/(pow(r->p[2], r->p[3])+pow(xt, r->p[3]));
1130 }
1131 *y=f;
1132 break;
1133 case 1811:
1134 xt=x-r->p[0]; if(xt<=0.0) {
1135 f=0.0;
1136 } else {
1137 a=pow(r->p[2], r->p[3])+pow(xt, r->p[3]);
1138 f=r->p[1]*r->p[3]*
1139 (pow(xt, r->p[3]-1.0)/a - pow(xt, 2.0*r->p[3]-1.0)/(a*a));
1140 }
1141 *y=f;
1142 break;
1143 case 1821:
1144 xt=x-r->p[0]; if(xt<=0.0) {
1145 f=0.0;
1146 } else {
1147 double a, c, b, k, xb, cb, cbxb, hill, hilld;
1148 a=r->p[1]; c=r->p[2]; b=r->p[3]; k=r->p[4];
1149 cb=pow(c, b); xb=pow(xt, b); cbxb=cb+xb;
1150 hilld= b*pow(xt, b-1.0)/cbxb - b*pow(xt, 2.0*b-1.0)/(cbxb*cbxb);
1151 hill=k*xb/cbxb;
1152 f=a*(hilld+hill);
1153 }
1154 *y=f;
1155 break;
1156 case 1833:
1157 xt=x; if(r->parNr>3) xt-=r->p[3];
1158 if(xt<=0.0) {
1159 f=0.0;
1160 } else {
1161 double c=0.0; if(r->parNr>2) c=r->p[2];
1162 double e=exp(-r->p[1]*xt);
1163 f=xt*e + (c/(r->p[1]*r->p[1]))*(1.0-(r->p[1]*xt+1.0)*e);
1164 f*=r->p[0];
1165 }
1166 *y=f;
1167 break;
1168 case 2233:
1169 if(x<=0.0) {
1170 f=1.0/(1.0-r->p[0]);
1171 } else {
1172 sqrx=x*x; cubx=sqrx*x;
1173 a=0.0+r->p[1]*x+r->p[3]*sqrx+r->p[5]*cubx;
1174 b=1.0+r->p[2]*x+r->p[4]*sqrx+r->p[6]*cubx;
1175 f=1.0/(1.0-r->p[0]*(1.0-(a/b)));
1176 }
1177 *y=f;
1178 break;
1179 case 2313:
1180 if(r->parNr<3) return(2);
1181 if(x<=0.0) {
1182 f=0.0;
1183 } else {
1184 double A1=r->p[1];
1185 double L1=r->p[2];
1186 double A2=0.0; if(r->parNr>3) A2=r->p[3];
1187 double L2=0.0; if(r->parNr>4) L2=r->p[4];
1188 double A3=0.0; if(r->parNr>5) A3=r->p[5];
1189 double L3=0.0; if(r->parNr>6) L3=r->p[6];
1190 f=(A1*x - A2 - A3)*exp(-L1*x) + A2*exp(-L2*x) + A3*exp(-L3*x);
1191 }
1192 *y=1.0/(1.0-r->p[0]*(1.0-f));
1193 break;
1194 case 2801:
1195 {
1196 double top, bottom, ec50, hillslope;
1197 top=r->p[0]; bottom=r->p[1]; ec50=r->p[2]; hillslope=r->p[3];
1198 if(x<=0.0) f=bottom;
1199 else f=bottom+(top-bottom)/(1.0+pow(ec50/x,hillslope));
1200 }
1201 *y=f;
1202 break;
1203 case 2802:
1204 {
1205 double top, bottom, logec50, hillslope;
1206 top=r->p[0]; bottom=r->p[1]; logec50=r->p[2]; hillslope=r->p[3];
1207 f=bottom+(top-bottom)/(1.0+pow(10.0, (logec50-x)*hillslope));
1208 }
1209 *y=f;
1210 break;
1211 case 2841:
1212 a=r->p[1]*pow(x, r->p[2]) / (pow(x, r->p[2]) + r->p[3]);
1213 f=1.0/(1.0-r->p[0]*(1.0-a));
1214 *y=f;
1215 break;
1216 case 3331: *y=nan(""); break; /* cannot be applied to one point only */
1217 case 9501: *y=nan(""); break; /* cannot be applied to one point only */
1218 case 9502: *y=nan(""); break; /* cannot be applied to one point only */
1219 case 9503: *y=nan(""); break; /* cannot be applied to one point only */
1220 case 9701: *y=nan(""); break; /* cannot be applied to one point only */
1221 default:
1222 *y=nan("");
1223 }
1224 if(isnan(*y)) return 1;
1225 return 0;
1226}
1227/*****************************************************************************/
1228
1229/*****************************************************************************/
1234int fitEvaltac(
1236 FitVOI *r,
1238 double *x,
1240 double *y,
1242 int dataNr
1243) {
1244 if(r==NULL || x==NULL || y==NULL || dataNr<1) return(1);
1245
1246 /* Special cases that can only be computed as TACs */
1247 if(r->type==3331) {
1248 double Ta=r->p[0];
1249 double Ti=r->p[1];
1250 double dT=r->p[r->parNr-2]; Ta+=dT;
1251 double tau=r->p[r->parNr-1];
1252 int n=(r->parNr-4)/2; if(n<1) return(2);
1253 for(int i=0; i<dataNr; i++) {
1254 if(x[i]<=Ta) y[i]=0.0;
1255 else if(x[i]<Ta+Ti) {
1256 double f=0.0;
1257 for(int j=0; j<n; j++) {
1258 double b;
1259 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];
1260 f+=b*(1.0-exp(-r->p[2*j+3]*(x[i]-Ta)));
1261 }
1262 if(Ti>0.0) f*=(1.0/Ti);
1263 y[i]=f;
1264 } else {
1265 double f=0.0;
1266 for(int j=0; j<n; j++) {
1267 double b;
1268 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];
1269 f+=b*(exp(-r->p[2*j+3]*(x[i]-Ta-Ti)) - exp(-r->p[2*j+3]*(x[i]-Ta)));
1270 }
1271 if(Ti>0.0) f*=(1.0/Ti);
1272 y[i]=f;
1273 }
1274 }
1275 if(simDispersion(x, y, dataNr, tau, 0.0, NULL)) return(2);
1276 return(0);
1277 }
1278
1279 /* Usual functions */
1280 for(int i=0; i<dataNr; i++) if(fitEval(r, x[i], y+i)) return(2);
1281 return 0;
1282}
1283/*****************************************************************************/
1284
1285/*****************************************************************************/
1290int fitIntegralEval(
1292 FitVOI *r,
1294 double x,
1296 double *yi
1297) {
1298 double a, f, xt, t;
1299
1300 if(r==NULL || yi==NULL) return 1;
1301 *yi=nan("");
1302 switch(r->type) {
1303 case 301:
1304 if(fabs(r->p[1])>1.0e-12) f=(r->p[0]/r->p[1])*(exp(r->p[1]*x)-1.0);
1305 else f=r->p[0]*x;
1306 *yi=f;
1307 break;
1308 case 302:
1309 f=0;
1310 if(fabs(r->p[1])>1.0e-12) a=(r->p[0]/r->p[1])*(exp(r->p[1]*x)-1.0);
1311 else a=r->p[0]*x;
1312 f+=a;
1313 if(fabs(r->p[3])>1.0e-12) a=(r->p[2]/r->p[3])*(exp(r->p[3]*x)-1.0);
1314 else a=r->p[2]*x;
1315 f+=a;
1316 *yi=f;
1317 break;
1318 case 303:
1319 f=0;
1320 if(fabs(r->p[1])>1.0e-12) a=(r->p[0]/r->p[1])*(exp(r->p[1]*x)-1.0);
1321 else a=r->p[0]*x;
1322 f+=a;
1323 if(fabs(r->p[3])>1.0e-12) a=(r->p[2]/r->p[3])*(exp(r->p[3]*x)-1.0);
1324 else a=r->p[2]*x;
1325 f+=a;
1326 if(fabs(r->p[5])>1.0e-12) a=(r->p[4]/r->p[5])*(exp(r->p[5]*x)-1.0);
1327 else a=r->p[4]*x;
1328 f+=a;
1329 *yi=f;
1330 break;
1331 case 304:
1332 f=0;
1333 if(fabs(r->p[1])>1.0e-12) a=(r->p[0]/r->p[1])*(exp(r->p[1]*x)-1.0);
1334 else a=r->p[0]*x;
1335 f+=a;
1336 if(fabs(r->p[3])>1.0e-12) a=(r->p[2]/r->p[3])*(exp(r->p[3]*x)-1.0);
1337 else a=r->p[2]*x;
1338 f+=a;
1339 if(fabs(r->p[5])>1.0e-12) a=(r->p[4]/r->p[5])*(exp(r->p[5]*x)-1.0);
1340 else a=r->p[4]*x;
1341 f+=a;
1342 if(fabs(r->p[7])>1.0e-12) a=(r->p[6]/r->p[7])*(exp(r->p[7]*x)-1.0);
1343 else a=r->p[6]*x;
1344 f+=a;
1345 *yi=f;
1346 break;
1347 case 305:
1348 f=0;
1349 if(fabs(r->p[1])>1.0e-12) a=(r->p[0]/r->p[1])*(exp(r->p[1]*x)-1.0);
1350 else a=r->p[0]*x;
1351 f+=a;
1352 if(fabs(r->p[3])>1.0e-12) a=(r->p[2]/r->p[3])*(exp(r->p[3]*x)-1.0);
1353 else a=r->p[2]*x;
1354 f+=a;
1355 if(fabs(r->p[5])>1.0e-12) a=(r->p[4]/r->p[5])*(exp(r->p[5]*x)-1.0);
1356 else a=r->p[4]*x;
1357 f+=a;
1358 if(fabs(r->p[7])>1.0e-12) a=(r->p[6]/r->p[7])*(exp(r->p[7]*x)-1.0);
1359 else a=r->p[6]*x;
1360 f+=a;
1361 if(fabs(r->p[9])>1.0e-12) a=(r->p[8]/r->p[9])*(exp(r->p[8]*x)-1.0);
1362 else a=r->p[8]*x;
1363 f+=a;
1364 *yi=f;
1365 break;
1366 case 321:
1367 f=(r->p[0]/r->p[1])*exp(r->p[1]*x)
1368 - r->p[0]*exp((r->p[1]+r->p[2])*x)/(r->p[1]+r->p[2]);
1369 *yi=f;
1370 break;
1371 case 322:
1372 f=0.0;
1373 a=(r->p[0]/r->p[1])*exp(r->p[1]*x)
1374 - r->p[0]*exp((r->p[1]+r->p[2])*x)/(r->p[1]+r->p[2]); f+=a;
1375 a=(r->p[3]/r->p[4])*exp(r->p[4]*x)
1376 - r->p[3]*exp((r->p[4]+r->p[5])*x)/(r->p[4]+r->p[5]); f+=a;
1377 *yi=f;
1378 break;
1379 case 323:
1380 f=0.0;
1381 a=(r->p[0]/r->p[1])*exp(r->p[1]*x)
1382 - r->p[0]*exp((r->p[1]+r->p[2])*x)/(r->p[1]+r->p[2]); f+=a;
1383 a=(r->p[3]/r->p[4])*exp(r->p[4]*x)
1384 - r->p[3]*exp((r->p[4]+r->p[5])*x)/(r->p[4]+r->p[5]); f+=a;
1385 a=(r->p[6]/r->p[7])*exp(r->p[7]*x)
1386 - r->p[6]*exp((r->p[7]+r->p[8])*x)/(r->p[7]+r->p[8]); f+=a;
1387 *yi=f;
1388 break;
1389 case 331: break; /* not yet applied */
1390 case 351: break; /* not yet applied */
1391 case 1312:
1392 t=r->p[4]; xt=x-t; f=0.0;
1393 if(xt>0.0) {
1394 double A1, A2, L1, L2;
1395 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3];
1396 if(L1!=0.0) {
1397 a=exp(L1*xt);
1398 f+=(A2/L1)*(1.0 - a);
1399 f+=(A1/(L1*L1))*(1.0 + a*(L1*xt - 1.0));
1400 }
1401 if(L2!=0.0) f+=(A2/L2)*(exp(L2*xt) - 1.0); else f+=A2*xt;
1402 }
1403 *yi=f;
1404 break;
1405 case 1313:
1406 t=r->p[6]; xt=x-t; f=0.0;
1407 if(xt>0.0) {
1408 double A1, A2, A3, L1, L2, L3;
1409 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3]; A3=r->p[4]; L3=r->p[5];
1410 if(L1!=0.0) {
1411 a=exp(L1*xt);
1412 f+=(A2/L1)*(1.0 - a);
1413 f+=(A3/L1)*(1.0 - a);
1414 f+=(A1/(L1*L1))*(1.0 + a*(L1*xt - 1.0));
1415 }
1416 if(L2!=0.0) f+=(A2/L2)*(exp(L2*xt) - 1.0); else f+=A2*xt;
1417 if(L3!=0.0) f+=(A3/L3)*(exp(L3*xt) - 1.0); else f+=A3*xt;
1418 }
1419 *yi=f;
1420 break;
1421 case 1314:
1422 t=r->p[8]; xt=x-t; f=0.0;
1423 if(xt>0.0) {
1424 double A1, A2, A3, A4, L1, L2, L3, L4;
1425 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3]; A3=r->p[4]; L3=r->p[5];
1426 A4=r->p[6]; L4=r->p[7];
1427 if(L1!=0.0) {
1428 a=exp(L1*xt);
1429 f+=(A2/L1)*(1.0 - a);
1430 f+=(A3/L1)*(1.0 - a);
1431 f+=(A4/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 if(L4!=0.0) f+=(A4/L4)*(exp(L4*xt) - 1.0); else f+=A4*xt;
1437 }
1438 *yi=f;
1439 break;
1440 case 1401: break; /* not yet applied */
1441 case 1402: break; /* not yet applied */
1442 case 1430:
1443 if(r->parNr<3) break;
1444 xt=x-r->p[0];
1445 f=0.0;
1446 if(xt>0.0) {
1447 for(int pi=2; pi<r->parNr; pi+=2) {
1448 double a, b;
1449 a=r->p[pi-1]; b=r->p[pi];
1450 f += a*(1.0-(b*xt+1.0)*exp(-b*xt))/(b*b);
1451 }
1452 }
1453 *yi=f;
1454 break;
1455 case 1431:
1456 f=0.0;
1457 if(x>0.0) f=r->p[0]*(1.0 - (r->p[1]*x + 1.0)*exp(-r->p[1]*x));
1458 *yi=f;
1459 break;
1460 case 1432:
1461 f=0.0;
1462 if(x>0.0) {
1463 a=r->p[0]/(r->p[1]*r->p[1]);
1464 f=a*(1.0 - (r->p[1]*x + 1.0)*exp(-r->p[1]*x));
1465 }
1466 *yi=f;
1467 break;
1468 case 1435:
1469 xt=x-r->p[0]; f=0.0;
1470 if(xt>0.0) {
1471 double a, b, m, n; a=r->p[1]; b=r->p[2]; m=r->p[3]; n=r->p[4];
1472 f+=m*((1.0-(n*a*xt+1.0)*exp(-n*a*xt)))/(n*n*a*a) + (1.0-exp(-a*xt))/a;
1473 f*=b;
1474 }
1475 *yi=f;
1476 break;
1477 case 1441: /* Probability density function of Erlang distribution, or actually,
1478 its integral is the cumulative distribution function of Erlang distribution */
1479 f=0.0;
1480 if(x>=0.0) {
1481 double A, k; A=r->p[0]; k=r->p[1]; if(!(k>=0.0)) return(2);
1482 int N=(int)round(r->p[2]); if(!(N>=0) || N>20) return(2); // Supports only N=0,1,2,..,20
1483 if(N==0) f=A;
1484 else if(N==1) f=A*(1.0-exp(-k*x));
1485 else if(N==2) f=A*(1.0 - exp(-k*x) - k*x*exp(-k*x));
1486 else {
1487 double s=1.0+k*x;
1488 for(int j=2; j<N; j++) s+=pow(k*x, j)/lfactorial(j);
1489 f=A*(1.0-s*exp(-k*x));
1490 }
1491 }
1492 *yi=f;
1493 break;
1494 case 2111: *yi=nan(""); break; /* no application */
1495 case 9501: *yi=nan(""); break; /* cannot be applied to one point only */
1496 case 9502: *yi=nan(""); break; /* cannot be applied to one point only */
1497 case 9503: *yi=nan(""); break; /* cannot be applied to one point only */
1498 case 9701: *yi=nan(""); break; /* cannot be applied to one point only */
1499 default:
1500 *yi=nan("");
1501 }
1502 if(isnan(*yi)) return(3);
1503 return(0);
1504}
1505/*****************************************************************************/
1506
1507/*****************************************************************************/
1512int fitIntegralEvaltac(
1514 FitVOI *r,
1516 double *x,
1518 double *yi,
1520 int dataNr
1521) {
1522 int i;
1523
1524 if(r==NULL || x==NULL || yi==NULL || dataNr<1) return(1);
1525 for(i=0; i<dataNr; i++) if(fitIntegralEval(r, x[i], yi+i)) return(2);
1526 return(0);
1527}
1528/*****************************************************************************/
1529
1530/*****************************************************************************/
1536int fitEvaltacframes(
1538 FitVOI *r,
1540 double *x1,
1542 double *x2,
1544 double *y,
1546 int dataNr
1547) {
1548 if(r==NULL || x1==NULL || x2==NULL || y==NULL || dataNr<1) return(1);
1549 for(int i=0; i<dataNr; i++) {
1550 double fd=x2[i]-x1[i]; // frame duration
1551 if(!(fd>=0.0)) return(1);
1552 if(fd<1.0E-025) { // if frame duration is about zero, then simply calculate f(x) for that
1553 if(fitEval(r, x2[i], y+i)) return(2);
1554 continue;
1555 }
1556 double v1, v2;
1557 if(fitIntegralEval(r, x1[i], &v1) || fitIntegralEval(r, x2[i], &v2)) return(2);
1558 y[i]=(v2-v1)/fd;
1559 }
1560 return(0);
1561}
1562/*****************************************************************************/
1563
1564/*****************************************************************************/
1569int fitDerivEval(
1571 FitVOI *r,
1573 double x,
1575 double *yd
1576) {
1577 double a, f, xt, t;
1578
1579 if(r==NULL || yd==NULL) return 1;
1580 *yd=nan("");
1581 switch(r->type) {
1582 case 301:
1583 break;
1584 case 302:
1585 break;
1586 case 303:
1587 break;
1588 case 304:
1589 break;
1590 case 305:
1591 break;
1592 case 321:
1593 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);
1594 *yd=f;
1595 break;
1596 case 322:
1597 f=0.0;
1598 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;
1599 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;
1600 *yd=f;
1601 break;
1602 case 323:
1603 f=0.0;
1604 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;
1605 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;
1606 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;
1607 *yd=f;
1608 break;
1609 case 331: break; /* not yet applied */
1610 case 351: break; /* not yet applied */
1611 case 1312:
1612 t=r->p[4]; xt=x-t; f=0.0;
1613 if(xt>0.0) {
1614 double A1, A2, L1, L2;
1615 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3];
1616 f= A1*exp(L1*xt) + A2*L2*exp(L2*xt) + (A1*xt -A2)*L1*exp(L1*xt);
1617 }
1618 *yd=f;
1619 break;
1620 case 1313:
1621 t=r->p[6]; xt=x-t; f=0.0;
1622 if(xt>0.0) {
1623 double A1, A2, A3, L1, L2, L3;
1624 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3]; A3=r->p[4]; L3=r->p[5];
1625 f= A1*exp(L1*xt) + A2*L2*exp(L2*xt) + A3*L3*exp(L3*xt) +
1626 (A1*xt -A2 -A3)*L1*exp(L1*xt);
1627 }
1628 *yd=f;
1629 break;
1630 case 1314: break;
1631 t=r->p[8]; xt=x-t; f=0.0;
1632 if(xt>0.0) {
1633 double A1, A2, A3, A4, L1, L2, L3, L4;
1634 A1=r->p[0]; L1=r->p[1]; A2=r->p[2]; L2=r->p[3]; A3=r->p[4]; L3=r->p[5];
1635 A4=r->p[6]; L4=r->p[7];
1636 f= A1*exp(L1*xt) + A2*L2*exp(L2*xt) + A3*L3*exp(L3*xt) +
1637 A4*L4*exp(L4*xt) + (A1*xt -A2 -A3 -A4)*L1*exp(L1*xt);
1638 }
1639 *yd=f;
1640 break;
1641 case 1401:
1642 case 1402:
1643 xt=x-r->p[3]; if(xt<=0.0 || r->p[2]==0.0) {
1644 f=0.0;
1645 } else {
1646 f=r->p[0]*pow(xt, r->p[1]-1.0)*exp(-xt/r->p[2])*(r->p[1]-(xt/r->p[2]));
1647 }
1648 *yd=f;
1649 break;
1650 case 1421:
1651 xt=x-r->p[0]; if(xt<=0.0) {
1652 *yd=0.0;
1653 } else {
1654 a=pow(xt/r->p[2], r->p[3]-1.0);
1655 *yd=r->p[1]*r->p[3]*a*exp(-a*xt/r->p[2])/r->p[2];
1656 }
1657 break;
1658 case 1801:
1659 xt=x-r->p[0]; if(xt<=0.0) {
1660 *yd=0.0;
1661 } else {
1662 a=pow(r->p[2], r->p[3])+pow(xt, r->p[3]);
1663 *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));
1664 }
1665 break;
1666 case 2111: *yd=nan(""); break; /* no application */
1667 case 9501: *yd=nan(""); break; /* cannot be applied to one point only */
1668 case 9502: *yd=nan(""); break; /* cannot be applied to one point only */
1669 case 9503: *yd=nan(""); break; /* cannot be applied to one point only */
1670 case 9701: *yd=nan(""); break; /* cannot be applied to one point only */
1671 default:
1672 *yd=nan("");
1673 }
1674 if(isnan(*yd)) return(3);
1675 return(0);
1676}
1677/*****************************************************************************/
1678
1679/*****************************************************************************/
1684int fitDerivEvaltac(
1686 FitVOI *r,
1688 double *x,
1690 double *yd,
1692 int dataNr
1693) {
1694 int i;
1695
1696 if(r==NULL || x==NULL || yd==NULL || dataNr<1) return(1);
1697 for(i=0; i<dataNr; i++) if(fitDerivEval(r, x[i], yd+i)) return(2);
1698 return(0);
1699}
1700/*****************************************************************************/
1701
1702/*****************************************************************************/
1707unsigned int factorial(
1711 unsigned int n
1712) {
1713 if(n<1) return(1);
1714 if(n>12) return(0);
1715 return(n*factorial(n-1));
1716}
1717/*****************************************************************************/
1718
1719/*****************************************************************************/
1724unsigned long long int lfactorial(
1728 unsigned long long int n
1729) {
1730 if(n<1) return(1);
1731 if(n>20) return(0);
1732 return(n*lfactorial(n-1));
1733}
1734/*****************************************************************************/
1735
1736/*****************************************************************************/
1744double igam(
1746 double a,
1748 double x
1749) {
1750 /* Check parameters */
1751 if(x==0.0) return(0.0);
1752 if((x<0.0) || (a<=0.0)) return(nan(""));
1753
1754 if((x>1.0) && (x>a)) return(1.0 - igamc(a, x));
1755
1756 /* Left tail of incomplete Gamma function: x^a * e^-x * Sum(k=0,Inf)(x^k / Gamma(a+k+1)) */
1757
1758 double ans, ax, c, r;
1759
1760 /* Compute x**a * exp(-x) / Gamma(a) */
1761 ax=a*log(x) - x - lgamma(a);
1762 if(ax<-DBL_MAX_10_EXP) return(0.0); // underflow
1763 ax=exp(ax);
1764
1765 /* power series */
1766 r=a; c=1.0; ans=1.0;
1767 do {
1768 r+=1.0;
1769 c*=x/r;
1770 ans+=c;
1771 } while(c/ans > DBL_EPSILON);
1772 return(ans*ax/a);
1773}
1774/*****************************************************************************/
1775
1776/*****************************************************************************/
1784double igamc(
1786 double a,
1788 double x
1789) {
1790 /* Check parameters */
1791 if((x<0.0) || (a<=0.0)) return(nan(""));
1792
1793 if((x<1.0) || (x<a)) return(1.0-igam(a, x));
1794
1795 double ans, ax, c, yc, r, t, y, z;
1796 double pk, pkm1, pkm2, qk, qkm1, qkm2;
1797 double big=4.503599627370496E+015;
1798 double biginv=2.22044604925031308085E-016;
1799
1800 ax = a*log(x) - x - lgamma(a);
1801 if(ax < -DBL_MAX_10_EXP) return(0.0); // underflow
1802 ax=exp(ax);
1803
1804 /* continued fraction */
1805 y=1.0-a; z=x+y+1.0; c=0.0;
1806 pkm2=1.0; qkm2=x; pkm1=x+1.0; qkm1=z*x;
1807 ans=pkm1/qkm1;
1808 do {
1809 c+=1.0; y+=1.0; z+=2.0;
1810 yc=y*c; pk=pkm1*z - pkm2*yc; qk=qkm1*z - qkm2*yc;
1811 if(qk!=0.0) {r=pk/qk; t=fabs((ans-r)/r); ans=r;}
1812 else t=1.0;
1813 pkm2=pkm1; pkm1=pk; qkm2=qkm1; qkm1=qk;
1814 if(fabs(pk)>big) {pkm2*=biginv; pkm1*=biginv; qkm2*=biginv; qkm1*=biginv;}
1815 } while(t>DBL_EPSILON);
1816 return(ans*ax);
1817}
1818/*****************************************************************************/
1819
1820/*****************************************************************************/
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.
MATHFUNC_TEST
int MATHFUNC_TEST
Definition mathfunc.c:5
fiterrmsg
char fiterrmsg[64]
Definition mathfunc.c:6
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:1744
lfactorial
unsigned long long int lfactorial(unsigned long long int n)
Definition mathfunc.c:1724
fitDerivEval
int fitDerivEval(FitVOI *r, double x, double *yd)
Definition mathfunc.c:1569
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:1234
fitEval
int fitEval(FitVOI *r, double x, double *y)
Definition mathfunc.c:618
igamc
double igamc(double a, double x)
Definition mathfunc.c:1784
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:1512
fitIntegralEval
int fitIntegralEval(FitVOI *r, double x, double *yi)
Definition mathfunc.c:1290
fitFunctionformat
int fitFunctionformat(int type, char *str)
Definition mathfunc.c:351
fitInit
void fitInit(FIT *fit)
Definition mathfunc.c:38
fitDerivEvaltac
int fitDerivEvaltac(FitVOI *r, double *x, double *yd, int dataNr)
Definition mathfunc.c:1684
factorial
unsigned int factorial(unsigned int n)
Definition mathfunc.c:1707
fitEvaltacframes
int fitEvaltacframes(FitVOI *r, double *x1, double *x2, double *y, int dataNr)
Definition mathfunc.c:1536
fitFunctionname
int fitFunctionname(int type, char *str)
Definition mathfunc.c:511
fitPrint
void fitPrint(FIT *fit)
Definition mathfunc.c:180
FIT::timeunit
int timeunit
Definition libtpccurveio.h:320
FIT::voiNr
int voiNr
Definition libtpccurveio.h:312
FIT::_voidataNr
int _voidataNr
Definition libtpccurveio.h:328
FIT::datafile
char datafile[FILENAME_MAX]
Definition libtpccurveio.h:314
FIT::program
char program[1024]
Definition libtpccurveio.h:326
FIT::unit
char unit[MAX_UNITS_LEN+1]
Definition libtpccurveio.h:318
FIT::studynr
char studynr[MAX_STUDYNR_LEN+1]
Definition libtpccurveio.h:316
FIT::voi
FitVOI * voi
Definition libtpccurveio.h:322
FIT::time
time_t time
Definition libtpccurveio.h:324
FitVOI::wss
double wss
Definition libtpccurveio.h:300
FitVOI::end
double end
Definition libtpccurveio.h:294
FitVOI::place
char place[MAX_REGIONSUBNAME_LEN+1]
Definition libtpccurveio.h:286
FitVOI::hemisphere
char hemisphere[MAX_REGIONSUBNAME_LEN+1]
Definition libtpccurveio.h:284
FitVOI::voiname
char voiname[MAX_REGIONSUBNAME_LEN+1]
Definition libtpccurveio.h:282
FitVOI::name
char name[MAX_REGIONNAME_LEN+1]
Definition libtpccurveio.h:280
FitVOI::start
double start
Definition libtpccurveio.h:292
FitVOI::p
double p[MAX_FITPARAMS]
Definition libtpccurveio.h:298
Generated on for TPCCLIB by doxygen 1.15.0