Browse code

tests passing - exchange needs to be fixed

sherman5 authored on 11/12/2017 23:32:14
Showing31 changed files

... ...
@@ -34,7 +34,7 @@ Imports:
34 34
 Suggests:
35 35
     testthat,
36 36
     lintr
37
-LinkingTo: Rcpp, BH, RcppArmadillo
37
+LinkingTo: Rcpp, BH
38 38
 License: GPL (==2)
39 39
 biocViews: GeneExpression, Transcription, GeneSetEnrichment,
40 40
     DifferentialExpression, Bayesian, Clustering, TimeCourse, RNASeq, Microarray,
... ...
@@ -1,8 +1,8 @@
1 1
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
2 2
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
3 3
 
4
-cogaps <- function(DMatrix, SMatrix, ABins, PBins, Config, ConfigNums, seed = -1L, messages = FALSE, singleCellRNASeq = FALSE) {
5
-    .Call('_CoGAPS_cogaps', PACKAGE = 'CoGAPS', DMatrix, SMatrix, ABins, PBins, Config, ConfigNums, seed, messages, singleCellRNASeq)
4
+cogaps <- function(DMatrix, SMatrix, nFactor, alphaA, alphaP, nEquil, nEquilCool, nSample, maxGibbsMassA, maxGibbsMassP, seed = -1L, messages = FALSE, singleCellRNASeq = FALSE) {
5
+    .Call('_CoGAPS_cogaps', PACKAGE = 'CoGAPS', DMatrix, SMatrix, nFactor, alphaA, alphaP, nEquil, nEquilCool, nSample, maxGibbsMassA, maxGibbsMassP, seed, messages, singleCellRNASeq)
6 6
 }
7 7
 
8 8
 run_catch_unit_tests <- function() {
... ...
@@ -64,123 +64,10 @@ gapsRun <- function(D, S, ABins = data.frame(), PBins = data.frame(),
64 64
                     alphaP = 0.01, nMaxP = 100000, max_gibbmass_paraP = 100.0,
65 65
                     seed=-1, messages=TRUE, singleCellRNASeq=FALSE)
66 66
 {
67
-  #Begin data type error checking code
68
-  charDataErrors = c(!is.character(simulation_id), !is.character(fixedDomain))
69
-  charCheck = c("simulation_id", "fixedDomain")
70
-
71
-  boolDataErrors = c(!is.logical(output_atomic), !is.logical(fixedBinProbs), !is.logical(sampleSnapshots))
72
-  boolCheck = c("output_atomic", "fixedBinProbs", "sampleSnapshots")
73
-
74
-  numericDataErrors = c(!is.numeric(nFactor), !is.numeric(nEquil), !is.numeric(nSample), !is.numeric(nOutR), !is.numeric(numSnapshots),
75
-                        !is.numeric(alphaA), !is.numeric(nMaxA), !is.numeric(max_gibbmass_paraA), !is.numeric(alphaP),
76
-                        !is.numeric(nMaxP), !is.numeric(max_gibbmass_paraP))
77
-  numericCheck = c("nFactor", "nEquil", "nSample", "nOutR", "numSnapshots", "alphaA", "nMaxA",
78
-                   "max_gibbmass_paraA", "alphaP",    "nMaxP", "max_gibbmass_paraP")
79
-
80
-  dataFrameErrors = c(!is.data.frame(D), !is.data.frame(S), !is.data.frame(ABins), !is.data.frame(PBins))
81
-  dataFrameCheck = c("D", "S", "ABins", "PBins")
82
-
83
-  matrixErrors = c(!is.matrix(D), !is.matrix(S), !is.matrix(ABins), !is.matrix(PBins))
84
-
85
-  if(any(charDataErrors))
86
-  {
87
-    #Check which of these is not a string
88
-    for(i in 1:length(charDataErrors))
89
-    {
90
-      if(charDataErrors[i] == TRUE)
91
-      {
92
-        stop(paste("Error in gapsRun: Argument",charCheck[i],"is of the incorrect type. Please see documentation for details."))
93
-      }
94
-    }
95
-  }
96
-
97
-  if(any(boolDataErrors))
98
-  {
99
-    #Check which of these is not a boolean
100
-    for(i in 1:length(boolDataErrors))
101
-    {
102
-      if(boolDataErrors[i] == TRUE)
103
-      {
104
-        stop(paste("Error in gapsRun: Argument",boolCheck[i],"is of the incorrect type. Please see documentation for details."))
105
-      }
106
-    }
107
-  }
108
-
109
-  if(any(numericDataErrors))
110
-  {
111
-    #Check which of these is not an integer/double
112
-    for(i in 1:length(numericDataErrors))
113
-    {
114
-      if(numericDataErrors[i] == TRUE)
115
-      {
116
-        stop(paste("Error in gapsRun: Argument",numericCheck[i],"is of the incorrect type. Please see documentation for details."))
117
-
118
-      }
119
-    }
120
-  }
121
-
122
-  #At least one of A, P, ABins, or PBins is not a matrix or data.frame
123
-  if(any((dataFrameErrors[1] && matrixErrors[1]), (dataFrameErrors[2] && matrixErrors[2]), (dataFrameErrors[3] && matrixErrors[3]), (dataFrameErrors[4] && matrixErrors[4])))
124
-  {
125
-    for(i in 1:length(dataFrameCheck))
126
-    {
127
-      if((dataFrameErrors[i] && matrixErrors[i]) == TRUE)
128
-      {
129
-        stop(paste("Error in gapsRun: Argument",dataFrameCheck[i],"is not a matrix or data.frame. Please see documentation for details."))
130
-
131
-      }
132
-    }
133
-  }
134
-
135 67
   #Floor the parameters that are integers to prevent allowing doubles.
136 68
   nFactor = floor(nFactor)
137 69
   nEquil = floor(nEquil)
138 70
   nSample = floor(nSample)
139
-  nOutR = floor(nOutR)
140
-  numSnapshots = floor(numSnapshots)
141
-  nMaxA = floor(nMaxA)
142
-  nMaxP = floor(nMaxP)
143
-
144
-
145
-  # pass all settings to C++ within a list
146
-  #    if (is.null(P)) {
147
-  Config = c(simulation_id, output_atomic, fixedBinProbs, fixedDomain, sampleSnapshots);
148
-
149
-  ConfigNums = c(nFactor, nEquil, nSample, nOutR, alphaA, nMaxA, max_gibbmass_paraA,
150
-                 alphaP, nMaxP, max_gibbmass_paraP, numSnapshots);
151
-
152
-  #Begin logic error checking code
153
-
154
-  #Check for negative or zero arguments
155
-  if(any(ConfigNums <= 0))
156
-  {
157
-    stop("Error in gapsRun: Numeric Arguments must be strictly positive.")
158
-  }
159
-
160
-  #Check for nonsensical inputs (such as numSnapshots < nEquil or nSample)
161
-  if((numSnapshots > nEquil) || (numSnapshots > nSample))
162
-  {
163
-    stop("Error in gapsRun: Cannot have more snapshots of A and P than equilibration and/or sampling iterations.")
164
-  }
165
-
166
-  if((nOutR > nEquil) || (nOutR > nSample))
167
-  {
168
-    stop("Error in gapsRun: Cannot have more output steps than equilibration and/or sampling iterations.")
169
-  }
170
-
171
-  if(nFactor > (ncol(D)))
172
-  {
173
-    warning("Warning in gapsRun: Number of requested patterns greater than columns of Data Matrix.")
174
-  }
175
-
176
-  #        P <- as.data.frame(matrix(nrow=1,c(1,1,1))) # make something to pass
177
-  #    } else {
178
-  #        Config = c(nFactor, simulation_id, nEquil, nSample, nOutR,
179
-  #        output_atomic, alphaA, nMaxA, max_gibbmass_paraA, lambdaA_scale_factor,
180
-  #        alphaP, nMaxP, max_gibbmass_paraP, lambdaP_scale_factor, 1)
181
-
182
-  #   }
183
-
184 71
 
185 72
   geneNames = rownames(D);
186 73
   sampleNames = colnames(D);
... ...
@@ -193,45 +80,45 @@ gapsRun <- function(D, S, ABins = data.frame(), PBins = data.frame(),
193 80
   }
194 81
 
195 82
   # call to C++ Rcpp code
196
-  cogapResult = cogaps(D, S, ABins, PBins, Config, ConfigNums, seed, messages,
197
-    singleCellRNASeq);
83
+  cogapResult = cogaps(as.matrix(D), as.matrix(S), nFactor, alphaA, alphaP, nEquil, 10, nSample, max_gibbmass_paraA,
84
+    max_gibbmass_paraP, seed, messages, singleCellRNASeq)
198 85
 
199 86
   # convert returned files to matrices to simplify visualization and processing
200
-  cogapResult$Amean = as.matrix(cogapResult$Amean);
201
-  cogapResult$Asd = as.matrix(cogapResult$Asd);
202
-  cogapResult$Pmean = as.matrix(cogapResult$Pmean);
203
-  cogapResult$Psd = as.matrix(cogapResult$Psd);
204
-
205
-  # label matrices
206
-  colnames(cogapResult$Amean) = patternNames;
207
-  rownames(cogapResult$Amean) = geneNames;
208
-  colnames(cogapResult$Asd) = patternNames;
209
-  rownames(cogapResult$Asd) = geneNames;
210
-  colnames(cogapResult$Pmean) = sampleNames;
211
-  rownames(cogapResult$Pmean) = patternNames;
212
-  colnames(cogapResult$Psd) = sampleNames;
213
-  rownames(cogapResult$Psd) = patternNames;
214
-
215
-  # calculate chi-squared of mean, this should be smaller than individual
216
-  # chi-squared sample values if sampling is good
217
-  calcChiSq = c(0);
218
-  MMatrix = (cogapResult$Amean %*% cogapResult$Pmean);
219
-
220
-
221
-  for(i in 1:(nrow(MMatrix)))
222
-  {
223
-    for(j in 1:(ncol(MMatrix)))
224
-    {
225
-      calcChiSq = (calcChiSq) + ((D[i,j] - MMatrix[i,j])/(S[i,j]))^(2);
226
-    }
227
-  }
228
-
229
-  cogapResult = c(cogapResult, calcChiSq);
230
-
231
-  names(cogapResult)[13] = "meanChi2";
232
-
233
-  if (messages)
234
-    message(paste("Chi-Squared of Mean:",calcChiSq))
87
+  #cogapResult$Amean = as.matrix(cogapResult$Amean);
88
+  #cogapResult$Asd = as.matrix(cogapResult$Asd);
89
+  #cogapResult$Pmean = as.matrix(cogapResult$Pmean);
90
+  #cogapResult$Psd = as.matrix(cogapResult$Psd);
91
+#
92
+  ## label matrices
93
+  #colnames(cogapResult$Amean) = patternNames;
94
+  #rownames(cogapResult$Amean) = geneNames;
95
+  #colnames(cogapResult$Asd) = patternNames;
96
+  #rownames(cogapResult$Asd) = geneNames;
97
+  #colnames(cogapResult$Pmean) = sampleNames;
98
+  #rownames(cogapResult$Pmean) = patternNames;
99
+  #colnames(cogapResult$Psd) = sampleNames;
100
+  #rownames(cogapResult$Psd) = patternNames;
101
+#
102
+  ## calculate chi-squared of mean, this should be smaller than individual
103
+  ## chi-squared sample values if sampling is good
104
+  #calcChiSq = c(0);
105
+  #MMatrix = (cogapResult$Amean %*% cogapResult$Pmean);
106
+#
107
+#
108
+  #for(i in 1:(nrow(MMatrix)))
109
+  #{
110
+  #  for(j in 1:(ncol(MMatrix)))
111
+  #  {
112
+  #    calcChiSq = (calcChiSq) + ((D[i,j] - MMatrix[i,j])/(S[i,j]))^(2);
113
+  #  }
114
+  #}
115
+#
116
+  #cogapResult = c(cogapResult, calcChiSq);
117
+#
118
+  #names(cogapResult)[13] = "meanChi2";
119
+#
120
+  #if (messages)
121
+  #  message(paste("Chi-Squared of Mean:",calcChiSq))
235 122
 
236 123
   return(cogapResult);
237 124
 }
238 125
new file mode 100644
... ...
@@ -0,0 +1,164 @@
1
+#include "Algorithms.h"
2
+#include "Matrix.h"
3
+
4
+#include <algorithm>
5
+
6
+#define GAPS_SQ(x) ((x) * (x))
7
+
8
+matrix_data_t gaps::algo::mean(const TwoWayMatrix &mat)
9
+{
10
+    return gaps::algo::sum(mat) / (mat.nRow() * mat.nCol());
11
+}
12
+
13
+matrix_data_t gaps::algo::sum(const Vector &vec)
14
+{
15
+    matrix_data_t sum = 0.0;
16
+    for (unsigned i = 0; i < vec.size(); ++i)
17
+    {
18
+        sum += vec(i);
19
+    }
20
+    return sum;
21
+}
22
+
23
+matrix_data_t gaps::algo::sum(const TwoWayMatrix &mat)
24
+{
25
+    matrix_data_t sum = 0.0;
26
+    for (unsigned r = 0; r < mat.nRow(); ++r)
27
+    {
28
+        sum += gaps::algo::sum(mat.getRow(r));
29
+    }
30
+    return sum;
31
+}
32
+
33
+matrix_data_t gaps::algo::nonZeroMean(const TwoWayMatrix &mat)
34
+{
35
+    matrix_data_t sum = 0.0;
36
+    unsigned int nNonZeros = 0;
37
+    for (unsigned int r = 0; r < mat.nRow(); ++r)
38
+    {
39
+        for (unsigned int c = 0; c < mat.nCol(); ++c)
40
+        {
41
+            if (mat(r,c) != 0.0)
42
+            {
43
+                nNonZeros++;
44
+                sum += mat(r,c);
45
+            }
46
+        }
47
+    }
48
+    return sum / nNonZeros;
49
+}
50
+
51
+bool gaps::algo::isRowZero(const RowMatrix &mat, unsigned row)
52
+{
53
+    return gaps::algo::sum(mat.getRow(row)) == 0;
54
+}
55
+    
56
+bool gaps::algo::isColZero(const ColMatrix &mat, unsigned col)
57
+{
58
+    return gaps::algo::sum(mat.getCol(col)) == 0;
59
+}
60
+
61
+static double deltaLL_A(unsigned row, unsigned col, double delta,
62
+const TwoWayMatrix &D, const TwoWayMatrix &S, const ColMatrix &A,
63
+const RowMatrix &P, const TwoWayMatrix &AP)
64
+{
65
+    double numer = 0.0, denom = 0.0, delLL = 0.0;
66
+    for (unsigned j = 0; j < D.nCol(); ++j)
67
+    {
68
+        numer = 2 * delta * (D(row,j) - AP(row,j)) * P(col,j) - GAPS_SQ(delta * P(col,j));
69
+        denom = 2 * GAPS_SQ(S(row,j));
70
+        delLL += numer / denom;
71
+    }
72
+    return delLL;
73
+}
74
+
75
+static double deltaLL_P(unsigned row, unsigned col, double delta,
76
+const TwoWayMatrix &D, const TwoWayMatrix &S, const ColMatrix &A,
77
+const RowMatrix &P, const TwoWayMatrix &AP)
78
+{
79
+    double numer = 0.0, denom = 0.0, delLL = 0.0;
80
+    for (unsigned i = 0; i < D.nRow(); ++i)
81
+    {
82
+        numer = 2 * delta * (D(i,col) - AP(i,col)) * A(i,row) - GAPS_SQ(delta * A(i,row));
83
+        denom = 2 * GAPS_SQ(S(i,col));
84
+        delLL += numer / denom;
85
+    }
86
+    return delLL;
87
+}
88
+
89
+double gaps::algo::deltaLL(const MatrixChange &ch, const TwoWayMatrix &D,
90
+const TwoWayMatrix &S, const ColMatrix &A, const RowMatrix &P,
91
+const TwoWayMatrix &AP)
92
+{
93
+    if (ch.nChanges == 1 && ch.label == 'A')
94
+    {
95
+        return deltaLL_A(ch.row1, ch.col1, ch.delta1, D, S, A, P, AP);
96
+    }
97
+    else if (ch.nChanges == 1 && ch.label == 'P')
98
+    {
99
+        return deltaLL_P(ch.row1, ch.col1, ch.delta1, D, S, A, P, AP);
100
+    }
101
+    else if (ch.nChanges == 2 && ch.label == 'A')
102
+    {
103
+        return deltaLL_A(ch.row1, ch.col1, ch.delta1, D, S, A, P, AP)
104
+            + deltaLL_A(ch.row2, ch.col2, ch.delta2, D, S, A, P, AP);
105
+    }
106
+    else
107
+    {
108
+        return deltaLL_P(ch.row1, ch.col1, ch.delta1, D, S, A, P, AP)
109
+            + deltaLL_P(ch.row2, ch.col2, ch.delta2, D, S, A, P, AP);
110
+    }
111
+}
112
+
113
+static AlphaParameters alphaParameters_A(unsigned row, unsigned col,
114
+const TwoWayMatrix &D, const TwoWayMatrix &S, const RowMatrix &P,
115
+const TwoWayMatrix &AP)
116
+{
117
+    double s = 0.0, su = 0.0, ratio = 0.0;
118
+    for (unsigned j = 0; j < D.nCol(); ++j)
119
+    {
120
+        ratio = P(col,j) / S(row,j);
121
+        s += GAPS_SQ(ratio);
122
+        su += P(col,j) * (D(row,j) - AP(row,j)) / GAPS_SQ(S(row,j));
123
+    }
124
+    return AlphaParameters(s, su);
125
+}
126
+
127
+static AlphaParameters alphaParameters_P(unsigned row, unsigned col,
128
+const TwoWayMatrix &D, const TwoWayMatrix &S, const ColMatrix &A,
129
+const TwoWayMatrix &AP)
130
+{
131
+    double s = 0.0, su = 0.0, ratio = 0.0;
132
+    for (unsigned i = 0; i < D.nRow(); ++i)
133
+    {
134
+        ratio = A(i,row) / S(i,col);
135
+        s += GAPS_SQ(ratio);
136
+        su += A(i,row) * (D(i,col) - AP(i,col)) / GAPS_SQ(S(i,col));
137
+    }
138
+    return AlphaParameters(s, su);
139
+}
140
+
141
+AlphaParameters gaps::algo::alphaParameters(const MatrixChange &ch,
142
+const TwoWayMatrix &D, const TwoWayMatrix &S, const ColMatrix &A,
143
+const RowMatrix &P, const TwoWayMatrix &AP)
144
+{
145
+    if (ch.nChanges == 1 && ch.label == 'A')
146
+    {
147
+        return alphaParameters_A(ch.row1, ch.col1, D, S, P, AP);
148
+    }
149
+    else if (ch.nChanges == 1 && ch.label == 'P')
150
+    {
151
+        return alphaParameters_P(ch.row1, ch.col1, D, S, A, AP);
152
+    }
153
+    else if (ch.nChanges == 2 && ch.label == 'A')
154
+    {
155
+        return alphaParameters_A(ch.row1, ch.col1, D, S, P, AP)
156
+            + alphaParameters_A(ch.row2, ch.col2, D, S, P, AP);
157
+    }
158
+    else
159
+    {
160
+        return alphaParameters_P(ch.row1, ch.col1, D, S, A, AP)
161
+            + alphaParameters_P(ch.row2, ch.col2, D, S, A, AP);
162
+    }
163
+}
164
+
0 165
new file mode 100644
... ...
@@ -0,0 +1,46 @@
1
+#ifndef __COGAPS_ALGORITHMS_H__
2
+#define __COGAPS_ALGORITHMS_H__
3
+
4
+#include "Matrix.h"
5
+
6
+struct AlphaParameters
7
+{
8
+    double s;
9
+    double su;
10
+    
11
+    AlphaParameters(double inS, double inSU)
12
+        : s(inS), su(inSU)
13
+    {}
14
+
15
+    AlphaParameters& operator+(const AlphaParameters &other)
16
+    {
17
+        s += other.s;
18
+        su -= other.su; // weird
19
+        return *this;
20
+    }
21
+};
22
+
23
+namespace gaps
24
+{
25
+
26
+namespace algo
27
+{
28
+    matrix_data_t sum(const Vector &vec);
29
+    matrix_data_t sum(const TwoWayMatrix &mat);
30
+    matrix_data_t mean(const TwoWayMatrix &mat);
31
+    matrix_data_t nonZeroMean(const TwoWayMatrix &mat);
32
+    bool isRowZero(const RowMatrix &mat, unsigned row);
33
+    bool isColZero(const ColMatrix &mat, unsigned col);
34
+
35
+    double deltaLL(const MatrixChange &ch, const TwoWayMatrix &D,
36
+        const TwoWayMatrix &S, const ColMatrix &A,
37
+        const RowMatrix &P, const TwoWayMatrix &AP);
38
+
39
+    AlphaParameters alphaParameters(const MatrixChange &ch,
40
+        const TwoWayMatrix &D, const TwoWayMatrix &S, const ColMatrix &A,
41
+        const RowMatrix &P, const TwoWayMatrix &AP);
42
+}
43
+
44
+}
45
+
46
+#endif
0 47
\ No newline at end of file
... ...
@@ -14,11 +14,13 @@
14 14
 
15 15
 #define EPSILON 1E-10
16 16
 
17
-AtomicSupport::AtomicSupport(uint64_t nrow, uint64_t ncol)
18
-    : mNumRows(nrow), mNumCols(ncol), mNumBins(nrow * ncol),
19
-      mNumAtoms(0), mTotalMass(0.0),
20
-      mMaxNumAtoms(std::numeric_limits<uint64_t>::max()),
21
-      mBinSize(std::numeric_limits<uint64_t>::max() / (nrow * ncol))
17
+AtomicSupport::AtomicSupport(char label, uint64_t nrow, uint64_t ncol,
18
+double alpha, double lambda)
19
+    :
20
+mNumRows(nrow), mNumCols(ncol), mNumBins(nrow * ncol),
21
+mNumAtoms(0), mTotalMass(0.0), mLabel(label), mAlpha(alpha), 
22
+mMaxNumAtoms(std::numeric_limits<uint64_t>::max()), mLambda(lambda),
23
+mBinSize(std::numeric_limits<uint64_t>::max() / (nrow * ncol))
22 24
 {}
23 25
 
24 26
 uint64_t AtomicSupport::getRow(uint64_t pos) const
... ...
@@ -55,14 +57,14 @@ AtomicProposal AtomicSupport::proposeDeath() const
55 57
 {
56 58
     uint64_t location = randomAtomPosition();
57 59
     uint64_t mass = mAtomicDomain.at(location);
58
-    return AtomicProposal('D', location, -mass);
60
+    return AtomicProposal(mLabel, 'D', location, -mass);
59 61
 }
60 62
 
61 63
 AtomicProposal AtomicSupport::proposeBirth() const
62 64
 {
63 65
     uint64_t location = randomFreePosition();
64 66
     uint64_t mass = gaps::random::exponential(mLambda);
65
-    return AtomicProposal('B', location, mass);
67
+    return AtomicProposal(mLabel, 'B', location, mass);
66 68
 }
67 69
 
68 70
 // move atom between adjacent atoms
... ...
@@ -82,7 +84,7 @@ AtomicProposal AtomicSupport::proposeMove() const
82 84
 
83 85
     uint64_t newLocation = gaps::random::uniform64(left, right);
84 86
 
85
-    return AtomicProposal('M', location, -mass, newLocation, mass);
87
+    return AtomicProposal(mLabel, 'M', location, -mass, newLocation, mass);
86 88
 }
87 89
 
88 90
 // exchange with adjacent atom to the right
... ...
@@ -111,8 +113,8 @@ AtomicProposal AtomicSupport::proposeExchange() const
111 113
 
112 114
     // preserve total mass
113 115
     return mass1 > mass2 ?
114
-        AtomicProposal('E', pos1, newMass - mass1, pos2, mass1 - newMass)
115
-        : AtomicProposal('E', pos1, mass2 - newMass, pos2, newMass - mass2);
116
+        AtomicProposal(mLabel, 'E', pos1, newMass - mass1, pos2, mass1 - newMass)
117
+        : AtomicProposal(mLabel, 'E', pos1, mass2 - newMass, pos2, newMass - mass2);
116 118
 }
