sherman5 authored on 12/12/2017 22:30:31
Showing16 changed files

... ...
@@ -1,5 +1,5 @@
1 1
 Package: CoGAPS
2
-Version: 2.6.0
2
+Version: 2.5.10
3 3
 Date: 2014-08-23
4 4
 Title: Coordinated Gene Activity in Pattern Sets
5 5
 Author: Wai-shing Lee, Conor Kelton, Ondrej Maxian, Jacob Carey, Genevieve
... ...
@@ -118,7 +118,7 @@ gapsRun <- function(D, S, ABins = data.frame(), PBins = data.frame(),
118 118
   names(cogapResult)[13] = "meanChi2";
119 119
 
120 120
   if (messages)
121
-    message(paste("Chi-Squared of Mean:",calcChiSq))
121
+    message(paste("Chi-Squared of Mean:", calcChiSq))
122 122
 
123 123
   return(cogapResult);
124 124
 }
... ...
@@ -1,8 +1,12 @@
1 1
 library(CoGAPS)
2 2
 
3
+print(packageVersion('CoGAPS'))
4
+
3 5
 data(SimpSim)
4 6
 nIter <- 1000
5
-eat <- gapsRun(SimpSim.D, SimpSim.S, nFactor=3, messages=FALSE, seed=12345)
7
+nFactor <- 7
8
+eat <- gapsRun(D=SimpSim.D, S=SimpSim.S, nEquil=nIter, nSample=nIter,
9
+    nFactor=nFactor, seed=12345)
6 10
 
7 11
 #data(GIST_TS_20084)
8 12
 #nIter <- 3000
... ...
@@ -1,6 +1,7 @@
1 1
 #!/bin/bash
2 2
 
3 3
 rm -f callgrind.out.*
4
-R -d "valgrind --tool=callgrind" -f profile_standard.R > valgrind_out.txt
4
+rm -f vgcore.*
5
+R -d "valgrind --tool=callgrind" -f profile_standard.R
5 6
 kcachegrind callgrind.out.*
6 7
 
7 8
deleted file mode 100644
... ...
@@ -1,20 +0,0 @@
1
-
2
-R version 3.4.1 (2017-06-30) -- "Single Candle"
3
-Copyright (C) 2017 The R Foundation for Statistical Computing
4
-Platform: x86_64-pc-linux-gnu (64-bit)
5
-
6
-R is free software and comes with ABSOLUTELY NO WARRANTY.
7
-You are welcome to redistribute it under certain conditions.
8
-Type 'license()' or 'licence()' for distribution details.
9
-
10
-  Natural language support but running in an English locale
11
-
12
-R is a collaborative project with many contributors.
13
-Type 'contributors()' for more information and
14
-'citation()' on how to cite R or R packages in publications.
15
-
16
-Type 'demo()' for some demos, 'help()' for on-line help, or
17
-'help.start()' for an HTML browser interface to help.
18
-Type 'q()' to quit R.
19
-
20
-> library(CoGAPS)
... ...
@@ -41,7 +41,7 @@ uint64_t AtomicSupport::randomFreePosition() const
41 41
     do
42 42
     {
43 43
         pos = gaps::random::uniform64();
44
-    } while (mAtomicDomain.count(pos));
44
+    } while (mAtomicDomain.count(pos) > 0);
45 45
     return pos;
46 46
 }
47 47
 
... ...
@@ -56,15 +56,15 @@ uint64_t AtomicSupport::randomAtomPosition() const
56 56
 AtomicProposal AtomicSupport::proposeDeath() const
57 57
 {
58 58
     uint64_t location = randomAtomPosition();
59
-    uint64_t mass = mAtomicDomain.at(location);
59
+    double mass = mAtomicDomain.at(location);
60 60
     return AtomicProposal(mLabel, 'D', location, -mass);
61 61
 }
62 62
 
63 63
 AtomicProposal AtomicSupport::proposeBirth() const
64 64
 {
65 65
     uint64_t location = randomFreePosition();
66
-    uint64_t mass = gaps::random::exponential(mLambda);
67
-    return AtomicProposal(mLabel, 'B', location, mass);
66
+    double mass = gaps::random::exponential(mLambda);
67
+    return AtomicProposal(mLabel, 'B', location, std::max(mass, EPSILON));
68 68
 }
