wolffd@0: function CPD = boolean_CPD(bnet, self, ftype, fname, pfail)
wolffd@0: % BOOLEAN_CPD Make a tabular CPD representing a (noisy) boolean function
wolffd@0: %
wolffd@0: % CPD = boolean_cpd(bnet, self, 'inline', f) uses the inline function f
wolffd@0: % to specify the CPT.
wolffd@0: % e.g., suppose X4 = X2 AND (NOT X3). Then we can write
wolffd@0: %    bnet.CPD{4} = boolean_CPD(bnet, 4, 'inline', inline('(x(1) & ~x(2)'));  
wolffd@0: % Note that x(1) refers pvals(1) = X2, and x(2) refers to pvals(2)=X3.
wolffd@0: %
wolffd@0: % CPD = boolean_cpd(bnet, self, 'named', f) assumes f is a function name.
wolffd@0: % f can be built-in to matlab, or a file.
wolffd@0: % e.g., If X4 = X2 AND X3, we can write
wolffd@0: %    bnet.CPD{4} = boolean_CPD(bnet, 4, 'named', 'and');
wolffd@0: % e.g., If X4 = X2 OR X3, we can write
wolffd@0: %    bnet.CPD{4} = boolean_CPD(bnet, 4, 'named', 'any');
wolffd@0: %
wolffd@0: % CPD = boolean_cpd(bnet, self, 'rnd') makes a random non-redundant bool fn.
wolffd@0: %
wolffd@0: % CPD = boolean_CPD(bnet, self, 'inline'/'named', f, pfail)
wolffd@0: % will put probability mass 1-pfail on f(parents), and put pfail on the other value.
wolffd@0: % This is useful for simulating noisy boolean functions.
wolffd@0: % If pfail is omitted, it is set to 0.
wolffd@0: % (Note that adding noise to a random (non-redundant) boolean function just creates a different
wolffd@0: % (potentially redundant) random boolean function.)
wolffd@0: %
wolffd@0: % Note: This cannot be used to simulate a noisy-OR gate.
wolffd@0: % Example: suppose C has parents A and B, and the
wolffd@0: % link of A->C fails with prob pA and the link B->C fails with pB.
wolffd@0: % Then the noisy-OR gate defines the following distribution
wolffd@0: %
wolffd@0: %  A  B  P(C=0)
wolffd@0: %  0  0  1.0
wolffd@0: %  1  0  pA
wolffd@0: %  0  1  pB
wolffd@0: %  1  1  pA * PB
wolffd@0: % 
wolffd@0: % By contrast, boolean_CPD(bnet, C, 'any', p) would define
wolffd@0: %
wolffd@0: %  A  B  P(C=0) 
wolffd@0: %  0  0  1-p    
wolffd@0: %  1  0  p      
wolffd@0: %  0  1  p
wolffd@0: %  1  1  p
wolffd@0: 
wolffd@0: 
wolffd@0: if nargin==0
wolffd@0:   % This occurs if we are trying to load an object from a file.
wolffd@0:   CPD = tabular_CPD(bnet, self);
wolffd@0:   return;
wolffd@0: elseif isa(bnet, 'boolean_CPD')
wolffd@0:   % This might occur if we are copying an object.
wolffd@0:   CPD = bnet;
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin < 5, pfail = 0; end
wolffd@0: 
wolffd@0: ps = parents(bnet.dag, self);
wolffd@0: ns = bnet.node_sizes;
wolffd@0: psizes = ns(ps);
wolffd@0: self_size = ns(self);
wolffd@0: 
wolffd@0: psucc = 1-pfail;
wolffd@0: 
wolffd@0: k = length(ps);
wolffd@0: switch ftype
wolffd@0:  case 'inline', f = eval_bool_fn(fname, k);
wolffd@0:  case 'named',  f = eval_bool_fn(fname, k);
wolffd@0:  case 'rnd',    f = mk_rnd_bool_fn(k);
wolffd@0:  otherwise,     error(['unknown function type ' ftype]);
wolffd@0: end
wolffd@0: 
wolffd@0: CPT = zeros(prod(psizes), self_size);
wolffd@0: ndx = find(f==0);
wolffd@0: CPT(ndx, 1) = psucc;
wolffd@0: CPT(ndx, 2) = pfail;
wolffd@0: ndx = find(f==1);
wolffd@0: CPT(ndx, 2) = psucc;
wolffd@0: CPT(ndx, 1) = pfail;
wolffd@0: if k > 0
wolffd@0:   CPT = reshape(CPT, [psizes self_size]);  
wolffd@0: end
wolffd@0: 
wolffd@0: clamp = 1;
wolffd@0: CPD = tabular_CPD(bnet, self, CPT, [], clamp);
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%%%%
wolffd@0: 
wolffd@0: function f = eval_bool_fn(fname, n)
wolffd@0: % EVAL_BOOL_FN Evaluate a boolean function on all bit vectors of length n
wolffd@0: % f = eval_bool_fn(fname, n)
wolffd@0: %
wolffd@0: % e.g. f = eval_bool_fn(inline('x(1) & x(3)'), 3)
wolffd@0: % returns   0     0     0     0     0     1     0     1
wolffd@0: 
wolffd@0: ns = 2*ones(1, n);
wolffd@0: f = zeros(1, 2^n);
wolffd@0: bits = ind2subv(ns, 1:2^n);
wolffd@0: for i=1:2^n
wolffd@0:   f(i) = feval(fname, bits(i,:)-1);
wolffd@0: end
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%
wolffd@0: 
wolffd@0: function f = mk_rnd_bool_fn(n)
wolffd@0: % MK_RND_BOOL_FN Make a random bit vector of length n that encodes a non-redundant boolean function
wolffd@0: % f = mk_rnd_bool_fn(n)
wolffd@0: 
wolffd@0: red = 1;
wolffd@0: while red
wolffd@0:   f = sample_discrete([0.5 0.5], 2^n, 1)-1;
wolffd@0:   red = redundant_bool_fn(f);
wolffd@0: end
wolffd@0: 
wolffd@0: %%%%%%%%
wolffd@0: 
wolffd@0: 
wolffd@0: function red = redundant_bool_fn(f)
wolffd@0: % REDUNDANT_BOOL_FN Does a boolean function depend on all its input values?
wolffd@0: % r = redundant_bool_fn(f)
wolffd@0: %
wolffd@0: % f is a vector of length 2^n, representing the output for each bit vector.
wolffd@0: % An input is redundant if there is no assignment to the other bits
wolffd@0: % which changes the output e.g., input 1 is redundant if u(2:n) s.t.,
wolffd@0: % f([0 u(2:n)]) <> f([1 u(2:n)]). 
wolffd@0: % A function is redundant it it has any redundant inputs.
wolffd@0: 
wolffd@0: n = log2(length(f));
wolffd@0: ns = 2*ones(1,n);
wolffd@0: red = 0;
wolffd@0: for i=1:n
wolffd@0:   ens = ns;
wolffd@0:   ens(i) = 1;
wolffd@0:   U = ind2subv(ens, 1:2^(n-1));
wolffd@0:   U(:,i) = 1;
wolffd@0:   f1 = f(subv2ind(ns, U));
wolffd@0:   U(:,i) = 2;
wolffd@0:   f2 = f(subv2ind(ns, U));
wolffd@0:   if isequal(f1, f2)
wolffd@0:     red = 1;
wolffd@0:     return;
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%%
wolffd@0: 
wolffd@0: function [b, iter] = rnd_truth_table(N)
wolffd@0: % RND_TRUTH_TABLE Construct the output of a random truth table s.t. each input is non-redundant
wolffd@0: % b = rnd_truth_table(N)
wolffd@0: %
wolffd@0: % N is the number of inputs. 
wolffd@0: % b is a random bit string of length N, representing the output of the truth table.
wolffd@0: % Non-redundant means that, for each input position k,
wolffd@0: % there are at least two bit patterns, u and v, that differ only in the k'th position,
wolffd@0: % s.t., f(u) ~= f(v), where f is the function represented by b.
wolffd@0: % We use rejection sampling to ensure non-redundancy.
wolffd@0: %
wolffd@0: % Example: b = [0 0 0 1  0 0 0 1] is indep of 3rd input (AND of inputs 1 and 2)
wolffd@0: 
wolffd@0: bits = ind2subv(2*ones(1,N), 1:2^N)-1;
wolffd@0: redundant = 1;
wolffd@0: iter = 0;
wolffd@0: while redundant & (iter < 4)
wolffd@0:   iter = iter + 1;
wolffd@0:   b = sample_discrete([0.5 0.5], 1, 2^N)-1;
wolffd@0:   redundant = 0;
wolffd@0:   for i=1:N
wolffd@0:     on = find(bits(:,i)==1);
wolffd@0:     off = find(bits(:,i)==0);
wolffd@0:     if isequal(b(on), b(off))
wolffd@0:       redundant = 1;
wolffd@0:       break;
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: