From: ryabin Date: Fri, 27 Mar 2009 19:25:25 +0000 (+0000) Subject: added the Calculators Thomas made in the fall. Added parameter and command error... X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=commitdiff_plain;h=eb1c88346fb246e95a6b38935b103f95e38b82ca added the Calculators Thomas made in the fall. Added parameter and command error checking. --- 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; };