Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

added source

  • Loading branch information...
commit 7b146cd1d21bce18763f0afdb34ff8ad1f8db96b 1 parent d71cb49
@jacobeisenstein authored
Showing with 1,064 additions and 0 deletions.
  1. +12 −0 bounds/scoreDoc.m
  2. +4 −0 bounds/scoreWords.m
  3. +16 −0 computePerplexity.m
  4. BIN  data/20news.mat
  5. BIN  data/nips.mat
  6. +24 −0 estimate/computeBetaDirichlet.m
  7. +31 −0 estimate/computeBetaSparseVariational.m
  8. +20 −0 estimate/evalLogNormal.m
  9. +27 −0 estimate/evalQTauLogAMinus1.m
  10. +10 −0 makeTopicReport.m
  11. +25 −0 preprocess.m
  12. +32 −0 runLDA.m
  13. +275 −0 sparseTAM.m
  14. +4 −0 startup.m
  15. +23 −0 tamEStep.m
  16. +272 −0 tamEStepInnerMex.c
  17. +54 −0 tamEStepNoMex.m
  18. +22 −0 utils/bayesian/fitDirichletPrior.m
  19. +5 −0 utils/bayesian/kldirichlet.m
  20. +15 −0 utils/holdoutWords.m
  21. +7 −0 utils/iterator/newDeltaIterator.m
  22. +10 −0 utils/iterator/newIterator.m
  23. +13 −0 utils/iterator/updateDeltaIterator.m
  24. +16 −0 utils/iterator/updateIterator.m
  25. +7 −0 utils/logNormalizeRows.m
  26. +4 −0 utils/logNormalizeVec.m
  27. +15 −0 utils/logToSimplex2.m
  28. +8 −0 utils/makeLogProbs.m
  29. +27 −0 utils/newton.m
  30. +10 −0 utils/normalize_rows.m
  31. +3 −0  utils/normalize_vec.m
  32. +8 −0 utils/randsample.m
  33. +15 −0 utils/sortidxs.m
  34. +3 −0  utils/spdiag.m
  35. +11 −0 utils/termCountsToWordList.m
  36. +36 −0 utils/termCountsToWordListMex.c
  37. BIN  utils/termCountsToWordListMex.mexglx
