/*
  gbget  (ver. 5.6) -- Basic data extraction and manipulation tool
  Copyright (C) 2011-2018 Giulio Bottazzi

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
  (version 2) as published by the Free Software Foundation;
  
  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 Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

#include "tools.h"

#if defined HAVE_LIBGSL
#include "assert.h"
#endif

#if defined HAVE_LIBGSL
#include "matheval.h"
#endif


/* input management */

/*

  parse_spec: parse the input line identifying the different pieces:
  file name, block spec, slice spec and transformation

  parse_slice: used to extract information about which blocks, columns
  or rows to pick from the command line specification

  new_matrix: build the new matrix from the read-in block, according
  to columns and rows specs.

  add_matrix: add the matrix to the output data

 */

void parse_spec(const char* specstring,char ** filename,char **blockslice,
		char **columnslice,char **rowslice,char **transform){

  char *piece=strdup (specstring);
  size_t idx1,idx2;
  char *pointer=piece;

  /* initial checks */
  if( (index(specstring,'[')!=NULL && index(specstring,']')==NULL ) 
      || (index(specstring,']')!=NULL && index(specstring,'[')==NULL ) ){
    fprintf(stderr,"ERROR (%s): error in block specification\n",GB_PROGNAME);
    exit(-1);
  }

  if( (index(specstring,'(')!=NULL && index(specstring,')')==NULL) 
      || (index(specstring,')')!=NULL && index(specstring,'(')==NULL)  ){
    fprintf(stderr,"ERROR (%s): error in slice specification\n",GB_PROGNAME);
    exit(-1);
  }


  /* filename */
  free(*filename);
  idx1=strcspn(pointer,"[(");
  if(idx1==0)
    *filename  = strdup("");
  else if(idx1<strlen(pointer))
    *filename  = my_strndup(pointer,idx1);
  else
    *filename  = strdup(pointer);

/*   fprintf(stderr,"file:-%s-\n",*filename); */

  /* point after the filename */
  pointer=piece+idx1;

  /* blocks */
  if(index(pointer,'[')!=NULL){
    idx1=strcspn(pointer,"[");
    idx2=strcspn(pointer+idx1,"]");
    *blockslice = my_strndup(pointer+idx1+1,idx2-1);
  }
  else
    *blockslice = strdup("");

/*   fprintf(stderr,"block:-%s-\n",*blockslice);  */

  /* define slice specification */
  if(index(pointer,'(')!=NULL){
    char *stmp3;
    idx1=strcspn(pointer,"(");
    idx2=strcspn(pointer+idx1,")");

    stmp3= my_strndup(pointer+idx1+1,idx2-1);

    if(index(stmp3,',')!=NULL){
      idx1=strcspn(stmp3,",");
      *columnslice=my_strndup(stmp3,idx1);
      *rowslice=strdup (stmp3+idx1+1);
    }
    else{
      *columnslice = strdup(stmp3);
      *rowslice=strdup ("");
    }
    free(stmp3);
  }
  else{
      *columnslice=strdup ("");
      *rowslice=strdup ("");
  }

/*   fprintf(stderr,"col:-%s-\n",*columnslice);  */
/*   fprintf(stderr,"row:-%s-\n",*rowslice);  */


  /* define trasformations */
  idx1=0; 
  if(index(pointer,'(')!=NULL) 
    idx1=strcspn(pointer,")");    
  
  idx2=0; 
  if(index(pointer,'[')!=NULL) 
    idx2=strcspn(pointer,"]"); 
  
  if(idx1==0 && idx2==0) 
    *transform = strdup(""); 
  else{ 
    if(idx2>idx1) 
      idx1=idx2; 
    *transform = strdup(pointer+idx1+1); 
  }
  
/*   fprintf(stderr,"pointer:%s len:%d idx1:%d idx2:%d trasf:%s\n", */
/* 	  pointer,strlen(pointer),idx1,idx2,*transform);  */


  free(piece);
/*   exit(0); */

  /* final checks */
  if(strspn(*blockslice,"-1234567890:;")<strlen(*blockslice)){
    fprintf(stderr,"ERROR (%s): wrong blocks specification: %s\n",GB_PROGNAME,*blockslice);
    exit(-1);
  }
  if(strspn(*columnslice,"-1234567890:;")<strlen(*columnslice)){
    fprintf(stderr,"ERROR (%s): wrong columns specification: %s\n",GB_PROGNAME,*columnslice);
    exit(-1);
  }
  if(strspn(*rowslice,"-1234567890:;")<strlen(*rowslice)){
    fprintf(stderr,"ERROR (%s): wrong rows specification: %s\n",GB_PROGNAME,*rowslice);
    exit(-1);
  }

}

/* split the spec according to ';' and build arrays of initial, final
   and skip entries */
void parse_slice(const char * slice,const char * typename,const int max,
		 size_t * piecenum, size_t **indexini, size_t **indexfin,int **indexskip){
	
  char *pieces=strdup (slice);
  char *stmp1=pieces;
  char *stmp2;

  *piecenum=0;

  /* split the ';' pieces */
  while( (stmp2=strsep (&stmp1,";")) != NULL)
    {

      /* default values */
      int ini=1;
      int fin=max;
      int skip=1;

      /* increase the number of pieces */
      (*piecenum)++;

      if(strcmp(stmp2,"")!=0){

	char *stmp3 = strdup(stmp2);
	char *stmp4 = stmp3;
	char *stmp5;

	/* read ini and fin */
	if( (stmp5=strsep (&stmp4,":")) != NULL && strlen(stmp5)>0 )
	  fin = ini = atoi(stmp5);
	if( (stmp5=strsep (&stmp4,":")) != NULL){
	  if(strlen(stmp5)>0) fin = atoi(stmp5); 
	  else fin=max;
	}

	/* treat negative elements */
	if(ini<0) ini=max+ini+1;
	if(fin<0) fin=max+fin+1;

	/* eventually set skip */
	if( (stmp5=strsep (&stmp4,",")) != NULL && strlen(stmp5)>0 )
	  skip= atoi(stmp5);
	else if(ini>fin) skip=-1;

	free(stmp3);

	/*check the boundaries and skip*/
	if(ini<1){
	  fprintf(stderr,
		  "WARNING (%s): requested %s %d out of range: set to %d\n",GB_PROGNAME,typename,ini,1);
	  ini=1;
	}
	else if(ini>max){
	  fprintf(stderr,"WARNING (%s): requested %s %d out of range: set to %d\n",GB_PROGNAME,typename,ini,max);
	  ini=max;
	}
	if(fin<1){
	  fprintf(stderr,"WARNING (%s): requested %s %d out of range: set to %d\n",GB_PROGNAME,typename,fin,1);
	  fin=1;
	}
	else if(fin>max){
	  fprintf(stderr,"WARNING (%s): requested %s %d out of range: set to %d\n",GB_PROGNAME,typename,fin,max);
	  fin=max;
	}
	if( ( (ini > fin) && (skip>0) ) || ((ini < fin) && (skip<0)) ){
	  fprintf(stderr,"WARNING (%s): requested %s out of range: increment reversed\n",
		  GB_PROGNAME,typename);
	  skip=-skip;
	}
      }

      /* assign values */
      *indexini  = (size_t *) my_realloc((void *) *indexini,(*piecenum)*sizeof(size_t));
      *indexfin  = (size_t *) my_realloc((void *) *indexfin,(*piecenum)*sizeof(size_t));
      *indexskip = (int *) my_realloc((void *) *indexskip,(*piecenum)*sizeof(int));
      
      (*indexini)[(*piecenum)-1]=ini;
      (*indexfin)[(*piecenum)-1]=fin;
      (*indexskip)[(*piecenum)-1]=skip;
      
    }
  
  
  free(pieces);

}


