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