function [opt_val, mapper] = gaussian_optimal(mu_src,mu_tgt,sigma_src,sigma_tgt)

	% compute objective (Computational OT book)
	sigma_src_half = sqrtm(sigma_src);
	temp = sigma_src_half*sigma_tgt*sigma_src_half;
    temp = (temp + temp')/2;
	temp1 = sqrtm(temp);
	opt_val = norm(mu_src-mu_tgt)^2 + trace(sigma_tgt+sigma_src) -2*trace(temp1);
	%fprintf('Gaussian optimal transport optimal objective: %g\n',opt_val);

	% % compute mapper (Computational OT book)
	% m = size(Xstr,1);
	%temp2 = inv(sqrtm(sigma_src));
	temp2 = inv(sigma_src_half);
	mapper = temp2*temp1*temp2;
	% Ypred = mu_tgt + temp3*( Xstr' - mu_src*ones(1,m) );
	% Ypred = Ypred'; % same orientation as Xstr

end