Browse code

still unsynced

sherman5 authored on 17/04/2018 23:22:57
Showing4 changed files

... ...
@@ -30,7 +30,7 @@ const std::string &cptFile, unsigned pumpThreshold, unsigned nPumpSamples)
30 30
     // create internal state from parameters and run from there
31 31
     GapsRunner runner(D, S, nFactor, nEquil, nEquilCool, nSample,
32 32
         nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seedUsed,
33
-        messages, singleCellRNASeq,  checkpointInterval, cptFile,
33
+        messages, singleCellRNASeq, checkpointInterval, cptFile,
34 34
         whichMatrixFixed, FP);
35 35
     return runner.run();
36 36
 }
... ...
@@ -138,7 +138,7 @@ void GapsRunner::updateSampler()
138 138
     mASampler.update(mIterA);
139 139
     mPSampler.sync(mASampler);
140 140
 
141
-    mNumUpdatesA += mIterP;
141
+    mNumUpdatesP += mIterP;
142 142
     mPSampler.update(mIterP);
143 143
     mASampler.sync(mPSampler);
144 144
 }
... ...
@@ -234,7 +234,8 @@ void GibbsSampler<T, MatA, MatB>::death(uint64_t pos, float mass, unsigned row,
234 234
 unsigned col)
235 235
 {
236 236
     removeMass(pos, mass, row, col);
237
-    mass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col, -mass) : mass;
237
+    float newMass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col, -mass) : 0.f;
238
+    mass = newMass < gaps::algo::epsilon ? mass : newMass;
238 239
     float deltaLL = impl()->computeDeltaLL(row, col, mass);
239 240
     if (deltaLL * mAnnealingTemp > std::log(gaps::random::uniform()))
240 241
     {
... ...
@@ -19,10 +19,15 @@ void ProposalQueue::setAlpha(float alpha)
19 19
 
20 20
 float ProposalQueue::deathProb(unsigned nAtoms) const
21 21
 {
22
-    double size = static_cast<double>(mDomainSize);
23
-    double term1 = (size - static_cast<double>(nAtoms)) / size;
24
-    double term2 = mAlpha * static_cast<double>(mNumBins) * term1;
25
-    return static_cast<double>(nAtoms) / (static_cast<double>(nAtoms) + term2);
22
+    //double size = static_cast<double>(mDomainSize);
23
+    //double term1 = (size - static_cast<double>(nAtoms)) / size;
24
+    //double term2 = mAlpha * static_cast<double>(mNumBins) * term1;
25
+    //return static_cast<double>(nAtoms) / (static_cast<double>(nAtoms) + term2);
26
+    float dMax = (float)mDomainSize;
27
+    float dNum = (float)nAtoms;
28
+    float maxTerm = (dMax - dNum) / dMax;
29
+
30
+    return dNum / (dNum + mAlpha * (float)mNumBins * maxTerm);
26 31
 }
27 32
 
28 33
 AtomicProposal ProposalQueue::makeProposal(const AtomicDomain &domain)
... ...
@@ -33,13 +38,13 @@ AtomicProposal ProposalQueue::makeProposal(const AtomicDomain &domain)
33 38
         return birth(domain);
34 39
     }
35 40
 
36
-    float bdProb = domain.size() < 2 ? 0.6667f : 0.5f;
37
-    float u = gaps::random::uniform();
38
-    if (u <= bdProb)
39
-    {
41
+    //float bdProb = domain.size() < 2 ? 0.6667f : 0.5f;
42
+    //float u = gaps::random::uniform();
43
+    //if (u <= bdProb)
44
+    //{
40 45
         return gaps::random::uniform() < deathProb(domain.size()) ? 
41 46
             death(domain) : birth(domain);
42
-    }
47
+    /*}
43 48
     else if (u < 0.75f || domain.size() < 2)
44 49
     {
45 50
         return move(domain);
... ...
@@ -47,7 +52,7 @@ AtomicProposal ProposalQueue::makeProposal(const AtomicDomain &domain)
47 52
     else
48 53
     {
49 54
         return exchange(domain);
50
-    }
55
+    }*/
51 56
 }
52 57
 
53 58
 AtomicProposal ProposalQueue::birth(const AtomicDomain &domain)