Repository: bab2min/tomotopy (language: python)

Issue:
`Corpus.extract_ngrams` crashes when empty doc is inserted
`Corpus.extract_ngrams` crashes when empty doc is inserted. 
Also, it would be nice if there is an option to use NPMI instead of PMI at `extract_ngrams` and `label.Candidate` provides candidates' count information (frequency).


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/generate_documentation.yml b/.github/workflows/generate_documentation.yml
index ec78e95..e62a819 100644
--- a/.github/workflows/generate_documentation.yml
+++ b/.github/workflows/generate_documentation.yml
@@ -39,7 +39,7 @@ jobs:
         mv variant/include/mapbox include/
     - name: build
       run: |
-        python3 -m pip install pdoc3 numpy
+        python3 -m pip install pdoc3==0.8.4
         export TOMOTOPY_LANG=${{ matrix.language }}
         python3 setup.py install
     - name: gen doc
diff --git a/.gitignore b/.gitignore
index 69d5d5f..e40a588 100644
--- a/.gitignore
+++ b/.gitignore
@@ -13,5 +13,7 @@ enwiki-1000.txt
 /html/
 .vscode/
 *.pyd
+*.log
 examples/*
 !examples/*.py
+enwiki-16000.txt
diff --git a/examples/hlda_basic.py b/examples/hlda_basic.py
new file mode 100644
index 0000000..bdf1a4e
--- /dev/null
+++ b/examples/hlda_basic.py
@@ -0,0 +1,54 @@
+import sys
+import tomotopy as tp
+import numpy as np
+
+def hlda_example(input_file, save_path):
+    from nltk.stem.porter import PorterStemmer
+    from nltk.corpus import stopwords
+    try:
+        cps = tp.utils.Corpus.load(input_file + '.cached.cps')
+    except IOError:
+        stemmer = PorterStemmer()
+        stops = set(stopwords.words('english'))
+        cps = tp.utils.Corpus(
+            tokenizer=tp.utils.SimpleTokenizer(stemmer=stemmer.stem), 
+            stopwords=lambda x: len(x) <= 2 or x in stops
+        )
+        cps.process(open(input_file, encoding='utf-8'))
+        cps.save(input_file + '.cached.cps')
+    
+    np.random.seed(42)
+    ridcs = np.random.permutation(len(cps))
+    test_idcs = ridcs[:20]
+    train_idcs = ridcs[20:]
+
+    test_cps = cps[test_idcs]
+    train_cps = cps[train_idcs]
+    
+    mdl = tp.HLDAModel(tw=tp.TermWeight.ONE, min_df=10, depth=4, rm_top=10, corpus=train_cps)
+    mdl.train(0)
+    print('Num docs:', len(mdl.docs), ', Vocab size:', len(mdl.used_vocabs), ', Num words:', mdl.num_words)
+    print('Removed top words:', mdl.removed_top_words)
+    print('Training...', file=sys.stderr, flush=True)
+    for _ in range(0, 1000, 10):
+        mdl.train(7)
+        mdl.train(3, freeze_topics=True)
+        print('Iteration: {:05}\tll per word: {:.5f}\tNum. of topics: {}'.format(mdl.global_step, mdl.ll_per_word, mdl.live_k))
+
+    for _ in range(0, 100, 10):
+        mdl.train(10, freeze_topics=True)
+        print('Iteration: {:05}\tll per word: {:.5f}\tNum. of topics: {}'.format(mdl.global_step, mdl.ll_per_word, mdl.live_k))
+
+    mdl.summary()
+    print('Saving...', file=sys.stderr, flush=True)
+    mdl.save(save_path, True)
+
+    test_result_cps, ll = mdl.infer(test_cps)
+    for doc in test_result_cps:
+        print(doc.path, doc.get_words(top_n=10))
+
+# You can get the sample data file 'enwiki-16000.txt'
+# at https://drive.google.com/file/d/1OfyJ9TqaMiqzO6Qw-c_jXL-pmSIZf5Xt/view?usp=sharing
+
+if __name__ == '__main__':
+    hlda_example('enwiki-16000.txt', 'test.hlda.tmm')
diff --git a/src/Labeling/FoRelevance.cpp b/src/Labeling/FoRelevance.cpp
index f941b6a..374c9d1 100644
--- a/src/Labeling/FoRelevance.cpp
+++ b/src/Labeling/FoRelevance.cpp
@@ -2,6 +2,7 @@
 #include <numeric>
 
 #include "FoRelevance.h"
+#include "Phraser.hpp"
 
 using namespace tomoto::label;
 
@@ -23,6 +24,26 @@ class DocWrapper
 	{
 		return doc->words[doc->wOrder.empty() ? idx : doc->wOrder[idx]];
 	}
+
+	auto begin() const -> decltype(doc->words.begin())
+	{
+		return doc->words.begin();
+	}
+
+	auto end() const -> decltype(doc->words.end())
+	{
+		return doc->words.end();
+	}
+
+	auto rbegin() const -> decltype(doc->words.rbegin())
+	{
+		return doc->words.rbegin();
+	}
+
+	auto rend() const -> decltype(doc->words.rend())
+	{
+		return doc->words.rend();
+	}
 };
 
 class DocIterator
@@ -61,9 +82,10 @@ std::vector<Candidate> PMIExtractor::extract(const tomoto::ITopicModel* tm) cons
 {
 	auto& vocabFreqs = tm->getVocabCf();
 	auto& vocabDf = tm->getVocabDf();
-	auto candidates = extractPMINgrams(DocIterator{ tm, 0 }, DocIterator{ tm, tm->getNumDocs() }, 
+	auto candidates = phraser::extractPMINgrams(DocIterator{ tm, 0 }, DocIterator{ tm, tm->getNumDocs() }, 
 		vocabFreqs, vocabDf,
-		candMinCnt, candMinDf, minLabelLen, maxLabelLen, maxCandidates, -99999.f
+		candMinCnt, candMinDf, minLabelLen, maxLabelLen, maxCandidates, 0.f,
+		normalized
 	);
 	if (minLabelLen <= 1)
 	{
@@ -77,6 +99,29 @@ std::vector<Candidate> PMIExtractor::extract(const tomoto::ITopicModel* tm) cons
 	return candidates;
 }
 
+
+std::vector<Candidate> tomoto::label::PMIBEExtractor::extract(const ITopicModel* tm) const
+{
+	auto& vocabFreqs = tm->getVocabCf();
+	auto& vocabDf = tm->getVocabDf();
+	auto candidates = phraser::extractPMIBENgrams(DocIterator{ tm, 0 }, DocIterator{ tm, tm->getNumDocs() },
+		vocabFreqs, vocabDf,
+		candMinCnt, candMinDf, minLabelLen, maxLabelLen, maxCandidates, 
+		0.f, 0.f
+	);
+	if (minLabelLen <= 1)
+	{
+		for (size_t i = 0; i < vocabDf.size(); ++i)
+		{
+			if (vocabFreqs[i] < candMinCnt) continue;
+			if (vocabDf[i] < candMinDf) continue;
+			candidates.emplace_back(0.f, i);
+		}
+	}
+	return candidates;
+}
+
+
 template<bool _lock>
 const Eigen::ArrayXi& FoRelevance::updateContext(size_t docId, const tomoto::DocumentBase* doc, const tomoto::Trie<tomoto::Vid, size_t>* root)
 {
diff --git a/src/Labeling/FoRelevance.h b/src/Labeling/FoRelevance.h
index 60c1f4b..481de5f 100644
--- a/src/Labeling/FoRelevance.h
+++ b/src/Labeling/FoRelevance.h
@@ -4,6 +4,7 @@
 #include "Labeler.h"
 #include "../Utils/EigenAddonOps.hpp"
 #include "../Utils/Trie.hpp"
+#include "../Utils/ThreadPool.hpp"
 
 /*
 Implementation of First-order Relevance for topic labeling by bab2min
@@ -16,166 +17,35 @@ namespace tomoto
 {
 	namespace label
 	{
-		template<typename _DocIter, typename _Freqs>
-		std::vector<Candidate> extractPMINgrams(_DocIter docBegin, _DocIter docEnd, 
-			_Freqs&& vocabFreqs, _Freqs&& vocabDf,
-			size_t candMinCnt, size_t candMinDf, size_t minNgrams, size_t maxNgrams, size_t maxCandidates, float minScore)
+		class PMIExtractor : public IExtractor
 		{
-			struct vvhash
-			{
-				size_t operator()(const std::pair<Vid, Vid>& k) const
-				{
-					return std::hash<Vid>{}(k.first) ^ std::hash<Vid>{}(k.second);
-				}
-			};
-
-			// counting unigrams & bigrams
-			std::unordered_map<std::pair<Vid, Vid>, size_t, vvhash> bigramCnt, bigramDf;
-
-			for(auto docIt = docBegin; docIt != docEnd; ++docIt)
-			{
-				std::unordered_set<std::pair<Vid, Vid>, vvhash> uniqBigram;
-				auto doc = *docIt;
-				Vid prevWord = doc[0];
-				for (size_t j = 1; j < doc.size(); ++j)
-				{
-					Vid curWord = doc[j];
-					if (curWord != non_vocab_id && vocabFreqs[curWord] >= candMinCnt && vocabDf[curWord] >= candMinDf)
-					{
-						if (prevWord != non_vocab_id && vocabFreqs[prevWord] >= candMinCnt && vocabDf[prevWord] >= candMinDf)
-						{
-							bigramCnt[std::make_pair(prevWord, curWord)]++;
-							uniqBigram.emplace(prevWord, curWord);
-						}
-					}
-					prevWord = curWord;
-				}
-
-				for (auto& p : uniqBigram) bigramDf[p]++;
-			}
-
-
-			// counting ngrams
-			std::vector<TrieEx<Vid, size_t>> trieNodes;
-
-			if (maxNgrams > 2)
-			{
-				std::unordered_set<std::pair<Vid, Vid>, vvhash> validPair;
-				for (auto& p : bigramCnt)
-				{
-					if (p.second >= candMinCnt) validPair.emplace(p.first);
-				}
-
-				trieNodes.resize(1);
-				auto allocNode = [&]() { return trieNodes.emplace_back(), & trieNodes.back(); };
-
-				for (auto docIt = docBegin; docIt != docEnd; ++docIt)
-				{
-					auto doc = *docIt;
-					if (trieNodes.capacity() < trieNodes.size() + doc.size() * maxNgrams)
-					{
-						trieNodes.reserve(std::max(trieNodes.size() + doc.size() * maxNgrams, trieNodes.capacity() * 2));
-					}
-
-					Vid prevWord = doc[0];
-					size_t labelLen = 0;
-					auto node = &trieNodes[0];
-					if (prevWord != non_vocab_id && vocabFreqs[prevWord] >= candMinCnt)
-					{
-						node = trieNodes[0].makeNext(prevWord, allocNode);
-						node->val++;
-						labelLen = 1;
-					}
-
-					for (size_t j = 1; j < doc.size(); ++j)
-					{
-						Vid curWord = doc[j];
-
-						if (curWord != non_vocab_id && vocabFreqs[curWord] < candMinCnt)
-						{
-							node = &trieNodes[0];
-							labelLen = 0;
-						}
-						else
-						{
-							if (labelLen >= maxNgrams)
-							{
-								node = node->getFail();
-								labelLen--;
-							}
-
-							if (validPair.count(std::make_pair(prevWord, curWord)))
-							{
-								auto nnode = node->makeNext(curWord, allocNode);
-								node = nnode;
-								do
-								{
-									nnode->val++;
-								} while (nnode = nnode->getFail());
-								labelLen++;
-							}
-							else
-							{
-								node = trieNodes[0].makeNext(curWord, allocNode);
-								node->val++;
-								labelLen = 1;
-							}
-						}
-						prevWord = curWord;
-					}
-				}
-			}
-
-			float totN = std::accumulate(vocabFreqs.begin(), vocabFreqs.end(), (size_t)0);
-
-			// calculating PMIs
-			std::vector<Candidate> candidates;
-			for (auto& p : bigramCnt)
-			{
-				auto& bigram = p.first;
-				if (p.second < candMinCnt) continue;
-				if (bigramDf[bigram] < candMinDf) continue;
-				auto pmi = std::log(p.second * totN
-					/ vocabFreqs[bigram.first] / vocabFreqs[bigram.second]);
-				if (pmi <= 0) continue;
-				candidates.emplace_back(pmi, bigram.first, bigram.second);
-			}
-
-			if (maxNgrams > 2)
+			size_t candMinCnt, candMinDf, minLabelLen, maxLabelLen, maxCandidates;
+			bool normalized;
+		public:
+			PMIExtractor(size_t _candMinCnt = 10, size_t _candMinDf = 2, 
+				size_t _minLabelLen = 1, size_t _maxLabelLen = 5, size_t _maxCandidates = 1000,
+				bool _normalized = false
+			)
+				: candMinCnt{ _candMinCnt }, candMinDf{ _candMinDf }, 
+				minLabelLen{ _minLabelLen }, maxLabelLen{ _maxLabelLen }, 
+				maxCandidates{ _maxCandidates }, normalized{ _normalized }
 			{
-				std::vector<Vid> rkeys;
-				trieNodes[0].traverse_with_keys([&](const TrieEx<Vid, size_t>* node, const std::vector<Vid>& rkeys)
-				{
-					if (rkeys.size() <= 2 || rkeys.size() < minNgrams || node->val < candMinCnt) return;
-					auto pmi = node->val / totN;
-					for (auto k : rkeys)
-					{
-						pmi *= totN / vocabFreqs[k];
-					}
-					pmi = std::log(pmi);
-					if (pmi < minScore) return;
-					candidates.emplace_back(pmi, rkeys);
-				}, rkeys);
 			}
 
-			std::sort(candidates.begin(), candidates.end(), [](const Candidate& a, const Candidate& b)
-			{
-				return a.score > b.score;
-			});
-			if (candidates.size() > maxCandidates) candidates.erase(candidates.begin() + maxCandidates, candidates.end());
-			return candidates;
-		}
-
+			std::vector<Candidate> extract(const ITopicModel* tm) const override;
+		};
 
-		class PMIExtractor : public IExtractor
+		class PMIBEExtractor : public IExtractor
 		{
 			size_t candMinCnt, candMinDf, minLabelLen, maxLabelLen, maxCandidates;
 		public:
-			PMIExtractor(size_t _candMinCnt = 10, size_t _candMinDf = 2, size_t _minLabelLen = 1, size_t _maxLabelLen = 5, size_t _maxCandidates = 1000)
-				: candMinCnt{ _candMinCnt }, candMinDf{ _candMinDf }, minLabelLen{ _minLabelLen}, maxLabelLen{ _maxLabelLen }, maxCandidates{ _maxCandidates }
+			PMIBEExtractor(size_t _candMinCnt = 10, size_t _candMinDf = 2, 
+				size_t _minLabelLen = 1, size_t _maxLabelLen = 5, size_t _maxCandidates = 1000
+			)
+				: candMinCnt{ _candMinCnt }, candMinDf{ _candMinDf }, minLabelLen{ _minLabelLen }, maxLabelLen{ _maxLabelLen }, maxCandidates{ _maxCandidates }
 			{
 			}
-			
+
 			std::vector<Candidate> extract(const ITopicModel* tm) const override;
 		};
 
@@ -212,7 +82,7 @@ namespace tomoto
 
 		public:
 			template<typename _Iter>
-			FoRelevance(const ITopicModel* _tm, 
+			FoRelevance(const ITopicModel* _tm,
 				_Iter candFirst, _Iter candEnd,
 				size_t _candMinDf = 2, float _smoothing = 0.1f, float _lambda = 0.1f, float _mu = 0.1f,
 				size_t _windowSize = (size_t)-1,
diff --git a/src/Labeling/Labeler.h b/src/Labeling/Labeler.h
index ab1ff36..6c2a1b4 100644
--- a/src/Labeling/Labeler.h
+++ b/src/Labeling/Labeler.h
@@ -10,6 +10,7 @@ namespace tomoto
 		struct Candidate
 		{
 			float score = 0;
+			size_t cf = 0, df = 0;
 			std::vector<Vid> w;
 			std::string name;
 
@@ -18,17 +19,17 @@ namespace tomoto
 			}
 
 			Candidate(float _score, Vid w1)
-				: w{ w1 }, score{ _score }
+				: score{ _score }, w{ w1 }
 			{
 			}
 
 			Candidate(float _score, Vid w1, Vid w2)
-				: w{ w1, w2 }, score{ _score }
+				: score{ _score }, w{ w1, w2 }
 			{
 			}
 
 			Candidate(float _score, const std::vector<Vid>& _w)
-				: w{ _w }, score{ _score }
+				: score{ _score }, w{ _w }
 			{
 			}
 		};
@@ -36,6 +37,7 @@ namespace tomoto
 		class IExtractor
 		{
 		public:
+			
 			virtual std::vector<Candidate> extract(const ITopicModel* tm) const = 0;
 			virtual ~IExtractor() {}
 		};
diff --git a/src/Labeling/Phraser.hpp b/src/Labeling/Phraser.hpp
new file mode 100644
index 0000000..e1fe91d
--- /dev/null
+++ b/src/Labeling/Phraser.hpp
@@ -0,0 +1,518 @@
+#pragma once
+
+#include <vector>
+#include <unordered_map>
+#include "Labeler.h"
+#include "../Utils/Trie.hpp"
+
+namespace tomoto
+{
+	namespace phraser
+	{
+		template<typename _DocIter>
+		void countUnigrams(std::vector<size_t>& unigramCf, std::vector<size_t>& unigramDf,
+			_DocIter docBegin, _DocIter docEnd
+		)
+		{
+			for (auto docIt = docBegin; docIt != docEnd; ++docIt)
+			{
+				auto doc = *docIt;
+				if (!doc.size()) continue;
+				std::unordered_set<Vid> uniqs;
+				for (size_t i = 0; i < doc.size(); ++i)
+				{
+					if (doc[i] == non_vocab_id) continue;
+					unigramCf[doc[i]]++;
+					uniqs.emplace(doc[i]);
+				}
+
+				for (auto w : uniqs) unigramDf[w]++;
+			}
+		}
+
+		template<typename _DocIter, typename _VvHash, typename _Freqs>
+		void countBigrams(std::unordered_map<std::pair<Vid, Vid>, size_t, _VvHash>& bigramCf,
+			std::unordered_map<std::pair<Vid, Vid>, size_t, _VvHash>& bigramDf,
+			_DocIter docBegin, _DocIter docEnd,
+			_Freqs&& vocabFreqs, _Freqs&& vocabDf,
+			size_t candMinCnt, size_t candMinDf
+		)
+		{
+			for (auto docIt = docBegin; docIt != docEnd; ++docIt)
+			{
+				std::unordered_set<std::pair<Vid, Vid>, _VvHash> uniqBigram;
+				auto doc = *docIt;
+				if (!doc.size()) continue;
+				Vid prevWord = doc[0];
+				for (size_t j = 1; j < doc.size(); ++j)
+				{
+					Vid curWord = doc[j];
+					if (curWord != non_vocab_id && vocabFreqs[curWord] >= candMinCnt && vocabDf[curWord] >= candMinDf)
+					{
+						if (prevWord != non_vocab_id && vocabFreqs[prevWord] >= candMinCnt && vocabDf[prevWord] >= candMinDf)
+						{
+							bigramCf[std::make_pair(prevWord, curWord)]++;
+							uniqBigram.emplace(prevWord, curWord);
+						}
+					}
+					prevWord = curWord;
+				}
+
+				for (auto& p : uniqBigram) bigramDf[p]++;
+			}
+		}
+
+		template<bool _reverse, typename _DocIter, typename _Freqs, typename _BigramPairs>
+		void countNgrams(std::vector<TrieEx<Vid, size_t>>& dest,
+			_DocIter docBegin, _DocIter docEnd,
+			_Freqs&& vocabFreqs, _Freqs&& vocabDf, _BigramPairs&& validPairs,
+			size_t candMinCnt, size_t candMinDf, size_t maxNgrams
+		)
+		{
+			if (dest.empty())
+			{
+				dest.resize(1);
+				dest.reserve(1024);
+			}
+			auto allocNode = [&]() { return dest.emplace_back(), & dest.back(); };
+
+			for (auto docIt = docBegin; docIt != docEnd; ++docIt)
+			{
+				auto doc = *docIt;
+				if (!doc.size()) continue;
+				if (dest.capacity() < dest.size() + doc.size() * maxNgrams)
+				{
+					dest.reserve(std::max(dest.size() + doc.size() * maxNgrams, dest.capacity() * 2));
+				}
+
+				Vid prevWord = _reverse ? *doc.rbegin() : *doc.begin();
+				size_t labelLen = 0;
+				auto node = &dest[0];
+				if (prevWord != non_vocab_id && vocabFreqs[prevWord] >= candMinCnt && vocabDf[prevWord] >= candMinDf)
+				{
+					node = dest[0].makeNext(prevWord, allocNode);
+					node->val++;
+					labelLen = 1;
+				}
+
+				const auto func = [&](Vid curWord)
+				{
+					if (curWord != non_vocab_id && (vocabFreqs[curWord] < candMinCnt || vocabDf[curWord] < candMinDf))
+					{
+						node = &dest[0];
+						labelLen = 0;
+					}
+					else
+					{
+						if (labelLen >= maxNgrams)
+						{
+							node = node->getFail();
+							labelLen--;
+						}
+
+						if (validPairs.count(_reverse ? std::make_pair(curWord, prevWord) : std::make_pair(prevWord, curWord)))
+						{
+							auto nnode = node->makeNext(curWord, allocNode);
+							node = nnode;
+							do
+							{
+								nnode->val++;
+							} while ((nnode = nnode->getFail()));
+							labelLen++;
+						}
+						else
+						{
+							node = dest[0].makeNext(curWord, allocNode);
+							node->val++;
+							labelLen = 1;
+						}
+					}
+					prevWord = curWord;
+				};
+
+				if (_reverse) std::for_each(doc.rbegin() + 1, doc.rend(), func);
+				else std::for_each(doc.begin() + 1, doc.end(), func);
+			}
+		}
+
+		inline void mergeNgramCounts(std::vector<TrieEx<Vid, size_t>>& dest, std::vector<TrieEx<Vid, size_t>>&& src)
+		{
+			if (src.empty()) return;
+			if (dest.empty()) dest.resize(1);
+
+			auto allocNode = [&]() { return dest.emplace_back(), & dest.back(); };
+
+			std::vector<Vid> rkeys;
+			src[0].traverse_with_keys([&](const TrieEx<Vid, size_t>* node, const std::vector<Vid>& rkeys)
+			{
+				if (dest.capacity() < dest.size() + rkeys.size() * rkeys.size())
+				{
+					dest.reserve(std::max(dest.size() + rkeys.size() * rkeys.size(), dest.capacity() * 2));
+				}
+				dest[0].build(rkeys.begin(), rkeys.end(), 0, allocNode)->val += node->val;
+			}, rkeys);
+		}
+
+		inline float branchingEntropy(const TrieEx<Vid, size_t>* node, size_t minCnt)
+		{
+			float entropy = 0;
+			size_t rest = node->val;
+			for (auto n : *node)
+			{
+				float p = n.second->val / (float)node->val;
+				entropy -= p * std::log(p);
+				rest -= n.second->val;
+			}
+			if (rest > 0)
+			{
+				float p = rest / (float)node->val;
+				entropy -= p * std::log(std::min(std::max(minCnt, (size_t)1), (size_t)rest) / (float)node->val);
+			}
+			return entropy;
+		}
+
+		template<typename _LocalData, typename _ReduceFn>
+		_LocalData parallelReduce(std::vector<_LocalData>&& data, _ReduceFn&& fn, ThreadPool* pool = nullptr)
+		{
+			if (pool)
+			{
+				for (size_t s = data.size(); s > 1; s = (s + 1) / 2)
+				{
+					std::vector<std::future<void>> futures;
+					size_t h = (s + 1) / 2;
+					for (size_t i = h; i < s; ++i)
+					{
+						futures.emplace_back(pool->enqueue([&, i, h](size_t)
+						{
+							_LocalData d = std::move(data[i]);
+							fn(data[i - h], std::move(d));
+						}));
+					}
+					for (auto& f : futures) f.get();
+				}
+			}
+			else
+			{
+				for (size_t i = 1; i < data.size(); ++i)
+				{
+					_LocalData d = std::move(data[i]);
+					fn(data[0], std::move(d));
+				}
+			}
+			return std::move(data[0]);
+		}
+
+		namespace detail
+		{
+			struct vvhash
+			{
+				size_t operator()(const std::pair<Vid, Vid>& k) const
+				{
+					return std::hash<Vid>{}(k.first) ^ std::hash<Vid>{}(k.second);
+				}
+			};
+		}
+
+		template<typename _DocIter, typename _Freqs>
+		std::vector<label::Candidate> extractPMINgrams(_DocIter docBegin, _DocIter docEnd,
+			_Freqs&& vocabFreqs, _Freqs&& vocabDf,
+			size_t candMinCnt, size_t candMinDf, size_t minNgrams, size_t maxNgrams, size_t maxCandidates,
+			float minScore, bool normalized = false,
+			ThreadPool* pool = nullptr)
+		{
+			// counting unigrams & bigrams
+			std::unordered_map<std::pair<Vid, Vid>, size_t, detail::vvhash> bigramCnt, bigramDf;
+
+			if (pool && pool->getNumWorkers() > 1)
+			{
+				using LocalCfDf = std::pair<
+					std::unordered_map<std::pair<Vid, Vid>, size_t, detail::vvhash>,
+					std::unordered_map<std::pair<Vid, Vid>, size_t, detail::vvhash>
+				>;
+				std::vector<LocalCfDf> localdata(pool->getNumWorkers());
+				std::vector<std::future<void>> futures;
+				const size_t stride = pool->getNumWorkers() * 8;
+				auto docIt = docBegin;
+				for (size_t i = 0; i < stride && docIt != docEnd; ++i, ++docIt)
+				{
+					futures.emplace_back(pool->enqueue([&, docIt, stride](size_t tid)
+					{
+						countBigrams(localdata[tid].first, localdata[tid].second, makeStrideIter(docIt, stride, docEnd), makeStrideIter(docEnd, stride, docEnd), vocabFreqs, vocabDf, candMinCnt, candMinDf);
+					}));
+				}
+
+				for (auto& f : futures) f.get();
+
+				auto r = parallelReduce(std::move(localdata), [](LocalCfDf& dest, LocalCfDf&& src)
+				{
+					for (auto& p : src.first) dest.first[p.first] += p.second;
+					for (auto& p : src.second) dest.second[p.first] += p.second;
+				}, pool);
+
+				bigramCnt = std::move(r.first);
+				bigramDf = std::move(r.second);
+			}
+			else
+			{
+				countBigrams(bigramCnt, bigramDf, docBegin, docEnd, vocabFreqs, vocabDf, candMinCnt, candMinDf);
+			}
+
+			// counting ngrams
+			std::vector<TrieEx<Vid, size_t>> trieNodes;
+			if (maxNgrams > 2)
+			{
+				std::unordered_set<std::pair<Vid, Vid>, detail::vvhash> validPairs;
+				for (auto& p : bigramCnt)
+				{
+					if (p.second >= candMinCnt && bigramDf[p.first] >= candMinDf) validPairs.emplace(p.first);
+				}
+
+				if (pool && pool->getNumWorkers() > 1)
+				{
+					using LocalFwBw = std::vector<TrieEx<Vid, size_t>>;
+					std::vector<LocalFwBw> localdata(pool->getNumWorkers());
+					std::vector<std::future<void>> futures;
+					const size_t stride = pool->getNumWorkers() * 8;
+					auto docIt = docBegin;
+					for (size_t i = 0; i < stride && docIt != docEnd; ++i, ++docIt)
+					{
+						futures.emplace_back(pool->enqueue([&, docIt, stride](size_t tid)
+						{
+							countNgrams<false>(localdata[tid],
+								makeStrideIter(docIt, stride, docEnd),
+								makeStrideIter(docEnd, stride, docEnd),
+								vocabFreqs, vocabDf, validPairs, candMinCnt, candMinDf, maxNgrams
+							);
+						}));
+					}
+
+					for (auto& f : futures) f.get();
+
+					auto r = parallelReduce(std::move(localdata), [&](LocalFwBw& dest, LocalFwBw&& src)
+					{
+						mergeNgramCounts(dest, std::move(src));
+					}, pool);
+
+					trieNodes = std::move(r);
+				}
+				else
+				{
+					countNgrams<false>(trieNodes,
+						docBegin, docEnd,
+						vocabFreqs, vocabDf, validPairs, candMinCnt, candMinDf, maxNgrams
+					);
+				}
+			}
+
+			float totN = std::accumulate(vocabFreqs.begin(), vocabFreqs.end(), (size_t)0);
+			const float logTotN = std::log(totN);
+
+			// calculating PMIs
+			std::vector<label::Candidate> candidates;
+			for (auto& p : bigramCnt)
+			{
+				auto& bigram = p.first;
+				if (p.second < candMinCnt) continue;
+				if (bigramDf[bigram] < candMinDf) continue;
+				auto pmi = std::log(p.second * totN
+					/ vocabFreqs[bigram.first] / vocabFreqs[bigram.second]);
+				if (normalized)
+				{
+					pmi /= std::log(totN / p.second);
+				}
+				if (pmi < minScore) continue;
+				candidates.emplace_back(pmi, bigram.first, bigram.second);
+				candidates.back().cf = p.second;
+				candidates.back().df = bigramDf[bigram];
+			}
+
+			if (maxNgrams > 2)
+			{
+				std::vector<Vid> rkeys;
+				trieNodes[0].traverse_with_keys([&](const TrieEx<Vid, size_t>* node, const std::vector<Vid>& rkeys)
+				{
+					if (rkeys.size() <= 2 || rkeys.size() < minNgrams || rkeys.size() > maxNgrams || node->val < candMinCnt) return;
+					auto pmi = std::log((float)node->val) - logTotN;
+					for (auto k : rkeys)
+					{
+						pmi += logTotN - std::log((float)vocabFreqs[k]);
+					}
+					if (normalized)
+					{
+						pmi /= (logTotN - std::log((float)node->val)) * (rkeys.size() - 1);
+					}
+					if (pmi < minScore) return;
+					candidates.emplace_back(pmi, rkeys);
+					candidates.back().cf = node->val;
+				}, rkeys);
+			}
+
+			std::sort(candidates.begin(), candidates.end(), [](const label::Candidate& a, const label::Candidate& b)
+			{
+				return a.score > b.score;
+			});
+			if (candidates.size() > maxCandidates) candidates.erase(candidates.begin() + maxCandidates, candidates.end());
+			return candidates;
+		}
+
+		template<typename _DocIter, typename _Freqs>
+		std::vector<label::Candidate> extractPMIBENgrams(_DocIter docBegin, _DocIter docEnd,
+			_Freqs&& vocabFreqs, _Freqs&& vocabDf,
+			size_t candMinCnt, size_t candMinDf, size_t minNgrams, size_t maxNgrams, size_t maxCandidates,
+			float minNPMI = 0, float minNBE = 0,
+			ThreadPool* pool = nullptr)
+		{
+			// counting unigrams & bigrams
+			std::unordered_map<std::pair<Vid, Vid>, size_t, detail::vvhash> bigramCnt, bigramDf;
+
+			if (pool && pool->getNumWorkers() > 1)
+			{
+				using LocalCfDf = std::pair<
+					std::unordered_map<std::pair<Vid, Vid>, size_t, detail::vvhash>,
+					std::unordered_map<std::pair<Vid, Vid>, size_t, detail::vvhash>
+				>;
+				std::vector<LocalCfDf> localdata(pool->getNumWorkers());
+				std::vector<std::future<void>> futures;
+				const size_t stride = pool->getNumWorkers() * 8;
+				auto docIt = docBegin;
+				for (size_t i = 0; i < stride && docIt != docEnd; ++i, ++docIt)
+				{
+					futures.emplace_back(pool->enqueue([&, docIt, stride](size_t tid)
+					{
+						countBigrams(localdata[tid].first, localdata[tid].second, makeStrideIter(docIt, stride, docEnd), makeStrideIter(docEnd, stride, docEnd), vocabFreqs, vocabDf, candMinCnt, candMinDf);
+					}));
+				}
+
+				for (auto& f : futures) f.get();
+
+				auto r = parallelReduce(std::move(localdata), [](LocalCfDf& dest, LocalCfDf&& src)
+				{
+					for (auto& p : src.first) dest.first[p.first] += p.second;
+					for (auto& p : src.second) dest.second[p.first] += p.second;
+				}, pool);
+
+				bigramCnt = std::move(r.first);
+				bigramDf = std::move(r.second);
+			}
+			else
+			{
+				countBigrams(bigramCnt, bigramDf, docBegin, docEnd, vocabFreqs, vocabDf, candMinCnt, candMinDf);
+			}
+
+			// counting ngrams
+			std::vector<TrieEx<Vid, size_t>> trieNodes, trieNodesBw;
+			if (maxNgrams > 2)
+			{
+				std::unordered_set<std::pair<Vid, Vid>, detail::vvhash> validPairs;
+				for (auto& p : bigramCnt)
+				{
+					if (p.second >= candMinCnt && bigramDf[p.first] >= candMinDf) validPairs.emplace(p.first);
+				}
+
+				if (pool && pool->getNumWorkers() > 1)
+				{
+					using LocalFwBw = std::pair<
+						std::vector<TrieEx<Vid, size_t>>,
+						std::vector<TrieEx<Vid, size_t>>
+					>;
+					std::vector<LocalFwBw> localdata(pool->getNumWorkers());
+					std::vector<std::future<void>> futures;
+					const size_t stride = pool->getNumWorkers() * 8;
+					auto docIt = docBegin;
+					for (size_t i = 0; i < stride && docIt != docEnd; ++i, ++docIt)
+					{
+						futures.emplace_back(pool->enqueue([&, docIt, stride](size_t tid)
+						{
+							countNgrams<false>(localdata[tid].first, 
+								makeStrideIter(docIt, stride, docEnd), 
+								makeStrideIter(docEnd, stride, docEnd), 
+								vocabFreqs, vocabDf, validPairs, candMinCnt, candMinDf, maxNgrams + 1
+							);
+							countNgrams<true>(localdata[tid].second, 
+								makeStrideIter(docIt, stride, docEnd), 
+								makeStrideIter(docEnd, stride, docEnd), 
+								vocabFreqs, vocabDf, validPairs, candMinCnt, candMinDf, maxNgrams + 1
+							);
+						}));
+					}
+
+					for (auto& f : futures) f.get();
+
+					auto r = parallelReduce(std::move(localdata), [&](LocalFwBw& dest, LocalFwBw&& src)
+					{
+						mergeNgramCounts(dest.first, std::move(src.first));
+						mergeNgramCounts(dest.second, std::move(src.second));
+					}, pool);
+
+					trieNodes = std::move(r.first);
+					trieNodesBw = std::move(r.second);
+				}
+				else
+				{
+					countNgrams<false>(trieNodes, 
+						docBegin, docEnd, 
+						vocabFreqs, vocabDf, validPairs, candMinCnt, candMinDf, maxNgrams + 1
+					);
+					countNgrams<true>(trieNodesBw, 
+						docBegin, docEnd, 
+						vocabFreqs, vocabDf, validPairs, candMinCnt, candMinDf, maxNgrams + 1
+					);
+				}
+			}
+
+			float totN = std::accumulate(vocabFreqs.begin(), vocabFreqs.end(), (size_t)0);
+			const float logTotN = std::log(totN);
+
+			// calculating PMIs
+			std::vector<label::Candidate> candidates;
+			for (auto& p : bigramCnt)
+			{
+				auto& bigram = p.first;
+				if (p.second < candMinCnt) continue;
+				if (bigramDf[bigram] < candMinDf) continue;
+				float npmi = std::log(p.second * totN
+					/ vocabFreqs[bigram.first] / vocabFreqs[bigram.second]);
+				npmi /= std::log(totN / p.second);
+				if (npmi < minNPMI) continue;
+
+				float rbe = branchingEntropy(trieNodes[0].getNext(bigram.first)->getNext(bigram.second), candMinCnt);
+				float lbe = branchingEntropy(trieNodesBw[0].getNext(bigram.second)->getNext(bigram.first), candMinCnt);
+				float nbe = std::sqrt(rbe * lbe) / std::log(p.second);
+				if (nbe < minNBE) continue;
+				candidates.emplace_back(npmi * nbe, bigram.first, bigram.second);
+				candidates.back().cf = p.second;
+				candidates.back().df = bigramDf[bigram];
+			}
+
+			if (maxNgrams > 2)
+			{
+				std::vector<Vid> rkeys;
+				trieNodes[0].traverse_with_keys([&](const TrieEx<Vid, size_t>* node, const std::vector<Vid>& rkeys)
+				{
+					if (rkeys.size() <= 2 || rkeys.size() < minNgrams || rkeys.size() > maxNgrams || node->val < candMinCnt) return;
+					auto npmi = std::log((float)node->val) - logTotN;
+					for (auto k : rkeys)
+					{
+						npmi += logTotN - std::log((float)vocabFreqs[k]);
+					}
+					npmi /= (logTotN - std::log((float)node->val)) * (rkeys.size() - 1);
+					if (npmi < minNPMI) return;
+
+					float rbe = branchingEntropy(node, candMinCnt);
+					float lbe = branchingEntropy(trieNodesBw[0].findNode(rkeys.rbegin(), rkeys.rend()), candMinCnt);
+					float nbe = std::sqrt(rbe * lbe) / std::log(node->val);
+					if (nbe < minNBE) return;
+					candidates.emplace_back(npmi * nbe, rkeys);
+					candidates.back().cf = node->val;
+				}, rkeys);
+			}
+
+			std::sort(candidates.begin(), candidates.end(), [](const label::Candidate& a, const label::Candidate& b)
+			{
+				return a.score > b.score;
+			});
+			if (candidates.size() > maxCandidates) candidates.erase(candidates.begin() + maxCandidates, candidates.end());
+			return candidates;
+		}
+	}
+}
diff --git a/src/TopicModel/CTModel.hpp b/src/TopicModel/CTModel.hpp
index 8e64068..aff3a3b 100644
--- a/src/TopicModel/CTModel.hpp
+++ b/src/TopicModel/CTModel.hpp
@@ -120,8 +120,8 @@ namespace tomoto
 			}*/
 		}
 
