Repository: bab2min/tomotopy (language: python)

Issue:
Error in CTM
I am working on a topic modelling project and want to use the Tomotopy implementation of CTM. So far I have succesfully implemented LDA and Pachinko in the same pipeline as CTM but for CTM I am getting the following error. I am uncertain of why this error arises or how to resolve this. We have tried changing the number of topics, the values for eta and introduced 'model.num_beta_sample=5'. I have around 4620 documents for training. Please find the relevant bit of the code below the error.

![image](https://user-images.githubusercontent.com/56274607/97922379-a78a7d80-1d5c-11eb-8756-bdcbd1e4bcdb.png)

`def train_CTM(documents, k, min_cf=0, rm_top=0, eta=0.01, model_burn_in=100, #smoothing_alpha=0.1,
              train_updates = 1000, train_iter = 10):
    
    # instantiate
    model = tp.CTModel(tw=tp.TermWeight.ONE, min_cf=min_cf, rm_top=rm_top, k=k, #smoothing_alpha = smoothing_alpha,
                      eta = eta)

    model.num_beta_sample = 5

    # add documents to model
    for doc in documents: model.add_doc(doc)
    
    # training**
    model.burn_in = model_burn_in
    # initialising 
    model.train(iter=0)
    print('Num docs:', len(model.docs), ', Vocab size:', len(model.used_vocabs), ', Num words:', model.num_words)
    print('Removed top words:', model.removed_top_words)
    print('Training...')
    # actual training 
    for i in range(0, train_updates, train_iter):
        model.train(train_iter)
        if i%100==0:print('Iteration: {}'.format(i))
    return model

def CTM_search(train_docs, test_docs):
    """ Wrapper function for CTM grid search """
    print("Starting CTM grid search ...")

    w = get_num_tokens(train_docs)

    # The GRID
    ks = [50, 100, 150, 200, 300]
    etas = [1/w, 10/w, 0.1/w]
    

    grid_results = pd.DataFrame(columns=["K","eta","perplexity"])

    start = time.time()

    for k in ks:
        startk = time.time()
        for eta in etas:
            print(" ------------------------------------------- ")
            print("| K = " + str(k) + "\t| eta=" + str(eta))
            try:
                model= train_CTM(train_docs, k,
                                eta=eta, 
                                model_burn_in=MODEL_BURN_IN,
                                train_updates = TRAIN_UPDATES, 
                                train_iter = TRAIN_ITERS)
            except Exception as e:
                print("The following error occurred:")
                print(e)
                print("The grid search will skip the current values.")
                grid_results = grid_results.append({"K":k, "eta":eta,"perplexity":np.nan}, 
                                                ignore_index=True)
                continue
            

            ll = get_test_LL(test_docs, model)
            pp = compute_test_pp(ll, test_docs)
            grid_results = grid_results.append({"K":k, "eta":eta,"perplexity":pp}, 
                                                ignore_index=True)
            print("Test perplexity = "+str(pp))

        endk = time.time()
        print("Time elapsed for k="+str(k)+" : " + str(round(endk - startk, 1)) + " s")

    end = time.time()
    print("Time elapsed: " + str(round(end - start, 1)) + " s")
    
    return grid_results`


Would you have any indication of where this may originate from or how to fix it?


Many thanks in advance for any insight you can provide.


Discussion / hints:
Is there any news on this one?
Thank @danielgarel for reporting the bug and @GiuliaLanzillotta for reminding it.
Umm... it seems there are some potential bugs in the implementation of Truncated Multivariate Normal distribution.

https://github.com/bab2min/tomotopy/blob/e95b656c46c0d486c09f64a7d7718ce0195ed326/src/Utils/TruncMultiNormal.hpp#L10-L69

Normally, truncation domain `[A, B]` satisfies `A <= B`, but sometimes it seems to be broken.
I'll test it with another implementation of Truncated Multivariate Normal function.

Fix patch:
diff --git a/README.kr.rst b/README.kr.rst
index 88b8b33..8e18e26 100644
--- a/README.kr.rst
+++ b/README.kr.rst
@@ -245,7 +245,7 @@ LDA모델로 1000회 iteration을 수행시 걸리는 시간을 초 단위로 
 
 예제 코드
 ---------
-tomotopy의 Python3 예제 코드는 https://github.com/bab2min/tomotopy/blob/master/examples/ 를 확인하시길 바랍니다.
+tomotopy의 Python3 예제 코드는 https://github.com/bab2min/tomotopy/blob/main/examples/ 를 확인하시길 바랍니다.
 
 예제 코드에서 사용했던 데이터 파일은 https://drive.google.com/file/d/18OpNijd4iwPyYZ2O7pQoPyeTAKEXa71J/view 에서 다운받을 수 있습니다.
 
diff --git a/README.rst b/README.rst
index 326a02a..6f21764 100644
--- a/README.rst
+++ b/README.rst
@@ -250,7 +250,7 @@ See `word_prior_example` in `example.py` for more details.
 
 Examples
 --------
-You can find an example python code of tomotopy at https://github.com/bab2min/tomotopy/blob/master/examples/ .
+You can find an example python code of tomotopy at https://github.com/bab2min/tomotopy/blob/main/examples/ .
 
 You can also get the data file used in the example code at https://drive.google.com/file/d/18OpNijd4iwPyYZ2O7pQoPyeTAKEXa71J/view .
 
diff --git a/src/TopicModel/CTModel.hpp b/src/TopicModel/CTModel.hpp
index aff3a3b..d862290 100644
--- a/src/TopicModel/CTModel.hpp
+++ b/src/TopicModel/CTModel.hpp
@@ -65,6 +65,7 @@ namespace tomoto
 			{
 				if (i == 0) pbeta = Eigen::Matrix<Float, -1, 1>::Ones(this->K);
 				else pbeta = doc.beta.col(i % numBetaSample).array().exp();
+
 				Float betaESum = pbeta.sum() + 1;
 				pbeta /= betaESum;
 				for (size_t k = 0; k < this->K; ++k)
@@ -78,7 +79,9 @@ namespace tomoto
 
 					Float c = betaESum * (1 - pbeta[k]);
 					lowerBound[k] = log(c * max_uk / (1 - max_uk));
-					upperBound[k] = log(c * min_unk / (1 - min_unk));
+					lowerBound[k] = std::max(std::min(lowerBound[k], (Float)100), (Float)-100);
+					upperBound[k] = log(c * min_unk / (1 - min_unk + epsilon));
+					upperBound[k] = std::max(std::min(upperBound[k], (Float)100), (Float)-100);
 					if (lowerBound[k] > upperBound[k])
 					{
 						THROW_ERROR_WITH_INFO(exception::TrainingError,
diff --git a/src/Utils/MultiNormalDistribution.hpp b/src/Utils/MultiNormalDistribution.hpp
index 7ef247a..f2d55b0 100644
--- a/src/Utils/MultiNormalDistribution.hpp
+++ b/src/Utils/MultiNormalDistribution.hpp
@@ -51,8 +51,9 @@ namespace tomoto
 					}
 					if (len > 1) newDist.cov /= len - 1;
 				}
-				newDist.l = newDist.cov.llt().matrixL();
-				newDist.logDet = newDist.l.diagonal().array().log().sum();
+				Eigen::MatrixXd l = newDist.cov.template cast<double>().llt().matrixL();
+				newDist.l = l.template cast<float>();
+				newDist.logDet = l.diagonal().array().log().sum();
 				return newDist;
 			}
 
diff --git a/src/Utils/TruncMultiNormal.hpp b/src/Utils/TruncMultiNormal.hpp
index 641918b..e6eb307 100644
--- a/src/Utils/TruncMultiNormal.hpp
+++ b/src/Utils/TruncMultiNormal.hpp
@@ -14,57 +14,55 @@ namespace tomoto
 			const Eigen::Matrix<_Ty, -1, 1>& lowerBound,
 			const Eigen::Matrix<_Ty, -1, 1>& upperBound,
 			_Rng& rng,
-			size_t iteration)
+			size_t burnIn
+		)
 		{
-			constexpr _Ty epsilon = 1e-6;
 			const size_t K = ret.size();
-			Eigen::Matrix<_Ty, -1, 1> bias = Eigen::Matrix<_Ty, -1, 1>::Zero(K), lowers, uppers;
-			auto& l = multiNormal.getCovL();
-			ret.setZero();
-
-			std::vector<size_t> ks(K);
-			std::iota(ks.begin(), ks.end(), 0);
-			for (size_t i = 0; i < iteration; ++i)
+			Eigen::Matrix<_Ty, -1, -1> l = multiNormal.getCovL();
+			ret = (lowerBound + upperBound) / 2;
+			Eigen::Matrix<_Ty, -1, 1> z = l.template triangularView<Eigen::Lower>().solve(ret - multiNormal.mean),
+				a = lowerBound - multiNormal.mean,
+				b = upperBound - multiNormal.mean,
+				t, at, bt;
+			for (size_t i = 0; i < burnIn; ++i)
 			{
-				// shuffle sampling orders except during initialization
-				if (i) std::shuffle(ks.begin(), ks.end(), rng);
-				for (size_t kx = 0; kx < K; ++kx)
+				for (size_t j = 0; j < K; ++j)
 				{
-					size_t k = ks[kx];
-					ret[k] = 0;
-					//bias = multiNormal.mean + l * ret;
-					//bias.tail(K - k) = multiNormal.mean.tail(K - k) + l.block(k, 0, K - k, K) * ret;
-					bias.tail(K - k) = multiNormal.mean.tail(K - k);
-					bias.tail(K - k).noalias() += l.block(k, 0, K - k, K) * ret;
-					lowers = (lowerBound - bias).tail(K - k).array() / l.col(k).tail(K - k).array();
-					uppers = (upperBound - bias).tail(K - k).array() / l.col(k).tail(K - k).array();
-					_Ty nLower = lowers[0], nUpper = uppers[0];
-					if (l(k, k) < 0) std::swap(nLower, nUpper);
-					if (i)
+					auto lj = l.col(j);
+					z[j] = 0;
+					t = l * z;
+					_Ty lower_pos = -INFINITY, upper_pos = INFINITY,
+						lower_neg = -INFINITY, upper_neg = INFINITY;
+					at = ((a - t).array() / lj.array()).matrix();
+					bt = ((b - t).array() / lj.array()).matrix();
+					for (size_t k = 0; k < K; ++k)
 					{
-						for (size_t j = 1; j < lowers.size(); ++j)
+						if (lj[k] > 0)
+						{
+							lower_pos = std::max(lower_pos, at[k]);
+							upper_pos = std::min(upper_pos, bt[k]);
+						}
+						else if (lj[k] < 0)
 						{
-							if (l.col(k)(j + k) > epsilon)
-							{
-								if (lowers[j] > nLower) nLower = lowers[j];
-								if (uppers[j] < nUpper) nUpper = uppers[j];
-							}
-							else if (l.col(k)(j + k) < -epsilon)
-							{
-								if (uppers[j] > nLower) nLower = uppers[j];
-								if (lowers[j] < nUpper) nUpper = lowers[j];
-							}
+							lower_neg = std::max(lower_neg, bt[k]);
+							upper_neg = std::min(upper_neg, at[k]);
 						}
 					}
-					if (abs(nLower - nUpper) <= 1e-4) ret[k] = (nLower + nUpper) / 2;
+					lower_pos = std::max(lower_pos, lower_neg);
+					upper_pos = std::min(upper_pos, upper_neg);
+					// this is due to numerical instability
+					if (lower_pos >= upper_pos)
+					{
+						std::cerr << __FILE__ << "(" << __LINE__ << "): wrong truncation range [" << lower_pos << ", " << upper_pos << "]" << std::endl;
+						z[j] = (lower_pos + upper_pos) / 2;
+					}
 					else
 					{
-						ret[k] = rtnorm::rtnorm(rng, nLower, nUpper);
+						z[j] = rtnorm::rtnorm(rng, lower_pos, upper_pos);
 					}
 				}
 			}
-			ret = l * ret;
-			ret += multiNormal.mean;
+			ret = (l * z) + multiNormal.mean;
 			return ret;
 		}
 
