TPCCLIB
Loading...
Searching...
No Matches
Functions | Variables
mathfunc.c File Reference

IO for FIT files and calculating function values. More...

#include "libtpccurveio.h"
#include <unistd.h>
#include "libtpcmodel.h"

Go to the source code of this file.

Functions

void fitEmpty (FIT *fit)
 
void fitInit (FIT *fit)
 
int fitWrite (FIT *fit, char *filename)
 
int fitSetmem (FIT *fit, int voiNr)
 
void fitPrint (FIT *fit)
 
int fitRead (char *filename, FIT *fit, int verbose)
 
int fitFunctionformat (int type, char *str)
 
int fitFunctionname (int type, char *str)
 
int fitEval (FitVOI *r, double x, double *y)
 
int fitEvaltac (FitVOI *r, double *x, double *y, int dataNr)
 
int fitIntegralEval (FitVOI *r, double x, double *yi)
 
int fitIntegralEvaltac (FitVOI *r, double *x, double *yi, int dataNr)
 
int fitDerivEval (FitVOI *r, double x, double *yd)
 
int fitDerivEvaltac (FitVOI *r, double *x, double *yd, int dataNr)
 
unsigned int factorial (unsigned int n)
 
unsigned long long int lfactorial (unsigned long long int n)
 
double igam (double a, double x)
 
double igamc (double a, double x)
 

Variables

int MATHFUNC_TEST
 
char fiterrmsg [64]
 

Detailed Description

IO for FIT files and calculating function values.

Author
Vesa Oikonen

Definition in file mathfunc.c.

Function Documentation

◆ factorial()

unsigned int factorial ( unsigned int n)

Calculate factorial of given number.

See also
lfactorial, igam, igamc
Returns
Returns factorial n!, or zero if factorial would cause wrap-around.
Parameters
nInteger n, from which the factorial is calculated.
Note
Wrap-around will occur if n>12.

Definition at line 1678 of file mathfunc.c.

1683 {
1684 if(n<1) return(1);
1685 if(n>12) return(0);
1686 return(n*factorial(n-1));
1687}
factorial
unsigned int factorial(unsigned int n)
Definition mathfunc.c:1678

Referenced by factorial().

◆ fitDerivEval()

int fitDerivEval ( FitVOI * r,
double x,
double * yd )

Evaluates yd=Df(x).

See also
fitEval, fitIntegralEval
Returns
Returns 0, if ok.
Parameters
rFit parameters of a single region
xTime where to evaluate the derivative of the function
ydThe derivative of the function is returned here

Definition at line 1540 of file mathfunc.c.

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}
FitVOI::p
double p[MAX_FITPARAMS]
Definition libtpccurveio.h:297

Referenced by fitDerivEvaltac().

◆ fitDerivEvaltac()

int fitDerivEvaltac ( FitVOI * r,
double * x,
double * yd,
int dataNr )

Evaluates an array yd[i]=Df(x[i]).

See also
fitEvaltac, fitIntegralEvaltac, fitDerivEvaltac
Returns
Returns 0, if ok.
Parameters
rFit parameters of a single region
xTimes where to evaluate the function derivatives
ydArray for the function derivatives
dataNrNr of (x,yd) data

Definition at line 1655 of file mathfunc.c.

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}
fitDerivEval
int fitDerivEval(FitVOI *r, double x, double *yd)
Definition mathfunc.c:1540

◆ fitEmpty()

void fitEmpty ( FIT * fit)

Free memory allocated for FIT. All contents are cleared.

See also
fitInit, fitRead, fitWrite
Parameters
fitPointer to FIT struct.

Definition at line 18 of file mathfunc.c.

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}
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

Referenced by fitRead(), and fitSetmem().

◆ fitEval()

int fitEval ( FitVOI * r,
double x,
double * y )

Evaluate y=f(x).

See also
fitEvaltac, fitIntegralEval, fitDerivEval
Returns
Returns 0, if ok.
Parameters
rFit parameters of a single region
xTime where to evaluate the function
yThe value of the function is returned here

