parfit.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 "tpcmodels.h"
#include "tpccsv.h"
#include "tpcift.h"
/*****************************************************************************/
#include "tpcpar.h"
/*****************************************************************************/
 
/*****************************************************************************/
int parWriteFIT(
  PAR *par,
  FILE *fp,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(fp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
  if(verbose>0) printf("%s():\n", __func__);
  if(par==NULL || par->tacNr<1 || par->parNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  int i, j, n;
  char tmp[256];
 
  /* Write magic string and fit program name */
  i=iftFindKey(&par->h, "program", 0);
  if(i<0 || !strnlen(par->h.item[i].value, 10)) strcpy(tmp, "libtpcpar");
  else strlcpy(tmp, par->h.item[i].value, 256);
  n=fprintf(fp, "%-11.11s %s\n", "FIT1", tmp);
  if(n<5) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /* Write fit date and time */
  i=iftFindKey(&par->h, "analysis_time", 0);
  if(i>=0) strlcpy(tmp, par->h.item[i].value, 256); else strcpy(tmp, "");
  fprintf(fp, "Date:\t%s\n", tmp);
 
  /* Write the name of the original datafile */
  i=iftFindKey(&par->h, "datafile", 0);
  if(i>=0) fprintf(fp, "Data file:\t%s\n", par->h.item[i].value);
  else fprintf(fp, "Data file:\t\n");
 
  /* Write the study number */
  i=iftFindKey(&par->h, "studynr", 0);
  if(i<0) i=iftFindKey(&par->h, "study_number", 0);
  if(i>=0) fprintf(fp, "Study number:\t%s\n", par->h.item[i].value);
 
  /* Write the concentration (y) unit of the original data */
  i=iftFindKey(&par->h, "unit", 0);
  if(i<0) i=iftFindKey(&par->h, "calibration_unit", 0);
  if(i>=0) fprintf(fp, "Data unit:\t%s\n", par->h.item[i].value);
  else fprintf(fp, "Data unit:\t\n");
 
  /* Write the time (x) unit of the original data */
  i=iftFindKey(&par->h, "timeunit", 0);
  if(i<0) i=iftFindKey(&par->h, "time_unit", 0);
  if(i>=0) fprintf(fp, "Time unit:\t%s\n", par->h.item[i].value);
  else fprintf(fp, "Time unit:\t\n");
 
  /* Write the tacNr to be saved */
  fprintf(fp, "Nr of VOIs:\t%d\n", par->tacNr);
 
  /* Write the Fit title */
  fprintf(fp, "%s %s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n",
    "Region", "Plane", "Start", "End", "dataNr", "WSS", "parNr", "Type", "Parameters");
 
  /* Write the individual fits */
  for(i=0; i<par->tacNr; i++) {
    /* Tac id */
    if(strlen(par->r[i].name)>0) {
      if(!roinameSubpart(par->r[i].name, "_- ", 0, tmp, 7)) strcpy(tmp, ".");
      fprintf(fp, "%s ", tmp);
      if(!roinameSubpart(par->r[i].name, "_- ", 1, tmp, 7)) strcpy(tmp, ".");
      fprintf(fp, "%s ", tmp);
      if(!roinameSubpart(par->r[i].name, "_- ", 2, tmp, 7)) strcpy(tmp, ".");
      fprintf(fp, "%s", tmp);
    } else {
      fprintf(fp, "tac%d . .", 1+i);
    }
    /* Time range, sample nr, (W)SS, parameter nr, function id */
    fprintf(fp, "\t%.3f\t%.3f\t%d\t%.2E\t%d\t%04d",
       par->r[i].start, par->r[i].end, par->r[i].dataNr,
       par->r[i].wss, par->parNr, modelOldId(par->r[i].model) );
    /* Function parameters */
    for(j=0; j<par->parNr; j++) fprintf(fp, "\t%.6E", par->r[i].p[j]);
    n=fprintf(fp, "\n");
  } // next fit
  if(n==0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
  
  /* Quit */
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int parReadFIT(
  PAR *par,
  CSV *csv,
  IFT *ift,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(par==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  parFree(par);
  
  if(csv==NULL || ift==NULL || csv->row_nr<1 || csv->col_nr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }  
  if(verbose>0) printf("%s()\n", __func__);
 
  if(verbose>10) {
    printf("\n CSV contents shown with separator ; \n\n");
    char sep=csv->separator;
    csv->separator=';'; csvWrite(csv, 0, stdout, NULL); csv->separator=sep;
    printf("\n IFT contents \n\n");
    iftWrite(ift, stdout, NULL);
    printf("\n");
  }
 
  int i, n, ret, titlerow, tunit=UNIT_UNKNOWN;
  char *cptr, tmp[256];
 
  /* TAC nr */
  i=iftFindKey(ift, "Nr of VOIs", 0);
  if(i<0) i=iftFindKey(ift, "voiNr", 0);
  if(i<0) i=iftFindKey(ift, "tacNr", 0);
  if(i<0) {
    parFree(par);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return TPCERROR_UNSUPPORTED;
  }
  int tacNr=atoi(ift->item[i].value);
 
  /* Find (and ignore) the parameter title line */
  i=0; if(verbose>2) printf("estimating parNr\n"); 
  while(i<csv->nr) {
    i=csvSearchField(csv, "Region", i); if(i<0) break; 
    if(i>=0 && csv->c[i].col==0) break;
    i++;
  }
  if(i<0 || i==csv->nr) {
    parFree(par);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  }
  titlerow=csv->c[i].row; if(verbose>3) printf("titlerow := %d\n", titlerow);
#if(1)
  /* Parameter nr from the row with most columns */
  int parNr=0;
  for(int r=0; r<tacNr; r++) {
    int n=csvRowLength(csv, r+titlerow+1); if(verbose>5) printf("n=%d at row %d\n", n, r+titlerow+1);
    if(n>parNr) parNr=n;
  }
  if(verbose>4) printf("  initial_parNr := %d\n", parNr);
#else
  /* Parameter nr from the next row */
  int parNr=csvRowLength(csv, titlerow+1);
#endif
  if(csv->separator!=' ') parNr-=7; else parNr-=9;
  if(verbose>2) {
    printf("  parNr := %d\n", parNr);
    printf("  tacNr := %d\n", tacNr);
  }
  if(parNr<1 || tacNr<1) {
    parFree(par);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  }
  /* Allocate space for parameters */
  ret=parAllocate(par, parNr, tacNr); if(ret!=TPCERROR_OK) {
    parFree(par);
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return ret;
  }
  par->parNr=parNr; par->tacNr=tacNr;
 
  /* Check magic string and get program name from the first non-comment IFT field */
  i=0; 
  while(ift->item[i].comment && i<ift->keyNr && ift->item[i].key!=NULL) i++;
  if(verbose>3) {
    printf("ift->item[%d].value := '%s'\n", i, ift->item[i].value);
  }
  if(i==ift->keyNr || strncasecmp(ift->item[i].value, "FIT1", 4)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return TPCERROR_UNSUPPORTED;
  }
  cptr=ift->item[i].value+4;
  if(strncpyCleanSpaces(tmp, cptr, 256)<1) {
    parFree(par);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return TPCERROR_UNSUPPORTED;
  }
  iftPut(&par->h, "program", tmp, 0, NULL);
 
  /* Get analysis time from IFT struct */
  i=iftFindKey(ift, "Date", 0);
  if(i>=0) iftPut(&par->h, "analysis_time", ift->item[i].value, 0, NULL);
 
  /* Filename of the original data */
  i=iftFindKey(ift, "Data file", 0);
  if(i>=0) iftPut(&par->h, "datafile", ift->item[i].value, 0, NULL);
 
  /* Study number */
  i=iftFindKey(ift, "Study number", 0);
  if(i<0) i=iftFindKey(ift, "studynr", 0);
  if(i>=0) iftPut(&par->h, "studynr", ift->item[i].value, 0, NULL);
 
  /* Units of the original data */
  i=iftFindKey(ift, "Data unit", 0);
  if(i>=0) iftPut(&par->h, "unit", ift->item[i].value, 0, NULL);
  i=iftFindKey(ift, "Time unit", 0);
  if(i>=0) {
    iftPut(&par->h, "timeunit", ift->item[i].value, 0, NULL);
    tunit=unitIdentify(ift->item[i].value);
  }
 
  /* Set parameter names */
  for(i=0; i<parNr; i++) sprintf(par->n[i].name, "p%d", 1+i);
 
  /* Read TAC parameters */
  ret=0; if(verbose>2) {printf("reading individual parameter values...\n"); fflush(stdout);}
  for(int j=0; j<tacNr; j++) {
    n=csvRowLength(csv, 1+j+titlerow); if(csv->separator!=' ') n-=7; else n-=9;
    if(n>parNr) {ret++; break;}
    i=0; n=1+j+titlerow;
    /* copy tac name */
    strlcpy(par->r[j].name, csvCell(csv, n, i), MAX_TACNAME_LEN);
    if(csv->separator==' ') {
      i++; 
      strlcat(par->r[j].name, " ", MAX_TACNAME_LEN);
      strlcat(par->r[j].name, cptr=csvCell(csv, n, i), MAX_TACNAME_LEN);
      i++;
      strlcat(par->r[j].name, " ", MAX_TACNAME_LEN);
      strlcat(par->r[j].name, cptr=csvCell(csv, n, i), MAX_TACNAME_LEN);
    }
    /* fit start and end times */
    i++; par->r[j].start=atofVerified(csvCell(csv, n, i));
    i++; par->r[j].end=atofVerified(csvCell(csv, n, i));
    if(tunit==UNIT_SEC) {par->r[j].start/=60.; par->r[j].end/=60.;}
    /* nr of fitted samples */
    i++; par->r[j].dataNr=atoi(csvCell(csv, n, i));
    /* (W)SS */
    i++; par->r[j].wss=atofVerified(csvCell(csv, n, i));
    /* nr of parameters for this TAC; this is not the same as the nr of fitted (free) 
       parameters that would be stored in PAR struct */
    i++; int tacParNr=atoi(csvCell(csv, n, i));
    par->r[j].fitNr=0; 
    /* function id */
    i++; par->r[j].model=modelOld2New(atoi(csvCell(csv, n, i)));
    /* function parameter values */
#if(1)
    int k; // read the ones that actually are saved... and then fill the rest with NaNs
    for(k=0; k<tacParNr; k++) {
      i++; par->r[j].p[k]=atofVerified(csvCell(csv, n, i));
    }
    for(; k<parNr; k++) par->r[j].p[k]=nan("");
#else
    for(int k=0; k<parNr; k++) {
      i++; par->r[j].p[k]=atofVerified(csvCell(csv, n, i));}
#endif
  }
  if(ret>0) {
    parFree(par);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  }
 
  /* Fix region names (mainly, remove dots) */
  {
    int i, tac;
    char fnsp[3][MAX_TACNAME_LEN+1];
    for(tac=0; tac<par->tacNr; tac++) {
      for(i=0; i<3; i++) {
        if(roinameSubpart(par->r[tac].name, "_ ", i, fnsp[i], MAX_TACNAME_LEN)==NULL ||
           strcmp(fnsp[i], ".")==0)
          strcpy(fnsp[i], "");
      }
      strcpy(par->r[tac].name, fnsp[0]);
      if(strlen(fnsp[1]) || strlen(fnsp[2])) strcat(par->r[tac].name, "_");
      if(strlen(fnsp[1])) strcat(par->r[tac].name, fnsp[1]);
      if(strlen(fnsp[2])) {strcat(par->r[tac].name, "_"); strcat(par->r[tac].name, fnsp[2]);}
    }
  }
 
  /* Set format in struct */
  par->format=PAR_FORMAT_FIT;
  
  /* Quit */
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/