69 69
 
70 70
 // move atom between adjacent atoms
... ...
@@ -72,7 +72,7 @@ AtomicProposal AtomicSupport::proposeMove() const
72 72
 {
73 73
     // get random atom
74 74
     uint64_t location = randomAtomPosition();
75
-    uint64_t mass = mAtomicDomain.at(location);
75
+    double mass = mAtomicDomain.at(location);
76 76
 
77 77
     // get an iterator to this atom
78 78
     std::map<uint64_t, double>::const_iterator it;
... ...
@@ -130,7 +130,7 @@ void AtomicSupport::updateAtomMass(uint64_t pos, double delta)
130 130
         mTotalMass += delta;
131 131
         mNumAtoms++;
132 132
     }
133
-    if (mAtomicDomain.at(pos) <= EPSILON) // check if atom is too small
133
+    if (mAtomicDomain.at(pos) < EPSILON) // check if atom is too small
134 134
     {
135 135
         mTotalMass -= mAtomicDomain.at(pos);
136 136
         mAtomicDomain.erase(pos);
... ...
@@ -138,8 +138,6 @@ void AtomicSupport::updateAtomMass(uint64_t pos, double delta)
138 138
     }
139 139
 }
140 140
 
141
-#define GEOM_F(x) (((double) (x)) / ((double) (x + 1)))
142
-#define POISS_F(x,y) ((double)((x) * (y)) / (double)((x) - (y)))
143 141
 AtomicProposal AtomicSupport::makeProposal() const
144 142
 {
145 143
     if (mNumAtoms == 0)
... ...
@@ -148,7 +146,7 @@ AtomicProposal AtomicSupport::makeProposal() const
148 146
     }
149 147
 
150 148
     double unif = gaps::random::uniform();
151
-    if (mNumAtoms < 2 && unif <= 0.6667 || unif <= 0.5) // birth/death
149
+    if ((mNumAtoms < 2 && unif <= 0.6667) || unif <= 0.5) // birth/death
152 150
     {
153 151
         if (mNumAtoms >= mMaxNumAtoms)
154 152
         {
... ...
@@ -158,15 +156,20 @@ AtomicProposal AtomicSupport::makeProposal() const
158 156
         double pDelete = 0.5; // default uniform prior for mAlpha == 0
159 157
         if (mAlpha > 0) // poisson prior
160 158
         {
161
-            pDelete = 1 + POISS_F(mMaxNumAtoms, mNumAtoms) / (mAlpha * mNumBins);
159
+            double dMax = (double)mMaxNumAtoms;
160
+            double dNum = (double)mNumAtoms;
161
+            double maxTerm = (dMax - dNum) / dMax;
162
+
163
+            pDelete = dNum / (dNum + mAlpha * (double)mNumBins * maxTerm);
162 164
         }
163 165
         else if (mAlpha < 0) // geometric prior
164 166
         {
165
-            pDelete = 1 + GEOM_F(mNumAtoms) / GEOM_F(-mAlpha * mNumBins);
167
+            double c = -mAlpha * mNumBins / (-mAlpha * mNumBins + 1.0);
168
+            pDelete = (double)mNumAtoms / ((mNumAtoms + 1) * c + (double)mNumAtoms);
166 169
         }
167 170
         return gaps::random::uniform() < pDelete ? proposeDeath() : proposeBirth();
168 171
     }
169
-    return mNumAtoms < 2 || unif <= 0.75 ? proposeMove() : proposeExchange();
172
+    return (mNumAtoms < 2 || unif >= 0.75) ? proposeMove() : proposeExchange();
170 173
 }
171 174
 
172 175
 void AtomicSupport::acceptProposal(const AtomicProposal &prop)
... ...
@@ -100,7 +100,7 @@ public:
100 100
     // setters
101 101
     void setAlpha(double alpha) {mAlpha = alpha;}
102 102
     void setLambda(double lambda) {mLambda = lambda;}
103
-    void setMaxNumAtoms(uint64_t max) {mMaxNumAtoms = max;}
103
+    //void setMaxNumAtoms(uint64_t max) {mMaxNumAtoms = max;}
104 104
 };