-		template<typename _DocIter>
-		void sampleGlobalLevel(ThreadPool* pool, _ModelState* localData, _RandGen* rgs, _DocIter first, _DocIter last) const
+		template<GlobalSampler _gs, typename _DocIter>
+		void sampleGlobalLevel(ThreadPool* pool, _ModelState*, _RandGen* rgs, _DocIter first, _DocIter last) const
 		{
 			if (this->globalStep < this->burnIn || !this->optimInterval || (this->globalStep + 1) % this->optimInterval != 0) return;
 
diff --git a/src/TopicModel/DT.h b/src/TopicModel/DT.h
index b5995c8..b623007 100644
--- a/src/TopicModel/DT.h
+++ b/src/TopicModel/DT.h
@@ -11,7 +11,7 @@ namespace tomoto
 		using DocumentLDA<_tw>::DocumentLDA;
 
 		uint64_t timepoint = 0;
-		ShareableVector<Float> eta;
+		ShareableMatrix<Float, -1, 1> eta;
 		sample::AliasMethod<> aliasTable;
 
 		DEFINE_SERIALIZER_AFTER_BASE_WITH_VERSION(BaseDocument, 0, timepoint);
diff --git a/src/TopicModel/DTModel.hpp b/src/TopicModel/DTModel.hpp
index 7cd5951..40f59d4 100644
--- a/src/TopicModel/DTModel.hpp
+++ b/src/TopicModel/DTModel.hpp
@@ -20,6 +20,7 @@ namespace tomoto
 
 		Eigen::Matrix<WeightType, -1, -1> numByTopic; // Dim: (Topic, Time)
 		Eigen::Matrix<WeightType, -1, -1> numByTopicWord; // Dim: (Topic * Time, Vocabs)
+		//ShareableMatrix<WeightType, -1, -1> numByTopicWord; // Dim: (Topic * Time, Vocabs)
 		DEFINE_SERIALIZER(numByTopic, numByTopicWord);
 	};
 
@@ -139,8 +140,6 @@ namespace tomoto
 		template<ParallelScheme _ps, typename _ExtraDocData>
 		void mergeState(ThreadPool& pool, _ModelState& globalState, _ModelState& tState, _ModelState* localData, _RandGen*, const _ExtraDocData& edd) const
 		{
-			std::vector<std::future<void>> res;
-
 			if (_ps == ParallelScheme::copy_merge)
 			{
 				tState = globalState;
@@ -157,17 +156,10 @@ namespace tomoto
 				}
 				Eigen::Map<Eigen::Matrix<WeightType, -1, 1>>{ globalState.numByTopic.data(), globalState.numByTopic.size() }
 					= globalState.numByTopicWord.rowwise().sum();
-
-				for (size_t i = 0; i < pool.getNumWorkers(); ++i)
-				{
-					res.emplace_back(pool.enqueue([&, i](size_t)
-					{
-						localData[i] = globalState;
-					}));
-				}
 			}
 			else if (_ps == ParallelScheme::partition)
 			{
+				std::vector<std::future<void>> res;
 				res = pool.enqueueToAll([&](size_t partitionId)
 				{
 					size_t b = partitionId ? edd.vChunkOffset[partitionId - 1] : 0,
@@ -175,7 +167,6 @@ namespace tomoto
 					globalState.numByTopicWord.block(0, b, globalState.numByTopicWord.rows(), e - b) = localData[partitionId].numByTopicWord;
 				});
 				for (auto& r : res) r.get();
-				res.clear();
 
 				// make all count being positive
 				if (_tw != TermWeight::one)
@@ -184,17 +175,11 @@ namespace tomoto
 				}
 				Eigen::Map<Eigen::Matrix<WeightType, -1, 1>>{ globalState.numByTopic.data(), globalState.numByTopic.size() }
 					= globalState.numByTopicWord.rowwise().sum();
-
-				res = pool.enqueueToAll([&](size_t threadId)
-				{
-					localData[threadId].numByTopic = globalState.numByTopic;
-				});
 			}
-			for (auto& r : res) r.get();
 		}
 
 		template<typename _DocIter>
-		void sampleGlobalLevel(ThreadPool* pool, _ModelState* localData, _RandGen* rgs, _DocIter first, _DocIter last)
+		void _sampleGlobalLevel(ThreadPool* pool, _ModelState*, _RandGen* rgs, _DocIter first, _DocIter last)
 		{
 			const auto K = this->K;
 			const Float eps = shapeA * (std::pow(shapeB + 1 + this->globalStep, -shapeC));
@@ -313,10 +298,10 @@ namespace tomoto
 			alphas = newAlphas;
 		}
 
-		template<typename _DocIter>
+		template<GlobalSampler _gs, typename _DocIter>
 		void sampleGlobalLevel(ThreadPool* pool, _ModelState* localData, _RandGen* rgs, _DocIter first, _DocIter last) const
 		{
-			// do nothing
+			if (_gs != GlobalSampler::inference) return const_cast<DerivedClass*>(this)->_sampleGlobalLevel(pool, localData, rgs, first, last);
 		}
 
 		void optimizeParameters(ThreadPool& pool, _ModelState* localData, _RandGen* rgs)
@@ -343,11 +328,11 @@ namespace tomoto
 			BaseClass::prepareDoc(doc, docId, wordSize);
 			if (docId == (size_t)-1)
 			{
-				doc.eta.init(nullptr, this->K);
+				doc.eta.init(nullptr, this->K, 1);
 			}
 			else
 			{
-				doc.eta.init((Float*)etaByDoc.col(docId).data(), this->K);
+				doc.eta.init((Float*)etaByDoc.col(docId).data(), this->K, 1);
 			}
 		}
 
@@ -427,7 +412,7 @@ namespace tomoto
 				numDocsByTime[doc.timepoint]++;
 				if (!initDocs)
 				{
-					doc.eta.init((Float*)etaByDoc.col(docId++).data(), this->K);
+					doc.eta.init((Float*)etaByDoc.col(docId++).data(), this->K, 1);
 				}
 			}
 
diff --git a/src/TopicModel/HDPModel.hpp b/src/TopicModel/HDPModel.hpp
index afd87ce..65935c3 100644
--- a/src/TopicModel/HDPModel.hpp
+++ b/src/TopicModel/HDPModel.hpp
@@ -96,7 +96,7 @@ namespace tomoto
 				ld.numTableByTopic.tail(newSize - oldSize).setZero();
 				ld.numByTopic.conservativeResize(newSize);
 				ld.numByTopic.tail(newSize - oldSize).setZero();
-				ld.numByTopicWord.conservativeResize(newSize, Eigen::NoChange);
+				ld.numByTopicWord.conservativeResize(newSize, V);
 				ld.numByTopicWord.block(oldSize, 0, newSize - oldSize, V).setZero();
 			}
 			else
