Mercurial > repos > xuebing > sharplabtool

comparison tools/plotting/histogram.py @ 0:9071e359b9a3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9071e359b9a3
1 #!/usr/bin/env python
2 #Greg Von Kuster
3
4 import sys
5 from rpy import *
6
7 assert sys.version_info[:2] >= ( 2, 4 )
8
9 def stop_err(msg):
10 sys.stderr.write(msg)
11 sys.exit()
12
13 def main():
14
15 # Handle input params
16 in_fname = sys.argv[1]
17 out_fname = sys.argv[2]
18 try:
19 column = int( sys.argv[3] ) - 1
20 except:
21 stop_err( "Column not specified, your query does not contain a column of numerical data." )
22 title = sys.argv[4]
23 xlab = sys.argv[5]
24 breaks = int( sys.argv[6] )
25 if breaks == 0:
26 breaks = "Sturges"
27 if sys.argv[7] == "true":
28 density = True
29 else: density = False
30 if len( sys.argv ) >= 9 and sys.argv[8] == "true":
31 frequency = True
32 else: frequency = False
33
34 matrix = []
35 skipped_lines = 0
36 first_invalid_line = 0
37 invalid_value = ''
38 i = 0
39 for i, line in enumerate( file( in_fname ) ):
40 valid = True
41 line = line.rstrip('\r\n')
42 # Skip comments
43 if line and not line.startswith( '#' ):
44 # Extract values and convert to floats
45 row = []
46 try:
47 fields = line.split( "\t" )
48 val = fields[column]
49 if val.lower() == "na":
50 row.append( float( "nan" ) )
51 except:
52 valid = False
53 skipped_lines += 1
54 if not first_invalid_line:
55 first_invalid_line = i+1
56 else:
57 try:
58 row.append( float( val ) )
59 except ValueError:
60 valid = False
61 skipped_lines += 1
62 if not first_invalid_line:
63 first_invalid_line = i+1
64 invalid_value = fields[column]
65 else:
66 valid = False
67 skipped_lines += 1
68 if not first_invalid_line:
69 first_invalid_line = i+1
70
71 if valid:
72 matrix += row
73
74 if skipped_lines < i:
75 try:
76 a = r.array( matrix )
77 r.pdf( out_fname, 8, 8 )
78 histogram = r.hist( a, probability=not frequency, main=title, xlab=xlab, breaks=breaks )
79 if density:
80 density = r.density( a )
81 if frequency:
82 scale_factor = len( matrix ) * ( histogram['mids'][1] - histogram['mids'][0] ) #uniform bandwidth taken from first 2 midpoints
83 density[ 'y' ] = map( lambda x: x * scale_factor, density[ 'y' ] )
84 r.lines( density )
85 r.dev_off()
86 except Exception, exc:
87 stop_err( "%s" %str( exc ) )
88 else:
89 if i == 0:
90 stop_err("Input dataset is empty.")
91 else:
92 stop_err( "All values in column %s are non-numeric." %sys.argv[3] )
93
94 print "Histogram of column %s. " %sys.argv[3]
95 if skipped_lines > 0:
96 print "Skipped %d invalid lines starting with line #%d, '%s'." % ( skipped_lines, first_invalid_line, invalid_value )
97
98 r.quit( save="no" )
99
100 if __name__ == "__main__":
101 main()