function [reconstructed_feas1, reconstructed_uct1] = bcmi_v1(feas1, mask1, gmm_mu1, gmm_sigma1, gmm_wei1, feamin1)

% BCMI_V1 Bounded conditional mean imputation with an approximate posterior
%
% Input:
%
% feas1      compressed spectral features (channels x frames)
% mask1      missing-data mask (channels x frames) reliable/unreliable 1/0
% gmm_mu1    GMM component means (dim x num components)
% gmm_sigma1 GMM component covariances (dim x dim x num components)
% gmm_wei1   GMM component weights (num components x 1)
% feamin1    minimum feature value
%
% Output:
%
% reconstructed_feas1  truncated posterior mean (channels x frames)
% reconstructed_uct1   truncated posterior variance (channels x frames)


% Copyright (C) 2014 Ulpu Remes
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy of
% this software and associated documentation files (the "Software"), to deal in
% the Software without restriction, including without limitation the rights to
% use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
% the Software, and to permit persons to whom the Software is furnished to do so,
% subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.


[feadim1, N1]=size(feas1);
[modeldim1,M1]=size(gmm_mu1);

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

% initialise variables

pos1=zeros(1,M1);
mixture_feas1=zeros(feadim1,M1);
mixture_uct1=zeros(feadim1,M1);
reconstructed_feas1=zeros(feadim1,N1);
reconstructed_uct1=zeros(feadim1,N1);

for n=1:N1
   
   unreliable_feas1 = find(mask1(:,n)==0);
   reliable_feas1 = find(mask1(:,n)==1);
  
   % reliable frames need not be reconstructed:
   
   if length(reliable_feas1)==feadim1
       reconstructed_feas1(:,n)=feas1(:,n);
       reconstructed_uct1(:,n)=zeros(feadim1,1);
       continue;
   end
   
   % reorganisation matrix for mixed frames
   
   eye1 = eye(feadim1);
   eye1 = eye1([reliable_feas1; unreliable_feas1],:);
   
    
    % FEATURE RECONSTRUCTION

    for m = 1:M1 

        if length(unreliable_feas1) < feadim1

        reorganised_feature=eye1*feas1(:,n);

        % compute mean and cov for unreliable data

        sigma1 = gmm_sigma1(:,:,m);
        sigma1 = eye1*sigma1*eye1';

        mu1 = gmm_mu1(:,m);
        mu1 = eye1*mu1;

        % define sr1, su1, sru1 etc.

        num_reliable1 = length(reliable_feas1);

        sr1 = sigma1(1:num_reliable1,1:num_reliable1);
        sm1 = sigma1(1:num_reliable1,num_reliable1+1:feadim1);
        sm2 = sigma1(num_reliable1+1:feadim1,1:num_reliable1);
        su1 = sigma1(num_reliable1+1:feadim1,num_reliable1+1:feadim1);

        mr2 = mu1(1:num_reliable1);
        mu2 = mu1(num_reliable1+1:feadim1);

        % marginal distribution for the unreliable components

        su3 = su1-sm2*pinv(sr1)*sm1;
        mu3 = mu2+sm2*pinv(sr1)*(reorganised_feature(1:num_reliable1)-mr2);

        %%%% feature estimate

        reorganised_feature(num_reliable1+1:feadim1)=mu3;
        mixture_feas1(:,m) = eye1'*reorganised_feature;
        mixture_uct1(:,m)=eye1'*[zeros(num_reliable1,1); max(diag(su3),0)];

        %%%% cluster posterior
        pos1(m)=gmm_wei1(m)*mvnpdf(reorganised_feature(1:num_reliable1),mr2,sr1)+exp(-700);

        else

        mixture_feas1(:,m) = gmm_mu1(:,m);
        mixture_uct1(:,m) = diag(gmm_sigma1(:,:,m));

        pos1(m)=gmm_wei1(m);
        end

    end

    reconstructed_feas1(:,n)=sum(bsxfun(@times,mixture_feas1,pos1/sum(pos1)),2);
    reconstructed_uct1(:,n)=sum(bsxfun(@times,mixture_uct1+power(bsxfun(@minus,mixture_feas1,reconstructed_feas1(:,n)),2),pos1/sum(pos1)),2);

end

% truncate posterior

lbm = (feamin1-reconstructed_feas1)./max(sqrt(reconstructed_uct1),exp(-700));
lbm = max(lbm,-exp(100));
ubm = (feas1-reconstructed_feas1)./max(sqrt(reconstructed_uct1),exp(-700));
denominator1 = max(normcdf(ubm)-normcdf(lbm),exp(-700));
normlbm = normpdf(lbm)./denominator1;
normubm = normpdf(ubm)./denominator1;



reconstructed_feas1 = reconstructed_feas1+(normlbm-normubm).*sqrt(reconstructed_uct1);
reconstructed_uct1 = reconstructed_uct1.*(1+lbm.*normlbm-ubm.*normubm-power(normlbm-normubm,2));
