]> git.donarmstrong.com Git - mothur.git/blob - kruskalwalliscommand.cpp
added kruskalwallis command
[mothur.git] / kruskalwalliscommand.cpp
1 /* 
2  * File:   kruskalwalliscommand.cpp
3  * Author: kiverson
4  *
5  * Created on June 26, 2012, 11:06 AM
6  */
7 #include "kruskalwalliscommand.h"
8
9 //**********************************************************************************************************************
10 class groupRank {
11 public:
12     string group;
13     double value;
14     double rank;
15 };
16 //**********************************************************************************************************************
17 vector<string> KruskalWallisCommand::setParameters(){   
18         try {
19                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
20                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
21                 
22                 vector<string> myArray;
23                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
24                 return myArray;
25         }
26         catch(exception& e) {
27                 m->errorOut(e, "KruskalWallisCommand", "setParameters");
28                 exit(1);
29         }
30 }
31 //**********************************************************************************************************************
32 string KruskalWallisCommand::getHelpString(){   
33         try {
34                 string helpString = "";
35                 helpString += "The kruskalwallis command parameter options are \n";
36         helpString += "Kruskal–Wallis one-way analysis of variance is a non-parametric method for testing whether samples originate from the same distribution.";
37                 return helpString;
38         }
39         catch(exception& e) {
40                 m->errorOut(e, "KruskalWallisCommand", "getHelpString");
41                 exit(1);
42         }
43 }
44 //**********************************************************************************************************************
45 string KruskalWallisCommand::getOutputFileNameTag(string type, string inputName=""){    
46         try {
47         string outputFileName = "";
48                 map<string, vector<string> >::iterator it;
49         
50         //is this a type this command creates
51         it = outputTypes.find(type);
52         if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
53         else {
54             if (type == "summary") {  outputFileName =  "cooccurence.summary"; }
55             else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
56         }
57         return outputFileName;
58         }
59         catch(exception& e) {
60                 m->errorOut(e, "KruskalWallisCommand", "getOutputFileNameTag");
61                 exit(1);
62         }
63 }
64 //**********************************************************************************************************************
65 KruskalWallisCommand::KruskalWallisCommand(){   
66         try {
67                 abort = true; calledHelp = true; 
68                 setParameters();
69         vector<string> tempOutNames;
70                 outputTypes["summary"] = tempOutNames;
71
72         }
73         catch(exception& e) {
74                 m->errorOut(e, "KruskalWallisCommand", "KruskalWallisCommand");
75                 exit(1);
76         }
77 }
78 //**********************************************************************************************************************
79 KruskalWallisCommand::KruskalWallisCommand(string option) {
80         try {
81                 abort = false; calledHelp = false;   
82                                 
83                 //allow user to run help
84                 if(option == "help") { help(); abort = true; calledHelp = true; }
85                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
86                 
87                 else {
88                         vector<string> myArray = setParameters();
89                         
90                         OptionParser parser(option);
91                         map<string,string> parameters = parser.getParameters();
92                         map<string,string>::iterator it;
93                         
94                         ValidParameters validParameter;
95                         
96                         //check to make sure all parameters are valid for command
97                         for (it = parameters.begin(); it != parameters.end(); it++) { 
98                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
99                         }
100
101                         
102                         //if the user changes the input directory command factory will send this info to us in the output parameter 
103                         string inputDir = validParameter.validFile(parameters, "inputdir", false);              
104                         if (inputDir == "not found"){   inputDir = "";          }
105                         else {
106                                 string path;
107                                 it = parameters.find("shared");
108                                 //user has given a template file
109                                 if(it != parameters.end()){ 
110                                         path = m->hasPath(it->second);
111                                         //if the user has not given a path then, add inputdir. else leave path alone.
112                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
113                                 }
114                         }
115                 
116             vector<string> tempOutNames;
117             outputTypes["summary"] = tempOutNames;
118
119
120                 }
121
122         }
123         catch(exception& e) {
124                 m->errorOut(e, "KruskalWallisCommand", "KruskalWallisCommand");
125                 exit(1);
126         }
127 }
128 //**********************************************************************************************************************
129 int KruskalWallisCommand::execute(){
130         try {
131                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
132         
133         //math goes here
134         
135         int N;
136         double ss, H;
137         double tmp = 0.0;
138                 
139         //merge all groups into a vector
140         //rank function here
141         
142         //test statistic
143         for (int i=0;i<N;i++) { tmp = tmp + (pow(rankSums[i],2) / counts[i]); }
144         
145         H = (12 / (N*(N+1))) * tmp - (3*(N+1));
146         
147         //ss = tmp - pow(accumulate(rankSums.begin(), rankSums.end(), 0), 2);
148         
149         //H = ss / ( (N * (N + 1))/12 );
150         
151         //correction for ties?
152         
153         //p-value calculation
154         
155         
156         
157                 return 0;
158         }
159         catch(exception& e) {
160                 m->errorOut(e, "KruskalWallisCommand", "execute");
161                 exit(1);
162         }
163 }
164 //**********************************************************************************************************************
165 multimap<double,double> KruskalWallisCommand::getRank(vector<groupRank> vec) {
166     try {
167         multimap<double,double> rankMap;
168         double rank = 1;
169         double previous;
170         double tie = 0.0;
171         int tiecount = 0;
172
173         sort (vec.begin(), vec.end());
174
175         for (int i=0;i<vec.size();i++) {
176             if (vec[i] != previous) { rankMap[rank] = vec[i]; }
177             else {tie = tie + rank; tiecount++;}
178             rank++;
179             previous = vec[i];
180         }
181     }
182     catch(exception& e) {
183                 m->errorOut(e, "KruskalWallisCommand", "getRank");
184                 exit(1);
185         }
186     
187 }
188 //**********************************************************************************************************************
189
190 //**********************************************************************************************************************
191 //**********************************************************************************************************************
192 //**********************************************************************************************************************