diff --git a/src/python/py_utils.cpp b/src/python/py_utils.cpp
index f39a42f..50b464a 100644
--- a/src/python/py_utils.cpp
+++ b/src/python/py_utils.cpp
@@ -480,6 +480,8 @@ PyObject* CorpusObject::addDoc(CorpusObject* self, PyObject* args, PyObject* kwa
 					throw runtime_error{ string{ "there is a document with uid = " } + PyUnicode_AsUTF8(repr) + " already." };
 				}
 				self->invmap.emplace(suid, self->docs.size());
+				doc.docUid = tomoto::SharedString{ uid };
+				continue;
 			}
 
 			Py_INCREF(value);
@@ -1506,9 +1508,11 @@ vector<size_t> insertCorpus(TopicModelObject* self, PyObject* _corpus, PyObject*
 
 		for (size_t i = 0; i < rdoc.words.size(); ++i)
 		{
-			if(rdoc.words[i] == tomoto::non_vocab_id) continue;
+			if (rdoc.words[i] == tomoto::non_vocab_id) continue;
+			
 			if (insert_into_empty) doc.words.emplace_back(rdoc.words[i]);
 			else doc.words.emplace_back(corpus->getVocabDict().mapToNewDict(rdoc.words[i], self->inst->getVocabDict()));
+			
 			if (!doc.rawStr.empty())
 			{
 				doc.origWordPos.emplace_back(rdoc.origWordPos[i]);
diff --git a/tomotopy/documentation.kr.rst b/tomotopy/documentation.kr.rst
index 09c213e..aeb9a9b 100644
--- a/tomotopy/documentation.kr.rst
+++ b/tomotopy/documentation.kr.rst
@@ -326,7 +326,7 @@ LDA모델로 1000회 iteration을 수행시 걸리는 시간을 초 단위로 
 
 예제 코드
 ---------
-tomotopy의 Python3 예제 코드는 https://github.com/bab2min/tomotopy/blob/master/examples/ 를 확인하시길 바랍니다.
+tomotopy의 Python3 예제 코드는 https://github.com/bab2min/tomotopy/blob/main/examples/ 를 확인하시길 바랍니다.
 
 예제 코드에서 사용했던 데이터 파일은 https://drive.google.com/file/d/18OpNijd4iwPyYZ2O7pQoPyeTAKEXa71J/view 에서 다운받을 수 있습니다.
 
diff --git a/tomotopy/documentation.rst b/tomotopy/documentation.rst
index a27d74d..2163ed9 100644
--- a/tomotopy/documentation.rst
+++ b/tomotopy/documentation.rst
@@ -329,7 +329,7 @@ See `word_prior_example` in `example.py` for more details.
 
 Examples
 --------
-You can find an example python code of tomotopy at https://github.com/bab2min/tomotopy/blob/master/examples/ .
+You can find an example python code of tomotopy at https://github.com/bab2min/tomotopy/blob/main/examples/ .
 
 You can also get the data file used in the example code at https://drive.google.com/file/d/18OpNijd4iwPyYZ2O7pQoPyeTAKEXa71J/view .
 


Test patch:
diff --git a/.github/workflows/pull_request_test.yml b/.github/workflows/pull_request_test.yml
index 0376c9b..b7c5f2d 100644
--- a/.github/workflows/pull_request_test.yml
+++ b/.github/workflows/pull_request_test.yml
@@ -14,7 +14,7 @@ jobs:
     strategy:
       max-parallel: 4
       matrix:
-        cp: [cp35-cp35m, cp36-cp36m, cp37-cp37m, cp38-cp38, cp39-cp39]
+        cp: [cp36-cp36m, cp37-cp37m, cp38-cp38, cp39-cp39]
 
     steps:        
     - uses: actions/checkout@v1
@@ -54,7 +54,7 @@ jobs:
     strategy:
       max-parallel: 4
       matrix:
-        python-version: [3.5, 3.6, 3.7, 3.8, 3.9]
+        python-version: [3.6, 3.7, 3.8, 3.9]
 
     steps:
     - uses: actions/checkout@v1
@@ -90,7 +90,7 @@ jobs:
     strategy:
       max-parallel: 4
       matrix:
-        python-version: [3.5, 3.6, 3.7, 3.8, 3.9]
+        python-version: [3.6, 3.7, 3.8, 3.9]
         architecture: [x86, x64]
 
     steps:
@@ -117,6 +117,13 @@ jobs:
       run: |
         python -m pip install numpy
         python setup.py build install
+    - name: Archive binary
+      uses: actions/upload-artifact@v2
+      with:
+        name: pyd ${{ matrix.python-version }} ${{ matrix.architecture }}
+        path: |
+          build/**/*.pyd
+          build/**/*.pdb
     - name: Test
       run: |
         python -m pip install pytest nltk