wolffd@0: function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
wolffd@0: % CPD_TO_LAMBDA_MSG Compute lambda message (gmux)
wolffd@0: % lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
wolffd@0: % Pearl p183 eq 4.52
wolffd@0: 
wolffd@0: % Let Y be this node, X1..Xn be the cts parents and M the discrete switch node.
wolffd@0: % e.g., for n=3, M=1
wolffd@0: %
wolffd@0: %  X1 X2 X3 M
wolffd@0: %   \
wolffd@0: %    \
wolffd@0: %      Y
wolffd@0: %
wolffd@0: % So the only case in which we send an informative message is if p=1=M.
wolffd@0: % To the other cts parents, we send the "know nothing" message.
wolffd@0: 
wolffd@0: switch msg_type
wolffd@0:  case 'd',
wolffd@0:   error('gaussian_CPD can''t create discrete msgs')
wolffd@0:  case 'g',
wolffd@0:   cps = ps(CPD.cps);
wolffd@0:   cpsizes = CPD.sizes(CPD.cps);
wolffd@0:   self_size = CPD.sizes(end);
wolffd@0:   i = find_equiv_posns(p, cps); % p is n's i'th cts parent
wolffd@0:   psz = cpsizes(i);
wolffd@0:   dps = ps(CPD.dps);
wolffd@0:   M = evidence{dps};
wolffd@0:   if isempty(M)
wolffd@0:     error('gmux node must have observed discrete parent')
wolffd@0:   end
wolffd@0:   P = msg{n}.lambda.precision;
wolffd@0:   if all(P == 0) | (cps(M) ~= p) % if we know nothing, or are sending to a disconnected parent
wolffd@0:     lam_msg.precision = zeros(psz, psz);
wolffd@0:     lam_msg.info_state = zeros(psz, 1);
wolffd@0:     return;
wolffd@0:   end
wolffd@0:   % We are sending a message to the only effectively connected parent.
wolffd@0:   % There are no other incoming pi messages.
wolffd@0:   Bmu = CPD.mean(:,M);
wolffd@0:   BSigma = CPD.cov(:,:,M);
wolffd@0:   Bi = CPD.weights(:,:,M);
wolffd@0:   if (det(P) > 0) | isinf(P) 
wolffd@0:     if isinf(P) % Y is observed
wolffd@0:       Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance
wolffd@0:       mu_lambda = msg{n}.lambda.mu; % observed_value;
wolffd@0:     else
wolffd@0:       Sigma_lambda = inv(P);
wolffd@0:       mu_lambda = Sigma_lambda * msg{n}.lambda.info_state;
wolffd@0:     end
wolffd@0:     C = inv(Sigma_lambda + BSigma);
wolffd@0:     lam_msg.precision = Bi' * C * Bi;
wolffd@0:     lam_msg.info_state = Bi' * C * (mu_lambda - Bmu);
wolffd@0:   else
wolffd@0:     % method that uses matrix inversion lemma
wolffd@0:     A = inv(P + inv(BSigma));
wolffd@0:     C = P - P*A*P;
wolffd@0:     lam_msg.precision = Bi' * C * Bi;
wolffd@0:     D = eye(self_size) - P*A;
wolffd@0:     z = msg{n}.lambda.info_state;
wolffd@0:     lam_msg.info_state = Bi' * (D*z - D*P*Bmu);
wolffd@0:   end
wolffd@0: end