Definition at line 645 of file mathfunc.c.

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}
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
MATHFUNC_TEST
int MATHFUNC_TEST
Definition mathfunc.c:5

Referenced by fitEvaltac().

◆ fitEvaltac()

int fitEvaltac ( FitVOI * r,
double * x,
double * y,
int dataNr )

Evaluates an array y[i]=f(x[i]).

See also
fitEval, fitIntegralEvaltac, fitDerivEvaltac
Returns
Returns 0, if ok.
Parameters
rFit parameters of a single region
xTimes where to evaluate the function
yArray for the function values
dataNrNr of (x,y) data

Definition at line 1252 of file mathfunc.c.

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}
simDispersion
int simDispersion(double *x, double *y, int n, double tau1, double tau2, double *tmp)
Definition simulate.c:1911
fitEval
int fitEval(FitVOI *r, double x, double *y)
Definition mathfunc.c:645

◆ fitFunctionformat()

int fitFunctionformat ( int type,
char * str )

Copies the description of a function type to the specified string which must have space for >=128 characters.

Parameters
typeThe number of function
strRepresentation of the format of the function

Definition at line 381 of file mathfunc.c.

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}

◆ fitFunctionname()

int fitFunctionname ( int type,
char * str )

Copies the name of the function to the specified string which must have space for >=128 characters.

Parameters
typeThe number of function
strName of the function

Definition at line 539 of file mathfunc.c.

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}

◆ fitInit()

void fitInit ( FIT * fit)

Initiate FIT structure. Call this once before first use.

See also
fitEmpty, fitRead, fitWrite

Definition at line 38 of file mathfunc.c.

40 {
41 if(fit==NULL) return;
42 memset(fit, 0, sizeof(FIT));
43 fit->_voidataNr=0; fit->voiNr=0;
44}

◆ fitIntegralEval()

int fitIntegralEval ( FitVOI * r,
double x,
double * yi )

Evaluates yi=Integral of f(x) between 0 and x.

See also
fitEval, fitDerivEval, fitIntegralEvaltac
Returns
Returns 0, if ok.
Parameters
rFit parameters of a single region
xTime where to evaluate integral of the function
yiThe integral value of the function is returned here

Definition at line 1308 of file mathfunc.c.

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}

Referenced by fitIntegralEvaltac().

◆ fitIntegralEvaltac()

int fitIntegralEvaltac ( FitVOI * r,
double * x,
double * yi,
int dataNr )

Evaluate an array yi[i]=Integral of f(x[i]) between 0 and x.

See also
fitEvaltac, fitIntegralEvaltac, fitDerivEvaltac
Returns
Returns 0, if ok.
Parameters
rFit parameters of a single region
xTimes where to evaluate the function integrals
yiArray for the function integral values
dataNrNr of (x,yi) data

Definition at line 1517 of file mathfunc.c.

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}
fitIntegralEval
int fitIntegralEval(FitVOI *r, double x, double *yi)
Definition mathfunc.c:1308

◆ fitPrint()

void fitPrint ( FIT * fit)

Print to stdout the contents of FIT data structure.

Mainly for testing purposes.

Parameters
fitPointer to FIT struct.

Definition at line 180 of file mathfunc.c.

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}
fitWrite
int fitWrite(FIT *fit, char *filename)
Definition mathfunc.c:54

◆ fitRead()

int fitRead ( char * filename,
FIT * fit,
int verbose )

Read FIT file contents to the specified data structure, emptying its old contents.

Returns
In case of an error, >0 is returned, and a description is written in fiterrmsg.
See also
fitInit, fitPrint, fitWrite
Parameters
filenamePointer to file name.
fitPointer to initiated FIT struct.
verboseVerbose level; if <=0, then nothing is printed into stdout.

