Initial python files, non-working

runstuff_BNPcovreg.m → matlab/runstuff_BNPcovreg.m

@@ -1,7 +1,7 @@
 %% SCRIPT FOR RUNNING BASIC SAMPLER ON SIMULATED DATA EXAMPLE %%
 
 % Load simulated data:
-load simData.mat
+load simData
 
 % y is a p x N matrix specifying the N observations of the p-dimensional data
 [p N] = size(y);

python/BNP_covreg_varinds.py

@@ -0,0 +1,324 @@
+# Autogenerated with SMOP version 
+# /usr/local/bin/smop BNP_covreg_varinds.m -o ../python/BNP_covreg_varinds.py
+
+from __future__ import division
+try:
+    from runtime import *
+except ImportError:
+    from smop.runtime import *
+
+def BNP_covreg_varinds_(y=None,prior_params=None,settings=None,restart=None,true_params=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 5-[y,prior_params,settings,restart,true_params].count(None)+len(args)
+
+    p,N=size_(y,nargout=2)
+    if exist_(char('true_params'),char('var')):
+        cov_true_diag=zeros_(p,N)
+        for tt in arange_(1,N).reshape(-1):
+            cov_true_diag[:,tt]=diag_(true_params.cov_true(arange_(),arange_(),tt))
+    if not exist_(char('restart'),char('var')):
+        restart=0
+    sample_K_flag=settings.sample_K_flag
+    latent_mean=settings.latent_mean
+    inds_y=settings.inds_y
+    y[not inds_y]=0
+    k=settings.k
+    L=settings.L
+    Niter=settings.Niter
+    trial=settings.trial
+    if not restart:
+        Stats[1:settings.saveEvery / settings.storeEvery]=struct_(char('zeta'),zeros_(L,k,N),char('psi'),zeros_(k,N),char('invSig_vec'),zeros_(1,p),char('theta'),zeros_(p,L),char('eta'),zeros_(k,N),char('phi'),zeros_(p,L),char('tau'),zeros_(1,L),char('K_ind'),0)
+        store_counter=1
+        delta=zeros_(1,L)
+        delta[1]=gamrnd_(prior_params.hypers.a1,1)
+        delta[2:L]=gamrnd_(prior_params.hypers.a2 * ones_(1,L - 1),1)
+        tau=exp_(cumsum_(log_(delta)))
+        phi=gamrnd_(prior_params.hypers.a_phi * ones_(p,L),1) / prior_params.hypers.b_phi
+        theta=zeros_(p,L)
+        for pp in arange_(1,p).reshape(-1):
+            theta[pp,:]=chol_(diag_(1.0 / (phi[pp,:].dot(tau)))).T * randn_(L,1)
+        xi=randn_(k,N)
+        psi=zeros_(k,N)
+        eta=psi + xi
+        invSig_vec=gamrnd_(prior_params.sig.a_sig * ones_(1,p),1) / prior_params.sig.b_sig
+        if sample_K_flag == 1 or sample_K_flag == 2:
+            Pk=cumsum_(prior_params.K.c_prior)
+            K_ind=1 + sum_(Pk[end()] * rand_(1) > Pk)
+        else:
+            K_ind=1
+        zeta=zeros_(L,k,N)
+        zeta=sample_zeta_(y,theta,eta,invSig_vec,zeta,prior_params.K.invK(arange_(),arange_(),K_ind),inds_y)
+        if not exist_(settings.saveDir,char('file')):
+            mkdir_(settings.saveDir)
+        if isfield_(settings,char('filename')):
+            settings_filename=strcat_(settings.saveDir,char('/'),settings.filename,char('_info4trial'),num2str_(trial))
+            init_stats_filename=strcat_(settings.saveDir,char('/'),settings.filename,char('initialStats_trial'),num2str_(trial))
+        else:
+            settings_filename=strcat_(settings.saveDir,char('/info4trial'),num2str_(trial))
+            init_stats_filename=strcat_(settings.saveDir,char('/initialStats_trial'),num2str_(trial))
+        if nargin > 4:
+            save_(settings_filename,char('y'),char('settings'),char('prior_params'),char('true_params'))
+        else:
+            save_(settings_filename,char('y'),char('settings'),char('prior_params'))
+        save_(init_stats_filename,char('zeta'),char('psi'),char('eta'),char('theta'),char('invSig_vec'),char('phi'),char('tau'),char('K_ind'))
+        num_iters=1
+        nstart=1
+    else:
+        load_([settings.saveDir,char('/BNP_covreg_stats'),char('iter'),num2str_(settings.lastIter),char('trial'),num2str_(settings.trial)])
+        theta=Stats_(end()).theta
+        eta=Stats_(end()).eta
+        zeta=Stats_(end()).zeta
+        phi=Stats_(end()).phi
+        tau=Stats_(end()).tau
+        psi=Stats_(end()).psi
+        nstart=settings.lastIter + 1
+        num_iters=1
+        store_counter=1
+    for nn in arange_(nstart,Niter).reshape(-1):
+        invSig_vec=sample_sig_(y,theta,eta,zeta,prior_params.sig,inds_y)
+        phi,tau=sample_hypers_(theta,phi,tau,prior_params.hypers,nargout=2)
+        theta=sample_theta_(y,eta,invSig_vec,zeta,phi,tau,inds_y)
+        if sample_K_flag == 1:
+            error_(char('need to code up K sampling for inds_y'))
+        else:
+            if sample_K_flag == 2:
+                error_(char('need to code up K sampling for inds_y'))
+            else:
+                K_ind=1
+        if latent_mean:
+            psi=sample_psi_margxi_(y,theta,invSig_vec,zeta,psi,prior_params.K.invK(arange_(),arange_(),K_ind),inds_y)
+        else:
+            psi=zeros_(size_(xi))
+        xi=sample_xi_(y,theta,invSig_vec,zeta,psi,inds_y)
+        eta=psi + xi
+        for ii in arange_(1,num_iters).reshape(-1):
+            zeta=sample_zeta_(y,theta,eta,invSig_vec,zeta,prior_params.K.invK(arange_(),arange_(),K_ind),inds_y)
+        num_iters=1
+        if rem_(nn,settings.storeEvery) == 0 and nn >= settings.saveMin:
+            Stats[store_counter].zeta=zeta
+            Stats[store_counter].eta=eta
+            Stats[store_counter].theta=theta
+            Stats[store_counter].invSig_vec=invSig_vec
+            Stats[store_counter].psi=psi
+            Stats[store_counter].phi=phi
+            Stats[store_counter].tau=tau
+            Stats[store_counter].K_ind=K_ind
+            store_counter=store_counter + 1
+        if rem_(nn,settings.saveEvery) == 0:
+            if isfield_(settings,char('filename')):
+                filename=strcat_(settings.saveDir,char('/'),settings.filename,char('iter'),num2str_(nn),char('trial'),num2str_(settings.trial))
+            else:
+                filename=strcat_(settings.saveDir,char('/BNP_covreg_stats'),char('iter'),num2str_(nn),char('trial'),num2str_(settings.trial))
+            save_(filename,char('Stats'))
+            store_counter=1
+        if not rem_(nn,100):
+            display_([char('Iter: '),num2str_(nn),char(', K_ind: '),num2str_(K_ind)])
+            if exist_(char('true_params'),char('var')):
+                cov_est=zeros_(p,N)
+                for tt in arange_(1,N).reshape(-1):
+                    cov_est[:,tt]=diag_(theta * zeta[:,:,tt] * zeta[:,:,tt].T * theta.T + diag_(1.0 / invSig_vec))
+                plot_(cov_true_diag.T,char('LineWidth'),2)
+                hold_(char('on'))
+                plot_(cov_est.T,char('--'),char('LineWidth'),2)
+                hold_(char('off'))
+                drawnow
+    return
+def sample_zeta_(y=None,theta=None,eta=None,invSig_vec=None,zeta=None,invK=None,inds_y=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 7-[y,theta,eta,invSig_vec,zeta,invK,inds_y].count(None)+len(args)
+
+    p,L=size_(theta,nargout=2)
+    k,N=size_(eta,nargout=2)
+    mu_tot=zeros_(p,N)
+    error_nn=zeros_(L,N)
+    for nn in arange_(1,N).reshape(-1):
+        mu_tot[:,nn]=theta * zeta[:,:,nn] * eta[:,nn]
+        error_nn[:,nn]=(theta ** 2).T * (invSig_vec.T.dot((1 - inds_y[:,nn])))
+    numiter=1
+    for nn in arange_(1,numiter).reshape(-1):
+        for ll in arange_(1,L).reshape(-1):
+            theta_ll=theta[:,ll]
+            for kk in randperm_(k).reshape(-1):
+                eta_kk=eta[kk,:]
+                zeta_ll_kk=squeeze_(zeta[ll,kk,:]).T
+                mu_tot=mu_tot - theta_ll[:,ones_(1,N)].dot(eta_kk[ones_(p,1),:]).dot(zeta_ll_kk[ones_(p,1),:])
+                A_lk_invSig_A_lk=(eta[kk,:] ** 2).dot(((theta[:,ll] ** 2).T * invSig_vec.T - error_nn[ll,:]))
+                theta_tmp=theta[:,ll].T.dot(invSig_vec)
+                ytilde=(y - mu_tot).dot(inds_y)
+                theta_lk=eta[kk,:].T.dot((theta_tmp * ytilde).T)
+                cholSig_lk_trans=numpy.linalg.solve(chol_(invK + diag_(A_lk_invSig_A_lk)),eye_(N))
+                m_lk=cholSig_lk_trans * (cholSig_lk_trans.T * theta_lk)
+                zeta[ll,kk,:]=m_lk + cholSig_lk_trans * randn_(N,1)
+                zeta_ll_kk=squeeze_(zeta[ll,kk,:]).T
+                mu_tot=mu_tot + theta_ll[:,ones_(1,N)].dot(eta_kk[ones_(p,1),:]).dot(zeta_ll_kk[ones_(p,1),:])
+    return zeta
+def sample_psi_margxi_(y=None,theta=None,invSig_vec=None,zeta=None,psi=None,invK=None,inds_y=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 7-[y,theta,invSig_vec,zeta,psi,invK,inds_y].count(None)+len(args)
+
+    p,L=size_(theta,nargout=2)
+    k,N=size_(psi,nargout=2)
+    Sigma_0=diag_(1.0 / invSig_vec)
+    mu_tot=zeros_(p,N)
+    Omega=zeros_(p,k,N)
+    OmegaInvOmegaOmegaSigma0=zeros_(k,p,N)
+    for nn in arange_(1,N).reshape(-1):
+        Omega[inds_y[:,nn],:,nn]=theta[inds_y[:,nn],:] * zeta[:,:,nn]
+        temp=numpy.linalg.solve((Omega[inds_y[:,nn],:,nn] * Omega[inds_y[:,nn],:,nn].T + Sigma_0[inds_y[:,nn],inds_y[:,nn]]),eye_(sum_(inds_y[:,nn])))
+        OmegaInvOmegaOmegaSigma0[:,inds_y[:,nn],nn]=Omega[inds_y[:,nn],:,nn].T * temp
+        mu_tot[:,nn]=Omega[:,:,nn] * psi[:,nn]
+    if sum_(sum_(mu_tot)) == 0:
+        numTotIters=50
+    else:
+        numTotIters=5
+    for numIter in arange_(1,numTotIters).reshape(-1):
+        for kk in randperm_(k).reshape(-1):
+            Omega_kk=squeeze_(Omega[:,kk,:])
+            psi_kk=psi[kk,:]
+            mu_tot=mu_tot - Omega_kk.dot(psi_kk[ones_(p,1),:])
+            theta_k=diag_(squeeze_(OmegaInvOmegaOmegaSigma0[kk,:,:]).T * (y - mu_tot))
+            Ak_invSig_Ak=diag_(squeeze_(OmegaInvOmegaOmegaSigma0[kk,:,:]).T * Omega_kk)
+            cholSig_k_trans=numpy.linalg.solve(chol_(invK + diag_(Ak_invSig_Ak)),eye_(N))
+            m_k=cholSig_k_trans * (cholSig_k_trans.T * theta_k)
+            psi[kk,:]=m_k + cholSig_k_trans * randn_(N,1)
+            psi_kk=psi[kk,:]
+            mu_tot=mu_tot + Omega_kk.dot(psi_kk[ones_(p,1),:])
+    return psi
+def sample_xi_(y=None,theta=None,invSig_vec=None,zeta=None,psi=None,inds_y=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 6-[y,theta,invSig_vec,zeta,psi,inds_y].count(None)+len(args)
+
+    p,N=size_(y,nargout=2)
+    L,k=size_(zeta[:,:,1],nargout=2)
+    xi=zeros_(k,N)
+    for nn in arange_(1,N).reshape(-1):
+        theta_zeta_n=theta * zeta[:,:,nn]
+        y_tilde_n=y[:,nn] - theta_zeta_n * psi[:,nn]
+        invSigMat=invSig_vec.dot(inds_y[:,nn].T)
+        invSigMat=invSigMat[ones_(k,1),:]
+        zeta_theta_invSig=theta_zeta_n.T.dot(invSigMat)
+        cholSig_xi_n_trans=numpy.linalg.solve(chol_(eye_(k) + zeta_theta_invSig * theta_zeta_n),eye_(k))
+        m_xi_n=cholSig_xi_n_trans * (cholSig_xi_n_trans.T * (zeta_theta_invSig * y_tilde_n))
+        xi[:,nn]=m_xi_n + cholSig_xi_n_trans * randn_(k,1)
+    return xi
+def sample_theta_(y=None,eta=None,invSig_vec=None,zeta=None,phi=None,tau=None,inds_y=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 7-[y,eta,invSig_vec,zeta,phi,tau,inds_y].count(None)+len(args)
+
+    p,N=size_(y,nargout=2)
+    L=size_(zeta,1)
+    theta=zeros_(p,L)
+    eta_tilde=zeros_(L,N)
+    for nn in arange_(1,N).reshape(-1):
+        eta_tilde[:,nn]=zeta[:,:,nn] * eta[:,nn]
+    eta_tilde=eta_tilde.T
+    for pp in arange_(1,p).reshape(-1):
+        inds_y_p=inds_y[pp,:].T
+        eta_tilde_p=eta_tilde.dot(inds_y_p[:,ones_(1,L)])
+        chol_Sig_theta_p_trans=numpy.linalg.solve(chol_(diag_(phi[pp,:].dot(tau)) + invSig_vec[pp] * (eta_tilde_p.T * eta_tilde_p)),eye_(L))
+        m_theta_p=invSig_vec[pp] * (chol_Sig_theta_p_trans * chol_Sig_theta_p_trans.T) * (eta_tilde_p.T * y[pp,:].T)
+        theta[pp,:]=m_theta_p + chol_Sig_theta_p_trans * randn_(L,1)
+    return theta
+def sample_sig_(y=None,theta=None,eta=None,zeta=None,prior_params=None,inds_y=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 6-[y,theta,eta,zeta,prior_params,inds_y].count(None)+len(args)
+
+    p,N=size_(y,nargout=2)
+    a_sig=prior_params.a_sig
+    b_sig=prior_params.b_sig
+    inds_vec=matlabarray([arange_(1,N)])
+    invSig_vec=zeros_(1,p)
+    for pp in arange_(1,p).reshape(-1):
+        sq_err=0
+        for nn in inds_vec[inds_y[pp,:]].reshape(-1):
+            sq_err=sq_err + (y[pp,nn] - theta[pp,:] * zeta[:,:,nn] * eta[:,nn]) ** 2
+        invSig_vec[pp]=gamrnd_(a_sig + 0.5 * sum_(inds_y[pp,:]),1) / (b_sig + 0.5 * sq_err)
+    return invSig_vec
+def sample_hypers_(theta=None,phi=None,tau=None,prior_params=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 4-[theta,phi,tau,prior_params].count(None)+len(args)
+
+    p,L=size_(theta,nargout=2)
+    a1=prior_params.a1
+    a2=prior_params.a2
+    a_phi=prior_params.a_phi
+    b_phi=prior_params.b_phi
+    a=matlabarray([a1,a2 * ones_(1,L - 1)])
+    delta=exp_([log_(tau[1]),diff_(log_(tau))])
+    for numIter in arange_(1,50).reshape(-1):
+        phi=gamrnd_((a_phi + 0.5) * ones_(p,L),1) / (b_phi + 0.5 * tau[ones_(1,p),:].dot((theta ** 2)))
+        sum_phi_theta=sum_(phi.dot((theta ** 2)),1)
+        for hh in arange_(1,L).reshape(-1):
+            tau_hh=exp_(cumsum_(log_(delta))).dot([zeros_(1,hh - 1),ones_(1,L - hh + 1) / delta[hh]])
+            delta[hh]=gamrnd_(a[hh] + 0.5 * p * (L - hh + 1),1) / (1 + 0.5 * sum_(tau_hh.dot(sum_phi_theta)))
+        tau=exp_(cumsum_(log_(delta)))
+    return phi,tau
+def sample_K_marg_zeta_(y=None,theta=None,eta=None,invSig_vec=None,prior_params=None,K_ind=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 6-[y,theta,eta,invSig_vec,prior_params,K_ind].count(None)+len(args)
+
+    k,N=size_(eta,nargout=2)
+    p,L=size_(theta,nargout=2)
+    c_prior=prior_params.c_prior
+    K=prior_params.K
+    grid_size=length_(c_prior)
+    Pk=- Inf * ones_(1,grid_size)
+    nbhd=copy_(Inf)
+    nbhd_vec=matlabarray([arange_(max_(K_ind - nbhd,1),min_(K_ind + nbhd,grid_size))])
+    tmp_mat_init=diag_(repmat_(1.0 / invSig_vec,1,N))
+    if p < sqrt_(N):
+        for cc in nbhd_vec.reshape(-1):
+            Kc=K[:,:,cc]
+            tmp_mat=copy_(tmp_mat_init)
+            for ll in arange_(1,L).reshape(-1):
+                theta_ll=theta[:,ll]
+                for kk in arange_(1,k).reshape(-1):
+                    eta_kk=eta[kk,:]
+                    tmp_mat=tmp_mat + kron_((eta_kk * eta_kk.T).dot(Kc),theta_ll * theta_ll.T)
+            Pk[cc]=normpdfln_(y[:],zeros_(p * N,1),tmp_mat)
+    else:
+        for cc in nbhd_vec.reshape(-1):
+            Kc=K[:,:,cc]
+            Kc=sparse_(Kc)
+            tmp_mat=copy_(tmp_mat_init)
+            for ll in arange_(1,L).reshape(-1):
+                theta_ll=theta[:,ll]
+                for kk in arange_(1,k).reshape(-1):
+                    eta_kk=eta[kk,:]
+                    eta_mat_kk=diag_(eta_kk)
+                    eta_theta_ll_kk=kron_(eta_mat_kk,theta_ll)
+                    tmp_mat=tmp_mat + eta_theta_ll_kk * Kc * eta_theta_ll_kk.T
+            Pk[cc]=normpdfln_(y[:],zeros_(p * N,1),tmp_mat)
+    nbhd_mask=zeros_(1,length_(prior_params.c_prior))
+    nbhd_mask[nbhd_vec]=1
+    Pk=Pk + log_(prior_params.c_prior).dot(nbhd_mask)
+    Pk=cumsum_(exp_(Pk - max_(Pk)))
+    K_ind=1 + sum_(Pk[end()] * rand_(1) > Pk)
+    return K_ind
+def sample_K_cond_zeta_(zeta=None,prior_params=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 2-[zeta,prior_params].count(None)+len(args)
+
+    L,k,N=size_(zeta,nargout=3)
+    c_prior=prior_params.c_prior
+    invK=prior_params.invK
+    zeta_tmp=reshape_(zeta,[L * k,N]).T
+    Pk=zeros_(1,length_(c_prior))
+    for ii in arange_(1,length_(c_prior)).reshape(-1):
+        Pk[ii]=- sum_(diag_(zeta_tmp.T * (invK[:,:,ii] * zeta_tmp)))
+    Pk=Pk - 0.5 * (L * k) * prior_params.logdetK + log_(prior_params.c_prior)
+    Pk=cumsum_(exp_(Pk - max_(Pk)))
+    K_ind=1 + sum_(Pk[end()] * rand_(1) > Pk)
+    return K_ind
+def normpdfln_(x=None,mu=None,Sigma=None,iSigma=None,*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = 4-[x,mu,Sigma,iSigma].count(None)+len(args)
+
+    if nargin < 4:
+        q=size_(Sigma,1)
+        iSigma=numpy.linalg.solve(Sigma,eye_(q))
+        logdet_iSigma=- 2 * sum_(log_(diag_(chol_(Sigma))))
+    else:
+        logdet_iSigma=2 * sum_(log_(diag_(chol_(iSigma))))
+    logp=0.5 * logdet_iSigma - 0.5 * (x - mu).T * (iSigma * (x - mu))
+    return logp