void new_matrix(double ***newdata,size_t *columnlength,size_t *rowlength,
		 double ** const vals,const size_t columns,const size_t rows,
		const char *columnslice,const char *rowslice){


   size_t i;
   
   size_t cslicenum=0;
   size_t *cini=NULL;
   size_t *cfin=NULL;
   int *cskip=NULL;

   size_t rslicenum=0;
   size_t *rini=NULL;
   size_t *rfin=NULL;
   int *rskip=NULL;

   /* read the specs */
   parse_slice(columnslice,"column",columns,&cslicenum,&cini,&cfin,&cskip);
   parse_slice(rowslice,"row",rows,&rslicenum,&rini,&rfin,&rskip);

   /* decide new matrix dimensions */
   *columnlength=0;
   for(i=0;i<cslicenum;i++)
     *columnlength += 
       (cfin[i]>cini[i] ? cfin[i]-cini[i] : cini[i]-cfin[i])/(cskip[i]>0?cskip[i]:-cskip[i])+1;

   *rowlength=0;
   for(i=0;i<rslicenum;i++)
     *rowlength += 
       (rfin[i]>rini[i] ? rfin[i]-rini[i] : rini[i]-rfin[i])/(rskip[i]>0?rskip[i]:-rskip[i])+1;

   
/*     for(i=0;i<cslicenum;i++)  */
/*       fprintf(stderr,"col[%zd] %zd:%zd:%d\n",i,cini[i],cfin[i],cskip[i]);  */
/*     for(i=0;i<rslicenum;i++)  */
/*       fprintf(stderr,"row[%zd] %zd:%zd:%d\n",i,rini[i],rfin[i],rskip[i]);  */
/*     fprintf(stderr,"clen=%zd rlen=%zd\n",*columnlength,*rowlength);  */

   
   /* build the new matrix */
   {
     size_t cs,rs;
     int c,r;
     size_t column=0,row=0;

     /* allocate the new matrix */
     /*    if(*newdata != NULL) free(*newdata); */
     double **data = (double **) my_alloc((*columnlength)*sizeof(double *));
     for(i=0;i<*columnlength;i++)
       data[i] = (double *) my_alloc((*rowlength)*sizeof(double));

     column=0;
     for(cs=0;cs<cslicenum;cs++){
       for(c=cini[cs]-1; ( cskip[cs]>0 ? c < cfin[cs] : c+2 > cfin[cs] ) ;c+=cskip[cs]){
	 row=0;
	 for(rs=0;rs<rslicenum;rs++){
	   for(r=rini[rs]-1; ( rskip[rs]>0 ? r < rfin[rs] : r+2 > rfin[rs] ) ;r+=rskip[rs]){
/*   	     fprintf(stderr,"val[%zd][%zd]=%f -> ",   */
/*   		     c,r,vals[c][r]);   */
	     data[column][row]=vals[c][r];
/*   	     fprintf(stderr,"data[%zd][%zd]=%f\n",   */
/*   		     column,row,data[column][row]);   */
	     row++;
	   }
	 }
	 column++;
       }
     }

     *newdata=data;
   }


   free(cini);
   free(cfin);
   free(cskip);

   free(rini);
   free(rfin);
   free(rskip);

}



void data_transpose(double *** vals,size_t *rows,size_t *columns){

  size_t i,j;

  /* allocate space */
  double **newvals = (double **) my_alloc((*rows)*sizeof(double *));
  for(i=0;i<*rows;i++){
    newvals[i] = (double *) my_alloc((*columns)*sizeof(double));
  }

  /* fill space */
  for(i=0;i<*rows;i++)
    for(j=0;j<*columns;j++)
      newvals[i][j] = (*vals)[j][i];

  /* deallocate old space */
  if(*columns>0)
    for(i=0;i<*columns;i++)
      free((*vals)[i]);
  free(*vals);
  *vals = newvals;

  /* swap axes */
  j=*columns;
  *columns=*rows;
  *rows=j;

}


void data_flat(double *** vals,size_t *rows,size_t *columns){

  size_t i,j,index;

  /* allocate space */
  double **newvals = (double **) my_alloc(sizeof(double *));
  newvals[0] = (double *) my_alloc((*columns)*(*rows)*sizeof(double));

  /* fill space */
  index=0;
  for(i=0;i<*columns;i++)
    for(j=0;j<*rows;j++)
      newvals[0][index++] = (*vals)[i][j];

  /* deallocate old space */
  for(i=0;i<*columns;i++)
    free((*vals)[i]);
  free(*vals);
  *vals = newvals;

  /* new axes */
  *rows=(*columns)*(*rows);
  *columns=1;
}

void data_diff(double *** vals,size_t *rows,size_t *columns){

  size_t i,j;

  /* allocate space */
  double **newvals = (double **) my_alloc((*columns-1)*sizeof(double *));
  for(i=0;i<*columns-1;i++)
    newvals[i] = (double *) my_alloc((*rows)*sizeof(double));

  /* fill space */
  for(i=0;i<*columns-1;i++)
    for(j=0;j<*rows;j++)
      newvals[i][j] = (*vals)[i+1][j]-(*vals)[i][j];

  /* deallocate old space */
  if(*columns>0)
    for(i=0;i<*columns;i++)
      free((*vals)[i]);
  free(*vals);
  *vals = newvals;

  /* new dimensions */
  *columns=*columns-1;

}

#if defined HAVE_LIBGSL

/* handy structure used in computation */
typedef struct function {
  void *f;
  size_t argnum;
/*   char *args[] = { "x" }; */
  char **args;			/* list of arguments names */
  size_t *indeces;	        /* the index (shift) for each variable */
  size_t *lags;	                /* the lag for each variable */
  double *values;		/* values used in computation */
  size_t maxindex;		/* the largest index */
  size_t maxlag;		/* the largest lag */
} Function;

