Mercurial > repos > devteam > canonical_correlation_analysis

comparison cca.py @ 0:9bc0c48a027f draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Imported from capsule None
author devteam
date Mon, 19 May 2014 12:34:48 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9bc0c48a027f
1 #!/usr/bin/env python
2
3 import sys, string
4 from rpy import *
5 import numpy
6
7 def stop_err(msg):
8 sys.stderr.write(msg)
9 sys.exit()
10
11 infile = sys.argv[1]
12 x_cols = sys.argv[2].split(',')
13 y_cols = sys.argv[3].split(',')
14
15 x_scale = x_center = "FALSE"
16 if sys.argv[4] == 'both':
17 x_scale = x_center = "TRUE"
18 elif sys.argv[4] == 'center':
19 x_center = "TRUE"
20 elif sys.argv[4] == 'scale':
21 x_scale = "TRUE"
22
23 y_scale = y_center = "FALSE"
24 if sys.argv[5] == 'both':
25 y_scale = y_center = "TRUE"
26 elif sys.argv[5] == 'center':
27 y_center = "TRUE"
28 elif sys.argv[5] == 'scale':
29 y_scale = "TRUE"
30
31 std_scores = "FALSE"
32 if sys.argv[6] == "yes":
33 std_scores = "TRUE"
34
35 outfile = sys.argv[7]
36 outfile2 = sys.argv[8]
37
38 fout = open(outfile,'w')
39 elems = []
40 for i, line in enumerate( file ( infile )):
41 line = line.rstrip('\r\n')
42 if len( line )>0 and not line.startswith( '#' ):
43 elems = line.split( '\t' )
44 break
45 if i == 30:
46 break # Hopefully we'll never get here...
47
48 if len( elems )<1:
49 stop_err( "The data in your input dataset is either missing or not formatted properly." )
50
51 x_vals = []
52
53 for k,col in enumerate(x_cols):
54 x_cols[k] = int(col)-1
55 x_vals.append([])
56
57 y_vals = []
58
59 for k,col in enumerate(y_cols):
60 y_cols[k] = int(col)-1
61 y_vals.append([])
62
63 skipped = 0
64 for ind,line in enumerate( file( infile )):
65 if line and not line.startswith( '#' ):
66 try:
67 fields = line.strip().split("\t")
68 valid_line = True
69 for col in x_cols+y_cols:
70 try:
71 assert float(fields[col])
72 except:
73 skipped += 1
74 valid_line = False
75 break
76 if valid_line:
77 for k,col in enumerate(x_cols):
78 try:
79 xval = float(fields[col])
80 except:
81 xval = NaN#
82 x_vals[k].append(xval)
83 for k,col in enumerate(y_cols):
84 try:
85 yval = float(fields[col])
86 except:
87 yval = NaN#
88 y_vals[k].append(yval)
89 except:
90 skipped += 1
91
92 x_vals1 = numpy.asarray(x_vals).transpose()
93 y_vals1 = numpy.asarray(y_vals).transpose()
94
95 x_dat= r.list(array(x_vals1))
96 y_dat= r.list(array(y_vals1))
97
98 try:
99 r.suppressWarnings(r.library("yacca"))
100 except:
101 stop_err("Missing R library yacca.")
102
103 set_default_mode(NO_CONVERSION)
104 try:
105 xcolnames = ["c%d" %(el+1) for el in x_cols]
106 ycolnames = ["c%d" %(el+1) for el in y_cols]
107 cc = r.cca(x=x_dat, y=y_dat, xlab=xcolnames, ylab=ycolnames, xcenter=r(x_center), ycenter=r(y_center), xscale=r(x_scale), yscale=r(y_scale), standardize_scores=r(std_scores))
108 ftest = r.F_test_cca(cc)
109 except RException, rex:
110 stop_err("Encountered error while performing CCA on the input data: %s" %(rex))
111
112 set_default_mode(BASIC_CONVERSION)
113 summary = r.summary(cc)
114
115 ncomps = len(summary['corr'])
116 comps = summary['corr'].keys()
117 corr = summary['corr'].values()
118 xlab = summary['xlab']
119 ylab = summary['ylab']
120
121 for i in range(ncomps):
122 corr[comps.index('CV %s' %(i+1))] = summary['corr'].values()[i]
123
124 ftest=ftest.as_py()
125 print >>fout, "#Component\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
126 print >>fout, "#Correlation\t%s" %("\t".join(["%.4g" % el for el in corr]))
127 print >>fout, "#F-statistic\t%s" %("\t".join(["%.4g" % el for el in ftest['statistic']]))
128 print >>fout, "#p-value\t%s" %("\t".join(["%.4g" % el for el in ftest['p.value']]))
129
130 print >>fout, "#X-Coefficients\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
131 for i,val in enumerate(summary['xcoef']):
132 print >>fout, "%s\t%s" %(xlab[i], "\t".join(["%.4g" % el for el in val]))
133
134 print >>fout, "#Y-Coefficients\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
135 for i,val in enumerate(summary['ycoef']):
136 print >>fout, "%s\t%s" %(ylab[i], "\t".join(["%.4g" % el for el in val]))
137
138 print >>fout, "#X-Loadings\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
139 for i,val in enumerate(summary['xstructcorr']):
140 print >>fout, "%s\t%s" %(xlab[i], "\t".join(["%.4g" % el for el in val]))
141
142 print >>fout, "#Y-Loadings\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
143 for i,val in enumerate(summary['ystructcorr']):
144 print >>fout, "%s\t%s" %(ylab[i], "\t".join(["%.4g" % el for el in val]))
145
146 print >>fout, "#X-CrossLoadings\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
147 for i,val in enumerate(summary['xcrosscorr']):
148 print >>fout, "%s\t%s" %(xlab[i], "\t".join(["%.4g" % el for el in val]))
149
150 print >>fout, "#Y-CrossLoadings\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
151 for i,val in enumerate(summary['ycrosscorr']):
152 print >>fout, "%s\t%s" %(ylab[i], "\t".join(["%.4g" % el for el in val]))
153
154 r.pdf( outfile2, 8, 8 )
155 #r.plot(cc)
156 for i in range(ncomps):
157 r.helio_plot(cc, cv = i+1, main = r.paste("Explained Variance for CV",i+1), type = "variance")
158 r.dev_off()