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 fitEvaltacframes (FitVOI *r, double *x1, double *x2, double *y, 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 1707 of file mathfunc.c.

1712 {
1713 if(n<1) return(1);
1714 if(n>12) return(0);
1715 return(n*factorial(n-1));
1716}
factorial
unsigned int factorial(unsigned int n)
Definition mathfunc.c:1707

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 1569 of file mathfunc.c.

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

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 1684 of file mathfunc.c.

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

◆ 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: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

Referenced by fitRead(), and fitSetmem().

◆ fitEval()

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

Evaluate y=f(x).

See also
fitEvaltac, fitEvaltacframes, 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 618 of file mathfunc.c.

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

Referenced by fitEvaltac(), and fitEvaltacframes().

◆ fitEvaltac()

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

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

See also
fitEvaltacframes, 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 1234 of file mathfunc.c.

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

◆ fitEvaltacframes()

int fitEvaltacframes ( FitVOI * r,
double * x1,
double * x2,
double * y,
int dataNr )

Evaluate an array y[i] of average f(x) between frame start (x1[i]) and end times (x2[i]). Many functions cannot currently be evaluated using this.

See also
fitEvaltac, fitEval, fitIntegralEvaltac
Returns
Returns 0, if ok.
Parameters
rFit parameters of a single region.
x1Frame start times.
x2Frame end times.
yArray for the function value means between x1 and x2.
dataNrNr of (x,yi) data.

Definition at line 1536 of file mathfunc.c.

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

◆ 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 351 of file mathfunc.c.

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}

◆ 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 511 of file mathfunc.c.

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}

◆ 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 1290 of file mathfunc.c.

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}

Referenced by fitEvaltacframes(), and 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, fitIntegralEval, 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 1512 of file mathfunc.c.

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}

◆ 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 /* 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}
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
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
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
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

◆ 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 1744 of file mathfunc.c.

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

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 1784 of file mathfunc.c.

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}

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 1724 of file mathfunc.c.

1729 {
1730 if(n<1) return(1);
1731 if(n>20) return(0);
1732 return(n*lfactorial(n-1));
1733}

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 for TPCCLIB by doxygen 1.15.0