void data_applyfuns(double *** myvals,size_t *myrows,size_t *mycolumns,const char *funslist){

  Function *F=NULL;
  size_t Fnum=0;

  char *funspecs=strdup(funslist);
  char *specstart=funspecs,*funspec;

  /* traslate to more handy names */
  double **data=*myvals;
  size_t columns=*mycolumns;
  size_t rows=*myrows;

/*   fprintf(stderr,"func spec=%s\n",funspecs); */


  /* parse line for functions specification */
  while( (funspec=strsep (&funspecs,";")) != NULL ){
    char *piece=strdup (funspec);
    char *stmp1=piece;
    char *stmp2;
    
    while( (stmp2=strsep (&stmp1,",")) != NULL ){/*take a piece*/
      size_t j;
      char *stmp3 = strdup(stmp2);
      /* 	fprintf(stderr,"token:%s\n",stmp3); */
      
      /* allocate new function */
      Fnum++;
      F=(Function *) my_realloc((void *) F,Fnum*sizeof(Function));
      
      /* generate the formal expression */
      F[Fnum-1].f = evaluator_create (stmp3);
      assert(F[Fnum-1].f);
      free(stmp3);
      
      /* retrive list of argument */
      {
	/* hack necessary due to libmatheval propotypes */
	int argnum;
	evaluator_get_variables (F[Fnum-1].f, &(F[Fnum-1].args), &argnum);
	F[Fnum-1].argnum = (size_t) argnum;
      }
      /* compute list of indeces */
      F[Fnum-1].indeces = (size_t *) my_alloc(F[Fnum-1].argnum*sizeof(size_t));
      F[Fnum-1].lags = (size_t *) my_alloc(F[Fnum-1].argnum*sizeof(size_t));
      F[Fnum-1].maxindex=0;
      F[Fnum-1].maxlag=0;
      for(j=0;j<F[Fnum-1].argnum;j++){
	stmp3 = (F[Fnum-1].args)[j];

	if(*stmp3 != 'x'){
	  fprintf(stderr,"ERROR (%s): variable name %s should begin with x\n",
		  GB_PROGNAME,stmp3);
	  exit(-1);
	}

	if( *(stmp3+1) == '\0' ){
	  (F[Fnum-1].indeces)[j] = 1;
	  (F[Fnum-1].lags)[j] = 0;

	}
	else {
	  char *endptr;

	  if( strchr(stmp3,'l') == NULL ){
	  
	    (F[Fnum-1].indeces)[j] = strtol(stmp3+1,&endptr,10);
	    (F[Fnum-1].lags)[j] = 0;
	    
	    if(*endptr != '\0'){
	      fprintf(stderr,"ERROR (%s): variable name %s must be of the form x#l#\n",
		    GB_PROGNAME,stmp3);
	      exit(-1);
	    }

	  }
	  else{
	    (F[Fnum-1].indeces)[j] = strtol(stmp3+1,&endptr,10);
	    (F[Fnum-1].lags)[j] = strtol(endptr+1,&endptr,10);

	    if(*endptr != '\0'){
	      fprintf(stderr,"ERROR (%s): variable name %s must be of the form x#l#\n",
		    GB_PROGNAME,stmp3);
	      exit(-1);
	    }

	  }
	  
	}

	/* possibly update the value of the maximum index */
	if( (F[Fnum-1].indeces)[j] > F[Fnum-1].maxindex)
	  F[Fnum-1].maxindex = (F[Fnum-1].indeces)[j];

	/* possibly update the value of the maximum lag */
	if( (F[Fnum-1].lags)[j] > F[Fnum-1].maxlag)
	  F[Fnum-1].maxlag = (F[Fnum-1].lags)[j];

	/* check index  */
	if((F[Fnum-1].indeces)[j]> columns){
	  fprintf(stderr,"ERROR (%s): column %zd not present in data\n",
		  GB_PROGNAME,(F[Fnum-1].indeces)[j]);
	  exit(-1);
	}

	/* check lag */
	if((F[Fnum-1].lags)[j]> rows){
	  fprintf(stderr,"ERROR (%s): lag %zd is longer then the data length\n",
		  GB_PROGNAME,(F[Fnum-1].lags)[j]);
	  exit(-1);
	}

      }

      /* allocate list of values */
      F[Fnum-1].values = (double *) my_alloc(F[Fnum-1].argnum*sizeof(double));

    }
    free(piece);
  }
  
  free(specstart);

  /* generate new matrix */
  {

    /* placeholder for new data */
    double **newvals;

    size_t i,j,h;
    /*       double **output = (double **) my_alloc(Fnum*sizeof(double *));       */
    /*       for(i=0;i<Fnum;i++) output[i] = (double *) my_alloc(rows*sizeof(double)); */
    
    size_t firstrow=0;

    /* compute the first row to print */
    for(j=0;j<Fnum;j++)
      if(F[j].maxlag>firstrow) firstrow=F[j].maxlag;

    /* allocate space */
    newvals = (double **) my_alloc(Fnum*sizeof(double *));
    for(i=0;i<Fnum;i++)
      newvals[i] = (double *) my_alloc((rows-firstrow)*sizeof(double));

    /* fill with new values */
    for(i=firstrow;i<rows;i++){
      for(j=0;j<Fnum;j++){
	for(h=0;h<F[j].argnum;h++)
	  (F[j].values)[h] = (F[j].indeces[h]>0 ?  data[F[j].indeces[h]-1][i-F[j].lags[h]] : i+1-F[j].lags[h]);
	newvals[j][i-firstrow]=evaluator_evaluate (F[j].f,F[j].argnum,F[j].args,F[j].values);
      }
    }

    /* free previous allocated space */
    for(i=0;i<*mycolumns;i++)
      free( (*myvals)[i]);
    free((*myvals));

    /* return newly allocated matrix */
    *mycolumns=Fnum;
    *myrows=rows-firstrow;
    *myvals = newvals;
  }
  

}

