function_parser.c

Go to the documentation of this file.

/*
 *function_parser
 *parses the wavefunctions from Craig Pryor's dotcode.
n 
 Copyright (C) 2006 Joseph Pingenot

 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 2 of the License, or
 (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

*/

#include <sys/stat.h>
#include <stdlib.h>
#include <fcntl.h>
#include <unistd.h>
#include <errno.h>
#include <string.h>
#include <stdio.h>
#include "function_parser.h"
#include "config.h"
#include "common.h"

#define FUNCTION_PARSER_GERROR_DOMAIN 5031

/*If the "File" is a directory, we will tack this on to the end.*/
#define FILE_TO_APPEND ".psi.dat"

/*For efficiency, we read in some number of points at a time.
 *4 sets of 8 bands
 *NOTE: the buffer size MUST end up greater than or equal to the header size!
 */
#define BUFFERED_POINTS 32
/*NOTE: this should have 8, but may in the future have more or less.
 *  We use this number elsewhere, so it's size of a single band at a single
 *  point.
 */
#define POINT_SIZE (sizeof(double)*2)
#define HEADER_SIZE (12+1+3+4*sizeof(int)+4*sizeof(double))
/*Don't forget we have 2 doubles: real and im parts of wavefunc*/
/*Assume 8 points.*/
#define DATA_BUFFER_SIZE (BUFFERED_POINTS*POINT_SIZE*8)

/*Helper functions: update band, x, y, z if needed.*/
void increment_y(int *y, int *z, struct wavefunc_header *header) {
  (*y)++;
  if((*y) >= header->Ny) {
    (*y)=0;
    (*z)++;
    /*We might want to warn *at* Nz at some point in the future.
     *Right now, we increment before we hit solid EOF (0 bytes returned)
     *  so we will end up getting to z = Nz.  More than that, though, means
     *  the file's broke.
     */
    if((*z) > header->Nz) {
      fprintf(stderr, "WARNING: increment_y: incremented z to %d, but Nz is %d!  This is a corrupted file!\n", *z, header->Nz);
    }
  }else{
  }
}
void increment_x(int *x, int *y, int *z, struct wavefunc_header *header) {
  (*x)++;
  if((*x) >= header->Nx) {
    (*x) = 0;
    increment_y(y, z, header);
  }else{
  }
}
void increment_band(int *band, int *x, int *y, int *z, struct wavefunc_header *header) {
  (*band)++;
  if((*band) >= header->Nc) {
    (*band) = 0;
    increment_x(x, y, z, header);
  }else{
  }
}


/*Documented in function_parser.h*/
int read_wavefunction(char* filename, struct wavefunc_header *header, double **matrix, int *err) {
  fprintf(stderr, "Got filename %s\n", filename);
  GFile* f = g_file_new_for_commandline_arg(filename);
  GError *e = NULL;
  int arr;
  arr = read_wavefunction_gfile(f, header, matrix, &e);
  if(e != NULL) {
    *err = e->code;
    fprintf(stderr, "function_parser::read_wavefunction(char*)<legacy>: Got error %d; message=%s\n", e->code, e->message);
    free(e);
  }
  return arr;
}

