function reconstructed_feas1 = mapi(feas1, mask1, gmm_model1, feamin1, Hm, fm)

% MAPI GMM-MAP-based missing feature reconstruction (Raj et al. 2004)
% 
% Input:
%
% feas1      compressed spectral features (channels x frames)
% mask1      missing-data mask (channels x frames) reliable/unreliable 1/0
% feamin1    minimum feature value (1x1)
% Hm         (optional)
% fm         (optional)
%
% Clean speech model:
%
% gmm_model1.mu (dim x num components)
% gmm_model1.sigma (dim x dim x num components)
% gmm_model1.weight (num components x 1)
%
% Output: 
% reconstructed_feas1 (channels x frames)

% (C) 2014 Ulpu Remes
% MIT license
% For license terms and references, see README.txt


% Check model and feature dimensions

[feadim1, N1]=size(feas1);
[modeldim1, M1]=size(gmm_model1.mu); 

assert(modeldim1==feadim1,'Feature dimension does not match model dimension');

if ~exist('Hm','var')||~exist('fm','var') % these variables are independent from the test data and can therefore be precomputed

    % calculate the quadratic and linear terms in 1/2*x'*H*x + f'*x

    Hm=zeros(modeldim1,modeldim1,M1);
    fm=zeros(modeldim1,M1);

    for m=1:M1

        Hm(:,:,m)=pinv(gmm_model1.sigma(:,:,m));
        fm(:,m)=-1*Hm(:,:,m)'*gmm_model1.mu(:,m);
    end 
end

% initialise variables

lik1=zeros(feadim1,1);
pos1=zeros(1,M1);
reconstructed1=zeros(feadim1,M1);
reconstructed_feas1=zeros(feadim1,N1);

for n=1:N1
    
   if sum(mask1(:,n))==feadim1 continue; end
   
   lbm=~mask1(:,n).*(feamin1*ones(feadim1,1))+mask1(:,n).*feas1(:,n);
   ubm=feas1(:,n);

   for m=1:M1
  
       % 2.1) calculate approx observation likelihoods and P(m|observation)
       
       for k = 1:feadim1
           if mask1(k,n)==0
           lik1(k)=0.5*erf((feas1(k,n)-gmm_model1.mu(k,m))/sqrt(2*gmm_model1.sigma(k,k,m)))+0.5;
           else
           lik1(k)=normpdf(feas1(k,n), gmm_model1.mu(k,m), sqrt(gmm_model1.sigma(k,k,m)));
           end
       end
       
       pos1(m) = gmm_model1.weight(m)*prod(lik1)+exp(-700);
       
       % 2.2) bounded optimisation
       
       reconstructed1(:,m)=qps_as(Hm(:,:,m),fm(:,m),lbm,ubm);
       
   end
   
   reconstructed_feas1(:,n)=sum(bsxfun(@times,reconstructed1,pos1/sum(pos1)),2);
end

% do not replace the reliable components
reconstructed_feas1 = mask1.*feas1+~mask1.*reconstructed_feas1;