void data_applyfuns_recursive(double *** myvals,size_t *myrows,size_t *mycolumns,const char *funslist){

  Function *F=NULL;
  size_t Fnum=0;

  char *funspecs=strdup(funslist);
  char *specstart=funspecs,*funspec;

  /* traslate to more handy names */
  double **data=*myvals;
  size_t columns=*mycolumns;
  size_t rows=*myrows;

/*   fprintf(stderr,"func spec=%s\n",funspecs); */


  /* parse line for functions specification */
  while( (funspec=strsep (&funspecs,";")) != NULL ){
    char *piece=strdup (funspec);
    char *stmp1=piece;
    char *stmp2;
    
    while( (stmp2=strsep (&stmp1,",")) != NULL ){/*take a piece*/
      size_t j;
      char *stmp3 = strdup(stmp2);
      /* 	fprintf(stderr,"token:%s\n",stmp3); */
      
      /* allocate new function */
      Fnum++;
      F=(Function *) my_realloc((void *) F,Fnum*sizeof(Function));
      
      /* generate the formal expression */
      F[Fnum-1].f = evaluator_create (stmp3);
      assert(F[Fnum-1].f);
      free(stmp3);
      
      /* retrive list of argument */
      {
	/* hack necessary due to libmatheval propotypes */
	int argnum;
	evaluator_get_variables (F[Fnum-1].f, &(F[Fnum-1].args), &argnum);
	F[Fnum-1].argnum = (size_t) argnum;
      }
      /* compute list of indeces */
      F[Fnum-1].indeces = (size_t *) my_alloc(F[Fnum-1].argnum*sizeof(size_t));
      F[Fnum-1].lags = (size_t *) my_alloc(F[Fnum-1].argnum*sizeof(size_t));
      F[Fnum-1].maxindex=0;
      F[Fnum-1].maxlag=0;
      for(j=0;j<F[Fnum-1].argnum;j++){
	stmp3 = (F[Fnum-1].args)[j];

	if(*stmp3 != 'x'){
	  fprintf(stderr,"ERROR (%s): variable name %s should begin with x\n",
		  GB_PROGNAME,stmp3);
	  exit(-1);
	}

	if( *(stmp3+1) == '\0' ){
	  (F[Fnum-1].indeces)[j] = 1;
	  (F[Fnum-1].lags)[j] = 0;

	}
	else {
	  char *endptr;

	  if( strchr(stmp3,'l') == NULL ){
	  
	    (F[Fnum-1].indeces)[j] = strtol(stmp3+1,&endptr,10);
	    (F[Fnum-1].lags)[j] = 0;
	    
	    if(*endptr != '\0'){
	      fprintf(stderr,"ERROR (%s): variable name %s must be of the form x#l#\n",
		    GB_PROGNAME,stmp3);
	      exit(-1);
	    }

	  }
	  else{
	    (F[Fnum-1].indeces)[j] = strtol(stmp3+1,&endptr,10);
	    (F[Fnum-1].lags)[j] = strtol(endptr+1,&endptr,10);

	    if(*endptr != '\0'){
	      fprintf(stderr,"ERROR (%s): variable name %s must be of the form x#l#\n",
		    GB_PROGNAME,stmp3);
	      exit(-1);
	    }

	  }

	}

	/* possibly update the value of the maximum index */
	if( (F[Fnum-1].indeces)[j] > F[Fnum-1].maxindex)
	  F[Fnum-1].maxindex = (F[Fnum-1].indeces)[j];

	/* possibly update the value of the maximum lag */
	if( (F[Fnum-1].lags)[j] > F[Fnum-1].maxlag)
	  F[Fnum-1].maxlag = (F[Fnum-1].lags)[j];

	/* check index  */
	if((F[Fnum-1].indeces)[j]> columns){
	  fprintf(stderr,"ERROR (%s): column %zd not present in data\n",
		  GB_PROGNAME,(F[Fnum-1].indeces)[j]);
	  exit(-1);
	}

	/* check lag */
	if((F[Fnum-1].lags)[j]> rows){
	  fprintf(stderr,"ERROR (%s): lag %zd is longer then the data length\n",
		  GB_PROGNAME,(F[Fnum-1].lags)[j]);
	  exit(-1);
	}

      }

      /* allocate list of values */
      F[Fnum-1].values = (double *) my_alloc(F[Fnum-1].argnum*sizeof(double));

    }
    free(piece);
  }
  
  free(specstart);

  /* generate new matrix */
  {

    /* placeholder for new data */
    double **newvals;

    size_t i,j,l,h;
    /*       double **output = (double **) my_alloc(Fnum*sizeof(double *));       */
    /*       for(i=0;i<Fnum;i++) output[i] = (double *) my_alloc(rows*sizeof(double)); */
    
    size_t firstrow=0;
    size_t totalcols=0;

    /* compute the first row to print */
    for(j=0;j<Fnum;j++)
      if(F[j].maxlag>firstrow) firstrow=F[j].maxlag;

    /* compute the total number of columns */
    for(j=0;j<Fnum;j++)
      totalcols += columns-F[j].maxindex+1;

    /* allocate space */
    newvals = (double **) my_alloc(totalcols*sizeof(double *));
    for(i=0;i<totalcols;i++)
      newvals[i] = (double *) my_alloc((rows-firstrow)*sizeof(double));

    /* fill with new values */
    for(i=firstrow;i<rows;i++){
      size_t col=0;
      for(j=0;j<Fnum;j++){
	for(l=0;l<columns-F[j].maxindex+1;l++){
	  for(h=0;h<F[j].argnum;h++)
	    (F[j].values)[h] = (F[j].indeces[h]>0 ?  data[l+F[j].indeces[h]-1][i-F[j].lags[h]] : i+1-F[j].lags[h]);
	  newvals[col+l][i-firstrow]=evaluator_evaluate (F[j].f,F[j].argnum,F[j].args,F[j].values);
	}
	col+=columns-F[j].maxindex+1;
      }
    }

    /* free previous allocated space */
    for(i=0;i<*mycolumns;i++)
      free( (*myvals)[i]);
    free((*myvals));

    /* allocate new space */
    *mycolumns=totalcols;
    *myrows=rows-firstrow;
    *myvals=newvals;

  }
  

}



