R/gllim_inverse_map.R

In epertham/xLLiM: High Dimensional Locally-Linear Mapping

gllim_inverse_map = function(y,theta,verb=0){
# %%%%%%%%%%%%%%%%% Inverse Mapping from Gllim Parameters %%%%%%%%%%%%%%%%%%%
# %%%% Author: Antoine Deleforge (July 2012) - antoine.deleforge@inria.fr %%%
# %% Converted to R: Emeline Perthame (2016) - perthame.emeline@gmail.com %%%
# % Description: Map N observations y using the inverse conditional
# % expectation E[x|y;theta] of the gllim model with parameters theta.
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %%%% Input %%%%
# % - y (DxN)               % Input observations to map
# % - theta  (struct)       % Gllim model parameters
# %   - theta.c (LxK)       % Gaussian means of X's prior
# %   - theta.Gamma (LxLxK) % Gaussian covariances of X's prior
# %   - theta.pi (1xK)      % Gaussian weights of X's prior
# %   - theta.A (DxLxK)     % Affine transformation matrices
# %   - theta.b (DxK)       % Affine transformation vectors
# %   - theta.Sigma (DxDxK) % Error covariances
# % - verb {0,1,2}          % Verbosity (def 1)
# %%%% Output %%%%
# % - x_exp (LxN)           % Posterior mean estimates E[xn|yn;theta]
# % - alpha (NxK)           % Weights of the posterior GMMs
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	D = nrow(y) ; N = ncol(y) ;
	L = nrow(theta$c) ; K = ncol(theta$c) ;
	
# % ======================Inverse density parameters=========================
	if(verb>=1) print('Compute K projections to X space and weights')
	proj=array(NaN,dim=c(L,N,K));     
	logalpha=matrix(0,N,K); 
	for (k in 1:K){
    if(verb>=2) print(paste("k=",k,sep=""));
    if(verb>=2) print('AbcG ');
    if (L==1) {Ak=theta$A[,,k,drop=FALSE];} else {Ak=matrix(theta$A[,,k],ncol=L,nrow=D);} # % DxL
    bk=theta$b[,k]; # % Dx1
    Sigmak=theta$Sigma[,,k]; # %DxD  ## OK 
    if (L==1) {ck=theta$c[,k,drop=FALSE];} else {ck=theta$c[,k];} # % Lx1
    Gammak=theta$Gamma[,,k]; # % LxL ## OK 
    
    if(verb>=2) print('cks ');
    cks=Ak%*%ck+bk; ## OK 
    
    if(verb>=2) print('Gks ');
    Gammaks=Sigmak+Ak%*%tcrossprod(Gammak,Ak); ## OK 
	
    if(verb>=2) print('iSks ');
    iSk = solve(Sigmak) #diag(1/diag(Sigmak))
    invSigmaks2=diag(L)+Gammak%*%crossprod(Ak, iSk)%*%Ak; 

    if(verb>=2) print('Aks ');
    Aks=solve(invSigmaks2,Gammak)%*%crossprod(Ak, iSk)      

    if(verb>=2) print('bks ');
    bks= solve(invSigmaks2) %*% (ck-Gammak%*%crossprod(Ak, iSk%*%bk))
   
    if(verb>=2) print('projections '); 
    proj[,,k]=sweep(Aks%*%y,1,bks,"+"); 

    if(verb>=2) print('logalpha '); 
    logalpha[,k]=log(theta$pi[k])+loggausspdf(y,cks,Gammaks);
    if(verb>=2) print(''); 
	}
	den=logsumexp(logalpha,2);
	logalpha= sweep(logalpha,1,den,"-") 
	alpha=exp(logalpha); 
	
	x_exp = lapply(1:L,function(l) proj[l,,]*alpha)
	x_exp = lapply(x_exp,rowSums)
	x_exp = do.call(rbind, x_exp)

return(list(x_exp=x_exp,alpha=alpha))
}