View
12 bounds/scoreDoc.m
@@ -0,0 +1,12 @@
+function [bound word_score lv_score] = scoreDoc(counts,beta,phi,x,sigma,alpha,e_log_theta)
+word_score = scoreWords(counts,beta);
+if nargin < 7,
+ e_log_theta = digamma(sigma) - digamma(sum(sigma));
+end
+if numel(alpha) == 1 && numel(sigma) > 1, alpha = repmat(alpha,1,numel(sigma)); end
+lv_score = e_log_theta * phi' * x';% ... %E[log p(z | theta)]
+lv_score = lv_score - x * sum(log(phi).*phi,2); %E[log q(phi)]
+lv_score = lv_score + gammaln(sum(alpha)) - sum(gammaln(alpha)) + e_log_theta * (alpha' - 1); % ... %E[log p(theta | alpha)]
+lv_score = lv_score - gammaln(sum(sigma)) + sum(gammaln(sigma)) - e_log_theta * (sigma' - 1); %E[log q(theta)]
+bound = word_score + lv_score;
+end
View
4 bounds/scoreWords.m
@@ -0,0 +1,4 @@
+function score = scoreWords(counts,beta)
+%function score = scoreWords(counts,beta)
+ score = sum(sum(counts .* beta));% - logsumexp(beta') * sum(counts,2);
+end
View
16 computePerplexity.m
@@ -0,0 +1,16 @@
+function [perplex ll ll_per_word] = computePerplexity(docs,topics,alpha,varargin)
+%function [perplex ll ll_per_word] = computePerplexity(docs,topics,alpha,varargin)
+%calls out to slow 3rd-party code (Wallach, Murray, Salakhutdinov, Mimno)
+%for computing perplexity, the right way
+
+[num_its ] =process_options(varargin,'num-its',100);
+doc_log_prob = 0;
+fprintf('computing perplexity');
+for i = 1:size(docs,1)
+ doc_log_prob = doc_log_prob + ldae_chibms(termCountsToWordList(docs(i,:)),topics,alpha,num_its);
+ if rem(i,10)==0, fprintf('='); end
+end
+fprintf('\n');
+ll_per_word = doc_log_prob / sum(sum(docs));
+perplex = exp(-ll_per_word);
+end
View
BIN  data/20news.mat
Binary file not shown
View
BIN  data/nips.mat
Binary file not shown
View
24 estimate/computeBetaDirichlet.m
@@ -0,0 +1,24 @@
+function [e_log_beta score lv_score word_score] = computeBetaDirichlet(ecounts,eta)
+% function beta = computeBetaBoring(ecounts)
+%
+% standard variational Bayesian treatment of E[log beta], where
+% beta ~ Dirichlet(eta)
+%
+% E[log(beta)] = digamma(counts + eta) - digamma(sum(counts + eta))
+% accepts multiple rows of ecounts
+% eta can be a vector or a scalar (indicating a symmetric prior)
+if nargin == 1, eta= 0; end
+[K W] = size(ecounts);
+if isscalar(eta), eta = repmat(eta,1,W); end
+
+word_score = 0; lv_score = 0;
+e_log_beta = zeros(size(ecounts));
+for k = 1:K
+ e_log_beta(k,:) = digamma(ecounts(k,:) + eta) - digamma(sum(ecounts(k,:)+eta));
+ if nargout > 1
+ word_score = word_score + ecounts(k,:)*e_log_beta(k,:)';
+ lv_score = lv_score - kldirichlet(ecounts(k,:)+eta,eta);
+ end
+end
+score = lv_score + word_score;
+end
View
31 estimate/computeBetaSparseVariational.m
@@ -0,0 +1,31 @@
+function eta = computeBetaSparseVariational(ecounts,eq_m,varargin)
+%function [eta bound] =
+%computeBetaSparseVariational(ecounts,eq_m,varargin)
+% newton optimization, variational EM for tau.
+[max_its verbose init_eta min_eta max_inv_tau] = ...
+ process_options(varargin,'max-its',1,...
+ 'verbose',false,'init-eta',[],'min-eta',1e-20,...
+ 'max-inv-tau',1e5);
+[W K] = size(ecounts); %eta = zeros(size(ecounts));
+
+if isempty(init_eta),
+ eta = zeros(W,1);
+ eq_inv_tau = ones(size(eta));
+else
+ eta = init_eta;
+ eta(eta.^2<min_eta.^2) = sign(eta(eta.^2<min_eta^2))*min_eta;
+ eq_inv_tau = 1./(eta.^2);
+end
+
+if ~verbose, fprintf('.'); end
+
+em_iter = newDeltaIterator(max_its,'debug',verbose,'thresh',1e-4);
+
+exp_eq_m = exp(eq_m);
+while ~(em_iter.done)
+ eta = newton(@evalLogNormal,eta,{ecounts,exp_eq_m,eq_inv_tau},'debug',verbose==1,'alpha',.1,'max-its',10000);
+ eq_inv_tau = 1./(eta.^2);
+ eq_inv_tau(eq_inv_tau >= max_inv_tau) = max_inv_tau;
+ em_iter = updateDeltaIterator(em_iter,eta);
+end
+end
View
20 estimate/evalLogNormal.m
@@ -0,0 +1,20 @@
+function [l g step] = evalLogNormal(eta,counts,exp_eq_m,invsigsq)
+Ck = sum(counts,1); C = sum(Ck); [W K] = size(counts);
+denom = repmat(exp(eta),1,K).*exp_eq_m;
+l = -(sum(eta' * counts) - Ck * log ( sum(denom) )' - 0.5 * trace(eta' * spdiag(invsigsq) * eta));
+if nargout > 1
+ beta = Ck * normalize_rows(denom') / (C + 1e-10); %expected beta
+ g = -(sum(counts,2) - C * beta' - invsigsq .* eta);
+ if nargout > 2
+ avec = -1./ (C * beta' + invsigsq);
+ a_times_g = (-g) .* avec;
+ c_times_a_times_beta = C * avec .* beta';
+ %step = -a_times_g + c_times_a_times_beta ./ (1 + beta * c_times_a_times_beta) * (beta * a_times_g);
+ step = -a_times_g + c_times_a_times_beta .* (beta * a_times_g ./ (1 + beta * c_times_a_times_beta));
+ end
+end
+end
+
+%these are slower
+%l = -(sum(eta' * counts) - Ck * log ( exp(eta)' * exp(eq_m) )' - 0.5 * sum(eta.*eta.*invsigsq));
+%l = -(sum(eta' * counts) - Ck * log ( exp(eta)' * exp(eq_m) )' - 0.5 * (eta.^2)'*invsigsq);
View
27 estimate/evalQTauLogAMinus1.m
@@ -0,0 +1,27 @@
+function [l g step] = evalQTauLogAMinus1(log_a_minus_1,b,etasq)
+%function [l g step] = evalQTauA(a,b,etasq)
+%
+%Q(tau) = Gamma(a,b)
+%eta \sim N(0,tau)
+%tau \sim 1/tau
+%
+%a is Wx1
+%b is Wx1
+%eta is Wx1
+
+a = exp(log_a_minus_1) + 1;
+
+%l = -.5 E[log tau] - .5 eta.^eta.^ E[1/tau] - E[log tau] - E[log Q(tau)]
+log_b = log(b);
+l = -(a + .5).*(digamma(a) + log_b) - .5 * etasq ./ ((a-1).*b) + a + gammaln(a) + a .* log_b;
+g = -(a + .5).*(trigamma(a)) + .5 * etasq ./ (b .* (a-1) .* (a-1)) + 1;
+
+l = -sum(l);
+g = -g .* (a - 1);
+
+step = 0; %do this later
+if nargout > 2
+ error('newton step size not yet implemented');
+end
+
+end
View
10 makeTopicReport.m
@@ -0,0 +1,10 @@
+function makeTopicReport(eta,vocab,varargin)
+%function makeTopicReport(eta,vocab,varargin)
+ [N background] = process_options(varargin,'N',10,'background',mean(eta));
+ for i = 1:size(eta,1)
+ if (size(eta,1)>1), fprintf('%d. ',i); end
+ num_legit = sum(abs(eta(i,:) - background)>1e-3);
+ fprintf('%s ',vocab{sortidxs(eta(i,:) - background,2,'descend',min(N,num_legit))})
+ fprintf('\n');
+ end
+end
View
25 preprocess.m
@@ -0,0 +1,25 @@
+function [tr_words te_words widx tr_idx te_idx] = preprocess(counts,varargin)
+%function [tr_words te_words widx tr_idx te_idx] = preprocess(counts,varargin)
+[debug max_words num_folds fold holdout] = process_options(varargin,'debug',0,'max-words',5000,'num-folds',5,'fold',1,'holdout',1);
+[N W] = size(counts);
+
+counts = holdoutWords(counts,holdout);
+
+te_idx = fold:num_folds:N;
+tr_idx = setdiff_sorted(1:N,te_idx);
+if debug
+ tr_idx = tr_idx(randsample(numel(tr_idx),min(1000,numel(tr_idx))));
+ te_idx = te_idx(randsample(numel(te_idx),min(100,numel(te_idx))));
+end
+count_sums = sum(counts(tr_idx,:)>0);
+max_words = min(max_words,sum(count_sums>0));
+
+widx = sortidxs(count_sums,2,'descend',max_words);
+tr_idx = tr_idx(sum(counts(tr_idx,widx),2)>0);
+te_idx = te_idx(sum(counts(te_idx,widx),2)>0);
+tr_words = [counts(tr_idx,widx)]; %sum(ap_full(tr_idx,nonwords),2)];
+te_words = [counts(te_idx,widx)]; %sum(ap_full(te_idx,nonwords),2)];
+
+%tr_words = holdoutWords(tr_words,holdout);
+%te_words = holdoutWords(te_words,holdout);
+end
View
32 runLDA.m
@@ -0,0 +1,32 @@
+function runLDA(K,seed,W,varargin)
+%function runLDA(K,seed,W,varargin)
+[debug dataset do_non_sparse do_sparse] = process_options(varargin,'debug',false,'dataset','20news','do-non-sparse',1,'do-sparse',1);
+
+directory = sprintf('traces.lda.%s',dataset);
+load(sprintf('data/%s.mat',dataset));
+
+if strcmp(dataset,'20news')
+ words = tr_data;
+ [tr_words te_words widx] = preprocess(words,'max-words',W,'debug',debug,'num-folds',10);
+else
+ [tr_words te_words widx] = preprocess(counts','max-words',W,'holdout',0.1,'debug',debug,'num-folds',50);
+ vocab = words;
+end
+
+if ~exist(directory,'dir')
+ mkdir(directory);
+end
+
+options = {'seed',seed,'te-x',te_words,'vocab',vocab(widx),'max-mstep-its',1000};
+if (debug), options = cat(2,options,'max-its',10); end
+
+basename = sprintf('%s/out.%d.%d',directory,K,seed);
+
+if do_non_sparse
+[ig theta_lda eta_lda] = sparseTAM(tr_words,K,'sparse',0,'compute-perplexity',1000,options{:});
+end
+if do_sparse
+[ig theta_sage eta_sage] = sparseTAM(tr_words,K,'sparse',1,'compute-perplexity',50,options{:});
+end
+save(sprintf('%s.final.mat',basename));
+end
View
275 sparseTAM.m
@@ -0,0 +1,275 @@
+function [acc theta eta_t eta_a] = sparseTAM(x,K,varargin)
+[sparse aspects alpha gamma_dirichlet te_x te_aspects vocab ...
+ seed hyperprior_update_interval compute_perplexity ...
+ topic_aspects verbose init_eta_t init_eta_a max_its max_mstep_its ...
+ save_prefix] = ...
+ process_options(varargin,'sparse',1,'aspects',ones(size(x,1),1),'alpha',1/K,...
+ 'init-gamma-dirichlet',1,'te-x',[],'te-aspects',[],...
+ 'vocab',[],'seed',[],'hyperprior-update-interval',1,...
+ 'compute-perplexity',1,'topic-aspects',0,'verbose',1,'init-eta-t',[],...
+ 'init-eta-a',[],'max-its',100,'max-mstep-its',100,'save-prefix',[]);
+disp(varargin);
+
+%screen words that have zero counts, because they mess stuff up
+word_counts = sum(x); x = x(:,word_counts>0);
+if ~isempty(te_x), te_x = te_x(:,word_counts>0); end
+if ~isempty(vocab), vocab = vocab(word_counts > 0); end
+
+alpha = alpha * ones(1,K); acc = 0;
+[D W] = size(x); theta = zeros(D,K); A = max(aspects);
+if A > 1, compute_perplexity = 0; end
+
+%% corpus initialization
+if ~isempty(seed), rand('seed',seed); end
+
+%mean
+m = makeLogProbs(sum(x))';
+%eta_t (topics)
+if sparse
+ if ~isempty(init_eta_t),eta_t = init_eta_t; else
+ eta_t = zeros(W,K);
+ if K>1, for k = 1:K,
+ eta_t(:,k) = computeBetaSparseVariational(full(sum(x(randsample(D,10),:))'),m,'max-its',max_mstep_its,'verbose',0);
+ end; end
+ end
+ %eta_a (aspects)
+ if ~isempty(init_eta_a), eta_a = init_eta_a; else
+ eta_a = zeros(W,A);
+ if A>1, for j = 1:A
+ if sum(aspects==j)>0, eta_a(:,j) = eta_a(:,j) + mean(x(aspects==j,:))'; end
+ end; end
+ end
+ %eta_ta (interactions)
+ eta_ta = zeros(W,K,A); if topic_aspects && A>1, for j = 1:A, for k = 1:K, eta_ta(:,k,j) = 0.1 * (eta_t(:,k) + eta_a(:,j)); end; end; end
+ %sums
+ eta_sum = makeEtaSum(m,eta_a,eta_t,eta_ta);
+else
+ assert(A==1); %can't do non-sparse aspect models now
+ for k = 1:K
+ eta_sum(:,k) = computeBetaDirichlet(full(sum(x(randsample(D,10),:))),gamma_dirichlet);
+ end
+ gamma_dirichlet = gamma_dirichlet * ones(1,W);
+end
+
+sigma = repmat(sum(x,2)/K,1,K) + repmat(alpha,D,1);
+q_a = zeros(D+1,max(aspects)); q_a(end,:) = 1/max(aspects); %prior
+
+eta_lv_score = 0; eta_ta_lv_score = 0; prior_prob = 0;
+%iter = newIterator(max_its,'debug',true,'thresh',1e-5);
+iter = newDeltaIterator(max_its,'debug',true,'thresh',1e-2);
+ecounts = zeros(W,K,A);
+while ~iter.done
+
+ %% E-step
+ word_score = 0; estep_lv_score = 0;
+ old_ecounts = ecounts;
+ ecounts = zeros(W,K,A);
+ for i = 1:D
+ if (rem(i,100)==0), fprintf('+'); end
+ old_sigma = sigma(i,:);
+ log_p_a = digamma(sum(q_a)) - digamma(sum(sum(q_a)));
+
+ [theta(i,:) q_a(i,:) new_counts sigma(i,:) score doc_lv_score doc_word_score] = tamEStep(x(i,:),eta_sum,alpha,log_p_a,old_sigma);
+ ecounts(:,:,aspects(i)) = ecounts(:,:,aspects(i)) + full(new_counts);
+ doc_lv_score = doc_lv_score + digamma(sum(q_a(:,aspects(i)))) - digamma(sum(sum(q_a))); %E[log P(a_i)] - E[log Q(a_i)]
+
+ word_score = word_score + doc_word_score;
+ estep_lv_score = estep_lv_score + doc_lv_score;
+ end
+ fprintf('fro norm of counts: %.3f\n',norm(reshape(ecounts,W,K*A)-reshape(old_ecounts,W,K*A),'fro'));
+
+ estep_lv_score = estep_lv_score - kldirichlet(sum(q_a),q_a(end,:));
+ fprintf('\n');
+ computeScore(word_score,estep_lv_score,eta_lv_score,prior_prob,'print',1);
+
+ %% M-step.
+ % we do not initialize the etas from their previous values, because
+ % otherwise we get locked into the initialization.
+ % however, this means we can't guarantee that the bound always
+ % improves. thus, we use a deltaiterator
+ if A>1 && K > 1, max_its = 1; else max_its = 0; end
+ mstep_iter = newIterator(max_its,'thresh',1e-5,'debug',true);
+ %mstep_iter = newDeltaIterator(max_its,'thresh',1e-2,'debug',true);
+ if sparse
+ eta_t_lv_score = zeros(K,1); eta_a_lv_score = zeros(A,1);
+ old_eta_t = eta_t;
+ while ~mstep_iter.done
+ if K > 1,
+ for k = 1:K
+ eq_m = logNormalizeRows(reshape(eta_sum(:,k,:),W,A)' - repmat(eta_t(:,k),1,A)');
+ eta_t(:,k) = computeBetaSparseVariational(reshape(ecounts(:,k,:),W,A),eq_m','max-its',max_mstep_its);
+ eta_sum = makeEtaSum(m,eta_a,eta_t,eta_ta);
+ end
+ end
+ fprintf(' ');
+ if A > 1,
+ for j = 1:A
+ if sparse
+ eq_m = logNormalizeRows(eta_sum(:,:,j)' - repmat(eta_a(:,j),1,K)');
+ eta_a(:,j) = computeBetaSparseVariational(ecounts(:,:,j),eq_m','max-its',max_mstep_its);
+ eta_sum = makeEtaSum(m,eta_a,eta_t,eta_ta);
+ else
+ assert(K==1);
+ eta_a(:,j) = computeBetaDirichlet(ecounts(:,1,j)');
+ end
+ end
+ end
+ fprintf(' ');
+ if topic_aspects, for k = 1:K, for j=1:A,
+ assert(sparse);
+ eq_m = logNormalizeRows(eta_sum(:,k,j)' - eta_ta(:,k,j)');
+ eta_ta(:,k,j) = computeBetaSparseVariational(ecounts(:,k,j),eq_m','max-its',max_mstep_its);
+ end; end;
+ eta_sum = makeEtaSum(m,eta_a,eta_t,eta_ta);
+ end
+ word_score = 0;
+ for j = 1:A, for k = 1:K
+ word_score = word_score + scoreWords(ecounts(:,k,j),eta_sum(:,k,j));
+ end; end
+ eta_lv_score = sum(eta_t_lv_score) + sum(eta_a_lv_score) + sum(sum(eta_ta_lv_score));
+ fprintf('\n');
+ total_score = computeScore(word_score,estep_lv_score,eta_lv_score,prior_prob,'print',0);
+ density_t = sum(sum(abs(eta_t) > 1e-5)) / numel(eta_t);
+ density_a = sum(sum(abs(eta_a) > 1e-5)) / numel(eta_a);
+ density_ta = sum(sum(sum(abs(eta_ta)>1e-5))) / numel(eta_ta);
+ fprintf('%.3f\tT=%.3f A=%.3f TA=%.3f norm diff=%.3f\n',total_score,density_t,density_a,density_ta,norm(eta_t-old_eta_t,'fro'));
+ %mstep_iter = updateIterator(mstep_iter,word_score + eta_lv_score);
+ mstep_iter = updateDeltaIterator(mstep_iter,reshape(eta_sum,W,A*K));
+ end
+ else
+ if A == 1
+ eta_sum(:,:,1) = computeBetaDirichlet(ecounts(:,:,1)',gamma_dirichlet)';
+ elseif K == 1
+ for a = 1:A
+ eta_sum(:,1,a) = computeBetaDirichlet(ecounts(:,1,a)',gamma_dirichlet);
+ end
+ else
+ error('in non-sparse model, either K or A must be 1');
+ end
+ end
+
+ %% status
+ %iter = updateIterator(iter,total_score);
+ iter = updateDeltaIterator(iter,[reshape(theta,D*K,1); reshape(q_a(1:end-1),D*A,1)]);
+ if ~isempty(vocab)
+ if sparse
+ if K > 1
+ fprintf(' ----- topics ----- \n');
+ if topic_aspects %print interaction terms
+ for k = 1:K
+ for j = 1:A
+ fprintf('K%dA%d\t ',k,j);
+ makeTopicReport(eta_sum(1:numel(vocab),k,j)',vocab,'N',10,'background',m'); %i'm not sure supplying the background is a good idea
+ end
+ fprintf('\n');
+ end
+ end
+ %same as above (for topics)
+ makeTopicReport(tprod(mean(q_a),[-3],eta_sum,[1 2 -3])',vocab,'N',20);
+ end
+ if A > 1
+ fprintf(' ----- aspects ----- \n');
+ makeTopicReport(eta_a(1:numel(vocab),:)',vocab,'N',20,'background',zeros(1,W));
+ end
+ else
+ if A > K
+ for k = 1:K,
+ if K>1, fprintf(' ----- topic %d ----- ',k); end
+ makeTopicReport(reshape(eta_sum(1:numel(vocab),k,:),W,A)',vocab,'N',20);
+ end
+ else
+ for a = 1:A
+ if A > 1, fprintf(' ----- aspect %d ----- \n', a); end
+ makeTopicReport(eta_sum(1:numel(vocab),:,a)',vocab,'N',20);
+ end
+ end
+ end
+ end
+
+ %% hyperpriors
+ prior_prob = 0;
+ if rem(iter.its,hyperprior_update_interval) == 0
+ %topics might be too sparse? consider removing this.
+ if K > 1
+ e_log_theta = digamma(sigma) - repmat(digamma(sum(sigma,2)),1,K);
+ alpha = fitDirichletPrior(e_log_theta);
+ end
+ if verbose > 0, fprintf('new alpha: %s\n',sprintf('%.2f ',alpha)); end
+ if ~sparse
+ %gamma = fitDirichletPrior(reshape(eta_sum,W,K*A)');
+ %here is a newton step that will increase your bound...
+ for i = 1:10
+ g_gamma = K * A * W * (digamma(W*gamma_dirichlet(1)) - digamma(gamma_dirichlet(1))) + sum(sum(sum(eta_sum)));
+ g_gamma = g_gamma - 2/gamma_dirichlet(1) + 1/(gamma_dirichlet(1)^2); %from IG(1,1) prior
+ h_gamma = K * A * W * (W * trigamma(W*gamma_dirichlet(1)) - trigamma(gamma_dirichlet(1)));
+ h_gamma = h_gamma + 2/(gamma_dirichlet(1)^2) - 2/(gamma_dirichlet(1)^3); %from IG(1,1) prior
+ gamma_dirichlet = gamma_dirichlet - 0.2 * g_gamma / h_gamma;
+ end
+ prior_prob = prior_prob - 2 * log(gamma_dirichlet(1)) - 1 / gamma_dirichlet(1); %from IG(1,1) prior
+ if verbose > 0, fprintf('new mean gamma = %.5f\n',mean(gamma_dirichlet)); end
+ end
+ end
+
+ %% TEST SET
+ if ~isempty(te_x)
+ theta_te = zeros(size(te_x,1),K);
+ qa_te = zeros(size(te_x,1),max(aspects));
+ e_log_a = digamma(sum(q_a)) - digamma(sum(sum(q_a)));
+
+ %use the Chib estimator from Wallach et al
+ if rem(iter.its,compute_perplexity)==0
+ if sparse, pred_topics = logToSimplex2(eta_sum');
+ else, pred_topics = logToSimplex2(computeBetaDirichlet(ecounts',gamma_dirichlet)); end
+ fprintf('PERPLEX %d: %.5f\n',iter.its,computePerplexity(te_x,pred_topics,alpha));
+ end
+ if A > 1 && rem(iter.its,5)==1
+ for i = 1:100:size(te_x,1)
+ [theta_te(i,:) qa_te(i,:)] = tamEStep(te_x(i,:),eta_sum(:,:,1:max(aspects)),alpha,e_log_a);
+ end
+ [ig preds] = max(qa_te(1:100:end,:)');
+ acc = mean(preds==te_aspects(1:100:end));
+ fprintf('ACC %d: = %.3f\n',iter.its,acc);
+ end
+ end
+ if ~isempty(save_prefix) && rem(iter.its,5)==1
+ save(sprintf('%s.%d.mat',save_prefix,iter.its));
+ end
+end
+if ~isempty(te_x)
+ if compute_perplexity
+ if sparse, pred_topics = normalize_rows(exp(eta_sum'));
+ else, pred_topics = logToSimplex2(computeBetaDirichlet(ecounts',gamma_dirichlet)); end
+ fprintf('FINAL PERPLEX\t%.5f\t%.5f\n',computePerplexity(te_x,pred_topics,alpha),iter.prev(1));
+ end
+ if A>1 && ~isempty(te_aspects )
+ for i = 1:size(te_x,1)
+ [theta_te(i,:) qa_te(i,:)] = tamEStep(te_x(i,:),eta_sum(:,:,1:max(aspects)),alpha,e_log_a);
+ end
+ [ig preds] = max(qa_te');
+ acc = mean(preds==te_aspects);
+ fprintf('FINAL ACC\t%.5f\t%.5f\n',acc,iter.prev(1));
+ end
+end
+if ~sparse %can't remember why i was doing this...
+ eta_t = eta_sum(:,:,1) - repmat(m,1,K);
+ eta_a = eta_sum(:,1,:) - repmat(m,1,A);
+end
+end
+
+function eta_sum = makeEtaSum(m,eta_a,eta_t,eta_ta)
+[W K A] = size(eta_ta);
+eta_sum = zeros(W,K,A);
+for j = 1:A
+ for k = 1:K
+ eta_sum(:,k,j) = logNormalizeVec(m + eta_a(:,j) + eta_t(:,k) + eta_ta(:,k));
+ end
+end
+end
+
+function total_score = computeScore(word_score, estep_lv_score, eta_lv_score, prior_score, varargin)
+[print] = process_options(varargin,'print',0);
+total_score = word_score + estep_lv_score + eta_lv_score + prior_score;
+if print,
+ fprintf('%.3f\t=\t%.3f\t+\t%.3f\t+\t%.3f\t+\t%.3f\n',total_score,word_score,estep_lv_score,eta_lv_score,prior_score);
+end
+end
View
4 startup.m
@@ -0,0 +1,4 @@
+addpath(genpath('3rd-party'));
+addpath(genpath('utils'));
+addpath(genpath('estimate'));
+addpath(genpath('bounds'));
View
23 tamEStep.m
@@ -0,0 +1,23 @@
+function [theta q_a out_counts sigma score lv_score word_score] = tamEStep(x,beta,alpha,a_log_prior,sigma)
+[W K A] = size(beta);
+mydudes = x>0;
+myx = x(mydudes);
+mybeta = beta(mydudes,:,:);
+
+if nargin == 4, sigma = ones(1,K) * (sum(x) / K) + alpha; end
+if isscalar(alpha) alpha = alpha * ones(1,K); end
+
+ran_mex_successfully = false;
+if exist('tamEStepInnerMex') == 3 %if we have this mex file
+ [sigma q_a ecounts score word_score lv_score] = tamEStepInnerMex(full(myx)',mybeta,alpha,a_log_prior,sigma);
+ %sometimes the mex version returns NaN, i haven't been able to debug
+ %this. so when it happens, i just back off to the non-mex version
+ ran_mex_successfully = (max(max(max(isnan(ecounts))))==0);
+end
+if ~ran_mex_successfully
+ [sigma q_a ecounts score word_score lv_score] = tamEStepNoMex(myx,mybeta,alpha,a_log_prior,sigma);
+end
+out_counts = zeros(W,K,A);
+out_counts(mydudes,:,:) = ecounts;
+theta = normalize_rows(sigma);
+end
View
272 tamEStepInnerMex.c
@@ -0,0 +1,272 @@
+#include <mex.h>
+#include <assert.h>
+#include <math.h>
+#include "util.h"
+
+/*function [theta q_a score lv_score word_score sigma] = tamEStep(x,beta,alpha,a_log_prior,sigma)
+ my_counts = [N x 1] term counts for all unique terms in the doc
+ my_beta = [N x K x A] log topic-word distribs for all unique terms in the doc
+ alpha = [1 x K] doc-topic prior
+ a_log_prior = [1 x A] doc-aspect prior
+
+ output:
+ sigma = [1 x K] variational posterior distrib over theta
+ q_a = [1 x A] variational posterior distrib over aspect
+ phi = [N x K] variational posterior over z (per token topic)
+
+ todo: include additional terms from KL(P(theta) || Q(theta)) that involve
+ alpha, so that you can do alpha updating. also, optimize further using fast exp or something.
+*/
+
+#define THRESH 0.00001
+#define MAX_ITS 25
+#define FAST_EXP 0
+
+int DEBUG=0;
+
+static union
+{
+ double d;
+ struct {
+#ifdef LITTLE_ENDIAN
+ int j,i;
+#else
+ int i,j;
+#endif
+ } n;
+} _eco;
+
+#define EXP_A (1048576/0.69314718055994530942)
+#define EXP_C 60801
+#define EXP(y) (_eco.n.i = EXP_A*(y) + (1072693248 - EXP_C), _eco.d)
+
+void my_assert(int condition){
+ if (!condition)
+ mexErrMsgTxt("assertion failed");
+}
+
+void logToSimplex4(const double* source, double* target, int K, int fast_exp){
+ double log_sum_exp,max_val;
+ int k;
+ max_val = source[0];
+ for (k = 1; k < K; k++)
+ if (source[k] > max_val) max_val = source[k];
+
+ log_sum_exp = 0;
+ for (k = 0; k < K; k++){
+ if (fast_exp) log_sum_exp += EXP(source[k] - max_val);
+ else log_sum_exp += exp(source[k] - max_val);
+ }
+ log_sum_exp = log(log_sum_exp);
+ for (k = 0; k < K; k++){
+ if (fast_exp) target[k] = EXP(source[k] - max_val - log_sum_exp);
+ else target[k] = exp(source[k] - log_sum_exp - max_val);
+ }
+}
+
+void logToSimplex(const double* source, double* target, int K){
+ logToSimplex4(source,target,K,FAST_EXP);
+}
+
+int getIdx(int i, int j, int k, int I, int J, int K){
+ return i + j * I + k * I * J;
+}
+
+double get(double* arr, int i, int j, int k, int I, int J, int K){
+ /* if (i % 10 == 2 && j % 10 ==3 && k % 10 == 4)*/
+ /* printf("%.3f %i %i/%i %i/%i %i/%i\n",arr[idx],idx,i,I,j,J,k,K); */
+ return arr[getIdx(i,j,k,I,J,K)];
+}
+
+double set(double* arr, double val, int i, int j, int k, int I, int J, int K){
+ arr[getIdx(i,j,k,I,J,K)] = val;
+}
+
+/*
+ function [sigma q_a out_counts score lv_score word_score]
+ = tamEStep(x,beta,alpha,a_log_prior,sigma)
+*/
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[]
+ )
+{
+ if (nrhs < 4 || nrhs > 5) {
+ mexErrMsgTxt("function [sigma q_a phi] = tamEStepInner(c,beta,alpha,a_log_prior,[sigma])");
+ return;
+ }
+ mwSize N,A,K;
+ int n,a,k,it;
+ double word_score, lv_score, log_p_k_theta_denom, log_p_w_a, old_score, change;
+ mwSize* beta_dims;
+ beta_dims = mxGetDimensions(prhs[1]);
+
+ N = mxGetM(prhs[0]);
+
+ /* if (N == 399)
+ DEBUG=1;
+ else
+ DEBUG=0;
+*/
+
+ K = beta_dims[1];
+
+ if (mxGetNumberOfDimensions(prhs[1]) > 2)
+ A = beta_dims[2];
+ else A = 1;
+
+ mxAssert (beta_dims[0] == N,"beta dims must match counts");
+ mxAssert (K == mxGetN(prhs[2]),"beta dims must match K");
+ mxAssert (A == mxGetN(prhs[3]),"beta dims must match A");
+ if (DEBUG)
+ printf("dimensions N=%i K=%i A=%i\n",N,K,A);
+
+ double* counts = mxGetPr(prhs[0]);
+ double* beta = mxGetPr(prhs[1]);
+ double* alpha = mxGetPr(prhs[2]);
+ double* log_p_a = mxGetPr(prhs[3]);
+
+ double* sigma;
+ if (nrhs == 5) {
+ sigma = mxGetPr(prhs[4]);
+ if (DEBUG){
+ printf("reading sigma:\n");
+ for (k = 0; k < K; k++){
+ printf("%.3f ",sigma[k]);
+ }
+ printf("\n");
+ }
+ }
+ else {
+ sigma = mxCalloc(K,sizeof(double));
+ for (k = 0; k < K; k++){
+ sigma[k] = alpha[k];
+ for (n = 0; n < N; n++){
+ sigma[k] += counts[n] / (double) K;
+ }
+ }
+ }
+
+ double* sigma_buffer = mxCalloc(K,sizeof(double));
+ double* phi = mxCalloc(K,sizeof(double));
+ double* phi_buffer = mxCalloc(K,sizeof(double));
+ double* q_a = mxCalloc(A,sizeof(double));
+ double* a_buffer = mxCalloc(A,sizeof(double));
+ double* log_p_k_theta = mxCalloc(K,sizeof(double));
+ double* out; /* for output counts */
+
+ logToSimplex(log_p_a,q_a,A);
+
+ old_score = 0;
+ for (it = 0; it < MAX_ITS; it++){
+ word_score = 0.0; lv_score = 0.0;
+
+ for (k = 0; k < K; k++) {
+ sigma_buffer[k] = alpha[k]; /*this is for storing sigma later */
+ log_p_k_theta_denom += sigma[k];
+ }
+ log_p_k_theta_denom = digamma(log_p_k_theta_denom);
+
+ /* parts of KL( P(theta) || Q(theta) )
+ Note we are leaving out constant parts:
+ -gammaln(sum(sigma)) + gammaln(sum(alpha) - sum(gammaln(alpha))
+ */
+ for (k = 0; k < K; k++){
+ log_p_k_theta[k] = digamma(sigma[k]) - log_p_k_theta_denom;
+ lv_score += log_p_k_theta[k] * (alpha[k] - sigma[k]);
+ lv_score += gammaln(sigma[k]);
+ }
+ for (a = 0; a < A; a++){
+ if (q_a[a] > 0)
+ lv_score += q_a[a] * (log_p_a[a] - log(q_a[a])); /* E[log P(a | var_theta)] - E[log q(a)]*/
+ a_buffer[a] = log_p_a[a]; /* for updating q(a) */
+ }
+
+ for (n = 0; n < N; n++){
+ /* compute phi */
+ for (k = 0; k < K; k++){
+ phi_buffer[k] = log_p_k_theta[k];
+ for (a = 0; a < A; a++){
+ phi_buffer[k] += q_a[a] * get(beta,n,k,a,N,K,A);
+ }
+ }
+ logToSimplex(phi_buffer,phi,K);
+
+ /* E[log p(z | theta)] - E[log Q(z)] */
+ for (k = 0; k < K; k++){
+ if (phi[k] > 0){
+ lv_score += counts[n] * phi[k] * (log_p_k_theta[k] - log(phi[k]));
+ phi[k] *= counts[n]; /* this is for efficiency later */
+ sigma_buffer[k] += phi[k];
+ for (a = 0; a < A; a++){
+ log_p_w_a = get(beta,n,k,a,N,K,A);
+ a_buffer[a] += phi[k] * log_p_w_a;
+ word_score += q_a[a] * phi[k] * log_p_w_a;
+ }
+ }
+ }
+ }
+ logToSimplex(a_buffer,q_a,A);
+ for (k = 0; k < K; k++) sigma[k] = sigma_buffer[k];
+
+ /* print status */
+ if (DEBUG){
+ for (a = 0; a < A; a++) printf("%.3e ",log(q_a[a]));
+ printf("\n");
+ for (k = 0; k < K; k++) printf("%.3f ",sigma[k]);
+ printf("\n");
+ }
+
+ change = (word_score + lv_score - old_score) / abs(old_score);
+ if (DEBUG) printf("Iteration %i (%.3e): %.6e = %.3e + %.3e\n",it,change,word_score + lv_score,word_score,lv_score);
+ if (change < THRESH && it > 0){
+ break;
+ }
+ old_score = word_score + lv_score;
+ }
+
+ /************************* OUTPUT ********************/
+ /* theta */
+ if (nlhs > 0){
+ plhs[0] = mxCreateDoubleMatrix(1,K,mxREAL);
+ out = mxGetPr(plhs[0]);
+ for (k = 0; k < K; k++) out[k] = sigma[k];
+ }
+ if (nlhs > 1){
+ plhs[1] = mxCreateDoubleMatrix(1,A,mxREAL);
+ out = mxGetPr(plhs[1]);
+ for (a = 0; a < A; a++) out[a] = q_a[a];
+ }
+
+ if (nlhs > 2){
+ int beta_dims_out[3];
+
+ beta_dims_out[0] = N;
+ beta_dims_out[1] = K;
+ beta_dims_out[2] = A;
+
+ plhs[2] = mxCreateNumericArray(3,beta_dims_out,mxDOUBLE_CLASS,mxREAL);
+ out = mxGetPr(plhs[2]);
+ for (n = 0; n < N; n++){
+ for (k = 0; k < K; k++){
+ phi_buffer[k] = log_p_k_theta[k];
+ for (a = 0; a < A; a++)
+ phi_buffer[k] += q_a[a] * get(beta,n,k,a,N,K,A);
+ }
+ logToSimplex(phi_buffer,phi,K);
+ for (k = 0; k < K; k++){
+ for (a = 0; a < A; a++)
+ set(out,counts[n] * q_a[a] * phi[k],n,k,a,N,K,A);
+ }
+ }
+ }
+ if (nlhs > 3)
+ plhs[3] = mxCreateDoubleScalar(word_score + lv_score);
+ if (nlhs > 4)
+ plhs[4] = mxCreateDoubleScalar(word_score);
+ if (nlhs > 5)
+ plhs[5] = mxCreateDoubleScalar(lv_score);
+
+}
+
View
54 tamEStepNoMex.m
@@ -0,0 +1,54 @@
+function [sigma q_a out_counts score lv_score word_score sigma] = tamEStepNoMex(x,beta,alpha,a_log_prior,sigma)
+%function [theta q_a out_counts score lv_score word_score sigma] = tamEStepNoMex(x,beta,alpha,a_log_prior,sigma)
+%
+%this contains a lot of optimizations that make it kind of complicated.
+%maybe i should release a simplified version?
+[W K A] = size(beta);
+if nargin < 5
+ if isscalar(alpha), sigma = ones(1,K) * (numel(x) / K + alpha);
+ else sigma = ones(1,K) * numel(x) / K + alpha; end
+end
+mydudes = x>0;
+q_a = normalize_vec(ones(1,A));
+myphi = 1/K * ones(sum(mydudes),K);
+myx = x(mydudes);
+mybeta = beta(mydudes,:,:);
+
+%new_counts = spalloc(K,W,sum(mydudes)*K);
+xrep = full(repmat(myx',1,K));
+iter = newIterator(5,'thresh',1e-4,'debug',false);
+%iter_delta = newDeltaIterator(1000,'thresh',1e-3);
+iter_ctr = 0;
+while ~iter.done
+ sigma = (myx * myphi + alpha);
+ dig_sig = digamma(sigma) - digamma(sum(sigma));
+ warning off;
+ myphi = logToSimplex2(tprod(mybeta,[1 2 -1],q_a,[3 -1],'n') + repmat(dig_sig,size(myx,2),1));
+ [q_a log_q_a] = logToSimplex2(squeeze(tprod(mybeta,[-1 -2 1],myphi,[-1 -2],'n'))' + a_log_prior);
+ warning on;
+ %score it. this is slow, so we do it less often.
+ if rem(iter_ctr,5) == 0
+ word_score = 0;
+ %new_counts(:,mydudes) = transpose(myphi.*xrep);
+ my_new_counts = transpose(myphi.*xrep);
+ %for a = 1:A
+ %new_counts = new_counts * q_a(a);
+ % word_score = word_score + q_a(a) * scoreWords(my_new_counts',mybeta(:,:,a));
+ %end
+ %same as above, but faster
+ word_score = q_a * tprod(my_new_counts,[-1 -2],mybeta,[-2 -1 1]);
+ %using mybeta(:,:,1) looks weird, but actually this part is being ignored, we're only using the latent variable score. so this is legit.
+ [ig ig2 lv_score] = scoreDoc(my_new_counts',mybeta(:,:,1),myphi,myx,sigma,alpha,dig_sig);
+ lv_score = lv_score + q_a * (a_log_prior - log_q_a)';
+ score = word_score + lv_score;
+
+ iter = updateIterator(iter,score);
+ end
+ iter_ctr = iter_ctr + 1;
+end
+%compute output counts for real
+if nargout >= 3
+ out_counts = zeros(W,K,A);
+ for k=1:K, for a=1:A, out_counts(mydudes,k,a) = myphi(:,k) .* myx' .* q_a(a); end; end
+end
+end
View
22 utils/bayesian/fitDirichletPrior.m
@@ -0,0 +1,22 @@
+function [ alpha ] = fitDirichletPrior( e_log_theta, alpha )
+%function [ alpha ] = fitDirichletPrior( e_log_theta, alpha )
+%maximize E[log P(theta | alpha)], for
+% theta | alpha ~ Dirichlet(alpha)
+
+if ~exist('alpha','var')
+ alpha = normalize_vec(sum(exp(e_log_theta)));
+end
+[solution fX] = minimize(log(alpha)',@ll_log_alpha,100,e_log_theta);
+alpha = exp(solution)';
+
+ function [l g] = ll_log_alpha(log_alpha,e_log_theta)
+ N = size(e_log_theta,1);
+ alpha = exp(log_alpha');
+ l = N * (gammaln(sum(alpha)) - sum(gammaln(alpha))) + (alpha - 1) * sum(e_log_theta)';
+ g = N * (digamma(sum(alpha)) - digamma(alpha)) + sum(e_log_theta);
+ g = transpose(-g .* alpha);
+ l = -l;
+ end
+
+end
+
View
5 utils/bayesian/kldirichlet.m
@@ -0,0 +1,5 @@
+function y = kldirichlet(p,q)
+ sumq = sum(q); sump = sum(p);
+ y = gammaln(sump) - gammaln(sumq) + sum(gammaln(q)-gammaln(p)) + ...
+ (p - q) * (digamma(p) - digamma(sump))';
+end
View
15 utils/holdoutWords.m
@@ -0,0 +1,15 @@
+function [train test] = holdoutWords(data,proportion)
+%function [train test] = holdoutWords(data)
+%takes sparse term-count representation, holds out 50% of words for training, 50% for test
+
+if nargin == 1, proportion = 0.5; end
+
+[D W] = size(data);
+
+[rows cols vals] = find(data);
+beta_param = proportion * ones(size(vals));
+proportions = randbeta(beta_param, 1- beta_param);
+train = sparse(rows,cols,round(vals .* proportions));
+test = sparse(rows,cols,round(vals .* (1-proportions)));
+
+end
View
7 utils/iterator/newDeltaIterator.m
@@ -0,0 +1,7 @@
+function iter = newDeltaIterator(max_its,varargin)
+ iter.max_its = max_its;
+ [iter.thresh iter.min_its iter.debug] = process_options(varargin,'thresh',1e-6,'min_its',0,'debug',false);
+ iter.prev = [];
+ iter.done = false;
+ iter.its = 0;
+end
View
10 utils/iterator/newIterator.m
@@ -0,0 +1,10 @@
+function iter = newIterator(max_its,varargin)
+%function iter = newIterator(max_its,varargin)
+%
+% create new iterator. assumes the value is always increasing
+ iter.max_its = max_its;
+ [iter.thresh iter.min_its iter.debug] = process_options(varargin,'thresh',0,'min-its',0,'debug',false);
+ iter.prev = -inf;
+ iter.done = false;
+ iter.its = 0;
+end
View
13 utils/iterator/updateDeltaIterator.m
@@ -0,0 +1,13 @@
+function iter = updateDeltaIterator(iter,val)
+ iter.its = iter.its + 1;
+ if iter.its > iter.max_its, iter.done = true; end
+ if isempty(iter.prev), iter.prev = val; return; end
+ delta_norm = norm(iter.prev - val,'fro');
+ if iter.debug,
+ fprintf('%d. %.3e\n',iter.its,delta_norm);
+ end
+ if delta_norm < iter.thresh
+ iter.done = true;
+ end
+ iter.prev = val;
+end
View
16 utils/iterator/updateIterator.m
@@ -0,0 +1,16 @@
+function iter = updateIterator( iter, val )
+%function [ iter ] = updateIterator( iter, val )
+iter.its = iter.its + 1;
+
+if isnan(val)
+ error('failure!');
+end
+
+change = (val - iter.prev) / abs(iter.prev);
+iter.prev = val;
+if iter.debug, fprintf('it: %d\tscore: %.3f\tchange: %.7f\n',iter.its,val,change); end
+if (iter.its > iter.max_its || (iter.its > iter.min_its && change < iter.thresh))
+ iter.done = true;
+end
+end
+
View
7 utils/logNormalizeRows.m
@@ -0,0 +1,7 @@
+function logp = logNormalizeRows(logp)
+%function logp = logNormalizeRows(logp)
+%logp_old = log(normalize_rows(exp(logp)));
+for k = 1:size(logp,1),
+ logp(k,:) = logp(k,:) - logsumexp(logp(k,:)');
+end
+end
View
4 utils/logNormalizeVec.m
@@ -0,0 +1,4 @@
+function logp = logNormalizeVec(logp)
+%function logp = logNormalizeRows(logp)
+logp = log(normalize_vec(exp(logp)));
+end
View
15 utils/logToSimplex2.m
@@ -0,0 +1,15 @@
+function [phi logphi] = logToSimplex2(logphi,minlog)
+%function [phi logphi] = logToSimplex2(logphi,minlog)
+%given log(x), return exp(log(x)) / sum(exp(log(x))).
+warning off
+lognormalizer = logsumexp(logphi,2);
+warning on
+logphi = logphi - repmat(lognormalizer,1,size(logphi,2));
+phi = exp(logphi);
+
+%sketchy... is this getting used?
+if nargin > 1,
+ logphi(isinf(logphi)) = minlog;
+end
+
+end
View
8 utils/makeLogProbs.m
@@ -0,0 +1,8 @@
+function logProbs = makeLogProbs(counts,varargin)
+%function logProbs = makeLogProbs(counts,varargin)
+%each row is a set of counts. convert to normalized log probabilities.
+[min_log_prob] = process_options(varargin,'min-log-prob',-200);
+logProbs = log(full(normalize_rows(counts)));
+logProbs(logProbs < min_log_prob) = min_log_prob;
+logProbs = log(normalize_rows(exp(logProbs)));
+end
View
27 utils/newton.m
@@ -0,0 +1,27 @@
+function [x fx] = newton(func,init,arguments,varargin)
+%function [x fx] = newton(func,init,arguments,varargin)
+% do a newton optimization
+% func must return the likelihood and the value of the update -H^{-1}g
+
+[max_its init_alpha thresh debug min_alpha] = process_options(varargin,'max-its',100,'alpha',0.1,'thresh',1e-4,'debug',0,'min-alpha',1e-5);
+its = 1;
+x = init;
+alpha = init_alpha;
+while its < max_its
+ [prev_score gradient step] = feval(func,x,arguments{:});
+ new_score = feval(func,x + alpha * step,arguments{:});
+ if debug, fprintf('%d: %.3f (alpha = %.3f)\n',its,new_score,alpha); end
+ if new_score < prev_score, alpha = alpha * 1.5; end
+ while (new_score > prev_score && alpha > min_alpha)
+ alpha = 0.5 * alpha;
+ new_score = feval(func,x + alpha * step,arguments{:});
+ if debug, fprintf('%d: %.3f (alpha = %.3f)\n',its,new_score,alpha); end
+ end
+ x = x + alpha * step;
+ fx(its) = new_score;
+ its = its + 1;
+ if prev_score - new_score < thresh
+ break;
+ end
+end
+end
View
10 utils/normalize_rows.m
@@ -0,0 +1,10 @@
+function out = normalize_rows(in)
+ out = scale_rows(in,1./row_sum(in));
+% rowSums = sum(x,2);
+% [I,J] = size(x);
+% p = zeros(I,J);
+% for i=1:I
+% p(i,:) = x(i,:) / rowSums(i);
+% end
+%
+end
View
3  utils/normalize_vec.m
@@ -0,0 +1,3 @@
+function [ p ] = normalize_vec( x )
+p = x / sum(x);
+end
View
8 utils/randsample.m
@@ -0,0 +1,8 @@
+function y = randsample(n,k)
+%function y = randsample(n,k)
+%
+%sample without replacement k integers from 1:n
+
+[ig idx] = sort(rand(n,1));
+y = idx(1:k);
+end
View
15 utils/sortidxs.m
@@ -0,0 +1,15 @@
+function [ idx val ] = sortidxs( x, dim, mode, topk )
+%function [ idx val ] = sortidxs( x, dim, mode, topk )
+if nargin < 2,
+ dim = 2;
+end
+if nargin < 3
+ mode = 'descend';
+end
+[val idx] = sort(x,dim,mode);
+if nargin == 4
+ idx = idx(1:topk);
+ val = val(1:topk);
+end
+end
+
View
3  utils/spdiag.m
@@ -0,0 +1,3 @@
+function y = spdiag(x)
+y = sparse(1:numel(x),1:numel(x),x);
+end
View
11 utils/termCountsToWordList.m
@@ -0,0 +1,11 @@
+function list = termCountsToWordList(term_counts)
+%function word_list = termCountsToWordList(term_counts)
+% takes a (possibly sparse) term counts matrix and returns a word list
+% calls a cute mex function
+for i = 1:size(term_counts,1)
+ [ig j s] = find(term_counts(i,:));
+ list{i} = termCountsToWordListMex(j,s);
+end
+if size(term_counts,1) == 1, list = list{1}; end
+end
+
View
36 utils/termCountsToWordListMex.c
@@ -0,0 +1,36 @@
+#include "mex.h"
+
+/**
+ here's how this works
+ prhs[0] = columns [row vector]
+ prhs[1] = counts [row vector]
+
+**/
+
+void mexFunction(
+ int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[]
+ )
+{
+ if (nrhs != 2 || nlhs != 1){ mexErrMsgTxt("out = func(cols,counts)"); }
+
+ double* cols, *counts, *out;
+ int i,j,ctr,tot_counts,num_cols;
+
+ tot_counts = 0;
+ cols = mxGetPr(prhs[0]);
+ counts = mxGetPr(prhs[1]);
+ num_cols = mxGetN(prhs[0]);
+
+ for (i = 0; i < num_cols; i++) tot_counts += counts[i];
+
+ plhs[0] = mxCreateDoubleMatrix(1,tot_counts,mxREAL);
+ out = mxGetPr(plhs[0]);
+
+ ctr = 0;
+ for (i = 0; i < num_cols ;i++){
+ for (j = 0; j < counts[i]; j++){
+ out[ctr++] = cols[i];
+ }
+ }
+}
View
BIN  utils/termCountsToWordListMex.mexglx
Binary file not shown
Please sign in to comment.
Something went wrong with that request. Please try again.