117 119
 
118 120
 void AtomicSupport::updateAtomMass(uint64_t pos, double delta)
... ...
@@ -180,13 +182,12 @@ MatrixChange AtomicSupport::getMatrixChange(const AtomicProposal &prop) const
180 182
 {
181 183
     if (prop.nChanges > 1)
182 184
     {
183
-        return MatrixChange(prop.type, getRow(prop.pos1), getCol(prop.pos1),
184
-            prop.delta1, getRow(prop.pos2), getCol(prop.pos2), prop.delta2);
185
+        return MatrixChange(prop.label, getRow(prop.pos1), getCol(prop.pos1), prop.delta1,
186
+            getRow(prop.pos2), getCol(prop.pos2), prop.delta2);
185 187
     }
186 188
     else
187 189
     {   
188
-        return MatrixChange(prop.type, getRow(prop.pos1), getCol(prop.pos1),
189
-            prop.delta1);
190
+        return MatrixChange(prop.label, getRow(prop.pos1), getCol(prop.pos1), prop.delta1);
190 191
     }
191 192
 }
192 193
 
... ...
@@ -11,6 +11,7 @@
11 11
 
12 12
 struct AtomicProposal
13 13
 {
14
+    char label;
14 15
     char type;
15 16
     unsigned nChanges;
16 17
 
... ...
@@ -20,12 +21,12 @@ struct AtomicProposal
20 21
     uint64_t pos2;
21 22
     double delta2;
22 23
 
23
-    AtomicProposal(char t, uint64_t p1, double m1)
24
-        : type(t), nChanges(1), pos1(p1), delta1(m1), pos2(0), delta2(0.0)
24
+    AtomicProposal(char l, char t, uint64_t p1, double m1)
25
+        : label(l), type(t), nChanges(1), pos1(p1), delta1(m1), pos2(0), delta2(0.0)
25 26
     {}
26 27
 
27
-    AtomicProposal(char t, uint64_t p1, double m1, uint64_t p2, double m2)
28
-        : type(t), nChanges(2), pos1(p1), delta1(m1), pos2(p2), delta2(m2)
28
+    AtomicProposal(char l, char t, uint64_t p1, double m1, uint64_t p2, double m2)
29
+        : label(l), type(t), nChanges(2), pos1(p1), delta1(m1), pos2(p2), delta2(m2)
29 30
     {}
30 31
 };