int read_wavefunction_gfile(GFile* f, struct wavefunc_header *header, double **matrix, GError **err) {
  GError *e=NULL;
  GFile *psifile = g_file_get_child(f, FILE_TO_APPEND);
  GFileInputStream *infile = g_file_read(psifile, NULL, &e);
  if(e!= NULL) {*err = e; return -3;}
  int count;
  GError *errn;
  char buffer[DATA_BUFFER_SIZE];
  int n_pts;
  int band;
  /*We allocate the matrix locally, so that we don't stomp on anything
   *unnecessarily. Kinda overkill, but whatever.*/
  double *mtrix;
  int pt;
  int ierr;
  char *bptr;
  int x, y, z;
  if((header == NULL) || (f == NULL)) {
    *err = g_error_new(FUNCTION_PARSER_GERROR_DOMAIN, 0, "function_parser::read_wavefunction_gfile: header and/or filename and/or matrix was NULL!");
    if(header == NULL)    return -1;
    if(f == NULL)  return -2;
    if(matrix == NULL)  return -5;
  }
  if((count = resilient_read(infile, buffer, HEADER_SIZE, &errn)) < HEADER_SIZE) {
#ifdef DEBUG3
    fprintf(stderr, "INFO: read_wavefunction: read in only %d bytes, not %lu (when reading header)\n", count, (long unsigned int)HEADER_SIZE);
#endif
    *err = errn;
    return -4;
  }
  if((ierr = parse_header(header, buffer)) < 0) {
    *err = g_error_new(FUNCTION_PARSER_GERROR_DOMAIN, ierr, "function_parser::read_wavefunction_gfile: Failed to parse header: %d", ierr);
    return -7;
  }
  
  /*Allocate the matrix storage*/
  if((mtrix = malloc(2*header->Nc*header->Nx*header->Ny*header->Nz*sizeof(double))) == NULL) {
    *err = g_error_new(FUNCTION_PARSER_GERROR_DOMAIN, 0, "function_parser::read_wavefunction_gfile: Failed to allocate storage for matrix");
    return -8;
  }
  /*FROM HERE ON OUT, WE MUST FREE MATRIX STORAGE IF WE'RE FAILING!*/
  band=x=y=z=0;
  while(((count = resilient_read(infile, buffer, DATA_BUFFER_SIZE, &errn)) != 0) && (count != -1)) {
#ifdef DEBUG3
    fprintf(stderr, "read in %d bytes, err=%d (%s)\n", count, (errn!=NULL)?(errn->code):0, (errn!=NULL)?(errn->message):"");
#endif
    /*We read in buffers!*/
    if((count % (POINT_SIZE*header->Nc)) != 0) {
      fprintf(stderr, "WARNING: read_wavefunction: did not read in an full point for some reason!  Ignoring extras\n");
    }
    n_pts = (count - (count%(POINT_SIZE*header->Nc)))/(POINT_SIZE*header->Nc);
#ifdef DEBUG3
    fprintf(stderr, "buffer holds %d bytes, n_pts=%d (n_pts*(POINT_SIZE=%lu)*(Nc=%d) = %lu)\n", count, n_pts, (long unsigned int)POINT_SIZE, header->Nc, (long unsigned int)n_pts*POINT_SIZE*header->Nc);
#endif
    /*point bptr at the start of the buffer.*/
    bptr = buffer;
    /*8 complex numbers per point (8 bands)*/
    n_pts *= header->Nc;
    for(pt=0; pt<n_pts; ++pt) {
      /*Re*/
      mtrix[band*2 + x*header->Nc*2 + y*header->Nx*header->Nc*2 + z*header->Ny*header->Nx*header->Nc*2] = get_binary_double(header->little_endian, bptr);
      /*update bptr.*/
      bptr += sizeof(double);
      /*Im*/
      mtrix[1+band*2 + x*header->Nc*2 + y*header->Nx*header->Nc*2 + z*header->Ny*header->Nx*header->Nc*2] = get_binary_double(header->little_endian, bptr);
      bptr += sizeof(double);
#ifdef DEBUG3
      fprintf(stderr, "band=%d x=%d y=%d z=%d: Re=%g Im=%g\n", band, x, y, z, mtrix[band*2 + x*header->Nc*2 + y*header->Nx*header->Nc*2 + z*header->Ny*header->Nx*header->Nc*2], mtrix[1+band*2 + x*header->Nc*2 + y*header->Nx*header->Nc*2 + z*header->Ny*header->Nx*header->Nc*2]);
#endif
      /*Update band, x, y, z
       *NOTE that increment_band increments x, y, and z if needed.
       */
      increment_band(&band, &x, &y, &z, header);
    }
  }
#ifdef DEBUG3
  fprintf(stderr, "read in %d bytes, err=%d (%s)\n", count, (errn!=NULL)?(errn->code):0, (errn!=NULL)?(errn->message):"");
#endif
  /*Set matrix to the temporary matrix.*/
  *matrix = mtrix;
  /*Return.  Success!*/
  *err = 0;
  return 0;
}

/*Reads in all of the data, despite interruptions and junk.
 *ARGS:
 * if       GFileInputStream* input file stream
 * buffer   place to store final data read in
 * count    number of chars to read
 * err      pointer to pointer to GError, for permanent errors
 *RETURNS:
 * if successful, number of bytes read and err is set to 0.
 * if failed, -1 and err is set to errno.
 * if partially successful, number of bytes read and err is set to errno.
 *NOTES:
 * * "partially successful" means that some data was read before a fatal error
 *    occurred.
 * * You almost certainly want to check errno if you care that a fatal error
 *    occurred.
 * * Will set *err to NULL if we have success.
 */
int resilient_read(GFileInputStream *is, char *buffer, int count, GError **err) {
  int in;
  /*Start off at the start of the buffer.*/
  char *buff=buffer;
  int to_read=count;
  GError *e = NULL;
  /*Read it right into buffer.  We're just gonna copy it over anyway.*/
  while((to_read > 0) && ((in = g_input_stream_read((GInputStream*)is, buff, to_read, NULL, &e)) < to_read)) {
    if(in > 0) {
      /*We read in data; just update where we are and how much we need to read.*/
      buff += in;
      to_read -= in;
      continue;
    }
    if(e != NULL) {
      /*Permanent error.*/
      *err=e;
      if((count - to_read) > 0) {
  /*We actaully read in data before error cropped up.*/
  return count - to_read;
      }else{
  /*No data, just error.*/
  return -1;
      }
    }
    if(in == 0) {
      /*we hit EOF.*/
      *err = NULL;
      /*If we started at EOF, to_read will be count, so we'll return 0.
       *Else, we will return the number of bytes read before EOF.*/
      return count - to_read;
    }
    e = NULL;
  }
  /*If we made it here, we just read in all of the data in one go.*/
  *err = NULL;
  return count;
}

