dftio.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include "tpcift.h"
#include "tpcisotope.h"
#include "tpccsv.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpctac.h"
/*****************************************************************************/
 
/*****************************************************************************/
int tacNameSplit(
  const char *rname,
  char *name1,
  char *name2,
  char *name3,
  unsigned int max_name_len
) {
  /* Initiate subnames */
  if(name1!=NULL) strcpy(name1, "");
  if(name2!=NULL) strcpy(name2, "");
  if(name3!=NULL) strcpy(name3, "");
 
  /* Check input */
  if(max_name_len<1) return 0;
  if(rname==NULL || *rname=='\0') return 0;
 
  /* Work with a copy of TAC name */
  char *temp=strdup(rname);
 
  /* Remove any quotes and trailing and initial spaces */
  if(strClean(temp)) {free(temp); return 0;}
  //printf("temp='%s'\n", temp);
 
  /* Allocate memory for local subnames */
  int fullLen=(int)strlen(rname);
  char sname[3][fullLen+1];
 
  /* Determine the subpart delimiters so that subpart number would be optimal */
  char space[32];
  strcpy(space, " \t"); 
  int n=strTokenNr(temp, space);
  if(n<3) {strcat(space, "_"); n=strTokenNr(temp, space);}
  if(n<3) {strcat(space, "-"); n=strTokenNr(temp, space);}
  //printf("n=%d space='%s'\n", n, space);
 
  /* get subparts */
  char *cptr=temp; int i;
  n=0; i=0;
  while(n<3 && *cptr!='\0') {
    //printf("  cptr='%s' n=%d i=%d\n", cptr, n, i);
    /* If delimiter */
    if(strspn(cptr, space)>0) {
      /* then end this subpart */
      sname[n][i]='\0';
      /* and prepare for the next */ 
      n++; i=0;
      cptr++; continue;
    }
    sname[n][i]=*cptr;
    cptr++; i++;
  } 
  if(n<3) {sname[n][i]='\0'; n++;}
  for( ; n<3 ;n++) strcpy(sname[n], "");
 
  for(i=0; i<3; i++) {
    /* If name was '.' then replace it with empty string */
    if(strcmp(sname[i], ".")==0) strcpy(sname[i], "");
    //printf("         -> '%s'\n", sname[i]);
  }
 
  /* If subname2 is empty and subname1 is too long, then divide it */
  if(sname[1][0]=='\0' && strlen(sname[0])>max_name_len)
    strlcpy(sname[1], &sname[0][max_name_len], fullLen+1);
  /* If subname3 is empty and subname2 is too long, then divide it */
  if(sname[2][0]=='\0' && strlen(sname[1])>max_name_len)
    strlcpy(sname[2], &sname[1][max_name_len], fullLen+1);
 
  /* Fill the output strings */
  char *nptr;
  for(i=n=0; i<3; i++) {
    //printf("         -> '%s'\n", sname[i]);
    if(i==0) nptr=name1; else if(i==1) nptr=name2; else nptr=name3;
    if(strlen(sname[i])>0) n=i+1;
    if(nptr==NULL) continue;
    strlcpy(nptr, sname[i], max_name_len+1);
  }
 
  free(temp);
  return n;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacWriteDFT(
  TAC *tac,
  FILE *fp,
  int extra,
  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(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(verbose>0) printf("%s()\n", __func__);
 
  int ret, prec=6;
 
  /* Make sure that TAC names are available */
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  tacEnsureNames(tac);
 
  /* Create local copy of headers, which will be edited here and written
     into DFT file, to preserve the original headers with TAC struct */
  if(verbose>99) iftWrite(&tac->h, stdout, NULL);
  IFT localh; iftInit(&localh);
  ret=iftDuplicate(&tac->h, &localh);
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return ret;
  }
 
  /* Writing obligatory header lines */
  if(verbose>1) printf("writing obligatory title lines\n");
  int n;
  char subname[7];
 
  /* 1st line with filetype identification string and region names */
  n=fprintf(fp, "DFT1"); if(n<1) {
    iftFree(&localh);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
  for(int ri=0; ri<tac->tacNr; ri++) {
    n=tacNameSplit(tac->c[ri].name, subname, NULL, NULL, 6);
    fprintf(fp, "\t%s", subname);
  }
  if(tacIsWeighted(tac)) fprintf(fp, "\t%s", "weight");
  n=fprintf(fp, "\n"); if(n<1) {
    iftFree(&localh);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /* 2nd line with study identification and 2nd subname */
  if(verbose>2) printf("  writing 2nd title line\n");
  char tmp[128];
  n=tacGetHeaderStudynr(&localh, tmp, status);
  if(n!=TPCERROR_OK) {
    strcpy(tmp, "."); if(verbose>5) printf("no studynr\n");
  } else tacSetHeaderStudynr(&localh, ""); // delete studynr from header
  n=fprintf(fp, "%s", tmp); if(n<1) {
    iftFree(&localh);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
  for(int ri=0; ri<tac->tacNr; ri++) {
    n=tacNameSplit(tac->c[ri].name, NULL, subname, NULL, 6);
    if(n<2 || !subname[0]) strcpy(subname, ".");
    fprintf(fp, "\t%s", subname);
  }
  if(tacIsWeighted(tac)) fprintf(fp, "\t%s", ".");
  n=fprintf(fp, "\n"); if(n<1) {
    iftFree(&localh);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /* 3rd line with calibration unit and 3rd subname */
  if(verbose>2) printf("  writing 3rd title line\n");
  tacSetHeaderUnit(&localh, UNIT_UNKNOWN); // delete any unit(s) from header
  if(tac->cunit==UNIT_UNKNOWN) strcpy(tmp, ".");
  else strcpy(tmp, unitName(tac->cunit));
  n=fprintf(fp, "%s", tmp); if(n<1) {
    iftFree(&localh);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
  for(int ri=0; ri<tac->tacNr; ri++) {
    n=tacNameSplit(tac->c[ri].name, NULL, NULL, subname, 6);
    if(n<3 || !subname[0]) strcpy(subname, ".");
    fprintf(fp, "\t%s", subname);
  }
  if(tacIsWeighted(tac)) fprintf(fp, "\t%s", ".");
  n=fprintf(fp, "\n"); if(n<1) {
    iftFree(&localh);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /* 4th line with x (time) type & unit and region volumes */
  if(verbose>2) printf("  writing 4th title line\n");
  if(unitIsDistance(tac->tunit)) strcpy(tmp, "Distance");
  else strcpy(tmp, "Time");
  if(tac->isframe!=0) strcat(tmp, "s");
  n=fprintf(fp, "%s (%s)", tmp, unitName(tac->tunit)); if(n<1) {
    iftFree(&localh);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
  for(int ri=0; ri<tac->tacNr; ri++) {
    if(!isnan(tac->c[ri].size)) sprintf(tmp, "%.*e", prec, tac->c[ri].size);
    else strcpy(tmp, ".");
    fprintf(fp, "\t%s", tmp);
  }
  if(tacIsWeighted(tac)) fprintf(fp,"\t%s", ".");
  n=fprintf(fp, "\n"); if(n<1) {
    iftFree(&localh);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /* Write data */
  if(verbose>1) printf("writing data table\n");
  double v;
  n=10;
  for(int fi=0; fi<tac->sampleNr; fi++) {
    /* x or  x1 and x2) */
    if(tac->isframe==0) v=tac->x[fi]; else v=tac->x1[fi];
    if(isnan(v)) n=fprintf(fp, "."); else n=fprintf(fp, "%.5f", v);
    if(n<1) break;
    if(tac->isframe!=0) {
      v=tac->x2[fi];
      if(isnan(v)) n=fprintf(fp, "\t."); else n=fprintf(fp, "\t%.5f", v);
      if(n<1) break;
    }
    /* Concentrations (y values) */
    for(int ri=0; ri<tac->tacNr; ri++) {
      if(isnan(tac->c[ri].y[fi])) n=fprintf(fp, "\t.");
      else n=fprintf(fp, "\t%.*e", prec, tac->c[ri].y[fi]);
      if(n<1) break;
    }
    if(n<1) break;
    /* Weight */
    if(tacIsWeighted(tac)) {
      if(isnan(tac->w[fi])) n=fprintf(fp, "\t.");
      else n=fprintf(fp, "\t%.*e", prec, tac->w[fi]);
      if(n<1) break;
    }
    if(n<1) break;
    n=fprintf(fp, "\n"); if(n<1) break;
  } // next sample
  if(n<1) {
    iftFree(&localh);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
  
  /* Write extra header, if requested */
  if(extra && localh.keyNr>0) {
    int ret=iftWrite(&localh, fp, status);
    if(ret!=TPCERROR_OK) {iftFree(&localh); return ret;}
  }
 
  iftFree(&localh);
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacReadDFT(
  TAC *tac,
  CSV *csv,
  IFT *hdr,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(tac==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  tacFree(tac);
  if(csv==NULL || csv->row_nr<5 || csv->col_nr<2) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }  
  if(verbose>0) printf("%s()\n", __func__);
 
  /* Verify the magic number first */
  if(strncasecmp(csv->c[0].content, "DFT", 3)!=0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  }
 
  int ret, i, j, m, n, isframe=0, is2cunit=0, is2tunit=0, tunit;
  char *cptr;
 
  /* Get some important information from title lines first */
  for(i=0; i<csv->nr; i++) if(csv->c[i].col==0 && csv->c[i].row==3) break;
  if(i==csv->nr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  }
  /* Do we have one or two x columns */
  // note that time unit may or may not be included in this field, depending
  // on field separator
  if(strncasecmp(csv->c[i].content, "Times", 5)==0) isframe=1;
  else if(strncasecmp(csv->c[i].content, "Time", 4)==0) isframe=0;
  else if(strncasecmp(csv->c[i].content, "Distances", 9)==0) isframe=1;
  else if(strncasecmp(csv->c[i].content, "Distance", 8)==0) isframe=0;
  else {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  }
  if(verbose>2) printf("  isframe=%d\n", isframe);
  /* Get time unit from first or second field, inside parens */
  cptr=strchr(csv->c[i].content, '(');
  if(cptr==NULL) cptr=strchr(csv->c[i].content, '[');
  if(cptr==NULL) {cptr=strchr(csv->c[i+1].content, '('); is2tunit=1;}
  if(cptr==NULL) {cptr=strchr(csv->c[i+1].content, '['); is2tunit=1;}
  if(cptr==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  } else {
    char *tmp=strdup(cptr); strCleanPars(tmp); tunit=unitIdentify(tmp);
    free(tmp);
  }
 
  /* Allocate memory for TAC data */
  m=csv->row_nr-4;
  n=csv->col_nr-1; if(isframe) n--;
  if(n<1 || m<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  }
  ret=tacAllocate(tac, m, n);
  statusSet(status, __func__, __FILE__, __LINE__, ret);
  if(ret!=TPCERROR_OK) return ret;
  tac->tacNr=n; tac->sampleNr=m; // these may be refined later
  tac->isframe=isframe;
  tac->tunit=tunit;
 
  /* Copy header to TAC struct */
  iftDuplicate(hdr, &tac->h);
 
  /* Read first title line, jumping over magic number */
  n=0;
  for(i=1; i<csv->nr; i++) {
    if(csv->c[i].row!=0 || n>=tac->tacNr) break;
    /* Copy TAC name */
    strlcpy(tac->c[n].name, csv->c[i].content, MAX_TACNAME_LEN+1);
    if(strcmp(tac->c[n].name, ".")==0) strcpy(tac->c[n].name, "");
    n++;
  }
  if(n<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  }
  tac->tacNr=n;
  if(verbose>2)
    printf("  tacNr=%d\n  sampleNr=%d\n", tac->tacNr, tac->sampleNr);
 
  /* Get study number */
  for(i=0; i<csv->nr; i++) if(csv->c[i].row==1 && csv->c[i].col==0) {
    if(strlen(csv->c[i].content)>0 && strcmp(csv->c[i].content, ".")!=0)
      tacSetHeaderStudynr(&tac->h, csv->c[i].content);
    break;
  }
  
  /* Get concentration unit from the first item on line 3 */
  for(i=0; i<csv->nr; i++) if(csv->c[i].row==2 && csv->c[i].col==0) break;
  tac->cunit=unitIdentify(csv->c[i].content);
  /* Carimas may (erroneously) write conc unit in two columns, with actual 
     unit in 2nd column inside parenthesis */
  if(strInPars(csv->c[i+1].content)
     && csvRowLength(csv, 2)==csvRowLength(csv, 1)+1)
  {
    if(verbose>2) printf("  concentration unit spans two columns\n");
    is2cunit=1;
    if(tac->cunit==UNIT_UNKNOWN) {
      char *tmp=strdup(csv->c[i+1].content);
      strCleanPars(tmp); tac->cunit=unitIdentify(tmp);
      free(tmp);
    }
  }
 
  
  /* Catenate additional TAC names if available on lines 2 and 3 */
  for(int ri=0; ri<tac->tacNr; ri++) {
    for(i=0; i<csv->nr; i++) if(csv->c[i].row==1 && csv->c[i].col==ri+1) break;
    for(j=i; j<csv->nr; j++) if(csv->c[j].row==2 && csv->c[j].col==ri+1) break;
    if(i==csv->nr || j==csv->nr) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
      return TPCERROR_INVALID_FORMAT;
    }
    if(is2cunit) j++; // add one in case cunit spans two columns
    m=strlen(csv->c[i].content); if(strcmp(csv->c[i].content, ".")==0) m=0;
    n=strlen(csv->c[j].content); if(strcmp(csv->c[j].content, ".")==0) n=0;
    if((m+n)>0 && (strlen(tac->c[ri].name)<MAX_TACNAME_LEN-1))
      strcat(tac->c[ri].name, "_");
    if(m>0 && (m+strlen(tac->c[ri].name)<MAX_TACNAME_LEN))
      strcat(tac->c[ri].name, csv->c[i].content);
    if(n>0 && (strlen(tac->c[ri].name)<MAX_TACNAME_LEN-1))
      strcat(tac->c[ri].name, "_");
    if(n>0 && (n+strlen(tac->c[ri].name)<MAX_TACNAME_LEN))
      strcat(tac->c[ri].name, csv->c[j].content);
  }
 
  /* Read sizes from 4th line */
  n=csvRowLength(csv, 3);
  for(i=0; i<csv->nr; i++) if(csv->c[i].row==3) {
    m=csv->c[i].col-1; if(is2tunit) m--; 
    if(m>=0 && m<tac->tacNr) tac->c[m].size=atofVerified(csv->c[i].content);
  }
#if(0)
  for(i=0, m=-1, n=0; i<csv->nr; i++) if(csv->c[i].row==3) {n++; if(m<0) m=i;}
  if(n==tac->tacNr+1) m+=1;
  else if(n==tac->tacNr+2) m+=2;
  else {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  }
  for(int ri=0; ri<tac->tacNr; ri++)
    tac->c[ri].size=atofVerified(csv->c[m++].content);
#endif
 
  /* Copy x and y data from CSV into TAC struct */
  int fi=0, ri=0, oknr=0;
  ret=0;
  double v;
  for(i=0; i<csv->nr; i++) if(csv->c[i].row==4) break;
  for(; i<csv->nr; i++) {
    if(verbose>10)
      printf("i=%d\trow=%d\tcol=%d\n", i, csv->c[i].row, csv->c[i].col);
    fi=csv->c[i].row-4; if(fi<0 || fi>=tac->sampleNr) {ret++; continue;}
    if(csv->c[i].col<0) {ret++; continue;}
    ri=csv->c[i].col-1; 
    if(tac->isframe) ri--; 
    if(ri>=tac->tacNr) {ret++; continue;}
    v=atofVerified(csv->c[i].content); if(!isnan(v)) oknr++;
    if(verbose>10) printf("   -> fi=%d\tri=%d\t=\t%g\n", fi, ri, v);
    if(ri<0) {
      if(tac->isframe) {
        if(ri==-2) tac->x1[fi]=v; else if(ri==-1) tac->x2[fi]=v;
      } else tac->x[fi]=v;
    } else tac->c[ri].y[fi]=v;
  }
  if(verbose>0 && ret>0)
    printf("%d error(s) in reading DFT file format.\n", ret);
  if(oknr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
    return TPCERROR_INVALID_FORMAT;
  }
  /* Set middle sample times when necessary */
  if(tac->isframe)
    for(fi=0; fi<tac->sampleNr; fi++)
      tac->x[fi]=0.5*(tac->x1[fi]+tac->x2[fi]);
 
  /* If units are not known, try to read units from header */
  if(tac->cunit==UNIT_UNKNOWN) tacGetHeaderUnit(tac, NULL);
  else tacSetHeaderUnit(&tac->h, UNIT_UNKNOWN);
  if(tac->tunit==UNIT_UNKNOWN) tacGetHeaderTimeunit(tac, NULL);
  else tacSetHeaderTimeunit(&tac->h, UNIT_UNKNOWN);
 
  /* Move column containing weights to its correct place in the struct */
  ret=tacWMove(tac, 1, status);
  if(ret!=TPCERROR_OK && ret!=TPCERROR_NO_WEIGHTS) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  if(tac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  tac->format=TAC_FORMAT_DFT;
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/