From eb1c88346fb246e95a6b38935b103f95e38b82ca Mon Sep 17 00:00:00 2001 From: ryabin Date: Fri, 27 Mar 2009 19:25:25 +0000 Subject: [PATCH] added the Calculators Thomas made in the fall. Added parameter and command error checking. --- Mothur.xcodeproj/project.pbxproj | 139 +++---- bergerparker.cpp | 38 ++ bergerparker.h | 32 ++ bstick.cpp | 103 ++++++ bstick.h | 33 ++ calculator.cpp | 343 +++++++---------- calculator.h | 24 +- collectcommand.cpp | 16 +- collectsharedcommand.cpp | 9 + errorchecking.cpp | 14 +- geom.cpp | 120 ++++++ geom.h | 39 ++ globaldata.cpp | 8 +- logsd.cpp | 106 ++++++ logsd.h | 38 ++ qstat.cpp | 76 ++++ qstat.h | 32 ++ rabundvector.cpp | 48 +++ rabundvector.hpp | 5 + sabundvector.cpp | 34 +- sabundvector.hpp | 7 +- sharedbdiversity.cpp | 126 +++++++ sharedbdiversity.h | 33 ++ sharedkstest.cpp | 79 ++++ sharedkstest.h | 29 ++ sharedordervector.h | 3 + sharedrabundvector.cpp | 43 ++- sharedrabundvector.h | 7 +- summarycommand.cpp | 22 +- summarysharedcommand.cpp | 10 +- validcalculator.cpp | 33 +- validparameter.cpp | 612 ++++++++++++++++--------------- validparameter.h | 28 +- 33 files changed, 1637 insertions(+), 652 deletions(-) create mode 100644 bergerparker.cpp create mode 100644 bergerparker.h create mode 100644 bstick.cpp create mode 100644 bstick.h create mode 100644 geom.cpp create mode 100644 geom.h create mode 100644 logsd.cpp create mode 100644 logsd.h create mode 100644 qstat.cpp create mode 100644 qstat.h create mode 100644 sharedbdiversity.cpp create mode 100644 sharedbdiversity.h create mode 100644 sharedkstest.cpp create mode 100644 sharedkstest.h diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj index 3c32966..b34805b 100644 --- a/Mothur.xcodeproj/project.pbxproj +++ b/Mothur.xcodeproj/project.pbxproj @@ -7,11 +7,6 @@ objects = { /* Begin PBXBuildFile section */ - 217595330F791EEE001CC3C6 /* sharedkulczynski.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 217595320F791EEE001CC3C6 /* sharedkulczynski.cpp */; }; - 2175953F0F792111001CC3C6 /* sharedkulczynskicody.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 2175953E0F792111001CC3C6 /* sharedkulczynskicody.cpp */; }; - 217595590F792490001CC3C6 /* sharedlennon.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 217595580F792490001CC3C6 /* sharedlennon.cpp */; }; - 217595640F7927BB001CC3C6 /* sharedmorisitahorn.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 217595630F7927BB001CC3C6 /* sharedmorisitahorn.cpp */; }; - 217595710F792A33001CC3C6 /* sharedbraycurtis.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 217595700F792A33001CC3C6 /* sharedbraycurtis.cpp */; }; 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 */; }; @@ -19,14 +14,6 @@ 3746107E0F4064D100460C57 /* weighted.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3746107D0F4064D100460C57 /* weighted.cpp */; }; 374610830F40652400460C57 /* unweighted.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 374610820F40652400460C57 /* unweighted.cpp */; }; 3746109D0F40657600460C57 /* unifracunweightedcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3746109C0F40657600460C57 /* unifracunweightedcommand.cpp */; }; - 374CD63F0F65832000D90B4A /* libshuffcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 374CD63E0F65832000D90B4A /* libshuffcommand.cpp */; }; - 374CD6F10F65A4C100D90B4A /* coverage.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 374CD6F00F65A4C100D90B4A /* coverage.cpp */; }; - 3765B47A0F77D15900C3EDC5 /* nocommands.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3765B4780F77D15900C3EDC5 /* nocommands.cpp */; }; - 3765B4E20F78055000C3EDC5 /* sharedochiai.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3765B4E10F78055000C3EDC5 /* sharedochiai.cpp */; }; - 3765B5160F780A7F00C3EDC5 /* sharedanderberg.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3765B5150F780A7F00C3EDC5 /* sharedanderberg.cpp */; }; - 3765B5380F7A5A9E00C3EDC5 /* heatmapcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3765B5370F7A5A9E00C3EDC5 /* heatmapcommand.cpp */; }; - 3765B5460F7A62A000C3EDC5 /* heatmap.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3765B5450F7A62A000C3EDC5 /* heatmap.cpp */; }; - 3782163D0F616079008E1F6D /* fullmatrix.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3782163C0F616079008E1F6D /* fullmatrix.cpp */; }; 379293C30F2DE73400B9034A /* treemap.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 379293C20F2DE73400B9034A /* treemap.cpp */; }; 379294700F2E191800B9034A /* parsimonycommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 3792946F0F2E191800B9034A /* parsimonycommand.cpp */; }; 3792948A0F2E258500B9034A /* parsimony.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 379294890F2E258500B9034A /* parsimony.cpp */; }; @@ -108,6 +95,14 @@ A70B53AA0F4CD7AD0064797E /* getgroupcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A70B53A40F4CD7AD0064797E /* getgroupcommand.cpp */; }; A70B53AB0F4CD7AD0064797E /* getlabelcommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A70B53A60F4CD7AD0064797E /* getlabelcommand.cpp */; }; A70B53AC0F4CD7AD0064797E /* getlinecommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A70B53A80F4CD7AD0064797E /* getlinecommand.cpp */; }; + EB1216880F619B83004A865F /* bergerparker.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB1216870F619B83004A865F /* bergerparker.cpp */; }; + EB1216E50F61ACFB004A865F /* bstick.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB1216E40F61ACFB004A865F /* bstick.cpp */; }; + EB1217230F61C9AC004A865F /* sharedkstest.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB1217220F61C9AC004A865F /* sharedkstest.cpp */; }; + EB6E68190F5F1C780044E17B /* qstat.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB6E68180F5F1C780044E17B /* qstat.cpp */; }; + EB6F015B0F6AC1670048081A /* sharedbdiversity.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB6F015A0F6AC1670048081A /* sharedbdiversity.cpp */; }; + EB9303EB0F534D9400E8EF26 /* logsd.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB9303EA0F534D9400E8EF26 /* logsd.cpp */; }; + EB9303F90F53517300E8EF26 /* geom.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EB9303F80F53517300E8EF26 /* geom.cpp */; }; + EBE2A41E0F7D5C6100A54BF8 /* jackknife2.cpp in Sources */ = {isa = PBXBuildFile; fileRef = EBE2A41D0F7D5C6100A54BF8 /* jackknife2.cpp */; }; /* End PBXBuildFile section */ /* Begin PBXCopyFilesBuildPhase section */ @@ -124,18 +119,6 @@ /* End PBXCopyFilesBuildPhase section */ /* Begin PBXFileReference section */ - 217595310F791EEE001CC3C6 /* sharedkulczynski.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedkulczynski.h; sourceTree = ""; }; - 217595320F791EEE001CC3C6 /* sharedkulczynski.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedkulczynski.cpp; sourceTree = ""; }; - 2175953D0F792111001CC3C6 /* sharedkulczynskicody.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedkulczynskicody.h; sourceTree = ""; }; - 2175953E0F792111001CC3C6 /* sharedkulczynskicody.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedkulczynskicody.cpp; sourceTree = ""; }; - 217595570F792490001CC3C6 /* sharedlennon.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedlennon.h; sourceTree = ""; }; - 217595580F792490001CC3C6 /* sharedlennon.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedlennon.cpp; sourceTree = ""; }; - 217595620F7927BB001CC3C6 /* sharedmorisitahorn.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedmorisitahorn.h; sourceTree = ""; }; - 217595630F7927BB001CC3C6 /* sharedmorisitahorn.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedmorisitahorn.cpp; sourceTree = ""; }; - 2175956F0F792A33001CC3C6 /* sharedbraycurtis.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedbraycurtis.h; sourceTree = ""; }; - 217595700F792A33001CC3C6 /* sharedbraycurtis.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedbraycurtis.cpp; sourceTree = ""; }; - 370936DE0F6E7A4A00EB4C2C /* nseqs.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = nseqs.h; sourceTree = ""; }; - 3709370A0F6E7FC100EB4C2C /* sharednseqs.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharednseqs.h; sourceTree = ""; }; 372E126E0F26365B0095CF7E /* readotucommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = readotucommand.h; sourceTree = ""; }; 372E126F0F26365B0095CF7E /* readotucommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = readotucommand.cpp; sourceTree = ""; }; 372E12940F263D5A0095CF7E /* readdistcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = readdistcommand.h; sourceTree = ""; }; @@ -149,22 +132,6 @@ 374610820F40652400460C57 /* unweighted.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = unweighted.cpp; sourceTree = ""; }; 3746109B0F40657600460C57 /* unifracunweightedcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = unifracunweightedcommand.h; sourceTree = ""; }; 3746109C0F40657600460C57 /* unifracunweightedcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = unifracunweightedcommand.cpp; sourceTree = ""; }; - 374CD63D0F65832000D90B4A /* libshuffcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = libshuffcommand.h; sourceTree = ""; }; - 374CD63E0F65832000D90B4A /* libshuffcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = libshuffcommand.cpp; sourceTree = ""; }; - 374CD6EF0F65A4C100D90B4A /* coverage.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = coverage.h; sourceTree = ""; }; - 374CD6F00F65A4C100D90B4A /* coverage.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = coverage.cpp; sourceTree = ""; }; - 3765B4780F77D15900C3EDC5 /* nocommands.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = nocommands.cpp; sourceTree = ""; }; - 3765B4790F77D15900C3EDC5 /* nocommands.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = nocommands.h; sourceTree = ""; }; - 3765B4E00F78055000C3EDC5 /* sharedochiai.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedochiai.h; sourceTree = ""; }; - 3765B4E10F78055000C3EDC5 /* sharedochiai.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedochiai.cpp; sourceTree = ""; }; - 3765B5140F780A7F00C3EDC5 /* sharedanderberg.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedanderberg.h; sourceTree = ""; }; - 3765B5150F780A7F00C3EDC5 /* sharedanderberg.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedanderberg.cpp; sourceTree = ""; }; - 3765B5360F7A5A9E00C3EDC5 /* heatmapcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = heatmapcommand.h; sourceTree = ""; }; - 3765B5370F7A5A9E00C3EDC5 /* heatmapcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = heatmapcommand.cpp; sourceTree = ""; }; - 3765B5440F7A62A000C3EDC5 /* heatmap.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = heatmap.h; sourceTree = ""; }; - 3765B5450F7A62A000C3EDC5 /* heatmap.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = heatmap.cpp; sourceTree = ""; }; - 3782163B0F616079008E1F6D /* fullmatrix.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = fullmatrix.h; sourceTree = ""; }; - 3782163C0F616079008E1F6D /* fullmatrix.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = fullmatrix.cpp; sourceTree = ""; }; 379293C10F2DE73400B9034A /* treemap.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = treemap.h; sourceTree = ""; }; 379293C20F2DE73400B9034A /* treemap.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = treemap.cpp; sourceTree = ""; }; 3792946E0F2E191800B9034A /* parsimonycommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = parsimonycommand.h; sourceTree = ""; }; @@ -337,6 +304,22 @@ A70B53A80F4CD7AD0064797E /* getlinecommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = getlinecommand.cpp; sourceTree = ""; }; A70B53A90F4CD7AD0064797E /* getlinecommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = getlinecommand.h; sourceTree = ""; }; C6859E8B029090EE04C91782 /* Mothur.1 */ = {isa = PBXFileReference; lastKnownFileType = text.man; path = Mothur.1; sourceTree = ""; }; + EB1216860F619B83004A865F /* bergerparker.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = bergerparker.h; sourceTree = ""; }; + EB1216870F619B83004A865F /* bergerparker.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = bergerparker.cpp; sourceTree = ""; }; + EB1216E30F61ACFB004A865F /* bstick.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = bstick.h; sourceTree = ""; }; + EB1216E40F61ACFB004A865F /* bstick.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = bstick.cpp; sourceTree = ""; }; + EB1217210F61C9AC004A865F /* sharedkstest.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedkstest.h; sourceTree = ""; }; + EB1217220F61C9AC004A865F /* sharedkstest.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedkstest.cpp; sourceTree = ""; }; + EB6E68170F5F1C780044E17B /* qstat.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = qstat.h; sourceTree = ""; }; + EB6E68180F5F1C780044E17B /* qstat.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = qstat.cpp; sourceTree = ""; }; + EB6F01590F6AC1670048081A /* sharedbdiversity.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sharedbdiversity.h; sourceTree = ""; }; + EB6F015A0F6AC1670048081A /* sharedbdiversity.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = sharedbdiversity.cpp; sourceTree = ""; }; + EB9303E90F534D9400E8EF26 /* logsd.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = logsd.h; sourceTree = ""; }; + EB9303EA0F534D9400E8EF26 /* logsd.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = logsd.cpp; sourceTree = ""; }; + EB9303F70F53517300E8EF26 /* geom.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = geom.h; sourceTree = ""; }; + EB9303F80F53517300E8EF26 /* geom.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = geom.cpp; sourceTree = ""; }; + EBE2A41C0F7D5C5C00A54BF8 /* jackknife2.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; name = jackknife2.h; path = ../jackknife2.h; sourceTree = SOURCE_ROOT; }; + EBE2A41D0F7D5C6100A54BF8 /* jackknife2.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; name = jackknife2.cpp; path = ../jackknife2.cpp; sourceTree = SOURCE_ROOT; }; /* End PBXFileReference section */ /* Begin PBXFrameworksBuildPhase section */ @@ -388,8 +371,6 @@ 37D927DD0F21331F001D4494 /* fileoutput.cpp */, 37D927E00F21331F001D4494 /* globaldata.hpp */, 37D927DF0F21331F001D4494 /* globaldata.cpp */, - 3765B5440F7A62A000C3EDC5 /* heatmap.h */, - 3765B5450F7A62A000C3EDC5 /* heatmap.cpp */, 37D927E60F21331F001D4494 /* inputdata.h */, 37D927E50F21331F001D4494 /* inputdata.cpp */, 37D927EA0F21331F001D4494 /* kmer.hpp */, @@ -417,6 +398,8 @@ 37D928210F21331F001D4494 /* shared.h */, 37D928200F21331F001D4494 /* shared.cpp */, 37D928420F21331F001D4494 /* singlelinkage.cpp */, + 37D928450F21331F001D4494 /* sparsematrix.hpp */, + 37D928440F21331F001D4494 /* sparsematrix.cpp */, 37D928480F21331F001D4494 /* summarydata.h */, 37D928490F21331F001D4494 /* summarydisplay.h */, 37D9284C0F21331F001D4494 /* utilities.hpp */, @@ -444,13 +427,26 @@ 37D927BB0F21331F001D4494 /* bootstrap.cpp */, 37D927C00F21331F001D4494 /* chao1.h */, 37D927BF0F21331F001D4494 /* chao1.cpp */, - 374CD6EF0F65A4C100D90B4A /* coverage.h */, - 374CD6F00F65A4C100D90B4A /* coverage.cpp */, + EB9303F70F53517300E8EF26 /* geom.h */, + EB9303F80F53517300E8EF26 /* geom.cpp */, + EB6E68170F5F1C780044E17B /* qstat.h */, + EB6E68180F5F1C780044E17B /* qstat.cpp */, + EB9303E90F534D9400E8EF26 /* logsd.h */, + EB9303EA0F534D9400E8EF26 /* logsd.cpp */, + EB1216860F619B83004A865F /* bergerparker.h */, + EB1216870F619B83004A865F /* bergerparker.cpp */, + EB1216E30F61ACFB004A865F /* bstick.h */, + EB1216E40F61ACFB004A865F /* bstick.cpp */, + EB1217210F61C9AC004A865F /* sharedkstest.h */, + EB1217220F61C9AC004A865F /* sharedkstest.cpp */, + EB6F01590F6AC1670048081A /* sharedbdiversity.h */, + EB6F015A0F6AC1670048081A /* sharedbdiversity.cpp */, + EBE2A41C0F7D5C5C00A54BF8 /* jackknife2.h */, + EBE2A41D0F7D5C6100A54BF8 /* jackknife2.cpp */, 37D927E80F21331F001D4494 /* jackknife.h */, 37D927E70F21331F001D4494 /* jackknife.cpp */, 37D927F50F21331F001D4494 /* npshannon.h */, 37D927F40F21331F001D4494 /* npshannon.cpp */, - 370936DE0F6E7A4A00EB4C2C /* nseqs.h */, 379294880F2E258500B9034A /* parsimony.h */, 379294890F2E258500B9034A /* parsimony.cpp */, 37D928020F21331F001D4494 /* rarecalc.h */, @@ -459,10 +455,6 @@ 37D9281E0F21331F001D4494 /* shannon.cpp */, 37D928230F21331F001D4494 /* sharedace.h */, 37D928220F21331F001D4494 /* sharedace.cpp */, - 3765B5140F780A7F00C3EDC5 /* sharedanderberg.h */, - 3765B5150F780A7F00C3EDC5 /* sharedanderberg.cpp */, - 2175956F0F792A33001CC3C6 /* sharedbraycurtis.h */, - 217595700F792A33001CC3C6 /* sharedbraycurtis.cpp */, 37D928250F21331F001D4494 /* sharedchao1.h */, 37D928240F21331F001D4494 /* sharedchao1.cpp */, 37D928290F21331F001D4494 /* sharedjabund.h */, @@ -471,17 +463,6 @@ 37D9282A0F21331F001D4494 /* sharedjclass.cpp */, 37D9282D0F21331F001D4494 /* sharedjest.h */, 37D9282C0F21331F001D4494 /* sharedjest.cpp */, - 217595310F791EEE001CC3C6 /* sharedkulczynski.h */, - 217595320F791EEE001CC3C6 /* sharedkulczynski.cpp */, - 2175953D0F792111001CC3C6 /* sharedkulczynskicody.h */, - 2175953E0F792111001CC3C6 /* sharedkulczynskicody.cpp */, - 217595570F792490001CC3C6 /* sharedlennon.h */, - 217595580F792490001CC3C6 /* sharedlennon.cpp */, - 217595620F7927BB001CC3C6 /* sharedmorisitahorn.h */, - 217595630F7927BB001CC3C6 /* sharedmorisitahorn.cpp */, - 3709370A0F6E7FC100EB4C2C /* sharednseqs.h */, - 3765B4E00F78055000C3EDC5 /* sharedochiai.h */, - 3765B4E10F78055000C3EDC5 /* sharedochiai.cpp */, 37D928350F21331F001D4494 /* sharedsobs.h */, 37D928340F21331F001D4494 /* sharedsobs.cpp */, 37AFC71D0F445386005F492D /* sharedsobscollectsummary.h */, @@ -527,14 +508,8 @@ A70B53A60F4CD7AD0064797E /* getlabelcommand.cpp */, A70B53A90F4CD7AD0064797E /* getlinecommand.h */, A70B53A80F4CD7AD0064797E /* getlinecommand.cpp */, - 3765B5360F7A5A9E00C3EDC5 /* heatmapcommand.h */, - 3765B5370F7A5A9E00C3EDC5 /* heatmapcommand.cpp */, 37D927E40F21331F001D4494 /* helpcommand.h */, 37D927E30F21331F001D4494 /* helpcommand.cpp */, - 374CD63D0F65832000D90B4A /* libshuffcommand.h */, - 374CD63E0F65832000D90B4A /* libshuffcommand.cpp */, - 3765B4790F77D15900C3EDC5 /* nocommands.h */, - 3765B4780F77D15900C3EDC5 /* nocommands.cpp */, 37D927FA0F21331F001D4494 /* parselistcommand.h */, 37D927F90F21331F001D4494 /* parselistcommand.cpp */, 3792946E0F2E191800B9034A /* parsimonycommand.h */, @@ -571,8 +546,6 @@ 37D927D50F21331F001D4494 /* datavector.hpp */, 37D927DC0F21331F001D4494 /* fastamap.h */, 37D927DB0F21331F001D4494 /* fastamap.cpp */, - 3782163B0F616079008E1F6D /* fullmatrix.h */, - 3782163C0F616079008E1F6D /* fullmatrix.cpp */, 37D927E20F21331F001D4494 /* groupmap.h */, 37D927E10F21331F001D4494 /* groupmap.cpp */, 37D927EE0F21331F001D4494 /* listvector.hpp */, @@ -591,8 +564,6 @@ 37D928300F21331F001D4494 /* sharedrabundvector.cpp */, 37D928330F21331F001D4494 /* sharedsabundvector.h */, 37D928320F21331F001D4494 /* sharedsabundvector.cpp */, - 37D928450F21331F001D4494 /* sparsematrix.hpp */, - 37D928440F21331F001D4494 /* sparsematrix.cpp */, 37AD4DB90F28E2FE00AA2D49 /* tree.h */, 37AD4DBA0F28E2FE00AA2D49 /* tree.cpp */, 379293C10F2DE73400B9034A /* treemap.h */, @@ -756,19 +727,14 @@ A70B53AA0F4CD7AD0064797E /* getgroupcommand.cpp in Sources */, A70B53AB0F4CD7AD0064797E /* getlabelcommand.cpp in Sources */, A70B53AC0F4CD7AD0064797E /* getlinecommand.cpp in Sources */, - 3782163D0F616079008E1F6D /* fullmatrix.cpp in Sources */, - 374CD63F0F65832000D90B4A /* libshuffcommand.cpp in Sources */, - 374CD6F10F65A4C100D90B4A /* coverage.cpp in Sources */, - 3765B47A0F77D15900C3EDC5 /* nocommands.cpp in Sources */, - 3765B4E20F78055000C3EDC5 /* sharedochiai.cpp in Sources */, - 3765B5160F780A7F00C3EDC5 /* sharedanderberg.cpp in Sources */, - 217595330F791EEE001CC3C6 /* sharedkulczynski.cpp in Sources */, - 2175953F0F792111001CC3C6 /* sharedkulczynskicody.cpp in Sources */, - 217595590F792490001CC3C6 /* sharedlennon.cpp in Sources */, - 217595640F7927BB001CC3C6 /* sharedmorisitahorn.cpp in Sources */, - 217595710F792A33001CC3C6 /* sharedbraycurtis.cpp in Sources */, - 3765B5380F7A5A9E00C3EDC5 /* heatmapcommand.cpp in Sources */, - 3765B5460F7A62A000C3EDC5 /* heatmap.cpp in Sources */, + EB9303EB0F534D9400E8EF26 /* logsd.cpp in Sources */, + EB9303F90F53517300E8EF26 /* geom.cpp in Sources */, + EB6E68190F5F1C780044E17B /* qstat.cpp in Sources */, + EB1216880F619B83004A865F /* bergerparker.cpp in Sources */, + EB1216E50F61ACFB004A865F /* bstick.cpp in Sources */, + EB1217230F61C9AC004A865F /* sharedkstest.cpp in Sources */, + EB6F015B0F6AC1670048081A /* sharedbdiversity.cpp in Sources */, + EBE2A41E0F7D5C6100A54BF8 /* jackknife2.cpp in Sources */, ); runOnlyForDeploymentPostprocessing = 0; }; @@ -782,7 +748,7 @@ GCC_DYNAMIC_NO_PIC = NO; GCC_ENABLE_FIX_AND_CONTINUE = YES; GCC_MODEL_TUNING = G5; - GCC_OPTIMIZATION_LEVEL = 3; + GCC_OPTIMIZATION_LEVEL = 0; GCC_PREPROCESSOR_DEFINITIONS = ( "_GLIBCXX_DEBUG=1", "_GLIBCXX_DEBUG_PEDANTIC=1", @@ -798,7 +764,6 @@ buildSettings = { DEBUG_INFORMATION_FORMAT = "dwarf-with-dsym"; GCC_MODEL_TUNING = G5; - GCC_OPTIMIZATION_LEVEL = 3; INSTALL_PATH = /usr/local/bin; PRODUCT_NAME = mothur; }; @@ -810,7 +775,6 @@ GCC_OPTIMIZATION_LEVEL = 3; GCC_WARN_ABOUT_RETURN_TYPE = YES; GCC_WARN_UNUSED_VARIABLE = YES; - INSTALL_MODE_FLAG = "a-w,a+rX"; PREBINDING = NO; SDKROOT = "$(DEVELOPER_SDK_DIR)/MacOSX10.5.sdk"; }; @@ -826,7 +790,6 @@ GCC_OPTIMIZATION_LEVEL = 0; GCC_WARN_ABOUT_RETURN_TYPE = YES; GCC_WARN_UNUSED_VARIABLE = YES; - INSTALL_MODE_FLAG = "a-w,a+rX"; PREBINDING = NO; SDKROOT = "$(DEVELOPER_SDK_DIR)/MacOSX10.5.sdk"; }; diff --git a/bergerparker.cpp b/bergerparker.cpp new file mode 100644 index 0000000..c8c25ea --- /dev/null +++ b/bergerparker.cpp @@ -0,0 +1,38 @@ +/* + * ssbp.cpp + * Mothur + * + * Created by Thomas Ryabin on 3/6/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ + +#include "bergerparker.h" +#include "calculator.h" + +/***************************************************************/ + +EstOutput BergerParker::getValues(SAbundVector* rank){ + try { + data.resize(1,0); + //Berger-Parker index + double BP = (double)rank->getNumSeqs()/(double)rank->getMaxRank(); + //cout << "BP index = " << 1/BP << "\n\n"; + + data[0] = BP; + if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } + + return data; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + +/***********************************************************************/ + diff --git a/bergerparker.h b/bergerparker.h new file mode 100644 index 0000000..f164ec8 --- /dev/null +++ b/bergerparker.h @@ -0,0 +1,32 @@ +#ifndef BERGERPARKER_H +#define BERGERPARKER_H +/* + * bergerparker.h + * Mothur + * + * Created by Thomas Ryabin on 3/6/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ + +#include "calculator.h" + +/*This class implements the SSBP estimator on single group. +It is a child of the calculator class.*/ + +/***********************************************************************/ + +class BergerParker : public Calculator { + +public: + BergerParker() : Calculator("bergerparker", 1) {}; + EstOutput getValues(SAbundVector*); + EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;}; + +private: +}; + +/***********************************************************************/ + +#endif + diff --git a/bstick.cpp b/bstick.cpp new file mode 100644 index 0000000..3a0e0e9 --- /dev/null +++ b/bstick.cpp @@ -0,0 +1,103 @@ +/* + * bstick.cpp + * Mothur + * + * Created by Thomas Ryabin on 3/6/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ + +#include "bstick.h" +#include "calculator.h" + +/***********************************************************************/ +double BStick::invSum(int index, double numSpec) +{ + double sum = 0; + for(int i = index; i <= numSpec; i++) + sum += 1/(double)i; + return sum; +} + +RAbundVector BStick::getRAbundVector(SAbundVector* rank){ + vector rData; + int mr = 1; + int nb = 0; + int ns = 0; + + for(int i = rank->size()-1; i > 0; i--) + { + int cur = rank->get(i); + if(mr == 1 && cur > 0) + mr = i; + nb += cur; + ns += i*cur; + for(int j = 0; j < cur; j++) + rData.push_back(i); + } + + RAbundVector rav = RAbundVector(rData, mr, nb, ns); + rav.setLabel(rank->getLabel()); + return rav; +} + +/***************************************************************************/ + +/***************************************************************************/ +EstOutput BStick::getValues(SAbundVector* rank){ + try { + data.resize(2,0); + rdata = getRAbundVector(rank); + double numInd = (double)rdata.getNumSeqs(); + double numSpec = (double)rdata.getNumBins(); + + double sumExp = 0; + double sumObs = 0; + double maxDiff = 0; + + for(int i = 0; i < rdata.size(); i++) + { + sumObs += rdata.get(i); + sumExp += numInd/numSpec*invSum(i+1,numSpec); + double diff = fabs(sumObs-sumExp); + if(diff > maxDiff) + maxDiff = diff; + } + + double DStatistic = maxDiff/numInd; + double critVal = 0; + /*cout << "BStick:\n"; + cout << "D-Statistic = " << DStatistic << "\n"; + cout << "Critical value for 95% confidence interval = ";*/ + if(rdata.size() > 20) + { + critVal = .886/sqrt(rdata.size()); + } + else + { + KOSTable table; + critVal = table.getConfLimit(numSpec); + } + /*cout << critVal << "\n"; + cout << "If D-Statistic is less than the critical value then the data fits the Broken Stick model w/ 95% confidence.\n\n";*/ + + data[0] = DStatistic; + data[1] = critVal; + if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } + if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; } + + return data; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + +/***********************************************************************/ + + diff --git a/bstick.h b/bstick.h new file mode 100644 index 0000000..f2107b1 --- /dev/null +++ b/bstick.h @@ -0,0 +1,33 @@ +#ifndef BSTICK_H +#define BSTICK_H +/* + * bstick.h + * Mothur + * + * Created by Thomas Ryabin on 3/6/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ +#include "calculator.h" + +/*This class implements the BStick estimator on single group. +It is a child of the calculator class.*/ + +/***********************************************************************/ + +class BStick : public Calculator { + +public: + BStick() : Calculator("bstick", 3) {}; + EstOutput getValues(SAbundVector*); + EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;}; + +private: + double invSum(int, double); + RAbundVector getRAbundVector(SAbundVector*); + RAbundVector rdata; +}; + +/***********************************************************************/ + +#endif diff --git a/calculator.cpp b/calculator.cpp index 6bac62d..ffd452a 100644 --- a/calculator.cpp +++ b/calculator.cpp @@ -82,6 +82,38 @@ double VecCalc::stError(vector vec){ exit(1); } } +/***********************************************************************/ +int VecCalc::sumElements(vector vec){ + try { + return sumElements(vec,0); + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the VecCalc class function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} +/***********************************************************************/ +int VecCalc::sumElements(vector vec, int index){ + try { + int sum = 0; + for(int i = index; i < vec.size(); i++) + sum += vec.at(i); + return sum; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the VecCalc class function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + /***********************************************************************/ double VecCalc::sumElements(vector vec){ try { @@ -135,6 +167,24 @@ double VecCalc::findMax(vector vec){ } } /***********************************************************************/ +int VecCalc::numNZ(vector vec){ + try { + int numNZ = 0; + for(int i = 0; i < vec.size(); i++) + if(vec.at(i) != 0) + numNZ++; + return numNZ; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function numNZ. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the VecCalc class function numNZ. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} +/***********************************************************************/ double VecCalc::numNZ(vector vec){ try { double numNZ = 0; @@ -510,181 +560,6 @@ vector VecCalc::getSData(char name[]){ } } /***********************************************************************/ -double BDiversity::getWhitt(vector vec1, vector vec2){ - try { - VecCalc vecCalc; - double numSpec = vecCalc.numNZ(vecCalc.addVecs(vec1,vec2)); - return 2*numSpec/(vecCalc.numNZ(vec1)+vecCalc.numNZ(vec2))-1; - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getWhitt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the BDiversity class function getWhitt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double BDiversity::getMS(vector vec1, vector vec2){ - try { - VecCalc vecCalc; - double a = vecCalc.numNZ(vecCalc.multVecs(vec1, vec2)); - double b = vecCalc.numPos(vecCalc.addVecs(vecCalc.addVecs(vec1, vecCalc.multiply(vec2, -1)), vecCalc.multiply(vecCalc.multVecs(vec1, vec2), -1))); - double c = vecCalc.numPos(vecCalc.addVecs(vecCalc.addVecs(vec2, vecCalc.multiply(vec1, -1)), vecCalc.multiply(vecCalc.multVecs(vec2, vec1), -1))); - return a/(a+b+c); - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getMS. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the BDiversity class function getMS. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double BDiversity::getSor(vector vec1, vector vec2){ - try { - double sum = 0; - double asum = 0; - double bsum = 0; - for(int i = 0; i < vec1.size(); i++) - { - asum += vec1.at(i); - bsum += vec2.at(i); - if(vec1.at(i) >= vec2.at(i)) - sum += vec2.at(i); - else - sum += vec1.at(i); - } - return 2*sum/(asum+bsum); - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getSor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the BDiversity class function getSor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -double BDiversity::getMor(vector vec1, vector vec2){ - try { - double sum = 0; - double asum = 0; - double bsum = 0; - double sum1 = 0; - double sum2 = 0; - for(int i = 0; i < vec1.size(); i++) - { - sum += vec1.at(i)*vec2.at(i); - asum += pow(vec1.at(i),2); - bsum += pow(vec2.at(i),2); - sum1 += vec1.at(i); - sum2 += vec2.at(i); - } - return 2*sum/((asum/pow(sum1,2)+bsum/pow(sum2,2))*(sum1*sum2)); - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getMor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the BDiversity class function getMor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -void BDiversity::printD(vector > columns, int type){ - try { - cout << " "; - for(int i = 0; i < columns.size(); i++) - cout << "Column:" << i << " "; - cout << "\n"; - for(int i = 0; i < columns.size(); i++) - { - cout << "Column " << i << ":"; - for(int j = 0; j < columns.size(); j++) - { - if(j > i) - { - double B; - if(type == 1) - B = getWhitt(columns.at(i), columns.at(j)); - else if(type == 2) - B = 1-getMS(columns.at(i), columns.at(j)); - else if(type == 3) - B = 1-getSor(columns.at(i), columns.at(j)); - else if(type == 4) - B = 1-getMor(columns.at(i), columns.at(j)); - - cout << B << " "; - } - else - cout << " "; - } - cout << "\n"; - } - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function printD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the BDiversity class function printD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} -/***********************************************************************/ -void BDiversity::doBD(vector vec, double cols)//vec = The data vector if the data table was read left-to-right, top-to-bottom. cols = the number of columns in the data table. -{ try { - VecCalc vecCalc; - vector > columns = vecCalc.gen2DVec(vec,cols,1); - - cout.setf(ios_base::fixed, ios_base::floatfield); - cout.precision(6);//This formats the data so the tables look pretty. - - //Whittaker's measure Bw (presence/absence data) - cout << "Whittaker's measure Bw (presence/absence data):\n"; - printD(columns, 1); - double sum = 0; - for(int i = 0; i < cols; i++) - sum += vecCalc.numNZ(columns.at(i)); - double meanRichness = sum/cols; - vector totVec = vecCalc.genTotVec(columns); - double totRichness = vecCalc.numNZ(totVec); - cout << "\nOverall B Diversity = " << totRichness/meanRichness << "\n\n\n";//The overall B Diversity - - //Marczewski-Steinhaus(MS) distance(Jaccard index) (presence/abscence data) - cout << "Marczewski-Steinhaus(MS) distance(Jaccard index) (presence/abscence data):\n"; - printD(columns, 2); - sum = 0; - for(int i = 0; i < cols; i++) - for(int j = 0; j < cols; j++) - if(j > i) - sum += vecCalc.numNZ(columns.at(i))+vecCalc.numNZ(columns.at(j)) - 2*vecCalc.numNZ(vecCalc.multVecs(columns.at(i), columns.at(j))); - cout << "\nOverall B Diversity = " << sum/cols << "\n\n\n";//The overall B Diversity - - //Sorensen quantitative index (abundance data) - cout << "Sorensen quantitative index (abundance data):\n"; - printD(columns, 3); - cout << "\n\n"; - - //Sorensen quantitative index (abundance data) - cout << "Morisita-Horn index (abundance data):\n"; - printD(columns,4); - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function doBD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the BDiversity class function doBD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } -} /***********************************************************************/ void BrokenStick::doBStick(vector vec)//vec = The data vector. @@ -731,7 +606,7 @@ double kEq(double k, double spec) return k/(1-k)*pow(1-k, spec)/(1-pow(1-k, spec)); } /***********************************************************************/ -void GeometricSeries::doGeomTest(vector vec)//vec = the data vector +/*void GeometricSeries::doGeomTest(vector vec)//vec = the data vector { try { VecCalc vecCalc; vec = vecCalc.quicksort(vec);//Sorts vec @@ -775,7 +650,7 @@ void GeometricSeries::doGeomTest(vector vec)//vec = the data vector cout << "An unknown error has occurred in the GeometricSeries class function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; exit(1); } -} +}*/ /***********************************************************************/ void Jackknifing::doJK(vector vec, double cols)//vec = the data vector if the data table was read left-to-right, top-to-bottom. cols = The number of columns in the data table. @@ -877,7 +752,7 @@ double logS(double num) } /***********************************************************************/ -void LogSD::doLogSD(vector indVec, vector specVec) //indVec = individuals vector, specVec = species vector +/*void LogSD::doLogSD(vector indVec, vector specVec) //indVec = individuals vector, specVec = species vector { try { VecCalc vecCalc; double numSpec = vecCalc.sumElements(specVec);//numSpec = The total number of species @@ -959,7 +834,7 @@ void LogSD::doLogSD(vector indVec, vector specVec) //indVec = in cout << "An unknown error has occurred in the LogSD class function doLogSD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; exit(1); } -} +}*/ /***********************************************************************/ void QStatistic::doQStat(vector vec)//vec = The data vector. { try { @@ -1074,40 +949,42 @@ void SSBPDiversityIndices::doSSBP(vector vec)//vec = The data vector. /***********************************************************************/ double TDTable::getConfLimit(int row, int col) //Rows are the degrees of freedom { try { - //Confidence Level .80 .90 .95 .98 .99 .998 .999 + //Found on http://www.vgtu.lt/leidiniai/elektroniniai/Probability.pdf/Table%203.pdf + + //Confidence Level .90 .95 .975 .99 .995 .999 .9995 double values[33][7] = {{3.078, 6.314, 12.706, 31.821, 63.656, 318.289, 636.578}, - {1.886, 2.920, 4.303, 6.965, 9.925, 22.328, 31.600}, - {1.638, 2.353, 3.182, 4.541, 5.841, 10.214, 12.924}, - {1.533, 2.132, 2.776, 3.747, 4.604, 7.173, 8.610}, - {1.476, 2.015, 2.571, 3.365, 4.032, 5.894, 6.869}, - {1.440, 1.943, 2.447, 3.143, 3.707, 5.208, 5.959}, - {1.415, 1.895, 2.365, 2.998, 3.499, 4.785, 5.408}, - {1.397, 1.860, 2.306, 2.896, 3.355, 4.501, 5.041}, - {1.383, 1.833, 2.262, 2.821, 3.250, 4.297, 4.781}, - {1.372, 1.812, 2.228, 2.764, 3.169, 4.144, 4.587}, - {1.363, 1.796, 2.201, 2.718, 3.106, 4.025, 4.437}, - {1.356, 1.782, 2.179, 2.681, 3.055, 3.930, 4.318}, - {1.350, 1.771, 2.160, 2.650, 3.012, 3.852, 4.221}, - {1.345, 1.761, 2.145, 2.624, 2.977, 3.787, 4.140}, - {1.341, 1.753, 2.131, 2.602, 2.947, 3.733, 4.073}, - {1.337, 1.746, 2.120, 2.583, 2.921, 3.686, 4.015}, - {1.333, 1.740, 2.110, 2.567, 2.898, 3.646, 3.965}, - {1.330, 1.734, 2.101, 2.552, 2.878, 3.610, 3.922}, - {1.328, 1.729, 2.093, 2.539, 2.861, 3.579, 3.883}, - {1.325, 1.725, 2.086, 2.528, 2.845, 3.552, 3.850}, - {1.323, 1.721, 2.080, 2.518, 2.831, 3.527, 3.819}, - {1.321, 1.717, 2.074, 2.508, 2.819, 3.505, 3.792}, - {1.319, 1.714, 2.069, 2.500, 2.807, 3.485, 3.768}, - {1.318, 1.711, 2.064, 2.492, 2.797, 3.467, 3.745}, - {1.316, 1.708, 2.060, 2.485, 2.787, 3.450, 3.725}, - {1.315, 1.706, 2.056, 2.479, 2.779, 3.435, 3.707}, - {1.314, 1.703, 2.052, 2.473, 2.771, 3.421, 3.689}, - {1.313, 1.701, 2.048, 2.467, 2.763, 3.408, 3.674}, - {1.311, 1.699, 2.045, 2.462, 2.756, 3.396, 3.660}, - {1.310, 1.697, 2.042, 2.457, 2.750, 3.385, 3.646}, - {1.296, 1.671, 2.000, 2.390, 2.660, 3.232, 3.460}, - {1.289, 1.658, 1.980, 2.358, 2.617, 3.160, 3.373}, - {1.282, 1.645, 1.960, 2.326, 2.576, 3.091, 3.291}}; + {1.886, 2.920, 4.303, 6.965, 9.925, 22.328, 31.600}, + {1.638, 2.353, 3.182, 4.541, 5.841, 10.214, 12.924}, + {1.533, 2.132, 2.776, 3.747, 4.604, 7.173, 8.610}, + {1.476, 2.015, 2.571, 3.365, 4.032, 5.894, 6.869}, + {1.440, 1.943, 2.447, 3.143, 3.707, 5.208, 5.959}, + {1.415, 1.895, 2.365, 2.998, 3.499, 4.785, 5.408}, + {1.397, 1.860, 2.306, 2.896, 3.355, 4.501, 5.041}, + {1.383, 1.833, 2.262, 2.821, 3.250, 4.297, 4.781}, + {1.372, 1.812, 2.228, 2.764, 3.169, 4.144, 4.587}, + {1.363, 1.796, 2.201, 2.718, 3.106, 4.025, 4.437}, + {1.356, 1.782, 2.179, 2.681, 3.055, 3.930, 4.318}, + {1.350, 1.771, 2.160, 2.650, 3.012, 3.852, 4.221}, + {1.345, 1.761, 2.145, 2.624, 2.977, 3.787, 4.140}, + {1.341, 1.753, 2.131, 2.602, 2.947, 3.733, 4.073}, + {1.337, 1.746, 2.120, 2.583, 2.921, 3.686, 4.015}, + {1.333, 1.740, 2.110, 2.567, 2.898, 3.646, 3.965}, + {1.330, 1.734, 2.101, 2.552, 2.878, 3.610, 3.922}, + {1.328, 1.729, 2.093, 2.539, 2.861, 3.579, 3.883}, + {1.325, 1.725, 2.086, 2.528, 2.845, 3.552, 3.850}, + {1.323, 1.721, 2.080, 2.518, 2.831, 3.527, 3.819}, + {1.321, 1.717, 2.074, 2.508, 2.819, 3.505, 3.792}, + {1.319, 1.714, 2.069, 2.500, 2.807, 3.485, 3.768}, + {1.318, 1.711, 2.064, 2.492, 2.797, 3.467, 3.745}, + {1.316, 1.708, 2.060, 2.485, 2.787, 3.450, 3.725}, + {1.315, 1.706, 2.056, 2.479, 2.779, 3.435, 3.707}, + {1.314, 1.703, 2.052, 2.473, 2.771, 3.421, 3.689}, + {1.313, 1.701, 2.048, 2.467, 2.763, 3.408, 3.674}, + {1.311, 1.699, 2.045, 2.462, 2.756, 3.396, 3.660}, + {1.310, 1.697, 2.042, 2.457, 2.750, 3.385, 3.646}, + {1.296, 1.671, 2.000, 2.390, 2.660, 3.232, 3.460}, + {1.289, 1.658, 1.980, 2.358, 2.617, 3.160, 3.373}, + {1.282, 1.645, 1.960, 2.326, 2.576, 3.091, 3.291}}; return values[row][col]; } catch(exception& e) { @@ -1119,6 +996,48 @@ double TDTable::getConfLimit(int row, int col) //Rows are the degrees of freedom exit(1); } } + +/***********************************************************************/ +double KOSTable::getConfLimit(int index) //Table of Critical values for N = 4-20 for One-Sample Kolmogorov-Smirnov Test +{ try { + //Confidence Level = .05 + //Sample size = 4-20. + //Found on http://www.utdallas.edu/~herve/MolinAbdi1998-LillieforsTechReport.pdf + + double values[21] = {.9011, + .6372, + .5203, + .4506, + 0.3754, + 0.3427, + 0.3245, + 0.3041, + 0.2875, + 0.2744, + 0.2616, + 0.2506, + 0.2426, + 0.2337, + 0.2257, + 0.2196, + 0.2128, + 0.2071, + 0.2018, + 0.1965, + 0.1920, + }; + return values[index]; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the TDTable class Function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the TDTable class function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + /*********************************************************************** void TrunLN::doTrunLN(vector indVec, vector specVec) //indVec = individuals vector, specVec = species vector { diff --git a/calculator.h b/calculator.h index 8c935af..c87d5c7 100644 --- a/calculator.h +++ b/calculator.h @@ -47,14 +47,18 @@ class VecCalc { // The methods seen in the order here is how they are ordered throughout the class. public: + VecCalc(){}; void printElements(vector); //This prints the values of the vector on one line with a space between each value. void printElements(vector); //This prints the values of the vector on one line with a space between each value. int findString(vector, string);//This returns the index of the given string in the given vector, if the string does not exist in the vector it returns -1. double mean(vector); //This returns the mean value of the vector. double stError(vector); //This returns the standard error of the vector. + int sumElements(vector, int); + int sumElements(vector); double sumElements(vector); //This returns the sum of all the values in the vector. double sumElements(vector, int); //This returns the sum of all the values in the vector excluding those whose index is before the given index. double findMax(vector); //This returns the maximum value in the vector. + int numNZ(vector); //This returns the number of non-zero values in the vector. double numNZ(vector); //This returns the number of non-zero values in the vector. double numPos(vector); //This returns the number of positive values in the vector. double findMaxDiff(vector, vector); //This returns the absolute value of the maximum difference between the two vectors. @@ -86,7 +90,7 @@ each combination. It also calculates the overall diversity for Whittaker's measu the Marczewski-Steinhaus distance.*/ -class BDiversity +/*class BDiversity { public: void doBD(vector, double);//Main method @@ -95,7 +99,7 @@ class BDiversity double getSor(vector, vector);//Sorensen quantitative index double getMor(vector, vector);//Morisita-Horn index void printD(vector >, int);//This prints a table that represents the given 2D vector, the second paramter specifies which method is to be used (1 for Whitt, 2 for MS, 3 for Sor, and 4 for Mor) -}; +};*/ /**************************************************************************************************/ @@ -115,11 +119,11 @@ class BrokenStick It prints the D-Statistic and the critical values for the Kolmogorov-Smirnov 1-sample test at the 95% confidence interval.*/ -class GeometricSeries +/*class GeometricSeries { public: void doGeomTest(vector); -}; +};*/ /**************************************************************************************************/ //This Class calculates the jackknifed estimate of the data and @@ -145,11 +149,11 @@ class KS2SampleTest It then generates a D-Statistic and prints the D-Statistic and the critical values for the Kolmogorov-Smirnov 1 sample test.*/ -class LogSD +/*class LogSD { public: void doLogSD(vector, vector); -}; +};*/ /**************************************************************************************************/ //This Class calculates and prints the Q-Statistic for the data. @@ -175,6 +179,14 @@ class TDTable double getConfLimit(int,int); }; +/**************************************************************************************************/ +//This Class stores the table of the confidence limits of the One-Sample Kolmogorov-Smirnov Test. +class KOSTable +{ + public: + double getConfLimit(int); +}; + /**************************************************************************************************/ /*This Class calculates the truncated lognormal for the data. It then prints the D-Statistic and the critical values for the diff --git a/collectcommand.cpp b/collectcommand.cpp index 000f85d..b3e8d8f 100644 --- a/collectcommand.cpp +++ b/collectcommand.cpp @@ -17,6 +17,11 @@ #include "npshannon.h" #include "shannon.h" #include "jackknife.h" +#include "geom.h" +#include "qstat.h" +#include "logsd.h" +#include "bergerparker.h" +#include "bstick.h" //********************************************************************************************************************** @@ -29,7 +34,6 @@ CollectCommand::CollectCommand(){ fileNameRoot = getRootName(globaldata->inputFileName); int i; validCalculator = new ValidCalculators(); - for (i=0; iEstimators.size(); i++) { if (validCalculator->isValidCalculator("single", globaldata->Estimators[i]) == true) { if (globaldata->Estimators[i] == "sobs") { @@ -51,6 +55,16 @@ CollectCommand::CollectCommand(){ cDisplays.push_back(new CollectDisplay(new Simpson(), new ThreeColumnFile(fileNameRoot+"simpson"))); }else if (globaldata->Estimators[i] == "bootstrap") { cDisplays.push_back(new CollectDisplay(new Bootstrap(), new OneColumnFile(fileNameRoot+"bootstrap"))); + }else if (globaldata->Estimators[i] == "geom") { + cDisplays.push_back(new CollectDisplay(new Geom(), new OneColumnFile(fileNameRoot+"geom"))); + }else if (globaldata->Estimators[i] == "qstat") { + cDisplays.push_back(new CollectDisplay(new QStat(), new OneColumnFile(fileNameRoot+"qstat"))); + }else if (globaldata->Estimators[i] == "logsd") { + cDisplays.push_back(new CollectDisplay(new LogSD(), new OneColumnFile(fileNameRoot+"logsd"))); + }else if (globaldata->Estimators[i] == "bergerparker") { + cDisplays.push_back(new CollectDisplay(new BergerParker(), new OneColumnFile(fileNameRoot+"bergerparker"))); + }else if (globaldata->Estimators[i] == "bstick") { + cDisplays.push_back(new CollectDisplay(new BStick(), new OneColumnFile(fileNameRoot+"bstick"))); } } } diff --git a/collectsharedcommand.cpp b/collectsharedcommand.cpp index 6ee2943..01a9343 100644 --- a/collectsharedcommand.cpp +++ b/collectsharedcommand.cpp @@ -19,6 +19,10 @@ #include "sharedsorest.h" #include "sharedthetayc.h" #include "sharedthetan.h" +<<<<<<< collectsharedcommand.cpp +#include "sharedkstest.h" +#include "sharedbdiversity.h" +======= #include "sharednseqs.h" #include "sharedochiai.h" #include "sharedanderberg.h" @@ -27,6 +31,7 @@ #include "sharedlennon.h" #include "sharedmorisitahorn.h" #include "sharedbraycurtis.h" +>>>>>>> 1.13 //********************************************************************************************************************** @@ -64,6 +69,10 @@ CollectSharedCommand::CollectSharedCommand(){ cDisplays.push_back(new CollectDisplay(new SharedThetaYC(), new SharedOneColumnFile(fileNameRoot+"shared.thetayc"))); }else if (globaldata->Estimators[i] == "sharedthetan") { cDisplays.push_back(new CollectDisplay(new SharedThetaN(), new SharedOneColumnFile(fileNameRoot+"shared.thetan"))); + }else if (globaldata->Estimators[i] == "sharedkstest") { + cDisplays.push_back(new CollectDisplay(new SharedKSTest(), new SharedOneColumnFile(fileNameRoot+"shared.kstest"))); + }else if (globaldata->Estimators[i] == "sharedbdiversity") { + cDisplays.push_back(new CollectDisplay(new SharedBDiversity(), new SharedOneColumnFile(fileNameRoot+"shared.bdiversity"))); }else if (globaldata->Estimators[i] == "sharednseqs") { cDisplays.push_back(new CollectDisplay(new SharedNSeqs(), new SharedOneColumnFile(fileNameRoot+"shared.nseqs"))); }else if (globaldata->Estimators[i] == "sharedochiai") { diff --git a/errorchecking.cpp b/errorchecking.cpp index b9b3100..59095f7 100644 --- a/errorchecking.cpp +++ b/errorchecking.cpp @@ -83,8 +83,10 @@ bool ErrorCheck::checkInput(string input) { splitAtComma(value, optionText); splitAtEquals(parameter, value); - //is it a valid parameter for the command - if (validParameter->isValidParameter(parameter, commandName) != true) { return false; } + //is it a valid parameter + if (validParameter->isValidParameter(parameter, commandName, value) != true) { return false; } + + if (parameter == "phylip" ) { phylipfile = value; } if (parameter == "column" ) { columnfile = value; } @@ -114,9 +116,9 @@ bool ErrorCheck::checkInput(string input) { if (errorFree) { //gets the last parameter and value value = optionText; splitAtEquals(parameter, value); - - //is it a valid parameter for the command - if (validParameter->isValidParameter(parameter, commandName) != true) { return false; } + //is it a valid parameter + if (validParameter->isValidParameter(parameter, commandName, value) != true) { return false; } + if (parameter == "phylip" ) { phylipfile = value; } if (parameter == "column" ) { columnfile = value; } @@ -442,7 +444,7 @@ void ErrorCheck::validateReadPhil() { }else if (sabundfile != "") { if ((listfile != "") || (rabundfile != "")) { cout << "When executing a read.otu you must enter ONLY ONE of the following: list, rabund or sabund." << endl; errorFree = false; } - }else if ((listfile == "") && (rabundfile == "") && (sabundfile == "")) { + }else if ((listfile == "") && (rabundfile == "") && (sabundfile == "") && (sharedfile == "")) { cout << "When executing a read.otu you must enter one of the following: list, rabund or sabund." << endl; errorFree = false; } diff --git a/geom.cpp b/geom.cpp new file mode 100644 index 0000000..8b484de --- /dev/null +++ b/geom.cpp @@ -0,0 +1,120 @@ +/* + * geom.cpp + * Mothur + * + * Created by Thomas Ryabin on 2/23/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ + +#include "geom.h" +#include "calculator.h" + +/***********************************************************************/ + +double Geom::kEq(double k, double spec){ + return k/(1-k)*pow(1-k, spec)/(1-pow(1-k, spec)); +} + +RAbundVector Geom::getRAbundVector(SAbundVector* rank){ + vector rData; + int mr = 1; + int nb = 0; + int ns = 0; + + for(int i = rank->size()-1; i > 0; i--) + { + int cur = rank->get(i); + if(mr == 1 && cur > 0) + mr = i; + nb += cur; + ns += i*cur; + for(int j = 0; j < cur; j++) + rData.push_back(i); + } + + RAbundVector rav = RAbundVector(rData, mr, nb, ns); + rav.setLabel(rank->getLabel()); + return rav; +} + +/***********************************************************************************/ + +/***********************************************************************************/ +EstOutput Geom::getValues(SAbundVector* rank){ + try { + data.resize(2,0); + + rdata = getRAbundVector(rank); + int numInd = rdata.getNumSeqs(); + int numSpec = rdata.getNumBins(); + int min = rdata.get(rdata.size()-1); + double k = .5; + double step = .49999; + + while(fabs(min - numInd*kEq(k, (double)numSpec)) > .0001) //This uses a binary search to find the value of k. + { + if(numInd*kEq(k, numSpec) > min) + k += step; + else + k -= step; + step /= 2; + } + + double cK = 1/(1-pow(1-k, numSpec)); + double sumExp = 0; + double sumObs = 0; + double maxDiff = 0; + + for(int i = 0; i < rdata.size(); i++) + { + sumObs += rdata.get(i); + sumExp += numInd*cK*k*pow(1-k, i); + double diff = fabs(sumObs-sumExp); + if(diff > maxDiff) + maxDiff = diff; + } + + double DStatistic = maxDiff/numInd; + double critVal = 0; + /*cout << "Geom:\n"; + cout << "D-Statistic = " << DStatistic << "\n"; + cout << "Critical value for 95% confidence interval = ";*/ + if(rdata.size() > 20) + { + critVal = .886/sqrt(rdata.size()); + } + else + { + KOSTable table; + critVal = table.getConfLimit(numSpec); + } + /*cout << critVal << "\n"; + cout << "If D-Statistic is less than the critical value then the data fits the Geometric Series model w/ 95% confidence.\n\n";*/ + + data[0] = DStatistic; + data[1] = critVal; + + if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } + if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; } + + return data; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + +/***********************************************************************/ + + + + + + + diff --git a/geom.h b/geom.h new file mode 100644 index 0000000..56bc112 --- /dev/null +++ b/geom.h @@ -0,0 +1,39 @@ +#ifndef GEOM_H +#define GEOM_H + +/* + * geom.h + * Mothur + * + * Created by Thomas Ryabin on 2/23/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ +#include "calculator.h" + +/* This class implements the LogSD estimator on single group. +It is a child of the calculator class. */ + +/***********************************************************************/ + +class Geom : public Calculator { + +public: + Geom() : Calculator("geom", 3) {}; + EstOutput getValues(SAbundVector*); + EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;}; + +private: + double kEq(double, double); + RAbundVector getRAbundVector(SAbundVector*); + RAbundVector rdata; +}; + +/***********************************************************************/ + +#endif + + + + + diff --git a/globaldata.cpp b/globaldata.cpp index a0ccdcb..cc102a5 100644 --- a/globaldata.cpp +++ b/globaldata.cpp @@ -155,7 +155,7 @@ void GlobalData::parseGlobalData(string commandString, string optionText){ //input defaults for calculators if (commandName == "collect.single") { - if ((calc == "default") || (calc == "")) { calc = "sobs-chao-ace-jack-bootstrap-shannon-npshannon-simpson"; } + if ((calc == "default") || (calc == "")) { calc = "sobs-chao-ace-jack-bootstrap-shannon-npshannon-simpson-geom-qstat-logsd-bergerparker-bstick"; } Estimators.clear(); splitAtDash(calc, Estimators); } @@ -165,17 +165,17 @@ void GlobalData::parseGlobalData(string commandString, string optionText){ splitAtDash(calc, Estimators); } if (commandName == "collect.shared") { - if ((calc == "default") || (calc == "")) { calc = "sharedsobs-sharedchao-sharedace-sharedjabund-sharedsorensonabund-sharedjclass-sharedsorclass-sharedjest-sharedsorest-sharedthetayc-sharedthetan"; } + if ((calc == "default") || (calc == "")) { calc = "sharedsobs-sharedchao-sharedace-sharedjabund-sharedsorensonabund-sharedjclass-sharedsorclass-sharedjest-sharedsorest-sharedthetayc-sharedthetan-sharedkstest-sharedbdiversity"; } Estimators.clear(); splitAtDash(calc, Estimators); } if (commandName == "summary.single") { - if ((calc == "default") || (calc == "")) { calc = "sobs-chao-ace-jack-bootstrap-shannon-npshannon-simpson"; } + if ((calc == "default") || (calc == "")) { calc = "sobs-chao-ace-jack-bootstrap-shannon-npshannon-simpson-geom-logsd-qstat-bergerparker-bstick"; } Estimators.clear(); splitAtDash(calc, Estimators); } if (commandName == "summary.shared") { - if ((calc == "default") || (calc == "")) { calc = "sharedsobs-sharedchao-sharedace-sharedjabund-sharedsorensonabund-sharedjclass-sharedsorclass-sharedjest-sharedsorest-sharedthetayc-sharedthetan"; } + if ((calc == "default") || (calc == "")) { calc = "sharedsobs-sharedchao-sharedace-sharedjabund-sharedsorensonabund-sharedjclass-sharedsorclass-sharedjest-sharedsorest-sharedthetayc-sharedthetan-sharedkstest-sharedbdiversity"; } Estimators.clear(); splitAtDash(calc, Estimators); } diff --git a/logsd.cpp b/logsd.cpp new file mode 100644 index 0000000..94fa459 --- /dev/null +++ b/logsd.cpp @@ -0,0 +1,106 @@ +/* + * logsd.cpp + * Mothur + * + * Created by Thomas Ryabin on 2/23/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ + +#include "logsd.h" +#include "calculator.h" + +/***********************************************************************/ + +double LogSD::logS(double x){ + return -(1-x)*log(1-x)/x; +} +EstOutput LogSD::getValues(SAbundVector* rank){ + try { + + /*test data VVV + int dstring[] = {0,37,22,12,12,11,11,6,4,3,5,2,4,2,3,2,2,4,2,0,4,4,1,1,0,1,0,0,2,2,0,0,0,2,2,0,0,0,1,1,3,0,2,0,0,0,0,0,2,0,0,1,1,1,0,0,0,0,1,0,0,1,0,0,2,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,1,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1}; + vector dvec; + for(int i = 0; i < 1800; i++) + dvec.push_back(dstring[i]); + int mr = 1799; + int nb = 197; + int ns = 6815; + SAbundVector rankw = SAbundVector(dvec, mr,nb,ns); + SAbundVector *rank = &rankw;*/ + + data.resize(2,0); + int numInd = rank->getNumSeqs(); + int numSpec = rank->getNumBins(); + double snRatio = (double)numSpec/numInd; + double x = .5; + double step = .4999999999; + + while(fabs(snRatio - logS(x)) > .00001) //This uses a binary search to find the value of x. + { + if(logS(x) > snRatio) + x += step; + else + x -= step; + step /= 2; + } + double alpha = numInd*(1-x)/x; + + int oct = 1; + double octSumObs = 0; + double sumObs = 0; + double octSumExp = 0; + double sumExp = 0; + double maxDiff = 0; + for(int y = 1; y < rank->size(); y++) + { + if(y - .5 < pow(2.0, oct)) + { + octSumObs += rank->get(y); + octSumExp += alpha*pow(x,y)/(y); + } + else + { + sumObs += octSumObs; + octSumObs = rank->get(y); + + sumExp += octSumExp; + octSumExp = alpha*pow(x,y)/(y); + oct++; + } + + if(y == rank->size()-1) + { + sumObs += octSumObs; + sumExp += octSumExp; + } + + double diff = fabs(sumObs - .5 - sumExp); + if(diff > maxDiff) + maxDiff = diff; + } + + double DStatistic = (maxDiff + .5)/numSpec; + /*cout << "LogSD:\n"; + cout << "D Test Statistic = " << DStatistic << "\n"; + cout << ".05 confidence value = " << .89196/sqrt(numSpec) << "\n"; + cout << "If D Test Statistic is greater than the critical value then the data fits the Log Series Distribution model w/ 95% confidence.\n\n";*/ + + data[0] = DStatistic; + data[1] = .89196/sqrt(numSpec); + if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } + if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; } + return data; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + +/***********************************************************************/ + diff --git a/logsd.h b/logsd.h new file mode 100644 index 0000000..b96115b --- /dev/null +++ b/logsd.h @@ -0,0 +1,38 @@ +#ifndef LOGSD_H +#define LOGSD_H + +/* + * logsd.h + * Mothur + * + * Created by Thomas Ryabin on 2/23/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ +#include "calculator.h" + +/*This class implements the LogSD estimator on single group. +It is a child of the calculator class.*/ + +/***********************************************************************/ + +class LogSD : public Calculator { + +public: + LogSD() : Calculator("logsd", 3) {}; + EstOutput getValues(SAbundVector*); + EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;}; + +private: + double logS(double); + RAbundVector rdata; +}; + +/***********************************************************************/ + +#endif + + + + + diff --git a/qstat.cpp b/qstat.cpp new file mode 100644 index 0000000..591c36a --- /dev/null +++ b/qstat.cpp @@ -0,0 +1,76 @@ +/* + * qstat.cpp + * Mothur + * + * Created by Thomas Ryabin on 3/4/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ + +#include "qstat.h" +#include "calculator.h" + +/***********************************************************************/ + +EstOutput QStat::getValues(SAbundVector* rank){ + try { + + /*test data VVV + int dstring[] = {0,0,1,4,2,0,2,1,1,1,1,1,0,1,1,2,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,2,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1}; + vector dvec; + for(int i = 0; i < 171; i++) + dvec.push_back(dstring[i]); + int mr = 170; + int nb = 29; + int ns = 884; + SAbundVector rankw = SAbundVector(dvec, mr,nb,ns); + SAbundVector *rank = &rankw;*/ + data.resize(1,0); + int numSpec = rank->getNumBins(); + int r1 = -1; + int r3 = -1; + int r1Ind = 0; + int r3Ind = 0; + int sumSpec = 0; + int iqSum = 0; + for(int i = 1; i < rank->size(); i++) + { + if(r1 != -1 && r3 != -1) + i = rank->size(); + + sumSpec += rank->get(i); + + if(r1 == -1 && sumSpec >= numSpec*.25) + { + r1 = rank->get(i); + r1Ind = i; + } + else if(r3 == -1 && sumSpec >= numSpec*.75) + { + r3 = rank->get(i); + r3Ind = i; + } + else if(sumSpec >= numSpec*.25 && sumSpec < numSpec*.75) + iqSum += rank->get(i); + } + + double qstat = (.5*r1 + iqSum + .5*r3)/log((double)r3Ind/r1Ind); + //cout << "QStat:\nQStatistic = " << qstat << "\n\n"; + + data[0] = qstat; + if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } + + return data; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + +/***********************************************************************/ + diff --git a/qstat.h b/qstat.h new file mode 100644 index 0000000..3911dc1 --- /dev/null +++ b/qstat.h @@ -0,0 +1,32 @@ +#ifndef QSTAT_H +#define QSTAT_H +/* + * qstat.h + * Mothur + * + * Created by Thomas Ryabin on 3/4/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ +#include "calculator.h" + +/*This class implements the LogSD estimator on single group. +It is a child of the calculator class.*/ + +/***********************************************************************/ + +class QStat : public Calculator { + +public: + QStat() : Calculator("qstat", 3) {}; + EstOutput getValues(SAbundVector*); + EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*) {return data;}; + +private: + RAbundVector rdata; +}; + +/***********************************************************************/ + +#endif + diff --git a/rabundvector.cpp b/rabundvector.cpp index d6b7d12..6b0a048 100644 --- a/rabundvector.cpp +++ b/rabundvector.cpp @@ -13,6 +13,7 @@ using namespace std; #include "rabundvector.hpp" #include "sabundvector.hpp" #include "ordervector.hpp" +#include "calculator.h" /***********************************************************************/ @@ -52,6 +53,26 @@ RAbundVector::RAbundVector(string id, vector rav) : DataVector(id), data(ra } } +/***********************************************************************/ + +RAbundVector::RAbundVector(vector rav, int mr, int nb, int ns) { + try { + numBins = nb; + maxRank = mr; + numSeqs = ns; + data = rav; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the RAbundVector class Function RAbundVector. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the RAbundVector class function RAbundVector. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + + /***********************************************************************/ @@ -160,6 +181,33 @@ int RAbundVector::size(){ /***********************************************************************/ +void RAbundVector::quicksort(){ + sort(data.rbegin(), data.rend()); +} + +/***********************************************************************/ + +int RAbundVector::sum(){ + VecCalc vecCalc; + return vecCalc.sumElements(data); +} + +/***********************************************************************/ + +int RAbundVector::sum(int index){ + VecCalc vecCalc; + return vecCalc.sumElements(data, index); +} + +/***********************************************************************/ + +int RAbundVector::numNZ(){ + VecCalc vecCalc; + return vecCalc.numNZ(data); +} + +/***********************************************************************/ + vector::reverse_iterator RAbundVector::rbegin(){ return data.rbegin(); } diff --git a/rabundvector.hpp b/rabundvector.hpp index 327b29b..2a866bc 100644 --- a/rabundvector.hpp +++ b/rabundvector.hpp @@ -17,6 +17,7 @@ class RAbundVector : public DataVector { public: RAbundVector(); RAbundVector(int); + RAbundVector(vector, int, int, int); // RAbundVector(const RAbundVector&); RAbundVector(string, vector); RAbundVector(const RAbundVector& bv) : DataVector(bv), data(bv.data), maxRank(bv.maxRank), numBins(bv.numBins), numSeqs(bv.numSeqs){}; @@ -33,6 +34,10 @@ public: void pop_back(); void resize(int); int size(); + void quicksort(); + int sum(); + int sum(int); + int numNZ(); vector::reverse_iterator rbegin(); vector::reverse_iterator rend(); diff --git a/sabundvector.cpp b/sabundvector.cpp index a785619..b293132 100644 --- a/sabundvector.cpp +++ b/sabundvector.cpp @@ -10,6 +10,7 @@ using namespace std; #include "sabundvector.hpp" #include "utilities.hpp" +#include "calculator.h" /***********************************************************************/ @@ -23,7 +24,6 @@ SAbundVector::SAbundVector(int size) : DataVector(), data(size, 0), maxRank(0), SAbundVector::SAbundVector(string id, vector sav) : DataVector(id), data(sav) { try { - for(int i=0;i sav) : DataVector(id), data(sa /***********************************************************************/ +SAbundVector::SAbundVector(vector dataVec, int mr, int nb, int ns) { + try { + data = dataVec; + maxRank = mr; + numBins = nb; + numSeqs = ns; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the SAbundVector class Function SAbundVector. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the SAbundVector class function SAbundVector. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} +/***********************************************************************/ + SAbundVector::SAbundVector(ifstream& f): DataVector(), maxRank(0), numBins(0), numSeqs(0) { try { int hold; @@ -65,6 +83,7 @@ SAbundVector::SAbundVector(ifstream& f): DataVector(), maxRank(0), numBins(0), n } } + /***********************************************************************/ void SAbundVector::set(int sabund, int abundance){ @@ -91,6 +110,7 @@ void SAbundVector::set(int sabund, int abundance){ } } + /***********************************************************************/ int SAbundVector::get(int index){ @@ -118,6 +138,18 @@ void SAbundVector::push_back(int abundance){ exit(1); } } +/***********************************************************************/ + +void SAbundVector::quicksort(){ + sort(data.rbegin(), data.rend()); +} + +/***********************************************************************/ + +int SAbundVector::sum(){ + VecCalc vecCalc; + return vecCalc.sumElements(data); +} /***********************************************************************/ diff --git a/sabundvector.hpp b/sabundvector.hpp index 3c68462..bd657aa 100644 --- a/sabundvector.hpp +++ b/sabundvector.hpp @@ -6,6 +6,7 @@ using namespace std; #include "datavector.hpp" #include "rabundvector.hpp" #include "ordervector.hpp" +#include "calculator.h" /* This class is a child to datavector. It represents OTU information at a certain distance. @@ -24,6 +25,7 @@ public: SAbundVector(); SAbundVector(int); // SAbundVector(const SAbundVector&); + SAbundVector(vector, int, int, int); SAbundVector(string, vector); SAbundVector(const SAbundVector& rv) : DataVector(rv.label), data(rv.data), maxRank(rv.maxRank), numBins(rv.numBins), numSeqs(rv.numSeqs){}; SAbundVector(ifstream&); @@ -36,13 +38,15 @@ public: void set(int, int); int get(int); void push_back(int); + void quicksort(); + int sum(); void resize(int); int size(); void print(ostream&); void print(string, ostream&); - RAbundVector getRAbundVector(); + RAbundVector getRAbundVector(); SAbundVector getSAbundVector(); OrderVector getOrderVector(map*); @@ -50,6 +54,7 @@ private: vector data; // bool needToUpdate; // void updateStats(); + int maxRank; int numBins; int numSeqs; diff --git a/sharedbdiversity.cpp b/sharedbdiversity.cpp new file mode 100644 index 0000000..d45a11f --- /dev/null +++ b/sharedbdiversity.cpp @@ -0,0 +1,126 @@ +/* + * bdiversity.cpp + * Mothur + * + * Created by Thomas Ryabin on 3/13/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ + +#include "sharedbdiversity.h" +#include "calculator.h" + +/***********************************************************************/ + +double SharedBDiversity::whitt(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ + double nz1 = (double)shared1->numNZ(); + double nz2 = (double)shared2->numNZ(); + double sum = nz1; + for(int i = 1; i < shared1->size(); i++) + { + if(shared1->get(i).abundance == 0 && shared2->get(i).abundance != 0) + sum++; + } + return 2*sum/(nz1+nz2)-1; +} + +/***********************************************************************/ + +double SharedBDiversity::ms(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ + double a = 0; + double b = 0; + double c = 0; + for(int i = 1; i < shared1->size(); i++) + { + int abund1 = shared1->get(i).abundance; + int abund2 = shared2->get(i).abundance; + + if(abund1 > 0 && abund2 > 0) + a++; + else if(abund1 > 0 && abund2 == 0) + b++; + else if(abund1 == 0 && abund2 > 0) + c++; + } + return (b+c)/(a+b+c); +} + +/***********************************************************************/ + +double SharedBDiversity::sor(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ + double sum = 0; + double asum = 0; + double bsum = 0; + for(int i = 1; i < shared1->size(); i++) + { + int abund1 = shared1->get(i).abundance; + int abund2 = shared2->get(i).abundance; + + asum += abund1; + bsum += abund2; + if(abund1 >= abund2) + sum += abund2; + else + sum += abund1; + } + return 2*sum/(asum+bsum); +} + +/***********************************************************************/ + +double SharedBDiversity::mor(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ + double multSum = 0; + double powSum1 = 0; + double powSum2 = 0; + double ind1 = 0; + double ind2 = 0; + for(int i = 1; i < shared1->size(); i++) + { + double abund1 = (double)shared1->get(i).abundance; + double abund2 = (double)shared2->get(i).abundance; + multSum += abund1*abund2; + powSum1 += pow(abund1, 2); + powSum2 += pow(abund2, 2); + ind1 += abund1; + ind2 += abund2; + } + return 2*multSum / ((powSum1/pow(ind1, 2) + powSum2/pow(ind2, 2)) * (ind1*ind2)); +} +/***********************************************************************/ + +EstOutput SharedBDiversity::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ + try { + data.resize(4,0); + + double whitt1 = whitt(shared1, shared2); + double jac1 = 1-ms(shared1, shared2); + double sor1 = sor(shared1, shared2); + double mor1 = mor(shared1, shared2); + /*cout << "Whittaker's Measure: " << whitt1 << "\n"; + cout << "Marczewski-Steinhaus distance: " << jac1 << "\n"; + cout << "Sorensen index: " << sor1 << "\n"; + cout << "Morisita-Horn index: " << mor1 << "\n";*/ + + data[0] = whitt1; + data[1] = jac1; + data[2] = sor1; + data[3] = mor1; + if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } + if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; } + if (isnan(data[2]) || isinf(data[0])) { data[2] = 0; } + if (isnan(data[3]) || isinf(data[1])) { data[3] = 0; } + + return data; + } + + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + +/***********************************************************************/ \ No newline at end of file diff --git a/sharedbdiversity.h b/sharedbdiversity.h new file mode 100644 index 0000000..b32d6d8 --- /dev/null +++ b/sharedbdiversity.h @@ -0,0 +1,33 @@ +#ifndef SHAREDBDIVERSITY_H +#define SHAREDBDIVERSITY_H +/* + * bdiversity.h + * Mothur + * + * Created by Thomas Ryabin on 3/13/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ +#include "calculator.h" + +/*This class implements the BDiversity estimator on 2 groups. +It is a child of the calculator class.*/ + +/***********************************************************************/ + +class SharedBDiversity : public Calculator { + +public: + SharedBDiversity() : Calculator("sharedbdiversity", 3) {}; + EstOutput getValues(SAbundVector*) {return data;}; + EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*); +private: + double whitt(SharedRAbundVector*, SharedRAbundVector*); + double ms(SharedRAbundVector*, SharedRAbundVector*); + double sor(SharedRAbundVector*, SharedRAbundVector*); + double mor(SharedRAbundVector*, SharedRAbundVector*); +}; + +/***********************************************************************/ + +#endif \ No newline at end of file diff --git a/sharedkstest.cpp b/sharedkstest.cpp new file mode 100644 index 0000000..214c4d8 --- /dev/null +++ b/sharedkstest.cpp @@ -0,0 +1,79 @@ +/* + * kstest.cpp + * Mothur + * + * Created by Thomas Ryabin on 3/6/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ + +#include "sharedkstest.h" +#include "calculator.h" +#include + +/***********************************************************************/ + +EstOutput SharedKSTest::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2){ + try { + data.resize(2,0); + + //Must return shared1 and shared2 to original order at conclusion of kstest + vector initData1 = shared1->getData(); + vector initData2 = shared2->getData(); + shared1->sortD(); + shared2->sortD(); + + int numNZ1 = shared1->numNZ(); + int numNZ2 = shared2->numNZ(); + double numInd1 = (double)shared1->getNumSeqs(); + double numInd2 = (double)shared2->getNumSeqs(); + + double maxDiff = -1; + double sum1 = 0; + double sum2 = 0; + for(int i = 1; i < shared1->getNumBins(); i++) + { + sum1 += shared1->get(i).abundance; + sum2 += shared2->get(i).abundance; + double diff = fabs((double)sum1/numInd1 - (double)sum2/numInd2); + if(diff > maxDiff) + maxDiff = diff; + } + + double DStatistic = maxDiff*numNZ1*numNZ2; + double a = pow((double)(numNZ1 + numNZ2)/(numNZ1*numNZ2),.5); + double pVal = exp(-2*pow(maxDiff/a,2)); + double critVal = 1.36*a*numNZ1*numNZ2; + + /*cout << "Kolmogorov-Smirnov 2-sample test:\n"; + if(numNZ1 > 25 && numNZ2 > 25) //If the sample sizes are both bigger than 25. + cout << "P-Value = " << pVal << "\nP-Value is the probability that the data sets are significantly different.\n"; + else + { + //cout << "90% Confidence Critical Value = " << 1.22*a*numNZ1*numNZ2 << "\n"; + cout << "D-Statistic = " << DStatistic << "\n"; + cout << "95% Confidence Critical Value = " << critVal << "\n"; + cout << "If D-Statistic is greater than the critical value then the data sets are significantly different at the 95% confidence level.\n\n"; + }*/ + + shared1->setData(initData1); + shared2->setData(initData2); + + data[0] = DStatistic; + data[1] = critVal; + if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } + if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; } + + return data; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + +/***********************************************************************/ diff --git a/sharedkstest.h b/sharedkstest.h new file mode 100644 index 0000000..dfe1c1a --- /dev/null +++ b/sharedkstest.h @@ -0,0 +1,29 @@ +#ifndef SHAREDKSTEST_H +#define SHAREDKSTEST_H +/* + * kstest.h + * Mothur + * + * Created by Thomas Ryabin on 3/6/09. + * Copyright 2009 __MyCompanyName__. All rights reserved. + * + */ +#include "calculator.h" + +/*This class implements the KSTest estimator on 2 groups. +It is a child of the calculator class.*/ + +/***********************************************************************/ + +class SharedKSTest : public Calculator { + +public: + SharedKSTest() : Calculator("sharedkstest", 3) {}; + EstOutput getValues(SAbundVector*) {return data;}; + EstOutput getValues(SharedRAbundVector*, SharedRAbundVector*); +private: +}; + +/***********************************************************************/ + +#endif diff --git a/sharedordervector.h b/sharedordervector.h index 94d91c5..a530ad6 100644 --- a/sharedordervector.h +++ b/sharedordervector.h @@ -23,6 +23,9 @@ struct individual { string group; int bin; int abundance; + bool operator()(const individual& i1, const individual& i2) { + return (i1.abundance > i2.abundance); + } }; #include "sabundvector.hpp" diff --git a/sharedrabundvector.cpp b/sharedrabundvector.cpp index b35b06d..846f9cd 100644 --- a/sharedrabundvector.cpp +++ b/sharedrabundvector.cpp @@ -14,6 +14,7 @@ using namespace std; #include "utilities.hpp" #include "sabundvector.hpp" #include "ordervector.hpp" +#include /***********************************************************************/ @@ -119,6 +120,11 @@ void SharedRAbundVector::set(int binNumber, int newBinSize, string groupname){ exit(1); } } +/***********************************************************************/ + +void SharedRAbundVector::setData(vector newData){ + data = newData; +} /***********************************************************************/ @@ -126,14 +132,32 @@ int SharedRAbundVector::getAbundance(int index){ return data[index].abundance; } +/***********************************************************************/ + +int SharedRAbundVector::numNZ(){ + int sum = 0; + for(int i = 1; i < numBins; i++) + if(data[i].abundance > 0) + sum++; + return sum; +} +/***********************************************************************/ +void SharedRAbundVector::sortD(){ + struct individual indObj; + sort(data.begin()+1, data.end(), indObj); +} /***********************************************************************/ individual SharedRAbundVector::get(int index){ return data[index]; } +/***********************************************************************/ +vector SharedRAbundVector::getData(){ + return data; +} /***********************************************************************/ void SharedRAbundVector::push_back(int binSize, int otu, string groupName){ @@ -241,7 +265,6 @@ int SharedRAbundVector::getNumSeqs(){ int SharedRAbundVector::getMaxRank(){ return maxRank; } - /***********************************************************************/ SharedRAbundVector SharedRAbundVector::getSharedRAbundVector(){ @@ -268,7 +291,25 @@ RAbundVector SharedRAbundVector::getRAbundVector() { exit(1); } } +/***********************************************************************/ +RAbundVector SharedRAbundVector::getRAbundVector2() { + try { + RAbundVector rav; + for(int i = 0; i < numBins; i++) + if(data[i].abundance != 0) + rav.push_back(data[i].abundance-1); + return rav; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the SharedRAbundVector class Function getRAbundVector. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the SharedRAbundVector class function getRAbundVector. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} /***********************************************************************/ SharedSAbundVector SharedRAbundVector::getSharedSAbundVector(){ diff --git a/sharedrabundvector.h b/sharedrabundvector.h index 41a5798..92f71cd 100644 --- a/sharedrabundvector.h +++ b/sharedrabundvector.h @@ -30,7 +30,7 @@ public: SharedRAbundVector(int); //SharedRAbundVector(string, vector); SharedRAbundVector(const SharedRAbundVector& bv) : DataVector(bv), data(bv.data), maxRank(bv.maxRank), numBins(bv.numBins), numSeqs(bv.numSeqs){}; -// SharedRAbundVector(ifstream&); + //SharedRAbundVector(ifstream&); ~SharedRAbundVector(); int getNumBins(); @@ -42,8 +42,12 @@ public: void setGroupIndex(int); void set(int, int, string); //OTU, abundance, groupname + void setData(vector ); individual get(int); + vector getData(); int getAbundance(int); + int numNZ(); + void sortD(); //Sorts the data in descending order. void push_back(int, int, string); //abundance, OTU, groupname void pop_back(); void resize(int); @@ -54,6 +58,7 @@ public: void print(ostream&); RAbundVector getRAbundVector(); + RAbundVector getRAbundVector2(); SAbundVector getSAbundVector(); OrderVector getOrderVector(map*); SharedOrderVector getSharedOrderVector(); diff --git a/summarycommand.cpp b/summarycommand.cpp index 5adf44f..8c8a2aa 100644 --- a/summarycommand.cpp +++ b/summarycommand.cpp @@ -17,6 +17,11 @@ #include "npshannon.h" #include "shannon.h" #include "jackknife.h" +#include "geom.h" +#include "logsd.h" +#include "qstat.h" +#include "bergerparker.h" +#include "bstick.h" //********************************************************************************************************************** @@ -28,10 +33,21 @@ SummaryCommand::SummaryCommand(){ for (i=0; iEstimators.size(); i++) { if (validCalculator->isValidCalculator("summary", globaldata->Estimators[i]) == true) { + if(globaldata->Estimators[i] == "sobs"){ sumCalculators.push_back(new Sobs()); }else if(globaldata->Estimators[i] == "chao"){ sumCalculators.push_back(new Chao1()); + }else if(globaldata->Estimators[i] == "geom"){ + sumCalculators.push_back(new Geom()); + }else if(globaldata->Estimators[i] == "logsd"){ + sumCalculators.push_back(new LogSD()); + }else if(globaldata->Estimators[i] == "qstat"){ + sumCalculators.push_back(new QStat()); + }else if(globaldata->Estimators[i] == "bergerparker"){ + sumCalculators.push_back(new BergerParker()); + }else if(globaldata->Estimators[i] == "bstick"){ + sumCalculators.push_back(new BStick()); }else if(globaldata->Estimators[i] == "ace"){ convert(globaldata->getAbund(), abund); if(abund < 5) @@ -107,16 +123,16 @@ int SummaryCommand::execute(){ if(globaldata->allLines == 1 || globaldata->lines.count(count) == 1 || globaldata->labels.count(sabund->getLabel()) == 1){ cout << sabund->getLabel() << '\t' << count << endl; - + outputFileHandle << sabund->getLabel(); for(int i=0;i data = sumCalculators[i]->getValues(sabund); outputFileHandle << '\t'; - sumCalculators[i]->print(outputFileHandle); + //sumCalculators[i]->print(outputFileHandle); } + outputFileHandle << endl; } - sabund = input->getSAbundVector(); count++; } diff --git a/summarysharedcommand.cpp b/summarysharedcommand.cpp index 44bd581..a7c38f7 100644 --- a/summarysharedcommand.cpp +++ b/summarysharedcommand.cpp @@ -20,6 +20,8 @@ #include "sharedsorest.h" #include "sharedthetayc.h" #include "sharedthetan.h" +#include "sharedkstest.h" +#include "sharedbdiversity.h" #include "sharedochiai.h" #include "sharedanderberg.h" #include "sharedkulczynski.h" @@ -64,6 +66,8 @@ SummarySharedCommand::SummarySharedCommand(){ sumCalculators.push_back(new SharedThetaYC()); }else if (globaldata->Estimators[i] == "sharedthetan") { sumCalculators.push_back(new SharedThetaN()); + }else if (globaldata->Estimators[i] == "sharedkstest") { + sumCalculators.push_back(new SharedKSTest()); }else if (globaldata->Estimators[i] == "sharednseqs") { sumCalculators.push_back(new SharedNSeqs()); }else if (globaldata->Estimators[i] == "sharedochiai") { @@ -81,9 +85,12 @@ SummarySharedCommand::SummarySharedCommand(){ }else if (globaldata->Estimators[i] == "sharedbraycurtis") { sumCalculators.push_back(new SharedBrayCurtis()); } + else if (globaldata->Estimators[i] == "sharedbdiversity") { + sumCalculators.push_back(new SharedBDiversity()); + } + } } - //reset calc for next command globaldata->setCalc(""); @@ -152,6 +159,7 @@ int SummarySharedCommand::execute(){ int n = 1; for (int k = 0; k < (lookup.size() - 1); k++) { // pass cdd each set of groups to commpare for (int l = n; l < lookup.size(); l++) { + outputFileHandle << order->getLabel() << '\t' << (lookup[k]->getGroup() + lookup[l]->getGroup()) << '\t' << '\t'; //print out label and group outputFileHandle << order->getLabel() << '\t'; //sort groups to be alphanumeric diff --git a/validcalculator.cpp b/validcalculator.cpp index 0ffa216..4c01445 100644 --- a/validcalculator.cpp +++ b/validcalculator.cpp @@ -126,16 +126,22 @@ bool ValidCalculators::isValidCalculator(string parameter, string calculator) { void ValidCalculators::initialSingle() { try { - single["sobs"] = "sobs"; - single["chao"] = "chao"; - single["ace"] = "ace"; - single["jack"] = "jack"; - single["shannon"] = "shannon"; - single["npshannon"] = "npshannon"; - single["simpson"] = "simpson"; - single["bootstrap"] = "bootstrap"; - single["default"] = "default"; + single["sobs"] = "sobs"; + single["chao"] = "chao"; + single["ace"] = "ace"; + single["jack"] = "jack"; + single["shannon"] = "shannon"; + single["npshannon"] = "npshannon"; + single["simpson"] = "simpson"; + single["bergerparker"] = "bergerparker"; + single["bootstrap"] = "bootstrap"; + single["geom"] = "geom"; + single["logsd"] = "logsd"; + single["qstat"] = "qstat"; + single["bstick"] = "bstick"; single["nseqs"] = "nseqs"; + single["default"] = "default"; + } catch(exception& e) { cout << "Standard Error: " << e.what() << " has occurred in the ValidCalculator class Function initialSingle. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; @@ -161,6 +167,8 @@ void ValidCalculators::initialShared() { shared["sharedsorest"] = "sharedsorest"; shared["sharedthetayc"] = "sharedthetayc"; shared["sharedthetan"] = "sharedthetan"; + shared["sharedkstest"] = "sharedkstest"; + shared["sharedbdiversity"] = "sharedbdiversity"; shared["sharednseqs"] = "sharednseqs"; shared["sharedochiai"] = "sharedochiai"; shared["sharedanderberg"] = "sharedanderberg"; @@ -216,7 +224,12 @@ void ValidCalculators::initialSummary() { summary["shannon"] = "shannon"; summary["npshannon"] = "npshannon"; summary["simpson"] = "simpson"; + summary["bergerparker"] = "bergerparker"; + summary["geom"] = "geom"; summary["bootstrap"] = "bootstrap"; + summary["logsd"] = "logsd"; + summary["qstat"] = "qstat"; + summary["bstick"] = "bstick"; summary["nseqs"] = "nseqs"; summary["default"] = "default"; } @@ -244,6 +257,8 @@ void ValidCalculators::initialSharedSummary() { sharedsummary["sharedsorest"] = "sharedsorest"; sharedsummary["sharedthetayc"] = "sharedthetayc"; sharedsummary["sharedthetan"] = "sharedthetan"; + sharedsummary["sharedkstest"] = "sharedkstest"; + sharedsummary["sharedbdiversity"] = "sharedbdiversity"; sharedsummary["sharednseqs"] = "sharednseqs"; sharedsummary["sharedochiai"] = "sharedochiai"; sharedsummary["sharedanderberg"] = "sharedanderberg"; diff --git a/validparameter.cpp b/validparameter.cpp index 6d325f7..50610f0 100644 --- a/validparameter.cpp +++ b/validparameter.cpp @@ -13,21 +13,38 @@ ValidParameters::ValidParameters() { try { - initialReaddist(); - initialReadotu(); - initialReadtree(); - initialCluster(); - initialDeconvolute(); - initialParsimony(); - initialCollectsingle(); - initialCollectshared(); - initialRarefactsingle(); - initialRarefactshared(); - initialSummarysingle(); - initialSummaryshared(); - initialUnifracweighted(); - initialUnifracunweighted(); - initialLibshuff(); + parameters["phylip"] = "phylip"; + parameters["column"] = "column"; + parameters["list"] = "list"; + parameters["rabund"] = "rabund"; + parameters["sabund"] = "sabund"; + parameters["name"] = "name"; + parameters["group"] = "group"; + parameters["order"] = "order"; + parameters["fasta"] = "fasta"; + parameters["tree"] = "tree"; + parameters["fileroot"] = "fileroot"; + parameters["cutoff"] = "cutoff"; + parameters["method"] = "method"; + parameters["format"] = "format"; + parameters["precision"] = "precision"; + parameters["label"] = "label"; + parameters["line"] = "line"; + parameters["iters"] = "iters"; + parameters["jumble"] = "jumble"; + parameters["freq"] = "freq"; + parameters["abund"] = "abund"; + parameters["random"] = "random"; + parameters["groups"] = "groups"; + parameters["calc"] = "calc"; + parameters["sharedrarefaction"] = "sharedrarefaction"; + parameters["sharedsummary"] = "sharedsummary"; + parameters["shared"] = "shared"; + parameters["single"] = "single"; + parameters["rarefaction"] = "rarefaction"; + + initCommandParameters(); + initParameterRanges(); } catch(exception& e) { cout << "Standard Error: " << e.what() << " has occurred in the ValidParameters class Function ValidParameters. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; @@ -43,198 +60,180 @@ ValidParameters::ValidParameters() { ValidParameters::~ValidParameters() {} + /***********************************************************************/ -bool ValidParameters::isValidParameter(string parameter, string command) { +bool ValidParameters::isValidParameter(string parameter, string command, string value) { try { + bool valid = false; + vector cParams = commandParameters[command]; + int numParams = cParams.size(); + for(int i = 0; i < numParams; i++) + { + if(cParams.at(i).compare(parameter) == 0) + { + valid = true; + i = numParams; + } + } + if(!valid) + { + cout << "'" << parameter << "' is not a valid parameter for the " << command << " command.\n"; + return false; + } + + if(parameterRanges.count(parameter) != 1) + return true; - if (command == "read.dist") { - //is it valid - if ((readdist.find(parameter)) != (readdist.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = readdist.begin(); it != readdist.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "read.otu") { - //is it valid - if ((readotu.find(parameter)) != (readotu.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = readotu.begin(); it != readotu.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - //are you looking for a calculator for a rarefaction parameter - }else if (command == "read.tree") { - //is it valid - if ((readtree.find(parameter)) != (readtree.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = readtree.begin(); it != readtree.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - //are you looking for a calculator for a summary parameter - }else if (command == "cluster") { - //is it valid - if ((cluster.find(parameter)) != (cluster.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = cluster.begin(); it != cluster.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } //are you looking for a calculator for a sharedsummary parameter - }else if (command == "deconvolute") { - //is it valid - if ((deconvolute.find(parameter)) != (deconvolute.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = deconvolute.begin(); it != deconvolute.end(); it++) { - cout << it->first; - } - cout << endl; - return false; } - }else if (command == "parsimony") { - //is it valid - if ((parsimony.find(parameter)) != (parsimony.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = parsimony.begin(); it != parsimony.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "collect.single") { - //is it valid - if ((collectsingle.find(parameter)) != (collectsingle.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = collectsingle.begin(); it != collectsingle.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "collect.shared") { - //is it valid - if ((collectshared.find(parameter)) != (collectshared.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = collectshared.begin(); it != collectshared.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "rarefaction.single") { - //is it valid - if ((rarefactsingle.find(parameter)) != (rarefactsingle.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = rarefactsingle.begin(); it != rarefactsingle.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "rarefaction.shared") { - //is it valid - if ((rarefactshared.find(parameter)) != (rarefactshared.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = rarefactshared.begin(); it != rarefactshared.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "summary.single") { - //is it valid - if ((summarysingle.find(parameter)) != (summarysingle.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = summarysingle.begin(); it != summarysingle.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "summary.shared") { - //is it valid - if ((summaryshared.find(parameter)) != (summaryshared.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = summaryshared.begin(); it != summaryshared.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "unifrac.weighted") { - //is it valid - if ((unifracweighted.find(parameter)) != (unifracweighted.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = unifracweighted.begin(); it != unifracweighted.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "unifrac.unweighted") { - //is it valid - if ((unifracunweighted.find(parameter)) != (unifracunweighted.end())) { - return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = unifracunweighted.begin(); it != unifracunweighted.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "libshuff") { - //is it valid - if ((libshuff.find(parameter)) != (libshuff.end())) { + int pVal; + double piSentinel = 3.14159; + vector range = parameterRanges[parameter]; + + valid = convertTest(value, pVal); + + if(!valid) + return false; + + + + /******************************************************************************************************** + Special Cases + *********************************************************************************************************/ + + if(parameter.compare("precision") == 0) + { + double logNum = log10((double)pVal); + double diff = (double)((int)logNum - logNum); + if(diff != 0) + { + cout << "The precision parameter can only take powers of 10 as a value (e.g. 10,1000,1000, etc.)\n"; + return false; + } + } + + /************************************************************************************************************/ + + + + double a,b,c,d,e; + + if(range.at(1).compare("NA") == 0) + a = piSentinel; + else + a = atoi(range.at(1).c_str()); + + if(range.at(3).compare("NA") == 0) + b = piSentinel; + else + b = atoi(range.at(3).c_str()); + + if(range.at(4).compare("between") == 0) + c = 0; + else if(range.at(4).compare("only") == 0) + c = 1; + else + { + cout << "The range can only be 'between' or 'only' the bounding numbers.\n"; + return false; + } + + if(range.at(0).compare(">") == 0) + d = 0; + else if(range.at(0).compare(">=") == 0 || range[3].compare("=>") == 0) + d = 1; + else + { + cout << "The parameter value can only be '>', '>=', or '=>' the lower bounding number.\n"; + return false; + } + + if(range.at(2).compare("<") == 0) + e = 0; + else if(range.at(2).compare("<=") == 0 || range[4].compare("=<") == 0) + e = 1; + else + { + cout << "The parameter value can only be '<', '<=', or '=<' the upper bounding number.\n"; + return false; + } + + bool a0 = pVal > a; + bool a1 = pVal >= a; + bool b0 = pVal < b; + bool b1 = pVal <= b; + + if(c != 1) + { + if(a == piSentinel && b == piSentinel) return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = libshuff.begin(); it != libshuff.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - }else if (command == "heatmap") { - //is it valid - if ((heatmap.find(parameter)) != (heatmap.end())) { + if(a != piSentinel && b == piSentinel) + { + if(d == 0) + valid = a0; + else + valid = a1; + } + else if(a == piSentinel && b != piSentinel) + { + if(e == 0) + valid = b0; + else + valid = b1; + } + else + { + if(d == 0 && e == 0) + valid = (a0 && b0); + else if(d == 0 && e == 1) + valid = (a0 && b1); + else if(d == 1 && e == 0) + valid = (a1 && b0); + else + valid = (a1 && b1); + } + } + else + { + if(a == piSentinel && b == piSentinel) return true; - }else { - cout << parameter << " is not a valid parameter for the " + command + " command. Valid parameters are "; - for (it = heatmap.begin(); it != heatmap.end(); it++) { - cout << it->first << ", "; - } - cout << endl; - return false; } - - //not a valid paramter - }else if (command == "help") { cout << parameter << " is not a valid parameter for the " + command + " command. There are no vaild parameters." << endl; - }else if (command == "quit") { cout << parameter << " is not a valid parameter for the " + command + " command. There are no vaild parameters." << endl; - }else if (command == "get.group") { cout << parameter << " is not a valid parameter for the " + command + " command. There are no vaild parameters." << endl; - }else if (command == "get.label") { cout << parameter << " is not a valid parameter for the " + command + " command. There are no vaild parameters." << endl; - }else if (command == "get.line") { cout << parameter << " is not a valid parameter for the " + command + " command. There are no vaild parameters." << endl; } + if(a != piSentinel && b == piSentinel) + valid = (pVal == a); + else if(a == piSentinel && b != piSentinel) + valid = (pVal == b); + else + valid = (pVal == a || pVal == b); + } - return false; + if(valid) + return true; + else + { + cout << "The '" << parameter << "' parameter needs to be "; + if(c == 1) + cout << "either '" << a << "' or '" << b << "'.\n"; + else + { + if(a != piSentinel) + { + cout << ">"; + if(d != 0) + cout << "="; + cout << " '" << a << "'"; + } + if(b == piSentinel) + cout << ".\n"; + else if(a != piSentinel) + cout << " and "; + if(b != piSentinel) + { + cout << "<"; + if(e != 0) + cout << "="; + cout << " '" << b << ".\n"; + } + } + return false; + } } catch(exception& e) { cout << "Standard Error: " << e.what() << " has occurred in the ValidParameters class Function isValidParameter. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; @@ -247,113 +246,128 @@ bool ValidParameters::isValidParameter(string parameter, string command) { } /***********************************************************************/ -void ValidParameters::initialReaddist() { - readdist["phylip"] = "phylip"; - readdist["column"] = "column"; - readdist["name"] = "name"; - readdist["group"] = "group"; - readdist["cutoff"] = "cutoff"; - readdist["precision"] = "precision"; -} -/***********************************************************************/ -void ValidParameters::initialReadotu() { - readotu["list"] = "list"; - readotu["rabund"] = "rabund"; - readotu["sabund"] = "sabund"; - readotu["shared"] = "shared"; - readotu["group"] = "group"; - readotu["order"] = "order"; - readotu["label"] = "label"; - readotu["line"] = "line"; -} -/***********************************************************************/ -void ValidParameters::initialReadtree() { - readtree["group"] = "group"; - readtree["tree"] = "tree"; -} -/***********************************************************************/ -void ValidParameters::initialCluster() { - cluster["cutoff"] = "cutoff"; - cluster["method"] = "method"; - cluster["precision"] = "precision"; -} -/***********************************************************************/ -void ValidParameters::initialDeconvolute() { - deconvolute["fasta"] = "fasta"; -} -/***********************************************************************/ -void ValidParameters::initialParsimony() { - parsimony["iters"] = "iters"; - parsimony["random"] = "random"; - parsimony["groups"] = "groups"; -} -/***********************************************************************/ -void ValidParameters::initialCollectsingle() { - collectsingle["label"] = "label"; - collectsingle["line"] = "line"; - collectsingle["freq"] = "freq"; - collectsingle["calc"] = "calc"; -} -/***********************************************************************/ -void ValidParameters::initialCollectshared() { - collectshared["label"] = "label"; - collectshared["line"] = "line"; - collectshared["freq"] = "freq"; - collectshared["calc"] = "calc"; - collectshared["jumble"] = "jumble"; -} -/***********************************************************************/ -void ValidParameters::initialRarefactsingle() { - rarefactsingle["label"] = "label"; - rarefactsingle["line"] = "line"; - rarefactsingle["freq"] = "freq"; - rarefactsingle["calc"] = "calc"; - rarefactsingle["iters"] = "iters"; -} -/***********************************************************************/ -void ValidParameters::initialRarefactshared() { - rarefactshared["label"] = "label"; - rarefactshared["line"] = "line"; - rarefactshared["jumble"] = "jumble"; - rarefactshared["calc"] = "calc"; - rarefactshared["iters"] = "iters"; -} -/***********************************************************************/ -void ValidParameters::initialSummarysingle() { - summarysingle["label"] = "label"; - summarysingle["line"] = "line"; - summarysingle["calc"] = "calc"; -} + /***********************************************************************/ -void ValidParameters::initialSummaryshared() { - summaryshared["label"] = "label"; - summaryshared["line"] = "line"; - summaryshared["calc"] = "calc"; - summaryshared["jumble"] = "jumble"; +void ValidParameters::initCommandParameters() { + try { + //{"parameter1","parameter2",...,"last parameter"}; + + string readdistArray[] = {"phylip","name","cutoff","precision"}; + commandParameters["read.dist"] = addParameters(readdistArray, sizeof(readdistArray)/sizeof(string)); + + string readotuArray[] = {"list","order","group","shared", "sabund"}; + commandParameters["read.otu"] = addParameters(readotuArray, sizeof(readotuArray)/sizeof(string)); + + string clusterArray[] = {"cutoff","precision","method"}; + commandParameters["cluster"] = addParameters(clusterArray, sizeof(clusterArray)/sizeof(string)); + + string deconvoluteArray[] = {"fasta"}; + commandParameters["deconvolute"] = addParameters(deconvoluteArray, sizeof(deconvoluteArray)/sizeof(string)); + + string collectsingleArray[] = {"freq","line","label","single","precision","abund"}; + commandParameters["collect.single"] = addParameters(collectsingleArray, sizeof(collectsingleArray)/sizeof(string)); + + string collectsharedArray[] = {"jumble","freq","line","label","shared","groups"}; + commandParameters["collect.shared"] = addParameters(collectsharedArray, sizeof(collectsharedArray)/sizeof(string)); + + string getgroupArray[] = {}; + commandParameters["get.group"] = addParameters(getgroupArray, sizeof(getgroupArray)/sizeof(string)); + + string getlabelArray[] = {}; + commandParameters["get.label"] = addParameters(getlabelArray, sizeof(getlabelArray)/sizeof(string)); + + string getlineArray[] = {}; + commandParameters["get.line"] = addParameters(getlineArray, sizeof(getlineArray)/sizeof(string)); + string rarefactionsingleArray[] = {"iters","freq","line","label","rarefaction","abund"}; + commandParameters["rarefaction.single"] = addParameters(rarefactionsingleArray, sizeof(rarefactionsingleArray)/sizeof(string)); + + string rarefactionsharedArray[] = {"iters","jumble","line","label","sharedrarefaction"}; + commandParameters["rarefaction.shared"] = addParameters(rarefactionsharedArray, sizeof(rarefactionsharedArray)/sizeof(string)); + + string summarysingleArray[] = {"line","label","summary","abund"}; + commandParameters["summary.single"] = addParameters(summarysingleArray, sizeof(summarysingleArray)/sizeof(string)); + + string summarysharedArray[] = {"jumble","line","label","sharedsummary"}; + commandParameters["summary.shared"] = addParameters(summarysharedArray, sizeof(summarysharedArray)/sizeof(string)); + + string quitArray[] = {}; + commandParameters["quit"] = addParameters(quitArray, sizeof(quitArray)/sizeof(string)); + + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the ValidParameters class Function isValidParameter. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the ValidParameters class function isValidParameter. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } } + /***********************************************************************/ -void ValidParameters::initialUnifracweighted() { - unifracweighted["iters"] = "iters"; - unifracweighted["groups"] = "groups"; -} + /***********************************************************************/ -void ValidParameters::initialUnifracunweighted() { - unifracunweighted["iters"] = "iters"; - unifracunweighted["groups"] = "groups"; +void ValidParameters::initParameterRanges() { + try { + int rangeSize = 5; + + /************************************************************************************************************** + {">=" or "=>" or ">" if the value should be greater than or equal to or just greater than the lower bound, + A number representing the lower bound ("NA" if there is no lower bound), + "<=" or "=<" or "<" if the value shoud be less than or equal to or just less than the upper bound, + A number representing the upper bound ("NA" if there is no lower bound), + "between" if between lower and upper bounds or "only" if exactly one of the bounds}; + + # = parameter + # (>, >=) lower bound, # (<, <=) upperbound, # should be (between, only) lower and upper bounds. + ***********************************************************************************************************/ + + string precisionArray[] = {">=","10", "<","NA", "between"}; + parameterRanges["precision"] = addParameters(precisionArray, rangeSize); + + string itersArray[] = {">=","10", "<","NA", "between"}; + parameterRanges["iters"] = addParameters(itersArray, rangeSize); + + string jumbleArray[] = {">","0", "<","1", "only"}; + parameterRanges["jumble"] = addParameters(jumbleArray, rangeSize); + + string freqArray[] = {">","1", "<","NA", "between"}; + parameterRanges["freq"] = addParameters(freqArray, rangeSize); + + string lineArray[] = {">=","1", "<","NA", "between"}; + parameterRanges["line"] = addParameters(lineArray, rangeSize); + + string abundArray[] = {">=","5", "<","NA", "between"}; + parameterRanges["abund"] = addParameters(abundArray, rangeSize); + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the ValidParameters class Function isValidParameter. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the ValidParameters class function isValidParameter. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } } + /***********************************************************************/ -void ValidParameters::initialLibshuff() { - libshuff["cutoff"] = "cutoff"; - libshuff["iters"] = "iters"; - libshuff["groups"] = "groups"; - libshuff["step"] = "step"; - libshuff["form"] = "form"; -} + /***********************************************************************/ -void ValidParameters::initialHeatmap() { - heatmap["label"] = "label"; - heatmap["line"] = "line"; +vector ValidParameters::addParameters(string parameters[], int size) { + try { + vector pVector (parameters, parameters+size); + return pVector; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the ValidParameters class Function isValidParameter. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the ValidParameters class function isValidParameter. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } } /***********************************************************************/ + diff --git a/validparameter.h b/validparameter.h index b3d56e5..61b1c5b 100644 --- a/validparameter.h +++ b/validparameter.h @@ -12,6 +12,7 @@ using namespace std; #include "mothur.h" +#include "utilities.hpp" //This class contains a list of all valid parameters in Mothur. //It has a function which will tell you if your parameter is valid. @@ -23,8 +24,12 @@ class ValidParameters { public: ValidParameters(); ~ValidParameters(); - bool isValidParameter(string, string); - + bool isValidParameter(string); + bool isValidParameter(string, string, string); + vector addParameters(string[], int); + void initCommandParameters(); + void initParameterRanges(); + private: map readdist; map readotu; @@ -44,23 +49,8 @@ class ValidParameters { map heatmap; map::iterator it; - - void initialReaddist(); - void initialReadotu(); - void initialReadtree(); - void initialCluster(); - void initialDeconvolute(); - void initialParsimony(); - void initialCollectsingle(); - void initialCollectshared(); - void initialRarefactsingle(); - void initialRarefactshared(); - void initialSummarysingle(); - void initialSummaryshared(); - void initialUnifracweighted(); - void initialUnifracunweighted(); - void initialLibshuff(); - void initialHeatmap(); + map > commandParameters; + map > parameterRanges; }; -- 2.39.2