annotate execute_dwt_var_perFeature.pl @ 0:d56c5d2e1a29 draft

Imported from capsule None
author devteam
date Mon, 27 Jan 2014 09:26:33 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
1 #!/usr/bin/perl -w
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
2 # Author: Erika Kvikstad
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
3
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
4 use warnings;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
5 use IO::Handle;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
6 use POSIX qw(floor ceil);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
7
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
8 $usage = "execute_dwt_var_perFeature.pl [TABULAR.in] [FEATURE] [ALPHA] [TABULAR.out] [PDF.out] \n";
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
9 die $usage unless @ARGV == 5;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
10
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
11 #get the input arguments
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
12 my $inputFile = $ARGV[0];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
13 my @features = split(/,/,$ARGV[1]);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
14 my $features_count = scalar(@features);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
15 my $alpha = $ARGV[2];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
16 my $outFile1 = $ARGV[3];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
17 my $outFile2 = $ARGV[4];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
18
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
19 open (INPUT, "<", $inputFile) || die("Could not open file $inputFile \n");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
20 open (OUTPUT2, ">", $outFile1) || die("Could not open file $outFile1 \n");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
21 open (OUTPUT3, ">", $outFile2) || die("Could not open file $outFile2 \n");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
22 #open (ERROR, ">", "error.txt") or die ("Could not open file error.txt \n");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
23
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
24 # choosing meaningful names for the output files
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
25 $pvalue = $outFile1;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
26 $pdf = $outFile2;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
27
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
28 # write R script
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
29 $r_script = "get_dwt_varPermut.r";
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
30
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
31 open(Rcmd, ">", "$r_script") or die "Cannot open $r_script \n\n";
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
32
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
33 print Rcmd "
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
34 ######################################################################
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
35 # plot multiscale wavelet variance
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
36 # create null bands by permuting the original data series
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
37 # generate plots and table of wavelet variance including p-values
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
38 ######################################################################
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
39 options(echo = FALSE)
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
40 #library(\"Rwave\");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
41 #library(\"wavethresh\");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
42 #library(\"waveslim\");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
43 # turn off diagnostics for de-bugging only, turn back on for functional tests on test
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
44 suppressMessages(require(\"Rwave\",quietly=TRUE,warn.conflicts = FALSE));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
45 suppressMessages(require(\"wavethresh\",quietly=TRUE,warn.conflicts = FALSE));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
46 suppressMessages(require(\"waveslim\",quietly=TRUE,warn.conflicts = FALSE));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
47 suppressMessages(require(\"bitops\",quietly=TRUE,warn.conflicts = FALSE));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
48
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
49 # to determine if data is properly formatted 2^N observations
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
50 is.power2<- function(x){x && !(bitAnd(x,x - 1));}
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
51
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
52 # dwt : discrete wavelet transform using Haar wavelet filter, simplest wavelet function but later can modify to let user-define the wavelet filter function
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
53 dwt_var_permut_getMax <- function(data, names, alpha, filter = 1,family=\"DaubExPhase\", bc = \"symmetric\", method = \"kendall\", wf = \"haar\", boundary = \"reflection\") {
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
54 max_var = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
55 matrix = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
56 title = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
57 final_pvalue = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
58 J = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
59 scale = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
60 out = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
61
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
62 print(class(data));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
63 print(names);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
64 print(alpha);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
65
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
66 par(mar=c(5,4,4,3),oma = c(4, 4, 3, 2), xaxt = \"s\", cex = 1, las = 1);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
67
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
68 title<-c(\"Wavelet\",\"Variance\",\"Pvalue\",\"Test\");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
69 print(title);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
70
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
71 for(i in 1:length(names)){
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
72 temp = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
73 results = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
74 wave1.dwt = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
75
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
76 # if data fails formatting check, do something
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
77
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
78 print(is.numeric(as.matrix(data)[, i]));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
79 if(!is.numeric(as.matrix(data)[, i]))
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
80 stop(\"data must be a numeric vector\");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
81
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
82 print(length(as.matrix(data)[, i]));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
83 print(is.power2(length(as.matrix(data)[, i])));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
84 if(!is.power2(length(as.matrix(data)[, i])))
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
85 stop(\"data length must be a power of two\");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
86
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
87
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
88 J <- wd(as.matrix(data)[, i], filter.number = filter, family=family, bc = bc)\$nlevels;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
89 print(J);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
90 temp <- vector(length = J);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
91 wave1.dwt <- dwt(as.matrix(data)[, i], wf = wf, J, boundary = boundary);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
92 #print(wave1.dwt);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
93
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
94 temp <- wave.variance(wave1.dwt)[-(J+1), 1];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
95 print(temp);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
96
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
97 #permutations code :
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
98 feature1 = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
99 null = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
100 var_lower=limit_lower=NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
101 var_upper=limit_upper=NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
102 med = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
103
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
104 limit_lower = alpha/2*1000;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
105 print(limit_lower);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
106 limit_upper = (1-alpha/2)*1000;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
107 print(limit_upper);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
108
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
109 feature1 = as.matrix(data)[,i];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
110 for (k in 1:1000) {
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
111 nk_1 = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
112 null.levels = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
113 var = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
114 null_wave1 = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
115
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
116 nk_1 = sample(feature1, length(feature1), replace = FALSE);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
117 null.levels <- wd(nk_1, filter.number = filter,family=family ,bc = bc)\$nlevels;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
118 var <- vector(length = length(null.levels));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
119 null_wave1 <- dwt(nk_1, wf = wf, J, boundary = boundary);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
120 var<- wave.variance(null_wave1)[-(null.levels+1), 1];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
121 null= rbind(null, var);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
122 }
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
123 null <- apply(null, 2, sort, na.last = TRUE);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
124 var_lower <- null[limit_lower, ];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
125 var_upper <- null[limit_upper, ];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
126 med <- (apply(null, 2, median, na.rm = TRUE));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
127
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
128 # plot
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
129 results <- cbind(temp, var_lower, var_upper);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
130 print(results);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
131 matplot(results, type = \"b\", pch = \"*\", lty = 1, col = c(1, 2, 2),xaxt='n',xlab=\"Wavelet Scale\",ylab=\"Wavelet variance\" );
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
132 mtext(names[i], side = 3, line = 0.5, cex = 1);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
133 axis(1, at = 1:J , labels=c(2^(0:(J-1))), las = 3, cex.axis = 1);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
134
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
135 # get pvalues by comparison to null distribution
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
136 #out <- (names[i]);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
137 for (m in 1:length(temp)){
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
138 print(paste(\"scale\", m, sep = \" \"));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
139 print(paste(\"var\", temp[m], sep = \" \"));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
140 print(paste(\"med\", med[m], sep = \" \"));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
141 pv = tail =scale = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
142 scale=2^(m-1);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
143 #out <- c(out, format(temp[m], digits = 3));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
144 if (temp[m] >= med[m]){
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
145 # R tail test
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
146 print(\"R\");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
147 tail <- \"R\";
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
148 pv <- (length(which(null[, m] >= temp[m])))/(length(na.exclude(null[, m])));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
149
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
150 } else {
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
151 if (temp[m] < med[m]){
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
152 # L tail test
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
153 print(\"L\");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
154 tail <- \"L\";
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
155 pv <- (length(which(null[, m] <= temp[m])))/(length(na.exclude(null[, m])));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
156 }
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
157 }
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
158 print(pv);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
159 out<-rbind(out,c(paste(\"Scale\", scale, sep=\"_\"),format(temp[m], digits = 3),pv,tail));
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
160 }
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
161 final_pvalue <-rbind(final_pvalue, out);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
162 }
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
163 colnames(final_pvalue) <- title;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
164 return(final_pvalue);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
165 }\n";
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
166
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
167 print Rcmd "
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
168 # execute
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
169 # read in data
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
170 data_test = final = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
171 sub = sub_names = NULL;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
172 data_test <- read.delim(\"$inputFile\",header=FALSE);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
173 pdf(file = \"$pdf\", width = 11, height = 8)\n";
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
174
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
175 for ($x=0;$x<$features_count;$x++){
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
176 $feature=$features[$x];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
177 print Rcmd "
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
178 if ($feature > ncol(data_test))
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
179 stop(\"column $feature doesn't exist\");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
180 sub<-data_test[,$feature];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
181 #sub_names <- colnames(data_test);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
182 sub_names<-colnames(data_test)[$feature];
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
183 final <- rbind(final,dwt_var_permut_getMax(sub, sub_names,$alpha));\n";
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
184 }
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
185
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
186 print Rcmd "
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
187
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
188 dev.off();
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
189 write.table(final, file = \"$pvalue\", sep = \"\\t\", quote = FALSE, row.names = FALSE);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
190
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
191 #eof\n";
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
192
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
193 close Rcmd;
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
194 system("R --no-restore --no-save --no-readline < $r_script > $r_script.out");
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
195
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
196 #close the input and output and error files
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
197 close(OUTPUT3);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
198 close(OUTPUT2);
d56c5d2e1a29 Imported from capsule None
devteam
parents:
diff changeset
199 close(INPUT);