31 32
 
... ...
@@ -33,6 +34,9 @@ class AtomicSupport
33 34
 {
34 35
 private:
35 36
 
37
+    // label of this atomic domain (A/P)
38
+    char mLabel;
39
+
36 40
     // storage of the atomic domain
37 41
     std::map<uint64_t, double> mAtomicDomain;
38 42
     uint64_t mNumAtoms;
... ...
@@ -47,10 +51,10 @@ private:
47 51
     uint64_t mNumBins;
48 52
     uint64_t mBinSize;
49 53
 
50
-    // average number of atoms per bin
54
+    // average number of atoms per bin (must be > 0)
51 55
     double mAlpha;
52 56
 
53
-    // expected magnitude of each atom
57
+    // expected magnitude of each atom (must be > 0)
54 58
     double mLambda;     
55 59
 
56 60
     // convert atomic position to row/col of the matrix
... ...
@@ -73,7 +77,8 @@ private:
73 77
 public:
74 78
 
75 79
     // constructor
76
-    AtomicSupport(uint64_t nrow, uint64_t ncol);
80
+    AtomicSupport(char label, uint64_t nrow, uint64_t ncol, double alpha=1.0,
81
+        double lambda=1.0);
77 82
 
78 83
     // create and accept a proposal
79 84
     AtomicProposal makeProposal() const;
... ...
@@ -95,9 +100,6 @@ public:
95 100
     void setAlpha(double alpha) {mAlpha = alpha;}
96 101
     void setLambda(double lambda) {mLambda = lambda;}
97 102
     void setMaxNumAtoms(uint64_t max) {mMaxNumAtoms = max;}
98
-
99
-    // TODO remove support for this
100
-    std::map<uint64_t, double> getDomain() {return mAtomicDomain;}
101 103
 };
102 104
 
103 105
 #endif
