function pot = mk_pot_from_logistic(bnode, bnet)
% MK_POT_FROM_LOGISTIC Make a potential from a logistic node.
% pot = mk_pot_from_logistic(bnode, bnet)

ns = bnet.node_sizes;
global xi


if ~subset(bnode.parents, bnet.observed_nodes)
  var = 1;
  xi = compute_var_param(bnode, bnet);
  %disp('hack mk_pot_from_logistic ********');
else
  var = 0;
end
  
switch bnet.jtree.pot_type
  case 'd',
    assert(subset(bnode.parents, bnet.observed_nodes));
    CPT = cond_dist_logistic(bnode, bnet);
    if ismember(bnode.self, bnet.observed_nodes)
      self_val = bnet.values(bnode.self);
      CPT = CPT(self_val);
    end
    ns(bnet.observed_nodes) = 1;
    ns(bnet.cts_nodes) = 1;
    fam = [bnode.parents bnode.self];
    pot = mk_table(fam, ns(fam));
    pot.T = myreshape(CPT, pot.sizes);
  case 'g',
    ns(bnet.observed_nodes) = 0;
    assert(ismember(bnode.self, bnet.observed_nodes));
    v = bnet.values(bnode.self);
    fam = [bnode.parents bnode.self];
    pot = mk_gpot(fam, ns(fam), 'c');
    if var
      %xi = compute_var_param(bnode, bnet, v-1);
      [g,h,K] = mk_var_approx(bnode.glim.w1, bnode.glim.b1, v-1, xi);
    else
      [g,h,K] = logistic_to_can_CPT(bnode, bnet, v);
    end
    pot.g = g; pot.h = h; pot.K = K;
  case 'cg',
    ns(intersect(bnet.observed_nodes, bnet.cts_nodes)) = 0;
    ns(intersect(bnet.observed_nodes, bnet.discrete_nodes)) = 1;
    pot = mk_cgpot(bnode.self, bnode.parents, ns, 'c');
    if ismember(bnode.self, bnet.observed_nodes)
      vs = bnet.values(bnode.self);
    else
      vs = 1:ns(bnode.self);
    end
    for j=1:length(vs)
      v = vs(j);
      if var
	%xi = compute_var_param(bnode, bnet, v-1);
	[pot.pot{j}.g, pot.pot{j}.h, pot.pot{j}.K] = mk_var_approx(bnode.glim.w1, bnode.glim.b1, v-1, xi);
      else
	[pot.pot{j}.g, pot.pot{j}.h, pot.pot{j}.K] = logistic_to_can_CPT(bnode, bnet, v);
      end
    end
  otherwise,
    error(['bad type ' bnet.jtree.pot_type]);
end


%%%%%
function [g,h,K] = logistic_to_can_CPT(bnode, bnet, self_val)

CPT = cond_dist_logistic(bnode, bnet);
prob = CPT(self_val);
g = log(prob);
n = bnode.parents_size;
h = zeros(n,1);
K = zeros(n,n);

%%%%%%%

function xi = compute_var_param(bnode, bnet)

% We use the variational lower bound in Jaakkola and Jordan, "A variational approach to Bayesian
% logistic regression models and their extensions", AI+Stats, 1996.
% See also Tipping, "Probabilistic visualisation of high-dimensional binary data", NIPS 98
% and Murphy, "A variational approximation for Bayesian networks with discrete and continuous
%   latent variables", UAI 99.

p = bnode.parents(1);
w = bnode.glim.w1;
b = bnode.glim.b1;

if ismember(p, bnet.observed_nodes)
  mu = bnet.values(p)
  xi = (w'*mu + b);
else
  mu = bnet.jtree.hybrid.marginal{p}.mu;
  S = bnet.jtree.hybrid.marginal{p}.Sigma;
  xi = sqrt(w'*(S+mu*mu')*w + 2*b*w'*mu + b^2);
end

% [g,h,K] are all symmetric in xi, so WLOG we can always use a positive coef for xi
% coef = 2r-1 = 1