105 105
 
106 106
 #endif
... ...
@@ -40,8 +40,6 @@ double maxGibbsMassP, int seed=-1, bool messages=false, bool singleCellRNASeq=fa
40 40
 
41 41
     chi2Vec.concat(chi2VecSample);
42 42
 
43
-    // compute statistics
44
-
45 43
     //Just leave the snapshots as empty lists
46 44
     return Rcpp::List::create(
47 45
         Rcpp::Named("Amean") = sampler.AMeanRMatrix(),
... ...
@@ -2,6 +2,7 @@
2 2
 #include "Algorithms.h"
3 3
 
4 4
 #include <Rcpp.h>
5
+#include <iostream>
5 6
 
6 7
 #define EPSILON 1E-10
7 8
 
... ...
@@ -25,6 +26,11 @@ mStatUpdates(0)
25 26
     mADomain.setLambda(alphaA * sqrt(nFactor / meanD));
26 27
     mPDomain.setAlpha(alphaP);
27 28
     mPDomain.setLambda(alphaP * sqrt(nFactor / meanD));
29
+
30
+    std::cout << "A lambda: " << mADomain.lambda() << "\n";
31
+    std::cout << "A alpha: " << mADomain.alpha() << "\n";
32
+    std::cout << "P lambda: " << mPDomain.lambda() << "\n";
33
+    std::cout << "P alpha: " << mPDomain.alpha() << "\n";
28 34
 }
29 35
 
30 36
 double GibbsSampler::getGibbsMass(const MatrixChange &change)
