wolffd@0: function [data, clamped] = mk_mutilated_samples(bnet, ncases, max_clamp, usecell)
wolffd@0: % GEN_MUTILATED_SAMPLES Do random interventions and then draw random samples
wolffd@0: % [data, clamped] = gen_mutilated_samples(bnet, ncases, max_clamp, usecell)
wolffd@0: %
wolffd@0: % At each step, we pick a random subset of size 0 .. max_clamp, and 
wolffd@0: % clamp these nodes to random values.
wolffd@0: %
wolffd@0: % data(i,m) is the value of node i in case m.
wolffd@0: % clamped(i,m) = 1 if node i in case m was set by intervention.
wolffd@0: 
wolffd@0: if nargin < 4, usecell = 1; end
wolffd@0: 
wolffd@0: ns = bnet.node_sizes;
wolffd@0: n = length(bnet.dag);
wolffd@0: if usecell
wolffd@0:   data = cell(n, ncases);
wolffd@0: else
wolffd@0:   data = zeros(n, ncases);
wolffd@0: end
wolffd@0: clamped = zeros(n, ncases);
wolffd@0: 
wolffd@0: csubsets = subsets(1:n, max_clamp, 0); % includes the empty set
wolffd@0: distrib_cset = normalise(ones(1, length(csubsets)));
wolffd@0: 
wolffd@0: for m=1:ncases
wolffd@0:   cset = csubsets{sample_discrete(distrib_cset)};
wolffd@0:   nvals = prod(ns(cset));
wolffd@0:   distrib_cvals = normalise(ones(1, nvals));
wolffd@0:   cvals = ind2subv(ns(cset), sample_discrete(distrib_cvals));
wolffd@0:   mutilated_bnet = do_intervention(bnet, cset, cvals);
wolffd@0:   ev = sample_bnet(mutilated_bnet);
wolffd@0:   if usecell
wolffd@0:     data(:,m) = ev;
wolffd@0:   else
wolffd@0:     data(:,m) = cell2num(ev);
wolffd@0:   end
wolffd@0:   clamped(cset,m) = 1;
wolffd@0: end