void data_applyselection(double *** myvals,size_t *myrows,size_t *mycolumns,const char *funslist){

  Function *F=NULL;
  size_t Fnum=0;

  char *funspecs=strdup(funslist);
  char *specstart=funspecs,*funspec;

  /* traslate to more handy names */
  double **data=*myvals;
  size_t columns=*mycolumns;
  size_t rows=*myrows;

  /* fprintf(stderr,"func spec=%s\n",funspecs);  */

  /* parse line for functions specification */
  while( (funspec=strsep (&funspecs,";")) != NULL ){
    char *piece=strdup (funspec);
    char *stmp1=piece;
    char *stmp2;
    
    while( (stmp2=strsep (&stmp1,",")) != NULL ){/*take a piece*/
      size_t j;
      char *stmp3 = strdup(stmp2);
      /* 	fprintf(stderr,"token:%s\n",stmp3); */
      
      /* allocate new function */
      Fnum++;
      F=(Function *) my_realloc((void *) F,Fnum*sizeof(Function));
      
      /* generate the formal expression */
      F[Fnum-1].f = evaluator_create (stmp3);
      assert(F[Fnum-1].f);
      free(stmp3);
      
      /* retrive list of argument */
      {
	/* hack necessary due to libmatheval propotypes */
	int argnum;
	evaluator_get_variables (F[Fnum-1].f, &(F[Fnum-1].args), &argnum);
	F[Fnum-1].argnum = (size_t) argnum;
      }
      /* compute list of indeces */
      F[Fnum-1].indeces = (size_t *) my_alloc(F[Fnum-1].argnum*sizeof(size_t));
      F[Fnum-1].lags = (size_t *) my_alloc(F[Fnum-1].argnum*sizeof(size_t));
      F[Fnum-1].maxindex=0;
      F[Fnum-1].maxlag=0;
      for(j=0;j<F[Fnum-1].argnum;j++){
	stmp3 = (F[Fnum-1].args)[j];

	if(*stmp3 != 'x'){
	  fprintf(stderr,"ERROR (%s): variable name %s should begin with x\n",
		  GB_PROGNAME,stmp3);
	  exit(-1);
	}

	if( *(stmp3+1) == '\0' ){
	  (F[Fnum-1].indeces)[j] = 1;
	  (F[Fnum-1].lags)[j] = 0;

	}
	else {
	  char *endptr;

	  if( strchr(stmp3,'l') == NULL ){
	  
	    (F[Fnum-1].indeces)[j] = strtol(stmp3+1,&endptr,10);
	    (F[Fnum-1].lags)[j] = 0;
	    
	    if(*endptr != '\0'){
	      fprintf(stderr,"ERROR (%s): variable name %s must be of the form x#l#\n",
		    GB_PROGNAME,stmp3);
	      exit(-1);
	    }

	  }
	  else{
	    (F[Fnum-1].indeces)[j] = strtol(stmp3+1,&endptr,10);
	    (F[Fnum-1].lags)[j] = strtol(endptr+1,&endptr,10);

	    if(*endptr != '\0'){
	      fprintf(stderr,"ERROR (%s): variable name %s must be of the form x#l#\n",
		    GB_PROGNAME,stmp3);
	      exit(-1);
	    }

	  }

	}

	/* possibly update the value of the maximum index */
	if( (F[Fnum-1].indeces)[j] > F[Fnum-1].maxindex)
	  F[Fnum-1].maxindex = (F[Fnum-1].indeces)[j];

	/* possibly update the value of the maximum lag */
	if( (F[Fnum-1].lags)[j] > F[Fnum-1].maxlag)
	  F[Fnum-1].maxlag = (F[Fnum-1].lags)[j];

	/* check index  */
	if((F[Fnum-1].indeces)[j]> columns){
	  fprintf(stderr,"ERROR (%s): column %zd not present in data\n",
		  GB_PROGNAME,(F[Fnum-1].indeces)[j]);
	  exit(-1);
	}

	/* check lag */
	if((F[Fnum-1].lags)[j]> rows){
	  fprintf(stderr,"ERROR (%s): lag %zd is longer then the data length\n",
		  GB_PROGNAME,(F[Fnum-1].lags)[j]);
	  exit(-1);
	}

      }

      /* allocate list of values */
      F[Fnum-1].values = (double *) my_alloc(F[Fnum-1].argnum*sizeof(double));

    }
    free(piece);
  }
  
  free(specstart);

  /* generate new matrix */
  {

    /* placeholder for new data */
    double **newvals;

    size_t newrows;

    size_t i,j,h;
    /*       double **output = (double **) my_alloc(Fnum*sizeof(double *));       */
    /*       for(i=0;i<Fnum;i++) output[i] = (double *) my_alloc(rows*sizeof(double)); */
    
    size_t firstrow=0;

    /* compute the first row to print */
    for(j=0;j<Fnum;j++)
      if(F[j].maxlag>firstrow) firstrow=F[j].maxlag;

    /* allocate space */
    newvals = (double **) my_alloc(columns*sizeof(double *));
    for(i=0;i<columns;i++)
      newvals[i] = (double *) my_alloc((rows-firstrow)*sizeof(double));
    
    /* fill with new values */
    newrows=0;
    for(i=firstrow;i<rows;i++){
      char valid=0;
      for(j=0;j<Fnum;j++){
	for(h=0;h<F[j].argnum;h++)	  
	  (F[j].values)[h] = (F[j].indeces[h]>0 ?  data[F[j].indeces[h]-1][i-F[j].lags[h]] : i+1);
	if(evaluator_evaluate (F[j].f,F[j].argnum,F[j].args,F[j].values) >= 0)
	  valid=1;
      }
      /* if passed all checks, add it */
      if(valid==1){
	for(j=0;j<columns;j++)
	  newvals[j][newrows]=data[j][i];
	newrows++;
      }
    }
    /* fprintf(stderr,"newrows %zd\n",newrows); */

    /* free previous allocated space */
    for(i=0;i<*mycolumns;i++)
      free( (*myvals)[i]);
    free((*myvals));

    /* return newly allocated matrix */
    *mycolumns=columns;
    *myrows=newrows;
    *myvals = newvals;

  }
  

}