@@ -155,7 +155,7 @@ namespace tomoto
 			if (_inc > 0 && tid >= doc.numByTopic.size())
 			{
 				size_t oldSize = doc.numByTopic.size();
-				doc.numByTopic.conservativeResize(tid + 1);
+				doc.numByTopic.conservativeResize(tid + 1, 1);
 				doc.numByTopic.tail(tid + 1 - oldSize).setZero();
 			}
 			constexpr bool _dec = _inc < 0 && _tw != TermWeight::one;
@@ -282,7 +282,7 @@ namespace tomoto
 						auto& doc = this->docs[j];
 						if (doc.numByTopic.size() >= K) continue;
 						size_t oldSize = doc.numByTopic.size();
-						doc.numByTopic.conservativeResize(K);
+						doc.numByTopic.conservativeResize(K, 1);
 						doc.numByTopic.tail(K - oldSize).setZero();
 					}
 				}, this->docs.size() * i / pool.getNumWorkers(), this->docs.size() * (i + 1) / pool.getNumWorkers()));
@@ -293,7 +293,6 @@ namespace tomoto
 		template<ParallelScheme _ps, typename _ExtraDocData>
 		void mergeState(ThreadPool& pool, _ModelState& globalState, _ModelState& tState, _ModelState* localData, _RandGen*, const _ExtraDocData& edd) const
 		{
-			std::vector<std::future<void>> res;
 			const size_t V = this->realV;
 			auto K = this->K;
 
@@ -303,7 +302,7 @@ namespace tomoto
 				globalState.numByTopic.conservativeResize(K);
 				globalState.numByTopic.tail(K - oldSize).setZero();
 				globalState.numTableByTopic.resize(K);
-				globalState.numByTopicWord.conservativeResize(K, Eigen::NoChange);
+				globalState.numByTopicWord.conservativeResize(K, V);
 				globalState.numByTopicWord.block(oldSize, 0, K - oldSize, V).setZero();
 			}
 
@@ -321,7 +320,7 @@ namespace tomoto
 			if (_tw != TermWeight::one)
 			{
 				globalState.numByTopic = globalState.numByTopic.cwiseMax(0);
-				globalState.numByTopicWord = globalState.numByTopicWord.cwiseMax(0);
+				globalState.numByTopicWord.matrix() = globalState.numByTopicWord.cwiseMax(0);
 			}
 
 
@@ -334,15 +333,6 @@ namespace tomoto
 				}
 			}
 			globalState.totalTable = globalState.numTableByTopic.sum();
-			
-			for (size_t i = 0; i < pool.getNumWorkers(); ++i)
-			{
-				res.emplace_back(pool.enqueue([&, this, i](size_t threadId)
-				{
-					localData[i] = globalState;
-				}));
-			}
-			for (auto& r : res) r.get();
 		}
 
 		/* this LL calculation is based on https://github.com/blei-lab/hdp/blob/master/hdp/state.cpp */
