Mercurial > repos > greg > ks_distribution

comparison ks_distribution.R @ 32:4c3fcfd55d94 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Mon, 26 Jun 2017 10:44:53 -0400
parents 812f09b96a62
children 5652711d5ae4
comparison
equal deleted inserted replaced
31:812f09b96a62 32:4c3fcfd55d94
13 13
14 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 14 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
15 args <- parse_args(parser, positional_arguments=TRUE) 15 args <- parse_args(parser, positional_arguments=TRUE)
16 opt <- args$options 16 opt <- args$options
17 17
18 set_component_colors = function(colors, number_comp)
19 {
20 # Handle colors for components.
21 if (is.na(colors))
22 {
23 # Randomly specify colors for components.
24 specified_colors <- c('red', 'yellow', 'green', 'black', 'blue', 'darkorange')
25 }
26 else
27 {
28 # Handle selected colors for components.
29 parser <- newJSONParser()
30 parser$addData(colors)
31 raw_colors <- parser$getObject()
32 specified_colors <- c()
33 for (raw_color in raw_colors)
34 {
35 specified_colors <- c(specified_colors, raw_color)
36 }
37 num_colors_specified = length(specified_colors)
38 if (num_colors_specified < number_comp)
39 {
40 for (i in num_colors_specified:number_comp)
41 {
42 if (!any(specified_colors=='red'))
43 {
44 specified_colors <- c(specified_colors, 'red')
45 }
46 else if (!any(specified_colors=='yellow'))
47 {
48 specified_colors <- c(specified_colors, 'yellow')
49 }
50 else if (!any(specified_colors=='green'))
51 {
52 specified_colors <- c(specified_colors, 'green')
53 }
54 else if (!any(specified_colors=='black'))
55 {
56 specified_colors <- c(specified_colors, 'black')
57 }
58 else if (!any(specified_colors=='blue'))
59 {
60 specified_colors <- c(specified_colors, 'blue')
61 }
62 else
63 {
64 specified_colors <- c(specified_colors, 'darkorange')
65 }
66 }
67 }
68 }
69 return specified_colors
70 }
71
18 get_pi_mu_var = function(components_data, number_comp) 72 get_pi_mu_var = function(components_data, number_comp)
19 { 73 {
20 # FixMe: enhance this to generically handle any integer value for number_comp. 74 # FixMe: enhance this to generically handle any integer value for number_comp.
21 if (number_comp == 1) 75 if (number_comp == 1)
22 { 76 {
56 } 110 }
57 results = c(pi, mu, var) 111 results = c(pi, mu, var)
58 return(results) 112 return(results)
59 } 113 }
60 114
61 plot_ks<-function(kaks_input, number_comp, colors, output, pi, mu, var) 115 plot_ks<-function(kaks_input, number_comp, specified_colors, output, pi, mu, var)
62 { 116 {
63 # Start PDF device driver to save charts to output. 117 # Start PDF device driver to save charts to output.
64 pdf(file=output, bg="white") 118 pdf(file=output, bg="white")
65 kaks <- read.table(file=kaks_input, header=T) 119 kaks <- read.table(file=kaks_input, header=T)
66 max_ks <- max(kaks$Ks, na.rm=TRUE) 120 max_ks <- max(kaks$Ks, na.rm=TRUE)
81 vx <- seq(1, 100) * (max_ks / 100) 135 vx <- seq(1, 100) * (max_ks / 100)
82 ymax <- max(nc) 136 ymax <- max(nc)
83 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3") 137 barplot(nc, space=0.25, offset=0, width=0.04, xlim=c(0, max_ks), ylim=c(0, ymax), col="lightpink1", border="lightpink3")
84 # Add x-axis. 138 # Add x-axis.
85 axis(1) 139 axis(1)
86 if (length(colors) == 0))
87 {
88 color <- c('red', 'yellow', 'green', 'black', 'blue', 'darkorange')
89 }
90 else
91 {
92 # Handle specified colors for components.
93 cStr <- unlist(colors)
94 color <- c()
95 items <- strsplit(cStr, ",")
96 for (item in items) {
97 color <- c(color, item)
98 }
99 num_colors_specified = length(color)
100 if (num_colors_specified < number_comp)
101 {
102 for (i in num_colors_specified:number_comp)
103 {
104 if (!any(color=='red'))
105 {
106 color <- c(color, 'red')
107 }
108 else if (!any(color=='yellow'))
109 {
110 color <- c(color, 'yellow')
111 }
112 else if (!any(color=='green'))
113 {
114 color <- c(color, 'green')
115 }
116 else if (!any(color=='black'))
117 {
118 color <- c(color, 'black')
119 }
120 else if (!any(color=='blue'))
121 {
122 color <- c(color, 'blue')
123 }
124 else
125 {
126 color <- c(color, 'darkorange')
127 }
128 }
129 }
130 }
131 for (i in 1:length(mu)) 140 for (i in 1:length(mu))
132 { 141 {
133 lines(vx, g[,i] * h, lwd=2, col=color[i]) 142 lines(vx, g[,i] * h, lwd=2, col=specified_colors[i])
134 } 143 }
135 } 144 }
136 145
137 calculate_fitted_density <- function(pi, mu, var, max_ks) 146 calculate_fitted_density <- function(pi, mu, var, max_ks)
138 { 147 {
151 } 160 }
152 } 161 }
153 return(fx) 162 return(fx)
154 } 163 }
155 164
156 # Handle colors for components.
157 if (is.na(opt$colors))
158 {
159 # Randomly specify colors for components.
160 specified_colors <- c('red', 'yellow', 'green', 'black', 'blue', 'darkorange')
161 }
162 else
163 {
164 # Handle selected colors for components.
165 parser <- newJSONParser()
166 parser$addData(opt$colors)
167 raw_colors <- parser$getObject()
168 specified_colors <- c()
169 for (raw_color in raw_colors)
170 {
171 specified_colors <- c(specified_colors, raw_color)
172 }
173 }
174
175 # Read in the components data. 165 # Read in the components data.
176 components_data <- read.delim(opt$components_input, header=TRUE) 166 components_data <- read.delim(opt$components_input, header=TRUE)
177 number_comp <- opt$number_comp 167 number_comp <- opt$number_comp
178 168
169 # Set component colors.
170 specified_colors = set_component_colors = function(opt$colors, number_comp)
171
179 # Set pi, mu, var. 172 # Set pi, mu, var.
180 items <- get_pi_mu_var(components_data, number_comp) 173 items <- get_pi_mu_var(components_data, number_comp)
181 if (number_comp == 1) 174 if (number_comp == 1)
182 { 175 {
183 pi <- items[1] 176 pi <- items[1]