2 * phylodiversitycommand.cpp
5 * Created by westcott on 4/30/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "phylodiversitycommand.h"
11 #include "treereader.h"
13 //**********************************************************************************************************************
14 vector<string> PhyloDiversityCommand::setParameters(){
17 CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree);
18 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pgroup);
19 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
20 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
21 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
22 CommandParameter pfreq("freq", "Number", "", "100", "", "", "",false,false); parameters.push_back(pfreq);
23 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
24 CommandParameter prarefy("rarefy", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(prarefy);
25 CommandParameter psummary("summary", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(psummary);
26 CommandParameter pcollect("collect", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pcollect);
27 CommandParameter pscale("scale", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pscale);
28 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
29 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
31 vector<string> myArray;
32 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
36 m->errorOut(e, "PhyloDiversityCommand", "setParameters");
40 //**********************************************************************************************************************
41 string PhyloDiversityCommand::getHelpString(){
43 string helpString = "";
44 helpString += "The phylo.diversity command parameters are tree, group, name, groups, iters, freq, processors, scale, rarefy, collect and summary. tree and group are required, unless you have valid current files.\n";
45 helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed. The group names are separated by dashes. By default all groups are used.\n";
46 helpString += "The iters parameter allows you to specify the number of randomizations to preform, by default iters=1000, if you set rarefy to true.\n";
47 helpString += "The freq parameter is used indicate when to output your data, by default it is set to 100. But you can set it to a percentage of the number of sequence. For example freq=0.10, means 10%. \n";
48 helpString += "The scale parameter is used indicate that you want your output scaled to the number of sequences sampled, default = false. \n";
49 helpString += "The rarefy parameter allows you to create a rarefaction curve. The default is false.\n";
50 helpString += "The collect parameter allows you to create a collectors curve. The default is false.\n";
51 helpString += "The summary parameter allows you to create a .summary file. The default is true.\n";
52 helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
53 helpString += "The phylo.diversity command should be in the following format: phylo.diversity(groups=yourGroups, rarefy=yourRarefy, iters=yourIters).\n";
54 helpString += "Example phylo.diversity(groups=A-B-C, rarefy=T, iters=500).\n";
55 helpString += "The phylo.diversity command output two files: .phylo.diversity and if rarefy=T, .rarefaction.\n";
56 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
60 m->errorOut(e, "PhyloDiversityCommand", "getHelpString");
66 //**********************************************************************************************************************
67 PhyloDiversityCommand::PhyloDiversityCommand(){
69 abort = true; calledHelp = true;
71 vector<string> tempOutNames;
72 outputTypes["phylodiv"] = tempOutNames;
73 outputTypes["rarefy"] = tempOutNames;
74 outputTypes["summary"] = tempOutNames;
77 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
81 //**********************************************************************************************************************
82 PhyloDiversityCommand::PhyloDiversityCommand(string option) {
84 abort = false; calledHelp = false;
86 //allow user to run help
87 if(option == "help") { help(); abort = true; calledHelp = true; }
88 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
91 vector<string> myArray = setParameters();;
93 OptionParser parser(option);
94 map<string,string> parameters = parser.getParameters();
95 map<string,string>::iterator it;
97 ValidParameters validParameter;
99 //check to make sure all parameters are valid for command
100 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
101 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
104 //initialize outputTypes
105 vector<string> tempOutNames;
106 outputTypes["phylodiv"] = tempOutNames;
107 outputTypes["rarefy"] = tempOutNames;
108 outputTypes["summary"] = tempOutNames;
110 //if the user changes the input directory command factory will send this info to us in the output parameter
111 string inputDir = validParameter.validFile(parameters, "inputdir", false);
112 if (inputDir == "not found"){ inputDir = ""; }
115 it = parameters.find("tree");
116 //user has given a template file
117 if(it != parameters.end()){
118 path = m->hasPath(it->second);
119 //if the user has not given a path then, add inputdir. else leave path alone.
120 if (path == "") { parameters["tree"] = inputDir + it->second; }
123 it = parameters.find("group");
124 //user has given a template file
125 if(it != parameters.end()){
126 path = m->hasPath(it->second);
127 //if the user has not given a path then, add inputdir. else leave path alone.
128 if (path == "") { parameters["group"] = inputDir + it->second; }
131 it = parameters.find("name");
132 //user has given a template file
133 if(it != parameters.end()){
134 path = m->hasPath(it->second);
135 //if the user has not given a path then, add inputdir. else leave path alone.
136 if (path == "") { parameters["name"] = inputDir + it->second; }
140 //check for required parameters
141 treefile = validParameter.validFile(parameters, "tree", true);
142 if (treefile == "not open") { treefile = ""; abort = true; }
143 else if (treefile == "not found") {
144 //if there is a current design file, use it
145 treefile = m->getTreeFile();
146 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
147 else { m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
148 }else { m->setTreeFile(treefile); }
150 //check for required parameters
151 groupfile = validParameter.validFile(parameters, "group", true);
152 if (groupfile == "not open") { groupfile = ""; abort = true; }
153 else if (groupfile == "not found") { groupfile = ""; }
154 else { m->setGroupFile(groupfile); }
156 namefile = validParameter.validFile(parameters, "name", true);
157 if (namefile == "not open") { namefile = ""; abort = true; }
158 else if (namefile == "not found") { namefile = ""; }
159 else { m->setNameFile(namefile); }
161 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(treefile); }
164 temp = validParameter.validFile(parameters, "freq", false); if (temp == "not found") { temp = "100"; }
165 m->mothurConvert(temp, freq);
167 temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
168 m->mothurConvert(temp, iters);
170 temp = validParameter.validFile(parameters, "rarefy", false); if (temp == "not found") { temp = "F"; }
171 rarefy = m->isTrue(temp);
172 if (!rarefy) { iters = 1; }
174 temp = validParameter.validFile(parameters, "summary", false); if (temp == "not found") { temp = "T"; }
175 summary = m->isTrue(temp);
177 temp = validParameter.validFile(parameters, "scale", false); if (temp == "not found") { temp = "F"; }
178 scale = m->isTrue(temp);
180 temp = validParameter.validFile(parameters, "collect", false); if (temp == "not found") { temp = "F"; }
181 collect = m->isTrue(temp);
183 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
184 m->setProcessors(temp);
185 m->mothurConvert(temp, processors);
187 groups = validParameter.validFile(parameters, "groups", false);
188 if (groups == "not found") { groups = ""; }
190 m->splitAtDash(groups, Groups);
191 m->setGroups(Groups);
194 if ((!collect) && (!rarefy) && (!summary)) { m->mothurOut("No outputs selected. You must set either collect, rarefy or summary to true, summary=T by default."); m->mothurOutEndLine(); abort=true; }
196 if (namefile == "") {
197 vector<string> files; files.push_back(treefile);
198 parser.getNameFile(files);
203 catch(exception& e) {
204 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
208 //**********************************************************************************************************************
210 int PhyloDiversityCommand::execute(){
213 if (abort == true) { if (calledHelp) { return 0; } return 2; }
215 m->setTreeFile(treefile);
216 TreeReader* reader = new TreeReader(treefile, groupfile, namefile);
217 vector<Tree*> trees = reader->getTrees();
218 tmap = trees[0]->getTreeMap();
222 vector<string> mGroups = m->getGroups();
223 vector<string> tGroups = tmap->getNamesOfGroups();
224 util.setGroups(mGroups, tGroups, "phylo.diversity"); //sets the groups the user wants to analyze
226 //incase the user had some mismatches between the tree and group files we don't want group xxx to be analyzed
227 for (int i = 0; i < mGroups.size(); i++) { if (mGroups[i] == "xxx") { mGroups.erase(mGroups.begin()+i); break; } }
228 m->setGroups(mGroups);
230 vector<string> outputNames;
232 //for each of the users trees
233 for(int i = 0; i < trees.size(); i++) {
235 if (m->control_pressed) { delete tmap; for (int j = 0; j < trees.size(); j++) { delete trees[j]; } for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; }
237 ofstream outSum, outRare, outCollect;
238 string outSumFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv.summary";
239 string outRareFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv.rarefaction";
240 string outCollectFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv";
242 if (summary) { m->openOutputFile(outSumFile, outSum); outputNames.push_back(outSumFile); outputTypes["summary"].push_back(outSumFile); }
243 if (rarefy) { m->openOutputFile(outRareFile, outRare); outputNames.push_back(outRareFile); outputTypes["rarefy"].push_back(outRareFile); }
244 if (collect) { m->openOutputFile(outCollectFile, outCollect); outputNames.push_back(outCollectFile); outputTypes["phylodiv"].push_back(outCollectFile); }
246 int numLeafNodes = trees[i]->getNumLeaves();
248 //create a vector containing indexes of leaf nodes, randomize it, select nodes to send to calculator
249 vector<int> randomLeaf;
250 for (int j = 0; j < numLeafNodes; j++) {
251 if (m->inUsersGroups(trees[i]->tree[j].getGroup(), mGroups) == true) { //is this a node from the group the user selected.
252 randomLeaf.push_back(j);
256 numLeafNodes = randomLeaf.size(); //reset the number of leaf nodes you are using
258 //each group, each sampling, if no rarefy iters = 1;
259 map<string, vector<float> > diversity;
261 //each group, each sampling, if no rarefy iters = 1;
262 map<string, vector<float> > sumDiversity;
264 //find largest group total
265 int largestGroup = 0;
266 for (int j = 0; j < mGroups.size(); j++) {
267 if (tmap->seqsPerGroup[mGroups[j]] > largestGroup) { largestGroup = tmap->seqsPerGroup[mGroups[j]]; }
269 //initialize diversity
270 diversity[mGroups[j]].resize(tmap->seqsPerGroup[mGroups[j]]+1, 0.0); //numSampled
273 //initialize sumDiversity
274 sumDiversity[mGroups[j]].resize(tmap->seqsPerGroup[mGroups[j]]+1, 0.0);
277 //convert freq percentage to number
279 if (freq < 1.0) { increment = largestGroup * freq;
280 }else { increment = freq; }
282 //initialize sampling spots
283 set<int> numSampledList;
284 for(int k = 1; k <= largestGroup; k++){ if((k == 1) || (k % increment == 0)){ numSampledList.insert(k); } }
285 if(largestGroup % increment != 0){ numSampledList.insert(largestGroup); }
287 //add other groups ending points
288 for (int j = 0; j < mGroups.size(); j++) {
289 if (numSampledList.count(diversity[mGroups[j]].size()-1) == 0) { numSampledList.insert(diversity[mGroups[j]].size()-1); }
292 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
294 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
297 vector<int> procIters;
299 int numItersPerProcessor = iters / processors;
301 //divide iters between processes
302 for (int h = 0; h < processors; h++) {
303 if(h == processors - 1){
304 numItersPerProcessor = iters - h * numItersPerProcessor;
306 procIters.push_back(numItersPerProcessor);
309 createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum);
311 }else{ //no need to paralellize if you dont want to rarefy
312 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
317 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
320 if (rarefy) { printData(numSampledList, sumDiversity, outRare, iters); }
324 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
326 m->mothurOutEndLine();
327 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
328 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
329 m->mothurOutEndLine();
334 catch(exception& e) {
335 m->errorOut(e, "PhyloDiversityCommand", "execute");
339 //**********************************************************************************************************************
340 int PhyloDiversityCommand::createProcesses(vector<int>& procIters, Tree* t, map< string, vector<float> >& div, map<string, vector<float> >& sumDiv, int numIters, int increment, vector<int>& randomLeaf, set<int>& numSampledList, ofstream& outCollect, ofstream& outSum){
342 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
345 vector<int> processIDS;
346 map< string, vector<float> >::iterator itSum;
348 //loop through and create all the processes you want
349 while (process != processors) {
353 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
356 driver(t, div, sumDiv, procIters[process], increment, randomLeaf, numSampledList, outCollect, outSum, false);
358 string outTemp = outputDir + toString(getpid()) + ".sumDiv.temp";
360 m->openOutputFile(outTemp, out);
362 //output the sumDIversity
363 for (itSum = sumDiv.begin(); itSum != sumDiv.end(); itSum++) {
364 out << itSum->first << '\t' << (itSum->second).size() << '\t';
365 for (int k = 0; k < (itSum->second).size(); k++) {
366 out << (itSum->second)[k] << '\t';
375 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
376 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
381 driver(t, div, sumDiv, procIters[0], increment, randomLeaf, numSampledList, outCollect, outSum, true);
383 //force parent to wait until all the processes are done
384 for (int i=0;i<(processors-1);i++) {
385 int temp = processIDS[i];
389 //get data created by processes
390 for (int i=0;i<(processors-1);i++) {
392 //input the sumDIversity
393 string inTemp = outputDir + toString(processIDS[i]) + ".sumDiv.temp";
395 m->openInputFile(inTemp, in);
397 //output the sumDIversity
398 for (int j = 0; j < sumDiv.size(); j++) {
402 in >> group >> size; m->gobble(in);
404 for (int k = 0; k < size; k++) {
408 sumDiv[group][k] += tempVal;
414 m->mothurRemove(inTemp);
422 catch(exception& e) {
423 m->errorOut(e, "PhyloDiversityCommand", "createProcesses");
427 //**********************************************************************************************************************
428 int PhyloDiversityCommand::driver(Tree* t, map< string, vector<float> >& div, map<string, vector<float> >& sumDiv, int numIters, int increment, vector<int>& randomLeaf, set<int>& numSampledList, ofstream& outCollect, ofstream& outSum, bool doSumCollect){
430 int numLeafNodes = randomLeaf.size();
431 vector<string> mGroups = m->getGroups();
433 for (int l = 0; l < numIters; l++) {
434 random_shuffle(randomLeaf.begin(), randomLeaf.end());
437 map<string, int> counts;
438 map< string, set<int> > countedBranch;
439 for (int j = 0; j < mGroups.size(); j++) { counts[mGroups[j]] = 0; countedBranch[mGroups[j]].insert(-2); } //add dummy index to initialize countedBranch sets
441 for(int k = 0; k < numLeafNodes; k++){
443 if (m->control_pressed) { return 0; }
445 //calc branch length of randomLeaf k
446 vector<float> br = calcBranchLength(t, randomLeaf[k], countedBranch);
448 //for each group in the groups update the total branch length accounting for the names file
449 vector<string> groups = t->tree[randomLeaf[k]].getGroup();
451 for (int j = 0; j < groups.size(); j++) {
452 int numSeqsInGroupJ = 0;
453 map<string, int>::iterator it;
454 it = t->tree[randomLeaf[k]].pcount.find(groups[j]);
455 if (it != t->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
456 numSeqsInGroupJ = it->second;
459 if (numSeqsInGroupJ != 0) { div[groups[j]][(counts[groups[j]]+1)] = div[groups[j]][counts[groups[j]]] + br[j]; }
461 for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
462 div[groups[j]][s] = div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths
464 counts[groups[j]] += numSeqsInGroupJ;
469 //add this diversity to the sum
470 for (int j = 0; j < mGroups.size(); j++) {
471 for (int g = 0; g < div[mGroups[j]].size(); g++) {
472 sumDiv[mGroups[j]][g] += div[mGroups[j]][g];
477 if ((collect) && (l == 0) && doSumCollect) { printData(numSampledList, div, outCollect, 1); }
478 if ((summary) && (l == 0) && doSumCollect) { printSumData(div, outSum, 1); }
484 catch(exception& e) {
485 m->errorOut(e, "PhyloDiversityCommand", "driver");
490 //**********************************************************************************************************************
492 void PhyloDiversityCommand::printSumData(map< string, vector<float> >& div, ofstream& out, int numIters){
495 out << "Groups\tnumSampled\tphyloDiversity" << endl;
497 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
499 vector<string> mGroups = m->getGroups();
500 for (int j = 0; j < mGroups.size(); j++) {
501 int numSampled = (div[mGroups[j]].size()-1);
502 out << mGroups[j] << '\t' << numSampled << '\t';
506 if (scale) { score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; }
507 else { score = div[mGroups[j]][numSampled] / (float)numIters; }
509 out << setprecision(4) << score << endl;
515 catch(exception& e) {
516 m->errorOut(e, "PhyloDiversityCommand", "printSumData");
520 //**********************************************************************************************************************
522 void PhyloDiversityCommand::printData(set<int>& num, map< string, vector<float> >& div, ofstream& out, int numIters){
525 out << "numSampled\t";
526 vector<string> mGroups = m->getGroups();
527 for (int i = 0; i < mGroups.size(); i++) { out << mGroups[i] << '\t'; }
530 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
532 for (set<int>::iterator it = num.begin(); it != num.end(); it++) {
533 int numSampled = *it;
535 out << numSampled << '\t';
537 for (int j = 0; j < mGroups.size(); j++) {
538 if (numSampled < div[mGroups[j]].size()) {
540 if (scale) { score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; }
541 else { score = div[mGroups[j]][numSampled] / (float)numIters; }
543 out << setprecision(4) << score << '\t';
544 }else { out << "NA" << '\t'; }
552 catch(exception& e) {
553 m->errorOut(e, "PhyloDiversityCommand", "printData");
557 //**********************************************************************************************************************
558 //need a vector of floats one branch length for every group the node represents.
559 vector<float> PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, map< string, set<int> >& counted){
562 //calc the branch length
563 //while you aren't at root
567 vector<string> groups = t->tree[leaf].getGroup();
568 sums.resize(groups.size(), 0.0);
570 map<string, map<int, double> > tempTotals; //maps node to total Branch Length
571 map< string, set<int> > tempCounted;
572 set<int>::iterator it;
575 if(t->tree[index].getBranchLength() != -1){
576 for (int k = 0; k < groups.size(); k++) {
577 sums[k] += abs(t->tree[index].getBranchLength());
578 counted[groups[k]].insert(index);
582 for (int k = 0; k < groups.size(); k++) {
583 tempTotals[groups[k]][index] = 0.0;
586 index = t->tree[index].getParent();
588 //while you aren't at root
589 while(t->tree[index].getParent() != -1){
591 if (m->control_pressed) { return sums; }
594 for (int k = 0; k < groups.size(); k++) {
595 map<string, int>::iterator itGroup = t->tree[index].pcount.find(groups[k]);
596 if (itGroup != t->tree[index].pcount.end()) { pcountSize++; }
598 //do both your chidren have have descendants from the users groups?
599 int lc = t->tree[index].getLChild();
600 int rc = t->tree[index].getRChild();
603 itGroup = t->tree[lc].pcount.find(groups[k]);
604 if (itGroup != t->tree[lc].pcount.end()) { LpcountSize++; }
607 itGroup = t->tree[rc].pcount.find(groups[k]);
608 if (itGroup != t->tree[rc].pcount.end()) { RpcountSize++; }
610 //if yes, add your childrens tempTotals
611 if ((LpcountSize != 0) && (RpcountSize != 0)) {
612 sums[k] += tempTotals[groups[k]][lc] + tempTotals[groups[k]][rc];
614 for (it = tempCounted[groups[k]].begin(); it != tempCounted[groups[k]].end(); it++) { counted[groups[k]].insert(*it); }
616 //cout << "added to total " << tempTotals[lc] << '\t' << tempTotals[rc] << endl;
617 if (t->tree[index].getBranchLength() != -1) {
618 if (counted[groups[k]].count(index) == 0) {
619 tempTotals[groups[k]][index] = abs(t->tree[index].getBranchLength());
620 tempCounted[groups[k]].insert(index);
622 tempTotals[groups[k]][index] = 0.0;
625 tempTotals[groups[k]][index] = 0.0;
627 }else { //if no, your tempTotal is your childrens temp totals + your branch length
628 tempTotals[groups[k]][index] = tempTotals[groups[k]][lc] + tempTotals[groups[k]][rc];
630 if (counted[groups[k]].count(index) == 0) {
631 tempTotals[groups[k]][index] += abs(t->tree[index].getBranchLength());
632 tempCounted[groups[k]].insert(index);
636 //cout << "temptotal = "<< tempTotals[i] << endl;
639 index = t->tree[index].getParent();
645 catch(exception& e) {
646 m->errorOut(e, "PhyloDiversityCommand", "calcBranchLength");
650 //**********************************************************************************************************************