<jupyter_start><jupyter_text>Topic Modeling on Roman Science*Before you run anything, write down your answer:* What do you expect are the subjects of science during the Roman Empire?In this notebook we will be using a **topic model** to explore a Roman encyclopedia.There are three cells with questions towards the end of the notebook.Once you are done, write your answer again (what is Roman science), having looked at Pliny's encyclopedia.<jupyter_code>import numpy as np
%load_ext cython

import re, sys, random, math
from collections import Counter
from timeit import default_timer as timer

from matplotlib import pyplot

word_pattern = re.compile("\w[\w\-\']*\w|\w")
source_directory = "../data/PlinyTheElder"

num_topics = 30
doc_smoothing = 0.5
word_smoothing = 0.01
%%cython

from cython.view cimport array as cvarray
import numpy as np
import random
from timeit import default_timer as timer

class Document:
    
    def __init__(self, long[:] doc_tokens, long[:] doc_topics, long[:] topic_changes, long[:] doc_topic_counts):
        self.doc_tokens = doc_tokens
        self.doc_topics = doc_topics
        self.topic_changes = topic_changes
        self.doc_topic_counts = doc_topic_counts

cdef class TopicModel:
    
    cdef long[:] topic_totals
    cdef long[:,:] word_topics
    cdef int num_topics
    cdef int vocab_size
    
    cdef double[:] topic_probs
    cdef double[:] topic_normalizers
    cdef float doc_smoothing
    cdef float word_smoothing
    cdef float smoothing_times_vocab_size
    
    documents = []
    vocabulary = []
    
    def __init__(self, num_topics, vocabulary, doc_smoothing, word_smoothing):
        self.num_topics = num_topics
        self.vocabulary.extend(vocabulary)
        self.vocab_size = len(vocabulary)
        
        self.doc_smoothing = doc_smoothing
        self.word_smoothing = word_smoothing
        self.smoothing_times_vocab_size = word_smoothing * self.vocab_size
        
        self.topic_totals = np.zeros(num_topics, dtype=int)
        self.word_topics = np.zeros((self.vocab_size, num_topics), dtype=int)
    
    def add_document(self, doc):
        cdef int word_id, topic
        
        self.documents.append(doc)
        
        for i in range(len(doc.doc_tokens)):
            word_id = doc.doc_tokens[i]
            topic = doc.doc_topics[i]
            
            self.word_topics[word_id,topic] += 1
            self.topic_totals[topic] += 1
            doc.doc_topic_counts[topic] += 1
            
    def sample(self, iterations):
        cdef int old_topic, new_topic, word_id, topic, i, doc_length
        cdef double sampling_sum = 0
        cdef double sample
        cdef long[:] word_topic_counts
        
        cdef long[:] doc_tokens
        cdef long[:] doc_topics
        cdef long[:] doc_topic_counts
        cdef long[:] topic_changes
        
        cdef double[:] uniform_variates
        cdef double[:] topic_probs = np.zeros(self.num_topics, dtype=float)
        cdef double[:] topic_normalizers = np.zeros(self.num_topics, dtype=float)
        
        for topic in range(self.num_topics):
            topic_normalizers[topic] = 1.0 / (self.topic_totals[topic] + self.smoothing_times_vocab_size)
        
        for iteration in range(iterations):
            for document in self.documents:
                doc_tokens = document.doc_tokens
                doc_topics = document.doc_topics
                doc_topic_counts = document.doc_topic_counts
                topic_changes = document.topic_changes
                
                doc_length = len(document.doc_tokens)
                uniform_variates = np.random.random_sample(doc_length)
                
                for i in range(doc_length):
                    word_id = doc_tokens[i]
                    old_topic = doc_topics[i]
                    word_topic_counts = self.word_topics[word_id,:]
        
                    ## erase the effect of this token
                    word_topic_counts[old_topic] -= 1
                    self.topic_totals[old_topic] -= 1
                    doc_topic_counts[old_topic] -= 1
        
                    topic_normalizers[old_topic] = 1.0 / (self.topic_totals[old_topic] + self.smoothing_times_vocab_size)
        
                    ###
                    ### SAMPLING DISTRIBUTION
                    ###
        
                    sampling_sum = 0.0
                    for topic in range(self.num_topics):
                        topic_probs[topic] = (doc_topic_counts[topic] + self.doc_smoothing) * (word_topic_counts[topic] + self.word_smoothing) * topic_normalizers[topic]
                        sampling_sum += topic_probs[topic]

                    sample = uniform_variates[i] * sampling_sum
        
                    new_topic = 0
                    while sample > topic_probs[new_topic]:
                        sample -= topic_probs[new_topic]
                        new_topic += 1
            
                    ## add the effect of this token back in
                    word_topic_counts[new_topic] += 1
                    self.topic_totals[new_topic] += 1
                    doc_topic_counts[new_topic] += 1
                    topic_normalizers[new_topic] = 1.0 / (self.topic_totals[new_topic] + self.smoothing_times_vocab_size)

                    doc_topics[i] = new_topic
        
                    if new_topic != old_topic:
                        #pass
                        topic_changes[i] += 1

    def topic_words(self, int topic, n_words=12):
        sorted_words = sorted(zip(self.word_topics[:,topic], self.vocabulary), reverse=True)
        return " ".join([w for x, w in sorted_words[:n_words]])

    def print_all_topics(self):
        for topic in range(self.num_topics):
            print(topic, self.topic_words(topic))
