tacw.c File Reference

Working with weights in TAC structure. More...

#include "tpcclibConfig.h"
#include "tpcift.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "tpctac.h"

Go to the source code of this file.

Functions
int	tacIsWeighted (TAC *tac)
int	tacWCopy (TAC *tac1, TAC *tac2, int i1, int i2)
int	tacWMove (TAC *tac, int ow, TPCSTATUS *status)
int	tacWByFreq (TAC *tac, isotope isot, TPCSTATUS *status)
unsigned int	tacWSampleNr (TAC *tac)
int	tacWeightNorm (TAC *tac, TPCSTATUS *status)
int	tacWeightModerate (TAC *tac, const double minprop, const int doZeroes, const int doNaNs, TPCSTATUS *status)
int	sifWeight (TAC *sif, isotope isot, TPCSTATUS *status)
int	tacSetWeights (TAC *tac, weights weightMethod, int weightNr, TPCSTATUS *status)

Detailed Description

Working with weights in TAC structure.

Definition in file tacw.c.

Function Documentation

sifWeight()

int sifWeight	(	TAC *	sif,
		isotope	isot,
		TPCSTATUS *	status )

Calculate weights for frames in SIF data based on true counts.

Weights are calculated from formula weight=(frame duration)^2 / (trues in a frame) which is applicable to non-decay corrected data. If weights are to be applied to decay-corrected data, provide function with isotope to use formula weight=(frame duration)^2 / (decay corrected trues in a frame) instead.

Reference: Mazoyer BM, Huesman RH, Budinger TF, Knittel BL. Dynamic PET data analysis. J Comput Assist Tomogr 1986; 10:645-653.

Formula is not applicable to initial frames where trues can be zero or even less, but weights for frames with no trues should still be high for fitting purposes; therefore, in those cases, weight is set to the maximum weight that was calculated from the data.

Weights are normalized to have an average of 1.0 for frames that have weight above zero.

See also

tacWCopy, tacWByFreq, tacWSampleNr, tacReadSIF, tacWeightModerate, tacWeightNorm, tacDecayCorrection, tacIsWeighted, tacSetWeights

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

sif	Pointer to SIF data, stored in TAC structure. Weights are added to the structure. Data must have at least two columns (TACs), which are assumed to represent prompts and randoms, to calculate trues as prompts - randoms. Counts in SIF are assumed to be not corrected for physical decay. Frame times are assumed to be in seconds.
isot	Isotope code, in case the weights are calculated for decay corrected data. Enter ISOTOPE_UNKNOWN if calculating weights for non-corrected data. SIF should contain isotope, too, but that is not used in this function.
status	Pointer to status data; enter NULL if not needed

