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