tacmodelinput.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcli.h"
#include "tpctac.h"
/*****************************************************************************/
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/
int tacFittime(
  TAC *d,
  double *startTime,
  double *endTime,
  int *first,
  int *last,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(first!=NULL) *first=0; 
  if(last!=NULL) *last=0;
  if(d==NULL || startTime==NULL || endTime==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(-1);
  }
 
  if(d->sampleNr<=0) {
    *startTime=*endTime=0.0; 
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(0);
  }
  if(*endTime<0.0) *endTime=1.0E+100;
 
  /* Change start and end times to seconds if necessary */
  double cf;
  cf=unitConversionFactor(UNIT_MIN, d->tunit);
  if(!isnan(cf)) {*startTime*=cf; *endTime*=cf;}
  if(verbose>1) {
    printf("startTime := %g\n", *startTime);
    printf("endTime := %g\n", *endTime);
  }
 
#if(0)
  /* Check that data range is not outside required range */
  {
    double s, e; 
    if(tacXRange(d, &s, &e) || e<*startTime || s>*endTime) {
      *startTime=*endTime=0.0; 
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X_RANGE);
      return(0);
    }
  }
#endif
 
  /* Get first and last data point inside the range */
  int s, e;
  {
    double f;
    s=e=-1;
    for(int i=0; i<d->sampleNr; i++) {
      if(d->isframe) f=0.5*(d->x1[i]+d->x2[i]); else f=d->x[i];
      if(verbose>3) printf("  f := %g\n", f);
      if(s<0 && f>=*startTime) s=i;
      if(f<=*endTime) e=i; else break;
    }
    if(s<0 || e<0) {
      *startTime=*endTime=0.0; 
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(0);
    }
    if(first!=NULL) *first=s; 
    if(last!=NULL) *last=e;
  }
 
  /* Correct fit range to frame start and end times */
  *startTime=(d->isframe ? d->x1[s] : d->x[s]);
  *endTime=  (d->isframe ? d->x2[e] : d->x[e]);
  /* Change start and end times back to minutes if necessary */
  if(!isnan(cf)) {*startTime/=cf; *endTime/=cf;}
  /* Return the sample number in the fit range */
  return(1+e-s);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacReadModelingData(
  const char *tissuefile,
  const char *inputfile1,
  const char *inputfile2,
  const char *inputfile3,
  double *fitdur,
  int cutInput,
  int *fitSampleNr,
  TAC *tis,
  TAC *inp,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tis==NULL || inp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(tissuefile==NULL || inputfile1==NULL || strnlen(inputfile1, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Check the function input */
  int input_nr=1;
  if(inputfile2!=NULL && strnlen(inputfile2, 1)>0) input_nr++;
  if(inputfile3!=NULL && strnlen(inputfile3, 1)>0) {
    if(input_nr<2) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
      return TPCERROR_NO_DATA;
    }
    input_nr++;
  }
  if(verbose>2) printf("input_nr := %d\n", input_nr);
 
  /* Delete any previous data */
  tacFree(tis); tacFree(inp);
  if(fitSampleNr!=NULL) *fitSampleNr=0;
 
  int ret;
  /*
   *  Read tissue data
   */
  if(verbose>1) printf("reading tissue data in %s\n", tissuefile);
  ret=tacRead(tis, tissuefile, status);
  if(ret!=TPCERROR_OK) return(ret);
 
  /* Do not check frame number; static scan is fine here */
  
  /* Check for NaN's */
  if(tacNaNs(tis)>0) {
    tacFree(tis);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  }
 
  /* Sort the tissue data by increasing sample times */
  ret=tacSortByTime(tis, status);
  if(ret!=TPCERROR_OK) {
    tacFree(tis);
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return ret;
  }
  
  /* Make sure that there is no overlap in sample frame times; samples must
     be sorted before this */
  ret=tacCorrectFrameOverlap(tis, status);
  if(ret!=TPCERROR_OK) {tacFree(tis); return ret;}
 
  /*
   *  Read first input data
   */
  if(verbose>1) printf("reading input data 1 in %s\n", inputfile1);
  ret=tacRead(inp, inputfile1, status);
  if(ret!=TPCERROR_OK) {tacFree(tis); return ret;}
  /* Check and correct the sample time unit */
  if(tis->tunit==UNIT_UNKNOWN) tis->tunit=inp->tunit;
  else if(inp->tunit==UNIT_UNKNOWN) inp->tunit=tis->tunit;
  if(inp->tunit==UNIT_UNKNOWN && verbose>0) {
    fprintf(stderr, "Warning: input sample time units not known.\n");}
  if(tis->tunit!=inp->tunit && tacXUnitConvert(inp, tis->tunit, status)) {
    tacFree(tis); tacFree(inp); return ret;
  }
  /* Check TAC nr */
  if(inp->tacNr>1) {
    if(verbose>0)
      fprintf(stderr, "Warning: using only first TAC in %s\n", inputfile1);
    inp->tacNr=1;
  }
  /* Check sample nr */
  if(inp->sampleNr<4) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    tacFree(tis); tacFree(inp); return TPCERROR_TOO_FEW;
  }
  /* Sort the data by increasing sample times */
  tacSortByTime(inp, NULL);
  /* If inp contains isotope and tis does not, then copy it */
  if(tacGetIsotope(tis)!=ISOTOPE_UNKNOWN) {
    int isotope=tacGetIsotope(inp);
    if(isotope!=ISOTOPE_UNKNOWN) tacSetIsotope(tis, isotope);
  }
 
  /*
   *  Read following input files, if available
   */
  char *fname;
  TAC tmptac; tacInit(&tmptac);
  for(int ii=2; ii<=input_nr; ii++) {
    if(ii==2) fname=(char*)inputfile2; else fname=(char*)inputfile3;
    if(verbose>1) printf("reading input data %d in %s\n", ii, fname);
    ret=tacRead(&tmptac, fname, status);
    if(ret!=TPCERROR_OK) {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;}
    /* Check TAC nr */
    if(tmptac.tacNr>1) {
      if(verbose>0)
        fprintf(stderr, "Warning: using only first TAC in %s\n", fname);
      tmptac.tacNr=1;
    }
    /* Check sample nr */
    if(tmptac.sampleNr<4) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return TPCERROR_TOO_FEW;
    }
    /* Sort the data by increasing sample times */
    tacSortByTime(&tmptac, NULL);
 
    /* Check and correct the sample time unit */
    if(tis->tunit==UNIT_UNKNOWN) tis->tunit=tmptac.tunit;
    else if(tmptac.tunit==UNIT_UNKNOWN) tmptac.tunit=tis->tunit;
    if(inp->tunit!=tmptac.tunit && tacXUnitConvert(&tmptac, inp->tunit, status))
    {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;
    }
 
    /* Check and correct the sample concentration unit */
    if(inp->cunit==UNIT_UNKNOWN) inp->cunit=tmptac.cunit;
    else if(tmptac.cunit==UNIT_UNKNOWN) tmptac.cunit=inp->cunit;
    if(inp->cunit!=tmptac.cunit && tacYUnitConvert(&tmptac, inp->cunit, status))
    {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;
    }
 
    /* Copy to input data */
    if(tacInterpolateInto(&tmptac, inp, NULL, NULL, status)!=TPCERROR_OK) {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;
    }
    
    tacFree(&tmptac);
  } // next input file
 
  /* Set time unit to min */
  ret=tacXUnitConvert(tis, UNIT_MIN, status);
  if(ret && verbose>0) {
    fprintf(stderr, "Warning: check that regional data times are in minutes.\n");
  }
  ret=tacXUnitConvert(inp, UNIT_MIN, status);
  if(ret && verbose>0) {
    fprintf(stderr, "Warning: check that input data times are in minutes.\n");
  }
  /* Check that input and tissue time ranges are about the same */
  double iend, tend;
  {
    ret=tacXRange(inp, NULL, &iend); if(ret==0) ret=tacXRange(tis, NULL, &tend);
    if(ret || iend<=0.0 || tend<=0.0) {
      tacFree(tis); tacFree(inp);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return TPCERROR_INVALID_XRANGE;
    }
    if(tend>10.0*iend || tend<0.10*iend) {
      if(verbose>0) fprintf(stderr, "Warning: check the sample time units.\n");
    }
  }
  ret=tacYUnitConvert(inp, tis->cunit, status);
  if(ret && verbose>0) {
    fprintf(stderr, "Warning: check the calibration units.\n");
  }
 
  /*
   *  Check and set fit time length
   */
  if(verbose>1) printf("checking and setting fit time length\n");
  /* Set fit duration */
  double starttime=0, endtime=*fitdur;
  int fnr=tacFittime(tis, &starttime, &endtime, NULL, NULL, status);
  if(verbose>3) {
    fprintf(stdout, "tis.sampleNr := %d\n", tis->sampleNr);
    fprintf(stdout, "starttime := %g\n", starttime);
    fprintf(stdout, "endtime := %g\n", endtime);
    //fprintf(stdout, "first := %d\n", first);
    //fprintf(stdout, "last := %d\n", last);
    fprintf(stdout, "fitSampleNr := %d\n", fnr);
  }
  *fitdur=endtime;
  if(fitSampleNr!=NULL) *fitSampleNr=fnr;
 
  /* Check that input data does not end much before fitdur */
  if(*fitdur>1.2*iend) {
    tacFree(tis); tacFree(inp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    return TPCERROR_TOO_FEW;
  }
 
  /* Cut off too many input samples to make calculation faster */
  if(cutInput && iend>*fitdur) {
    if(verbose>1) printf("cut off too many input samples\n");
    int i; double f;
    for(i=0; i<inp->sampleNr; i++) {
      if(inp->isframe) f=0.5*(inp->x1[i]+inp->x2[i]); else f=inp->x[i];
      if(f>(*fitdur)) break;
    }
    if(i<inp->sampleNr) i++; 
    inp->sampleNr=i;
    if(inp->sampleNr<4) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
      tacFree(tis); tacFree(inp); return TPCERROR_TOO_FEW;
    }
  }
  if(verbose>2) fprintf(stdout, "inp.sampleNr := %d\n", inp->sampleNr);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacReadReference(
  TAC *tis,
  const char *reference,
  TAC *ref,
  int *refIndex,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tis==NULL || reference==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return 0;
  }
  if(strnlen(reference, 1)<1 || tis->sampleNr<1 || !tacIsX(tis)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return 0;
  }
 
  /* First, try to read reference from file */
  TAC temp; tacInit(&temp);
  if(verbose>1) printf("trying to read file %s\n", reference);
  if(tacRead(&temp, reference, status)==TPCERROR_OK) {
    /* Check sample nr */
    if(temp.sampleNr<tis->sampleNr) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
      tacFree(&temp); return 0;
    }
    int n, ret=0, refi=0;
    /* Convert units */
    if(tacXUnitConvert(&temp, tis->tunit, status)) ret++;
    if(tacYUnitConvert(&temp, tis->cunit, status)) ret++;
    if(verbose>0 && ret>0)
      fprintf(stderr, "Warning: check the units of reference and tissue data.\n");
    /* Check that sample times do match */
    if(!tacXMatch(tis, &temp, verbose-1)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      tacFree(&temp); return 0;
    } else if(tis->isframe && !temp.isframe) {
      /* Copy frame start and end times, in case those had been lost */
      tacXCopy(tis, &temp, 0, tis->sampleNr-1);
      temp.isframe=1;
    }
    /* Select the best reference region */
    if(temp.tacNr==1) {
      temp.c[0].sw=2; refi=0; // just one
    } else {
      for(int i=0; i<temp.tacNr; i++) temp.c[i].sw=1;
      refi=tacSelectBestReference(&temp); if(refi<0) refi=0;
      temp.c[refi].sw=2;
    }
    if(verbose>1) printf("selected ref region := %s\n", temp.c[refi].name);
    /* If ref contains isotope and tis does not, then copy it */
    if(tacGetIsotope(tis)!=ISOTOPE_UNKNOWN) {
      int isotope=tacGetIsotope(&temp);
      if(isotope!=ISOTOPE_UNKNOWN) tacSetIsotope(tis, isotope);
    }
    /* Copy to output struct, if pointer is given; otherwise copy to tis */
    if(ref!=NULL) {
      if(verbose>2) printf("adding %d ref tacs to ref struct\n", temp.tacNr);
      ret=tacDuplicate(&temp, ref);
      if(ret!=TPCERROR_OK) {
        statusSet(status, __func__, __FILE__, __LINE__, ret);
        tacFree(&temp); return 0;
      }
      if(refIndex!=NULL) *refIndex=refi;
    } else {
      if(verbose>2) printf("adding %d ref tacs to tis struct\n", temp.tacNr);
      for(int i=0; i<tis->tacNr; i++) tis->c[i].sw=0;
      ret=tacAllocateMore(tis, temp.tacNr);
      for(int i=0; i<temp.tacNr && ret==TPCERROR_OK; i++)
        ret=tacCopyTacc(&temp.c[i], &tis->c[tis->tacNr+i], tis->sampleNr);
      if(ret!=TPCERROR_OK) {
        statusSet(status, __func__, __FILE__, __LINE__, ret);
        tacFree(&temp); return 0;
      }
      if(refIndex!=NULL) *refIndex=tis->tacNr+refi;
      tis->tacNr+=temp.tacNr;
    }
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    n=temp.tacNr; tacFree(&temp);
    if(verbose>1) printf("%d ref region(s) found in file.\n", n);
    return(n);
    // reference file was read and processed
  } else if(verbose>1) {
    printf("  '%s' was not a file\n", reference); fflush(stdout);
  }
 
  /* So 'reference' did not contain valid, accessible filename,
     but is it a region name or number? */
  if(verbose>1) printf("trying to find reference in TAC struct\n");
  int n, refi;
  n=tacSelectTACs(tis, reference, 1, status);
  if(verbose>1) {printf("%d region(s) matched.\n", n); fflush(stdout);}
  if(n<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_REFERENCE);
    return(0);
  }
  if(n==tis->tacNr) {
    if(verbose>2) {printf("... which means all region(s) matched.\n"); fflush(stdout);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(0);
  }
  refi=tacSelectBestReference(tis);
  if(refi<0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(0);
  }
  tis->c[refi].sw=2; if(refIndex!=NULL) *refIndex=refi;
 
  /* If output struct was not given, then we are ready */
  if(ref==NULL) {
    if(verbose>1) {printf("reference TAC tagged.\n"); fflush(stdout);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(n);
  }
 
  /* Move reference TACs to output struct */
  if(verbose>1) printf("moving reference TAC(s) to output TAC struct\n");
  int ret;
  ret=tacDuplicate(tis, ref);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return 0;
  }
  int i=ref->tacNr-1; ret=0;
  while(i>=0 && !ret) {
    if(ref->c[i].sw==0) ret=tacDeleteTACC(ref, i); 
    i--;
  }
  i=tis->tacNr-1;
  while(i>=0 && !ret) {
    if(tis->c[i].sw) ret=tacDeleteTACC(tis, i); 
    i--;
  }
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return 0;
  }
  if(refIndex!=NULL) { // set reference index to index inside ref struct
    *refIndex=-1;
    for(i=0; i<ref->tacNr; i++) if(ref->c[i].sw==2) {*refIndex=i; break;}
  }
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(n);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacReadModelingInput(
  const char *inputfile1,
  const char *inputfile2,
  const char *inputfile3,
  TAC *inp,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(inp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(inputfile1==NULL || strnlen(inputfile1, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }  
 
  /* Check the function input */
  int input_nr=1;
  if(inputfile2!=NULL && strnlen(inputfile2, 1)>0) input_nr++;
  if(inputfile3!=NULL && strnlen(inputfile3, 1)>0) {
    if(input_nr<2) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
      return TPCERROR_NO_DATA;
    }
    input_nr++;
  }
  if(verbose>2) {printf("input_nr := %d\n", input_nr); fflush(stdout);}
 
  /* Delete any previous data */
  tacFree(inp);
 
  int ret;
  /*
   *  Read first input data
   */
  if(verbose>1) {
    printf("reading input data 1 in %s\n", inputfile1); fflush(stdout);
  }
  ret=tacRead(inp, inputfile1, status);
  if(ret!=TPCERROR_OK) {return ret;}
  /* Check for NaN's */
  if(tacNaNs(inp)>0) {
    tacFree(inp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  }
  /* Sort the data by increasing sample times */
  ret=tacSortByTime(inp, status);
  if(ret!=TPCERROR_OK) {
    tacFree(inp); return ret;
  }
  /* Check TAC nr */
  if(inp->tacNr>1) {
    if(verbose>0)
      fprintf(stderr, "Warning: using only first TAC in %s\n", inputfile1);
    inp->tacNr=1;
  }
  /* Check sample nr */
  if(inp->sampleNr<4) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    tacFree(inp); return TPCERROR_TOO_FEW;
  }
 
  /*
   *  Read following input files, if available
   */
  char *fname;
  TAC tmptac; tacInit(&tmptac);
  for(int ii=2; ii<=input_nr; ii++) {
    if(ii==2) fname=(char*)inputfile2; else fname=(char*)inputfile3;
    if(verbose>1) {
      printf("reading input data %d in %s\n", ii, fname);
      fflush(stdout);
    }
    ret=tacRead(&tmptac, fname, status);
    if(ret!=TPCERROR_OK) {
      tacFree(inp); tacFree(&tmptac); return ret;}
    /* Check TAC nr */
    if(tmptac.tacNr>1) {
      if(verbose>0)
        fprintf(stderr, "Warning: using only first TAC in %s\n", fname);
      tmptac.tacNr=1;
    }
    /* Check sample nr */
    if(tmptac.sampleNr<4) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
      tacFree(inp); tacFree(&tmptac); return TPCERROR_TOO_FEW;
    }
    /* Sort the data by increasing sample times */
    tacSortByTime(&tmptac, status);
    if(ret!=TPCERROR_OK) {
      tacFree(inp); tacFree(&tmptac); return ret;
    }
    /* Check and correct the sample time unit */
    if(inp->tunit==UNIT_UNKNOWN) inp->tunit=tmptac.tunit;
    else if(tmptac.tunit==UNIT_UNKNOWN) tmptac.tunit=inp->tunit;
    if(inp->tunit!=tmptac.tunit && tacXUnitConvert(&tmptac, inp->tunit, status))
    {
      tacFree(inp); tacFree(&tmptac); return ret;
    }
    /* Check and correct the sample concentration unit */
    if(inp->cunit==UNIT_UNKNOWN) inp->cunit=tmptac.cunit;
    else if(tmptac.cunit==UNIT_UNKNOWN) tmptac.cunit=inp->cunit;
    if(inp->cunit!=tmptac.cunit && tacYUnitConvert(&tmptac, inp->cunit, status))
    {
      tacFree(inp); tacFree(&tmptac); return ret;
    }
    /* Copy to input data */
    if(tacInterpolateInto(&tmptac, inp, NULL, NULL, status)!=TPCERROR_OK) {
      tacFree(inp); tacFree(&tmptac); return ret;
    }
    tacFree(&tmptac);
  } // next input file
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/