Definition at line 363 of file tacw.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(%s)\n", __func__, isotopeName(isot));
  if(sif==NULL || sif->tacNr<1 || sif->sampleNr<1) {
    if(sif!=NULL && verbose>2) {
      printf("  sif.tacNr := %d\n  sif.sampleNr := %d\n", sif->tacNr, sif->sampleNr);
    }
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  sif->weighting=WEIGHTING_OFF;
  if(sif->tacNr<2 || !sif->isframe) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return TPCERROR_UNSUPPORTED;
  }
  /* Check times */
  double xmin, xmax;
  if(tacXRange(sif, &xmin, &xmax) || !(xmax>xmin)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
 
  /* Calculate trues */
  double trues[sif->sampleNr];
  for(int i=0; i<sif->sampleNr; i++) trues[i]=sif->c[0].y[i]-sif->c[1].y[i];
  if(verbose>2) {
    printf("\nTrues:\n");
    for(int i=0; i<sif->sampleNr; i++) printf("\t%g\n", trues[i]);
    fflush(stdout);
  }
  /* Correct for decay, if requested */
  if(isot!=ISOTOPE_UNKNOWN) {
    for(int i=0; i<sif->sampleNr; i++)
      trues[i]*=decayCorrectionFactorFromIsotope(isot, sif->x1[i]/60., (sif->x2[i]-sif->x1[i])/60.);
    if(verbose>2) {
      printf("\nTrues after decay correction:\n");
      for(int i=0; i<sif->sampleNr; i++) printf("\t%g\n", trues[i]);
      fflush(stdout);
    }
  }
  /* Calculate weights */
  double wMax=0.0;
  for(int i=0; i<sif->sampleNr; i++) {
    double fdur=sif->x2[i]-sif->x1[i]; if(fdur<0.0) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return TPCERROR_INVALID_X;
    }
    if(trues[i]>0.0) {
      sif->w[i]=fdur*fdur/trues[i];
      if(sif->w[i]>wMax) wMax=sif->w[i];
    } else {
      sif->w[i]=nan("");
    }
  }
  if(!(wMax>0.0)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  if(verbose>1) {printf("\twMax := %g\n", wMax); fflush(stdout);}
  /* Set maximum weight to frames which had no true counts */
  for(int i=0; i<sif->sampleNr; i++) if(!(trues[i]>0.0) && isnan(sif->w[i])) sif->w[i]=wMax;
  if(verbose>2) {
    printf("\nWeights before normalization:\n");
    for(int i=0; i<sif->sampleNr; i++) printf("\t%g\n", sif->w[i]);
    fflush(stdout);
  }
 
  /* Normalize weights */
  sif->weighting=WEIGHTING_ON_COUNTS;
  if(tacWeightNorm(sif, status)!=TPCERROR_OK) return(status->error);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by tacSetWeights().

tacIsWeighted()

int tacIsWeighted

(

TAC *

tac

)

Check if TAC contains weights as indicated by the 'weighting' field.

See also

tacSetWeights, tacWMove, sifWeight, tacWByFreq, imgHasWeights, tacIsX, tacIsSize

Returns

Returns 1 if weighted, and 0 if not or in case of an error.

Parameters

tac	Pointer to TAC structure.

Definition at line 24 of file tacw.c.

  {
  if(tac==NULL) return(0);
  if(tac->weighting==WEIGHTING_UNKNOWN || tac->weighting==WEIGHTING_OFF) return(0);
  if(tac->weighting<WEIGHTING_LAST) return(1);
  return(0);
}

Referenced by parAllocateWithTAC(), tacDecayCorrection(), tacSampleXRange(), tacSetWeights(), tacWeightModerate(), tacWeightNorm(), tacWMove(), tacWriteCSV(), tacWriteDFT(), tacWritePMOD(), tacWriteSheetIntoXML(), and tacWriteXML().

tacSetWeights()

int tacSetWeights	(	TAC *	tac,
		weights	weightMethod,
		int	weightNr,
		TPCSTATUS *	status )

Set TAC sample weights, based on given weighting scheme.

Weights are normalized to have average weight 1.0, with the sum of weights equals the weightNr or number of frames, whichever is smaller.

See also

tacIsWeighted, tacWSampleNr, tacWCopy, tacWMove, sifWeight, tacWByFreq, tacIsX

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

tac	Pointer to TAC structure. Data must be sorted by increasing sample times.
weightMethod	Weighting scheme: WEIGHTING_ON_COUNTS, WEIGHTING_ON_F, WEIGHTING_ON_FD, WEIGHTING_OFF. With WEIGHTING_UNKNOWN weightNr (below) is applied to weights and remaining weights are re-normalized.
weightNr	Number of samples (frames) to weight; if TAC contains more samples than this, then those are set to have zero weight. Enter a large number to give weight to all available samples.
status	Pointer to status data; enter NULL if not needed

Definition at line 462 of file tacw.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(tac, %d, %d)\n", __func__, weightMethod, weightNr); fflush(stdout);}
 
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(weightNr<1) weightMethod=WEIGHTING_OFF;
  if(weightNr>tac->sampleNr) weightNr=tac->sampleNr;
  int origSampleNr=tac->sampleNr;
 
  /* Get isotope from TAC data */
  isotope fisot=tacGetIsotope(tac);
  if(verbose>3) {printf("  isotope :=  %s\n", isotopeName(fisot)); fflush(stdout);}
  if(fisot==ISOTOPE_UNKNOWN) {
    if(weightMethod==WEIGHTING_ON_COUNTS || weightMethod==WEIGHTING_ON_FD) {
      tac->weighting=WEIGHTING_OFF;
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNKNOWN_ISOTOPE);
      return TPCERROR_UNKNOWN_ISOTOPE;
    }
  }
 
  if(tac->format==TAC_FORMAT_SIF && weightMethod==WEIGHTING_ON_GENERAL)
    weightMethod=WEIGHTING_ON_COUNTS;
 
  if(weightMethod==WEIGHTING_OFF) {
 
    if(verbose>2 && !tacIsWeighted(tac)) {
      printf("  Note: data did not contain weights.\n"); fflush(stdout);}
    tac->weighting=WEIGHTING_OFF;
    for(int i=0; i<tac->sampleNr; i++) tac->w[i]=0.0;
    for(int i=0; i<weightNr; i++) tac->w[i]=1.0;
    if(verbose>1) {printf("  weights removed.\n"); fflush(stdout);}
 
  } else if(weightMethod==WEIGHTING_ON_GENERAL) {
 
    if(verbose>0) {fprintf(stderr, "Error: cannot apply general weighting.\n"); fflush(stderr);}
    for(int i=0; i<tac->sampleNr; i++) tac->w[i]=0.0;
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
 
  } else if(weightMethod==WEIGHTING_UNKNOWN) {
 
    for(int i=weightNr; i<tac->sampleNr; i++) tac->w[i]=0.0;
 
    if(!tacIsWeighted(tac)) {
      if(verbose>2) {printf("  Note: data does not contain weights.\n"); fflush(stdout);}
      for(int i=0; i<weightNr; i++) tac->w[i]=1.0;
    } else {
      if(verbose>2) {printf("  Note: data contains weights.\n"); fflush(stdout);}
      tac->sampleNr=weightNr;
      int ret=tacWeightNorm(tac, status);
      tac->sampleNr=origSampleNr;
      if(ret!=TPCERROR_OK) return(ret);
    }
 
  } else if(weightMethod==WEIGHTING_ON_COUNTS) {
 
    if(verbose>2 && tacIsWeighted(tac)) {
      printf("  Note: data already contained weights.\n"); fflush(stdout);}
    tac->weighting=WEIGHTING_OFF;
    for(int i=0; i<tac->sampleNr; i++) tac->w[i]=0.0;
 
    if(tac->format==TAC_FORMAT_SIF) {
 
      /* TAC contains SIF data, just add weights */
      tac->sampleNr=weightNr;
      int ret=sifWeight(tac, fisot, status);
      tac->sampleNr=origSampleNr;
      if(ret!=TPCERROR_OK) return(ret);
      tac->weighting=WEIGHTING_ON_COUNTS;
      if(verbose>1) {printf("  weights added.\n"); fflush(stdout);}
 
    } else {
 
      double weight_moderate=100.0;
 
      /* Create temp SIF */
      if(!unitIsTime(tac->tunit) || !tacIsX(tac)) {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
        return TPCERROR_INVALID_X;
      }
      if(!unitIsRadioactivity(tac->cunit) && !unitIsRAConc(tac->cunit)) {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INCOMPATIBLE_UNIT);
        return TPCERROR_INCOMPATIBLE_UNIT;
      }
      TAC sif; tacInit(&sif);
      if(tacAllocate(&sif, weightNr, 3)!=TPCERROR_OK) {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INCOMPATIBLE_DATA);
        return TPCERROR_INCOMPATIBLE_DATA;
      }
      sif.sampleNr=weightNr;
      sif.tacNr=2; // 3rd contains trues that is not saved
      sif.format=TAC_FORMAT_SIF;
      sif.isframe=1; // SIF always contains frame start and end times
      sif.tunit=tac->tunit;
      sif.cunit=tac->cunit;
      tacSetIsotope(&sif, fisot);
      tacXCopy(tac, &sif, 0, sif.sampleNr-1);
      int ret=tacXUnitConvert(&sif, UNIT_SEC, status);
      if(ret!=TPCERROR_OK) {tacFree(&sif); return(ret);}
      for(int j=0; j<sif.sampleNr; j++) {
        double sum=0.0, wsum=0.0;
        for(int i=0; i<tac->tacNr; i++) {
          if(!isfinite(tac->c[i].y[j])) continue;
          double w=tac->c[i].size; if(!(w>0.0)) w=1.0;
          sum+=w*tac->c[i].y[j]; wsum+=w;
        }
        if(wsum>0.0) sif.c[0].y[j]=sum/wsum; else sif.c[0].y[j]=0.0;
      } // next time frame
      ret=tacDecayCorrection(&sif, fisot, 0, status);
      if(ret!=TPCERROR_OK) {tacFree(&sif); return(ret);}
      if(unitIsRadioactivity(sif.cunit)) {
        ret=tacYUnitConvert(&sif, UNIT_BQ, status);
      } else {
        ret=tacYUnitConvert(&sif, UNIT_BQ_PER_ML, status);
        for(int i=0; i<sif.sampleNr; i++) sif.c[0].y[i]*=1000.; // Counts from one litre
        sif.cunit=UNIT_BQ;
      }
      if(ret!=TPCERROR_OK) {tacFree(&sif); return(ret);}
      /* Multiply with frame duration (Bq/s -> Bq/frame) */
      for(int i=0; i<sif.sampleNr; i++) sif.c[0].y[i]*=(sif.x2[i]-sif.x1[i]);
      sif.cunit=UNIT_COUNTS;
      /* Set prompts, randoms, and trues */
      for(int i=0; i<sif.sampleNr; i++) {
        sif.c[2].y[i]=sif.c[0].y[i];
        sif.c[1].y[i]=0.0;
      }
      /* Calculate weights based on approximated SIF data */
      ret=sifWeight(&sif, fisot, status);
      if(ret==TPCERROR_OK) ret=tacWeightNorm(&sif, status);
      if(ret!=TPCERROR_OK) {tacFree(&sif); return(ret);}
 
      double prop=0.0; if(weight_moderate>1.0) prop=1.0/weight_moderate;
      ret=tacWeightModerate(&sif, prop, 0, 0, status);
      if(ret!=TPCERROR_OK) {tacFree(&sif); return(ret);}
 
      /* Copy weights into TAC */
      {
        int j;
        for(j=0; j<sif.sampleNr; j++) tac->w[j]=sif.w[j];
        for( ; j<tac->sampleNr; j++) tac->w[j]=0.0;
      }
      tacFree(&sif);
      tac->weighting=WEIGHTING_ON_COUNTS;
      if(verbose>1) {printf("  weights added.\n"); fflush(stdout);}
 
    }
 
  } else if(weightMethod==WEIGHTING_ON_F) {
 
    if(verbose>2 && tacIsWeighted(tac)) {
      printf("  Note: data already contained weights.\n"); fflush(stdout);}
    tac->weighting=WEIGHTING_OFF;
    for(int i=0; i<tac->sampleNr; i++) tac->w[i]=0.0;
    tac->sampleNr=weightNr;
    int ret=tacWByFreq(tac, ISOTOPE_UNKNOWN, status);
    if(ret==TPCERROR_OK) ret=tacWeightNorm(tac, status);
    tac->sampleNr=origSampleNr;
    if(ret!=TPCERROR_OK) return(ret);
    if(verbose>1) {printf("  weights added.\n"); fflush(stdout);}
 
  } else if(weightMethod==WEIGHTING_ON_FD) {
 
    if(verbose>2 && tacIsWeighted(tac)) {
      printf("  Note: data already contained weights.\n"); fflush(stdout);}
    tac->weighting=WEIGHTING_OFF;
    for(int i=0; i<tac->sampleNr; i++) tac->w[i]=0.0;
    tac->sampleNr=weightNr;
    int ret=tacWByFreq(tac, fisot, status);
    if(ret==TPCERROR_OK) ret=tacWeightNorm(tac, status);
    tac->sampleNr=origSampleNr;
    if(ret!=TPCERROR_OK) return(ret);
    if(verbose>1) {printf("  weights added.\n"); fflush(stdout);}
 
  } else {
    if(verbose>0) {fprintf(stderr, "Error: invalid weighting scheme.\n"); fflush(stderr);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacWByFreq()

int tacWByFreq	(	TAC *	tac,
		isotope	isot,
		TPCSTATUS *	status )

Add weight to TAC data based on sample frequency or time frame length.

Weights are scaled so that weight sum = sample number. Any existing weighting is overwritten.

See also

tacSetWeights, tacIsWeighted, sifWeight, tacWSampleNr, tacWCopy, tacWMove, tacWeightModerate

Returns

Return enum tpcerror (TPCERROR_OK when successful).

Parameters

tac	Pointer to TAC struct where weights will be added.

Note

Samples/frames must be sorted by sample time, but duplicate samples are allowed.

Parameters

isot	Isotope code, in case the weights are calculated for decay corrected data. Enter ISOTOPE_UNKNOWN if calculating weights for non-corrected data, or when decay correction is ignored.
status	Pointer to status data; enter NULL if not needed

Definition at line 134 of file tacw.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(%s)\n", __func__, isotopeName(isot));
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  
  tac->weighting=WEIGHTING_ON_GENERAL;
  if(tac->sampleNr==1) {tac->w[0]=1.0; return(TPCERROR_OK);}
 
  if(tac->isframe!=0) { /* weights based on frame lengths */
    tac->weighting=WEIGHTING_ON_F;
 
    for(int i=0; i<tac->sampleNr; i++) {
      tac->w[i]=tac->x2[i]-tac->x1[i]; 
    }
 
  } else { /* weights based on sample distance */
    tac->weighting=WEIGHTING_ON_F;
 
    int i, i1, i2;
    double t, t1, t2, f;
    for(i=0; i<tac->sampleNr; i++) {
      t=t1=t2=tac->x[i];
      /* Find the closest sample time before this one */
      for(i1=i; i1>=0; i1--) {t1=tac->x[i1]; if(t1<t) break;} 
      /* Find the closest sample time after this one */ 
      for(i2=i; i2<tac->sampleNr; i2++) {t2=tac->x[i2]; if(t2>t) break;} 
      /* Mean sample distance */
      f=0.0;
      if(t1<t) f+=t-t1; else f+=t2-t;
      if(t2>t) f+=t2-t; else f+=t-t1;
      f*=0.5; if(f<=0.0) f=1.0;
      tac->w[i]=f; 
    }
 
  }
 
  /* Account for decay correction, if requested */
  if(isot!=ISOTOPE_UNKNOWN) {
    tac->weighting=WEIGHTING_ON_FD;
    double lambda=lambdaFromIsotope(isot);
    if(tac->tunit==UNIT_SEC) lambda/=60.; else if(tac->tunit==UNIT_MSEC) lambda/=60000.;
    if(tac->isframe!=0) {
      for(int i=0; i<tac->sampleNr; i++)
        tac->w[i]*=decayCorrectionFactorFromLambda(-lambda, tac->x1[i], (tac->x2[i]-tac->x1[i])); 
    } else {
      for(int i=0; i<tac->sampleNr; i++)
        tac->w[i]*=decayCorrectionFactorFromLambda(-lambda, tac->x[i], 0.0); 
    }
  }
 
  /* Scale weights so that sum of weights equals sample number */
  double sumw=0.0;
  for(int i=0; i<tac->sampleNr; i++) if(isfinite(tac->w[i])) sumw+=tac->w[i];
  if(!(sumw>0.0)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  sumw/=(double)tac->sampleNr;
  for(int i=0; i<tac->sampleNr; i++) tac->w[i]/=sumw;
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by tacSetWeights().

tacWCopy()

int tacWCopy	(	TAC *	tac1,
		TAC *	tac2,
		int	i1,
		int	i2 )

Copy weights from one TAC structure to another.

See also

tacXCopy, tacWMove, sifWeight, tacWByFreq, tacIsWeighted, tacSetWeights

Author

Vesa Oikonen

Returns

Returns TPCERROR status.

Parameters

tac1	Pointer to source TAC structure.
tac2	Pointer to target TAC structure.
i1	Index of the first sample to copy.
i2	Index of the last sample to copy.

Definition at line 41 of file tacw.c.

  {
  if(tac1==NULL || tac2==NULL) return(TPCERROR_FAIL);
  if(i1<0 || i1>i2) return(TPCERROR_FAIL);
  if(i2>=tac1->_sampleNr || i2>=tac2->_sampleNr) return(TPCERROR_FAIL);
  
  for(int i=i1; i<=i2; i++) {
    tac2->w[i]=tac1->w[i];
  }
  tac2->weighting=tac1->weighting;
  return(TPCERROR_OK);
}

Referenced by tacDuplicate(), and tacExtract().

tacWeightModerate()

int tacWeightModerate	(	TAC *	tac,
		const double	minprop,
		const int	doZeroes,
		const int	doNaNs,
		TPCSTATUS *	status )

Moderate weights so that the min weight is not less than given proportion of the maximum weight.

After moderation, all weights are normalized to have and average of one.

See also

sifWeight, tacWeighNorm, tacWCopy, tacWByFreq, tacWSampleNr, tacReadSIF

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

tac	Pointer to TAC structure (which can contain SIF or TAC data).
minprop	The minimum proportion of frame weight to maximum frame weight. Must be less than one. If proportion is lower than this limit, then weight is set to proportion times maximum. Setting minimum proportion to zero or less does not change any weights, but weights are still normalized. By default, zero weights or NaNs are not changed.
doZeroes	Set zero weights, too, to the minimum proportion of maximum frame weight.
doNaNs	Set NaN weights, too, to the minimum proportion of maximum frame weight.
status	Pointer to status data; enter NULL if not needed

Definition at line 277 of file tacw.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(%g, %d, %d)\n", __func__, minprop, doZeroes, doNaNs);
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(!isfinite(minprop) || minprop>=1.0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  if(!tacIsWeighted(tac)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_WEIGHTS);
    return TPCERROR_NO_WEIGHTS;
  }
 
  double wMax=nan("");
  for(int i=0; i<tac->sampleNr; i++) if(isfinite(tac->w[i])) {
    if(!isfinite(wMax)) {wMax=tac->w[i]; continue;}
    if(tac->w[i]>wMax) wMax=tac->w[i];
  }
  if(verbose>2) {printf("weight_maximum := %g\n", wMax); fflush(stdout);}
  if(!isfinite(wMax) || wMax<=0.0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_WEIGHTS);
    return TPCERROR_NO_WEIGHTS;
  }
 
  double wLimit=minprop*wMax;
  if(verbose>2) {printf("weight_limit := %g\n", wLimit); fflush(stdout);}
 
  int n=0;
  if(wLimit>0.0) {
    for(int i=0; i<tac->sampleNr; i++) {
      if(tac->w[i]==0.0) {
        if(doZeroes) {tac->w[i]=wLimit; n++;}
        continue;
      }
      if(isnan(tac->w[i])) {
        if(doNaNs) {tac->w[i]=wLimit; n++;}
        continue;
      }
      if(tac->w[i]<wLimit) {tac->w[i]=wLimit; n++;}
    }
  }
  if(verbose>1 && minprop>0.0) {printf("  %d weight(s) set to %g.\n", n, wLimit); fflush(stdout);}
 
  return(tacWeightNorm(tac, status));
}

Referenced by tacSetWeights().

tacWeightNorm()

int tacWeightNorm	(	TAC *	tac,
		TPCSTATUS *	status )

Scale weights so that average weight is 1.0, and the sum of weights equals the nr of frames.

See also

tacWeightModerate,sifWeight, tacWCopy, tacWByFreq, tacWSampleNr, tacReadSIF

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

tac	Pointer to TAC structure (which can contain SIF or TAC data).
status	Pointer to status data; enter NULL if not needed

Definition at line 237 of file tacw.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(!tacIsWeighted(tac)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_WEIGHTS);
    return TPCERROR_NO_WEIGHTS;
  }
 
  double wSum=0.0; int wNr=0;
  for(int i=0; i<tac->sampleNr; i++)
    if(isfinite(tac->w[i]) && tac->w[i]>0.0) {wSum+=tac->w[i]; wNr++;}
  if(wNr==0 || !(wSum>0.0)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  double wMean=wSum/(double)wNr;
  for(int i=0; i<tac->sampleNr; i++) 
    if(isfinite(tac->w[i]) && tac->w[i]>0.0) tac->w[i]/=wMean;
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}

