src/tools/average.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdio.h>
  40 #include <math.h>
  41 #include "string2.h"
  42 #include "smalloc.h"
  43 #include "statutil.h"
  44 #include "vec.h"
  45
  46 static void my_usage(char *prog,char *arg)
  47 {
  48   fprintf(stderr," Usage: %s [-p] [-s] [-c #columns]"
  49           " naver < infile > outfile\n",prog);
  50   fprintf(stderr,"-p means picoseconds rather than nanoseconds\n");
  51   fprintf(stderr,"-s means silent rather than verbose\n");
  52   fprintf(stderr,"Don't ever use argument %s again!\n",arg);
  53   exit(1);
  54 }
  55
  56 void lsq_y_ax_b_double(int n, double x[], double y[], double *a, double *b)
  57 {
  58   int    i;
  59   double yx,xx,sx,sy;
  60
  61   yx=xx=sx=sy=0.0;
  62   for (i=0; (i < n); i++) {
  63     yx+=y[i]*x[i];
  64     xx+=x[i]*x[i];
  65     sx+=x[i];
  66     sy+=y[i];
  67   }
  68   *a=(n*yx-sy*sx)/(n*xx-sx*sx);
  69   *b=(sy-(*a)*sx)/n;
  70 }
  71
  72 int main(int argc, char *argv[])
  73 {
  74   double *x,**y,value=0.001;
  75   double *yav,yyy,yyy2,ymin,ymax,aver,var,sd;
  76   int    i,j,k,nav=100,ncol=1,MAX=50000;
  77   char   buf[STRLEN];
  78   bool   bSilent=FALSE,bVerySilent=FALSE;
  79   double   lsq_a,lsq_b,rms_res;
  80
  81   for(i=1; (i<argc); i++) {
  82     if (argv[i][0] == '-')
  83       switch (argv[i][1]) {
  84       case 'p':
  85         value=1.0;
  86         break;
  87       case 'm':
  88         MAX=iscan(argc,argv,&i);
  89         break;
  90       case 'c':
  91         ncol=iscan(argc,argv,&i);
  92         break;
  93       case 's':
  94         bSilent=TRUE;
  95         break;
  96       case 'S':
  97         bVerySilent=bSilent=TRUE;
  98         break;
  99       default:
 100         my_usage(argv[0],argv[i]);
 101       }
 102     else
 103       nav=strtol(argv[i],NULL,10);
 104   }
 105   if (!bSilent)
 106     fprintf(stderr,"Will average stdin with %d columns, over %d points\n",
 107             ncol,nav);
 108
 109
 110   snew(x,MAX);
 111   snew(y,ncol);
 112   for(i=0; (i<ncol); i++) {
 113     snew(y[i],MAX);
 114   }
 115   snew(yav,MAX);
 116   i=0;
 117   do {
 118     if (scanf("%s",buf) == 1) {
 119       if ((buf[0] != '@') && (buf[0] != '#')) {
 120         sscanf(buf,"%lf",&x[i]);
 121         for(k=0; (k<ncol); k++)
 122           scanf("%lf",&(y[k][i]));
 123         if (i >= nav) {
 124           if (!bVerySilent)
 125             printf("%10e",x[i-nav/2]*value);
 126           for(k=0; (k<ncol); k++) {
 127             yav[k]=0;
 128             for(j=i-nav+1; (j<=i); j++)
 129               yav[k]+=y[k][j];
 130             yav[k]/=nav;
 131             if (!bVerySilent)
 132               printf("  %10e",yav[k]);
 133           }
 134           if (!bVerySilent)
 135             printf("\n");
 136         }
 137         i++;
 138       }
 139       else {
 140         while (getc(stdin) != '\n')
 141           ;
 142       }
 143     }
 144     else
 145       break;
 146   } while (((int)strlen(buf) > 0) && (i<MAX));
 147
 148
 149   if (!bSilent)
 150     fprintf(stderr,"%3s  %12s  %12s  %12s  %12s  %12s  %12s  %12s\n",
 151             "i","min","aver","max","var","sd","drift","fluc");
 152   for(k=0; (k<ncol); k++) {
 153     aver=y[k][0];
 154     ymin=aver;
 155     ymax=aver;
 156     yyy2=aver*aver;
 157     for(j=1; (j<i); j++) {
 158       yyy=y[k][j];
 159       aver+=yyy;
 160       yyy2+=(yyy*yyy);
 161       if (yyy < ymin)
 162         ymin=yyy;
 163       else if (yyy > ymax)
 164         ymax=yyy;
 165     }
 166     aver/=i;
 167     var=yyy2/i-aver*aver;
 168     sd=sqrt(var);
 169
 170     lsq_y_ax_b_double(i,x,y[k],&lsq_a,&lsq_b);
 171     rms_res=0;
 172     for(j=0; (j<i);j++)
 173       rms_res+=sqr(y[k][j]-(lsq_a*x[j]+lsq_b));
 174     rms_res=sqrt(rms_res/i);
 175
 176     fprintf(stderr,"%3d  %12g  %12g  %12g  %12g  %12g  %12g  %12g\n",
 177             k,ymin,aver,ymax,var,sd,lsq_a,rms_res);
 178   }
 179   return 0;
 180 }
 181
 182