@@ -400,13 +390,14 @@ namespace tomoto
 			{
 				this->globalState.numByTopic = Eigen::Matrix<WeightType, -1, 1>::Zero(K);
 				this->globalState.numTableByTopic = Eigen::Matrix<int32_t, -1, 1>::Zero(K);
-				this->globalState.numByTopicWord = Eigen::Matrix<WeightType, -1, -1>::Zero(K, V);
+				//this->globalState.numByTopicWord = Eigen::Matrix<WeightType, -1, -1>::Zero(K, V);
+				this->globalState.numByTopicWord.init(nullptr, K, V);
 			}
 		}
 
 		void prepareDoc(_DocType& doc, size_t docId, size_t wordSize) const
 		{
-			doc.numByTopic.init(nullptr, this->K);
+			doc.numByTopic.init(nullptr, this->K, 1);
 			doc.numTopicByTable.clear();
 			doc.Zs = tvector<Tid>(wordSize);
 			if (_tw != TermWeight::one) doc.wordWeights.resize(wordSize);
@@ -577,7 +568,7 @@ namespace tomoto
 	template<typename _TopicModel>
 	void DocumentHDP<_tw>::update(WeightType * ptr, const _TopicModel & mdl)
 	{
-		this->numByTopic.init(ptr, mdl.getK());
+		this->numByTopic.init(ptr, mdl.getK(), 1);
 		for (size_t i = 0; i < this->Zs.size(); ++i)
 		{
 			if (this->words[i] >= mdl.getV()) continue;
diff --git a/src/TopicModel/HLDAModel.hpp b/src/TopicModel/HLDAModel.hpp
index f838b6a..6736ae1 100644
--- a/src/TopicModel/HLDAModel.hpp
+++ b/src/TopicModel/HLDAModel.hpp
@@ -119,19 +119,26 @@ namespace tomoto
 
 			DEFINE_SERIALIZER(nodes, levelBlocks);
 
-			template<bool _MakeNewPath = true>
+			template<bool _makeNewPath = true>
 			void calcNodeLikelihood(Float gamma, size_t levelDepth)
 			{
 				nodeLikelihoods.resize(nodes.size());
 				nodeLikelihoods.array() = -INFINITY;
-				updateNodeLikelihood<_MakeNewPath>(gamma, levelDepth, &nodes[0]);
+				updateNodeLikelihood<_makeNewPath>(gamma, levelDepth, &nodes[0]);
+				if (!_makeNewPath)
+				{
+					for (size_t i = 0; i < levelBlocks.size(); ++i)
+					{
+						if (levelBlocks[i] < levelDepth - 1) nodeLikelihoods.segment((i + 1) * blockSize, blockSize).array() = -INFINITY;
+					}
+				}
 			}
 
-			template<bool _MakeNewPath = true>
+			template<bool _makeNewPath = true>
 			void updateNodeLikelihood(Float gamma, size_t levelDepth, NCRPNode* node, Float weight = 0)
 			{
 				size_t idx = node - nodes.data();
-				const Float pNewNode = _MakeNewPath ? log(gamma / (node->numCustomers + gamma)) : -INFINITY;
+				const Float pNewNode = _makeNewPath ? log(gamma / (node->numCustomers + gamma)) : -INFINITY;
 				nodeLikelihoods[idx] = weight + (((size_t)node->level < levelDepth - 1) ? pNewNode : 0);
 				for(auto * child = node->getChild(); child; child = child->getSibling())
 				{
@@ -187,7 +194,7 @@ namespace tomoto
 				std::vector<std::future<void>> futures;
 				futures.reserve(levelBlocks.size());
 
-				auto calc = [this, eta, realV, &doc, &ld](size_t threadId, size_t b)
+				auto calc = [&, eta, realV](size_t threadId, size_t b)
 				{
 					Float cnt = 0;
 					Vid prevWord = -1;
@@ -284,7 +291,7 @@ namespace tomoto
 					size_t oldSize = ld.numByTopic.rows();
 					size_t newSize = std::max(nodes.size(), ((oldSize + oldSize / 2 + 7) / 8) * 8);
 					ld.numByTopic.conservativeResize(newSize);
-					ld.numByTopicWord.conservativeResize(newSize, Eigen::NoChange);
+					ld.numByTopicWord.conservativeResize(newSize, ld.numByTopicWord.cols());
 					ld.numByTopic.segment(oldSize, newSize - oldSize).setZero();
 					ld.numByTopicWord.block(oldSize, 0, newSize - oldSize, ld.numByTopicWord.cols()).setZero();
 				}
@@ -317,13 +324,13 @@ namespace tomoto
 		typename _Derived = void,
 		typename _DocType = DocumentHLDA<_tw>,
 		typename _ModelState = ModelStateHLDA<_tw>>
-	class HLDAModel : public LDAModel<_tw, _RandGen, flags::shared_state, _Interface,
+	class HLDAModel : public LDAModel<_tw, _RandGen, flags::partitioned_multisampling, _Interface,
 		typename std::conditional<std::is_same<_Derived, void>::value, HLDAModel<_tw, _RandGen>, _Derived>::type,
 		_DocType, _ModelState>
 	{
 	protected:
 		using DerivedClass = typename std::conditional<std::is_same<_Derived, void>::value, HLDAModel<_tw, _RandGen>, _Derived>::type;
-		using BaseClass = LDAModel<_tw, _RandGen, flags::shared_state, _Interface, DerivedClass, _DocType, _ModelState>;
+		using BaseClass = LDAModel<_tw, _RandGen, flags::partitioned_multisampling, _Interface, DerivedClass, _DocType, _ModelState>;
 		friend BaseClass;
 		friend typename BaseClass::BaseClass;
 		using WeightType = typename BaseClass::WeightType;
@@ -341,11 +348,11 @@ namespace tomoto
 		}
 
 		// Words of all documents should be sorted by ascending order.
-		template<bool _MakeNewPath = true>
+		template<GlobalSampler _gs>
 		void samplePathes(_DocType& doc, ThreadPool* pool, _ModelState& ld, _RandGen& rgs) const
 		{
-			if(_MakeNewPath) ld.nt->nodes[doc.path.back()].dropPathOne();
-			ld.nt->template calcNodeLikelihood<_MakeNewPath>(gamma, this->K);
+			if(_gs != GlobalSampler::inference) ld.nt->nodes[doc.path.back()].dropPathOne();
+			ld.nt->template calcNodeLikelihood<_gs == GlobalSampler::train>(gamma, this->K);
 
 			std::vector<Float> newTopicWeights(this->K - 1);
 			std::vector<WeightType> cntByLevel(this->K);
@@ -355,7 +362,7 @@ namespace tomoto
 				if (doc.words[w] >= this->realV) break;
 				addWordToOnlyLocal<-1>(ld, doc, w, doc.words[w], doc.Zs[w]);
 
-				if (_MakeNewPath)
+				if (_gs == GlobalSampler::train)
 				{
 					if (doc.words[w] != prevWord)
 					{
@@ -371,7 +378,7 @@ namespace tomoto
 				}
 			}
 
-			if (_MakeNewPath)
+			if (_gs == GlobalSampler::train)
 			{
 				for (size_t l = 1; l < this->K; ++l)
 				{
@@ -386,7 +393,7 @@ namespace tomoto
 			size_t newPath = sample::sampleFromDiscreteAcc(ld.nt->nodeLikelihoods.data(),
 				ld.nt->nodeLikelihoods.data() + ld.nt->nodeLikelihoods.size(), rgs);
 
-			if(_MakeNewPath) newPath = ld.nt->template generateLeafNode<_tw>(newPath, this->K, ld);
+			if(_gs == GlobalSampler::train) newPath = ld.nt->template generateLeafNode<_tw>(newPath, this->K, ld);
 			doc.path.back() = newPath;
 			for (size_t l = this->K - 2; l > 0; --l)
 			{
@@ -398,7 +405,7 @@ namespace tomoto
 				if (doc.words[w] >= this->realV) break;
 				addWordToOnlyLocal<1>(ld, doc, w, doc.words[w], doc.Zs[w]);
 			}
-			if (_MakeNewPath) ld.nt->nodes[doc.path.back()].addPathOne();
+			if (_gs != GlobalSampler::inference) ld.nt->nodes[doc.path.back()].addPathOne();
 		}
 
 		template<int _inc>
@@ -426,6 +433,7 @@ namespace tomoto
 		template<bool _asymEta>
 		Float* getZLikelihoods(_ModelState& ld, const _DocType& doc, size_t docId, size_t vid) const
 		{
+			if (_asymEta) THROW_ERROR_WITH_INFO(exception::Unimplemented, "Unimplemented features");
 			const size_t V = this->realV;
 			assert(vid < V);
 			auto& zLikelihood = ld.zLikelihood;
@@ -439,50 +447,14 @@ namespace tomoto
 			return &zLikelihood[0];
 		}
 
-		void sampleTopics(_DocType& doc, size_t docId, _ModelState& ld, _RandGen& rgs) const
-		{
-			for (size_t w = 0; w < doc.words.size(); ++w)
-			{
-				if (doc.words[w] >= this->realV) continue;
-				addWordTo<-1>(ld, doc, w, doc.words[w], doc.Zs[w]);
-				Float* dist;
-				if (this->etaByTopicWord.size())
-				{
-					THROW_ERROR_WITH_INFO(exception::Unimplemented, "Unimplemented features");
-				}
-				else
-				{
-					dist = static_cast<const DerivedClass*>(this)->template
-						getZLikelihoods<false>(ld, doc, docId, doc.words[w]);
-				}
-				doc.Zs[w] = sample::sampleFromDiscreteAcc(dist, dist + this->K, rgs);
-				addWordTo<1>(ld, doc, w, doc.words[w], doc.Zs[w]);
-			}
-		}
-
-		template<ParallelScheme _ps, bool _infer, typename _ExtraDocData>
-		void sampleDocument(_DocType& doc, const _ExtraDocData& edd, size_t docId, _ModelState& ld, _RandGen& rgs, size_t iterationCnt, size_t partitionId = 0) const
-		{
-			sampleTopics(doc, docId, ld, rgs);
-		}
-
-		template<typename _DocIter>
-		void sampleGlobalLevel(ThreadPool* pool, _ModelState* localData, _RandGen* rgs, _DocIter first, _DocIter last)
+		template<GlobalSampler _gs, typename _DocIter>
+		void sampleGlobalLevel(ThreadPool* pool, _ModelState* globalData, _RandGen* rgs, _DocIter first, _DocIter last) const
 		{
 			for (auto doc = first; doc != last; ++doc)
 			{
-				samplePathes<>(*doc, pool, *localData, rgs[0]);
-			}
-			localData->nt->markEmptyBlocks();
-		}
-
-		template<typename _DocIter>
-		void sampleGlobalLevel(ThreadPool* pool, _ModelState* localData, _RandGen* rgs, _DocIter first, _DocIter last) const
-		{
-			for (auto doc = first; doc != last; ++doc)
-			{
-				samplePathes<false>(*doc, pool, *localData, rgs[0]);
+				samplePathes<_gs>(*doc, pool, *globalData, rgs[0]);
 			}
+			if (_gs != GlobalSampler::inference) globalData->nt->markEmptyBlocks();
 		}
 
 		template<typename _DocIter>
@@ -539,7 +511,8 @@ namespace tomoto
 			if (initDocs)
 			{
 				this->globalState.numByTopic = Eigen::Matrix<WeightType, -1, 1>::Zero(this->K);
-				this->globalState.numByTopicWord = Eigen::Matrix<WeightType, -1, -1>::Zero(this->K, V);
+				//this->globalState.numByTopicWord = Eigen::Matrix<WeightType, -1, -1>::Zero(this->K, V);
+				this->globalState.numByTopicWord.init(nullptr, this->K, V);
 				this->globalState.nt->nodes.resize(detail::NodeTrees::blockSize);
 			}
 		}
@@ -547,7 +520,7 @@ namespace tomoto
 		void prepareDoc(_DocType& doc, size_t docId, size_t wordSize) const
 		{
 			sortAndWriteOrder(doc.words, doc.wOrder);
-			doc.numByTopic.init(nullptr, this->K);
+			doc.numByTopic.init(nullptr, this->K, 1);
 			doc.Zs = tvector<Tid>(wordSize);
 			doc.path.resize(this->K);
 			for (size_t l = 0; l < this->K; ++l) doc.path[l] = l;
@@ -595,6 +568,31 @@ namespace tomoto
 			return cnt;
 		}
 
+		template<ParallelScheme _ps>
+		void distributeMergedState(ThreadPool& pool, _ModelState& globalState, _ModelState* localData) const
+		{
+			std::vector<std::future<void>> res;
+			if (_ps == ParallelScheme::copy_merge)
+			{
+				for (size_t i = 0; i < pool.getNumWorkers(); ++i)
+				{
+					res.emplace_back(pool.enqueue([&, i](size_t)
+					{
+						localData[i] = globalState;
+					}));
+				}
+			}
+			else if (_ps == ParallelScheme::partition)
+			{
+				res = pool.enqueueToAll([&](size_t threadId)
+				{
+					localData[threadId].numByTopicWord.init((WeightType*)globalState.numByTopicWord.data(), globalState.numByTopicWord.rows(), globalState.numByTopicWord.cols());
+					localData[threadId].numByTopic = globalState.numByTopic;
+				});
+			}
+			for (auto& r : res) r.get();
+		}
+
 	public:
 		DEFINE_SERIALIZER_AFTER_BASE_WITH_VERSION(BaseClass, 0, gamma);
 		DEFINE_TAGGED_SERIALIZER_AFTER_BASE_WITH_VERSION(BaseClass, 1, 0x00010001, gamma);
@@ -671,7 +669,7 @@ namespace tomoto
 	template<typename _TopicModel>
 	inline void DocumentHLDA<_tw>::update(WeightType * ptr, const _TopicModel & mdl)
 	{
-		this->numByTopic.init(ptr, mdl.getLevelDepth());
+		this->numByTopic.init(ptr, mdl.getLevelDepth(), 1);
 		for (size_t i = 0; i < this->Zs.size(); ++i)
 		{
 			if (this->words[i] >= mdl.getV()) continue;
diff --git a/src/TopicModel/HPAModel.hpp b/src/TopicModel/HPAModel.hpp
index 1e31e88..0974e57 100644
--- a/src/TopicModel/HPAModel.hpp
+++ b/src/TopicModel/HPAModel.hpp
@@ -28,7 +28,8 @@ namespace tomoto
 		typename _Interface = IHPAModel,
 		typename _Derived = void,
 		typename _DocType = DocumentHPA<_tw>,
-		typename _ModelState = ModelStateHPA<_tw>>
+		typename _ModelState = ModelStateHPA<_tw>
+	>
 	class HPAModel : public LDAModel<_tw, _RandGen, 0, _Interface,
 		typename std::conditional<std::is_same<_Derived, void>::value, HPAModel<_tw, _RandGen, _Exclusive>, _Derived>::type,
 		_DocType, _ModelState>
@@ -250,8 +251,6 @@ namespace tomoto
 		template<ParallelScheme _ps, typename _ExtraDocData>
 		void mergeState(ThreadPool& pool, _ModelState& globalState, _ModelState& tState, _ModelState* localData, _RandGen*, const _ExtraDocData& edd) const
 		{
-			std::vector<std::future<void>> res;
-
 			tState = globalState;
 			globalState = localData[0];
 			for (size_t i = 1; i < pool.getNumWorkers(); ++i)
@@ -276,15 +275,6 @@ namespace tomoto
 				globalState.numByTopicWord[1] = globalState.numByTopicWord[1].cwiseMax(0);
 				globalState.numByTopicWord[2] = globalState.numByTopicWord[2].cwiseMax(0);
 			}
-
-			for (size_t i = 0; i < pool.getNumWorkers(); ++i)
-			{
-				res.emplace_back(pool.enqueue([&, this, i](size_t threadId)
-				{
-					localData[i] = globalState;
-				}));
-			}
-			for (auto& r : res) r.get();
 		}
 
 		std::vector<uint64_t> _getTopicsCount() const
@@ -379,7 +369,7 @@ namespace tomoto
 
 		void prepareDoc(_DocType& doc, size_t docId, size_t wordSize) const
 		{
-			doc.numByTopic.init(nullptr, this->K + 1);
+			doc.numByTopic.init(nullptr, this->K + 1, 1);
 			doc.numByTopic1_2 = Eigen::Matrix<WeightType, -1, -1>::Zero(this->K, K2 + 1);
 			doc.Zs = tvector<Tid>(wordSize);
 			doc.Z2s = tvector<Tid>(wordSize);
@@ -575,7 +565,7 @@ namespace tomoto
 	template<typename _TopicModel>
 	void DocumentHPA<_tw>::update(WeightType * ptr, const _TopicModel & mdl)
 	{
-		this->numByTopic.init(ptr, mdl.getK() + 1);
+		this->numByTopic.init(ptr, mdl.getK() + 1, 1);
 		this->numByTopic1_2 = Eigen::Matrix<WeightType, -1, -1>::Zero(mdl.getK(), mdl.getK2() + 1);
 		for (size_t i = 0; i < this->Zs.size(); ++i)
 		{
diff --git a/src/TopicModel/LDA.h b/src/TopicModel/LDA.h
index f62264d..718a63b 100644
--- a/src/TopicModel/LDA.h
+++ b/src/TopicModel/LDA.h
@@ -5,32 +5,67 @@ namespace tomoto
 {
     enum class TermWeight { one, idf, pmi, size };
 
-	template<typename _Scalar>
-	struct ShareableVector : Eigen::Map<Eigen::Matrix<_Scalar, -1, 1>>
+	template<typename _Scalar, Eigen::Index _rows, Eigen::Index _cols>
+	struct ShareableMatrix : Eigen::Map<Eigen::Matrix<_Scalar, _rows, _cols>>
 	{
-		Eigen::Matrix<_Scalar, -1, 1> ownData;
-		ShareableVector(_Scalar* ptr = nullptr, Eigen::Index len = 0) 
-			: Eigen::Map<Eigen::Matrix<_Scalar, -1, 1>>(nullptr, 0)
+		using BaseType = Eigen::Map<Eigen::Matrix<_Scalar, _rows, _cols>>;
+		Eigen::Matrix<_Scalar, _rows, _cols> ownData;
+
+		ShareableMatrix(_Scalar* ptr = nullptr, Eigen::Index rows = 0, Eigen::Index cols = 0) 
+			: BaseType(nullptr, _rows != -1 ? _rows : 0, _cols != -1 ? _cols : 0)
 		{
-			init(ptr, len);
+			init(ptr, rows, cols);
 		}
 
-		void init(_Scalar* ptr, Eigen::Index len)
+		ShareableMatrix(const ShareableMatrix& o)
+			: BaseType(nullptr, _rows != -1 ? _rows : 0, _cols != -1 ? _cols : 0), ownData{ o.ownData }
 		{
-			if (!ptr && len)
+			if (o.ownData.data())
+			{
+				new (this) BaseType(ownData.data(), ownData.rows(), ownData.cols());
+			}
+			else
 			{
-				ownData = Eigen::Matrix<_Scalar, -1, 1>::Zero(len);
+				new (this) BaseType((_Scalar*)o.data(), o.rows(), o.cols());
+			}
+		}
+
+		ShareableMatrix(ShareableMatrix&& o) = default;
+
+		ShareableMatrix& operator=(const ShareableMatrix& o)
+		{
+			if (o.ownData.data())
+			{
+				ownData = o.ownData;
+				new (this) BaseType(ownData.data(), ownData.rows(), ownData.cols());
+			}
+			else
+			{
+				new (this) BaseType((_Scalar*)o.data(), o.rows(), o.cols());
+			}
+			return *this;
+		}
+
+		ShareableMatrix& operator=(ShareableMatrix&& o) = default;
+
+		void init(_Scalar* ptr, Eigen::Index rows, Eigen::Index cols)
+		{
+			if (!ptr && rows && cols)
+			{
+				ownData = Eigen::Matrix<_Scalar, _rows, _cols>::Zero(_rows != -1 ? _rows : rows, _cols != -1 ? _cols : cols);
 				ptr = ownData.data();
 			}
-			// is this the best way??
-			this->m_data = ptr;
-			((Eigen::internal::variable_if_dynamic<Eigen::Index, -1>*)&this->m_rows)->setValue(len);
+			else
+			{
+				ownData = Eigen::Matrix<_Scalar, _rows, _cols>{};
+			}
+			new (this) BaseType(ptr, _rows != -1 ? _rows : rows, _cols != -1 ? _cols : cols);
 		}
 
-		void conservativeResize(size_t newSize)
+		void conservativeResize(size_t newRows, size_t newCols)
 		{
-			ownData.conservativeResize(newSize);
-			init(ownData.data(), ownData.size());
+			ownData.conservativeResize(_rows != -1 ? _rows : newRows, _cols != -1 ? _cols : newCols);
+			new (this) BaseType(ownData.data(), ownData.rows(), ownData.cols());
 		}
 
 		void becomeOwner()
@@ -38,9 +73,26 @@ namespace tomoto
 			if (ownData.data() != this->m_data)
 			{
 				ownData = *this;
-				init(ownData.data(), ownData.size());
+				new (this) BaseType(ownData.data(), ownData.rows(), ownData.cols());
 			}
 		}
+
+		void serializerRead(std::istream& istr)
+		{
+			uint32_t rows = serializer::readFromStream<uint32_t>(istr);
+			uint32_t cols = serializer::readFromStream<uint32_t>(istr);
+			init(nullptr, rows, cols);
+			if (!istr.read((char*)this->data(), sizeof(_Scalar) * this->size()))
+				throw std::ios_base::failure(std::string("reading type '") + typeid(_Scalar).name() + std::string("' is failed"));
+		}
+
+		void serializerWrite(std::ostream& ostr) const
+		{
+			serializer::writeToStream<uint32_t>(ostr, (uint32_t)this->rows());
+			serializer::writeToStream<uint32_t>(ostr, (uint32_t)this->cols());
+			if (!ostr.write((const char*)this->data(), sizeof(_Scalar) * this->size()))
+				throw std::ios_base::failure(std::string("writing type '") + typeid(_Scalar).name() + std::string("' is failed"));
+		}
 	};
 
 	template<typename _Base, TermWeight _tw>
@@ -85,7 +137,7 @@ namespace tomoto
 
 		tvector<Tid> Zs;
 		tvector<Float> wordWeights;
-		ShareableVector<WeightType> numByTopic;
+		ShareableMatrix<WeightType, -1, 1> numByTopic;
 
 		DEFINE_SERIALIZER_AFTER_BASE_WITH_VERSION(DocumentBase, 0, Zs, wordWeights);
 		DEFINE_TAGGED_SERIALIZER_AFTER_BASE_WITH_VERSION(DocumentBase, 1, 0x00010001, Zs, wordWeights);
diff --git a/src/TopicModel/LDACVB0Model.hpp b/src/TopicModel/LDACVB0Model.hpp
index 6b55e87..403b4a6 100644
--- a/src/TopicModel/LDACVB0Model.hpp
+++ b/src/TopicModel/LDACVB0Model.hpp
@@ -163,7 +163,7 @@ namespace tomoto
 			{
 				res.emplace_back(pool.enqueue([&, this, ch, chStride](size_t threadId)
 				{
-					forRandom((this->docs.size() - 1 - ch) / chStride + 1, rgs[threadId](), [&, this](size_t id)
+					forShuffled((this->docs.size() - 1 - ch) / chStride + 1, rgs[threadId](), [&, this](size_t id)
 					{
 						static_cast<DerivedClass*>(this)->template sampleDocument<ParallelScheme::copy_merge>(
 							this->docs[id * chStride + ch], 0, id * chStride + ch,
diff --git a/src/TopicModel/LDAModel.hpp b/src/TopicModel/LDAModel.hpp
index 27ddc17..760c84a 100644
--- a/src/TopicModel/LDAModel.hpp
+++ b/src/TopicModel/LDAModel.hpp
@@ -58,7 +58,8 @@ namespace tomoto
 
 		Eigen::Matrix<Float, -1, 1> zLikelihood;
 		Eigen::Matrix<WeightType, -1, 1> numByTopic; // Dim: (Topic, 1)
-		Eigen::Matrix<WeightType, -1, -1> numByTopicWord; // Dim: (Topic, Vocabs)
+		//Eigen::Matrix<WeightType, -1, -1> numByTopicWord; // Dim: (Topic, Vocabs)
+		ShareableMatrix<WeightType, -1, -1> numByTopicWord; // Dim: (Topic, Vocabs)
 		DEFINE_SERIALIZER(numByTopic, numByTopicWord);
 	};
 
@@ -137,7 +138,8 @@ namespace tomoto
 		typename _Interface,
 		typename _Derived,
 		typename _DocType,
-		typename _ModelState>
+		typename _ModelState
+	>
 	class HDPModel;
 
 	template<TermWeight _tw, typename _RandGen,
@@ -145,7 +147,8 @@ namespace tomoto
 		typename _Interface = ILDAModel,
 		typename _Derived = void, 
 		typename _DocType = DocumentLDA<_tw>,
-		typename _ModelState = ModelStateLDA<_tw>>
+		typename _ModelState = ModelStateLDA<_tw>
+	>
 	class LDAModel : public TopicModel<_RandGen, _Flags, _Interface,
 		typename std::conditional<std::is_same<_Derived, void>::value, LDAModel<_tw, _RandGen, _Flags>, _Derived>::type,
 		_DocType, _ModelState>,
@@ -306,25 +309,23 @@ namespace tomoto
 				e = edd.chunkOffsetByDoc(partitionId + 1, docId);
 			}
 
-			size_t vOffset = (_ps == ParallelScheme::partition && partitionId) ? edd.vChunkOffset[partitionId - 1] : 0;
-
 			for (size_t w = b; w < e; ++w)
 			{
 				if (doc.words[w] >= this->realV) continue;
-				addWordTo<-1>(ld, doc, w, doc.words[w] - vOffset, doc.Zs[w]);
+				static_cast<const DerivedClass*>(this)->template addWordTo<-1>(ld, doc, w, doc.words[w], doc.Zs[w]);
 				Float* dist;
 				if (etaByTopicWord.size())
 				{
 					dist = static_cast<const DerivedClass*>(this)->template
-						getZLikelihoods<true>(ld, doc, docId, doc.words[w] - vOffset);
+						getZLikelihoods<true>(ld, doc, docId, doc.words[w]);
 				}
 				else
 				{
 					dist = static_cast<const DerivedClass*>(this)->template
-						getZLikelihoods<false>(ld, doc, docId, doc.words[w] - vOffset);
+						getZLikelihoods<false>(ld, doc, docId, doc.words[w]);
 				}
 				doc.Zs[w] = sample::sampleFromDiscreteAcc(dist, dist + K, rgs);
-				addWordTo<1>(ld, doc, w, doc.words[w] - vOffset, doc.Zs[w]);
+				static_cast<const DerivedClass*>(this)->template addWordTo<1>(ld, doc, w, doc.words[w], doc.Zs[w]);
 			}
 		}
 
@@ -335,7 +336,7 @@ namespace tomoto
 			// single-threaded sampling
 			if (_ps == ParallelScheme::none)
 			{
-				forRandom((size_t)std::distance(docFirst, docLast), rgs[0](), [&](size_t id)
+				forShuffled((size_t)std::distance(docFirst, docLast), rgs[0](), [&](size_t id)
 				{
 					static_cast<const DerivedClass*>(this)->presampleDocument(docFirst[id], id, *localData, *rgs, this->globalStep);
 					static_cast<const DerivedClass*>(this)->template sampleDocument<_ps, _infer>(
@@ -344,7 +345,7 @@ namespace tomoto
 
 				});
 			}
-			// multi-threaded sampling on partition ad update into global
+			// multi-threaded sampling on partition and update into global
 			else if (_ps == ParallelScheme::partition)
 			{
 				const size_t chStride = pool.getNumWorkers();
@@ -353,7 +354,7 @@ namespace tomoto
 					res = pool.enqueueToAll([&, i, chStride](size_t partitionId)
 					{
 						size_t didx = (i + partitionId) % chStride;
-						forRandom(((size_t)std::distance(docFirst, docLast) + (chStride - 1) - didx) / chStride, rgs[partitionId](), [&](size_t id)
+						forShuffled(((size_t)std::distance(docFirst, docLast) + (chStride - 1) - didx) / chStride, rgs[partitionId](), [&](size_t id)
 						{
 							if (i == 0)
 							{
@@ -380,7 +381,7 @@ namespace tomoto
 				{
 					res.emplace_back(pool.enqueue([&, ch, chStride](size_t threadId)
 					{
-						forRandom(((size_t)std::distance(docFirst, docLast) + (chStride - 1) - ch) / chStride, rgs[threadId](), [&](size_t id)
+						forShuffled(((size_t)std::distance(docFirst, docLast) + (chStride - 1) - ch) / chStride, rgs[threadId](), [&](size_t id)
 						{
 							static_cast<const DerivedClass*>(this)->presampleDocument(
 								docFirst[id * chStride + ch], id * chStride + ch, 
@@ -396,6 +397,16 @@ namespace tomoto
 				for (auto& r : res) r.get();
 				res.clear();
 			}
+			else
+			{
+				throw std::runtime_error{ "Unsupported ParallelScheme" };
+			}
+		}
+
+		template<ParallelScheme _ps, bool _infer, typename _DocIter>
+		void performSamplingGlobal(ThreadPool* pool, _ModelState& globalState, _RandGen* rgs, 
+			_DocIter docFirst, _DocIter docLast) const
+		{
 		}
 
 		template<typename _DocIter, typename _ExtraDocData>
@@ -444,7 +455,8 @@ namespace tomoto
 				size_t b = partitionId ? edd.vChunkOffset[partitionId - 1] : 0,
 					e = edd.vChunkOffset[partitionId];
 
-				localData[partitionId].numByTopicWord = globalState.numByTopicWord.block(0, b, globalState.numByTopicWord.rows(), e - b);
+				//localData[partitionId].numByTopicWord.matrix() = globalState.numByTopicWord.block(0, b, globalState.numByTopicWord.rows(), e - b);
+				localData[partitionId].numByTopicWord.init((WeightType*)globalState.numByTopicWord.data(), globalState.numByTopicWord.rows(), globalState.numByTopicWord.cols());
 				localData[partitionId].numByTopic = globalState.numByTopic;
 				if (!localData[partitionId].zLikelihood.size()) localData[partitionId].zLikelihood = globalState.zLikelihood;
 			});
@@ -467,16 +479,29 @@ namespace tomoto
 		}
 
 		template<ParallelScheme _ps>
-		void trainOne(ThreadPool& pool, _ModelState* localData, _RandGen* rgs)
+		void trainOne(ThreadPool& pool, _ModelState* localData, _RandGen* rgs, bool freeze_topics = false)
 		{
 			std::vector<std::future<void>> res;
 			try
 			{
-				performSampling<_ps, false>(pool, localData, rgs, res, 
-					this->docs.begin(), this->docs.end(), eddTrain);
+				static_cast<DerivedClass*>(this)->template performSampling<_ps, false>(pool, localData, rgs, res,
+					this->docs.begin(), this->docs.end(), eddTrain
+				);
 				static_cast<DerivedClass*>(this)->updateGlobalInfo(pool, localData);
 				static_cast<DerivedClass*>(this)->template mergeState<_ps>(pool, this->globalState, this->tState, localData, rgs, eddTrain);
-				static_cast<DerivedClass*>(this)->template sampleGlobalLevel<>(&pool, localData, rgs, this->docs.begin(), this->docs.end());
+				static_cast<DerivedClass*>(this)->template performSamplingGlobal<_ps, false>(&pool, this->globalState, rgs, 
+					this->docs.begin(), this->docs.end()
+				);
+				
+				if(freeze_topics) static_cast<DerivedClass*>(this)->template sampleGlobalLevel<GlobalSampler::freeze_topics>(
+					&pool, &this->globalState, rgs, this->docs.begin(), this->docs.end()
+				);
+				else static_cast<DerivedClass*>(this)->template sampleGlobalLevel<GlobalSampler::train>(
+					&pool, &this->globalState, rgs, this->docs.begin(), this->docs.end()
+				);
+
+				static_cast<DerivedClass*>(this)->template distributeMergedState<_ps>(pool, this->globalState, localData);
+				
 				if (this->globalStep >= this->burnIn && optimInterval && (this->globalStep + 1) % optimInterval == 0)
 				{
 					static_cast<DerivedClass*>(this)->optimizeParameters(pool, localData, rgs);
@@ -503,8 +528,6 @@ namespace tomoto
 		template<ParallelScheme _ps, typename _ExtraDocData>
 		void mergeState(ThreadPool& pool, _ModelState& globalState, _ModelState& tState, _ModelState* localData, _RandGen*, const _ExtraDocData& edd) const
 		{
-			std::vector<std::future<void>> res;
-
 			if (_ps == ParallelScheme::copy_merge)
 			{
 				tState = globalState;
@@ -517,10 +540,27 @@ namespace tomoto
 				// make all count being positive
 				if (_tw != TermWeight::one)
 				{
-					globalState.numByTopicWord = globalState.numByTopicWord.cwiseMax(0);
+					globalState.numByTopicWord.matrix() = globalState.numByTopicWord.cwiseMax(0);
+				}
+				globalState.numByTopic = globalState.numByTopicWord.rowwise().sum();
+			}
+			else if (_ps == ParallelScheme::partition)
+			{
+				// make all count being positive
+				if (_tw != TermWeight::one)
+				{
+					globalState.numByTopicWord.matrix() = globalState.numByTopicWord.cwiseMax(0);
 				}
 				globalState.numByTopic = globalState.numByTopicWord.rowwise().sum();
+			}
+		}
 
+		template<ParallelScheme _ps>
+		void distributeMergedState(ThreadPool& pool, _ModelState& globalState, _ModelState* localData) const
+		{
+			std::vector<std::future<void>> res;
+			if (_ps == ParallelScheme::copy_merge)
+			{
 				for (size_t i = 0; i < pool.getNumWorkers(); ++i)
 				{
 					res.emplace_back(pool.enqueue([&, i](size_t)
@@ -531,22 +571,6 @@ namespace tomoto
 			}
 			else if (_ps == ParallelScheme::partition)
 			{
-				res = pool.enqueueToAll([&](size_t partitionId)
-				{
-					size_t b = partitionId ? edd.vChunkOffset[partitionId - 1] : 0,
-						e = edd.vChunkOffset[partitionId];
-					globalState.numByTopicWord.block(0, b, globalState.numByTopicWord.rows(), e - b) = localData[partitionId].numByTopicWord;
-				});
-				for (auto& r : res) r.get();
-				res.clear();
-
-				// make all count being positive
-				if (_tw != TermWeight::one)
-				{
-					globalState.numByTopicWord = globalState.numByTopicWord.cwiseMax(0);
-				}
-				globalState.numByTopic = globalState.numByTopicWord.rowwise().sum();
-
 				res = pool.enqueueToAll([&](size_t threadId)
 				{
 					localData[threadId].numByTopic = globalState.numByTopic;
@@ -560,16 +584,11 @@ namespace tomoto
 		ex) document pathing at hLDA model
 		* if pool is nullptr, workers has been already pooled and cannot branch works more.
 		*/
-		template<typename _DocIter>
+		template<GlobalSampler _gs, typename _DocIter>
 		void sampleGlobalLevel(ThreadPool* pool, _ModelState* localData, _RandGen* rgs, _DocIter first, _DocIter last) const
 		{
 		}
 
-		template<typename _DocIter>
-		void sampleGlobalLevel(ThreadPool* pool, _ModelState* localData, _RandGen* rgs, _DocIter first, _DocIter last)
-		{
-		}
-
 		template<typename _DocIter>
 		double getLLDocs(_DocIter _first, _DocIter _last) const
 		{
@@ -592,16 +611,33 @@ namespace tomoto
 			double ll = 0;
 			const size_t V = this->realV;
 			// topic-word distribution
-			auto lgammaEta = math::lgammaT(eta);
-			ll += math::lgammaT(V*eta) * K;
-			for (Tid k = 0; k < K; ++k)
+			if (etaByTopicWord.size())
 			{
-				ll -= math::lgammaT(ld.numByTopic[k] + V * eta);
-				for (Vid v = 0; v < V; ++v)
+				for (Tid k = 0; k < K; ++k)
 				{
-					if (!ld.numByTopicWord(k, v)) continue;
-					ll += math::lgammaT(ld.numByTopicWord(k, v) + eta) - lgammaEta;
-					assert(std::isfinite(ll));
+					Float etasum = etaByTopicWord.row(k).sum();
+					ll += math::lgammaT(etasum) - math::lgammaT(ld.numByTopic[k] + etasum);
+					for (Vid v = 0; v < V; ++v)
+					{
+						if (!ld.numByTopicWord(k, v)) continue;
+						ll += math::lgammaT(ld.numByTopicWord(k, v) + etaByTopicWord(v, k)) - math::lgammaT(etaByTopicWord(v, k));
+						assert(std::isfinite(ll));
+					}
+				}
+			}
+			else
+			{
+				auto lgammaEta = math::lgammaT(eta);
+				ll += math::lgammaT(V * eta) * K;
+				for (Tid k = 0; k < K; ++k)
+				{
+					ll -= math::lgammaT(ld.numByTopic[k] + V * eta);
+					for (Vid v = 0; v < V; ++v)
+					{
+						if (!ld.numByTopicWord(k, v)) continue;
+						ll += math::lgammaT(ld.numByTopicWord(k, v) + eta) - lgammaEta;
+						assert(std::isfinite(ll));
+					}
 				}
 			}
 			return ll;
@@ -637,9 +673,9 @@ namespace tomoto
 		void prepareDoc(_DocType& doc, size_t docId, size_t wordSize) const
 		{
 			sortAndWriteOrder(doc.words, doc.wOrder);
-			doc.numByTopic.init(getTopicDocPtr(docId), K);
+			doc.numByTopic.init(getTopicDocPtr(docId), K, 1);
 			doc.Zs = tvector<Tid>(wordSize);
-			if(_tw != TermWeight::one) doc.wordWeights.resize(wordSize, 1);
+			if(_tw != TermWeight::one) doc.wordWeights.resize(wordSize);
 		}
 
 		void prepareWordPriors()
@@ -664,7 +700,8 @@ namespace tomoto
 			if (initDocs)
 			{
 				this->globalState.numByTopic = Eigen::Matrix<WeightType, -1, 1>::Zero(K);
-				this->globalState.numByTopicWord = Eigen::Matrix<WeightType, -1, -1>::Zero(K, V);
+				//this->globalState.numByTopicWord = Eigen::Matrix<WeightType, -1, -1>::Zero(K, V);
+				this->globalState.numByTopicWord.init(nullptr, K, V);
 			}
 			if(m_flags & flags::continuous_doc_data) numByTopicDoc = Eigen::Matrix<WeightType, -1, -1>::Zero(K, this->docs.size());
 		}
@@ -791,12 +828,18 @@ namespace tomoto
 				for (size_t i = 0; i < maxIter; ++i)
 				{
 					std::vector<std::future<void>> res;
-					performSampling<_ps, true>(pool,
+					static_cast<const DerivedClass*>(this)->template performSampling<_ps, true>(pool,
 						(m_flags & flags::shared_state) ? &tmpState : localData.data(), rgs.data(), res,
-						docFirst, docLast, edd);
+						docFirst, docLast, edd
+					);
 					static_cast<const DerivedClass*>(this)->template mergeState<_ps>(pool, tmpState, tState, localData.data(), rgs.data(), edd);
-					static_cast<const DerivedClass*>(this)->template sampleGlobalLevel<>(
-						&pool, (m_flags & flags::shared_state) ? &tmpState : localData.data(), rgs.data(), docFirst, docLast);
+					static_cast<const DerivedClass*>(this)->template performSamplingGlobal<_ps, true>(&pool, tmpState, rgs.data(),
+						docFirst, docLast
+					);
+					static_cast<const DerivedClass*>(this)->template sampleGlobalLevel<GlobalSampler::inference>(
+						&pool, (m_flags & flags::shared_state) ? &tmpState : localData.data(), rgs.data(), docFirst, docLast
+					);
+					static_cast<const DerivedClass*>(this)->template distributeMergedState<_ps>(pool, tmpState, localData.data());
 				}
 				double ll = static_cast<const DerivedClass*>(this)->getLLRest(tmpState) - static_cast<const DerivedClass*>(this)->getLLRest(this->globalState);
 				ll += static_cast<const DerivedClass*>(this)->template getLLDocs<>(docFirst, docLast);
@@ -817,7 +860,9 @@ namespace tomoto
 					{
 						static_cast<const DerivedClass*>(this)->presampleDocument(*d, -1, tmpState, rgc, i);
 						static_cast<const DerivedClass*>(this)->template sampleDocument<ParallelScheme::none, true>(*d, edd, -1, tmpState, rgc, i);
-						static_cast<const DerivedClass*>(this)->template sampleGlobalLevel<>(
+						static_cast<const DerivedClass*>(this)->template performSamplingGlobal<_ps, true>(&pool, tmpState, &rgc,
+							&*d, &*d + 1);
+						static_cast<const DerivedClass*>(this)->template sampleGlobalLevel<GlobalSampler::inference>(
 							&pool, &tmpState, &rgc, &*d, &*d + 1);
 					}
 					double ll = static_cast<const DerivedClass*>(this)->getLLRest(tmpState) - gllRest;
@@ -845,7 +890,9 @@ namespace tomoto
 							static_cast<const DerivedClass*>(this)->template sampleDocument<ParallelScheme::none, true>(
 								*d, edd, -1, tmpState, rgc, i
 							);
-							static_cast<const DerivedClass*>(this)->template sampleGlobalLevel<>(
+							static_cast<const DerivedClass*>(this)->template performSamplingGlobal<_ps, true>(nullptr, tmpState, &rgc,
+								&*d, &*d + 1);
+							static_cast<const DerivedClass*>(this)->template sampleGlobalLevel<GlobalSampler::inference>(
 								nullptr, &tmpState, &rgc, &*d, &*d + 1
 							);
 						}
@@ -1036,7 +1083,7 @@ namespace tomoto
 	template<typename _TopicModel>
 	void DocumentLDA<_tw>::update(WeightType* ptr, const _TopicModel& mdl)
 	{
-		numByTopic.init(ptr, mdl.getK());
+		numByTopic.init(ptr, mdl.getK(), 1);
 		for (size_t i = 0; i < Zs.size(); ++i)
 		{
 			if (this->words[i] >= mdl.getV()) continue;
diff --git a/src/TopicModel/MGLDAModel.hpp b/src/TopicModel/MGLDAModel.hpp
index 414e965..f20d36b 100644
--- a/src/TopicModel/MGLDAModel.hpp
+++ b/src/TopicModel/MGLDAModel.hpp
@@ -107,18 +107,16 @@ namespace tomoto
 				e = edd.chunkOffsetByDoc(partitionId + 1, docId);
 			}
 
-			size_t vOffset = (_ps == ParallelScheme::partition && partitionId) ? edd.vChunkOffset[partitionId - 1] : 0;
-
 			const auto K = this->K;
 			for (size_t w = b; w < e; ++w)
 			{
 				if (doc.words[w] >= this->realV) continue;
-				addWordTo<-1>(ld, doc, w, doc.words[w] - vOffset, doc.Zs[w] - (doc.Zs[w] < K ? 0 : K), doc.sents[w], doc.Vs[w], doc.Zs[w] < K ? 0 : 1);
-				auto dist = getVZLikelihoods(ld, doc, doc.words[w] - vOffset, doc.sents[w]);
+				addWordTo<-1>(ld, doc, w, doc.words[w], doc.Zs[w] - (doc.Zs[w] < K ? 0 : K), doc.sents[w], doc.Vs[w], doc.Zs[w] < K ? 0 : 1);
+				auto dist = getVZLikelihoods(ld, doc, doc.words[w], doc.sents[w]);
 				auto vz = sample::sampleFromDiscreteAcc(dist, dist + T * (K + KL), rgs);
 				doc.Vs[w] = vz / (K + KL);
 				doc.Zs[w] = vz % (K + KL);
-				addWordTo<1>(ld, doc, w, doc.words[w] - vOffset, doc.Zs[w] - (doc.Zs[w] < K ? 0 : K), doc.sents[w], doc.Vs[w], doc.Zs[w] < K ? 0 : 1);
+				addWordTo<1>(ld, doc, w, doc.words[w], doc.Zs[w] - (doc.Zs[w] < K ? 0 : K), doc.sents[w], doc.Vs[w], doc.Zs[w] < K ? 0 : 1);
 			}
 		}
 
@@ -294,7 +292,7 @@ namespace tomoto
 			doc.Zs = tvector<Tid>(wordSize);
 			doc.Vs.resize(wordSize);
 			if (_tw != TermWeight::one) doc.wordWeights.resize(wordSize);
-			doc.numByTopic.init(nullptr, this->K + KL);
+			doc.numByTopic.init(nullptr, this->K + KL, 1);
 			doc.numBySentWin = Eigen::Matrix<WeightType, -1, -1>::Zero(S, T);
 			doc.numByWin = Eigen::Matrix<WeightType, -1, 1>::Zero(S + T - 1);
 			doc.numByWinL = Eigen::Matrix<WeightType, -1, 1>::Zero(S + T - 1);
@@ -308,7 +306,8 @@ namespace tomoto
 			if (initDocs)
 			{
 				this->globalState.numByTopic = Eigen::Matrix<WeightType, -1, 1>::Zero(this->K + KL);
-				this->globalState.numByTopicWord = Eigen::Matrix<WeightType, -1, -1>::Zero(this->K + KL, V);
+				//this->globalState.numByTopicWord = Eigen::Matrix<WeightType, -1, -1>::Zero(this->K + KL, V);
+				this->globalState.numByTopicWord.init(nullptr, this->K + KL, V);
 			}
 		}
 
@@ -533,7 +532,7 @@ namespace tomoto
 	template<typename _TopicModel>
 	void DocumentMGLDA<_tw>::update(WeightType * ptr, const _TopicModel & mdl)
 	{
-		this->numByTopic.init(ptr, mdl.getK() + mdl.getKL());
+		this->numByTopic.init(ptr, mdl.getK() + mdl.getKL(), 1);
 		numBySent.resize(*std::max_element(sents.begin(), sents.end()) + 1);
 		for (size_t i = 0; i < this->Zs.size(); ++i)
 		{
diff --git a/src/TopicModel/PAModel.hpp b/src/TopicModel/PAModel.hpp
index d29d3d8..ec8ff6a 100644
--- a/src/TopicModel/PAModel.hpp
+++ b/src/TopicModel/PAModel.hpp
@@ -144,7 +144,7 @@ namespace tomoto
 				size_t b = partitionId ? edd.vChunkOffset[partitionId - 1] : 0,
 					e = edd.vChunkOffset[partitionId];
 
-				localData[partitionId].numByTopicWord = globalState.numByTopicWord.block(0, b, globalState.numByTopicWord.rows(), e - b);
+				localData[partitionId].numByTopicWord.matrix() = globalState.numByTopicWord.block(0, b, globalState.numByTopicWord.rows(), e - b);
 				localData[partitionId].numByTopic = globalState.numByTopic;
 				localData[partitionId].numByTopic1_2 = globalState.numByTopic1_2;
 				localData[partitionId].numByTopic2 = globalState.numByTopic2;
@@ -157,8 +157,6 @@ namespace tomoto
 		template<ParallelScheme _ps, typename _ExtraDocData>
 		void mergeState(ThreadPool& pool, _ModelState& globalState, _ModelState& tState, _ModelState* localData, _RandGen*, const _ExtraDocData& edd) const
 		{
-			std::vector<std::future<void>> res;
-
 			if (_ps == ParallelScheme::copy_merge)
 			{
 				tState = globalState;
@@ -177,19 +175,12 @@ namespace tomoto
 					globalState.numByTopic = globalState.numByTopic.cwiseMax(0);
 					globalState.numByTopic1_2 = globalState.numByTopic1_2.cwiseMax(0);
 					globalState.numByTopic2 = globalState.numByTopic2.cwiseMax(0);
-					globalState.numByTopicWord = globalState.numByTopicWord.cwiseMax(0);
-				}
-
-				for (size_t i = 0; i < pool.getNumWorkers(); ++i)
-				{
-					res.emplace_back(pool.enqueue([&, this, i](size_t threadId)
-					{
-						localData[i] = globalState;
-					}));
+					globalState.numByTopicWord.matrix() = globalState.numByTopicWord.cwiseMax(0);
 				}
 			}
 			else if (_ps == ParallelScheme::partition)
 			{
+				std::vector<std::future<void>> res;
 				res = pool.enqueueToAll([&](size_t partitionId)
 				{
 					size_t b = partitionId ? edd.vChunkOffset[partitionId - 1] : 0,
@@ -197,7 +188,6 @@ namespace tomoto
 					globalState.numByTopicWord.block(0, b, globalState.numByTopicWord.rows(), e - b) = localData[partitionId].numByTopicWord;
 				});
 				for (auto& r : res) r.get();
-				res.clear();
 
 				tState.numByTopic1_2 = globalState.numByTopic1_2;
 				globalState.numByTopic1_2 = localData[0].numByTopic1_2;
@@ -209,11 +199,31 @@ namespace tomoto
 				// make all count being positive
 				if (_tw != TermWeight::one)
 				{
-					globalState.numByTopicWord = globalState.numByTopicWord.cwiseMax(0);
+					globalState.numByTopicWord.matrix() = globalState.numByTopicWord.cwiseMax(0);
 				}
 				globalState.numByTopic = globalState.numByTopic1_2.rowwise().sum();
 				globalState.numByTopic2 = globalState.numByTopicWord.rowwise().sum();
 
+			}
+		}
+
+
+		template<ParallelScheme _ps>
+		void distributeMergedState(ThreadPool& pool, _ModelState& globalState, _ModelState* localData) const
+		{
+			std::vector<std::future<void>> res;
+			if (_ps == ParallelScheme::copy_merge)
+			{
+				for (size_t i = 0; i < pool.getNumWorkers(); ++i)
+				{
+					res.emplace_back(pool.enqueue([&, i](size_t)
+					{
+						localData[i] = globalState;
+					}));
+				}
+			}
+			else if (_ps == ParallelScheme::partition)
+			{
 				res = pool.enqueueToAll([&](size_t threadId)
 				{
 					localData[threadId].numByTopic = globalState.numByTopic;
@@ -221,7 +231,6 @@ namespace tomoto
 					localData[threadId].numByTopic2 = globalState.numByTopic2;
 				});
 			}
-
 			for (auto& r : res) r.get();
 		}
 
@@ -304,7 +313,8 @@ namespace tomoto
 				this->globalState.numByTopic = Eigen::Matrix<WeightType, -1, 1>::Zero(this->K);
 				this->globalState.numByTopic2 = Eigen::Matrix<WeightType, -1, 1>::Zero(K2);
 				this->globalState.numByTopic1_2 = Eigen::Matrix<WeightType, -1, -1>::Zero(this->K, K2);
-				this->globalState.numByTopicWord = Eigen::Matrix<WeightType, -1, -1>::Zero(K2, V);
+				//this->globalState.numByTopicWord = Eigen::Matrix<WeightType, -1, -1>::Zero(K2, V);
+				this->globalState.numByTopicWord.init(nullptr, K2, V);
 			}
 		}
 
diff --git a/src/TopicModel/PT.h b/src/TopicModel/PT.h
new file mode 100644
index 0000000..cee1fd0
--- /dev/null
+++ b/src/TopicModel/PT.h
@@ -0,0 +1,27 @@
+#pragma once
+#include "LDA.h"
+
+namespace tomoto
+{
+	template<TermWeight _tw>
+	struct DocumentPTM : public DocumentLDA<_tw>
+	{
+		using BaseDocument = DocumentLDA<_tw>;
+		using DocumentLDA<_tw>::DocumentLDA;
+		using WeightType = typename DocumentLDA<_tw>::WeightType;
+		
+		uint64_t pseudoDoc = 0;
+
+		DEFINE_SERIALIZER_AFTER_BASE_WITH_VERSION(BaseDocument, 0, pseudoDoc);
+		DEFINE_TAGGED_SERIALIZER_AFTER_BASE_WITH_VERSION(BaseDocument, 1, 0x00010001, pseudoDoc);
+	};
+
+	class IPTModel : public ILDAModel
+	{
+	public:
+		using DefaultDocType = DocumentPTM<TermWeight::one>;
+		static IPTModel* create(TermWeight _weight, size_t _K = 1, size_t _P = 100, 
+			Float alpha = 0.1, Float eta = 0.01, Float lambda = 0.01, size_t seed = std::random_device{}(),
+			bool scalarRng = false);
+	};
+}
\ No newline at end of file
diff --git a/src/TopicModel/PTModel.cpp b/src/TopicModel/PTModel.cpp
new file mode 100644
index 0000000..1d62f77
--- /dev/null
+++ b/src/TopicModel/PTModel.cpp
@@ -0,0 +1,10 @@
+#include "PTModel.hpp"
+
+namespace tomoto
+{
+
+	IPTModel* IPTModel::create(TermWeight _weight, size_t _K, size_t _P, Float _alpha, Float _eta, Float _lambda, size_t seed, bool scalarRng)
+	{
+		TMT_SWITCH_TW(_weight, scalarRng, PTModel, _K, _P, _alpha, _eta, _lambda, seed);
+	}
+}
\ No newline at end of file
diff --git a/src/TopicModel/PTModel.hpp b/src/TopicModel/PTModel.hpp
new file mode 100644
index 0000000..669095c
--- /dev/null
+++ b/src/TopicModel/PTModel.hpp
@@ -0,0 +1,273 @@
+#pragma once
+#include "LDAModel.hpp"
+#include "PT.h"
+
+/*
+Implementation of Pseudo-document topic model using Gibbs sampling by bab2min
+
+Zuo, Y., Wu, J., Zhang, H., Lin, H., Wang, F., Xu, K., & Xiong, H. (2016, August). Topic modeling of short texts: A pseudo-document view. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining (pp. 2105-2114).
+*/
+
+namespace tomoto
+{
+	template<TermWeight _tw>
+	struct ModelStatePTM : public ModelStateLDA<_tw>
+	{
+		using WeightType = typename ModelStateLDA<_tw>::WeightType;
+
+		Eigen::Array<Float, -1, 1> pLikelihood;
+		Eigen::ArrayXi numDocsByPDoc;
+		Eigen::Matrix<WeightType, -1, -1> numByTopicPDoc;
+
+		//DEFINE_SERIALIZER_AFTER_BASE(ModelStateLDA<_tw>);
+	};
+
+	template<TermWeight _tw, typename _RandGen,
+		typename _Interface = IPTModel,
+		typename _Derived = void,
+		typename _DocType = DocumentPTM<_tw>,
+		typename _ModelState = ModelStatePTM<_tw>>
+	class PTModel : public LDAModel<_tw, _RandGen, flags::continuous_doc_data | flags::partitioned_multisampling, _Interface,
+		typename std::conditional<std::is_same<_Derived, void>::value, PTModel<_tw, _RandGen>, _Derived>::type,
+		_DocType, _ModelState>
+	{
+	protected:
+		using DerivedClass = typename std::conditional<std::is_same<_Derived, void>::value, PTModel<_tw, _RandGen>, _Derived>::type;
+		using BaseClass = LDAModel<_tw, _RandGen, flags::continuous_doc_data | flags::partitioned_multisampling, _Interface, DerivedClass, _DocType, _ModelState>;
+		friend BaseClass;
+		friend typename BaseClass::BaseClass;
+		using WeightType = typename BaseClass::WeightType;
+
+		static constexpr char TMID[] = "PTM";
+
+		uint64_t numPDocs;
+		Float lambda;
+		uint32_t pseudoDocSamplingInterval = 10;
+
+		void optimizeParameters(ThreadPool& pool, _ModelState* localData, _RandGen* rgs)
+		{
+			const auto K = this->K;
+			for (size_t i = 0; i < 10; ++i)
+			{
+				Float denom = this->calcDigammaSum(&pool, [&](size_t i) { return this->globalState.numByTopicPDoc.col(i).sum(); }, numPDocs, this->alphas.sum());
+				for (size_t k = 0; k < K; ++k)
+				{
+					Float nom = this->calcDigammaSum(&pool, [&](size_t i) { return this->globalState.numByTopicPDoc(k, i);}, numPDocs, this->alphas(k));
+					this->alphas(k) = std::max(nom / denom * this->alphas(k), 1e-5f);
+				}
+			}
+		}
+
+		void samplePseudoDoc(ThreadPool* pool, _ModelState& ld, _RandGen& rgs, _DocType& doc) const
+		{
+			if (doc.getSumWordWeight() == 0) return;
+			Eigen::Array<WeightType, -1, 1> docTopicDist = Eigen::Array<WeightType, -1, 1>::Zero(this->K);
+			for (size_t i = 0; i < doc.words.size(); ++i)
+			{
+				if (doc.words[i] >= this->realV) continue;
+				this->template addWordTo<-1>(ld, doc, i, doc.words[i], doc.Zs[i]);
+				typename std::conditional<_tw != TermWeight::one, float, int32_t>::type weight
+					= _tw != TermWeight::one ? doc.wordWeights[i] : 1;
+				docTopicDist[doc.Zs[i]] += weight;
+			}
+			--ld.numDocsByPDoc[doc.pseudoDoc];
+			
+			if (pool)
+			{
+				std::vector<std::future<void>> futures;
+				for (size_t w = 0; w < pool->getNumWorkers(); ++w)
+				{
+					futures.emplace_back(pool->enqueue([&](size_t, size_t w)
+					{
+						for (size_t p = w; p < numPDocs; p += pool->getNumWorkers())
+						{
+							Float ax = math::lgammaSubt(ld.numByTopicPDoc.col(p).array().template cast<Float>() + this->alphas.array(), docTopicDist.template cast<Float>()).sum();
+							Float bx = math::lgammaSubt(ld.numByTopicPDoc.col(p).sum() + this->alphas.sum(), docTopicDist.sum());
+							ld.pLikelihood[p] = ax - bx;
+						}
+					}, w));
+				}
+				for (auto& f : futures) f.get();
+			}
+			else
+			{
+				for (size_t p = 0; p < numPDocs; ++p)
+				{
+					Float ax = math::lgammaSubt(ld.numByTopicPDoc.col(p).array().template cast<Float>() + this->alphas.array(), docTopicDist.template cast<Float>()).sum();
+					Float bx = math::lgammaSubt(ld.numByTopicPDoc.col(p).sum() + this->alphas.sum(), docTopicDist.sum());
+					ld.pLikelihood[p] = ax - bx;
+				}
+			}
+			ld.pLikelihood = (ld.pLikelihood - ld.pLikelihood.maxCoeff()).exp();
+			ld.pLikelihood *= ld.numDocsByPDoc.template cast<Float>() + lambda;
+
+			sample::prefixSum(ld.pLikelihood.data(), numPDocs);
+			doc.pseudoDoc = sample::sampleFromDiscreteAcc(ld.pLikelihood.data(), ld.pLikelihood.data() + numPDocs, rgs);
+
+			++ld.numDocsByPDoc[doc.pseudoDoc];
+			doc.numByTopic.init(ld.numByTopicPDoc.col(doc.pseudoDoc).data(), this->K, 1);
+			for (size_t i = 0; i < doc.words.size(); ++i)
+			{
+				if (doc.words[i] >= this->realV) continue;
+				this->template addWordTo<1>(ld, doc, i, doc.words[i], doc.Zs[i]);
+			}
+		}
+
+		template<ParallelScheme _ps, bool _infer, typename _DocIter>
+		void performSamplingGlobal(ThreadPool* pool, _ModelState& globalState, _RandGen* rgs,
+			_DocIter docFirst, _DocIter docLast) const
+		{
+			if (this->globalStep % pseudoDocSamplingInterval) return;
+			for (; docFirst != docLast; ++docFirst)
+			{
+				samplePseudoDoc(pool, globalState, rgs[0], *docFirst);
+			}
+		}
+		
+		template<typename _DocIter>
+		double getLLDocs(_DocIter _first, _DocIter _last) const
+		{
+			double ll = 0;
+			// doc-topic distribution
+			for (; _first != _last; ++_first)
+			{
+				auto& doc = *_first;
+			}
+			return ll;
+		}
+
+		double getLLRest(const _ModelState& ld) const
+		{
+			double ll = BaseClass::getLLRest(ld);
+			const size_t V = this->realV;
+			ll -= math::lgammaT(ld.numDocsByPDoc.sum() + lambda * numPDocs) - math::lgammaT(lambda * numPDocs);
+			// pseudo_doc-topic distribution
+			for (size_t p = 0; p < numPDocs; ++p)
+			{
+				ll += math::lgammaT(ld.numDocsByPDoc[p] + lambda) - math::lgammaT(lambda);
+				ll -= math::lgammaT(ld.numByTopicPDoc.col(p).sum() + this->alphas.sum()) - math::lgammaT(this->alphas.sum());
+				for (Tid k = 0; k < this->K; ++k)
+				{
+					ll += math::lgammaT(ld.numByTopicPDoc(k, p) + this->alphas[k]) - math::lgammaT(this->alphas[k]);
+				}
+			}
+			return ll;
+		}
+
+		void prepareDoc(_DocType& doc, size_t docId, size_t wordSize) const
+		{
+			sortAndWriteOrder(doc.words, doc.wOrder);
+			doc.numByTopic.init((WeightType*)this->globalState.numByTopicPDoc.col(0).data(), this->K, 1);
+			doc.Zs = tvector<Tid>(wordSize);
+			if (_tw != TermWeight::one) doc.wordWeights.resize(wordSize);
+		}
+
+		void initGlobalState(bool initDocs)
+		{
+			this->alphas.resize(this->K);
+			this->alphas.array() = this->alpha;
+			this->globalState.pLikelihood = Eigen::Matrix<Float, -1, 1>::Zero(numPDocs);
+			this->globalState.numDocsByPDoc = Eigen::ArrayXi::Zero(numPDocs);
+			this->globalState.numByTopicPDoc = Eigen::Matrix<WeightType, -1, -1>::Zero(this->K, numPDocs);
+			BaseClass::initGlobalState(initDocs);
+		}
+
+		struct Generator
+		{
+			std::uniform_int_distribution<uint64_t> psi;
+			std::uniform_int_distribution<Tid> theta;
+		};
+
+		Generator makeGeneratorForInit(const _DocType*) const
+		{
+			return Generator{ 
+				std::uniform_int_distribution<uint64_t>{0, numPDocs - 1},
+				std::uniform_int_distribution<Tid>{0, (Tid)(this->K - 1)}
+			};
+		}
+
+		template<bool _Infer>
+		void updateStateWithDoc(Generator& g, _ModelState& ld, _RandGen& rgs, _DocType& doc, size_t i) const
+		{
+			if (i == 0)
+			{
+				doc.pseudoDoc = g.psi(rgs);
+				++ld.numDocsByPDoc[doc.pseudoDoc];
+				doc.numByTopic.init(ld.numByTopicPDoc.col(doc.pseudoDoc).data(), this->K, 1);
+			}
+			auto& z = doc.Zs[i];
+			auto w = doc.words[i];
+			if (this->etaByTopicWord.size())
+			{
+				auto col = this->etaByTopicWord.col(w);
+				z = sample::sampleFromDiscrete(col.data(), col.data() + col.size(), rgs);
+			}
+			else
+			{
+				z = g.theta(rgs);
+			}
+			this->template addWordTo<1>(ld, doc, i, w, z);
+		}
+
+		template<ParallelScheme _ps, bool _infer, typename _DocIter, typename _ExtraDocData>
+		void performSampling(ThreadPool& pool, _ModelState* localData, _RandGen* rgs, std::vector<std::future<void>>& res,
+			_DocIter docFirst, _DocIter docLast, const _ExtraDocData& edd) const
+		{
+			// single-threaded sampling
+			if (_ps == ParallelScheme::none)
+			{
+				forShuffled((size_t)std::distance(docFirst, docLast), rgs[0](), [&](size_t id)
+				{
+					static_cast<const DerivedClass*>(this)->presampleDocument(docFirst[id], id, *localData, *rgs, this->globalStep);
+					static_cast<const DerivedClass*>(this)->template sampleDocument<_ps, _infer>(
+						docFirst[id], edd, id,
+						*localData, *rgs, this->globalStep, 0);
+
+				});
+			}
+			// multi-threaded sampling on partition and update into global
+			else if (_ps == ParallelScheme::partition)
+			{
+				const size_t chStride = pool.getNumWorkers();
+				for (size_t i = 0; i < chStride; ++i)
+				{
+					res = pool.enqueueToAll([&, i, chStride](size_t partitionId)
+					{
+						forShuffled((size_t)std::distance(docFirst, docLast), rgs[partitionId](), [&](size_t id)
+						{
+							if ((docFirst[id].pseudoDoc + partitionId) % chStride != i) return;
+							static_cast<const DerivedClass*>(this)->template sampleDocument<_ps, _infer>(
+								docFirst[id], edd, id,
+								localData[partitionId], rgs[partitionId], this->globalStep, partitionId
+							);
+						});
+					});
+					for (auto& r : res) r.get();
+					res.clear();
+				}
+			}
+			else
+			{
+				throw std::runtime_error{ "Unsupported ParallelScheme" };
+			}
+		}
+
+	public:
+		DEFINE_SERIALIZER_AFTER_BASE_WITH_VERSION(BaseClass, 0, numPDocs, lambda);
+		DEFINE_TAGGED_SERIALIZER_AFTER_BASE_WITH_VERSION(BaseClass, 1, 0x00010001, numPDocs, lambda);
+
+		PTModel(size_t _K = 1, size_t _P = 100, Float _alpha = 1.0, Float _eta = 0.01, Float _lambda = 0.01,
+			size_t _rg = std::random_device{}())
+			: BaseClass(_K, _alpha, _eta, _rg), numPDocs(_P), lambda(_lambda)
+		{
+		}
+
+		void updateDocs()
+		{
+			for (auto& doc : this->docs)
+			{
+				doc.template update<>(this->getTopicDocPtr(doc.pseudoDoc), *static_cast<DerivedClass*>(this));
+			}
+		}
+	};
+}
diff --git a/src/TopicModel/TopicModel.hpp b/src/TopicModel/TopicModel.hpp
index 64182c8..719e3cd 100644
--- a/src/TopicModel/TopicModel.hpp
+++ b/src/TopicModel/TopicModel.hpp
@@ -121,6 +121,7 @@ namespace tomoto
 	};
 
 	enum class ParallelScheme { default_, none, copy_merge, partition, size };
+	enum class GlobalSampler { train, freeze_topics, inference, size };
 
 	inline const char* toString(ParallelScheme ps)
 	{
@@ -236,7 +237,7 @@ namespace tomoto
 		virtual const std::vector<uint64_t>& getVocabCf() const = 0;
 		virtual const std::vector<uint64_t>& getVocabDf() const = 0;
 
-		virtual int train(size_t iteration, size_t numWorkers, ParallelScheme ps = ParallelScheme::default_) = 0;
+		virtual int train(size_t iteration, size_t numWorkers, ParallelScheme ps = ParallelScheme::default_, bool freeze_topics = false) = 0;
 		virtual size_t getGlobalStep() const = 0;
 		virtual void prepare(bool initDocs = true, size_t minWordCnt = 0, size_t minWordDf = 0, size_t removeTopN = 0) = 0;
 		
@@ -588,7 +589,7 @@ namespace tomoto
 			return ps;
 		}
 
-		int train(size_t iteration, size_t numWorkers, ParallelScheme ps) override
+		int train(size_t iteration, size_t numWorkers, ParallelScheme ps, bool freeze_topics = false) override
 		{
 			if (!numWorkers) numWorkers = std::thread::hardware_concurrency();
 			ps = getRealScheme(ps);
@@ -606,16 +607,20 @@ namespace tomoto
 				localRG.emplace_back(rg());
 			}
 
-			for (size_t i = 0; i < numWorkers; ++i)
+			if (ps == ParallelScheme::copy_merge)
 			{
-				if(ps == ParallelScheme::copy_merge) localData.emplace_back(static_cast<_Derived*>(this)->globalState);
+				for (size_t i = 0; i < numWorkers; ++i)
+				{
+					localData.emplace_back(static_cast<_Derived*>(this)->globalState);
+				}
 			}
-
-			if (ps == ParallelScheme::partition)
+			else if (ps == ParallelScheme::partition)
 			{
 				localData.resize(numWorkers);
-				static_cast<_Derived*>(this)->updatePartition(*cachedPool, globalState, localData.data(), docs.begin(), docs.end(), 
-					static_cast<_Derived*>(this)->eddTrain);
+				static_cast<_Derived*>(this)->updatePartition(
+					*cachedPool, globalState, localData.data(), docs.begin(), docs.end(), 
+					static_cast<_Derived*>(this)->eddTrain
+				);
 			}
 
 			auto state = ps == ParallelScheme::none ? &globalState : localData.data();
@@ -629,15 +634,15 @@ namespace tomoto
 						{
 						case ParallelScheme::none:
 							static_cast<_Derived*>(this)->template trainOne<ParallelScheme::none>(
-								*cachedPool, state, localRG.data());
+								*cachedPool, state, localRG.data(), freeze_topics);
 							break;
 						case ParallelScheme::copy_merge:
 							static_cast<_Derived*>(this)->template trainOne<ParallelScheme::copy_merge>(
-								*cachedPool, state, localRG.data());
+								*cachedPool, state, localRG.data(), freeze_topics);
 							break;
 						case ParallelScheme::partition:
 							static_cast<_Derived*>(this)->template trainOne<ParallelScheme::partition>(
-								*cachedPool, state, localRG.data());
+								*cachedPool, state, localRG.data(), freeze_topics);
 							break;
 						}
 						break;
diff --git a/src/Utils/Trie.hpp b/src/Utils/Trie.hpp
index 6d8a757..df480e9 100644
--- a/src/Utils/Trie.hpp
+++ b/src/Utils/Trie.hpp
@@ -26,10 +26,31 @@ namespace tomoto
 		}
 	};
 
+	template<class _Map, class _Node>
+	class TrieIterator : public _Map::const_iterator
+	{
+		using Base = typename _Map::const_iterator;
+		using Key = typename _Map::key_type;
+		const _Node* base = nullptr;
+	public:
+
+		TrieIterator(const Base& it, const _Node* _base)
+			: Base(it), base(_base)
+		{
+		}
+
+		std::pair<const Key, const _Node*> operator*() const
+		{
+			auto p = Base::operator*();
+			return std::make_pair(p.first, base + p.second);
+		}
+	};
+
 	template<class _Key, class _Value, class _KeyStore = ConstAccess<std::map<_Key, int32_t>>, class _Trie = void>
 	struct Trie
 	{
 		using Node = typename std::conditional<std::is_same<_Trie, void>::value, Trie, _Trie>::type;
+		using iterator = TrieIterator<_KeyStore, Node>;
 		_KeyStore next = {};
 		int32_t fail = 0;
 		_Value val = {};
@@ -47,13 +68,23 @@ namespace tomoto
 			return fail ? (Node*)this + fail : nullptr;
 		}
 
+		iterator begin() const
+		{
+			return { next.begin(), (const Node*)this };
+		}
+
+		iterator end() const
+		{
+			return { next.end(), (const Node*)this };
+		}
+
 		template<typename _TyIter, typename _FnAlloc>
-		void build(_TyIter first, _TyIter last, const _Value& _val, _FnAlloc&& alloc)
+		Node* build(_TyIter first, _TyIter last, const _Value& _val, _FnAlloc&& alloc)
 		{
 			if (first == last)
 			{
 				if (!val) val = _val;
-				return;
+				return (Node*)this;
 			}
 
 			auto v = *first;
@@ -61,13 +92,13 @@ namespace tomoto
 			{
 				next[v] = alloc() - this;
 			}
-			getNext(v)->build(++first, last, _val, alloc);
+			return getNext(v)->build(++first, last, _val, alloc);
 		}
 
 		template<typename _TyIter>
 		Node* findNode(_TyIter begin, _TyIter end)
 		{
-			if (begin == end) return this;
+			if (begin == end) return (Node*)this;
 			auto n = getNext(*begin);
 			if (n) return n->findNode(++begin, end);
 			return nullptr;
@@ -173,21 +204,21 @@ namespace tomoto
 		int32_t parent = 0;
 
 		template<typename _TyIter, typename _FnAlloc>
-		void build(_TyIter first, _TyIter last, const _Value& _val, _FnAlloc&& alloc)
+		TrieEx* build(_TyIter first, _TyIter last, const _Value& _val, _FnAlloc&& alloc)
 		{
 			if (first == last)
 			{
 				if (!this->val) this->val = _val;
-				return;
+				return this;
 			}
 
 			auto v = *first;
-			if (!getNext(v))
+			if (!this->getNext(v))
 			{
 				this->next[v] = alloc() - this;
 				this->getNext(v)->parent = -this->next[v];
 			}
-			this->getNext(v)->build(++first, last, _val, alloc);
+			return this->getNext(v)->build(++first, last, _val, alloc);
 		}
 
 		template<typename _FnAlloc>
diff --git a/src/Utils/Utils.hpp b/src/Utils/Utils.hpp
index f7cdbc8..d3da7f2 100644
--- a/src/Utils/Utils.hpp
+++ b/src/Utils/Utils.hpp
@@ -70,7 +70,7 @@ namespace tomoto
 	}
 
 	template<class UnaryFunction>
-	UnaryFunction forRandom(size_t N, size_t seed, UnaryFunction f)
+	UnaryFunction forShuffled(size_t N, size_t seed, UnaryFunction f)
 	{
 		static size_t primes[16] = {
 			65537, 65539, 65543, 65551, 65557, 65563,
@@ -206,132 +206,137 @@ namespace tomoto
 	}
 
 	template <typename _UnaryFunc, typename _Iterator>
-	class TransformIter
+	class TransformIter : public _Iterator
 	{
 	private:
-		_Iterator i;
 		_UnaryFunc f;
 	public:
 		using reference = typename std::result_of<
 			const _UnaryFunc(typename std::iterator_traits<_Iterator>::reference)
 		>::type;
 		using value_type = reference;
-
-		using pointer = void;
-		using iterator_category = typename std::iterator_traits<_Iterator>::iterator_category;
-		using difference_type = typename std::iterator_traits<_Iterator>::difference_type;
-
+		
 		TransformIter(const _Iterator& _iter = {}, _UnaryFunc _f = {})
-			: i(_iter), f(_f)
+			: _Iterator(_iter), f(_f)
 		{}
 		
 		reference operator*()
 		{
-			return f(*i);
+			return f(_Iterator::operator*());
 		}
 
 		const reference operator*() const
 		{
-			return f(*i);
+			return f(_Iterator::operator*());
 		}
 
 		reference operator[](std::size_t idx)
 		{
-			return f(i[idx]);
+			return f(_Iterator::operator[](idx));
 		}
 
 		const reference operator[](std::size_t idx) const
 		{
-			return f(i[idx]);
+			return f(_Iterator::operator[](idx));
 		}
 
 		TransformIter& operator++()
 		{
-			++i;
+			_Iterator::operator++();
 			return *this;
 		}
 
-		TransformIter& operator++(int)
+		TransformIter operator++(int)
 		{
 			auto c = *this;
-			++i;
+			_Iterator::operator++();
 			return c;
 		}
 
 		TransformIter& operator--()
 		{
-			--i;
+			_Iterator::operator--();
 			return *this;
 		}
 
-		TransformIter& operator--(int)
+		TransformIter operator--(int)
 		{
 			auto c = *this;
-			--i;
+			_Iterator::operator--();
 			return c;
 		}
 
 		TransformIter operator+(int n) const
 		{
-			return { f, i + n };
+			return { _Iterator::operator+(n), f };
 		}
 
 		TransformIter operator-(int n) const
 		{
-			return { f, i - n };
+			return { _Iterator::operator-(n), f };
 		}
 
 		TransformIter& operator+=(int n)
 		{
-			i += n;
+			_Iterator::operator+=(n);
 			return *this;
 		}
 
 		TransformIter& operator-=(int n)
 		{
-			i -= n;
+			_Iterator::operator-=(n);
 			return *this;
 		}
 
 		typename std::iterator_traits<_Iterator>::difference_type operator-(const TransformIter& o) const
 		{
-			return i - o.i;
+			return (const _Iterator&)*this - (const _Iterator&)o;
 		}
 
-		bool operator==(const TransformIter& o) const
-		{
-			return i == o.i;
-		}
+	};
 
-		bool operator!=(const TransformIter& o) const
-		{
-			return i != o.i;
-		}
+	template <typename _UnaryFunc, typename _Iterator> 
+	TransformIter<_UnaryFunc, _Iterator> makeTransformIter(const _Iterator& iter, _UnaryFunc f)
+	{
+		return { iter, f };
+	}
 
-		bool operator<(const TransformIter& o) const
+	template <typename _Iterator>
+	class StrideIter : public _Iterator
+	{
+		size_t stride;
+		const _Iterator end;
+	public:
+		StrideIter(const _Iterator& iter, size_t _stride = 1, const _Iterator& _end = {})
+			: _Iterator{ iter }, stride{ _stride }, end{ _end }
 		{
-			return i < o.i;
 		}
 
-		bool operator>(const TransformIter& o) const
-		{
-			return i > o.i;
-		}
+		StrideIter(const StrideIter&) = default;
+		StrideIter(StrideIter&&) = default;
 
-		bool operator<=(const TransformIter& o) const
+		StrideIter& operator++()
 		{
-			return i <= o.i;
+			for (size_t i = 0; i < stride && *this != end; ++i)
+			{
+				_Iterator::operator++();
+			}
+			return *this;
 		}
 
-		bool operator>=(const TransformIter& o) const
+		StrideIter& operator--()
 		{
-			return i >= o.i;
+			for (size_t i = 0; i < stride && *this != end; ++i)
+			{
+				_Iterator::operator--();
+			}
+			return *this;
 		}
 	};
 
-	template <typename _UnaryFunc, typename _Iterator> 
-	TransformIter<_UnaryFunc, _Iterator> makeTransformIter(const _Iterator& iter, _UnaryFunc f)
+	template <typename _Iterator>
+	StrideIter<_Iterator> makeStrideIter(const _Iterator& iter, size_t stride, const _Iterator& end = {})
 	{
-		return { iter, f };
+		return { iter, stride, end };
 	}
 }
diff --git a/src/Utils/math.h b/src/Utils/math.h
index 27a9714..363e7ce 100644
--- a/src/Utils/math.h
+++ b/src/Utils/math.h
@@ -153,28 +153,32 @@ namespace tomoto
 		inline float lgammaT(float x) { return detail::LUT_lgamma::get(x); }
 		inline float digammaT(float x) { return detail::LUT_digamma::get(x); }
 
+		// approximation : lgamma(z) ~= (z+2.5)ln(z+3) - z - 3 + 0.5 ln (2pi) + 1/12/(z + 3) - ln (z(z+1)(z+2))
 		template<class _T>
 		inline auto lgammaApprox(_T z) -> decltype((z + 2.5)* log(z + 3) - (z + 3) + 0.91893853 + 1. / 12. / (z + 3) - log(z * (z + 1) * (z + 2)))
 		{
-			// approximation : lgamma(z) ~= (z+2.5)ln(z+3) - z - 3 + 0.5 ln (2pi) + 1/12/(z + 3) - ln (z(z+1)(z+2))
 			return (z + 2.5) * log(z + 3) - (z + 3) + 0.91893853 + 1. / 12. / (z + 3) - log(z * (z + 1) * (z + 2));
 		}
 
+		// calc lgamma(z + a) - lgamma(z)
 		template<class _T, class _U>
-		inline _T lgammaSubt(_T z, _U a) // calc lgamma(z + a) - lgamma(z)
+		inline auto lgammaSubt(_T z, _U a) -> decltype((z + a + 1.5)* log(z + a + 2) - (z + 1.5) * log(z + 2) - a + (1. / (z + a + 2) - 1. / (z + 2)) / 12. - log(((z + a) * (z + a + 1)) / (z * (z + 1))))
 		{
 			return (z + a + 1.5) * log(z + a + 2) - (z + 1.5) * log(z + 2) - a + (1. / (z + a + 2) - 1. / (z + 2)) / 12. - log(((z + a) * (z + a + 1)) / (z * (z + 1)));
 		}
 
+		// approximation : digamma(z) ~= ln(z+4) - 1/2/(z+4) - 1/12/(z+4)^2 - 1/z - 1/(z+1) - 1/(z+2) - 1/(z+3)
 		template<class _T>
 		inline auto digammaApprox(_T z) -> decltype(log(z + 4) - 1. / 2. / (z + 4) - 1. / 12. / ((z + 4) * (z + 4)) - 1. / z - 1. / (z + 1) - 1. / (z + 2) - 1. / (z + 3))
 		{
-			// approximation : digamma(z) ~= ln(z+4) - 1/2/(z+4) - 1/12/(z+4)^2 - 1/z - 1/(z+1) - 1/(z+2) - 1/(z+3)
 			return log(z + 4) - 1. / 2. / (z + 4) - 1. / 12. / ((z + 4) * (z + 4)) - 1. / z - 1. / (z + 1) - 1. / (z + 2) - 1. / (z + 3);
 		}
 
+		// calc digamma(z + a) - digamma(z)
 		template<class _T, class _U>
-		inline _T digammaSubt(_T z, _U a) // calc digamma(z + a) - digamma(z)
+		inline auto digammaSubt(_T z, _U a) -> decltype(log((z + a + 2) / (z + 2)) - (1 / (z + a + 2) - 1 / (z + 2)) / 2 - (1 / (z + a + 2) / (z + a + 2) - 1 / (z + 2) / (z + 2)) / 12
+			- 1. / (z + a) - 1. / (z + a + 1)
+			- 1. / z - 1. / (z + 1))
 		{
 			return log((z + a + 2) / (z + 2)) - (1 / (z + a + 2) - 1 / (z + 2)) / 2 - (1 / (z + a + 2) / (z + a + 2) - 1 / (z + 2) / (z + 2)) / 12
 				- 1. / (z + a) - 1. / (z + a + 1)
diff --git a/src/Utils/tvector.hpp b/src/Utils/tvector.hpp
index 2e59b50..52db212 100644
--- a/src/Utils/tvector.hpp
+++ b/src/Utils/tvector.hpp
@@ -172,10 +172,14 @@ namespace tomoto
 
 		reverse_iterator rbegin() noexcept { return reverse_iterator(_last); }
 
+		const_reverse_iterator rbegin() const noexcept { return reverse_iterator(_last); }
+
 		const_reverse_iterator crbegin() const noexcept { return reverse_iterator(_last); }
 
 		reverse_iterator rend() noexcept { return reverse_iterator(_first); }
 
+		const_reverse_iterator rend() const noexcept { return reverse_iterator(_first); }
+
 		const_reverse_iterator crend() const noexcept { return reverse_iterator(_first); }
 
 		// 23.3.11.3, capacity:
diff --git a/src/python/PyUtils.h b/src/python/PyUtils.h
index 53c9253..36a6a96 100644
--- a/src/python/PyUtils.h
+++ b/src/python/PyUtils.h
@@ -408,12 +408,12 @@ namespace py
 		};
 	}
 
-	struct cast_to_signed_t
-	{
-	};
-
+	struct cast_to_signed_t {};
 	static constexpr cast_to_signed_t cast_to_signed{};
 
+	struct force_list_t {};
+	static constexpr force_list_t force_list{};
+
 	template<typename _Ty>
 	struct ValueBuilder<std::vector<_Ty>, 
 		typename std::enable_if<std::is_arithmetic<_Ty>::value>::type>
@@ -542,6 +542,17 @@ namespace py
 		return ret;
 	}
 
+	template<typename _Ty>
+	inline PyObject* buildPyValue(const std::vector<_Ty>& v, force_list_t)
+	{
+		PyObject* ret = PyList_New(v.size());
+		for (size_t i = 0; i < v.size(); ++i)
+		{
+			PyList_SetItem(ret, i, buildPyValue(v[i]));
+		}
+		return ret;
+	}
+
 	template<typename _Ty, typename _Tx>
 	inline typename std::enable_if<
 		!std::is_arithmetic<
@@ -725,6 +736,12 @@ namespace py
 			}
 		}
 	};
+
+	class wrong_type_error : public std::runtime_error
+	{
+	public:
+		using std::runtime_error::runtime_error;
+	};
 }
 
 #define PRINT_WARN(msg) do{ py::WarningLog::get().print(std::cerr, msg); } while(0)
diff --git a/src/python/docs.h b/src/python/docs.h
index 1410749..d3848eb 100644
--- a/src/python/docs.h
+++ b/src/python/docs.h
@@ -356,7 +356,7 @@ word : str
 )"");
 
 DOC_SIGNATURE_EN_KO(LDA_train__doc__,
-	"train(self, iter=10, workers=0, parallel=0)",
+	"train(self, iter=10, workers=0, parallel=0, freeze_topics=False)",
 	u8R""(Train the model using Gibbs-sampling with `iter` iterations. Return `None`. 
 After calling this method, you cannot `tomotopy.LDAModel.add_doc` or `tomotopy.LDAModel.set_word_prior` more.
 
@@ -370,7 +370,11 @@ workers : int
 parallel : Union[int, tomotopy.ParallelScheme]
     .. versionadded:: 0.5.0
     
-    the parallelism scheme for training. the default value is ParallelScheme.DEFAULT which means that tomotopy selects the best scheme by model.
+    the parallelism scheme for training. the default value is `tomotopy.ParallelScheme.DEFAULT` which means that tomotopy selects the best scheme by model.
+freeze_topics : bool
+    .. versionadded:: 0.10.1
+
+    prevents to create a new topic when training. Only valid for `tomotopy.HLDAModel`
 )"",
 u8R""(깁스 샘플링을 `iter` 회 반복하여 현재 모델을 학습시킵니다. 반환값은 `None`입니다. 
 이 메소드가 호출된 이후에는 더 이상 `tomotopy.LDAModel.add_doc`로 현재 모델에 새로운 학습 문헌을 추가시킬 수 없습니다.
@@ -386,6 +390,10 @@ parallel : Union[int, tomotopy.ParallelScheme]
     .. versionadded:: 0.5.0
 
     학습에 사용할 병렬화 방법. 기본값은 ParallelScheme.DEFAULT로 이는 모델에 따라 최적의 방법을 tomotopy가 알아서 선택하도록 합니다.
+freeze_topics : bool
+    .. versionadded:: 0.10.1
+
+    학습 시 새로운 토픽을 생성하지 못하도록 합니다. 이 파라미터는 오직 `tomotopy.HLDAModel`에만 유효합니다.
 )"");
 
 DOC_SIGNATURE_EN_KO(LDA_get_topic_words__doc__,
diff --git a/src/python/label_docs.h b/src/python/label_docs.h
index ac4ff84..6063c44 100644
--- a/src/python/label_docs.h
+++ b/src/python/label_docs.h
@@ -31,6 +31,14 @@ DOC_VARIABLE_EN_KO(Candidate_name__doc__,
 	u8R""(an actual name of the candidate for topic label)"",
 	u8R""(토픽 레이블의 실제 이름)"");
 
+DOC_VARIABLE_EN_KO(Candidate_cf__doc__,
+    u8R""(collection frequency of the candidate (read-only))"",
+    u8R""(후보의 장서빈도(읽기전용))"");
+
+DOC_VARIABLE_EN_KO(Candidate_df__doc__,
+    u8R""(document frequency of the candidate (read-only))"",
+    u8R""(후보의 문헌빈도(읽기전용))"");
+
 /*
 	class PMIExtractor
 */
diff --git a/src/python/py_LDA.cpp b/src/python/py_LDA.cpp
index 7e5be56..b5e7b6e 100644
--- a/src/python/py_LDA.cpp
+++ b/src/python/py_LDA.cpp
@@ -178,9 +178,9 @@ PyObject* LDA_getWordPrior(TopicModelObject* self, PyObject* args, PyObject *kwa
 
 static PyObject* LDA_train(TopicModelObject* self, PyObject* args, PyObject *kwargs)
 {
-	size_t iteration = 10, workers = 0, ps = 0;
-	static const char* kwlist[] = { "iter", "workers", "parallel", nullptr };
-	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|nnn", (char**)kwlist, &iteration, &workers, &ps)) return nullptr;
+	size_t iteration = 10, workers = 0, ps = 0, fixed = 0;
+	static const char* kwlist[] = { "iter", "workers", "parallel", "freeze_topics", nullptr };
+	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|nnnp", (char**)kwlist, &iteration, &workers, &ps, &fixed)) return nullptr;
 	try
 	{
 		if (!self->inst) throw runtime_error{ "inst is null" };
@@ -190,7 +190,7 @@ static PyObject* LDA_train(TopicModelObject* self, PyObject* args, PyObject *kwa
 			inst->prepare(true, self->minWordCnt, self->minWordDf, self->removeTopWord);
 			self->isPrepared = true;
 		}
-		inst->train(iteration, workers, (tomoto::ParallelScheme)ps);
+		inst->train(iteration, workers, (tomoto::ParallelScheme)ps, !!fixed);
 		Py_INCREF(Py_None);
 		return Py_None;
 	}
@@ -282,6 +282,22 @@ PyObject* LDA_infer(TopicModelObject* self, PyObject* args, PyObject *kwargs)
 			auto ll = self->inst->infer(docs, iteration, tolerance, workers, (tomoto::ParallelScheme)ps, !!together);
 			return py::buildPyTuple(py::UniqueObj{ (PyObject*)cps }, ll);
 		}
+		else if (PyObject_TypeCheck(argDoc, &UtilsDocument_type))
+		{
+			auto* doc = (DocumentObject*)argDoc;
+			if (doc->corpus->tm != self) throw runtime_error{ "`doc` was from another model, not fit to this model" };
+			if (doc->owner)
+			{
+				std::vector<tomoto::DocumentBase*> docs;
+				docs.emplace_back((tomoto::DocumentBase*)doc->getBoundDoc());
+				float ll = self->inst->infer(docs, iteration, tolerance, workers, (tomoto::ParallelScheme)ps, !!together)[0];
+				return Py_BuildValue("(Nf)", py::buildPyValue(self->inst->getTopicsByDoc(doc->getBoundDoc())), ll);
+			}
+			else
+			{
+				return Py_BuildValue("(Ns)", py::buildPyValue(self->inst->getTopicsByDoc(doc->getBoundDoc())), nullptr);
+			}
+		}
 		else if ((iter = py::UniqueObj{ PyObject_GetIter(argDoc) }) != nullptr)
 		{
 			std::vector<tomoto::DocumentBase*> docs;
@@ -312,21 +328,7 @@ PyObject* LDA_infer(TopicModelObject* self, PyObject* args, PyObject *kwargs)
 		}
 		else
 		{
-			PyErr_Clear();
-			if (!PyObject_TypeCheck(argDoc, &UtilsDocument_type)) throw runtime_error{ "`doc` must be tomotopy.Document type or list of tomotopy.Document" };
-			auto* doc = (DocumentObject*)argDoc;
-			if (doc->corpus->tm != self) throw runtime_error{ "`doc` was from another model, not fit to this model" };
-			if (doc->owner)
-			{
-				std::vector<tomoto::DocumentBase*> docs;
-				docs.emplace_back((tomoto::DocumentBase*)doc->getBoundDoc());
-				float ll = self->inst->infer(docs, iteration, tolerance, workers, (tomoto::ParallelScheme)ps, !!together)[0];
-				return Py_BuildValue("(Nf)", py::buildPyValue(self->inst->getTopicsByDoc(doc->getBoundDoc())), ll);
-			}
-			else
-			{
-				return Py_BuildValue("(Ns)", py::buildPyValue(self->inst->getTopicsByDoc(doc->getBoundDoc())), nullptr);
-			}
+			throw runtime_error{ "`doc` must be tomotopy.Document type or list of tomotopy.Document" };	
 		}
 	}
 	catch (const bad_exception&)
diff --git a/src/python/py_PA.cpp b/src/python/py_PA.cpp
index 77caa30..5987039 100644
--- a/src/python/py_PA.cpp
+++ b/src/python/py_PA.cpp
@@ -222,6 +222,24 @@ static PyObject* PA_infer(TopicModelObject* self, PyObject* args, PyObject *kwar
 			auto ll = self->inst->infer(docs, iteration, tolerance, workers, (tomoto::ParallelScheme)ps, !!together);
 			return py::buildPyTuple(py::UniqueObj{ (PyObject*)cps }, ll);
 		}
+		else if (PyObject_TypeCheck(argDoc, &UtilsDocument_type))
+		{
+			auto* doc = (DocumentObject*)argDoc;
+			if (doc->corpus->tm != self) throw runtime_error{ "'doc' was from another model, not fit to this model" };
+			if (doc->owner)
+			{
+				std::vector<tomoto::DocumentBase*> docs;
+				docs.emplace_back((tomoto::DocumentBase*)doc->doc);
+				double ll = self->inst->infer(docs, iteration, tolerance, workers, (tomoto::ParallelScheme)ps, !!together)[0];
+				return Py_BuildValue("((NN)f)", py::buildPyValue(inst->getTopicsByDoc(doc->getBoundDoc())),
+					py::buildPyValue(inst->getSubTopicsByDoc(doc->getBoundDoc())), ll);
+			}
+			else
+			{
+				return Py_BuildValue("((NN)s)", py::buildPyValue(inst->getTopicsByDoc(doc->getBoundDoc())),
+					py::buildPyValue(inst->getSubTopicsByDoc(doc->getBoundDoc())), nullptr);
+			}
+		}
 		else if ((iter = py::UniqueObj{ PyObject_GetIter(argDoc) }) != nullptr)
 		{
 			std::vector<tomoto::DocumentBase*> docs;
@@ -256,23 +274,7 @@ static PyObject* PA_infer(TopicModelObject* self, PyObject* args, PyObject *kwar
 		}
 		else
 		{
-			PyErr_Clear();
-			if (!PyObject_TypeCheck(argDoc, &UtilsDocument_type)) throw runtime_error{ "'doc' must be tomotopy.Document type or list of tomotopy.Document" };
-			auto* doc = (DocumentObject*)argDoc;
-			if (doc->corpus->tm != self) throw runtime_error{ "'doc' was from another model, not fit to this model" };
-			if (doc->owner)
-			{
-				std::vector<tomoto::DocumentBase*> docs;
-				docs.emplace_back((tomoto::DocumentBase*)doc->doc);
-				double ll = self->inst->infer(docs, iteration, tolerance, workers, (tomoto::ParallelScheme)ps, !!together)[0];
-				return Py_BuildValue("((NN)f)", py::buildPyValue(inst->getTopicsByDoc(doc->getBoundDoc())), 
-					py::buildPyValue(inst->getSubTopicsByDoc(doc->getBoundDoc())), ll);
-			}
-			else
-			{
-				return Py_BuildValue("((NN)s)", py::buildPyValue(inst->getTopicsByDoc(doc->getBoundDoc())),
-					py::buildPyValue(inst->getSubTopicsByDoc(doc->getBoundDoc())), nullptr);
-			}
+			throw runtime_error{ "'doc' must be tomotopy.Document type or list of tomotopy.Document" };
 		}
 	}
 	catch (const bad_exception&)
diff --git a/src/python/py_SLDA.cpp b/src/python/py_SLDA.cpp
index 609349b..4c94c1f 100644
--- a/src/python/py_SLDA.cpp
+++ b/src/python/py_SLDA.cpp
@@ -243,33 +243,34 @@ static PyObject* SLDA_estimateVars(TopicModelObject* self, PyObject* args, PyObj
 	{
 		if (!self->inst) throw runtime_error{ "inst is null" };
 		auto* inst = static_cast<tomoto::ISLDAModel*>(self->inst);
-		if (py::UniqueObj iter = py::UniqueObj{ PyObject_GetIter(argDoc) })
+		try 
 		{
-			py::UniqueObj nextDoc;
-			std::vector<const tomoto::DocumentBase*> docs;
-			while ((nextDoc = py::UniqueObj{ PyIter_Next(iter) }))
-			{
-				if (!PyObject_TypeCheck(nextDoc, &UtilsDocument_type)) throw runtime_error{ "`doc` must be tomotopy.Document or list of tomotopy.Document" };
-				auto* doc = (DocumentObject*)nextDoc.get();
-				if (doc->corpus->tm != self) throw runtime_error{ "`doc` was from another model, not fit to this model" };
-				docs.emplace_back(doc->getBoundDoc());
-			}
-			if (PyErr_Occurred()) return nullptr;
-			return py::buildPyValueTransform(docs.begin(), docs.end(), [&](const tomoto::DocumentBase* d)
-			{
-				return inst->estimateVars(d);
-			});
+			if (!PyObject_TypeCheck(argDoc, &UtilsDocument_type)) throw py::wrong_type_error{ "`doc` must be tomotopy.Document or list of tomotopy.Document" };
+			auto* doc = (DocumentObject*)argDoc;
+			if (doc->corpus->tm != self) throw py::wrong_type_error{ "`doc` was from another model, not fit to this model" };
+
+			return py::buildPyValue(inst->estimateVars(doc->getBoundDoc()));
 		}
-		else
+		catch (const py::wrong_type_error&)
 		{
 			PyErr_Clear();
 		}
-
-		if (!PyObject_TypeCheck(argDoc, &UtilsDocument_type)) throw runtime_error{ "`doc` must be tomotopy.Document or list of tomotopy.Document" };
-		auto* doc = (DocumentObject*)argDoc;
-		if (doc->corpus->tm != self) throw runtime_error{ "`doc` was from another model, not fit to this model" };
-
-		return py::buildPyValue(inst->estimateVars(doc->getBoundDoc()));
+		
+		py::UniqueObj iter = py::UniqueObj{ PyObject_GetIter(argDoc) };
+		py::UniqueObj nextDoc;
+		std::vector<const tomoto::DocumentBase*> docs;
+		while ((nextDoc = py::UniqueObj{ PyIter_Next(iter) }))
+		{
+			if (!PyObject_TypeCheck(nextDoc, &UtilsDocument_type)) throw py::wrong_type_error{ "`doc` must be tomotopy.Document or list of tomotopy.Document" };
+			auto* doc = (DocumentObject*)nextDoc.get();
+			if (doc->corpus->tm != self) throw py::wrong_type_error{ "`doc` was from another model, not fit to this model" };
+			docs.emplace_back(doc->getBoundDoc());
+		}
+		if (PyErr_Occurred()) return nullptr;
+		return py::buildPyValueTransform(docs.begin(), docs.end(), [&](const tomoto::DocumentBase* d)
+		{
+			return inst->estimateVars(d);
+		});
 	}
 	catch (const bad_exception&)
 	{
diff --git a/src/python/py_label.cpp b/src/python/py_label.cpp
index d6751fe..d8f2ce5 100644
--- a/src/python/py_label.cpp
+++ b/src/python/py_label.cpp
@@ -128,10 +128,46 @@ static PyObject* Candidate_getScore(CandidateObject* self, void* closure)
 	}
 }
 
+static PyObject* Candidate_getCf(CandidateObject* self, void* closure)
+{
+	try
+	{
+		return py::buildPyValue(self->cand.cf);
+	}
+	catch (const bad_exception&)
+	{
+		return nullptr;
+	}
+	catch (const exception& e)
+	{
+		PyErr_SetString(PyExc_Exception, e.what());
+		return nullptr;
+	}
+}
+
+static PyObject* Candidate_getDf(CandidateObject* self, void* closure)
+{
+	try
+	{
+		return py::buildPyValue(self->cand.df);
+	}
+	catch (const bad_exception&)
+	{
+		return nullptr;
+	}
+	catch (const exception& e)
+	{
+		PyErr_SetString(PyExc_Exception, e.what());
+		return nullptr;
+	}
+}
+
 static PyGetSetDef Candidate_getseters[] = {
 	{ (char*)"words", (getter)Candidate_getWords, nullptr, Candidate_words__doc__, nullptr },
 	{ (char*)"name", (getter)Candidate_getName, (setter)Candidate_setName, Candidate_name__doc__, nullptr },
 	{ (char*)"score", (getter)Candidate_getScore, nullptr, Candidate_score__doc__, nullptr },
+	{ (char*)"cf", (getter)Candidate_getCf, nullptr, Candidate_cf__doc__, nullptr },
+	{ (char*)"df", (getter)Candidate_getDf, nullptr, Candidate_df__doc__, nullptr },
 	{ nullptr },
 };
 
@@ -265,13 +301,13 @@ static PyMethodDef Labeler_methods[] =
 
 static int PMIExtractor_init(ExtractorObject *self, PyObject *args, PyObject *kwargs)
 {
-	size_t minCf = 10, minDf = 5, minLen = 1, maxLen = 5, maxCand = 5000;
-	static const char* kwlist[] = { "min_cf", "min_df", "min_len", "max_len", "max_cand", nullptr };
-	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|nnnnn", (char**)kwlist,
-		&minCf, &minDf, &minLen, &maxLen, &maxCand)) return -1;
+	size_t minCf = 10, minDf = 5, minLen = 1, maxLen = 5, maxCand = 5000, normalized = 0;
+	static const char* kwlist[] = { "min_cf", "min_df", "min_len", "max_len", "max_cand", "normalized", nullptr };
+	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|nnnnnp", (char**)kwlist,
+		&minCf, &minDf, &minLen, &maxLen, &maxCand, &normalized)) return -1;
 	try
 	{
-		self->inst = new tomoto::label::PMIExtractor{ minCf, minDf, minLen, maxLen, maxCand };
+		self->inst = new tomoto::label::PMIExtractor{ minCf, minDf, minLen, maxLen, maxCand, !!normalized };
 	}
 	catch (const exception& e)
 	{
@@ -350,7 +386,7 @@ static int FoRelevance_init(LabelerObject *self, PyObject *args, PyObject *kwarg
 				pcands.emplace_back(&((CandidateObject*)item.get())->cand);
 			}
 		}
-		auto deref = [](tomoto::label::Candidate* p) { return *p; };
+		auto deref = [](tomoto::label::Candidate* p)->tomoto::label::Candidate& { return *p; };
 		self->inst = new tomoto::label::FoRelevance{ 
 			tm->inst, 
 			tomoto::makeTransformIter(pcands.begin(), deref), 
diff --git a/src/python/py_main.cpp b/src/python/py_main.cpp
index 40888be..b56270e 100644
--- a/src/python/py_main.cpp
+++ b/src/python/py_main.cpp
@@ -12,8 +12,8 @@ void char2Byte(const char* strBegin, const char* strEnd, vector<uint32_t>& start
 {
 	if (strBegin == strEnd) return;
 	vector<size_t> charPos;
-	auto it = strBegin, end = strEnd;
-	for (; it != end; )
+	auto it = strBegin;
+	for (; it != strEnd; )
 	{
 		charPos.emplace_back(it - strBegin);
 		uint8_t c = *it;
@@ -35,7 +35,7 @@ void char2Byte(const char* strBegin, const char* strEnd, vector<uint32_t>& start
 		}
 		else it += 1;
 	}
-	charPos.emplace_back(std::distance(strBegin, strEnd));
+	charPos.emplace_back(strEnd - strBegin);
 
 	for (size_t i = 0; i < startPos.size(); ++i)
 	{
@@ -50,6 +50,11 @@ void char2Byte(const string& str, vector<uint32_t>& startPos, vector<uint16_t>&
 	return char2Byte(&str[0], &str[0] + str.size(), startPos, length);
 }
 
+void char2Byte(const tomoto::SharedString& str, vector<uint32_t>& startPos, vector<uint16_t>& length)
+{
+	return char2Byte(str.begin(), str.end(), startPos, length);
+}
+
 void TopicModelObject::dealloc(TopicModelObject* self)
 {
 	DEBUG_LOG("TopicModelObject Dealloc " << self);
diff --git a/src/python/py_utils.cpp b/src/python/py_utils.cpp
index ae3361a..f39a42f 100644
--- a/src/python/py_utils.cpp
+++ b/src/python/py_utils.cpp
@@ -2,6 +2,7 @@
 #include "module.h"
 #include "utils.h"
 #include "label.h"
+#include "../Labeling/Phraser.hpp"
 #include "../Labeling/FoRelevance.h"
 
 using namespace std;
@@ -23,10 +24,10 @@ namespace py
 		PyObject* operator()(const tomoto::RawDoc& v)
 		{
 			PyObject* ret = PyTuple_New(5);
-			PyTuple_SET_ITEM(ret, 0, buildPyValue(v.words));
+			PyTuple_SET_ITEM(ret, 0, buildPyValue(v.words, py::force_list));
 			PyTuple_SET_ITEM(ret, 1, buildPyValue(v.rawStr));
-			PyTuple_SET_ITEM(ret, 2, buildPyValue(v.origWordPos));
-			PyTuple_SET_ITEM(ret, 3, buildPyValue(v.origWordLen));
+			PyTuple_SET_ITEM(ret, 2, buildPyValue(v.origWordPos, py::force_list));
+			PyTuple_SET_ITEM(ret, 3, buildPyValue(v.origWordLen, py::force_list));
 			PyObject* dict = PyDict_New();
 			for (auto& p : v.misc)
 			{
@@ -221,6 +222,26 @@ class RawDocWrapper
 	{
 		return raw.words[idx];
 	}
+
+	auto begin() const -> decltype(raw.words.begin())
+	{
+		return raw.words.begin();
+	}
+
+	auto end() const -> decltype(raw.words.end())
+	{
+		return raw.words.end();
+	}
+
+	auto rbegin() const -> decltype(raw.words.rbegin())
+	{
+		return raw.words.rbegin();
+	}
+
+	auto rend() const -> decltype(raw.words.rend())
+	{
+		return raw.words.rend();
+	}
 };
 
 
@@ -264,8 +285,13 @@ int CorpusObject::init(CorpusObject* self, PyObject* args, PyObject* kwargs)
 
 	if (!dep)
 	{
-		dep = PyObject_CallObject((PyObject*)&UtilsVocab_type, nullptr);
+		PyObject* utils = PyImport_AddModule("tomotopy.utils");
+		py::UniqueObj classCorpus{ PyObject_GetAttrString(utils, "Corpus") };
+		py::UniqueObj classVocabDict{ PyObject_GetAttrString(classCorpus, "_VocabDict") };
+
+		dep = PyObject_CallObject(classVocabDict, nullptr);
 		((VocabObject*)dep)->vocabs = new tomoto::Dictionary;
+		((VocabObject*)dep)->size = -1;
 	}
 	else Py_INCREF(dep);
 
@@ -355,18 +381,17 @@ PyObject* CorpusObject::setstate(CorpusObject* self, PyObject* args)
 			self->docs.emplace_back(move(doc));
 		}
 		if (PyErr_Occurred()) throw bad_exception{};
+		Py_INCREF(Py_None);
+		return Py_None;
 	}
 	catch (const bad_exception&)
 	{
-		return nullptr;
 	}
 	catch (const exception& e)
 	{
 		PyErr_SetString(PyExc_Exception, e.what());
-		return nullptr;
 	}
-	Py_INCREF(Py_None);
-	return Py_None;
+	return nullptr;
 }
 
 PyObject* CorpusObject::addDoc(CorpusObject* self, PyObject* args, PyObject* kwargs)
@@ -481,9 +506,10 @@ PyObject* CorpusObject::extractNgrams(CorpusObject* self, PyObject* args, PyObje
 {
 	size_t minCf = 10, minDf = 5, maxLen = 5, maxCand = 5000;
 	float minScore = -INFINITY;
-	static const char* kwlist[] = { "min_cf", "min_df", "max_len", "max_cand", "min_score", nullptr };
-	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|nnnnf", (char**)kwlist,
-		&minCf, &minDf, &maxLen, &maxCand, &minScore)) return nullptr;
+	size_t normalized = 0, workers = 1;
+	static const char* kwlist[] = { "min_cf", "min_df", "max_len", "max_cand", "min_score", "normalized", "workers", nullptr };
+	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|nnnnfpn", (char**)kwlist,
+		&minCf, &minDf, &maxLen, &maxCand, &minScore, &normalized, &workers)) return nullptr;
 	try
 	{
 		if (!self->isIndependent())
@@ -511,9 +537,9 @@ PyObject* CorpusObject::extractNgrams(CorpusObject* self, PyObject* args, PyObje
 		};
 		auto docBegin = tomoto::makeTransformIter(self->docs.begin(), tx);
 		auto docEnd = tomoto::makeTransformIter(self->docs.end(), tx);
-		auto cands = tomoto::label::extractPMINgrams(docBegin, docEnd,
+		auto cands = tomoto::phraser::extractPMINgrams(docBegin, docEnd,
 			cf, df,
-			minCf, minDf, 2, maxLen, maxCand, minScore
+			minCf, minDf, 2, maxLen, maxCand, minScore, normalized
 		);
 
 		PyObject* ret = PyList_New(0);
@@ -538,6 +564,7 @@ PyObject* CorpusObject::extractNgrams(CorpusObject* self, PyObject* args, PyObje
 	}
 }
 
+// TODO: It loses some ngram patterns. Fix me!
 PyObject* CorpusObject::concatNgrams(CorpusObject* self, PyObject* args, PyObject* kwargs)
 {
 	PyObject* cands;
@@ -672,9 +699,9 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 	try 
 	{
 		// indexing by int
-		if (PyLong_Check(idx))
+		Py_ssize_t v = PyLong_AsLongLong(idx);
+		if (v != -1 || !(PyErr_Occurred() && (PyErr_Clear(), true)))
 		{
-			Py_ssize_t v = PyLong_AsLong(idx);
 			if(v >= len(self) || -v > len(self)) throw out_of_range{ "IndexError: " + to_string(v) };
 			auto doc = (DocumentObject*)PyObject_CallFunctionObjArgs((PyObject*)&UtilsDocument_type, (PyObject*)self, nullptr);
 			if (!doc) throw bad_exception{};
@@ -704,6 +731,7 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 			if (self->isIndependent())
 			{
 				auto ret = (CorpusObject*)PyObject_CallFunctionObjArgs((PyObject*)&UtilsCorpus_type, (PyObject*)self->vocab, nullptr);
+				if (!ret) throw bad_exception{};
 				for (Py_ssize_t i = start; i < end; i += step)
 				{
 					ret->docs.emplace_back(self->docs[i]);
@@ -713,6 +741,7 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 			else if (self->made)
 			{
 				auto ret = (CorpusObject*)PyObject_CallFunctionObjArgs((PyObject*)&UtilsCorpus_type, (PyObject*)self->tm, nullptr);
+				if (!ret) throw bad_exception{};
 				for (Py_ssize_t i = start; i < end; i += step)
 				{
 					ret->docsMade.emplace_back(self->docsMade[i]);
@@ -723,6 +752,7 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 			else if(self->isSubDocs())
 			{
 				auto ret = (CorpusObject*)PyObject_CallFunctionObjArgs((PyObject*)&UtilsCorpus_type, (PyObject*)self->tm, nullptr);
+				if (!ret) throw bad_exception{};
 				for (Py_ssize_t i = start; i < end; i += step)
 				{
 					ret->docIdcs.emplace_back(self->docIdcs[i]);
@@ -733,6 +763,7 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 			else 
 			{
 				auto ret = (CorpusObject*)PyObject_CallFunctionObjArgs((PyObject*)&UtilsCorpus_type, (PyObject*)self->tm, nullptr);
+				if (!ret) throw bad_exception{};
 				for (Py_ssize_t i = start; i < end; i += step)
 				{
 					ret->docIdcs.emplace_back(i);
@@ -742,14 +773,14 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 			}
 		}
 		// indexing by list of uid or int
-		else if (PyList_Check(idx))
+		else if (py::UniqueObj{ PyObject_GetIter(idx) })
 		{
 			vector<size_t> idcs;
 			py::foreach<PyObject*>(idx, [&](PyObject* o)
 			{
-				if (PyLong_Check(o))
+				Py_ssize_t v = PyLong_AsLongLong(o);
+				if (v != -1 || !(PyErr_Occurred() && (PyErr_Clear(), true)))
 				{
-					Py_ssize_t v = PyLong_AsLong(o);
 					if (v >= len(self) || -v > len(self))
 					{
 						throw out_of_range{ "IndexError. len = " + to_string(len(self)) + ", idx = " + to_string(v) };
@@ -766,7 +797,7 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 				}
 				else
 				{
-					py::UniqueObj ty{ PyObject_Type(idx) };
+					py::UniqueObj ty{ PyObject_Type(o) };
 					py::UniqueObj repr{ PyObject_Str(ty) };
 					throw runtime_error{ string{"Unsupported indexing type "} + PyUnicode_AsUTF8(repr) };
 				}
@@ -775,6 +806,7 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 			if (self->isIndependent())
 			{
 				auto ret = (CorpusObject*)PyObject_CallFunctionObjArgs((PyObject*)&UtilsCorpus_type, (PyObject*)self->vocab, nullptr);
+				if (!ret) throw bad_exception{};
 				for (auto i : idcs)
 				{
 					ret->docs.emplace_back(self->docs[i]);
@@ -784,6 +816,7 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 			else if (self->made)
 			{
 				auto ret = (CorpusObject*)PyObject_CallFunctionObjArgs((PyObject*)&UtilsCorpus_type, (PyObject*)self->tm, nullptr);
+				if (!ret) throw bad_exception{};
 				for (auto i : ret->docIdcs)
 				{
 					ret->docsMade.emplace_back(self->docsMade[i]);
@@ -794,6 +827,7 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 			else
 			{
 				auto ret = (CorpusObject*)PyObject_CallFunctionObjArgs((PyObject*)&UtilsCorpus_type, (PyObject*)self->tm, nullptr);
+				if (!ret) throw bad_exception{};
 				ret->docIdcs = move(idcs);
 				for (auto i : ret->docIdcs)
 				{
@@ -811,18 +845,16 @@ PyObject* CorpusObject::getitem(CorpusObject* self, PyObject* idx)
 	}
 	catch (const bad_exception&)
 	{
-		return nullptr;
 	}
 	catch (const out_of_range& e)
 	{
 		PyErr_SetString(PyExc_KeyError, e.what());
-		return nullptr;
 	}
 	catch (const exception& e)
 	{
 		PyErr_SetString(PyExc_Exception, e.what());
-		return nullptr;
 	}
+	return nullptr;
 }
 
 PyObject* CorpusObject::iter(CorpusObject* self)
diff --git a/src/python/utils.h b/src/python/utils.h
index ba74bcd..af8be13 100644
--- a/src/python/utils.h
+++ b/src/python/utils.h
@@ -57,7 +57,7 @@ struct CorpusObject
 
 	inline bool isIndependent() const
 	{
-		return !!PyObject_TypeCheck(vocab, &UtilsVocab_type);
+		return vocab && !!PyObject_TypeCheck(vocab, &UtilsVocab_type);
 	}
 
 	inline bool isSubDocs() const
diff --git a/tomotopy/_summary.py b/tomotopy/_summary.py
index 102e77f..7b6b77f 100644
--- a/tomotopy/_summary.py
+++ b/tomotopy/_summary.py
@@ -228,7 +228,7 @@ def topics_info_HLDAModel(mdl, file, topic_word_top_n):
 
     def print_hierarchical(k=0, level=0):
         words = ' '.join(w for w, _ in mdl.get_topic_words(k, top_n=topic_word_top_n))
-        print('| {}#{} ({}) : {}'.format('  ' * level, k, topic_cnt[k], words), file=file)
+        print('| {}#{} ({}, {}) : {}'.format('  ' * level, k, topic_cnt[k], mdl.num_docs_of_topic(k), words), file=file)
         for c in np.sort(mdl.children_topics(k)):
             print_hierarchical(c, level + 1)
 
diff --git a/tomotopy/_version.py b/tomotopy/_version.py
index fcca93b..0fb6f35 100644
--- a/tomotopy/_version.py
+++ b/tomotopy/_version.py
@@ -1,1 +1,1 @@
-__version__ = '0.10.0'
\ No newline at end of file
+__version__ = '0.10.1'
\ No newline at end of file
diff --git a/tomotopy/utils.py b/tomotopy/utils.py
index 77f5a84..ad29311 100644
--- a/tomotopy/utils.py
+++ b/tomotopy/utils.py
@@ -102,7 +102,7 @@ def process(self, data_feeder):
         
         return num
 
-    def save(self, filename:str):
+    def save(self, filename:str, protocol=0):
         '''Save the current instance into the file `filename`. 
 
 Parameters
@@ -132,7 +132,7 @@ def load(filename:str):
     def __len__(self):
         return super().__len__()
     
-    def extract_ngrams(self, min_cf=10, min_df=5, max_len=5, max_cand=5000, min_score=float('-inf')):
+    def extract_ngrams(self, min_cf=10, min_df=5, max_len=5, max_cand=5000, min_score=float('-inf'), normalized=False, workers=0):
         '''..versionadded:: 0.10.0
 
 Extract frequent n-grams using PMI score
@@ -149,13 +149,18 @@ def extract_ngrams(self, min_cf=10, min_df=5, max_len=5, max_cand=5000, min_scor
     Maximum number of n-grams to be extracted
 min_score : float
     Minium PMI score of n-grams to be extracted
+normalized : bool
+    whether to use Normalized PMI or just PMI
+workers : int
+    an integer indicating the number of workers to perform samplings. 
+    If `workers` is 0, the number of cores in the system will be used.
 
 Returns
 -------
 candidates : List[tomotopy.label.Candidate]
     The extracted n-gram candidates in `tomotopy.label.Candidate` type
         '''
-        return super().extract_ngrams(min_cf, min_df, max_len, max_cand, min_score)
+        return super().extract_ngrams(min_cf, min_df, max_len, max_cand, min_score, normalized, workers)
     
     def concat_ngrams(self, cands, delimiter='_'):
         '''..versionadded:: 0.10.0
@@ -196,13 +201,17 @@ def __init__(self, stemmer=None, pattern:str=None, lowercase=True):
 
     def __call__(self, raw:str, user_data=None):
         if self._stemmer:
-            for g in self._pat.finditer(raw.lower() if self._lowercase else raw):
+            for g in self._pat.finditer(raw if self._lowercase else raw):
                 start, end = g.span(0)
-                yield self._stemmer(g.group(0)), start, end - start
+                word = g.group(0)
+                if self._lowercase: word = word.lower()
+                yield self._stemmer(word), start, end - start
         else:
-            for g in self._pat.finditer(raw.lower() if self._lowercase else raw):
+            for g in self._pat.finditer(raw if self._lowercase else raw):
                 start, end = g.span(0)
-                yield g.group(0), start, end - start
+                word = g.group(0)
+                if self._lowercase: word = word.lower()
+                yield word, start, end - start
 
 import os
 if os.environ.get('TOMOTOPY_LANG') == 'kr':


Test patch:
diff --git a/test/unit_test.py b/test/unit_test.py
index ad7df85..dd749f5 100644
--- a/test/unit_test.py
+++ b/test/unit_test.py
@@ -1,48 +1,51 @@
+import sys, os
 import tomotopy as tp
-import sys
+
+curpath = os.path.dirname(os.path.realpath(__file__))
+print(curpath)
 
 model_cases = [
-    (tp.LDAModel, 'test/sample.txt', 0, None, {'k':40}, None),
-    (tp.LLDAModel, 'test/sample_with_md.txt', 1, lambda x:x, {'k':5}, None),
-    (tp.PLDAModel, 'test/sample_with_md.txt', 0, None, {'latent_topics':2, 'topics_per_label':2}, None),
-    (tp.PLDAModel, 'test/sample_with_md.txt', 1, lambda x:x, {'latent_topics':2, 'topics_per_label':2}, None),
-    (tp.HLDAModel, 'test/sample.txt', 0, None, {'depth':3}, [tp.ParallelScheme.NONE]),
-    (tp.CTModel, 'test/sample.txt', 0, None, {'k':10}, None),
-    (tp.HDPModel, 'test/sample.txt', 0, None, {'initial_k':10}, [tp.ParallelScheme.COPY_MERGE]),
-    (tp.MGLDAModel, 'test/sample.txt', 0, None, {'k_g':5, 'k_l':5}, None),
-    (tp.PAModel, 'test/sample.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
-    (tp.HPAModel, 'test/sample.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
-    (tp.DMRModel, 'test/sample_with_md.txt', 1, lambda x:'_'.join(x), {'k':10}, None),
-    (tp.SLDAModel, 'test/sample_with_md.txt', 1, lambda x:list(map(float, x)), {'k':10, 'vars':'b'}, None),
-    (tp.DTModel, 'test/sample_tp.txt', 1, lambda x:int(x[0]), {'k':10, 't':13}, None),
-    (tp.GDMRModel, 'test/sample_tp.txt', 1, lambda x:list(map(float, x)), {'k':10, 'degrees':[3]}, None),
+    (tp.LDAModel, curpath + '/sample.txt', 0, None, {'k':40}, None),
+    (tp.LLDAModel, curpath + '/sample_with_md.txt', 1, lambda x:x, {'k':5}, None),
+    (tp.PLDAModel, curpath + '/sample_with_md.txt', 0, None, {'latent_topics':2, 'topics_per_label':2}, None),
+    (tp.PLDAModel, curpath + '/sample_with_md.txt', 1, lambda x:x, {'latent_topics':2, 'topics_per_label':2}, None),
+    (tp.HLDAModel, curpath + '/sample.txt', 0, None, {'depth':3}, None),
+    (tp.CTModel, curpath + '/sample.txt', 0, None, {'k':10}, None),
+    (tp.HDPModel, curpath + '/sample.txt', 0, None, {'initial_k':10}, [tp.ParallelScheme.COPY_MERGE]),
+    (tp.MGLDAModel, curpath + '/sample.txt', 0, None, {'k_g':5, 'k_l':5}, None),
+    (tp.PAModel, curpath + '/sample.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
+    (tp.HPAModel, curpath + '/sample.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
+    (tp.DMRModel, curpath + '/sample_with_md.txt', 1, lambda x:'_'.join(x), {'k':10}, None),
+    (tp.SLDAModel, curpath + '/sample_with_md.txt', 1, lambda x:list(map(float, x)), {'k':10, 'vars':'b'}, None),
+    (tp.DTModel, curpath + '/sample_tp.txt', 1, lambda x:int(x[0]), {'k':10, 't':13}, None),
+    (tp.GDMRModel, curpath + '/sample_tp.txt', 1, lambda x:list(map(float, x)), {'k':10, 'degrees':[3]}, None),
 ]
 
 model_raw_cases = [
-    (tp.LDAModel, 'test/sample_raw.txt', 0, None, {'k':10}, None),
-    (tp.HLDAModel, 'test/sample_raw.txt', 0, None, {'depth':3}, [tp.ParallelScheme.NONE]),
-    (tp.CTModel, 'test/sample_raw.txt', 0, None, {'k':10}, None),
-    (tp.HDPModel, 'test/sample_raw.txt', 0, None, {'initial_k':10}, [tp.ParallelScheme.COPY_MERGE]),
-    (tp.MGLDAModel, 'test/sample_raw.txt', 0, None, {'k_g':5, 'k_l':5}, None),
-    (tp.PAModel, 'test/sample_raw.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
-    (tp.HPAModel, 'test/sample_raw.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
+    (tp.LDAModel, curpath + '/sample_raw.txt', 0, None, {'k':10}, None),
+    (tp.HLDAModel, curpath + '/sample_raw.txt', 0, None, {'depth':3}, None),
+    (tp.CTModel, curpath + '/sample_raw.txt', 0, None, {'k':10}, None),
+    (tp.HDPModel, curpath + '/sample_raw.txt', 0, None, {'initial_k':10}, [tp.ParallelScheme.COPY_MERGE]),
+    (tp.MGLDAModel, curpath + '/sample_raw.txt', 0, None, {'k_g':5, 'k_l':5}, None),
+    (tp.PAModel, curpath + '/sample_raw.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
+    (tp.HPAModel, curpath + '/sample_raw.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
 ]
 
 model_corpus_cases = [
-    (tp.LDAModel, 'test/sample.txt', 0, None, {'k':10}, None),
-    (tp.LLDAModel, 'test/sample_with_md.txt', 1, lambda x:{'labels':x}, {'k':5}, None),
-    (tp.PLDAModel, 'test/sample_with_md.txt', 0, None, {'latent_topics':2, 'topics_per_label':2}, None),
-    (tp.PLDAModel, 'test/sample_with_md.txt', 1, lambda x:{'labels':x}, {'latent_topics':2, 'topics_per_label':2}, None),
-    (tp.HLDAModel, 'test/sample.txt', 0, None, {'depth':3}, [tp.ParallelScheme.NONE]),
-    (tp.CTModel, 'test/sample.txt', 0, None, {'k':10}, None),
-    (tp.HDPModel, 'test/sample.txt', 0, None, {'initial_k':10}, [tp.ParallelScheme.COPY_MERGE]),
-    (tp.MGLDAModel, 'test/sample.txt', 0, None, {'k_g':5, 'k_l':5}, None),
-    (tp.PAModel, 'test/sample.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
-    (tp.HPAModel, 'test/sample.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
-    (tp.DMRModel, 'test/sample_with_md.txt', 1, lambda x:{'metadata':'_'.join(x)}, {'k':10}, None),
-    (tp.SLDAModel, 'test/sample_with_md.txt', 1, lambda x:{'y':list(map(float, x))}, {'k':10, 'vars':'b'}, None),
-    (tp.DTModel, 'test/sample_tp.txt', 1, lambda x:{'timepoint':int(x[0])}, {'k':10, 't':13}, None),
-    (tp.GDMRModel, 'test/sample_tp.txt', 1, lambda x:{'metadata':list(map(float, x))}, {'k':10, 'degrees':[3]}, None),
+    (tp.LDAModel, curpath + '/sample.txt', 0, None, {'k':10}, None),
+    (tp.LLDAModel, curpath + '/sample_with_md.txt', 1, lambda x:{'labels':x}, {'k':5}, None),
+    (tp.PLDAModel, curpath + '/sample_with_md.txt', 0, None, {'latent_topics':2, 'topics_per_label':2}, None),
+    (tp.PLDAModel, curpath + '/sample_with_md.txt', 1, lambda x:{'labels':x}, {'latent_topics':2, 'topics_per_label':2}, None),
+    (tp.HLDAModel, curpath + '/sample.txt', 0, None, {'depth':3}, None),
+    (tp.CTModel, curpath + '/sample.txt', 0, None, {'k':10}, None),
+    (tp.HDPModel, curpath + '/sample.txt', 0, None, {'initial_k':10}, [tp.ParallelScheme.COPY_MERGE]),
+    (tp.MGLDAModel, curpath + '/sample.txt', 0, None, {'k_g':5, 'k_l':5}, None),
+    (tp.PAModel, curpath + '/sample.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
+    (tp.HPAModel, curpath + '/sample.txt', 0, None, {'k1':5, 'k2':10}, [tp.ParallelScheme.COPY_MERGE]),
+    (tp.DMRModel, curpath + '/sample_with_md.txt', 1, lambda x:{'metadata':'_'.join(x)}, {'k':10}, None),
+    (tp.SLDAModel, curpath + '/sample_with_md.txt', 1, lambda x:{'y':list(map(float, x))}, {'k':10, 'vars':'b'}, None),
+    (tp.DTModel, curpath + '/sample_tp.txt', 1, lambda x:{'timepoint':int(x[0])}, {'k':10, 't':13}, None),
+    (tp.GDMRModel, curpath + '/sample_tp.txt', 1, lambda x:{'metadata':list(map(float, x))}, {'k':10, 'degrees':[3]}, None),
 ]
 
 def train1(cls, inputFile, mdFields, f, kargs, ps):
@@ -213,7 +216,7 @@ def train_multi_corpus(cls, inputFile, mdFields, f, kargs, ps):
     print('Corpus2')
     for d in tcorpus2[:10]: print(d.get_ll())
 
-def test_estimate_SLDA_PARTITION(cls=tp.SLDAModel, inputFile='test/sample_with_md.txt', mdFields=1, f=lambda x:list(map(float, x)), kargs={'k':10, 'vars':'b'}, ps=tp.ParallelScheme.PARTITION):
+def test_estimate_SLDA_PARTITION(cls=tp.SLDAModel, inputFile=curpath + '/sample_with_md.txt', mdFields=1, f=lambda x:list(map(float, x)), kargs={'k':10, 'vars':'b'}, ps=tp.ParallelScheme.PARTITION):
     print('Test estimate')
     tw = 0
     print('Initialize model %s with TW=%s ...' % (str(cls), ['one', 'idf', 'pmi'][tw]))
@@ -244,7 +247,7 @@ def test_auto_labeling():
     corpus = tp.utils.Corpus(tokenizer=tp.utils.SimpleTokenizer(stemmer=stemmer.stem), 
         stopwords=lambda x: len(x) <= 2 or x in stopwords)
     # data_feeder yields a tuple of (raw string, user data) or a str (raw string)
-    corpus.process(open('test/sample_raw.txt', encoding='utf-8'))
+    corpus.process(open(curpath + '/sample_raw.txt', encoding='utf-8'))
     
     ngrams = corpus.extract_ngrams(min_cf=5, min_df=3)
     for c in ngrams:
@@ -280,7 +283,7 @@ def test_docs():
     corpus = tp.utils.Corpus(tokenizer=tp.utils.SimpleTokenizer(stemmer=stemmer.stem), 
         stopwords=lambda x: len(x) <= 2 or x in stopwords)
     # data_feeder yields a tuple of (raw string, user data) or a str (raw string)
-    for i, line in enumerate(open('test/sample_raw.txt', encoding='utf-8')):
+    for i, line in enumerate(open(curpath + '/sample_raw.txt', encoding='utf-8')):
         corpus.add_doc(raw=line, uid='doc{:05}'.format(i), etc=len(line))
 
     def _test_doc(doc, etc=False):
@@ -318,7 +321,7 @@ def test_corpus_transform():
     corpus = tp.utils.Corpus(tokenizer=tp.utils.SimpleTokenizer(stemmer=stemmer.stem), 
         stopwords=lambda x: len(x) <= 2 or x in stopwords)
     # data_feeder yields a tuple of (raw string, user data) or a str (raw string)
-    for i, line in enumerate(open('test/sample_raw.txt', encoding='utf-8')):
+    for i, line in enumerate(open(curpath + '/sample_raw.txt', encoding='utf-8')):
         corpus.add_doc(raw=line, uid='doc{:05}'.format(i), metadata=line[0])
     
     def xform(misc):
@@ -330,7 +333,7 @@ def xform(misc):
 
 def test_hdp_to_lda():
     mdl = tp.HDPModel(tw=tp.TermWeight.ONE, min_df=5, rm_top=5, alpha=0.5, gamma=0.5, initial_k=5)
-    for n, line in enumerate(open('test/sample.txt', encoding='utf-8')):
+    for n, line in enumerate(open(curpath + '/sample.txt', encoding='utf-8')):
         ch = line.strip().split()
         mdl.add_doc(ch)
     mdl.burn_in = 100
@@ -354,7 +357,7 @@ def test_hdp_to_lda():
 
 def test_coherence():
     mdl = tp.LDAModel(tw=tp.TermWeight.ONE, k=20, min_df=5, rm_top=5)
-    for n, line in enumerate(open('test/sample.txt', encoding='utf-8')):
+    for n, line in enumerate(open(curpath + '/sample.txt', encoding='utf-8')):
         ch = line.strip().split()
         mdl.add_doc(ch)
     mdl.train(1000)