104 106
new file mode 100644
... ...
@@ -0,0 +1,66 @@
1
+#include "GibbsSampler.h"
2
+
3
+#include <Rcpp.h>
4
+#include <ctime>
5
+
6
+void runGibbsSampler(GibbsSampler &sampler, unsigned iterations,
7
+    bool updateTemp=false);
8
+
9
+// [[Rcpp::export]]
10
+Rcpp::List cogaps(Rcpp::NumericMatrix DMatrix, Rcpp::NumericMatrix SMatrix,
11
+unsigned nFactor, double alphaA, double alphaP, unsigned nEquil,
12
+unsigned nEquilCool, unsigned nSample, double maxGibbsMassA,
13
+double maxGibbsMassP, int seed=-1, bool messages=false, bool singleCellRNASeq=false)
14
+{
15
+    // set seed
16
+    uint32_t seedUsed = seed >= 0 ? static_cast<uint32_t>(seed)
17
+        : static_cast<uint32_t>(std::time(0));
18
+    gaps::random::setSeed(seedUsed);
19
+
20
+    // create the gibbs sampler
21
+    GibbsSampler sampler(DMatrix, SMatrix, nFactor, alphaA, alphaP,
22
+        maxGibbsMassA, maxGibbsMassP, singleCellRNASeq);
23
+
24
+    // run the sampler for each phase of the algorithm
25
+    runGibbsSampler(sampler, nEquil, true);
26
+    runGibbsSampler(sampler, nEquilCool);
27
+    runGibbsSampler(sampler, nSample);
28
+
29
+    //Just leave the snapshots as empty lists
30
+    return Rcpp::List::create(Rcpp::Named("randSeed") = seedUsed,
31
+                              Rcpp::Named("numAtomA") = sampler.totalNumAtoms('A'));
32
+}
33
+
34
+void runGibbsSampler(GibbsSampler &sampler, unsigned iterations, bool updateTemp)
35
+{
36
+    unsigned nIterA = 10;
37
+    unsigned nIterP = 10;
38
+        
39
+    double tempDenom = (double)iterations / 2.0;
40
+
41
+    for (unsigned i = 0; i < iterations; ++i)
42
+    {
43
+        if (updateTemp)
44
+        {
45
+            sampler.setAnnealingTemp(std::min((i + 1) / tempDenom, 1.0));
46
+        }
47
+
48
+        for (unsigned j = 0; j < nIterA; ++j)
49
+        {
50
+            sampler.update('A');
51
+        }
52
+
53
+        for (unsigned j = 0; j < nIterP; ++j)
54
+        {
55
+            sampler.update('P');
56
+        }
57
+
58
+        double tempChiSq = sampler.chi2();
59
+        uint64_t numAtomsA = sampler.totalNumAtoms('A');
60
+        uint64_t numAtomsP = sampler.totalNumAtoms('P');
61
+
62
+        nIterA = gaps::random::poisson(std::max(numAtomsA, (uint64_t)10));
63
+        nIterP = gaps::random::poisson(std::max(numAtomsP, (uint64_t)10));
64
+    }
65
+
66
+}
0 67
\ No newline at end of file
1 68
deleted file mode 100644
... ...
@@ -1,397 +0,0 @@
1
-// cogapsR.cpp
2
-
3
-// =============================================================================
4
-// This is the main code for Cogaps. (7th Sep, 2013)
5
-// This file also works as an interface to R via Rcpp (24th July, 2014)
6
-// =============================================================================
7
-
8
-#include <iostream>       // for use with standard I/O
9
-#include <string>         // for string processing
10
-#include <fstream>        // for output to files
11
-#include <limits>         // for extracting numerical limits of C++
12
-#include <ctime>
13
-
14
-
15
-// ------ incorporated to use Cogaps_options ------------
16
-#include <vector>
17
-#include <iomanip>
18
-#include <boost/algorithm/string.hpp>
19
-// ------------------------------------------------------
20
-#include "Random.h"   // for incorporating a random number generator.
21
-#include <RcppArmadillo.h>
22
-#include "Matrix.h"  // for incorporating a Matrix class
23
-#include "AtomicSupport.h"  // for incorporating an Atomic class
24
-#include "GAPSNorm.h"  // for incorporating calculation of statistics in cogaps.
25
-#include "GibbsSampler.h" // for incorporating the GibbsSampler which
26
-// does all the atomic space to matrix conversion
27
-// and sampling actions.
28
-
29
-// ------------------------------------------------------
30
-
31
-//namespace bpo = boost::program_options;
32
-using namespace std;
33
-using std::vector;
34
-
35
-// [[Rcpp::export]]
36
-Rcpp::List cogaps(Rcpp::NumericMatrix DMatrix, Rcpp::NumericMatrix SMatrix,
37
-Rcpp::NumericMatrix ABins, Rcpp::NumericMatrix PBins,
38
-Rcpp::CharacterVector Config, Rcpp::NumericVector ConfigNums, int seed=-1, 
39
-bool messages=false, bool singleCellRNASeq=false)
40
-{
41
-    // set seed
42
-    uint32_t seedUsed = seed >= 0 ? static_cast<uint32_t>(seed)
43
-        : static_cast<uint32_t>(std::time(0));
44
-    gaps::random::setSeed(seedUsed);
45
-
46
-    //---------------------
47
-    // ===========================================================================
48
-    // Part 1) Initialization:
49
-    // In this section, we read in the system parameters from the paremter file
50
-    // parameter.txt, and matrices D and S from datafile.txt.
51
-    // Then we initialize A and P in both their atomic domains and
52
-    // matrix forms.
53
-    // ===========================================================================
54
-//---------------------------------------------
55
-//CK CODE FOR CHANGING R VARIABLES OF THE CONFIG FILE INTO C VARIABLES
56
-//THIS REPLACES COGAPS_PROGAM_OPTIONS
57
-    string temp;
58
-    double tempNumInput;
59
-    tempNumInput = (ConfigNums[1]);
60
-    unsigned long nEquil = (tempNumInput);
61
-    //For equilibration
62
-    tempNumInput = (ConfigNums[2]);
63
-    unsigned long nSample = (tempNumInput);
64
-    // For Sampling
65
-    tempNumInput = (ConfigNums[3]);
66
-    unsigned long nObsR = tempNumInput;
67
-    //How many times we output the Chi-Sq
68
-    tempNumInput = (ConfigNums[0]);
69
-    unsigned int nFactor = tempNumInput;
70
-//Number of patterns
71
-    tempNumInput = (ConfigNums[4]);
72
-    double alphaA = tempNumInput;
73
-//scale parameter A
74
-    tempNumInput = (ConfigNums[7]);
75
-    double alphaP = tempNumInput;
76
-    //scale parameter P
77
-    tempNumInput = (ConfigNums[5]);
78
-    double nMaxA = tempNumInput;
79
-    //max atoms in A domain
80
-    tempNumInput = (ConfigNums[8]);
81
-    double nMaxP = tempNumInput;
82
-    //max atoms in P domain
83
-    string simulation_id = Rcpp::as<string>(Config[0]);
84
-    //Variable to name optional output files
85
-    tempNumInput = (ConfigNums[6]);
86
-    double max_gibbsmass_paraA = tempNumInput;
87
-//max atom mass in A domain
88
-    tempNumInput = (ConfigNums[9]);
89
-    double max_gibbsmass_paraP = tempNumInput;
90
-    //max atom mass in P domain
91
-    //Lambda Scale Factors removed and made to be constant. Can be hard coded from GibbsSampler.cpp
92
-    temp = Rcpp::as<string>(Config[1]);
93
-    bool Q_output_atomic;
94
-
95
-    if (temp == "TRUE" || temp == "true") {
96
-        Q_output_atomic = true;    // whether to output the atomic space info to a file specified by simulation ID
97
-
98
-    } else {
99
-        Q_output_atomic = false;
100
-    }
101
-
102
-    temp = Rcpp::as<string>(Config[2]);
103
-    bool fixBinProbs;
104
-
105
-    if (temp == "TRUE" || temp == "true") {
106
-        fixBinProbs = true;    // whether we are allowing variable bin sizes.
107
-
108
-    } else {
109
-        fixBinProbs = false;
110
-    }
111
-
112
-    temp = Rcpp::as<string>(Config[3]); //The Domain to allow variable bin sizes for (A, P, Both or Neither)
113
-    string fixedDomainStr = temp;
114
-    temp = Rcpp::as<string>(Config[4]); //Whether or not to run snapshots
115
-    bool SampleSnapshots;
116
-
117
-    if (temp == "TRUE" || temp == "true") {
118
-        SampleSnapshots = true;
119
-
120
-    } else {
121
-        SampleSnapshots = false;
122
-    }
123
-
124
-    tempNumInput = (ConfigNums[10]); //The number of snapshots saved during the Sampling (nSample/numSnapshots) increments
125
-    int numSnapshots = tempNumInput;
126
-    //Code to make the D and S matrices read from R into C++ vectors to make into Elana's Matrix Objects in Matrix.cpp
127
-    //Now also allows for variable bin sizes to be created
128
-    //Code to establish the sizes and initialize the C++ vectors to pass
129
-
130
-    //--------------------END CREATING D and S C++ VECTORS
131
-    // Parameters or structures to be calculated or constructed:
132
-    unsigned long nIterA = 10;    // initial inner loop iterations for A
133
-    unsigned long nIterP = 10;    // initial inner loop iterations for P
134
-    unsigned long atomicSize = 0; // number of atomic points
135
-
136
-    // Initialize the GibbsSampler.
137
-    GibbsSampler GibbsSamp(DMatrix, SMatrix, nFactor, alphaA, alphaP);
138
-
139
-    // ---------------------------------------------------------------------------
140
-    // Based on the information of D, construct and initialize for A and P both
141
-    // the matrices and atomic spaces.
142
-    //GibbsSamp.init_AMatrix_and_PMatrix(); // initialize A and P matrices
143
-    //This Section now is to handle the many possibilities for Variable Bin Sizes (Priors)
144
-    //A for variable A bins, P for variable P Bins, B for both and N for regular uniform bin sizes
145
-    char fixedDomain = fixedDomainStr[0];
146
-
147
-    // ===========================================================================
148
-    // Part 2) Equilibration:
149
-    // In this section, we let the system eqilibrate with nEquil outer loop
150
-    // iterations. Within each outer loop iteration, A is iterated nIterA times
151
-    // and P is iterated nIterP times. After equilibration, we update nIterA and
152
-    // nIterP according to the expected number of atoms in the atomic spaces
153
-    // of A and P respectively.
154
-    // ===========================================================================
155
-    double chi2;
156
-    double tempChiSq;
157
-    double tempAtomA;
158
-    double tempAtomP;
159
-    int outCount = 0;
160
-    int numOutputs = nObsR;                  //R'S OUTPUT CONTROL
161
-    int totalChiSize = nSample + nEquil;
162
-    Rcpp::NumericVector chiVect(totalChiSize); //INITIALIZE THE VECTOR HOLDING THE CHISQUARE.
163
-    Rcpp::NumericVector nAEquil(nEquil);       //INITIALIZE THE VECTOR HOLDING THE ATOMS FOR EACH MATRIX FOR EACH EQUIL/SAMP
164
-    Rcpp::NumericVector nASamp(nSample);
165
-    Rcpp::NumericVector nPEquil(nEquil);
166
-    Rcpp::NumericVector nPSamp(nSample);
167
-
168
-    for (unsigned long ext_iter = 1; ext_iter <= nEquil; ++ext_iter)
169
-    {
170
-        GibbsSamp.set_AnnealingTemperature();
171
-
172
-        for (unsigned long iterA = 1; iterA <= nIterA; ++iterA)
173
-        {
174
-            GibbsSamp.update('A');
175
-        }
176
-
177
-        for (unsigned long iterP = 1; iterP <= nIterP; ++iterP)
178
-        {
179
-            GibbsSamp.update('P');
180
-        }
181
-
182
-        //Finds the current ChiSq and places it into the vector to be returned to R (and output on occasion)
183
-        tempChiSq = GibbsSamp.get_sysChi2();
184
-        chiVect[(ext_iter) - 1] = tempChiSq;
185
-        // ----------- output computing info ---------
186
-        tempAtomA = GibbsSamp.getTotNumAtoms('A');
187
-        tempAtomP = GibbsSamp.getTotNumAtoms('P');
188
-        nAEquil[outCount] = tempAtomA;
189
-        nPEquil[outCount] = tempAtomP;
190
-        outCount++;
191
-
192
-        if (ext_iter % numOutputs == 0 && messages)
193
-        {
194
-            Rcpp::Rcout << "Equil:" << ext_iter << " of " << nEquil <<
195
-                        ", Atoms:" << tempAtomA << "("
196
-                        << tempAtomP << ")" <<
197
-                        "  Chi2 = " << tempChiSq << endl;
198
-        }
199
-
200
-        // -------------------------------------------
201
-        // re-calculate nIterA and nIterP to the expected number of atoms
202
-        nIterA = (unsigned long) Random::poisson(max((double) GibbsSamp.getTotNumAtoms('A'), 10.));
203
-        nIterP = (unsigned long) Random::poisson(max((double) GibbsSamp.getTotNumAtoms('P'), 10.));
204
-        // --------------------------------------------
205
-    }  // end of for-block for equilibration
206
-
207
-    // ===========================================================================
208
-    // Part 2.5) Equilibration Settling:
209
-    // Allow the Equilibration to settle for 10% of the set Equilibrations
210
-    // ===========================================================================
211
-    int nEquilCool = floor(.1 * nEquil);
212
-
213
-    for (unsigned long ext_iter = 1; ext_iter <= nEquilCool; ++ext_iter)
214
-    {
215
-        for (unsigned long iterA = 1; iterA <= nIterA; ++iterA)
216
-        {
217
-            GibbsSamp.update('A');
218
-        }
219
-
220
-        for (unsigned long iterP = 1; iterP <= nIterP; ++iterP)
221
-        {
222
-            GibbsSamp.update('P');
223
-        }
224
-    }
225
-
226
-    // ===========================================================================
227
-    // Part 3) Sampling:
228
-    // After the system equilibriates in Part 2, we sample the systems with an
229
-    // outer loop of nSample iterations. Within each outer loop iteration, A is
230
-    // iterated nIterA times and P is iterated nIterP times. After sampling,
231
-    // we update nIterA and nIterP according to the expected number of atoms in
232
-    // the atomic spaces of A and P respectively.
233
-    // ===========================================================================
234
-    // Initialize Snapshots of A and P
235
-    vector <vector <vector <double> > > ASnap;
236
-    vector <vector <vector <double> > > PSnap;
237
-    unsigned int statindx = 0;
238
-    outCount = 0;
239
-
240
-    for (unsigned long i = 1; i <= nSample; ++i) {
241
-        for (unsigned long iterA = 1; iterA <= nIterA; ++iterA) {
242
-            GibbsSamp.update('A');
243
-        }
244
-
245
-        GibbsSamp.check_atomic_matrix_consistency('A');
246
-
247
-        for (unsigned long iterP = 1; iterP <= nIterP; ++iterP) {
248
-            GibbsSamp.update('P');
249
-        }
250
-
251
-        GibbsSamp.check_atomic_matrix_consistency('P');
252
-
253
-        if (Q_output_atomic == true) {
254
-            GibbsSamp.output_atomicdomain('A', i);
255
-            GibbsSamp.output_atomicdomain('P', i);
256
-        }
257
-
258
-        statindx += 1;
259
-        GibbsSamp.compute_statistics_prepare_matrices(statindx);
260
-        //Do the same as above.
261
-        tempChiSq = GibbsSamp.get_sysChi2();
262
-        chiVect[(nEquil + i) - 1] = tempChiSq;
263
-        // ----------- output computing info ---------
264
-        tempAtomA = GibbsSamp.getTotNumAtoms('A');
265
-        tempAtomP = GibbsSamp.getTotNumAtoms('P');
266
-        nASamp[outCount] = tempAtomA;
267
-        nPSamp[outCount] = tempAtomP;
268
-        outCount++;
269
-
270
-        if (i % numOutputs == 0) {
271
-            if (messages) {
272
-                Rcpp::Rcout << "Samp: " << i << " of " << nSample <<
273
-                            ", Atoms:" << tempAtomA << "("
274
-                            << tempAtomP << ")" <<
275
-                            "  Chi2 = " << tempChiSq << endl;
276
-
277
-                if (i == nSample) {
278
-                    chi2 = 2.*GibbsSamp.cal_logLikelihood();
279
-                    Rcpp::Rcout << " *** Check value of final chi2: " << chi2 << " **** " << endl;
280
-                }
281
-            }
282
-        }
283
-
284
-        if (SampleSnapshots && (i % (nSample / numSnapshots) == 0)) {
285
-            vector <vector <vector <double> > > NormedMats = GibbsSamp.getNormedMatrices();
286
-            ASnap.push_back(NormedMats[0]);
287
-            PSnap.push_back(NormedMats[1]);
288
-        }
289
-
290
-        // -------------------------------------------
291
-        // re-calculate nIterA and nIterP to the expected number of atoms
292
-        nIterA = (unsigned long) Random::poisson(max((double) GibbsSamp.getTotNumAtoms('A'), 10.));
293
-        nIterP = (unsigned long) Random::poisson(max((double) GibbsSamp.getTotNumAtoms('P'), 10.));
294
-        // --------------------------------------------
295
-    }  // end of for-block for Sampling
296
-
297
-    // ===========================================================================
298
-    // Part 4) Calculate statistics:
299
-    // In this final section, we calculate all statistics pertaining to the final
300
-    // sample and check the results.
301
-    // ===========================================================================
302
-    vector<vector <double> > AMeanVector;
303
-    vector<vector <double> > AStdVector;
304
-    vector<vector <double> > PMeanVector;
305
-    vector<vector <double> > PStdVector;
306
-    GibbsSamp.compute_statistics(statindx,
307
-                                 AMeanVector, AStdVector, PMeanVector, PStdVector);          // compute statistics like mean and s.d.
308
-    //CODE FOR CONVERTING ALL THE VECTORS FOR THE DATAFILES INTO NUMERIC VECTORS OR LISTS
309
-    int numRow = AMeanVector.size();
310
-    int numCol = AMeanVector[0].size() ;
311
-    Rcpp::NumericMatrix AMeanMatrix(numRow, numCol) ;
312
-
313
-    for (int i = 0; i < numRow; i++) {
314
-        for (int j = 0; j < numCol; j++) {
315
-            AMeanMatrix(i, j) = AMeanVector[i][j] ;
316
-        }
317
-    }
318
-
319
-    numRow = AStdVector.size();
320
-    numCol = AStdVector[0].size() ;
321
-    Rcpp::NumericMatrix AStdMatrix(numRow, numCol) ;
322
-
323
-    for (int i = 0; i < numRow; i++) {
324
-        for (int j = 0; j < numCol; j++) {
325
-            AStdMatrix(i, j) = AStdVector[i][j] ;
326
-        }
327
-    }
328
-
329
-    numRow = PMeanVector.size();
330
-    numCol = PMeanVector[0].size() ;
331
-    Rcpp::NumericMatrix PMeanMatrix(numRow, numCol) ;
332
-
333
-    for (int i = 0; i < numRow; i++) {
334
-        for (int j = 0; j < numCol; j++) {
335
-            PMeanMatrix(i, j) = PMeanVector[i][j] ;
336
-        }
337
-    }
338
-
339
-    numRow = PStdVector.size();
340
-    numCol = PStdVector[0].size() ;
341
-    Rcpp::NumericMatrix PStdMatrix(numRow, numCol) ;
342
-
343
-    for (int i = 0; i < numRow; i++) {
344
-        for (int j = 0; j < numCol; j++) {
345
-            PStdMatrix(i, j) = PStdVector[i][j] ;
346
-        }
347
-    }
348
-
349
-    //Code for transferring Snapshots in R
350
-    int numSnaps = numSnapshots; //Arbitrary to keep convention
351
-
352
-    if (SampleSnapshots == true) {
353
-        numRow = AMeanVector.size();
354
-        numCol = AMeanVector[0].size() ;
355
-        arma::cube ASnapR(numRow, numCol, numSnaps);
356
-
357
-        for (int k = 0; k < numSnaps; k++) {
358
-            for (int i = 0; i < numRow; ++i) {
359
-                for (int j = 0; j < numCol; ++j) {
360
-                    ASnapR(i, j, k) = ASnap[k][i][j];
361
-                }
362
-            }
363
-        }
364
-
365
-        numRow = PMeanVector.size();
366
-        numCol = PMeanVector[0].size() ;
367
-        arma::cube PSnapR(numRow, numCol, numSnaps);
368
-
369
-        for (int k = 0; k < numSnaps; k++) {
370
-            for (int i = 0; i < numRow; ++i) {
371
-                for (int j = 0; j < numCol; ++j) {
372
-                    PSnapR(i, j, k) = PSnap[k][i][j];
373
-                }
374
-            }
375
-        }
376
-
377
-        Rcpp::List fileContainer =  Rcpp::List::create(Rcpp::Named("Amean") = AMeanMatrix,
378
-                                    Rcpp::Named("Asd") = AStdMatrix, Rcpp::Named("Pmean") = PMeanMatrix, Rcpp::Named("Psd") = PStdMatrix,
379
-                                    Rcpp::Named("ASnapshots") = ASnapR, Rcpp::Named("PSnapshots") = PSnapR,
380
-                                    Rcpp::Named("atomsAEquil") = nAEquil, Rcpp::Named("atomsASamp") = nASamp,
381
-                                    Rcpp::Named("atomsPEquil") = nPEquil, Rcpp::Named("atomsPSamp") = nPSamp, Rcpp::Named("chiSqValues") = chiVect,
382
-                                    Rcpp::Named("randSeed") = seedUsed);
383
-        return (fileContainer);
384
-
385
-    } else {
386
-        //Just leave the snapshots as empty lists
387
-        Rcpp::List ASnapR = Rcpp::List::create();
388
-        Rcpp::List PSnapR = Rcpp::List::create();
389
-        Rcpp::List fileContainer =  Rcpp::List::create(Rcpp::Named("Amean") = AMeanMatrix,
390
-                                    Rcpp::Named("Asd") = AStdMatrix, Rcpp::Named("Pmean") = PMeanMatrix, Rcpp::Named("Psd") = PStdMatrix,
391
-                                    Rcpp::Named("ASnapshots") = ASnapR, Rcpp::Named("PSnapshots") = PSnapR,
392
-                                    Rcpp::Named("atomsAEquil") = nAEquil, Rcpp::Named("atomsASamp") = nASamp,
393
-                                    Rcpp::Named("atomsPEquil") = nPEquil, Rcpp::Named("atomsPSamp") = nPSamp, Rcpp::Named("chiSqValues") = chiVect,
394
-                                    Rcpp::Named("randSeed") = seedUsed);
395
-        return (fileContainer);
396
-    }
397
-}
398 0
deleted file mode 100644
... ...
@@ -1,804 +0,0 @@
1
-
2
-// CoGAPS C++ Version
3
-//
4
-// Functions to compute matrix multiplication and
5
-// Likelihood function
6
-//
7
-// History: v 1.0  Jan 16, 2014
8
-//          Updated to include mapped methods August 7, 2014
9
-
10
-
11
-#include <iostream>
12
-#include <cmath>
13
-#include <stdexcept>
14
-#include "GAPSNorm.h"
15
-// -- checking
16
-#include <iomanip>
17
-// -----
18
-
19
-#include "Matrix.h"
20
-#include "MatAlgo.h"
21
-
22
-using std::vector;
23
-
24
-#define GAPS_SQ(x) ((x) * (x))
25
-
26
-double gaps::norm::calChi2(const Matrix &D, const Matrix &AP, const Matrix &S)
27
-{
28
-    double chi2 = 0;
29
-    for (unsigned i = 0; i < D.nRow(); ++i)
30
-    {
31
-        for (unsigned j = 0; j < D.nCol(); ++j)
32
-        {
33
-            chi2 += GAPS_SQ(D(i,j) - AP(i,j)) / GAPS_SQ(S(i,j));
34
-        }
35
-    }
36
-    return chi2;
37
-}
38
-
39
-// ---------------------------------------------------------------------------
40
-// Calcuation of the change in the log-likelihood when matrix A (or P) is
41
-// augmented with matrix ElemChange (~delA or delP) when there is only ONE
42
-// change. Mathematically, it calculates changes to the following:
43
-//
44
-// Let the only change be delA_{mn}, define M = D-AP;
45
-// delloglikelihood = \sum_j [2*M_{mj}*(delA_{mn}*P_{nj})-(delA_{mn}*P_{nj})^2]/2
46
-//                          / S_{mj}^2
47
-double GAPSNorm::calcDeltaLL1E(char matrix_label, const Matrix &D, const Matrix &S,
48
-const Matrix &A, const Matrix &P, const std::vector<ElementChange> ElemChange,
49
-unsigned int nFactor)
50
-{
51
-    // ------- Calculate changes in log-likelihood --------
52
-    // Extract where and how large the change is
53
-    unsigned int chRow, chCol;
54
-    double delelem;
55
-    double sqTerm = 0;
56
-    double mockTerm = 0.;
57
-    double delloglikelihood = 0.;
58
-    chRow = ElemChange[0].row;    // chRow = iRow that carries a change
59
-    chCol = ElemChange[0].col;    // chCol = iCol that carries a change
60
-    delelem = ElemChange[0].delta;  // delelem = change in the element
61
-
62
-    switch (matrix_label)
63
-    {
64
-        case 'A':
65
-        {
66
-            // ---- Form M = D - A*P, in particular, we need only M[chRow][]
67
-            Vector M(D.nCol());
68
-
69
-            for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
70
-            {
71
-                M(iCol) = D(chRow, iCol);
72
-                for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
73
-                {
74
-                    M(iCol) -= A(chRow,iPattern) * P(iPattern,iCol);
75
-                }
76
-                mockTerm = 2.0 * M(iCol) * delelem * P(chCol,iCol);
77
-                sqTerm = pow(delelem * P(chCol,iCol), 2);
78
-                delloglikelihood += (mockTerm - sqTerm) / 2.0 / pow(S(chRow,iCol), 2);
79
-            }
80
-            break;
81
-        }
82
-
83
-        case 'P':
84
-        {
85
-            // ---- Form M = D - A*P, in particular, we need only M[][chCol]
86
-            Vector M(D.nRow());
87
-
88
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
89
-            {
90
-                M(iRow) = D(iRow,chCol);
91
-                for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
92
-                {
93
-                    M(iRow) -= A(iRow,iPattern) * P(iPattern,chCol);
94
-                }
95
-                mockTerm = 2.0 * M(iRow) * delelem * A(iRow,chRow);
96
-                sqTerm = pow(delelem * A(iRow,chRow), 2);
97
-                delloglikelihood += (mockTerm - sqTerm) / 2.0 / pow(S(iRow,chCol), 2);
98
-            }
99
-            break;
100
-        }
101
-    }
102
-    return delloglikelihood;
103
-}
104
-
105
-
106
-// ---------------------------------------------------------------------------
107
-// Calcuation of the change in the log-likelihood when matrix A (or P) is
108
-// augmented with matrix ElemChange (~delA or delP) when there are TWO changes.
109
-// Mathematically, it calculates changes to the following:
110
-//
111
-// Let the only change be delA_{mn} and delA_{rs}, define M = D-AP;
112
-// Two cases:
113
-// I) m = r
114
-// let Del_j = delA_{mn}*P_{nj} + delA_{ms}*P_{sj}
115
-// delloglikelihood = \sum_j [2*M_{mj}*Del_j-Del_j^2] /2/S_{mj}^2
116
-// II) m != r
117
-// let dellog1 = \sum_j [2*M_{mj}*(delA_{mn}*P_{nj})-(delA_{mn}*P_{nj})^2] /2
118
-//          /S_{mj}^2
119
-//     dellog2 = \sum_j [2*M_{rj}*(delA_{rs}*P_{sj})-(delA_{rs}*P_{sj})^2] /2
120
-//          /S_{rj}^2
121
-//     delloglikelihood = dellog1 + dellog2
122
-
123
-double GAPSNorm::calcDeltaLL2E(char matrix_label, const Matrix &D, const Matrix &S,
124
-const Matrix &A, const Matrix &P, const std::vector<ElementChange> ElemChange,
125
-unsigned int nFactor)
126
-{
127
-    // ------- Calculate changes in log-likelihood --------
128
-    // Extract where and how much the change is.
129
-    unsigned int chRow[2], chCol[2];
130
-    double delelem[2];
131
-    chRow[0] = ElemChange[0].row;  // chRow[0] = iRow for the first change
132
-    chCol[0] = ElemChange[0].col;  // chCol[0] = iCol for the first change
133
-    delelem[0] = ElemChange[0].delta; // delelem[0] = first change
134
-    chRow[1] = ElemChange[1].row;  // chRow[1] = iRow for the second change
135
-    chCol[1] = ElemChange[1].col;  // chCol[1] = iCol for the second change
136
-    delelem[1] = ElemChange[1].delta; // delelem[1] = second change
137
-
138
-    // Calculate delloglikelihood.
139
-    double mockTerm0 = 0.; // temp variables
140
-    double sqTerm0 = 0;    // temp variables
141
-    double mockTerm1 = 0;  // temp variables
142
-    double sqTerm1 = 0;    // temp variables
143
-    unsigned int iRow, iCol, iPattern;  // loop counters
144
-    double delloglikelihood = 0.;    // target quantity to compute
145
-
146
-    switch (matrix_label)
147
-    {
148
-        case 'A':
149
-        {
150
-            // ---- Form M = D - A*P, in particular, we need only M[chRow[0]][]
151
-            // and M[chRow[1]][].
152
-            double M[2][D.nCol()];
153
-
154
-            for (iCol = 0; iCol < D.nCol(); ++ iCol)
155
-            {
156
-                M[0][iCol] = D(chRow[0],iCol);
157
-                M[1][iCol] = D(chRow[1],iCol);
158
-
159
-                for (iPattern = 0; iPattern < nFactor; ++iPattern)
160
-                {
161
-                    M[0][iCol] -= A(chRow[0],iPattern) * P(iPattern,iCol);
162
-                    M[1][iCol] -= A(chRow[1],iPattern) * P(iPattern,iCol);
163
-                }
164
-            }
165
-
166
-            // ---- Two conditions to calculate delloglikelihood.
167
-            if (chRow[0] == chRow[1])
168
-            {
169
-                for (iCol = 0; iCol < D.nCol(); ++iCol)
170
-                {
171
-                    mockTerm0 = 2.0 * M[0][iCol] * (delelem[0] * P(chCol[0],iCol) +
172
-                                                 delelem[1] * P(chCol[1],iCol));
173
-                    sqTerm0 = pow(delelem[0] * P(chCol[0],iCol) + delelem[1] * P(chCol[1],iCol), 2);
174
-                    delloglikelihood += (mockTerm0 - sqTerm0) / 2. / pow(S(chRow[0],iCol), 2);
175
-                }
176
-
177
-            }
178
-            else
179
-            {
180
-                for (iCol = 0; iCol < D.nCol(); ++iCol)
181
-                {
182
-                    mockTerm0 = 2.0 * M[0][iCol] * delelem[0] * P(chCol[0],iCol);
183
-                    sqTerm0 = pow(delelem[0] * P(chCol[0],iCol), 2);
184
-                    mockTerm1 = 2.0 * M[1][iCol] * delelem[1] * P(chCol[1],iCol);
185
-                    sqTerm1 = pow(delelem[1] * P(chCol[1],iCol), 2);
186
-                    delloglikelihood += (mockTerm0 - sqTerm0) / 2.0 / pow(S(chRow[0],iCol), 2) +
187
-                                        (mockTerm1 - sqTerm1) / 2.0 / pow(S(chRow[1],iCol), 2);
188
-                }
189
-            }
190
-            break;
191
-        }
192
-
193
-        case 'P':
194
-        {
195
-            // ---- Form M = D - A*P, in particular, we need only M[][chCol[0]]
196
-            // and M[][chCol[1]].
197
-            double M[D.nRow()][2];
198
-
199
-            for (iRow = 0; iRow < D.nRow(); ++ iRow)
200
-            {
201
-                M[iRow][0] = D(iRow,chCol[0]);
202
-                M[iRow][1] = D(iRow,chCol[1]);
203
-
204
-                for (iPattern = 0; iPattern < nFactor; ++iPattern)
205
-                {
206
-                    M[iRow][0] -= A(iRow,iPattern) * P(iPattern,chCol[0]);
207
-                    M[iRow][1] -= A(iRow,iPattern) * P(iPattern,chCol[1]);
208
-                }
209
-            }
210
-
211
-            // ------ Two conditions to calculate the delloglikelihood.
212
-            if (chCol[0] == chCol[1])
213
-            {
214
-                for (iRow = 0; iRow < D.nRow(); ++iRow)
215
-                {
216
-                    mockTerm0 = 2.0 * M[iRow][0] * (A(iRow,chRow[0]) * delelem[0] +
217
-                        A(iRow,chRow[1]) * delelem[1]);
218
-                    sqTerm0 = pow(A(iRow,chRow[0]) * delelem[0] + A(iRow,chRow[1]) * delelem[1], 2);
219
-                    delloglikelihood += (mockTerm0 - sqTerm0) / 2. / pow(S(iRow,chCol[0]), 2);
220
-                }
221
-
222
-            }
223
-            else
224
-            {
225
-                for (iRow = 0; iRow < D.nRow(); ++iRow)
226
-                {
227
-                    mockTerm0 = 2.0 * M[iRow][0] * A(iRow,chRow[0]) * delelem[0];
228
-                    sqTerm0 = pow(A(iRow,chRow[0]) * delelem[0], 2);
229
-                    mockTerm1 = 2.0 * M[iRow][1] * A(iRow,chRow[1]) * delelem[1];
230
-                    sqTerm1 = pow(A(iRow,chRow[1]) * delelem[1], 2);
231
-                    delloglikelihood += (mockTerm0 - sqTerm0) / 2.0 / pow(S(iRow,chCol[0]), 2) +
232
-                                        (mockTerm1 - sqTerm1) / 2.0 / pow(S(iRow,chCol[1]), 2);
233
-                }
234
-            }
235
-            break;
236
-        }
237
-    }
238
-    return delloglikelihood;
239
-}
240
-
241
-// ---------------------------------------------------------------------------
242
-// Calcuation of the change in the log-likelihood when matrix A (or P) is augmented
243
-// with matrix ElemChange (~delA or delP) generally.
244
-// Mathematically, it calculates changes to the following:
245
-//
246
-// Let the only change be delA, define M = D-AP;
247
-// delloglikelihood = \sum_{ij} [2*M_{ij}*(delA*P)_{ij}-(delA*P)_{ij}^2] /2/S_{ij}^2
248
-
249
-double GAPSNorm::calcDeltaLLGen(char matrix_label, const Matrix &D, const Matrix &S,
250
-const Matrix &A, const Matrix &P, const std::vector<ElementChange> ElemChange,
251
-unsigned int nFactor)
252
-{
253
-    double delloglikelihood = 0;
254
-    // ---- Form M = D - A*P -----
255
-    Matrix M(D.nRow(), D.nCol());
256
-
257
-    for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
258
-    {
259
-        for (unsigned int iCol = 0; iCol < D.nCol(); ++ iCol)
260
-        {
261
-            M(iRow,iCol) = D(iRow,iCol);
262
-            for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
263
-            {
264
-                M(iRow,iCol) -= A(iRow,iPattern) * P(iPattern,iCol);
265
-            }
266
-        }
267
-    }
268
-
269
-    switch (matrix_label)
270
-    {
271
-        case 'A':
272
-        {
273
-            unsigned int chRow, chCol;
274
-            double delelem;
275
-            unsigned int sizeChange = ElemChange.size();
276
-            // ---- Construct delA -----
277
-            double delA[D.nRow()][nFactor];
278
-
279
-            for (unsigned int m = 0; m < D.nRow(); m++)
280
-            {
281
-                for (unsigned int n = 0; n < nFactor ; n++)
282
-                {
283
-                    delA[m][n] = 0.0;
284
-                }
285
-            }
286
-
287
-            for (unsigned int m = 0; m < sizeChange; ++ m)
288
-            {
289
-                chRow = ElemChange[m].row;
290
-                chCol = ElemChange[m].col;
291
-                delelem = ElemChange[m].delta;
292
-                delA[chRow][chCol] += delelem;
293
-            }
294
-
295
-            // ----- Compute delA*P ----
296
-            double delAP[D.nRow()][D.nCol()];
297
-
298
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
299
-            {
300
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
301
-                {
302
-                    delAP[iRow][iCol] = 0.;
303
-                    for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
304
-                    {
305
-                        delAP[iRow][iCol] += delA[iRow][iPattern] * P(iPattern,iCol);
306
-                    }
307
-                }
308
-            }
309
-
310
-            // ------ Compute delloglikelihood -------
311
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
312
-            {
313
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
314
-                {
315
-                    delloglikelihood += (2.0 * M(iRow,iCol) * delAP[iRow][iCol] - pow(delAP[iRow][iCol], 2))
316
-                                        / 2.0 / pow(S(iRow,iCol), 2);
317
-                }
318
-            }
319
-            break;
320
-        } // end of the A sub-block
321
-
322
-        case 'P':
323
-        {
324
-            unsigned int chRow, chCol;
325
-            double delelem;
326
-            unsigned int sizeChange = ElemChange.size();
327
-            // ---- Construct delP -----
328
-            double delP[nFactor][D.nCol()];
329
-
330
-            for (unsigned int m = 0; m < nFactor; m++)
331
-            {
332
-                for (unsigned int n = 0; n < D.nCol() ; n++)
333
-                {
334
-                    delP[m][n] = 0. ;
335
-                }
336
-            }
337
-
338
-            for (unsigned int m = 0; m < sizeChange; ++ m)
339
-            {
340
-                chRow = ElemChange[m].row;
341
-                chCol = ElemChange[m].col;
342
-                delelem = ElemChange[m].delta;
343
-                delP[chRow][chCol] += delelem;
344
-            }
345
-
346
-            // ----- Compute A*delP ----
347
-            double AdelP[D.nRow()][D.nCol()];
348
-
349
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
350
-            {
351
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
352
-                {
353
-                    AdelP[iRow][iCol] = 0.;
354
-                    for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
355
-                    {
356
-                        AdelP[iRow][iCol] += A(iRow,iPattern) * delP[iPattern][iCol];
357
-                    }
358
-                }
359
-            }
360
-
361
-            // ------ Compute delloglikelihood -------
362
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
363
-            {
364
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
365
-                {
366
-                    delloglikelihood += (2.0 * M(iRow,iCol) * AdelP[iRow][iCol] - pow(AdelP[iRow][iCol], 2))
367
-                                        / 2.0 / pow(S(iRow,iCol), 2);
368
-                }
369
-            }
370
-            break;
371
-        } // end of the P sub-block
372
-    } // end of switch
373
-    return delloglikelihood;
374
-} // end of calcDeltaGen
375
-
376
-// ---------------------------------------------------------------------------
377
-// Calcuation of the change in the log-likelihood when matrix A (or P) is augmented
378
-// with a change that occurs in a fixed pattern. This means that an entire
379
-// column/row of A or P is changed.
380
-// Because changing an entire column of A or row of P completely changes the matrix
381
-// A*P, we recalculate A*P completely
382
-double GAPSNorm::calcDeltaLLMap(char matrix_label, const Matrix &D, const Matrix &S,
383
-const Matrix &A, const Matrix &P, std::vector<double> &newPat, unsigned int chPat,
384
-unsigned int nFactor)
385
-{
386
-    double delloglikelihood = 0;
387
-    // ---- Form M = D - A*P -----
388
-    Matrix M(D.nRow(), D.nCol());
389
-
390
-    for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
391
-    {
392
-        for (unsigned int iCol = 0; iCol < D.nCol(); ++ iCol)
393
-        {
394
-            M(iRow,iCol) = D(iRow,iCol);
395
-            for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
396
-            {
397
-                M(iRow,iCol) -= A(iRow,iPattern) * P(iPattern,iCol);
398
-            }
399
-        }
400
-    }
401
-
402
-    switch (matrix_label)
403
-    {
404
-        case 'A':
405
-        {
406
-            // ---- Construct delA -----
407
-            double delA[D.nRow()][nFactor];
408
-
409
-            for (unsigned int m = 0; m < D.nRow(); m++)
410
-            {
411
-                for (unsigned int n = 0; n < nFactor ; n++)
412
-                {
413
-                    delA[m][n] = 0. ;
414
-                }
415
-            }
416
-
417
-            // The change here is along one column, chPat:
418
-            // Remember that the passed in vector represents
419
-            // the actual new row, not the change, hence the
420
-            // adjustment.
421
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
422
-            {
423
-                delA[iRow][chPat] += (newPat.at(iRow) - A(iRow,chPat));
424
-            }
425
-
426
-            // ----- Compute delA*P ----
427
-            double delAP[D.nRow()][D.nCol()];
428
-
429
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
430
-            {
431
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
432
-                {
433
-                    delAP[iRow][iCol] = 0.;
434
-                    for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
435
-                    {
436
-                        delAP[iRow][iCol] += delA[iRow][iPattern] * P(iPattern,iCol);
437
-                    }
438
-                }
439
-            }
440
-
441
-            // ------ Compute delloglikelihood -------
442
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
443
-            {
444
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
445
-                {
446
-                    delloglikelihood += (2.0 * M(iRow,iCol) * delAP[iRow][iCol] - pow(delAP[iRow][iCol], 2))
447
-                                        / 2.0 / pow(S(iRow,iCol), 2);
448
-                }
449
-            }
450
-
451
-            break;
452
-        } // end of the A sub-block
453
-
454
-        case 'P':
455
-        {
456
-            // ---- Construct delP -----
457
-            double delP[nFactor][D.nCol()];
458
-
459
-            for (unsigned int m = 0; m < nFactor; m++)
460
-            {
461
-                for (unsigned int n = 0; n < D.nCol() ; n++)
462
-                {
463
-                    delP[m][n] = 0. ;
464
-                }
465
-            }
466
-
467
-            // The change here is along one row, chPat:
468
-            // Remember that the passed in vector represents
469
-            // the actual new row, not the change, hence the
470
-            // adjustment.
471
-            for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
472
-            {
473
-                delP[chPat][iCol] += (newPat.at(iCol) - P(chPat,iCol));
474
-            }
475
-
476
-            // ----- Compute A*delP ----
477
-            double AdelP[D.nRow()][D.nCol()];
478
-
479
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
480
-            {
481
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
482
-                {
483
-                    AdelP[iRow][iCol] = 0.;
484
-                    for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
485
-                    {
486
-                        AdelP[iRow][iCol] += A(iRow,iPattern) * delP[iPattern][iCol];
487
-                    }
488
-                }
489
-            }
490
-
491
-            // ------ Compute delloglikelihood -------
492
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
493
-            {
494
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
495
-                {
496
-                    delloglikelihood += (2.0 * M(iRow,iCol) * AdelP[iRow][iCol] - pow(AdelP[iRow][iCol], 2))
497
-                                        / 2. / pow(S(iRow,iCol), 2);
498
-                }
499
-            }
500
-            break;
501
-        } // end of the P sub-block
502
-    } // end of switch
503
-    return delloglikelihood;
504
-} // end of calcDeltaLLGenMap
505
-
506
-// For move exchange / multiple pattern changes
507
-double GAPSNorm::calcDeltaLL2Map(char matrix_label, const Matrix &D, const Matrix &S,
508
-const Matrix &A, const Matrix &P,
509
-std::vector<double> &newPat1, unsigned int chPat1,
510
-std::vector <double> &newPat2, unsigned int chPat2,
511
-unsigned int nFactor)
512
-{
513
-    double delloglikelihood = 0;
514
-    // ---- Form M = D - A*P -----
515
-    Matrix M(D.nRow(), D.nCol());
516
-
517
-    for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
518
-    {
519
-        for (unsigned int iCol = 0; iCol < D.nCol(); ++ iCol)
520
-        {
521
-            M(iRow,iCol) = D(iRow,iCol);
522
-            for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
523
-            {
524
-                M(iRow,iCol) -= A(iRow,iPattern) * P(iPattern,iCol);
525
-            }
526
-        }
527
-    }
528
-
529
-    switch (matrix_label)
530
-    {
531
-        case 'A':
532
-        {
533
-            // ---- Construct delA -----
534
-            double delA[D.nRow()][nFactor];
535
-
536
-            for (unsigned int m = 0; m < D.nRow(); m++)
537
-            {
538
-                for (unsigned int n = 0; n < nFactor ; n++)
539
-                {
540
-                    delA[m][n] = 0. ;
541
-                }
542
-            }
543
-
544
-            // The change here is along 2 columns, chPat1, and chPat2:
545
-            // Remember that the passed in vector represents
546
-            // the actual new row, not the change, hence the
547
-            // adjustment.
548
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
549
-            {
550
-                delA[iRow][chPat1] += (newPat1.at(iRow) - A(iRow,chPat1));
551
-                delA[iRow][chPat2] += (newPat2.at(iRow) - A(iRow,chPat2));
552
-            }
553
-
554
-            // ----- Compute delA*P ----
555
-            double delAP[D.nRow()][D.nCol()];
556
-
557
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
558
-            {
559
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
560
-                {
561
-                    delAP[iRow][iCol] = 0.;
562
-                    for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
563
-                    {
564
-                        delAP[iRow][iCol] += delA[iRow][iPattern] * P(iPattern,iCol);
565
-                    }
566
-                }
567
-            }
568
-
569
-            // ------ Compute delloglikelihood -------
570
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
571
-            {
572
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
573
-                {
574
-                    delloglikelihood += (2.0 * M(iRow,iCol) * delAP[iRow][iCol] - pow(delAP[iRow][iCol], 2))
575
-                                        / 2.0 / pow(S(iRow,iCol), 2);
576
-                }
577
-            }
578
-            break;
579
-        } // end of the A sub-block
580
-
581
-        case 'P':
582
-        {
583
-            // ---- Construct delP -----
584
-            double delP[nFactor][D.nCol()];
585
-
586
-            for (unsigned int m = 0; m < nFactor; m++)
587
-            {
588
-                for (unsigned int n = 0; n < D.nCol() ; n++)
589
-                {
590
-                    delP[m][n] = 0.0;
591
-                }
592
-            }
593
-
594
-            // The change here is along one row, chPat:
595
-            // Remember that the passed in vector represents
596
-            // the actual new row, not the change, hence the
597
-            // adjustment.
598
-            for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
599
-            {
600
-                delP[chPat1][iCol] += (newPat1.at(iCol) - P(chPat1,iCol));
601
-                delP[chPat2][iCol] += (newPat2.at(iCol) - P(chPat2,iCol));
602
-            }
603
-
604
-            // ----- Compute A*delP ----
605
-            double AdelP[D.nRow()][D.nCol()];
606
-
607
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
608
-            {
609
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
610
-                {
611
-                    AdelP[iRow][iCol] = 0.;
612
-                    for (unsigned int iPattern = 0; iPattern < nFactor; ++iPattern)
613
-                    {
614
-                        AdelP[iRow][iCol] += A(iRow,iPattern) * delP[iPattern][iCol];
615
-                    }
616
-                }
617
-            }
618
-
619
-            // ------ Compute delloglikelihood -------
620
-            for (unsigned int iRow = 0; iRow < D.nRow(); ++iRow)
621
-            {
622
-                for (unsigned int iCol = 0; iCol < D.nCol(); ++iCol)
623
-                {
624
-                    delloglikelihood += (2.0 * M(iRow,iCol) * AdelP[iRow][iCol] - pow(AdelP[iRow][iCol], 2))
625
-                                        / 2.0 / pow(S(iRow,iCol), 2);
626
-                }
627
-            }
628
-            break;
629
-        } // end of the P sub-block
630
-    } // end of switch
631
-    return delloglikelihood;
632
-} // end of calcDeltaLLGenMap
633
-
634
-
635
-// ---------------------------------------------------------------------------
636
-// Calcuation of the parameters s and su for exchange action
637
-// will be used GibbsSampler to determine whether or not Gibbs exchange
638
-// action is possible. If it is, s and su are adjusted for annealingtemp
639
-// and used to determine distribution on alpha. For full proof, see Fertig (2009)
640
-
641
-std::pair<double, double> GAPSNorm::calcAlphaParameters(char the_matrix_label, unsigned int nFactor,
642
-const Matrix &D, const Matrix &S, const Matrix &AOrig, const Matrix &POrig,
643
-unsigned int iGene1, unsigned int iPattern1,
644
-unsigned int iGene2, unsigned int iPattern2,
645
-unsigned int iSample1, unsigned int iSample2)
646
-{
647
-    double s = 0.0;
648
-    double su = 0.0;
649
-    double mock = 0.0;
650
-    double mock1 = 0.0;
651
-    double mock2 = 0.0;
652
-
653
-    switch (the_matrix_label)
654
-    {
655
-        case 'A': {
656
-            // ---------- EXCHANGE ACTION WITH A ----------------------------
657
-            double Aeff;
658
-
659
-            // compute the distribution parameters
660
-            if (iGene1 == iGene2)
661
-            {
662
-                for (int jSample = 0; jSample < D.nCol(); jSample++)
663
-                {
664
-                    // Calculate the mock term
665
-                    mock = D(iGene1,jSample);
666
-
667
-                    for (int jPattern = 0; jPattern < nFactor; jPattern++)
668
-                    {
669
-                        Aeff = AOrig(iGene1,jPattern);
670
-                        mock -= Aeff * POrig(jPattern,jSample);
671
-                    } // end of for-block that make changes to elements in A
672
-
673
-                    s += ((POrig(iPattern1,jSample) - POrig(iPattern2,jSample)) *
674
-                          (POrig(iPattern1,jSample) - POrig(iPattern2,jSample))) /
675
-                         pow(S(iGene1,jSample), 2);
676
-                    su += mock * (POrig(iPattern1,jSample) - POrig(iPattern2,jSample)) /
677
-                          pow(S(iGene1,jSample), 2);
678
-                }
679
-            }  // end of case iGene1 = iGene2
680
-
681
-            else // iGene1 != iGene2
682
-            {
683
-                for (int jSample = 0; jSample < D.nCol(); jSample++)
684
-                {
685
-                    mock1 = D(iGene1,jSample);
686
-                    mock2 = D(iGene2,jSample);
687
-
688
-                    for (int jPattern = 0; jPattern < nFactor; jPattern++)
689
-                    {
690
-                        Aeff = AOrig(iGene1,jPattern);
691
-                        mock1 -= Aeff * POrig(jPattern,jSample);
692
-                        Aeff = AOrig(iGene2,jPattern);
693
-                        mock2 -= Aeff * POrig(jPattern,jSample);
694
-                    } // end of for-block for adding changes to A
695
-
696
-                    s  += pow(POrig(iPattern1,jSample) / S(iGene1,jSample), 2) +
697
-                          pow(POrig(iPattern2,jSample) / S(iGene2,jSample), 2);
698
-                    su += mock1 * POrig(iPattern1,jSample) / pow(S(iGene1,jSample), 2) -
699
-                          mock2 * POrig(iPattern2,jSample) / pow(S(iGene2,jSample), 2);
700
-                }
701
-            } // end of if-block for calculating s and su for case 'A'
702
-            break;
703
-        } // end of switch block for EXCHANGE ACTION with A
704
-
705
-        case 'P':
706
-        {
707
-            // ---------- EXCHANGE ACTION WITH P ----------------------------
708
-            double Peff;
709
-
710
-            // compute the distribution parameters
711
-            if (iSample1 == iSample2)
712
-            {
713
-                for (int jGene = 0; jGene < D.nRow(); jGene++)
714
-                {
715
-                    // Calculate the mock term
716
-                    mock = D(jGene,iSample1);
717
-
718
-                    for (int jPattern = 0; jPattern < nFactor; jPattern++)
719
-                    {
720
-                        Peff = POrig(jPattern,iSample1);
721
-                        mock -=  AOrig(jGene,jPattern) * Peff;
722
-                    } // end of for-block that make changes to elements in P
723
-
724
-                    s += pow(((AOrig(jGene,iPattern1) - AOrig(jGene,iPattern2)) /
725
-                              S(jGene,iSample1)), 2);
726
-                    su += mock * (AOrig(jGene,iPattern1) - AOrig(jGene,iPattern2)) /
727
-                          pow(S(jGene,iSample1), 2);
728
-                }
729
-            }  // end of case iSample1 = iSample2
730
-
731
-            else // iSample1 != iSample2
732
-            {
733
-                for (int jGene = 0; jGene < D.nRow(); jGene++)
734
-                {
735
-                    mock1 = D(jGene,iSample1);
736
-                    mock2 = D(jGene,iSample2);
737
-
738
-                    for (int jPattern = 0; jPattern < nFactor; jPattern++)
739
-                    {
740
-                        Peff = POrig(jPattern,iSample1);
741
-                        mock1 -= AOrig(jGene,jPattern) * Peff;
742
-                        Peff = POrig(jPattern,iSample2);
743
-                        mock2 -= AOrig(jGene,jPattern) * Peff;
744
-                    } // end of for-block for adding changes to P
745
-
746
-                    s  += pow(AOrig(jGene,iPattern1) / S(jGene,iSample1), 2) +
747
-                          pow(AOrig(jGene,iPattern2) / S(jGene,iSample2), 2);
748
-                    su += mock1 * AOrig(jGene,iPattern1) / pow(S(jGene,iSample1), 2) -
749
-                          mock2 * AOrig(jGene,iPattern2) / pow(S(jGene,iSample2), 2);
750
-                }
751
-            } // end of if-block for calculating s and su for case 'P'
752
-
753
-            // end of compute distribution parameters for P
754
-            break;
755
-        } // end of switch block for EXCHANGE ACTION with P
756
-    } // end of switch block for EXCHANGE ACTION
757
-    return std::pair<double, double>(s, su);
758
-} // end of calcalphaparameters method
759
-
760
-vector <vector <vector <double> > > GibbsSampler::getNormedMatrices()
761
-{
762
-    vector <vector <double> > AMatrixNormed;
763
-    AMatrixNormed.resize(_nRow, vector<double>(_nFactor, 0.0));
764
-    vector <vector <double> > PMatrixNormed;
765
-    PMatrixNormed.resize(_nFactor, vector<double>(_nCol, 0.0));
766
-    vector<double> k(_nFactor);  // normalized vector
767
-
768
-    // compute the normalization vector
769
-    for (int m = 0; m < _nFactor; ++m)
770
-    {
771
-        k[m] = 0.;
772
-
773
-        for (int n = 0; n < _nCol; ++n)
774
-        {
775
-            k[m] += _PMatrix(m,n);
776
-        }
777
-
778
-        if (k[m] == 0) // when the whole row of P is zero, then don't do anything
779
-        {
780
-            k[m] = 1.0;
781
-        }
782
-    }
783
-
784
-    for (int m = 0; m < _nRow; ++m)
785
-    {
786
-        for (int n = 0; n < _nFactor; ++n)
787
-        {
788
-            AMatrixNormed[m][n] = _AMatrix(m,n) * k[n];
789
-        }
790
-    }
791
-
792
-    for (int m = 0; m < _nFactor; ++m)
793
-    {
794
-        for (int n = 0; n < _nCol; ++n)
795
-        {
796
-            PMatrixNormed[m][n] = _PMatrix(m,n) / k[m];
797
-        }
798
-    }
799
-
800
-    vector <vector <vector <double> > > NormedMatrices;
801
-    NormedMatrices.push_back(AMatrixNormed);
802
-    NormedMatrices.push_back(PMatrixNormed);
803
-    return NormedMatrices;
804
-}
805 0
\ No newline at end of file
806 1
deleted file mode 100644
... ...
@@ -1,38 +0,0 @@
1
-#ifndef __COGAPS_GAPSNORM_H__
2
-#define __COGAPS_GAPSNORM_H__
3
-
4
-#include "Matrix.h"
5
-
6
-#include <iostream>
7
-#include <vector>
8
-#include <utility>
9
-
10
-namespace gaps
11
-{
12
-
13
-namespace norm
14
-{
15
-    double calChi2(const Matrix &D, const Matrix &S, const Matrix &A,
16
-        const Matrix &P, unsigned int nFactor);
17
-
18
-    double calcDeltaLL1E(const Matrix &D, const Matrix &S, const Matrix &AP,
19
-        const MatrixChange &change);
20
-
21
-    double calcDeltaLL2E(const Matrix &D, const Matrix &S, const Matrix &AP,
22
-        const MatrixChange &change);
23
-                                
24
-    /** NEW METHOD
25
-    *   @short Calculate the parameters involved in exchange action for further use
26
-    *   @return a pair <s, su> of alpha distribution parameters for further calculations.
27
-    */
28
-    /*std::pair<double, double> calcAlphaParameters(char the_matrix_label, unsigned int nFactor,
29
-        const Matrix &D, const Matrix &S, const Matrix &AOrig, const Matrix &POrig,
30
-        unsigned int iGene1, unsigned int iPattern1,
31
-        unsigned int iGene2, unsigned int iPattern2,
32
-        unsigned int iSample1, unsigned int iSample2);*/
33
-
34
-} // namespace norm
35
-
36
-} // namespace gaps
37
-
38
-#endif
... ...
@@ -1,276 +1,244 @@
1 1
 #include "GibbsSampler.h"