#endif

 void matrix_transform(double ***matrix,size_t *columnlength,size_t *rowlength,
		       const char *transform,const int o_verbose){

   size_t itmp1;

   for(itmp1=0;itmp1<strlen(transform);itmp1++){
     switch(transform[itmp1]){
     case 't':
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"   %zdx%zd\t[->tran]\t",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       data_transpose(matrix,rowlength,columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"%zdx%zd\n",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       break;
     case 'f':
	/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
	  if(o_verbose==1)
	    fprintf(stderr,"   %zdx%zd\t[->flat]\t",*rowlength,*columnlength);
	  /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
	  data_flat(matrix,rowlength,columnlength);
	  /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
	  if(o_verbose==1)
	    fprintf(stderr,"%zdx%zd\n",*rowlength,*columnlength);
	  /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
	  break;
     case 'l':
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"   %zdx%zd\t[-> log]\t",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       {
	 size_t i,j;
	 for(i=0;i<*columnlength;i++)
	   for(j=0;j<*rowlength;j++)
	     (*matrix)[i][j]=log((*matrix)[i][j]);
       }
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"%zdx%zd\n",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       break;
     case 'd':
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"   %zdx%zd\t[->diff]\t",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       data_diff(matrix,rowlength,columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"%zdx%zd\n",*rowlength,*columnlength);
	  /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       break;
     case 'D':
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"   %zdx%zd\t[->denan]\t",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       denan_data(matrix,rowlength,columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"%zdx%zd\n",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       break;
     case 'w':
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"   %zdx%zd\t[->norm]\t",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       {
	 size_t i,j;
	 for(i=0;i<*columnlength;i++){
	   double sum=0;
	   for(j=0;j<*rowlength;j++)
	     if(finite( (*matrix)[i][j] ))
	       sum += (*matrix)[i][j];
	   for(j=0;j<*rowlength;j++)
	     (*matrix)[i][j]/=sum;
	 }
       }
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"%zdx%zd\n",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       break;
     case 'z':
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"   %zdx%zd\t[->detr]\t",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       {
	 size_t i,j;
	 for(i=0;i<*columnlength;i++){
	   double mean=0;
	   size_t num=0;
	   for(j=0;j<*rowlength;j++)
	     if(finite( (*matrix)[i][j] )){
	       mean += (*matrix)[i][j];
	       num++;
	     }
	   mean /= num;
	   for(j=0;j<*rowlength;j++)
	     (*matrix)[i][j]-=mean;
	 }
       }
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"%zdx%zd\n",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       break;
     case 'Z':
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"   %zdx%zd\t[->zscr]\t",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       {
	 size_t i,j;
	 for(i=0;i<*columnlength;i++){
	   double mean=0,stdev=0,error=0,dtmp1;
	   size_t num=0;
	   for(j=0;j<*rowlength;j++)
	     if(finite( (*matrix)[i][j] ))
	       {
		 mean  += (*matrix)[i][j];
		 num++;
	       }
	   mean /= num;
	   for(j=0;j<*rowlength;j++)
	     if(finite( (*matrix)[i][j] ))
	       {
		 error += (	dtmp1= (*matrix)[i][j] - mean);
		 stdev += dtmp1*dtmp1;
	       }
	   stdev = sqrt( (stdev-error*error/num)/(num-1));
	   for(j=0;j<*rowlength;j++)
	     (*matrix)[i][j]=((*matrix)[i][j]-mean)/stdev;
	 }
       }
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"%zdx%zd\n",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       break;
     case 'P':
       /* ignored; the data will be printed at the end of this slice */

       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"   %zdx%zd\t[print ]\t",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
       if(o_verbose==1)
	 fprintf(stderr,"%zdx%zd\n",*rowlength,*columnlength);
       /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

       break;
#if defined HAVE_LIBGSL
     case '<':
       {
	 size_t i;
	 char *funslist;
	 /* select the list of functions specification */
	 for (i=itmp1;i<strlen(transform) && transform[i] != '>';i++);
	 funslist=my_strndup(transform+itmp1+1,i-itmp1-1);

	 /* apply functions */
	 if(funslist[0]=='@'){	/* recursive application */
	   free(funslist);
	   funslist=my_strndup(transform+itmp1+2,i-itmp1-2);
	   data_applyfuns_recursive(matrix,rowlength,columnlength,funslist);
	 }
	 else 			/* non-recursive application */
	   data_applyfuns(matrix,rowlength,columnlength,funslist);

	 /* update the index */
	 itmp1=i;
       }
       break;
     case '{':
       {
	 size_t i;
	 char *funslist;
	 /* select the list of functions specification */
	 for (i=itmp1;i<strlen(transform) && transform[i] != '}';i++);
	 funslist=my_strndup(transform+itmp1+1,i-itmp1-1);

	 /* apply selection */
	 data_applyselection(matrix,rowlength,columnlength,funslist);

	 /* update the index */
	 itmp1=i;
       }
       break;

#endif
     default:
       fprintf(stderr,
	       "WARNING (%s): unknown transformation \"%c\"; request ignored\n",
	       GB_PROGNAME,transform[itmp1]);
     }
   }
 }

void add_matrix(double ***data, size_t **length, size_t *datanum,
		double **newdata,size_t columnlength,size_t rowlength){
  
  
  (*data)= (double **) my_realloc((void **) (*data),(*datanum+columnlength)*sizeof(double *));
  (*length)= (size_t *) my_realloc((void *) (*length),(*datanum+columnlength)*sizeof(size_t));
  
  {
    size_t i;
    for(i=0;i<columnlength;i++){
      (*data)[*datanum]=newdata[i];
      (*length)[*datanum]=rowlength;
      (*datanum)++;
    }
  }
}


/* check if the data set is a matrix */
int isamatrix(size_t datanum,size_t *length){
  
  size_t i,min,max;
  
  if(datanum==0) 
    return 1;
  
  min=max=length[0];
  for(i=1;i<datanum;i++){
    if(length[i]>max) max = length[i];
    if(length[i]<min) min = length[i];
  }
  
  if(max==min)
    return 1; 			/* true */
  else
    return 0;			/* false */

}


void print_data(size_t *mydatanum, size_t **mylength, double ***mydata, 
		const char* finaltransform, const int o_verbose){

  size_t i;
  size_t datanum=*mydatanum;
  double **data=*mydata;
  size_t *length=*mylength;

  if(datanum==0) return ;

  /*OUTPUT*/
  if(isamatrix(datanum,length)==1){/* final output has equally long columns */
    size_t rowlength = length[0];
    
    /* apply final transformations */
    if(strlen(finaltransform)>0)
      matrix_transform(&data,&datanum,&rowlength,
		       finaltransform,o_verbose);

    /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
    if(o_verbose==1)
      fprintf(stderr,"output => %zdx%zd\n",rowlength,datanum);
    /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

    /* print the result */
    PRINTMATRIXBYCOLS(stdout,data,rowlength,datanum);


  }
  else{/* final output has columns of different lengths */

   /* can't apply final transformations */
   if(strlen(finaltransform)>0)
     fprintf(stderr,
	     "WARNING (%s): cannot transform final output; request ignored\n",GB_PROGNAME);

    /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
   if(o_verbose==1){
     fprintf(stderr,"output (%zd columns of size) =>",datanum);
     for(i=0;i<datanum;i++)
       fprintf(stderr," %zd;",length[i]);
     fprintf(stderr,"\n");
   }
    /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

   PRINTNOMATRIXBYCOLS(stdout,data,length,datanum);
    
  }
  
  /* free printed data */
  free(*mylength);
  *mylength=NULL;
  if(datanum>0) for(i=0;i<datanum;i++) free(data[i]);
  *mydatanum = 0;
  free(data);
  *mydata=NULL;

}

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

  /* actual input data array and its dimensions */
  double ***vals=NULL;
  size_t blocks=0;
  size_t *rows=NULL,*columns=NULL;

  /* actual input specification */
  char *filename=strdup("");
  char *old_filename=strdup("");
  char *splitstring=strdup(" \t");

  /* list of specs */
  char **specs=NULL;
  size_t specsnum;
  size_t spec;

  /* output data arrays */
  double **data=NULL;
  size_t datanum = 0;
  size_t *length=NULL;

  /* final transformations */
  char *finaltransform=strdup("");


  /* OPTIONS */
  int o_binary_in=0;
  int o_binary_out=0;
  int o_verbose=0;
  
  /* ------------------------------------------ */
  int opt;
  /* ------------------------------------------ */

  /* COMMAND LINE PROCESSING */
      
  while((opt=getopt_long(argc,argv,"vhF:e:o:s:t:b:",gb_long_options, &gb_option_index))!=EOF){
    if(opt==0){
      gbutils_header(argv[0],stdout);
      exit(0);
    }
    else if(opt=='?'){
      fprintf(stderr,"option %c not recognized\n",optopt);
      exit(-1);
    }
    else if(opt=='v'){
      /*increase verbosity*/
      o_verbose=1;
    }
    else if(opt=='b'){
      /*set binary input or output*/

      switch(atoi(optarg)){
      case 0 :
	break;
      case 1 :
	o_binary_in=1;
	break;
      case 2 :
	o_binary_out=1;
	PRINTMATRIXBYCOLS=printmatrixbycols_bin;
	PRINTNOMATRIXBYCOLS=printnomatrixbycols_bin;
	break;
      case 3 :
	o_binary_in=1;
	o_binary_out=1;
	PRINTMATRIXBYCOLS=printmatrixbycols_bin;
	PRINTNOMATRIXBYCOLS=printnomatrixbycols_bin;
	break;
      default:
	fprintf(stderr,"unclear binary specification %d, option ignored.\n",atoi(optarg));
	break;
      }
    }
    else if(opt=='h'){
      /*print help*/
      fprintf(stdout,"Print slices of tabular data from files and apply transformations. Data\n");
      fprintf(stdout,"are read from text files with fields separated by space (use option -F\n");
      fprintf(stdout,"to specify a different separator). Inside data file, data-blocks are   \n");
      fprintf(stdout,"separated by two empty lines. File can be compressed with zlib (.gz).\n");
      fprintf(stdout,"Usage: %s [options] 'filename[index](C,R)trans'                      \n\n",argv[0]);
      fprintf(stdout," filename  is the input file. If not specified it default to stdin or  \n");
      fprintf(stdout,"           the last specified filename if any.\n");
      fprintf(stdout," index     stands for a data-block index.                              \n");
      fprintf(stdout," index     stands for a data-block index.                              \n");
      fprintf(stdout," C,R       stands for columns and rows spec given as \"min:max:skip\" to\n");
      fprintf(stdout,"           select from \"min\" to \"max\" every \"skip\" steps. If negative \n");
      fprintf(stdout,"           min and max are counted from the end. By default all data   \n");
      fprintf(stdout,"           are printed (\"1:-1:1\"). If min>max then count is reversed \n");
      fprintf(stdout,"           and skip must be negative (-1 by default). Different specs  \n");
      fprintf(stdout,"           are separated by semicolon ';' and considered sequentially. \n");
      fprintf(stdout," trans     is a list of transformations applied to selected data: 'd'  \n");
      fprintf(stdout,"           take the diff of subsequent columns; 'D' remove all rows with\n");
      fprintf(stdout,"           at least one Not-A-Number (NAN) entry; 'f' flatten the output\n");
      fprintf(stdout,"           piling all columns; 'l' take log of all entries, 'P' print all\n");
      fprintf(stdout,"           entries collected as a data-block; 't' transpose the matrix \n");
      fprintf(stdout,"           of data; 'z' subtract from the entries in each column their \n");
      fprintf(stdout,"           mean; 'Z' replace the entry in each column with their zscore;\n");
      fprintf(stdout,"           'w' divide the entry in each columns by their mean.          \n\n");
      fprintf(stdout,"           '<..;..>' functions separated by semicolons in angle brackets\n");
      fprintf(stdout,"           can be used for generic data transformation; the function is \n");
      fprintf(stdout,"           computed for each row of data. Variables names are 'x' followed\n");
      fprintf(stdout,"           by the number of the column and optionally by 'l' and the number\n");
      fprintf(stdout,"           of lags. For instance 'x2+x3l1' means the sum of the entries in\n");
      fprintf(stdout,"           the 2nd column plus the entries in the 3rd column in the previous\n");
      fprintf(stdout,"           row. 'x0' stands for the row number and 'x' is equal to 'x1'\n\n");
      fprintf(stdout,"           '<@..;..>' if the functions specification starts with a '@' the\n");
      fprintf(stdout,"           functions are computed recursively along the columns. In this\n");
      fprintf(stdout,"           case the number after the 'x' is the relative column counted  \n");
      fprintf(stdout,"           starting from the one considered at each step.              \n\n");
      fprintf(stdout,"           '{...}' a function in curly brackets can be use to select data:\n");
      fprintf(stdout,"           only rows that return a non-negative value are retained     \n");
      fprintf(stdout,"Options:                                                               \n");
      fprintf(stdout," -F  set the input fields separators (default ' \\t')                  \n");
      fprintf(stdout," -o  set the output format (default '%%12.6e')                         \n");
      fprintf(stdout," -e  set the output format for empty fields  (default '%%13s')         \n");
      fprintf(stdout," -s  set the output separation string  (default ' ')                   \n");
      fprintf(stdout," -t  define global transformations applied before each output (default '')\n");
      fprintf(stdout," -v  verbose mode\n");
      fprintf(stdout,"Examples:\n");
      fprintf(stdout," gbget 'file(1:3)ld'  select the first three columns in 'file', take the\n");
      fprintf(stdout,"                      log and the difference of successive columns;\n");
      fprintf(stdout," gbget 'file(2,-10:-1)  <x^2> select the last ten elements of the second'\n");
      fprintf(stdout,"                        of 'file' and print their squares\n");
      fprintf(stdout," gbget '[2]()' '[1]()' < ...  select the second and first data block from the\n");
      fprintf(stdout,"                              standard input.\n");
      fprintf(stdout," gbget 'file(1:3)<x1*x2-x3>'  select the first three columns in 'file' and\n");
      fprintf(stdout,"                              in each row multiply the first two entries and.\n");
      fprintf(stdout,"                              subtract the third.\n");
      fprintf(stdout," gbget 'file()<@x1+x2>'  print the sum of two subsequent columns\n");
      fprintf(stdout," gbget 'file(1:3){x2-2}'  select the first three columns in 'file' for the\n");
      fprintf(stdout,"                          rows whose second field is not lower then 2 \n");
      exit(0);
    }
    else if(opt=='F'){
      /*set the input fields separator string*/
      free(splitstring);
      splitstring = strdup(optarg);
    }
    else if(opt=='e'){
      /*set the format for empty fields*/
      EMPTY = strdup(optarg);
    }
    else if(opt=='o'){
      /*set the output string */
      FLOAT = strdup(optarg);
    }
    else if(opt=='s'){
      /*set the separator*/
      SEP = strdup(optarg);
    }
    else if(opt=='t'){
      /* set the final transformations */
      size_t len1 = strlen(optarg);
      size_t len2 = strlen(finaltransform);
      finaltransform  =  (char *) my_realloc((void *) finaltransform ,(len1+len2+1)*sizeof(char));
      strcpy(finaltransform+len2,optarg);
    }
  }
  /* END OF COMMAND LINE PROCESSING */

  /* initialize global variables */
  initialize_program(argv[0]);

  /* build the list of space and tabs separated specs */
  specsnum=0;
  for(spec=optind;spec<(size_t) argc;spec++){

    char *piece=strdup (argv[spec]);
    char *stmp1=piece,*stmp2;

    while( (stmp2=strsep (&stmp1," \t")) != NULL)
      if(strlen(stmp2) >0){
	specsnum++;
	specs  =  (char **) my_realloc((void *) specs,specsnum*sizeof(char *));
	specs[specsnum-1] = strdup(stmp2);
      }

    free(piece);
  }


  /* +++++++++++++++++++++++++++++++++++++++++++++ */
  if(o_verbose){
    fprintf(stderr,"Output format='%s' Separator='%s' Empty format='%s'  Empty string='%s'  ",
	    FLOAT,SEP,EMPTY," ");

    if(o_binary_in || o_binary_out){
      fprintf(stderr,"Binary:");
      if(o_binary_in)
	fprintf(stderr," in");
      if(o_binary_out)
	fprintf(stderr," out");
    }
    fprintf(stderr,"\n");
  }

  /* +++++++++++++++++++++++++++++++++++++++++++++ */


  /* parse line for data specification */
  /* --------------------------------- */    
  for(spec=0;spec<specsnum;spec++){

    /* added matrix */
    size_t rowlength,columnlength;
    double **newdata=NULL;

    /* different pieces of the spec */
    char *blockslice=NULL, *columnslice=NULL, *rowslice=NULL, *transform=NULL;
    
    parse_spec(specs[spec],&filename,&blockslice,&columnslice,&rowslice,&transform);

    /* set proper filename */
    if(strcmp(filename,"")==0 && strcmp(old_filename,"")!=0){
      free(filename);filename=strdup(old_filename);
    }

    /* fprintf(stderr,"old filename:%s, filename:%s, blockslice:%s, columnslice:%s, rowslice:%s, transform:%s,\n",old_filename,filename,blockslice,columnslice,rowslice,transform);   */ 
    /* return 0;  */

    /* load new data if the file changed */
    if( spec==0 || ( strcmp(filename,old_filename) !=0 )){
      int file;

      /*open file*/
      if(strcmp(filename,"")==0) file=0;
      else if( (file = open(filename,O_RDONLY)) == -1){
	fprintf(stderr,"WARNING (%s): failed to open file \"%s\"; request ignored\n",GB_PROGNAME,filename);
	continue;
      }

      
      /* deallocate if necessary */
      if(vals != NULL){
	size_t i,j ;
	for(i=0;i<blocks;i++){
	  for(j=0;j<columns[i];j++)
	    free(vals[i][j]);
	  free(vals[i]);
	}
	free(vals);
	free(rows);
	free(columns);
	blocks=0;
	rows=NULL;
	columns=NULL;
	vals=NULL;
      }

      /*read from file*/
      if(o_binary_in==1)
	loadblocks_bin(&vals,&blocks,&rows,&columns,file);
      else
	loadblocks(&vals,&blocks,&rows,&columns,file,splitstring);

      /* store values */
      free(old_filename);old_filename=strdup(filename);

      /* close file */
      if(file==0)
	lseek(file,0,SEEK_SET);
      close(file);

      /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
      if(o_verbose){
	size_t b;
	fprintf(stderr,"<- source:%s",(strcmp(filename,"")==0?"(stdin)":filename));
	fprintf(stderr,"\n");
	for(b=0;b<blocks;b++)
	  fprintf(stderr,"   block:%zd rows:%zd cols:%zd\n",b,rows[b],columns[b]);	
	fprintf(stderr,"\n");
      }
      /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

    }


    /* insert in turn any matrix specified in the block spec */
    {

      size_t bslicenum=0;
      size_t b,bs;
      size_t *bini=NULL;
      size_t *bfin=NULL;
      int *bskip=NULL;

      /* build array of block slice */
      parse_slice(blockslice,"block",blocks,&bslicenum,&bini,&bfin,&bskip);

      for(bs=0;bs<bslicenum;bs++){

	/* for(b=bini[bs]-1; bskipsign*(b- bfin[bs]+1) <= 0;b+=bskip[bs]){ */
	for(b=bini[bs]-1; ( bskip[bs]>0 ? b < bfin[bs] : b+2 > bfin[bs] );b+=bskip[bs]){

	  /* create the new matrix */
	  new_matrix(&newdata,&columnlength,&rowlength,
		     vals[b],columns[b],rows[b],columnslice,rowslice);


	  /* apply transformations to the new matrix */
	  if(strlen(transform)>0 && strcmp(transform,"P") )
	    matrix_transform(&newdata,&columnlength,&rowlength,
			     transform,o_verbose);

	  /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
	  if(o_verbose==1)
	    fprintf(stderr,"-> %zdx%zd\n",rowlength,columnlength);
	  /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
	  
	  /* add the new matrix to the data */
	  add_matrix(&data,&length,&datanum,
		     newdata,columnlength,rowlength);
	  
	  /* if P transformation was specified then print data collected so
	     far and restart to collect them */
	  if(index(transform,'P') != NULL){
	    
	    /* final transformation and print */
	    print_data(&datanum, &length, &data, finaltransform,o_verbose);
	    
	    /* add two empty lines in ASCII output */
	    if(o_binary_out==0)
	      printf("\n\n");
	    
	  }
	}
      }
      
      free(bini);
      free(bfin);
      free(bskip);

    }
    
    /* free everything */
    free(columnslice);
    free(rowslice);
    free(blockslice);
    free(transform);
    free(newdata);


  }/* END OF DATA CREATION */

  /* free some space */
  free(old_filename);
  free(filename);
  free(splitstring);

  if(specsnum>0){
    size_t i;
    for(i=0;i<specsnum;i++) free(specs[i]);
    free(specs);
  }

  if(vals != NULL){
    size_t i,j ;
    for(i=0;i<blocks;i++){
      for(j=0;j<columns[i];j++)
	free(vals[i][j]);
      free(vals[i]);
      }
    free(vals);
    free(rows);
    free(columns);
    blocks=0;
  }

  /* final transformation and print*/

  /* print */
  print_data(&datanum, &length, &data, finaltransform,o_verbose);


  {
    size_t i;

    for(i=0;i<datanum;i++) free(data[i]);
    free(length);
    free(data);
  }

  free(finaltransform);

  finalize_program();

  return 0 ;
}