src/tools/average.c

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