2
+#include "Algorithms.h"
2 3
 
3 4
 #include <Rcpp.h>
4 5
 
5
-#if 0
6
+#define EPSILON 1E-10
6 7
 
7 8
 GibbsSampler::GibbsSampler(Rcpp::NumericMatrix D, Rcpp::NumericMatrix S,
8
-unsigned int nFactor, double alphaA, double alphaP)
9
+unsigned int nFactor, double alphaA, double alphaP, double maxGibbsMassA,
10
+double maxGibbsMassP, bool singleCellRNASeq)
9 11
     :
10 12
 mDMatrix(D), mSMatrix(S), mAMatrix(D.nrow(), nFactor),
11 13
 mPMatrix(nFactor, D.ncol()), mAPMatrix(D.nrow(), D.ncol()),
12
-mADomain(D.nrow(), nFactor), mPDomain(nFactor, D.ncol())
14
+mADomain('A', D.nrow(), nFactor), mPDomain('P', nFactor, D.ncol()),
15
+mMaxGibbsMassA(maxGibbsMassA), mMaxGibbsMassP(maxGibbsMassP),
16
+mAnnealingTemp(0.0), mChi2(0.0), mSingleCellRNASeq(singleCellRNASeq)
13 17
 {
14
-    double meanD = mSingleCellRNASeq ? gaps::mat_algo::nonZeroMean(mDMatrix)
15
-        : gaps::mat_algo::mean(mDMatrix);
18
+    double meanD = mSingleCellRNASeq ? gaps::algo::nonZeroMean(mDMatrix)
19
+        : gaps::algo::mean(mDMatrix);
16 20
 
17 21
     mADomain.setAlpha(alphaA);
18
-    mADomain.setLambda(alphaA * sqrt(nFactor / meanD) * lambdaScaleFactorA);
22
+    mADomain.setLambda(alphaA * sqrt(nFactor / meanD));
19 23
     mPDomain.setAlpha(alphaP);
20
-    mPDomain.setLambda(alphaP * sqrt(nFactor / meanD) * lambdaScaleFactorP);
24
+    mPDomain.setLambda(alphaP * sqrt(nFactor / meanD));
21 25
 }
