GetSubSim.cpp

#include <cstdio>
#include <cstdlib>
#include <sys/stat.h>
#include <cstring>
#include <vector>
#include <math.h>

#include <string> //for str operations
using namespace std;

// a general matrix style structure
template <typename T>
struct matrix{
  int x;
  int y;
  T** data;
};

template <typename T>
struct array{
  int x;
  T* data;
};

template <typename T>
T *collapse(std::vector<T> &v){
  T *tmp = new T[v.size()];
  for(int i=0;i<v.size();i++)
    tmp[i] = v[i];
  return tmp;
}

//function to cleanup our generic matrix structure
template <typename T>
void cleanup(matrix<T> &m){//using a reference to avoid copying the data
  for(int i=0;i<m.x;i++)
    delete [] m.data[i];
  delete [] m.data;
}

// get the filesize of a file
size_t fsize(const char* fname){
  struct stat st ;
  stat(fname,&st);
  return st.st_size;
}

// find out if a file exists
int fexists(const char* str) {
  struct stat buffer ;
  return (stat(str, &buffer )==0 );
}

// a nice wrapper for getting files
FILE *getFILE(const char*fname,const char* mode) {
  int writeFile = 0;
  for(size_t i=0;i<strlen(mode);i++)
    if(mode[i]=='w')
      writeFile = 1;
  if(writeFile&&fexists(fname)){//DRAGON
    fprintf(stderr,"\t-> File exists: %s exiting...\n",fname);
    exit(0);
  }
  FILE *fp;
  if(NULL==(fp=fopen(fname,mode))){
    fprintf(stderr,"\t->Error opening FILE handle for file:%s exiting\n",fname);
    exit(0);
  }
  return fp;
}

// CHECK HERE FOR USAGE

// this program aims at selecting only a subset of individuals and sites from the file .glf.gz generated by ngsSim
// takes the first nsites_new and nind_new out of the original values

int main (int argc, char *argv[]) {

  char *infile=NULL;
  char *outfile=NULL;

  int argPos = 1, nind = 0, increment = 0, nind_new = 0, check=0, ncat=10;
  long int nsites = 0, nsites_new = 0;

  // CHECK HERE FOR USAGE

  while (argPos<argc) {
    increment = 0;
    if(strcmp(argv[argPos],"-infile")==0)
      infile = argv[argPos+1];
    else if(strcmp(argv[argPos],"-nind")==0)
      nind = atoi(argv[argPos+1]);
    else if(strcmp(argv[argPos],"-nsites")==0)
      nsites = atoi(argv[argPos+1]);
    else if(strcmp(argv[argPos],"-ncat")==0)
      ncat = atoi(argv[argPos+1]);
    else if(strcmp(argv[argPos],"-outfile")==0)
      outfile = argv[argPos+1];
    else if(strcmp(argv[argPos],"-nind_new")==0)
      nind_new = atoi(argv[argPos+1]);
    else if(strcmp(argv[argPos],"-nsites_new")==0)
      nsites_new = atoi(argv[argPos+1]);
    else if(strcmp(argv[argPos],"-check")==0)
      check = atoi(argv[argPos+1]);
    else {
      printf("\tUnknown arguments: %s\n",argv[argPos]);
      return 0;
    }
    argPos = argPos + 2 + increment;
  } 

  FILE *fin = getFILE(infile,"rb");

  FILE *fout = getFILE(outfile, "wb");

  size_t filesize = fsize(infile);
  if((filesize % (sizeof(double)*(nind*nsites*ncat)) )) {
    fprintf(stderr, "\n\t-> Possible error, binaryfiles might be broken. %d %lu\n", nind, nsites);
    exit(-1);
  }

  int conta_ind=1; // 1 based
  int conta_site=1; // 1 based

  long int ntot=0;

  for (int i=0; i<(nsites*nind); i++) {

    double *tmp = new double[ncat];
    fread(tmp,sizeof(double),ncat,fin);
    
    if (conta_site<=nsites_new && conta_ind<=nind_new) {
      fwrite(tmp, sizeof(double), ncat, fout);
      ntot++;
    }
 
    conta_ind++;
    if ((conta_ind-1)==nind) { conta_site++; conta_ind=1; }

    delete [] tmp;

  }

  fclose(fin);
  fclose(fout);

  fprintf(stderr, "\nExpected %lu lines and written %lu lines\n", nsites_new*nind_new, ntot);

  if (check) {
    FILE *fp = getFILE(outfile, "rb");
    size_t filesize =fsize(outfile);
    if((filesize %(sizeof(double)*(nind_new*nsites_new*ncat)) )) {
      fprintf(stderr,"\n\t-> Possible error, binaryfiles might be broken\n");
      fclose(fp);
      exit(-1);
    }
    fclose(fp);
  }

} // main