/*Fills in header struct from a buffer containing header info.
 *ARGS:
 *  header     pointer to header struct to fill in
 *  buff       buffer containing header data.
 *RETURNS:
 *  0          if successful
 * -1          if bad format string.
 * -2          if bad endian value (accepted is 'l' and 'b')
 * -3          if bad format string (i.e. doesn't start with 'v' and/or end
 *                  with '\n')
 * -4          if we don't understand this version.
 */
int parse_header(struct wavefunc_header *header, char *buff) {
  char *b=buff;
  /*First, we fill in the format and endianness info, since
   *   we will need for converting the ints and things.
   */
  memcpy(&(header->format[0]), b, 12);
  header->format[12]='\0';
  b += 12;
  header->endian = b[0];
  b++;
  header->version_string[0]=b[0];
  header->version_string[1]=b[1];
  header->version_string[2]=b[2];
  header->version_string[3]='\0';
  b += 3;
  /*Checks and fill in derive info (e.g. endianness int)*/
  if(strcmp(header->format, "CompWaveFxnb") != 0) return -1;
  if(header->endian == 'b') {
    header->little_endian=0;
  }else if(header->endian == 'l') {
    header->little_endian=1;
  }else{
    return -2;
  }
  if(header->version_string[0] != 'v') return -3;
  if(header->version_string[1] != '2') return -4;
  if(header->version_string[2] != '\n') return -3;
  else header->version = 2;

  /*Now fetch doubles and ints.*/
  header->Nc = get_binary_int(header->little_endian, b);
  b += FILE_SIZE_INT;
  header->Nx = get_binary_int(header->little_endian, b);
  b += FILE_SIZE_INT;
  header->Ny = get_binary_int(header->little_endian, b);
  b += FILE_SIZE_INT;
  header->Nz = get_binary_int(header->little_endian, b);
  b += FILE_SIZE_INT;
  
  header->dx = get_binary_double(header->little_endian, b);
  b += FILE_SIZE_DOUBLE;
  header->dy = get_binary_double(header->little_endian, b);
  b += FILE_SIZE_DOUBLE;
  header->dz = get_binary_double(header->little_endian, b);
  b += FILE_SIZE_DOUBLE;
  header->energy = get_binary_double(header->little_endian, b);
  b += FILE_SIZE_DOUBLE;

  return 0;
}

/*Reads in a binary double from a buffer
 *ARGS:
 *  le     buffer is in little endian format.
 *  buffer buffer containing the data to be read.
 *RETURNS:
 *  the number stored in buffer
 *NOTES:
 *  If you pass too few bytes, you're in trouble!
 */
double get_binary_double(int le, char *buffer) {
  double value;
  char *ptr = (char*)(&value);
  int i, j;
#ifdef DEBUG2
  fprintf(stderr, "INFO: get_binary_double: le=%d, WORDS_BIGENDIAN=%d\n", le,
#ifdef WORDS_BIGENDIAN
    1
#else
    0
#endif
    );
#endif
#ifdef WORDS_BIGENDIAN
  if(!le) {
#else
  if(le) {
#endif
    /*Same endian-ness.  No conversion needed.*/
#ifdef DEBUG2
    fprintf(stderr, "INFO: get_binary_double: no need to convert endianness.\n");
#endif
    memcpy(ptr, buffer, FILE_SIZE_DOUBLE);
  }else{
    /*They're different.  We must bit-flip.*/
#ifdef DEBUG2
    fprintf(stderr, "INFO: get_binary_double: MUST convert endianness.\n");
#endif
    for(i=0, j=(FILE_SIZE_DOUBLE-1); i<FILE_SIZE_DOUBLE; ++i, --j) {
      ptr[i]=buffer[j];
    }
  }
  return value;
}

/*Reads in a binary int from a buffer
 *ARGS:
 *  le     buffer is in little endian format.
 *  buffer buffer containing the data to be read.
 *RETURNS:
 *  the number stored in buffer
 *NOTES:
 *  If you pass too few bytes, you're in trouble!
 */
int get_binary_int(int le, char *buffer) {
  int value;
  char *ptr = (char*)(&value);
  int i, j;
#ifdef DEBUG2
  fprintf(stderr, "INFO: get_binary_int: le=%d, WORDS_BIGENDIAN=%d\n", le, 
    #ifdef WORDS_BIGENDIAN
    1
    #else
    0
    #endif
    );
#endif
#ifdef WORDS_BIGENDIAN
  if(!le) {
#else
  if(le) {
#endif
    /*Same endian-ness.  No conversion needed.*/
#ifdef DEBUG2
    fprintf(stderr, "INFO: get_binary_int: no need to convert endianness.\n");
#endif
    memcpy(ptr, buffer, FILE_SIZE_INT);
  }else{
    /*They're different.  We must bit-flip.*/
#ifdef DEBUG2
    fprintf(stderr, "INFO: get_binary_int: MUST convert endianness.\n");
#endif
    for(i=0, j=(FILE_SIZE_INT-1); i<FILE_SIZE_INT; ++i, --j) {
      ptr[i]=buffer[j];
    }
  }
  return value;
}