22 26
 
23
-double GibbsSampler::getGibbsMass(char matrixLabel, double origMass,
24
-unsigned int iRow, unsigned int iCol, const Matrix &otherMatrix,
25
-const Matrix &currentChainMatrix, const Matrix &D, const Matrix &S,
26
-double rng)
27
+double GibbsSampler::getGibbsMass(const MatrixChange &change)
27 28
 {
28
-    double s  = 0.0;
29
-    double su = 0.0;
30
-
31
-    if (matrixLabel == 'A')
32
-    {
33
-        double psq = 0.0, ssq = 0.0;
34
-        for (unsigned j = 0; j < mDMatrix.nCol(); ++j)
35
-        {
36
-            psq = mPMatrix(col,j) * mPMatrix(col,j);
37
-            ssq = mSMatrix(row,j) * mSMatrix(row,j);
38
-           
39
-            s += (mAnnealingTemp * psq / ssq) / 2;
40
-            su += s * (mDMatrix(row,j) - mAPMatrix(row,j));
41
-        }
42
-    }
43
-    else
44
-    {
45
-        double asq = 0.0, ssq = 0.0;
46
-        for (unsigned i = 0; i < mDMatrix.nRow(); ++i)
47
-        {
48
-            asq = mAMatrix(i,row) * mAMatrix(i,row);
49
-            ssq = mSMatrix(i,col) * mSMatrix(i,col);
50
-           
51
-            s += (mAnnealingTemp * asq / ssq) / 2;
52
-            su += s * (mDMatrix(i,col) - mAPMatrix(i,col));
53
-        }