delay.c File Reference

Time delay between PET tissue and blood curve. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "tpcextensions.h"
#include "tpctac.h"
#include "tpcli.h"
#include "tpclinopt.h"
#include "tpctacmod.h"

Go to the source code of this file.

Functions
int	tacDelay (TAC *tac, double dt, int ti, TPCSTATUS *status)
	Move TAC y values (concentrations) in time, keeping sample times (x values) intact.
void	initDelayCMFitData (DELAYCMFITDATA *d)
	Before first use, initiate the data structure for delay time estimation.
void	freeDelayCMFitData (DELAYCMFITDATA *d)
	After last use, free memory in the data structure for delay time estimation.
int	tacDelayCMFit (TAC *btac, TAC *ttac, int ci, double dtmin, double dtmax, double fitend, double dtstep, double *dt, int mode, int model, DELAYCMFITDATA *tdata, TPCSTATUS *status)
	Fit time delay between PET tissue and plasma or blood curve.

Detailed Description

Time delay between PET tissue and blood curve.

Definition in file delay.c.

Function Documentation

freeDelayCMFitData()

void freeDelayCMFitData

(

DELAYCMFITDATA *

d

)

After last use, free memory in the data structure for delay time estimation.

See also

tacDelayCMFit, initDelayCMFitData

Parameters

d	Pointer to the data structure to be initialized.

Definition at line 150 of file delay.c.

  {
  if(d==NULL) return;
  d->iNr=d->tNr=d->moveNr=0;
  tacFree(&d->dtac); tacFree(&d->ditac); tacFree(&d->diitac); tacFree(&d->ddtac); tacFree(&d->dddtac);
  d->mode=0; d->model=0;
  d->llsq_n=d->llsq_m=0;
  free(d->llsq_a); free(d->llsq_mat);
  free(d->llsq_b); free(d->llsq_x); free(d->llsq_wp); free(d->llsq_zz); 
  free(d->llsq_i);
  return;
}

Referenced by tacDelayCMFit().

initDelayCMFitData()

void initDelayCMFitData

(

DELAYCMFITDATA *

d

)

Before first use, initiate the data structure for delay time estimation.

See also

tacDelayCMFit, freeDelayCMFitData

Parameters

d	Pointer to the data structure to be initialized.

Definition at line 131 of file delay.c.

  {
  if(d==NULL) return;
  d->iNr=d->tNr=d->moveNr=0;
  tacInit(&d->dtac); tacInit(&d->ditac); tacInit(&d->diitac); tacInit(&d->ddtac); tacInit(&d->dddtac);
  d->mode=0; d->model=0;
  d->llsq_n=d->llsq_m=0;
  d->llsq_mat=NULL;
  d->llsq_a=NULL;
  d->llsq_b=d->llsq_x=d->llsq_wp=d->llsq_zz=NULL; 
  d->llsq_i=NULL;
  return;
}

Referenced by tacDelayCMFit().

tacDelay()

int tacDelay	(	TAC *	tac,
		double	dt,
		int	ti,
		TPCSTATUS *	status )

Move TAC y values (concentrations) in time, keeping sample times (x values) intact.

If TAC has frame start and end times, this function tries to retain the AUC with stepwise frame interpretation. If you wish to prevent this, then set tac->isframe=0 before calling this function and after it set back tac->isframe=1.

See also

tacDelayCMFit, tacInterpolate, tacFramesToSteps, simDispersion, liInterpolate, tacVb

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

tac	Pointer to TAC structure.

Precondition

Data must be sorted by increasing x. No gaps or overlap between frames allowed.

See also

tacSortByTime, tacSetXContiguous

Postcondition

Data y values are modified, x values will stay the same.

Parameters