## Read the stoplist file

stoplist = set()
with open("{}/stoplist.txt".format(source_directory), encoding="utf-8") as stop_reader:
    for line in stop_reader:
        line = line.rstrip()
        stoplist.add(line)


## Read the documents file
        
word_counts = Counter()
documents = []

for line in open("{}/documents.txt".format(source_directory), encoding="utf-8"):
    #line = line.lower()
    
    tokens = word_pattern.findall(line)
    
    ## remove stopwords, short words, and upper-cased words
    tokens = [w for w in tokens if not w in stoplist and len(w) >= 3 and not w[0].isupper()]
    word_counts.update(tokens)
    
    doc_topic_counts = np.zeros(num_topics, dtype=int)
    
    documents.append({ "original": line, "token_strings": tokens, "topic_counts": doc_topic_counts })


## Now that we're done reading from disk, we can count the total
##  number of words.

vocabulary = list(word_counts.keys())
vocabulary_size = len(vocabulary)
word_ids = { w: i for (i, w) in enumerate(vocabulary) }
smoothing_times_vocab_size = word_smoothing * vocabulary_size

#word_topics = np.zeros((len(vocabulary), num_topics), dtype=int)
#topic_totals = np.zeros(num_topics, dtype=int)

for document in documents:
    tokens = document["token_strings"]
    doc_topic_counts = document["topic_counts"]
    
    doc_tokens = np.ndarray(len(tokens), dtype=int)
    doc_topics = np.ndarray(len(tokens), dtype=int)
    topic_changes = np.zeros(len(tokens), dtype=int)
    
    for i, w in enumerate(tokens):
        word_id = word_ids[w]
        topic = random.randrange(num_topics)
        
        doc_tokens[i] = word_id
        doc_topics[i] = topic
        
        ## Update counts: 
        #word_topics[word_id][topic] += 1
        #topic_totals[topic] += 1
        #doc_topic_counts[topic] += 1
    
    document["doc_tokens"] = doc_tokens
    document["doc_topics"] = doc_topics
    document["topic_changes"] = topic_changes<jupyter_output><empty_output><jupyter_text>This cell actually runs the modelThis may take some time to run. It will print "Done!" at the end.* Do the topics you see make sense? Which ones are most meaningful, and which ones are least meaningful? Copy specific examples.* Are the topics consistent? Compare your results to others at your table. Write specific examples of similar and different (but recognizably related) versions of topics.[Answer below]<jupyter_code>model = TopicModel(num_topics, vocabulary, doc_smoothing, word_smoothing)

for document in documents:
    c_doc = Document(document["doc_tokens"], document["doc_topics"], document["topic_changes"], document["topic_counts"])
    model.add_document(c_doc)

sampling_dist = np.zeros(num_topics, dtype=float)
#topic_normalizers = np.zeros(num_topics, dtype=float)
#for topic in range(num_topics):
#    topic_normalizers[topic] = 1.0 / (topic_totals[topic] + smoothing_times_vocab_size)

doc_topic_probs = np.zeros((len(model.documents), num_topics))
word_topic_probs = np.zeros((len(vocabulary), num_topics))

# Initial burn-in iterations
for i in range(20):
    start = timer()
    model.sample(50)
    elapsed_time = timer() - start
    print("Iteration {}, {:.2f} seconds per iteration".format(i * 50, elapsed_time / 50))
    model.print_all_topics()
    print()

# Saved samples
for i in range(5):
    model.sample(10)
    
    for doc_id, doc in enumerate(model.documents):
        for word_id, topic in zip(doc.doc_tokens, doc.doc_topics):
            doc_topic_probs[doc_id,topic] += 1
            word_topic_probs[word_id,topic] += 1

            
print("Done!")
            
# Normalize
doc_row_sums = doc_topic_probs.sum(axis=1)
doc_topic_probs /= doc_row_sums[:,np.newaxis]

word_col_sums = word_topic_probs.sum(axis=0)
word_topic_probs /= word_col_sums[np.newaxis,:]<jupyter_output><empty_output><jupyter_text>This cell shows the prevalence of each topic from the beginning of the work to the end.* Which topics are most concentrated? Which topics are most dispersed through the work? How does this concentration/dispersion relate to their interpretability?* What can you infer about the structure of Pliny's encyclopedia?[Answer below]<jupyter_code>for topic in range(num_topics):
    print(topic, model.topic_words(topic, n_words=6))
    pyplot.plot(doc_topic_probs[:,topic])
    pyplot.show()<jupyter_output><empty_output><jupyter_text>This cell prints the documents with the largest proportion of a specified topic.* Does looking at example documents help explain what the model is finding? Provide positive and negative examples if possible. Add cells to look at at least five topics.[Answer below]<jupyter_code>def top_docs(topic, n_docs=10):
    for doc_id in np.argsort(-doc_topic_probs[:,topic])[:n_docs]:
        print(doc_id, doc_topic_probs[doc_id,topic], documents[doc_id]["original"])
top_docs(6)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09236684171042761, 'fraction_numerical': 0.005532633158289573, 'mean_word_length': 2.776558073654391, 'pattern_counts': {'":': 3, '<': 17, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '346108', 'license': 'Apache-2.0', 'n_tokens_mistral': 3400, 'n_tokens_neox': 3262, 'n_words': 807, 'path': 'week10/TopicModels.ipynb', 'repo_name': 'mimno/info-3350-fall-2019', 'size': '10664'}