... ...
@@ -130,7 +136,7 @@ const AtomicProposal &proposal, double rejectProb)
130 136
 {
131 137
     MatrixChange change = domain.getMatrixChange(proposal);
132 138
     double delLL = computeDeltaLL(change);
133
-    if (delLL * mAnnealingTemp > rejectProb)
139
+    if (delLL * mAnnealingTemp >= rejectProb)
134 140
     {
135 141
         domain.acceptProposal(proposal);
136 142
         change.label == 'A' ? mAMatrix.update(change) : mPMatrix.update(change);
... ...
@@ -323,6 +329,10 @@ void GibbsSampler::updateStatistics()
323 329
     for (unsigned r = 0; r < mPMatrix.nRow(); ++r)
324 330
     {
325 331
         normVec(r) = gaps::algo::sum(mPMatrix.getRow(r));
332
+        if (normVec(r) == 0)
333
+        {
334
+            normVec(r) = 1.0;
335
+        }
326 336
     }
327 337
 
328 338
     for (unsigned c = 0; c < mAMeanMatrix.nCol(); ++c)
... ...
@@ -11,4 +11,4 @@ OBJECTS =   Algorithms.o \
11 11
             cpp_tests/testMatrix.o \
12 12
             cpp_tests/testAlgorithms.o \
13 13
             cpp_tests/testAtomicSupport.o \
14
-            cpp_tests/testGibbsSampler.o
14
+            cpp_tests/testGibbsSampler.o
15 15
\ No newline at end of file
... ...
@@ -21,7 +21,8 @@ struct MatrixChange
21 21
     matrix_data_t delta2;
22 22
 
23 23
     MatrixChange(char l, unsigned r, unsigned c, double d)
24
-        : label(l), nChanges(1), row1(r), col1(c), delta1(d)
24
+        : label(l), nChanges(1), row1(r), col1(c), delta1(d), row2(0),
25
+        col2(0), delta2(0.0)
25 26
     {}
26 27
 
27 28
     MatrixChange(char l, unsigned r1, unsigned c1, double d1, unsigned r2,
... ...
@@ -4,12 +4,10 @@
4 4
 TEST_CASE("Test AtomicSupport.h")
5 5
 {
6 6
     unsigned nrow = 100, ncol = 500;
7
-    uint64_t maxAtoms = 100;
8 7
 
9 8
     SECTION("Make and Convert proposals")
10 9
     {
11 10
         AtomicSupport domain('A', nrow, ncol, 1.0, 0.5);
12
-        domain.setMaxNumAtoms(maxAtoms);
13 11
         REQUIRE(domain.alpha() == 1.0);
14 12
         REQUIRE(domain.lambda() == 0.5);
15 13
         
... ...
@@ -19,8 +17,6 @@ TEST_CASE("Test AtomicSupport.h")
19 17
         REQUIRE(prop.label == 'A');
20 18
         REQUIRE(prop.type == 'B');
21 19
         REQUIRE(prop.nChanges == 1);
22
-        REQUIRE(prop.pos1 == 0xf21db672ab2f52a4);
23
-        REQUIRE(prop.delta1 == 1.0);
24 20
         REQUIRE(prop.pos2 == 0);
25 21
         REQUIRE(prop.delta2 == 0.0);
26 22
         
... ...
@@ -55,10 +51,51 @@ TEST_CASE("Test AtomicSupport.h")
55 51
 
56 52
             double oldMass = domain.totalMass();
57 53
 
54
+            uint64_t nOld = domain.numAtoms();
55
+
58 56
             domain.acceptProposal(prop);
57
+
58
+            if (prop.type == 'B')
59
+            {
60
+                REQUIRE(domain.numAtoms() == nOld + 1);
61
+            }
62
+            else if (prop.type == 'D')
63
+            {
64
+                REQUIRE(domain.numAtoms() == nOld - 1);
65
+            }
66
+            /*else
67
+            {
68
+                REQUIRE(domain.numAtoms() == nOld);
69
+            }*/
59 70
         
60 71
             REQUIRE(domain.totalMass() == oldMass + prop.delta1 + prop.delta2);
61
-            REQUIRE(domain.numAtoms() <= maxAtoms);
62 72
         }
63 73
     }
74
+
75
+    SECTION("Proposal Distribution")
76
+    {
77
+        AtomicSupport domain('A', nrow, ncol, 0.01, 0.05);
78
+        
79
+        gaps::random::setSeed(1);
80
+
81
+        unsigned nB = 0, nD = 0, nM = 0, nE = 0;
82
+        for (unsigned i = 0; i < 5000; ++i)
83
+        {
84
+            AtomicProposal prop = domain.makeProposal();
85
+            domain.acceptProposal(prop);
86
+
87
+            switch(prop.type)
88
+            {
89
+                case 'B': nB++; break;
90
+                case 'D': nD++; break;
91
+                case 'M': nM++; break;
92
+                case 'E': nE++; break;
93
+            }
94
+        }
95
+        REQUIRE(domain.numAtoms() > 100);
96
+        REQUIRE(nB > 500);
97
+        //REQUIRE(nD > 500);
98
+        REQUIRE(nM > 500);
99
+        REQUIRE(nE > 500);
100
+    }
64 101
 }
65 102
\ No newline at end of file
... ...
@@ -22,7 +22,7 @@ TEST_CASE("Test GibbsSampler.h")
22 22
 
23 23
     SECTION("Create GibbsSampler")
24 24
     {
25
-        GibbsSampler sampler(rD, rS, 10, 1, 1, 5.0, 5.0, false);
25
+        GibbsSampler sampler(rD, rS, 10, 0.01, 0.01, 1.0, 1.0, false);
26 26
 
27 27
         REQUIRE(sampler.chi2() == 0.0);
28 28
         REQUIRE(sampler.totalNumAtoms('A') == 0);
... ...
@@ -31,12 +31,63 @@ TEST_CASE("Test GibbsSampler.h")
31 31
 
32 32
     SECTION("Update GibbsSampler")
33 33
     {
34
-        GibbsSampler sampler(rD, rS, 10, 1, 1, 5.0, 5.0, false);
34
+        GibbsSampler sampler(rD, rS, 10, 0.01, 0.01, 1.0, 1.0, false);
35 35
 
36
-        for (unsigned i = 0; i < 1000; ++i)
36
+        for (unsigned i = 0; i < 5000; ++i)
37 37
         {
38 38
             REQUIRE_NOTHROW(sampler.update('A'));
39 39
             REQUIRE_NOTHROW(sampler.update('P'));
40 40
         }
41
+        REQUIRE(sampler.totalNumAtoms('A') > 10);
42
+        REQUIRE(sampler.totalNumAtoms('P') > 10);
43
+    }
44
+
45
+    SECTION("GibbsSampler Statistics")
46
+    {
47
+        GibbsSampler sampler(rD, rS, 10, 0.01, 0.01, 1.0, 1.0, false);
48
+
49
+        for (unsigned i = 0; i < 1000; ++i)
50
+        {
51
+            for (unsigned j = 0; j < 10; ++j)
52
+            {
53
+                sampler.update('A');
54
+                sampler.update('P');
55
+            }
56
+            sampler.updateStatistics();
57
+        }
58
+        Rcpp::NumericMatrix AMean = sampler.AMeanRMatrix();
59
+        Rcpp::NumericMatrix AStd = sampler.AStdRMatrix();
60
+        Rcpp::NumericMatrix PMean = sampler.PMeanRMatrix();
61
+        Rcpp::NumericMatrix PStd = sampler.PStdRMatrix();
62
+
63
+        REQUIRE(AMean.nrow() == rD.nrow());
64
+        REQUIRE(AMean.ncol() == 10);
65
+
66
+        REQUIRE(AStd.nrow() == rD.nrow());
67
+        REQUIRE(AStd.ncol() == 10);
68
+
69
+        REQUIRE(PMean.nrow() == 10);
70
+        REQUIRE(PMean.ncol() == rD.ncol());
71
+
72
+        REQUIRE(PStd.nrow() == 10);
73
+        REQUIRE(PStd.ncol() == rD.ncol());
74
+
75
+        for (unsigned r = 0; r < AMean.nrow(); ++r)
76
+        {
77
+            for (unsigned c = 0; c < AMean.ncol(); ++c)
78
+            {
79
+                REQUIRE(AMean(r,c) >= 0.0);
80
+                REQUIRE(AStd(r,c) >= 0.0);
81
+            }
82
+        }
83
+
84
+        for (unsigned r = 0; r < PMean.nrow(); ++r)
85
+        {
86
+            for (unsigned c = 0; c < PMean.ncol(); ++c)
87
+            {
88
+                REQUIRE(PMean(r,c) >= 0.0);
89
+                REQUIRE(PStd(r,c) >= 0.0);
90
+            }
91
+        }
41 92
     }
42 93
 }
43 94
\ No newline at end of file
... ...
@@ -13,13 +13,19 @@ TEST_CASE("Test Random.h - Random Number Generation")
13 13
     SECTION("Test uniform distribution over unit interval")
14 14
     {
15 15
         double min = 1, max = 0;
16
-        for (unsigned i = 0; i < 10000; ++i)
16
+        double sum = 0.0;
17
+        unsigned N = 10000;
18
+        for (unsigned i = 0; i < N; ++i)
17 19
         {
18 20
             min = std::min(gaps::random::uniform(), min);
19 21
             max = std::max(gaps::random::uniform(), max);
22
+            sum += gaps::random::uniform();
20 23
         }
24
+        REQUIRE(sum / N == Approx(0.5).epsilon(0.01));
21 25
         REQUIRE(min >= 0);
26
+        REQUIRE(min < 0.01);
22 27
         REQUIRE(max <= 1);
28
+        REQUIRE(max > 0.99);
23 29
     }
24 30
 
25 31
     SECTION("Test uniform distribution over general interval")
... ...
@@ -1,6 +1,6 @@
1
-context("lints")
1
+#context("lints")
2 2
 
3
-test_that("Package Style", {
4
-    skip("Don't lint for now")
5
-    lintr::expect_lint_free()
6
-})
3
+#test_that("Package Style", {
4
+#    skip("Don't lint for now")
5
+#    lintr::expect_lint_free()
6
+#})
... ...
@@ -4,7 +4,9 @@ test_that("GAPS Simple Simulation",
4 4
 {
5 5
     data(SimpSim)
6 6
     nIter <- 1000
7
-    expect_error(gapsRun(SimpSim.D, SimpSim.S, nFactor=3, messages=FALSE), NA)
7
+    res <- gapsRun(SimpSim.D, SimpSim.S, nFactor=3, messages=FALSE)
8
+
9
+    expect_true(!is.na(res$meanChi2))
8 10
 })
9 11
 
10 12
 #test_that("GAPSmap Simple Simulation",