From: pschloss Date: Tue, 14 Apr 2009 14:14:38 +0000 (+0000) Subject: *** empty log message *** X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=commitdiff_plain;h=197c6d3303439582502840980d6a85cf3aab2314 *** empty log message *** --- diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj index 967e8e4..e9fbcfe 100644 --- a/Mothur.xcodeproj/project.pbxproj +++ b/Mothur.xcodeproj/project.pbxproj @@ -11,7 +11,6 @@ 372E12700F26365B0095CF7E /* readotucommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 372E126F0F26365B0095CF7E /* readotucommand.cpp */; }; 372E12960F263D5A0095CF7E /* readdistcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 372E12950F263D5A0095CF7E /* readdistcommand.cpp */; }; 372E12ED0F264D320095CF7E /* commandfactory.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 372E12EC0F264D320095CF7E /* commandfactory.cpp */; }; - 3731C1F30F8CFC0A0065A2AD /* treegroupscommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3731C1F20F8CFC0A0065A2AD /* treegroupscommand.cpp */; }; 374610780F40645300460C57 /* unifracweightedcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 374610770F40645300460C57 /* unifracweightedcommand.cpp */; }; 3746107E0F4064D100460C57 /* weighted.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3746107D0F4064D100460C57 /* weighted.cpp */; }; 374610830F40652400460C57 /* unweighted.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 374610820F40652400460C57 /* unweighted.cpp */; }; @@ -19,7 +18,6 @@ 37519A6B0F80E6EB00FED5E8 /* sharedanderbergs.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37519A6A0F80E6EB00FED5E8 /* sharedanderbergs.cpp */; }; 37519AA10F810D0200FED5E8 /* venncommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37519AA00F810D0200FED5E8 /* venncommand.cpp */; }; 37519AB50F810FAE00FED5E8 /* venn.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37519AB40F810FAE00FED5E8 /* venn.cpp */; }; - 375873E70F7D63E90040F377 /* coverage.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 375873E50F7D63E90040F377 /* coverage.cpp */; }; 375873EC0F7D64520040F377 /* fullmatrix.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 375873EB0F7D64520040F377 /* fullmatrix.cpp */; }; 375873EF0F7D646F0040F377 /* heatmap.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 375873EE0F7D646F0040F377 /* heatmap.cpp */; }; 375873F20F7D64800040F377 /* heatmapcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 375873F00F7D64800040F377 /* heatmapcommand.cpp */; }; @@ -40,7 +38,6 @@ 37AFC71F0F445386005F492D /* sharedsobscollectsummary.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37AFC71E0F445386005F492D /* sharedsobscollectsummary.cpp */; }; 37B28F680F27590100808A62 /* deconvolutecommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37B28F670F27590100808A62 /* deconvolutecommand.cpp */; }; 37C1D9730F86506E0059E3F0 /* binsequencecommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37C1D9720F86506E0059E3F0 /* binsequencecommand.cpp */; }; - 37C7F3A90F8E90AD00E91C2B /* sharedutilities.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37C7F3A80F8E90AD00E91C2B /* sharedutilities.cpp */; }; 37D928550F21331F001D4494 /* ace.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37D927B80F21331F001D4494 /* ace.cpp */; }; 37D928560F21331F001D4494 /* averagelinkage.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37D927BA0F21331F001D4494 /* averagelinkage.cpp */; }; 37D928570F21331F001D4494 /* bootstrap.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37D927BB0F21331F001D4494 /* bootstrap.cpp */; }; @@ -110,6 +107,9 @@ 37D9289F0F21331F001D4494 /* validparameter.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37D928530F21331F001D4494 /* validparameter.cpp */; }; 37E5F3E30F29FD4200F8D827 /* treenode.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37E5F3E20F29FD4200F8D827 /* treenode.cpp */; }; 37E5F4920F2A3DA800F8D827 /* readtreecommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 37E5F4910F2A3DA800F8D827 /* readtreecommand.cpp */; }; + 7E412F490F8D21B600381DD0 /* slibshuff.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 7E412F480F8D21B600381DD0 /* slibshuff.cpp */; }; + 7E412FEA0F8D3E2C00381DD0 /* libshuff.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 7E412FE90F8D3E2C00381DD0 /* libshuff.cpp */; }; + 7E4130F80F8E58FA00381DD0 /* dlibshuff.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 7E4130F60F8E58FA00381DD0 /* dlibshuff.cpp */; }; 8DD76F6A0486A84900D96B5E /* Mothur.1 in CopyFiles */ = {isa = PBXBuildFile; fileRef = C6859E8B029090EE04C91782 /* Mothur.1 */; }; A70B53AA0F4CD7AD0064797E /* getgroupcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A70B53A40F4CD7AD0064797E /* getgroupcommand.cpp */; }; A70B53AB0F4CD7AD0064797E /* getlabelcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A70B53A60F4CD7AD0064797E /* getlabelcommand.cpp */; }; @@ -144,8 +144,6 @@ 372E12940F263D5A0095CF7E /* readdistcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = readdistcommand.h; sourceTree = ""; }; 372E12950F263D5A0095CF7E /* readdistcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = readdistcommand.cpp; sourceTree = ""; }; 372E12EC0F264D320095CF7E /* commandfactory.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = commandfactory.cpp; sourceTree = ""; }; - 3731C1F10F8CFC0A0065A2AD /* treegroupscommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = treegroupscommand.h; sourceTree = ""; }; - 3731C1F20F8CFC0A0065A2AD /* treegroupscommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = treegroupscommand.cpp; sourceTree = ""; }; 374610760F40645300460C57 /* unifracweightedcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = unifracweightedcommand.h; sourceTree = ""; }; 374610770F40645300460C57 /* unifracweightedcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = unifracweightedcommand.cpp; sourceTree = ""; }; 3746107C0F4064D100460C57 /* weighted.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = weighted.h; sourceTree = ""; }; @@ -160,8 +158,6 @@ 37519AA00F810D0200FED5E8 /* venncommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = venncommand.cpp; sourceTree = ""; }; 37519AB30F810FAE00FED5E8 /* venn.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = venn.h; sourceTree = ""; }; 37519AB40F810FAE00FED5E8 /* venn.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = venn.cpp; sourceTree = ""; }; - 375873E50F7D63E90040F377 /* coverage.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = coverage.cpp; sourceTree = ""; }; - 375873E60F7D63E90040F377 /* coverage.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = coverage.h; sourceTree = ""; }; 375873EA0F7D64520040F377 /* fullmatrix.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = fullmatrix.h; sourceTree = ""; }; 375873EB0F7D64520040F377 /* fullmatrix.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = fullmatrix.cpp; sourceTree = ""; }; 375873ED0F7D646F0040F377 /* heatmap.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = heatmap.h; sourceTree = ""; }; @@ -205,8 +201,6 @@ 37B28F670F27590100808A62 /* deconvolutecommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = deconvolutecommand.cpp; sourceTree = ""; }; 37C1D9710F86506E0059E3F0 /* binsequencecommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = binsequencecommand.h; sourceTree = ""; }; 37C1D9720F86506E0059E3F0 /* binsequencecommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = binsequencecommand.cpp; sourceTree = ""; }; - 37C7F3A70F8E90AD00E91C2B /* sharedutilities.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedutilities.h; sourceTree = ""; }; - 37C7F3A80F8E90AD00E91C2B /* sharedutilities.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedutilities.cpp; sourceTree = ""; }; 37D927B80F21331F001D4494 /* ace.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = ace.cpp; sourceTree = ""; }; 37D927B90F21331F001D4494 /* ace.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = ace.h; sourceTree = ""; }; 37D927BA0F21331F001D4494 /* averagelinkage.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = averagelinkage.cpp; sourceTree = ""; }; @@ -340,6 +334,7 @@ 37D928490F21331F001D4494 /* summarydisplay.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = summarydisplay.h; sourceTree = ""; }; 37D9284A0F21331F001D4494 /* summarysharedcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = summarysharedcommand.cpp; sourceTree = ""; }; 37D9284B0F21331F001D4494 /* summarysharedcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = summarysharedcommand.h; sourceTree = ""; }; + 37D9284C0F21331F001D4494 /* utilities.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = utilities.hpp; sourceTree = ""; }; 37D9284D0F21331F001D4494 /* uvest.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = uvest.cpp; sourceTree = ""; }; 37D9284E0F21331F001D4494 /* uvest.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = uvest.h; sourceTree = ""; }; 37D9284F0F21331F001D4494 /* validcalculator.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = validcalculator.cpp; sourceTree = ""; }; @@ -353,6 +348,12 @@ 37E5F3E20F29FD4200F8D827 /* treenode.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = treenode.cpp; sourceTree = ""; }; 37E5F4900F2A3DA800F8D827 /* readtreecommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = readtreecommand.h; sourceTree = ""; }; 37E5F4910F2A3DA800F8D827 /* readtreecommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = readtreecommand.cpp; sourceTree = ""; }; + 7E412F420F8D213C00381DD0 /* libshuff.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = libshuff.h; sourceTree = ""; }; + 7E412F470F8D21B600381DD0 /* slibshuff.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = slibshuff.h; sourceTree = ""; }; + 7E412F480F8D21B600381DD0 /* slibshuff.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = slibshuff.cpp; sourceTree = ""; }; + 7E412FE90F8D3E2C00381DD0 /* libshuff.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = libshuff.cpp; sourceTree = ""; }; + 7E4130F60F8E58FA00381DD0 /* dlibshuff.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = dlibshuff.cpp; sourceTree = ""; }; + 7E4130F70F8E58FA00381DD0 /* dlibshuff.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = dlibshuff.h; sourceTree = ""; }; 8DD76F6C0486A84900D96B5E /* mothur */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = mothur; sourceTree = BUILT_PRODUCTS_DIR; }; A70B53A40F4CD7AD0064797E /* getgroupcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = getgroupcommand.cpp; sourceTree = ""; }; A70B53A50F4CD7AD0064797E /* getgroupcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = getgroupcommand.h; sourceTree = ""; }; @@ -401,10 +402,16 @@ 08FB7795FE84155DC02AAC07 /* Source */ = { isa = PBXGroup; children = ( + 7E4130F60F8E58FA00381DD0 /* dlibshuff.cpp */, + 7E4130F70F8E58FA00381DD0 /* dlibshuff.h */, 37D927BA0F21331F001D4494 /* averagelinkage.cpp */, 37D928A60F2133C0001D4494 /* calculators */, 37D927C20F21331F001D4494 /* cluster.hpp */, 37D927C10F21331F001D4494 /* cluster.cpp */, + 7E412F470F8D21B600381DD0 /* slibshuff.h */, + 7E412F480F8D21B600381DD0 /* slibshuff.cpp */, + 7E412F420F8D213C00381DD0 /* libshuff.h */, + 7E412FE90F8D3E2C00381DD0 /* libshuff.cpp */, 37D928A90F2133E5001D4494 /* commands */, 37D927C60F21331F001D4494 /* collect.h */, 37D927C50F21331F001D4494 /* collect.cpp */, @@ -454,11 +461,10 @@ 37D9281C0F21331F001D4494 /* sequence.cpp */, 37D928210F21331F001D4494 /* shared.h */, 37D928200F21331F001D4494 /* shared.cpp */, - 37C7F3A70F8E90AD00E91C2B /* sharedutilities.h */, - 37C7F3A80F8E90AD00E91C2B /* sharedutilities.cpp */, 37D928420F21331F001D4494 /* singlelinkage.cpp */, 37D928480F21331F001D4494 /* summarydata.h */, 37D928490F21331F001D4494 /* summarydisplay.h */, + 37D9284C0F21331F001D4494 /* utilities.hpp */, 37519AB30F810FAE00FED5E8 /* venn.h */, 37519AB40F810FAE00FED5E8 /* venn.cpp */, ); @@ -489,8 +495,6 @@ EB1216E40F61ACFB004A865F /* bstick.cpp */, 37D927C00F21331F001D4494 /* chao1.h */, 37D927BF0F21331F001D4494 /* chao1.cpp */, - 375873E60F7D63E90040F377 /* coverage.h */, - 375873E50F7D63E90040F377 /* coverage.cpp */, EB9303F70F53517300E8EF26 /* geom.h */, EB9303F80F53517300E8EF26 /* geom.cpp */, EB9303E90F534D9400E8EF26 /* logsd.h */, @@ -616,8 +620,6 @@ 37D928460F21331F001D4494 /* summarycommand.cpp */, 37D9284B0F21331F001D4494 /* summarysharedcommand.h */, 37D9284A0F21331F001D4494 /* summarysharedcommand.cpp */, - 3731C1F10F8CFC0A0065A2AD /* treegroupscommand.h */, - 3731C1F20F8CFC0A0065A2AD /* treegroupscommand.cpp */, 3746109B0F40657600460C57 /* unifracunweightedcommand.h */, 3746109C0F40657600460C57 /* unifracunweightedcommand.cpp */, 374610760F40645300460C57 /* unifracweightedcommand.h */, @@ -826,7 +828,6 @@ EB1216E50F61ACFB004A865F /* bstick.cpp in Sources */, EB1217230F61C9AC004A865F /* sharedkstest.cpp in Sources */, EB6F015B0F6AC1670048081A /* sharedbdiversity.cpp in Sources */, - 375873E70F7D63E90040F377 /* coverage.cpp in Sources */, 375873EC0F7D64520040F377 /* fullmatrix.cpp in Sources */, 375873EF0F7D646F0040F377 /* heatmap.cpp in Sources */, 375873F20F7D64800040F377 /* heatmapcommand.cpp in Sources */, @@ -843,8 +844,9 @@ 37519AB50F810FAE00FED5E8 /* venn.cpp in Sources */, 37C1D9730F86506E0059E3F0 /* binsequencecommand.cpp in Sources */, 370B88070F8A4EE4005AB382 /* getoturepcommand.cpp in Sources */, - 3731C1F30F8CFC0A0065A2AD /* treegroupscommand.cpp in Sources */, - 37C7F3A90F8E90AD00E91C2B /* sharedutilities.cpp in Sources */, + 7E412F490F8D21B600381DD0 /* slibshuff.cpp in Sources */, + 7E412FEA0F8D3E2C00381DD0 /* libshuff.cpp in Sources */, + 7E4130F80F8E58FA00381DD0 /* dlibshuff.cpp in Sources */, ); runOnlyForDeploymentPostprocessing = 0; }; @@ -897,7 +899,7 @@ ppc, i386, ); - GCC_OPTIMIZATION_LEVEL = 0; + GCC_OPTIMIZATION_LEVEL = 3; GCC_WARN_ABOUT_RETURN_TYPE = YES; GCC_WARN_UNUSED_VARIABLE = YES; PREBINDING = NO; diff --git a/coverage.cpp b/coverage.cpp deleted file mode 100644 index dea8dff..0000000 --- a/coverage.cpp +++ /dev/null @@ -1,185 +0,0 @@ -/* - * coverage.cpp - * Mothur - * - * Created by Sarah Westcott on 3/9/09. - * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved. - * - */ - -/* This class library coverage at the given distances of the Cramer-von Mises statistic. - you may refer to the "Integration of Microbial Ecology and Statistics: A Test To Compare Gene Libraries" - paper in Applied and Environmental Microbiology, Sept. 2004, p. 5485-5492 0099-2240/04/$8.00+0 - DOI: 10.1128/AEM.70.9.5485-5492.2004 Copyright 2004 American Society for Microbiology for more information. */ - - -#include "coverage.h" - -//********************************************************************************************************************** -Coverage::Coverage() { - globaldata = GlobalData::getInstance(); - numUserGroups = globaldata->Groups.size(); - numGroups = globaldata->gGroupmap->getNumGroups(); -} - -//********************************************************************************************************************** -void Coverage::getValues(FullMatrix* matrix, vector< vector< vector > >& data, vector dist) { - try { - vector min; - vector groups; - - //initialize data - data.resize(dist.size()); - for (int l = 0; l < data.size(); l++) { - data[l].resize(numGroups); - for (int k = 0; k < data[l].size(); k++) { - data[l][k].push_back(0.0); - } - } - - /**************************************/ - //get the minimums for each comparision - /**************************************/ - int count = 0; - - for (int i = 0; i < numGroups; i++) { - for (int j = 0; j < numGroups; j++) { - - //is this "box" one the user wants analyzed? - if ((inUsersGroups(globaldata->gGroupmap->namesOfGroups[i], globaldata->Groups) == true) && (inUsersGroups(globaldata->gGroupmap->namesOfGroups[j], globaldata->Groups) == true)) { - - min = matrix->getMins(count); //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB; - - //find the coverage at this distance - sort(min.begin(), min.end()); - - //loop through each distance and fill data - for (int k = 0; k < data.size(); k++) { - - int index = -1; - //find index in min where value is higher than d - for (int m = 0; m < min.size(); m++) { - if (min[m] > dist[k]) { index = m; break; } - } - - // if you don't find one than all the mins are less than d - if (index == -1) { index = min.size(); } - - //save value in data - data[k][i][j] = 1.0 - ((min.size()-index)/(float)min.size()); - - } - } - count++; - } - } - - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} - -//********************************************************************************************************************** -void Coverage::getValues(FullMatrix* matrix, vector< vector< vector > >& data, vector dist, string mode) { - try { - vector min1; - vector min2; - vector groups; - - //initialize data - data.resize(dist.size()); - for (int l = 0; l < data.size(); l++) { - data[l].resize(numGroups); - for (int k = 0; k < data[l].size(); k++) { - data[l][k].push_back(0.0); - } - } - - int count = 0; - int count2 = 0; - - //for each box - for (int i = 0; i < numGroups; i++) { - for (int j = 0; j < numGroups; j++) { - - if (i != j) { - //is this "box" one the user wants analyzed? - if ((inUsersGroups(globaldata->gGroupmap->namesOfGroups[i], globaldata->Groups) == true) && (inUsersGroups(globaldata->gGroupmap->namesOfGroups[j], globaldata->Groups) == true)) { - - matrix->shuffle(globaldata->gGroupmap->namesOfGroups[i], globaldata->gGroupmap->namesOfGroups[j]); - - min1 = matrix->getMins(count); //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB; - min2 = matrix->getMins(count2); //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB; - - //find the coverage at this distance - sort(min1.begin(), min1.end()); - - //find the coverage at this distance - sort(min2.begin(), min2.end()); - - float distDiff = 0; - - //loop through each distance and fill data - for (int k = 0; k < data.size(); k++) { - //****** coverage of AA **********// - int index = -1; - //find index in min where value is higher than d - for (int m = 0; m < min1.size(); m++) { - if (min1[m] > dist[k]) { index = m; break; } - } - - // if you don't find one than all the mins are less than d - if (index == -1) { index = min1.size(); } - - //****** coverage of AB **********// - int index2 = -1; - //find index in min where value is higher than d - for (int m = 0; m < min2.size(); m++) { - if (min2[m] > dist[k]) { index2 = m; break; } - } - - // if you don't find one than all the mins are less than d - if (index2 == -1) { index2 = min2.size(); } - - //coverage of ii - float covII = 1.0 - ((min1.size()-index)/(float)min1.size()); - - //coverage of ij - float covIJ = 1.0 - ((min2.size()-index2)/(float)min2.size()); - - //save value in data (Caa - Cab)^2 * distDiff - data[k][i][j] = ((covII-covIJ) * (covII-covIJ)) * distDiff; - - //update distDiff - if (k < data.size() - 1) { - distDiff = dist[k+1] - dist[k]; - } - } - - //put matrix back to original - matrix->restore(); - } - } - count2++; - } - count += numGroups+1; //go from AA to BB to CC - } - - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - -} - diff --git a/coverage.h b/coverage.h deleted file mode 100644 index b4dccf7..0000000 --- a/coverage.h +++ /dev/null @@ -1,40 +0,0 @@ -#ifndef COVERAGE_H -#define COVERAGE_H - -/* - * coverage.h - * Mothur - * - * Created by Sarah Westcott on 3/9/09. - * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved. - * - */ - -#include "mothur.h" -#include "fullmatrix.h" -#include "globaldata.hpp" - -using namespace std; - -/***********************************************************************/ - -class Coverage { - - public: - Coverage(); - ~Coverage(){}; - void getValues(FullMatrix*, vector< vector< vector > >&, vector, string); //matrix, container for results, vector of distances, mode - for random matrices - void getValues(FullMatrix*, vector< vector< vector > >&, vector); //for user matrix - - private: - GlobalData* globaldata; - int numGroups, numUserGroups; - -}; - - -/***********************************************************************/ - - - -#endif \ No newline at end of file diff --git a/dlibshuff.cpp b/dlibshuff.cpp new file mode 100644 index 0000000..9db8872 --- /dev/null +++ b/dlibshuff.cpp @@ -0,0 +1,85 @@ +/* + * DLibshuff.cpp + * Mothur + * + * Created by Pat Schloss on 4/8/09. + * Copyright 2009 Patrick D. Schloss. All rights reserved. + * + */ + +#include "dlibshuff.h" + +/***********************************************************************/ + +DLibshuff::DLibshuff(FullMatrix* D, int it, float step, float co) : Libshuff(D, it, step, co){ + + numDXs = int(cutOff / stepSize); + +} + +/***********************************************************************/ + +float DLibshuff::evaluatePair(int i, int j){ + return dCalculate(i,j); +} + +/***********************************************************************/ + +vector > DLibshuff::evaluateAll(){ + savedMins.resize(numGroups); + vector > dCXYValues(numGroups); + + for(int i=0;i nx = calcN(minX); + vector nxy = calcN(minXY); + + double sum = 0; + + for(int i=0;i DLibshuff::calcN(vector minVector){ + + vector counts(numDXs,0); + + int precision = int(1 / stepSize); + + for(int i=0;i > evaluateAll(); + float evaluatePair(int, int); + +private: + int numDXs; + double dCalculate(int, int); + vector calcN(vector); +}; + +#endif diff --git a/fullmatrix.cpp b/fullmatrix.cpp index f4acd78..bc53b7d 100644 --- a/fullmatrix.cpp +++ b/fullmatrix.cpp @@ -10,6 +10,7 @@ #include "fullmatrix.h" /**************************************************************************/ + //This constructor reads a distance matrix file and stores the data in the matrix. FullMatrix::FullMatrix(ifstream& filehandle) { try{ @@ -27,12 +28,14 @@ FullMatrix::FullMatrix(ifstream& filehandle) { group = groupmap->getGroup(name); if(group == "not found") { cout << "Error: Sequence '" << name << "' was not found in the group file, please correct." << endl; exit(1); } - index[0].groupname = group; + index.resize(numSeqs); + index[0].groupName = group; index[0].seqName = name; //determine if matrix is square or lower triangle //if it is square read the distances for the first sequence char d; + bool square; while((d=filehandle.get()) != EOF){ //is d a number meaning its square @@ -59,7 +62,20 @@ FullMatrix::FullMatrix(ifstream& filehandle) { //sort sequences so they are gathered in groups for processing sortGroups(0, numSeqs-1); - + + groups.push_back(index[0].groupName); + sizes.push_back(1); + int groupCount = 0; + + for(int i=1;i> name; group = groupmap->getGroup(name); - index[i].groupname = group; + index[i].groupName = group; index[i].seqName = name; if(group == "not found") { cout << "Error: Sequence '" << name << "' was not found in the group file, please correct." << endl; exit(1); } @@ -125,7 +141,7 @@ void FullMatrix::readLTMatrix(ifstream& filehandle) { filehandle >> name; group = groupmap->getGroup(name); - index[i].groupname = group; + index[i].groupName = group; index[i].seqName = name; if(group == "not found") { cout << "Error: Sequence '" << name << "' was not found in the group file, please correct." << endl; exit(1); } @@ -154,6 +170,7 @@ void FullMatrix::readLTMatrix(ifstream& filehandle) { } /**************************************************************************/ + void FullMatrix::sortGroups(int low, int high){ try{ @@ -164,15 +181,15 @@ void FullMatrix::sortGroups(int low, int high){ /* compare value */ //what group does this row belong to - string z = index[(low + high) / 2].groupname; + string z = index[(low + high) / 2].groupName; /* partition */ do { /* find member above ... */ - while(index[i].groupname < z) i++; + while(index[i].groupName < z) i++; /* find element below ... */ - while(index[j].groupname > z) j--; + while(index[j].groupName > z) j--; if(i <= j) { /* swap rows*/ @@ -190,9 +207,9 @@ void FullMatrix::sortGroups(int low, int high){ } //swap map elements - z = index[i].groupname; - index[i].groupname = index[j].groupname; - index[j].groupname = z; + z = index[i].groupName; + index[i].groupName = index[j].groupName; + index[j].groupName = z; name = index[i].seqName; index[i].seqName = index[j].seqName; @@ -225,345 +242,47 @@ void FullMatrix::sortGroups(int low, int high){ } /**************************************************************************/ -int FullMatrix::getNumSeqs(){ return numSeqs; } -/**************************************************************************/ -//print out matrix -void FullMatrix::printMatrix(ostream& out) { - try{ - for (int i = 0; i < numSeqs; i++) { - out << "row " << i << " group = " << index[i].groupname << " name = " << index[i].seqName << endl; - for (int j = 0; j < numSeqs; j++) { - out << matrix[i][j] << " "; - } - out << endl; - } - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function printMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the FullMatrix class function printMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - -} -/**************************************************************************/ -void FullMatrix::setBounds(){ - try{ - numGroups = globaldata->gGroupmap->namesOfGroups.size(); - - //sort globaldata->gGroupmap.namesOfGroups so that it will match the matrix - sort(globaldata->gGroupmap->namesOfGroups.begin(), globaldata->gGroupmap->namesOfGroups.end()); - - //one for each comparision - //minsForRows.resize(numGroups*numGroups); - - /*************************************************/ - //find where in matrix each group starts and stops - /*************************************************/ - bounds.resize(numGroups); - - bounds[0] = 0; - bounds[numGroups] = numSeqs; +float FullMatrix::get(int i, int j){ return matrix[i][j]; } - //for each group find bounds of subgroup/comparison - for (int i = 1; i < numGroups; i++) { - getBounds(bounds[i], globaldata->gGroupmap->namesOfGroups[i-1]); - } - - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function getMinsForRowsVectors. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the FullMatrix class function getMinsForRowsVectors. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } +/**************************************************************************/ -} -/**************************************************************************/ -vector FullMatrix::getMins(int x) { - try{ - //clear out old data - minsForRows.clear(); - - /************************************************************/ - //fill the minsForRows vector for the box the user wants - /************************************************************/ - int count = 0; - int lowBoundx = bounds[0]; //where first group starts - int lowBoundy = bounds[0]; - int highBoundx = bounds[1]; //where second group starts - int highBoundy = bounds[1]; - - int countx = 1; //index in bound - int county = 1; //index in bound - - //find the bounds for the box the user wants - for (int i = 0; i < (numGroups * numGroups); i++) { - - //are you at the box? - if (count == x) { break; } - else { count++; } - - //move to next box - if (county < numGroups) { - county++; - highBoundy = bounds[county]; - lowBoundy = bounds[county-1]; - }else{ //you are moving to a new row of "boxes" - county = 1; - countx++; - highBoundx = bounds[countx]; - lowBoundx = bounds[countx-1]; - highBoundy = bounds[county]; - lowBoundy = bounds[county-1]; - } - } - - //each row in the box - for (int x = lowBoundx; x < highBoundx; x++) { - float min4Row = 100000.0; - //each entry in that row - for (int y = lowBoundy; y < highBoundy; y++) { - //if you are not on the diagonal and you are less than previous minimum - if ((x != y) && (matrix[x][y] < min4Row)) { - min4Row = matrix[x][y]; - } - } - //save minimum value for that row in minsForRows vector of vectors - minsForRows.push_back(min4Row); - } - - return minsForRows; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function getMins. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the FullMatrix class function getMins. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/**************************************************************************/ -void FullMatrix::getBounds(int& higher, string group) { - try{ - bool gotLower = false; - - //for each group find bounds of subgroup/comparison - for (it = index.begin(); it != index.end(); it++) { - if (it->second.groupname == group) { - gotLower = true; - }else if ((gotLower == true) && (it->second.groupname != group)) { higher = it->first; break; } - } - - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function getBounds. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the FullMatrix class function getBounds. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } +vector FullMatrix::getGroups(){ return groups; } -} +/**************************************************************************/ -/**************************************************************************/ -//print out matrix -void FullMatrix::printMinsForRows(ostream& out) { - try{ - for (int j = 0; j < minsForRows.size(); j++) { - out << minsForRows[j] << " "; - } - out << endl; +vector FullMatrix::getSizes(){ return sizes; } - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function printMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the FullMatrix class function printMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } +/**************************************************************************/ -} +int FullMatrix::getNumGroups(){ return groups.size(); } -/**************************************************************************/ -//shuffles the sequences in the 2 groups passed in. -void FullMatrix::shuffle(string groupA, string groupB){ - try{ - vector rows2Swap; - vector shuffled; - float y = 0; - string name = ""; - - - /********************************/ - //save rows you want to randomize - /********************************/ - //go through the matrix map to find the rows from groups you want to randomize - for (it = index.begin(); it != index.end(); it++) { - //is this row from group A or B? - if ((it->second.groupname == groupA) || (it->second.groupname == groupB)) { - rows2Swap.push_back(it->first); - shuffled.push_back(it->first); - } - } - - //randomize rows to shuffle in shuffled - random_shuffle(shuffled.begin(), shuffled.end()); - - /***************************************/ - //swap rows and columns to randomize box - /***************************************/ - for (int i = 0; i < shuffled.size(); i++) { +/**************************************************************************/ - //record the swaps you are making so you can undo them in restore function - restoreIndex[i].a = shuffled[i]; - restoreIndex[i].b = rows2Swap[i]; - - /* swap rows*/ - for (int h = 0; h < numSeqs; h++) { - y = matrix[shuffled[i]][h]; - matrix[shuffled[i]][h] = matrix[rows2Swap[i]][h]; - matrix[rows2Swap[i]][h] = y; - } - - /* swap columns */ - for (int b = 0; b < numSeqs; b++) { - y = matrix[b][shuffled[i]]; - matrix[b][shuffled[i]] = matrix[b][rows2Swap[i]]; - matrix[b][rows2Swap[i]] = y; - } - - //swap map elements - name = index[shuffled[i]].seqName; - index[shuffled[i]].seqName = index[rows2Swap[i]].seqName; - index[rows2Swap[i]].seqName = name; +int FullMatrix::getNumSeqs(){ return numSeqs; } - } - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function shuffle. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the FullMatrix class function shuffle. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} /**************************************************************************/ -//unshuffles the matrix. -void FullMatrix::restore(){ - try{ - float y = 0; - string name = ""; - - //reverse iterate through swaps and undo them to restore original matrix and index map. - for(it2 = restoreIndex.rbegin(); it2 != restoreIndex.rend(); it2++) { - /* swap rows */ - - for (int h = 0; h < numSeqs; h++) { - y = matrix[it2->second.a][h]; - matrix[it2->second.a][h] = matrix[it2->second.b][h]; - matrix[it2->second.b][h] = y; - } - - /* swap columns */ - for (int b = 0; b < numSeqs; b++) { - y = matrix[b][it2->second.a]; - matrix[b][it2->second.a] = matrix[b][it2->second.b]; - matrix[b][it2->second.b] = y; - } - - - //swap map elements - name = index[it2->second.a].seqName; - index[it2->second.a].seqName = index[it2->second.b].seqName; - index[it2->second.b].seqName = name; - } - - //clear restore for next shuffle - restoreIndex.clear(); - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function restore. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the FullMatrix class function restore. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/**************************************************************************/ -void FullMatrix::getDist(vector& distances) { +void FullMatrix::printMatrix(ostream& out) { try{ - map dist; //holds the distances for the integral form - map::iterator it; - - /************************************************************/ - //fill the minsForRows vectors for each group the user wants - /************************************************************/ - int lowBoundx = bounds[0]; //where first group starts - int lowBoundy = bounds[0]; - int highBoundx = bounds[1]; //where second group starts - int highBoundy = bounds[1]; - - int countx = 1; //index in bound - int county = 1; //index in bound - - //go through each "box" in the matrix - for (int i = 0; i < (numGroups * numGroups); i++) { - //each row in the box - for (int x = lowBoundx; x < highBoundx; x++) { - float min4Row = 100000.0; - //each entry in that row - for (int y = lowBoundy; y < highBoundy; y++) { - //if you are not on the diagonal and you are less than previous minimum - if ((x != y) && (matrix[x][y] < min4Row)){ - min4Row = matrix[x][y]; - } - } - //save minimum value - dist[min4Row] = min4Row; - } - - //****** reset bounds to process next "box" ******** - //if you still have more "boxes" in that row - if (county < numGroups) { - county++; - highBoundy = bounds[county]; - lowBoundy = bounds[county-1]; - }else{ //you are moving to a new row of "boxes" - county = 1; - countx++; - highBoundx = bounds[countx]; - lowBoundx = bounds[countx-1]; - highBoundy = bounds[county]; - lowBoundy = bounds[county-1]; + for (int i = 0; i < numSeqs; i++) { + out << "row " << i << " group = " << index[i].groupName << " name = " << index[i].seqName << endl; + for (int j = 0; j < numSeqs; j++) { + out << matrix[i][j] << " "; } + out << endl; } - - //store distances in users vector - for (it = dist.begin(); it != dist.end(); it++) { - distances.push_back(it->first); - } - } catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function restore. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function printMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; exit(1); } catch(...) { - cout << "An unknown error has occurred in the FullMatrix class function restore. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + cout << "An unknown error has occurred in the FullMatrix class function printMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; exit(1); } + } +/**************************************************************************/ + diff --git a/fullmatrix.h b/fullmatrix.h index 078cdb7..5912b6c 100644 --- a/fullmatrix.h +++ b/fullmatrix.h @@ -17,52 +17,39 @@ using namespace std; struct Names { - string groupname; string seqName; + string groupName; }; -struct Swap { - int a; - int b; -}; - - class FullMatrix { - public: - FullMatrix(){}; - FullMatrix(ifstream&); - ~FullMatrix(){}; +public: + FullMatrix(){}; + FullMatrix(ifstream&); + ~FullMatrix(){}; - int getNumSeqs(); - void printMatrix(ostream&); - void setBounds(); //requires globaldata->Groups to be filled - vector getMins(int); //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB; - void getDist(vector&); //fills a vector with the valid distances for the integral form. - void shuffle(string, string); //shuffles the sequences in the groups passed in. - void restore(); //unshuffles the matrix. + int getNumSeqs(); + vector getSizes(); + vector getGroups(); + int getNumGroups(); + void printMatrix(ostream&); + float get(int, int); - private: - void sortGroups(int, int); //this function sorts the sequences within the matrix. - void getBounds(int&, string); - void readSquareMatrix(ifstream&); - void readLTMatrix(ifstream&); - void printMinsForRows(ostream&); - - map index; // row in vector, sequence group. need to know this so when we sort it can be updated. - map restoreIndex; //a map of the swaps made so you can undo them in restore. - map::iterator it; - map::reverse_iterator it2; - - vector< vector > matrix; //a 2D distance matrix of all the sequences and their distances to eachother. - vector minsForRows; //vector< minimum distance for that subrow> - one for each comparison. - vector bounds; //bounds[1] = starting row in matrix from group B, bounds[2] = starting row in matrix from group C, bounds[3] = no need to find upper bound of C because its numSeqs. - - GroupMap* groupmap; //maps sequences to groups they belong to. - GlobalData* globaldata; - int numSeqs, numGroups, numUserGroups; - bool square; - +private: + vector< vector > matrix; //a 2D distance matrix of all the sequences and their distances to eachother. + void readSquareMatrix(ifstream&); + void readLTMatrix(ifstream&); + vector index; // row in vector, sequence group. need to know this so when we sort it can be updated. + vector sizes; + vector groups; + void sortGroups(int, int); //this function sorts the sequences within the matrix. + + + GroupMap* groupmap; //maps sequences to groups they belong to. + int numSeqs; + int numGroups; + int numUserGroups; + GlobalData* globaldata; }; #endif \ No newline at end of file diff --git a/libshuff.cpp b/libshuff.cpp new file mode 100644 index 0000000..b7e5f27 --- /dev/null +++ b/libshuff.cpp @@ -0,0 +1,163 @@ +/* + * libshuffform.cpp + * Mothur + * + * Created by Pat Schloss on 4/8/09. + * Copyright 2009 Patrick D. Schloss. All rights reserved. + * + */ + +#include "libshuff.h" + +/***********************************************************************/ + +void swap(int& i,int& j){ int t = i; i = j; j = t; } + +/***********************************************************************/ + +Libshuff::Libshuff(FullMatrix* D, int it, float step, float co) : matrix(D), iters(it), stepSize(step), cutOff(co){ + try{ + groupNames = matrix->getGroups(); + groupSizes = matrix->getSizes(); + numGroups = matrix->getNumGroups(); + + initializeGroups(matrix); + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the Libshuff class Function Libshuff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Libshuff class function Libshuff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + +} + +/***********************************************************************/ + +void Libshuff::initializeGroups(FullMatrix* matrix){ + try{ + groups.resize(numGroups); + savedGroups.resize(numGroups); + + savedGroups.resize(numGroups); + for(int i=0;i > > Libshuff::getSavedMins(){ + return savedMins; +} + +/***********************************************************************/ + +vector Libshuff::getMinX(int x){ + try{ + vector minX(groupSizes[x], 0); + for(int i=0;i 1 ? (i==0 ? matrix->get(groups[x][0], groups[x][1]) : matrix->get(groups[x][i], groups[x][0])) : 0.0); + for(int j=0;jget(groups[x][i], groups[x][j]); + if(dx < minX[i]){ minX[i] = dx; } + } + } + } + return minX; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the Libshuff class Function getMinX. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Libshuff class function getMinX. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + +} + +/***********************************************************************/ + +vector Libshuff::getMinXY(int x, int y){ + try{ + vector minXY(groupSizes[x], 0); + + for(int i=0;iget(groups[x][i], groups[y][0]); + for(int j=0;jget(groups[x][i], groups[y][j]); + if(dxy v(nv); + + int index=0; + for(int k=0;k0;k--){ + int z = (int)(rand() % k); + swap(v[z],v[k]); + } + + index=0; + for(int k=0;k > evaluateAll() = 0; + virtual float evaluatePair(int,int) = 0; + void randomizeGroups(int, int); + void resetGroup(int); + vector > > getSavedMins(); + +protected: + void initializeGroups(FullMatrix*); + vector getMinX(int); + vector getMinXY(int, int); + + vector > > savedMins; + + + FullMatrix* matrix; + vector groupSizes; + vector groupNames; + vector > groups; + vector > savedGroups; + vector minX; + vector minXY; + float cutOff; + int iters; + float stepSize; + + int numGroups; +}; + +#endif diff --git a/libshuffcommand.cpp b/libshuffcommand.cpp index 33a1a3b..1e38226 100644 --- a/libshuffcommand.cpp +++ b/libshuffcommand.cpp @@ -14,46 +14,29 @@ #include "libshuffcommand.h" +#include "libshuff.h" +#include "slibshuff.h" +#include "dlibshuff.h" //********************************************************************************************************************** - LibShuffCommand::LibShuffCommand(){ try { - globaldata = GlobalData::getInstance(); - convert(globaldata->getCutOff(), cutOff); - convert(globaldata->getIters(), iters); - convert(globaldata->getStep(), step); - form = globaldata->getForm(); - matrix = globaldata->gMatrix; - coverageFile = getRootName(globaldata->getPhylipFile()) + "coverage"; - summaryFile = getRootName(globaldata->getPhylipFile()) + "slsummary"; - openOutputFile(coverageFile, out); - openOutputFile(summaryFile, outSum); + srand( (unsigned)time( NULL ) ); - //set the groups to be analyzed - setGroups(); + globaldata = GlobalData::getInstance(); + convert(globaldata->getCutOff(), cutOff); //get the cutoff + convert(globaldata->getIters(), iters); //get the number of iterations + convert(globaldata->getStep(), step); //get the step size for the discrete command + matrix = globaldata->gMatrix; //get the distance matrix + setGroups(); //set the groups to be analyzed - //file headers for coverage file - out << "D" << '\t'; - for (int i = 0; i < groupComb.size(); i++) { - out << "C" + groupComb[i] << '\t'; + if(globaldata->getForm() == "discrete"){ + form = new DLibshuff(matrix, iters, step, cutOff); } - - for (int i = 0; i < numGroups; i++) { - for (int j = 0; j < numGroups; j++) { - //don't output AA to AA - if (i != j) { - out << "Delta" + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t'; - } - } + else{ + form = new SLibshuff(matrix, iters, cutOff); } - out << endl; - - numComp = numGroups*numGroups; - - coverage = new Coverage(); - } catch(exception& e) { cout << "Standard Error: " << e.what() << " has occurred in the LibShuffCommand class Function LibShuffCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; @@ -68,134 +51,42 @@ LibShuffCommand::LibShuffCommand(){ //********************************************************************************************************************** -LibShuffCommand::~LibShuffCommand(){ - delete coverage; -} - -//********************************************************************************************************************** - int LibShuffCommand::execute(){ try { - //deltaValues[0] = scores for the difference between AA and AB. - //cValues[0][0][0] = AA at distance 0.0, cValues[0][0][1] = AB at distance 0.0, cValues[0][0][2] = AC at distance 0.0, cValues[0][1][0] = BA at distance 0.0, cValues[0][1][1] = BB... - Progress* reading; - reading = new Progress("Comparing to random:", iters); - - sumDelta.resize(numComp-numGroups, 0.0); - - matrix->setBounds(); - - //load distances - if (form != "discrete") { matrix->getDist(dist); } - else { - float f = 0.0; - while (f <= cutOff) { - dist.push_back(f); - f += step; - } - } - - /*****************************/ - //get values for users matrix - /*****************************/ - - //clear out old Values - deltaValues.clear(); - deltaValues.resize(dist.size()); - - coverage->getValues(matrix, cValues, dist); - - float distDiff = 0; - - //loop through each distance and load sumdelta - for (int p = 0; p < cValues.size(); p++) { - //find delta values - int count = 0; - for (int i = 0; i < numGroups; i++) { - for (int j = 0; j < numGroups; j++) { - //don't save AA to AA - if (i != j) { - //(Caa - Cab)^2 * distDiff - deltaValues[p].push_back(((cValues[p][i][i]-cValues[p][i][j]) * (cValues[p][i][i]-cValues[p][i][j])) * distDiff); //* distDiff - sumDelta[count] += deltaValues[p][count]; - count++; - } - } - } - if (p < cValues.size() - 1) { - distDiff = dist[p+1] - dist[p]; - } - } - - printCoverageFile(); - - /*******************************************************************************/ - //create and score random matrixes finding the sumDelta values for summary file - /******************************************************************************/ - //initialize rsumDelta - rsumDelta.resize(numComp-numGroups); - for (int l = 0; l < rsumDelta.size(); l++) { - for (int w = 0; w < iters; w++) { - rsumDelta[l].push_back(0.0); - } - } + savedDXYValues = form->evaluateAll(); + savedMinValues = form->getSavedMins(); + pValueCounts.assign(numGroups, 0); + for(int i=0;igetValues(matrix, cValues, dist, "random"); - - //loop through each distance and load rsumdelta - for (int p = 0; p < cValues.size(); p++) { - //find delta values - int count = 0; - for (int i = 0; i < numGroups; i++) { - for (int j = 0; j < numGroups; j++) { - //don't save AA to AA - if (i != j) { - //rsumDelta[3][500] = the sum of the delta scores for BB-BC for random matrix # 500. - rsumDelta[count][m] += cValues[p][i][j]; // where cValues[p][0][1] = delta value at distance p of AA-AB, cValues[p][1][2] = delta value at distance p of BB-BC. - count++; - } - } + for(int i=0;inewLine(groupNames[i]+'-'+groupNames[j], iters); + for(int p=0;prandomizeGroups(i,j); + if(form->evaluatePair(i,j) >= savedDXYValues[i][j]) { pValueCounts[i][j]++; } + if(form->evaluatePair(j,i) >= savedDXYValues[j][i]) { pValueCounts[j][i]++; } + reading->update(p); } + form->resetGroup(i); + form->resetGroup(j); } - - //clear out old Values - reading->update(m); - cValues.clear(); - } - reading->finish(); delete reading; - - /**********************************************************/ - //find the signifigance of the user matrix' sumdelta values - /**********************************************************/ - - for (int t = 0; t < rsumDelta.size(); t++) { - //sort rsumDelta t - sort(rsumDelta[t].begin(), rsumDelta[t].end()); - - //the index of the score higher than yours is returned - //so if you have 1000 random matrices the index returned is 100 - //then there are 900 matrices with a score greater then you. - //giving you a signifigance of 0.900 - int index = findIndex(sumDelta[t], t); - - //the signifigance is the number of trees with the users score or higher - sumDeltaSig.push_back((iters-index)/(float)iters); - } - + cout << endl; printSummaryFile(); + printCoverageFile(); //clear out users groups globaldata->Groups.clear(); + delete form; return 0; } @@ -208,27 +99,81 @@ int LibShuffCommand::execute(){ exit(1); } } + //********************************************************************************************************************** + void LibShuffCommand::printCoverageFile() { try { - //format output - out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); + + ofstream outCov; + summaryFile = getRootName(globaldata->getPhylipFile()) + "libshuff.coverage"; + openOutputFile(summaryFile, outCov); + outCov.setf(ios::fixed, ios::floatfield); outCov.setf(ios::showpoint); + cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint); + + map > allDistances; + map >::iterator it; + + vector > indices(numGroups); + int numIndices = numGroups * numGroups; - //loop through each distance - for (int p = 0; p < cValues.size(); p++) { - out << setprecision(6) << dist[p] << '\t'; - //print out coverage values - for (int i = 0; i < numGroups; i++) { - for (int j = 0; j < numGroups; j++) { - out << cValues[p][i][j] << '\t'; + int index = 0; + for(int i=0;i prevRow = it->second; + it++; + + for(it;it!=allDistances.end();it++){ + for(int i=0;isecond.size();i++){ + it->second[i] += prevRow[i]; + } + prevRow = it->second; + } + + vector lastRow = allDistances.rbegin()->second; + outCov << setprecision(8); + + outCov << "dist"; + for (int i = 0; i < numGroups; i++){ + outCov << '\t' << groupNames[i]; + } + for (int i=0;ifirst << '\t'; + for(int i=0;isecond[indices[i][i]]/(float)lastRow[indices[i][i]] << '\t'; + } + for(int i=0;isecond[indices[i][j]]/(float)lastRow[indices[i][j]] << '\t'; + outCov << it->second[indices[j][i]]/(float)lastRow[indices[j][i]] << '\t'; + } + } + outCov << endl; } } @@ -241,38 +186,45 @@ void LibShuffCommand::printCoverageFile() { exit(1); } } + //********************************************************************************************************************** + void LibShuffCommand::printSummaryFile() { try { - //format output + + ofstream outSum; + summaryFile = getRootName(globaldata->getPhylipFile()) + "libshuff.summary"; + openOutputFile(summaryFile, outSum); + outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint); + cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint); - for (int i = 0; i < numGroups; i++) { - for (int j = 0; j < numGroups; j++) { - //don't output AA to AA - if (i != j) { - outSum << "Delta " + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t'<< "DeltaSig " + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t'; - cout << "Delta " + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t'<< "DeltaSig " + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t'; + cout << setw(20) << left << "Comparison" << '\t' << setprecision(8) << "dCXYScore" << '\t' << "Significance" << endl; + outSum << setw(20) << left << "Comparison" << '\t' << setprecision(8) << "dCXYScore" << '\t' << "Significance" << endl; + + int precision = (int)log10(iters); + for(int i=0;i (1/(float)iters)) { - outSum << setprecision(6) << sumDelta[i] << '\t' << setprecision(globaldata->getIters().length()) << sumDeltaSig[i] << '\t'; - cout << setprecision(6) << sumDelta[i] << '\t' << setprecision(globaldata->getIters().length()) << sumDeltaSig[i] << '\t'; - }else { - outSum << setprecision(6) << sumDelta[i] << '\t' << setprecision(globaldata->getIters().length()) << "<" << (1/float(iters)) << '\t'; - cout << setprecision(6) << sumDelta[i] << '\t' << setprecision(globaldata->getIters().length()) << "<" << (1/float(iters)) << '\t'; - } - } - outSum << endl; - cout << endl; - } catch(exception& e) { cout << "Standard Error: " << e.what() << " has occurred in the LibShuffCommand class Function printSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; @@ -285,6 +237,7 @@ void LibShuffCommand::printSummaryFile() { } //********************************************************************************************************************** + void LibShuffCommand::setGroups() { try { //if the user has not entered specific groups to analyze then do them all @@ -293,7 +246,7 @@ void LibShuffCommand::setGroups() { for (int i=0; i < numGroups; i++) { globaldata->Groups.push_back(globaldata->gGroupmap->namesOfGroups[i]); } - }else { + } else { if (globaldata->getGroups() != "all") { //check that groups are valid for (int i = 0; i < globaldata->Groups.size(); i++) { @@ -311,8 +264,8 @@ void LibShuffCommand::setGroups() { globaldata->Groups.push_back(globaldata->gGroupmap->namesOfGroups[i]); } cout << "When using the groups parameter you must have at least 2 valid groups. I will run the command using all the groups in your groupfile." << endl; - }else { numGroups = globaldata->Groups.size(); } - }else { //users wants all groups + } else { numGroups = globaldata->Groups.size(); } + } else { //users wants all groups numGroups = globaldata->gGroupmap->getNumGroups(); globaldata->Groups.clear(); for (int i=0; i < numGroups; i++) { @@ -320,21 +273,22 @@ void LibShuffCommand::setGroups() { } } } - + //sort so labels match sort(globaldata->Groups.begin(), globaldata->Groups.end()); //sort sort(globaldata->gGroupmap->namesOfGroups.begin(), globaldata->gGroupmap->namesOfGroups.end()); - + groupNames = globaldata->Groups; + // number of comparisons i.e. with groups A,B,C = AA, AB, AC, BA, BB, BC...; - for (int i=0; iGroups[i] + "-" + globaldata->Groups[l]); - } - } +// for (int i=0; iGroups[i] + "-" + globaldata->Groups[l]); +// } +// } } catch(exception& e) { cout << "Standard Error: " << e.what() << " has occurred in the LibShuffCommand class Function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; @@ -345,23 +299,5 @@ void LibShuffCommand::setGroups() { exit(1); } } -/***********************************************************/ -int LibShuffCommand::findIndex(float score, int index) { - try{ - for (int i = 0; i < rsumDelta[index].size(); i++) { - if (rsumDelta[index][i] >= score) { return i; } - } - return rsumDelta[index].size(); - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the LibShuffCommand class Function findIndex. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the LibShuffCommand class function findIndex. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} /***********************************************************/ - diff --git a/libshuffcommand.h b/libshuffcommand.h index 950bb00..18bfade 100644 --- a/libshuffcommand.h +++ b/libshuffcommand.h @@ -11,8 +11,8 @@ */ #include "command.hpp" -#include "coverage.h" #include "fullmatrix.h" +#include "libshuff.h" using namespace std; @@ -22,33 +22,25 @@ class LibShuffCommand : public Command { public: LibShuffCommand(); - ~LibShuffCommand(); + ~LibShuffCommand(){}; int execute(); private: - vector< vector< vector > > cValues; // vector -one for each distance level. - vector< vector > deltaValues; // vector< vector of delta scores, one for each comparison.> -one at each distance - vector sumDelta; //sum of delta scores, one for each comparison. - vector sumDeltaSig; //number of random matrixes with that delta value or ?? - vector< vector > rsumDelta; //vector< vector > - vector groupComb; - vector dist; - + vector groupNames; void setGroups(); - int findIndex(float, int); void printCoverageFile(); void printSummaryFile(); GlobalData* globaldata; - Coverage* coverage; FullMatrix* matrix; + Libshuff* form; float cutOff, step; int numGroups, numComp, iters; - string coverageFile, summaryFile, form; - ofstream out, outSum; - - + string coverageFile, summaryFile; + vector > pValueCounts; + vector > savedDXYValues; + vector > > savedMinValues; }; #endif diff --git a/progress.cpp b/progress.cpp index 9f3327f..a01ae71 100644 --- a/progress.cpp +++ b/progress.cpp @@ -16,14 +16,33 @@ const int totalTicks = 50; const char marker = '|'; +/***********************************************************************/ + +Progress::Progress(){ + try { + cout << "********************#****#****#****#****#****#****#****#****#****#****#"; + + nTicks = 0; + finalPos = 0; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the Progress class Function Progress. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Progress class function Progress. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + /***********************************************************************/ Progress::Progress(string job, int end){ try { cout << "********************#****#****#****#****#****#****#****#****#****#****#\n"; - cout << job << marker; + cout << setw(20) << left << job << setw(1) << marker; cout.flush(); - + nTicks = 0; finalPos = end; } @@ -39,6 +58,27 @@ Progress::Progress(string job, int end){ /***********************************************************************/ +void Progress::newLine(string job, int end){ + try { + cout << endl; + cout << setw(20) << left << job << setw(1) << marker; + cout.flush(); + + nTicks = 0; + finalPos = end; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the Progress class Function newline. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Progress class function newline. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + +/***********************************************************************/ + void Progress::update(const int currentPos){ try { int ratio = int(totalTicks * (float)currentPos / finalPos); diff --git a/progress.hpp b/progress.hpp index ae0303c..a2ae7d1 100644 --- a/progress.hpp +++ b/progress.hpp @@ -8,8 +8,10 @@ using namespace std; class Progress { public: + Progress(); Progress(string, int); void update(int); + void newLine(string, int); void finish(); private: diff --git a/slibshuff.cpp b/slibshuff.cpp new file mode 100644 index 0000000..a9710de --- /dev/null +++ b/slibshuff.cpp @@ -0,0 +1,116 @@ +/* + * slibshuff.cpp + * Mothur + * + * Created by Pat Schloss on 4/8/09. + * Copyright 2009 Patrick D. Schloss. All rights reserved. + * + */ + +#include "slibshuff.h" + +/***********************************************************************/ + +SLibshuff::SLibshuff(FullMatrix* D, int it, float co) : Libshuff(D, it, 0, co){} + +/***********************************************************************/ + +float SLibshuff::evaluatePair(int i, int j){ + return sCalculate(i,j); +} + +/***********************************************************************/ + +vector > SLibshuff::evaluateAll(){ + try{ + savedMins.resize(numGroups); + + vector > dCXYValues(numGroups); + + for(int i=0;i > evaluateAll(); + float evaluatePair(int, int); + +private: + double sCalculate(int, int); +}; + +#endif diff --git a/validparameter.cpp b/validparameter.cpp index 41ea7cf..98dbde8 100644 --- a/validparameter.cpp +++ b/validparameter.cpp @@ -31,7 +31,6 @@ ValidParameters::ValidParameters() { ValidParameters::~ValidParameters() {} - /***********************************************************************/ bool ValidParameters::isValidParameter(string parameter, string command, string value) { try { @@ -229,7 +228,8 @@ void ValidParameters::initCommandParameters() { try { //{"parameter1","parameter2",...,"last parameter"}; - string readdistArray[] = {"phylip","column","name","cutoff","precision","group"}; + string readdistArray[] = {"phylip","column", "name","cutoff","precision", "group"}; + commandParameters["read.dist"] = addParameters(readdistArray, sizeof(readdistArray)/sizeof(string)); string readotuArray[] = {"list","order","shared", "line", "label","group","sabund", "rabund"};