Referenced by sifWeight(), tacDecayCorrection(), tacSetWeights(), and tacWeightModerate().

tacWMove()

int tacWMove	(	TAC *	tac,
		int	ow,
		TPCSTATUS *	status )

Identify column containing weights in TAC structure, and move weights to the correct place.

TAC name starting with string 'weight' is identified as weight column.

Sets tac->weighting to WEIGHTING_ON or WEIGHTING_OFF, depending on whether weight column was found or not.

See also

tacWCopy, tacSetWeights

Author

Vesa Oikonen

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

tac	Pointer to TAC structure.
ow	Overwrite (1) or do not overwrite (0) existing weights.
status	Pointer to status data; enter NULL if not needed

Definition at line 75 of file tacw.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  int ret;
 
  if(tac->weighting==WEIGHTING_UNKNOWN) tac->weighting=WEIGHTING_OFF;
 
  /* Search for the weight column */
  int wc=-1;
  for(int i=0; i<tac->tacNr; i++) {
    if(strncasecmp(tac->c[i].name, "WEIGHT", 6)!=0) continue;
    if(wc<0) wc=i;
    else if(ow==0) {
      /* Error if found more than one */
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_DUPLICATE_DATA);
      return TPCERROR_DUPLICATE_DATA;
    }
  }
  /* Error if not found */
  if(wc<0) {
    if(ow) tac->weighting=WEIGHTING_OFF;
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_WEIGHTS);
    return TPCERROR_NO_WEIGHTS;
  }
  /* Error if overwriting not allowed */
  if(ow==0 && tacIsWeighted(tac)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_DUPLICATE_DATA);
    return TPCERROR_DUPLICATE_DATA;
  }
  /* Copy weights to their right place */
  for(int i=0; i<tac->sampleNr; i++) tac->w[i]=tac->c[wc].y[i];
  tac->weighting=WEIGHTING_ON_GENERAL;
  /* Then delete the original column */
  ret=tacDeleteTACC(tac, wc);
  statusSet(status, __func__, __FILE__, __LINE__, ret);
  return(ret);
}

Referenced by tacReadCarimasTxt(), tacReadCSV(), tacReadDFT(), and tacReadPMOD().

tacWSampleNr()

unsigned int tacWSampleNr

(

TAC *

tac

)

Get the number of samples in TAC that have weight > 0. Missing (NaN) sample values are included as long as weight is not missing. If weights are not set, then nr of all samples is returned.

Returns

The number of samples with weight above zero.

See also

tacWByFreq, aicSS, parFreeNr, tacWeightModerate, tacIsWeighted, tacSetWeights

Parameters

tac	Pointer to the TAC struct.

Definition at line 219 of file tacw.c.

  {
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) return((unsigned int)0);
  if(tac->weighting==WEIGHTING_OFF) return((unsigned int)tac->sampleNr);
  unsigned int n=0;
  for(int i=0; i<tac->sampleNr; i++) if(tac->w[i]>0.0) n++;
  if(n==0 && tac->weighting==WEIGHTING_UNKNOWN) n=(unsigned int)tac->sampleNr;
  return(n);
}