dt	Time (in TAC tunit's) to move the TAC forward (negative time moves TAC backward).
ti	Index [0..tacNr-1] of the TAC to move, or <0 to move all TACs.
status	Pointer to status data; enter NULL if not needed.

Definition at line 29 of file delay.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(tac, %g, %d)\n", __func__, dt, ti);
  /* Check the function input */
  if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(ti>=tac->tacNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(fabs(dt)<1.0E-12) {
    if(verbose>1) printf("  requested delay time is zero; nothing done.\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK);
  }
  /* Check that input does not have any missing values */
  if(tacXNaNs(tac)>0 || tacYNaNs(tac, ti)>0) {
    if(verbose>1) printf("  data contains NaN(s).\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  } 
 
  /* Get and check the x range */
  double xmin=nan(""), xmax=nan("");
  if(tacXRange(tac, &xmin, &xmax)!=0 || dt>=xmax || dt<=-xmax) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  if(verbose>2) printf("  xrange: %g - %g\n", xmin, xmax);
 
 
  /* If TAC contains just frame mid times, then simple dot-to-dot interpolation */
  if(tac->isframe==0) {
    /* Make temporary x array */
    double tx[tac->sampleNr];
    for(int i=0; i<tac->sampleNr; i++) tx[i]=tac->x[i]+dt;
    /* One TAC at a time */
    for(int j=0; j<tac->tacNr; j++) if(ti<0 || ti==j) {
      /* Make temporary y array */
      double ty[tac->sampleNr];
      for(int i=0; i<tac->sampleNr; i++) ty[i]=tac->c[j].y[i];
      /* Interpolate temporary data to original times */
      if(liInterpolate(tx, ty, tac->sampleNr, tac->x, tac->c[j].y, NULL, NULL, tac->sampleNr, 0, 1, 
                       verbose-10)!=0)
      {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
        return TPCERROR_INVALID_VALUE;
      }
    }
    /* That's it then */
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK);
  }
 
  /* We have frame start and end times */
  /* Convert the data into stepwise dot-to-dot data */
  TAC stac; tacInit(&stac);
  {
    int ret=tacFramesToSteps(tac, &stac, NULL);
    if(ret!=TPCERROR_OK) {
      tacFree(&stac);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }
  }
  /* Change times of the stepwise data */
  for(int i=0; i<stac.sampleNr; i++) stac.x[i]+=dt;
  /* Interpolate stepwise data into original frame times, one TAC at a time */
  for(int j=0; j<tac->tacNr; j++) if(ti<0 || ti==j) {
    if(liInterpolateForPET(stac.x, stac.c[j].y, stac.sampleNr, tac->x1, tac->x2, tac->c[j].y, 
                           NULL, NULL, tac->sampleNr, 0, 1, verbose-10))
    {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return TPCERROR_INVALID_VALUE;
    }
  }
  tacFree(&stac);
 
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacDelayCMFit()

int tacDelayCMFit	(	TAC *	btac,
		TAC *	ttac,
		int	ci,
		double	dtmin,
		double	dtmax,
		double	fitend,
		double	dtstep,
		double *	dt,
		int	mode,
		int	model,
		DELAYCMFITDATA *	tdata,
		TPCSTATUS *	status )

Fit time delay between PET tissue and plasma or blood curve.

Initial part of tissue curve is fitted with a range of input curves moved in time, using linearised one or two-tissue compartment model with blood volume term. The time that provides the best fit is selected. It is assumed that tissue heterogeneity, dispersion of the blood curve, metabolites, and variable plasma-to-blood ratio will only minimally impact the time delay estimate, since only the initial part of the data is fitted, especially if reversible two-tissue model (R2tCM) is used.

See also

tacDelay, simDispersion

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

btac	Pointer to input data (blood or plasma curve). Units of input and tissue TAC must match. All but first curve in TAC are ignored. Enter NULL in subsequent runs
ttac	Pointer to tissue data. Units of input and tissue TAC must match. All but first curve in TAC are ignored.
ci	Index [0..ttac->tacNr-1] of tissue tac to use.
dtmin	Minimum delay to be tested, for example -10 s, in the same units as in TACs. Not used in subsequent runs with tdata.
dtmax	Maximum delay to be tested, for example +50 s, in the same units as in TACs. Not used in subsequent runs with tdata.
fitend	Fitting uses data from 0 to end time, for example 120 s, in the same units as in TACs. In case of short PET scan, end time may need to be reduced so that input TAC with negative delay extends to the fit end. Not used in subsequent runs with tdata.
dtstep	The step length for moving the input curve, for example 1 s, in the same units as in TACs. Not used in subsequent runs with tdata.
dt	Pointer for the time delay estimate. Positive delay time means that tissue curve is delayed as compared to the input curve, and vice versa. Thus, input curve needs to be moved by the delay time to match the tissue curve.
mode	Mode for multilinear equations: 0=integrals and 2nd integrals, 1=derivates and 2nd derivates (experimental), 2=integrals and derivates (experimental).
model	Model: 0=R1TCM, 1=I2TCM, 2=R2TCM.
tdata	Pointer to temporary data structure, which quickens subsequent runs with other tissue TACs; enter NULL if not needed.

See also

initDelayCMFitData, freeDelayCMFitData

Parameters

status

Pointer to status data; enter NULL if not needed.

Definition at line 179 of file delay.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>2) {
    printf("%s(inp, tis, %d, %g, %g, %g, %g, %d, ...)\n", __func__, ci, dtmin, dtmax, fitend, dtstep, mode);
    fflush(stdout);
  }
 
 
  /* Check provided data */
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
  if(ttac==NULL || (btac==NULL && tdata==NULL)) return(TPCERROR_FAIL);
  if(btac!=NULL && (btac->tacNr<1 || btac->sampleNr<4)) return(TPCERROR_FAIL);
  if(ci<0 || ci>=ttac->tacNr || ttac->sampleNr<4) return(TPCERROR_FAIL);
  if(model<0 || model>2) return(TPCERROR_FAIL);
  if(dt!=NULL) *dt=nan("");
 
  /* Initiate (possibly reused) data structure */
  DELAYCMFITDATA *b, ltdata; initDelayCMFitData(&ltdata);
  if(tdata!=NULL) {
    b=tdata;
    if(verbose>3) {
      if(b->moveNr==0) printf("  given temp data will be filled\n"); 
      else printf("  filled temp data will be used\n");
      fflush(stdout);
    }
  } else {
    b=&ltdata;
    if(verbose>4) {printf("  local temp data created\n"); fflush(stdout);}
  }
 
  /* Check and set temp data when run first time with the given TTAC */
  if(btac!=NULL || b->moveNr==0) {
    /* Delay range and steps*/
    if(!(dtmax>dtmin) || !(dtstep>0.001)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    b->dtmin=dtmin; b->dtmax=dtmax; b->dtstep=dtstep;
    /* How many delayed curves needed */
    b->moveNr=1+(b->dtmax-b->dtmin)/b->dtstep;
    if(verbose>3) printf("  moveNr := %d\n", b->moveNr);
    if(b->moveNr<5 || b->moveNr>5000) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    /* Get x range */
    if(tacXRange(ttac, &b->ttac_xmin, &b->ttac_xmax)!=0 || tacXRange(btac, &b->btac_xmin, &b->btac_xmax)!=0) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    if(verbose>3) printf("  TTAC range := %g - %g\n  BTAC range := %g - %g\n",
                         b->ttac_xmin, b->btac_xmax, b->ttac_xmin, b->btac_xmax);
    /* Check that time units do match */
    if(btac->tunit!=ttac->tunit) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return(TPCERROR_INVALID_X);
    }
    /* Check that input curve is long enough to fit duration */
    if(!(fitend>0.0)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    b->fitend=fitend;
    if(b->fitend>b->ttac_xmax) b->fitend=b->ttac_xmax;
    if(b->fitend>b->btac_xmax) b->fitend=b->btac_xmax;
    if(b->dtmin<0.0 && b->btac_xmax<(b->fitend+b->dtmin))
      b->fitend+=b->dtmin;
    if(verbose>3) printf("fitend=%g -> %g\n", fitend, b->fitend);
    b->iNr=0; for(unsigned int i=0; i<(unsigned int)btac->sampleNr; i++) if(btac->x[i]<=b->fitend) b->iNr++; else break;
    b->tNr=0; for(unsigned int i=0; i<(unsigned int)ttac->sampleNr; i++) if(ttac->x[i]<=b->fitend) b->tNr++; else break;
    //printf("iNr=%u tNr=%u\n", b->iNr, b->tNr);
    if(b->iNr<5 || b->tNr<5) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
 
    /* Further, when run the first time, make input TACs moved in time, interpolated and integrated 
       to tissue times */
    if(verbose>5) {printf("  preparing delayed input data\n"); fflush(stdout);}
    b->mode=mode;
    b->model=model;
    int ret;
    ret=tacAllocate(&b->dtac, b->tNr, b->moveNr+1); // additional place for tissue curve
    if(ret!=TPCERROR_OK) {
      freeDelayCMFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(TPCERROR_OUT_OF_MEMORY);
    }
    b->dtac.sampleNr=b->tNr;
    b->dtac.tacNr=b->moveNr+1;  // tissue as the last curve
    b->dtac.isframe=ttac->isframe;
    for(int i=0; i<b->dtac.sampleNr; i++) {
      b->dtac.x1[i]=ttac->x1[i]; b->dtac.x[i]=ttac->x[i]; b->dtac.x2[i]=ttac->x2[i];}
    ret=tacDuplicate(&b->dtac, &b->ditac);
    if(ret==TPCERROR_OK) ret=tacDuplicate(&b->dtac, &b->diitac); 
    if(ret==TPCERROR_OK) ret=tacDuplicate(&b->dtac, &b->ddtac); 
    if(ret==TPCERROR_OK) ret=tacDuplicate(&b->dtac, &b->dddtac); 
    if(ret!=TPCERROR_OK) {
      freeDelayCMFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }
    /* Integrate/derivate input data, change times, and interpolate to tissue data */
    double bd[btac->sampleNr], bdd[btac->sampleNr];
    double bi[btac->sampleNr], bii[btac->sampleNr];
    if(liDerivate3(btac->x, btac->c[0].y, btac->sampleNr, bd, bdd, 0)) {
      freeDelayCMFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return(TPCERROR_INVALID_VALUE);
    }
    if(btac->isframe==0) {
      if(liIntegrate(btac->x, btac->c[0].y, btac->sampleNr, bi, 3, 0) ||
         liIntegrate(btac->x, bi, btac->sampleNr, bii, 3, 0))
        ret=TPCERROR_INVALID_VALUE;
    } else {
      if(liIntegrateFE(btac->x1, btac->x2, btac->c[0].y, btac->sampleNr, bi, bii, 0))
        ret=TPCERROR_INVALID_VALUE;
    }
    if(ret) {
      freeDelayCMFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }
    double dx[btac->sampleNr], dx2[btac->sampleNr], *dxp;
    if(btac->isframe==0) dxp=dx; else dxp=dx2;
    for(int di=0; di<b->dtac.tacNr-1; di++) {
      for(int i=0; i<btac->sampleNr; i++) {
        dx[i]=btac->x[i]+b->dtmin+di*b->dtstep;
        dx2[i]=btac->x2[i]+b->dtmin+di*b->dtstep;
      }
      if(b->dtac.isframe==0) {
        ret=liInterpolate(dx, btac->c[0].y, btac->sampleNr, b->dtac.x, b->dtac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
    
        if(!ret) ret=liInterpolate(dxp, bi, btac->sampleNr, b->dtac.x, b->ditac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolate(dxp, bii, btac->sampleNr, b->dtac.x, b->diitac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolate(dx, bd, btac->sampleNr, b->dtac.x, b->ddtac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolate(dx, bdd, btac->sampleNr, b->dtac.x, b->dddtac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
      } else {
        ret=liInterpolateForPET(dx, btac->c[0].y, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->dtac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dxp, bi, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->ditac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dxp, bii, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->diitac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dx, bd, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->ddtac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dx, bdd, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->dddtac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
      }
      if(ret) {
        freeDelayCMFitData(&ltdata);
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
        return(TPCERROR_INVALID_XRANGE);
      }
    }
    if(verbose>6) {
      printf("integrated input curve 1:\n");
      for(int i=0; i<b->dtac.sampleNr; i++) 
        printf("%g  %g %g %g\n", b->dtac.x[i], b->dtac.c[0].y[i], b->ditac.c[0].y[i], b->diitac.c[0].y[i]);
      fflush(stdout);
    }
 
 
    /* Allocate memory for linear LSQ */
    if(verbose>5) {printf("  preparing data for LLSQ\n"); fflush(stdout);}
    b->llsq_n=3+b->model; 
    b->llsq_m=b->dtac.sampleNr;
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    b->llsq_mat=(double*)malloc((b->llsq_n*b->llsq_m)*sizeof(double));
    if(b->llsq_mat==NULL) {freeDelayCMFitData(&ltdata); return(TPCERROR_OUT_OF_MEMORY);}
    b->llsq_a=(double**)malloc(b->llsq_n*sizeof(double*));
    if(b->llsq_a==NULL) {freeDelayCMFitData(&ltdata); return(TPCERROR_OUT_OF_MEMORY);}
    for(int ni=0; ni<b->llsq_n; ni++) b->llsq_a[ni]=b->llsq_mat+ni*b->llsq_m;
    b->llsq_b=(double*)malloc(b->llsq_m*sizeof(double));
    b->llsq_x=(double*)malloc(b->llsq_n*sizeof(double));
    b->llsq_wp=(double*)malloc(b->llsq_n*sizeof(double));
    b->llsq_zz=(double*)malloc(b->llsq_m*sizeof(double));
    b->llsq_i=(int*)malloc(b->llsq_n*sizeof(int));
    if(b->llsq_b==NULL || b->llsq_x==NULL || b->llsq_wp==NULL || b->llsq_zz==NULL || b->llsq_i==NULL) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
      freeDelayCMFitData(&ltdata); return(TPCERROR_OUT_OF_MEMORY);}
  }
 
  /* Integrate tissue curve */
  {
    if(verbose>5) {printf("  preparing tissue data\n"); fflush(stdout);}
    int ti=b->dtac.tacNr-1; // index of tissue curve after the input curves
    for(int i=0; i<b->dtac.sampleNr; i++) b->dtac.c[ti].y[i]=ttac->c[ci].y[i];
 
    int ret=liDerivate(b->dtac.x, b->dtac.c[ti].y, b->dtac.sampleNr, b->ddtac.c[ti].y, b->dddtac.c[ti].y, 0);
    if(b->dtac.isframe==0) {
      if(!ret) ret=liIntegrate(b->dtac.x, b->dtac.c[ti].y, b->dtac.sampleNr, b->ditac.c[ti].y, 3, 0);
      if(!ret) ret=liIntegrate(b->dtac.x, b->ditac.c[ti].y, b->dtac.sampleNr, b->diitac.c[ti].y, 3, 0);
    } else {
      if(!ret) ret=liIntegrateFE(b->dtac.x1, b->dtac.x2, b->dtac.c[ti].y, b->dtac.sampleNr, 
                                 b->ditac.c[ti].y, b->diitac.c[ti].y, 0);
    }
    if(ret) {
      freeDelayCMFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
 
    if(verbose>10) {
      printf("integrated tissue curves:\n");
      for(int i=0; i<b->dtac.sampleNr; i++) 
        printf("%g  %g %g %g\n", b->dtac.x[i], b->dtac.c[ti].y[i], b->ditac.c[ti].y[i], b->diitac.c[ti].y[i]);
      fflush(stdout);
    }
  }
 
 
  /* Multilinear fitting with each delayed input */
  if(verbose>5) {printf("  LLSQ\n"); fflush(stdout);}
  double ss, ssmin=nan("");
  int mini=-1;
  for(int di=0; di<b->dtac.tacNr-1; di++) {
    /* Setup data matrix A and vector B */
    for(int mi=0; mi<b->llsq_m; mi++) {
      b->llsq_b[mi]=b->dtac.c[b->dtac.tacNr-1].y[mi]; // TTAC
      b->llsq_mat[mi]=b->dtac.c[di].y[mi]; // BTAC
      if(mode==1) {
        b->llsq_mat[mi+1*b->llsq_m]=b->ddtac.c[di].y[mi]; // BTAC derivate
        b->llsq_mat[mi+2*b->llsq_m]=-b->ddtac.c[b->dtac.tacNr-1].y[mi]; // TTAC derivate
        if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->dddtac.c[di].y[mi]; // BTAC derivate
        if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->dddtac.c[b->dtac.tacNr-1].y[mi]; // TTAC derivate
      } else if(mode==2) {
        b->llsq_mat[mi+1*b->llsq_m]=b->ditac.c[di].y[mi]; // BTAC integral
        b->llsq_mat[mi+2*b->llsq_m]=-b->ditac.c[b->dtac.tacNr-1].y[mi]; // TTAC integral
        if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->ddtac.c[di].y[mi]; // BTAC derivate
        if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->ddtac.c[b->dtac.tacNr-1].y[mi]; // TTAC derivate
      } else { // mode 0
        b->llsq_mat[mi+1*b->llsq_m]=b->ditac.c[di].y[mi]; // BTAC integral
        b->llsq_mat[mi+2*b->llsq_m]=-b->ditac.c[b->dtac.tacNr-1].y[mi]; // TTAC integral
        if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->diitac.c[di].y[mi]; // BTAC 2nd integral
        if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->diitac.c[b->dtac.tacNr-1].y[mi]; // TTAC 2nd integral
      }
      b->llsq_mat[mi]=b->dtac.c[di].y[mi]; // BTAC
      b->llsq_mat[mi+1*b->llsq_m]=b->ditac.c[di].y[mi]; // BTAC integral
      b->llsq_mat[mi+2*b->llsq_m]=-b->ditac.c[b->dtac.tacNr-1].y[mi]; // TTAC integral
      if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->diitac.c[di].y[mi]; // BTAC 2nd integral
      if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->diitac.c[b->dtac.tacNr-1].y[mi]; // TTAC 2nd integral
    }
    if(verbose>3) {
      int failnr=0;
      for(int mi=0; mi<b->llsq_m; mi++) {
        if(!isfinite(b->llsq_b[mi])) failnr++;
        if(!isfinite(b->llsq_mat[mi])) failnr++;
        if(!isfinite(b->llsq_mat[mi+1*b->llsq_m])) failnr++;
        if(!isfinite(b->llsq_mat[mi+2*b->llsq_m])) failnr++;
        if(!isfinite(b->llsq_mat[mi+3*b->llsq_m])) failnr++;
        if(!isfinite(b->llsq_mat[mi+4*b->llsq_m])) failnr++;
      }
      if(failnr>0) printf("%d bad values\n", failnr);
    }
    if(verbose>9 && b->llsq_m<100) {
      printf("---- B x A ----\n");
      for(int mi=0; mi<b->llsq_m; mi++) {
        printf("  %g ", b->llsq_b[mi]);
        for(int ni=0; ni<b->llsq_n; ni++)
          printf(" %g", b->llsq_mat[mi+ni*b->llsq_m]);
         printf("\n");
      }
    }
 
    int ret=nnls(b->llsq_a, b->llsq_m, b->llsq_n, b->llsq_b, b->llsq_x, &ss, b->llsq_wp, b->llsq_zz, b->llsq_i);
    if(ret<2 && isfinite(ss)) {
      if(!isfinite(ssmin) || ssmin>ss) {ssmin=ss; mini=di;}
      if(verbose>2) printf("   %d %d %g %g\n", di, mini, ss, ssmin);
    } else {
      if(verbose>1) printf("  failed NNLS\n");
    }
  }
  if(mini<0) {
    if(verbose>1) fprintf(stderr, "Error: all fits failed.\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
    freeDelayCMFitData(&ltdata); return(TPCERROR_BAD_FIT);
  }
  if(dt!=NULL && mini>=0) *dt=b->dtmin+mini*b->dtstep;
 
  freeDelayCMFitData(&ltdata); // the local one, not the one possibly given by caller 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}