Definition at line 196 of file mathfunc.c.

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}
get_datetime
int get_datetime(char *str, struct tm *date, int verbose)
Definition datetime.c:322
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
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
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
fitSetmem
int fitSetmem(FIT *fit, int voiNr)
Definition mathfunc.c:154
fitEmpty
void fitEmpty(FIT *fit)
Definition mathfunc.c:18
fiterrmsg
char fiterrmsg[64]
Definition mathfunc.c:6
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

◆ fitSetmem()

int fitSetmem ( FIT * fit,
int voiNr )

Allocate memory for FIT data. Any previous contents are destroyed.

Returns
Returns 0 when successful, and >0 in case of an error.
See also
fitInit, fitEmpty, fitRead
Parameters
fitPointer to FIT struct.
voiNrNr of TACs to allocate.

Definition at line 154 of file mathfunc.c.

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}

Referenced by fit_allocate_with_dft(), and fitRead().

◆ fitWrite()

int fitWrite ( FIT * fit,
char * filename )

Write function parameters in FIT into specified file.

If necessary, a backup file (+BACKUP_EXTENSION) is created.

Returns
In case of an error, >0 is returned, and a description is written in fiterrmsg.
See also
fitPrint, fitRead, fitInit
Parameters
fitPointer to FIT struct.
filenameFilename.

Definition at line 54 of file mathfunc.c.

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}
backupExistingFile
int backupExistingFile(char *filename, char *backup_ext, char *status)
Definition backup.c:14
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
petTunit
char * petTunit(int tunit)
Definition petunits.c:226

Referenced by fitPrint().

◆ igam()

double igam ( double a,
double x )

Cumulative gamma distribution, or Regularized gamma function, more specifically, lower incomplete gamma function divided by gamma function.

Standard gamma distribution is assumed (Beta=1). f(a,x) = (1/Gamma(a)) * Integral(0,x)(e^-t * t^(a-1))dt

Returns
Returns the value of regularized gamma function, or NaN in case of an error.
See also
igamc, inverfc, factorial
Parameters
aShape parameter alpha; must be > 0.
xIntegral from 0 to x; must be >= 0.

Definition at line 1715 of file mathfunc.c.

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}
igamc
double igamc(double a, double x)
Definition mathfunc.c:1756

Referenced by fitEval(), and igamc().

◆ igamc()

double igamc ( double a,
double x )

Regularized gamma function, more specifically, upper incomplete gamma function divided by gamma function.

f(a,x) = (1/Gamma(a)) * Integral(x,Inf)(e^-t * t^(a-1))dt Standard gamma distribution is assumed (Beta=1).

Returns
Returns the value of regularized gamma function, or NaN in case of an error.
See also
igam, inverfc, factorial
Parameters
aShape parameter alpha; must be > 0.
xIntegral from x to infinity; must be >= 0

Definition at line 1756 of file mathfunc.c.

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}

Referenced by igam().

◆ lfactorial()

unsigned long long int lfactorial ( unsigned long long int n)

Calculate factorial of given number.

See also
factorial, igam, igamc
Returns
Returns factorial n!, or zero if factorial would cause wrap-around.
Parameters
nInteger n, from which the factorial is calculated.
Note
Wrap-around will occur if n>20.

Definition at line 1695 of file mathfunc.c.

1700 {
1701 if(n<1) return(1);
1702 if(n>20) return(0);
1703 return(n*lfactorial(n-1));
1704}

Referenced by fitEval(), fitIntegralEval(), and lfactorial().

Variable Documentation

◆ fiterrmsg

char fiterrmsg[64]

Error message from FIT functions

Definition at line 6 of file mathfunc.c.

Referenced by fitRead(), and fitWrite().

◆ MATHFUNC_TEST

int MATHFUNC_TEST

Verbose prints from FIT functions

Definition at line 5 of file mathfunc.c.

Referenced by fitEval(), fitPrint(), and fitWrite().

Generated on Fri Sep 19 2025 13:32:06 for TPCCLIB by doxygen 1.12.0