#!/usr/bin/python # -*- coding: utf-8 -*- # #------------------------------------------------------------------------------ #anacolt.py, Copyright Bjoern Olausson #------------------------------------------------------------------------------ #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. # #To view the license visit #http://www.gnu.org/licenses/old-licenses/gpl-2.0.html #or write to #Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA #02110-1301 USA #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ # #This script is intended to pars NAMD colvar trajectory files and plot/extract #data from any collumn you like # __version__ = '1.2.0' __version_info__ = tuple([int(num) for num in __version__.split('.')]) import sys import re import glob import itertools import pylab import math import numpy as np from smooth import smooth from optparse import OptionParser, OptionGroup usage = '''Run without -x and -y to see all possible values for the input file.''' opa = OptionParser(usage=usage, version=__version__) opa.add_option("-f", action="append", dest="coltf", metavar="FILE", help="NAMD colvar trajectory file") opa.add_option("-o", action="store", dest="saveplot", metavar="FILE", help="Filename to which we save the plot (.png) and RAW data (.dat)") opa.add_option("-c", action="store", dest="concat", metavar="FILE", help="Concatenate all X/Y data in one file") opa.add_option("-s", action="store_true", dest="singleplot", help="Plot each X/Y pair to a standallone image") opa.add_option("-m", action="store_true", dest="singlefigure", help="Plot each X/Y pair to single figure") opa.add_option("-j", action="store", dest="sm", type="int", help="Smooth data over x points") group = OptionGroup(opa, "You can choose freely from the presented keywords \ which will occure on either X or Y Axis") group.add_option("-x", action="append", dest="xax", help="Set collumn header which data to appear on the X-axis.") group.add_option("-y", action="append", dest="yax", help="Set collumn header which data to appear on the Y-axis. Use \"*\" to \ plot all collumns against a _single_ X collumn. If you want to plot multible Y\ values, you have to asign every Y value a X value and vice versa. For example:\ anaout.py -x step -y colvar_name1 -y step -y colvar_name2") #group.add_option("-l", action="store", dest="yaxl", #help="Use Python list expression to select collumns for the Y-Axis. E. g.\ # Plot every third collumn starting from the second to the end: 0:,1::3") opa.add_option_group(group) # instruct optparse to parse the program's command line (options, args) = opa.parse_args() # Checking for required otions if not options.coltf: opa.print_help() opa.error("-f must be set") coltf = options.coltf xax = options.xax yax = options.yax sm = options.sm concat = options.concat saveplot = options.saveplot singleplot = options.singleplot singlefigure = options.singlefigure # Function for natural sorting! def sort_nicely(l): """ Sort the given list in the way that humans expect. http://www.codinghorror.com/blog/2007/12/ sorting-for-humans-natural-sort-order.html""" convert = lambda text: int(text) if text.isdigit() else text alphanum_key = lambda key: [convert(c) for c in re.split('([0-9]+)', key)] l.sort(key=alphanum_key) return l filelist_unsorted = [] for args in coltf: expanded = glob.glob(args) filelist_unsorted.extend(expanded) filelist = sort_nicely(filelist_unsorted) f = filelist[0] # Extract the collumn header fo = open(f) header = fo.readline().strip("#").split() fo.close() if not xax or not yax: print "Available collumns" print "------------------------------------------" print ", ".join(header) print "------------------------------------------" sys.exit(0) if yax[0] == '*': xindex = header.index(xax[0]) yc = range(len(header)) del yc[xindex] xc = [xindex] * len(yc) else: xc = [] for item in xax: xc.append(header.index(item)) yc = [] for item in yax: yc.append(header.index(item)) print "Loading %(file)s" % {"file": f} a = np.loadtxt(f, dtype='float', comments='#') #pylab.plot(a[0:,0], a[yaxl]) #pylab.show() if len(filelist) > 1: for f in filelist[1:]: print "Loading %(file)s" % {"file": f} b = np.loadtxt(f, dtype='float', comments='#') c = b[1:, :] a = np.vstack((a, c)) if sm: shp = a.shape for i in range(shp[1]): collumn = a[:, i] smoothed = np.array(smooth(collumn, sm)) try: asm = np.vstack((asm, smoothed)) except NameError, e: asm = smoothed a = asm.transpose() # Print the average for the selected values print "Some stats for the selected values:" for xAxis, yAxis in itertools.izip(xc, yc): X = a[:, xAxis] Y = a[:, yAxis] Xname = header[xAxis] Yname = header[yAxis] print Yname print("MIN: %(MIN)s@%(FMIN)s \ Max: %(MAX)s@%(FMAX)s \ MEAN: %(MEAN)s \ STD: %(STD)s | \ SUM: %(SUM)s" % {"MIN": Y.min(), "FMIN": X[Y.argmin()], "MAX": Y.max(), "FMAX": X[Y.argmax()], "MEAN": Y.mean(), "STD": Y.std(), "SUM": Y.sum() } ) if concat: for xAxis, yAxis in itertools.izip(xc, yc): X = a[:, xAxis].reshape(-1,1) Y = a[:, yAxis].reshape(-1,1) Xname = header[xAxis] Yname = header[yAxis] XY = np.hstack((X, Y)) try: cc = np.hstack((cc, XY)) except NameError, e: cc = XY try: cn.extend([Xname, Yname]) except NameError, e: cn = [Xname, Yname] fout = open(concat, "w") fout.writelines("%s\n" % (" ".join(cn))) for i in cc: j = [str(s) for s in i.tolist()] l = " ".join(j) fout.writelines("%s\n" % (l)) fout.close() def SaveData(Xvalues, Yvalues, Xname, Yname, saveplot): if saveplot: savedata = str(Xname) + "-" + str(Yname) + "_" + str(saveplot) + ".dat" fout = open(savedata, "w") for i, j in itertools.izip(Xvalues, Yvalues): fout.writelines("%s %s\n" % (i, j)) fout.close() def SaveImage(Xname, Yname, saveplot): if saveplot: saveimg_png = str(Xname) + "-" + str(Yname) + "_" + str(saveplot) + ".png" saveimg_svg = str(Xname) + "-" + str(Yname) + "_" + str(saveplot) + ".svg" pylab.savefig(saveimg_png) pylab.savefig(saveimg_svg) def DynaPltSingle(Xvalues, Yvalues, Xname, Yname): SaveData(Xvalues, Yvalues, Xname, Yname, saveplot) pylab.plot(Xvalues, Yvalues) Xlegend.append(Xname) Ylegend.append(Yname) def DynaPltMulti(Xvalues, Yvalues, Xname, Yname, figure): pylab.figure(figure) pylab.plot(Xvalues, Yvalues) pylab.xlabel(Xname) pylab.ylabel(Yname) SaveImage(Xname, Yname, saveplot) SaveData(Xvalues, Yvalues, Xname, Yname, saveplot) def DynaPltFig(Xvalues, Yvalues, Xname, Yname, ax): ax.plot(Xvalues, Yvalues) pylab.xlabel(Xname) pylab.ylabel(Yname) SaveData(Xvalues, Yvalues, Xname, Yname, saveplot) fign = 0 Xlegend = [] Xlabel = [] Ylegend = [] Ylabel = [] if singleplot: for xAxis, yAxis in itertools.izip(xc, yc): DynaPltSingle(a[:, xAxis], a[:, yAxis], header[xAxis], header[yAxis]) pylab.legend((Ylegend), shadow=True) pylab.ylabel(Ylegend) pylab.xlabel(Xlegend) SaveImage(xAxis, 'MULTI', saveplot) elif singlefigure: fig = pylab.figure() figurex = 1 nplots = len(yc) nsubp = int(math.ceil(math.sqrt(nplots))) for xAxis, yAxis in itertools.izip(xc, yc): ax = fig.add_subplot(nsubp, nsubp, figurex) DynaPltFig(a[:, xAxis], a[:, yAxis], header[xAxis], header[yAxis], ax) figurex += 1 SaveImage(xAxis, 'MULTI', saveplot) else: figure = 0 for xAxis, yAxis in itertools.izip(xc, yc): DynaPltMulti(a[:, xAxis], a[:, yAxis], header[xAxis], header[yAxis], figure) figure += 1 pylab.show()