Mauch MIREX 2010

syntax: glob

_beattracker

_chordtools

_chroma

_chromadata

_dbn

_FullBNT

_logfiles

_misc

_parameters

_segmentation

_song

_writetools

output

output_*

/assocarray.m/1.1.1.1/Wed May 29 15:59:52 2002//

/subsref.m/1.1.1.1/Wed Aug  4 19:36:30 2004//

D

FullBNT/BNT/@assocarray

:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt

function A = assocarray(keys, vals)

% ASSOCARRAY Make an associative array

% function A = assocarray(keys, vals)

%

% keys{i} is the i'th string, vals{i} is the i'th value.

% After construction, A('foo') will return the value associated with foo.

A.keys = keys;

A.vals = vals;

A = class(A, 'assocarray');

function val = subsref(A, S)

% SUBSREF Subscript reference for an associative array

% A('foo') will return the value associated with foo.

% If there are multiple identicaly keys, the first match is returned.

% Currently the search is sequential.

i = 1;

while i <= length(A.keys)

  if strcmp(S.subs{1}, A.keys{i})

    val = A.vals{i};

    return;

  end

  i = i + 1;

end

error(['can''t find ' S.subs{1}])

/boolean_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002//

D

FullBNT/BNT/CPDs/@boolean_CPD

:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt

function CPD = boolean_CPD(bnet, self, ftype, fname, pfail)

% BOOLEAN_CPD Make a tabular CPD representing a (noisy) boolean function

%

% CPD = boolean_cpd(bnet, self, 'inline', f) uses the inline function f

% to specify the CPT.

% e.g., suppose X4 = X2 AND (NOT X3). Then we can write

%    bnet.CPD{4} = boolean_CPD(bnet, 4, 'inline', inline('(x(1) & ~x(2)'));  

% Note that x(1) refers pvals(1) = X2, and x(2) refers to pvals(2)=X3.

%

% CPD = boolean_cpd(bnet, self, 'named', f) assumes f is a function name.

% f can be built-in to matlab, or a file.

% e.g., If X4 = X2 AND X3, we can write

%    bnet.CPD{4} = boolean_CPD(bnet, 4, 'named', 'and');

% e.g., If X4 = X2 OR X3, we can write

%    bnet.CPD{4} = boolean_CPD(bnet, 4, 'named', 'any');

%

% CPD = boolean_cpd(bnet, self, 'rnd') makes a random non-redundant bool fn.

%

% CPD = boolean_CPD(bnet, self, 'inline'/'named', f, pfail)

% will put probability mass 1-pfail on f(parents), and put pfail on the other value.

% This is useful for simulating noisy boolean functions.

% If pfail is omitted, it is set to 0.

% (Note that adding noise to a random (non-redundant) boolean function just creates a different

% (potentially redundant) random boolean function.)

%

% Note: This cannot be used to simulate a noisy-OR gate.

% Example: suppose C has parents A and B, and the

% link of A->C fails with prob pA and the link B->C fails with pB.

% Then the noisy-OR gate defines the following distribution

%

%  A  B  P(C=0)

%  0  0  1.0

%  1  0  pA

%  0  1  pB

%  1  1  pA * PB

% 

% By contrast, boolean_CPD(bnet, C, 'any', p) would define

%

%  A  B  P(C=0) 

%  0  0  1-p    

%  1  0  p      

%  0  1  p

%  1  1  p

if nargin==0

  % This occurs if we are trying to load an object from a file.

  CPD = tabular_CPD(bnet, self);

  return;

elseif isa(bnet, 'boolean_CPD')

  % This might occur if we are copying an object.

  CPD = bnet;

  return;

end

if nargin < 5, pfail = 0; end

ps = parents(bnet.dag, self);

ns = bnet.node_sizes;

psizes = ns(ps);

self_size = ns(self);

psucc = 1-pfail;

k = length(ps);

switch ftype

 case 'inline', f = eval_bool_fn(fname, k);

 case 'named',  f = eval_bool_fn(fname, k);

 case 'rnd',    f = mk_rnd_bool_fn(k);

 otherwise,     error(['unknown function type ' ftype]);

end

CPT = zeros(prod(psizes), self_size);

ndx = find(f==0);

CPT(ndx, 1) = psucc;

CPT(ndx, 2) = pfail;

ndx = find(f==1);

CPT(ndx, 2) = psucc;

CPT(ndx, 1) = pfail;

if k > 0

  CPT = reshape(CPT, [psizes self_size]);  

end

clamp = 1;

CPD = tabular_CPD(bnet, self, CPT, [], clamp);

%%%%%%%%%%%%

function f = eval_bool_fn(fname, n)

% EVAL_BOOL_FN Evaluate a boolean function on all bit vectors of length n

% f = eval_bool_fn(fname, n)

%

% e.g. f = eval_bool_fn(inline('x(1) & x(3)'), 3)

% returns   0     0     0     0     0     1     0     1

ns = 2*ones(1, n);

f = zeros(1, 2^n);

bits = ind2subv(ns, 1:2^n);

for i=1:2^n

  f(i) = feval(fname, bits(i,:)-1);

end

%%%%%%%%%%%%%%%

function f = mk_rnd_bool_fn(n)

% MK_RND_BOOL_FN Make a random bit vector of length n that encodes a non-redundant boolean function

% f = mk_rnd_bool_fn(n)

red = 1;

while red

  f = sample_discrete([0.5 0.5], 2^n, 1)-1;

  red = redundant_bool_fn(f);

end

%%%%%%%%

function red = redundant_bool_fn(f)

% REDUNDANT_BOOL_FN Does a boolean function depend on all its input values?

% r = redundant_bool_fn(f)

%

% f is a vector of length 2^n, representing the output for each bit vector.

% An input is redundant if there is no assignment to the other bits

% which changes the output e.g., input 1 is redundant if u(2:n) s.t.,

% f([0 u(2:n)]) <> f([1 u(2:n)]). 

% A function is redundant it it has any redundant inputs.

n = log2(length(f));

ns = 2*ones(1,n);

red = 0;

for i=1:n

  ens = ns;

  ens(i) = 1;

  U = ind2subv(ens, 1:2^(n-1));

  U(:,i) = 1;

  f1 = f(subv2ind(ns, U));

  U(:,i) = 2;

  f2 = f(subv2ind(ns, U));

  if isequal(f1, f2)

    red = 1;

    return;

  end

end

%%%%%%%%%%

function [b, iter] = rnd_truth_table(N)

% RND_TRUTH_TABLE Construct the output of a random truth table s.t. each input is non-redundant

% b = rnd_truth_table(N)

%

% N is the number of inputs. 

% b is a random bit string of length N, representing the output of the truth table.

% Non-redundant means that, for each input position k,

% there are at least two bit patterns, u and v, that differ only in the k'th position,

% s.t., f(u) ~= f(v), where f is the function represented by b.

% We use rejection sampling to ensure non-redundancy.

%

% Example: b = [0 0 0 1  0 0 0 1] is indep of 3rd input (AND of inputs 1 and 2)

bits = ind2subv(2*ones(1,N), 1:2^N)-1;

redundant = 1;

iter = 0;

while redundant & (iter < 4)

  iter = iter + 1;

  b = sample_discrete([0.5 0.5], 1, 2^N)-1;

  redundant = 0;

  for i=1:N

    on = find(bits(:,i)==1);

    off = find(bits(:,i)==0);

    if isequal(b(on), b(off))

      redundant = 1;

      break;

    end

  end

end

/deterministic_CPD.m/1.1.1.1/Mon Oct  7 13:26:36 2002//

D

FullBNT/BNT/CPDs/@deterministic_CPD

:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt

function CPD = deterministic_CPD(bnet, self, fname, pfail)

% DETERMINISTIC_CPD Make a tabular CPD representing a (noisy) deterministic function

%

% CPD = deterministic_CPD(bnet, self, fname)

% This calls feval(fname, pvals) for each possible vector of parent values.

% e.g., suppose there are 2 ternary parents, then pvals = 

%  [1 1], [2 1], [3 1],   [1 2], [2 2], [3 2],   [1 3], [2 3], [3 3]

% If v = feval(fname, pvals(i)), then

%  CPD(x | parents=pvals(i)) = 1 if x==v, and = 0 if x<>v

% e.g., suppose X4 = X2 AND (NOT X3). Then

%    bnet.CPD{4} = deterministic_CPD(bnet, 4, inline('((x(1)-1) & ~(x(2)-1)) + 1'));  

% Note that x(1) refers pvals(1) = X2, and x(2) refers to pvals(2)=X3

% See also boolean_CPD.

%

% CPD = deterministic_CPD(bnet, self, fname, pfail)

% will put probability mass 1-pfail on f(parents), and distribute pfail over the other values.

% This is useful for simulating noisy deterministic functions.

% If pfail is omitted, it is set to 0.

%

if nargin==0

  % This occurs if we are trying to load an object from a file.

  CPD = tabular_CPD(bnet, self);

  return;

elseif isa(bnet, 'deterministic_CPD')

  % This might occur if we are copying an object.

  CPD = bnet;

  return;

end

if nargin < 4, pfail = 0; end

ps = parents(bnet.dag, self);

ns = bnet.node_sizes;

psizes = ns(ps);

self_size = ns(self);

psucc = 1-pfail;

CPT = zeros(prod(psizes), self_size);

pvals = zeros(1, length(ps));

for i=1:prod(psizes)

  pvals = ind2subv(psizes, i);

  x = feval(fname, pvals);

  %fprintf('%d ', [pvals x]); fprintf('\n');

  if psucc == 1

    CPT(i, x) = 1;

  else

    CPT(i, x) = psucc;

    rest = mysetdiff(1:self_size, x);

    CPT(i, rest) = pfail/length(rest);

  end

end

CPT = reshape(CPT, [psizes self_size]);  

CPD = tabular_CPD(bnet, self, 'CPT',CPT, 'clamped',1);

function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)

% CPD_TO_LAMBDA_MSG Compute lambda message (discrete)

% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)

% Pearl p183 eq 4.52

switch msg_type

  case 'd',

   T = prod_CPT_and_pi_msgs(CPD, n, ps, msg, p);

   mysize = length(msg{n}.lambda);

   lambda = dpot(n, mysize, msg{n}.lambda);

   T = multiply_by_pot(T, lambda);

   lam_msg = pot_to_marginal(marginalize_pot(T, p));

   lam_msg = lam_msg.T;           

 case 'g',

  error('discrete_CPD can''t create Gaussian msgs')

end

function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence)

% COMPUTE_PI Compute pi vector (discrete) 

% pi = compute_pi(CPD, msg_type, n, ps, msg, evidence)

% Pearl p183 eq 4.51

switch msg_type

  case 'd',

   T = prod_CPT_and_pi_msgs(CPD, n, ps, msg);

   pi = pot_to_marginal(marginalize_pot(T, n));

   pi = pi.T(:);                   

 case 'g', 

  error('can only convert discrete CPD to Gaussian pi if observed')

end

function pot = CPD_to_scgpot(CPD, domain, ns, cnodes, evidence)

% CPD_TO_SCGPOT Convert a CPD to a CG potential, incorporating any evidence (discrete)

% pot = CPD_to_scgpot(CPD, domain, ns, cnodes, evidence)

%

% domain is the domain of CPD.

% node_sizes(i) is the size of node i.

% cnodes

% evidence{i} is the evidence on the i'th node.

%odom = domain(~isemptycell(evidence(domain)));

%vals = cat(1, evidence{odom});

%map = find_equiv_posns(odom, domain);

%index = mk_multi_index(length(domain), map, vals);

CPT = CPD_to_CPT(CPD);

%CPT = CPT(index{:});

CPT = CPT(:);

%ns(odom) = 1;

potarray = cell(1, length(CPT));

for i=1:length(CPT)

  %p = CPT(i);

  potarray{i} = scgcpot(0, 0, CPT(i));

  %scpot{i} = scpot(0, 0);

end

pot = scgpot(domain, [], [], ns, potarray);

/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002//

/CPD_to_pi.m/1.1.1.1/Wed May 29 15:59:52 2002//

/CPD_to_scgpot.m/1.1.1.1/Wed May 29 15:59:52 2002//

/README/1.1.1.1/Wed May 29 15:59:52 2002//

/convert_CPD_to_table_hidden_ps.m/1.1.1.1/Wed May 29 15:59:52 2002//

/convert_obs_CPD_to_table.m/1.1.1.1/Wed May 29 15:59:52 2002//

/convert_to_pot.m/1.1.1.1/Fri Feb 20 22:00:38 2004//

/convert_to_sparse_table.c/1.1.1.1/Wed May 29 15:59:52 2002//

/convert_to_table.m/1.1.1.1/Wed May 29 15:59:52 2002//

/discrete_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002//

/dom_sizes.m/1.1.1.1/Wed May 29 15:59:52 2002//

/log_prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002//

/prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002//

/sample_node.m/1.1.1.1/Wed May 29 15:59:52 2002//

D

A D/Old////

A D/private////

FullBNT/BNT/CPDs/@discrete_CPD

:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt

/convert_to_pot.m/1.1.1.1/Wed May 29 15:59:52 2002//

/convert_to_table.m/1.1.1.1/Wed May 29 15:59:52 2002//

/prob_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002//

/prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002//

D

FullBNT/BNT/CPDs/@discrete_CPD/Old

:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt

function pot = convert_to_pot(CPD, pot_type, domain, evidence)

% CONVERT_TO_POT Convert a tabular CPD to one or more potentials

% pots = convert_to_pot(CPD, pot_type, domain, evidence)

%

% pots{i} = CPD evaluated using evidence(domain(:,i))

% If 'domains' is a single row vector, pots will be an object, not a cell array.

ncases = size(domain,2);

assert(ncases==1); % not yet vectorized

sz = dom_sizes(CPD);

ns = zeros(1, max(domain));

ns(domain) = sz;

local_ev = evidence(domain);

obs_bitv = ~isemptycell(local_ev);

odom = domain(obs_bitv);

T = convert_to_table(CPD, domain, local_ev, obs_bitv);

switch pot_type

 case 'u',

  pot = upot(domain, sz, T, 0*myones(sz));  

 case 'd',

  ns(odom) = 1;

  pot = dpot(domain, ns(domain), T);          

 case {'c','g'},

  % Since we want the output to be a Gaussian, the whole family must be observed.

  % In other words, the potential is really just a constant.

  p = T;

  %p = prob_node(CPD, evidence(domain(end)), evidence(domain(1:end-1)));

  ns(domain) = 0;

  pot = cpot(domain, ns(domain), log(p));       

 case 'cg',

  T = T(:);

  ns(odom) = 1;

  can = cell(1, length(T));

  for i=1:length(T)

    can{i} = cpot([], [], log(T(i)));

  end

  pot = cgpot(domain, [], ns, can);   

 otherwise,

  error(['unrecognized pot type ' pot_type])

end

function T = convert_to_table(CPD, domain, local_ev, obs_bitv)

% CONVERT_TO_TABLE Convert a discrete CPD to a table

% function T = convert_to_table(CPD, domain, local_ev, obs_bitv)

%

% We convert the CPD to a CPT, and then lookup the evidence on the discrete parents.

% The resulting table can easily be converted to a potential.

CPT = CPD_to_CPT(CPD);

obs_child_only = ~any(obs_bitv(1:end-1)) & obs_bitv(end);

if obs_child_only

  sz = size(CPT);

  CPT = reshape(CPT, prod(sz(1:end-1)), sz(end));

  o = local_ev{end};

  T = CPT(:, o);

else

  odom = domain(obs_bitv);  

  vals = cat(1, local_ev{find(obs_bitv)}); % undo cell array

  map = find_equiv_posns(odom, domain);

  index = mk_multi_index(length(domain), map, vals);

  T = CPT(index{:});

end

function p = prob_CPD(CPD, domain, ns, cnodes, evidence)

% PROB_CPD Compute prob of a node given evidence on the parents (discrete)

% p = prob_CPD(CPD, domain, ns, cnodes, evidence)

%

% domain is the domain of CPD.

% node_sizes(i) is the size of node i.

% cnodes = all the cts nodes

% evidence{i} is the evidence on the i'th node.

ps = domain(1:end-1);

self = domain(end);

CPT = CPD_to_CPT(CPD);

if isempty(ps)

  T = CPT;

else

  assert(~any(isemptycell(evidence(ps))));

  pvals = cat(1, evidence{ps});

  i = subv2ind(ns(ps), pvals(:)');

  T = reshape(CPT, [prod(ns(ps)) ns(self)]);

  T = T(i,:);

end

p = T(evidence{self});

function [P, p] = prob_node(CPD, self_ev, pev)

% PROB_NODE Compute prod_m P(x(i,m)| x(pi_i,m), theta_i) for node i (discrete)

% [P, p] = prob_node(CPD, self_ev, pev)

%

% self_ev(m) is the evidence on this node in case m.

% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).

% (These may also be cell arrays.)

%

% p(m) = P(x(i,m)| x(pi_i,m), theta_i) 

% P = prod p(m)

if iscell(self_ev), usecell = 1; else usecell = 0; end

ncases = length(self_ev);

sz = dom_sizes(CPD);

nparents = length(sz)-1;

if nparents == 0

  assert(isempty(pev));

else

  assert(isequal(size(pev), [nparents ncases]));

end

n = length(sz);

dom = 1:n;

p = zeros(1, ncases);

if nparents == 0

  for m=1:ncases

    if usecell

      evidence = {self_ev{m}};

    else

      evidence = num2cell(self_ev(m));

    end

    T = convert_to_table(CPD, dom, evidence);

    p(m) = T;

  end

else

  for m=1:ncases

    if usecell

      evidence = cell(1,n);

      evidence(1:n-1) = pev(:,m);

      evidence(n) = self_ev(m);

    else

      evidence = num2cell([pev(:,m)', self_ev(m)]);

    end

    T = convert_to_table(CPD, dom, evidence);

    p(m) = T;

  end

end

P = prod(p);

Any CPD on a discrete child with discrete parents

can be represented as a table (although this might be quite big).

discrete_CPD uses this tabular representation to implement various

functions. Subtypes are free to implement more efficient versions.

function T = convert_CPD_to_table_hidden_ps(CPD, child_obs)

% CONVERT_CPD_TO_TABLE_HIDDEN_PS Convert a discrete CPD to a table

% T = convert_CPD_to_table_hidden_ps(CPD, child_obs)

%

% This is like convert_to_table, except that we are guaranteed that

% none of the parents have evidence on them.

% child_obs may be an integer (1,2,...) or [].

CPT = CPD_to_CPT(CPD);

if isempty(child_obs)

  T = CPT(:);

else

  sz = dom_sizes(CPD);

  if length(sz)==1 % no parents

    T = CPT(child_obs);

  else

    CPT = reshape(CPT, prod(sz(1:end-1)), sz(end));

    T = CPT(:, child_obs);

  end

end

function T = convert_to_table(CPD, domain, evidence)

% CONVERT_TO_TABLE Convert a discrete CPD to a table

% T = convert_to_table(CPD, domain, evidence)

%

% We convert the CPD to a CPT, and then lookup the evidence on the discrete parents.

% The resulting table can easily be converted to a potential.

CPT = CPD_to_CPT(CPD);

odom = domain(~isemptycell(evidence(domain)));

vals = cat(1, evidence{odom});

map = find_equiv_posns(odom, domain);

index = mk_multi_index(length(domain), map, vals);

T = CPT(index{:});

function pot = convert_to_pot(CPD, pot_type, domain, evidence)

% CONVERT_TO_POT Convert a discrete CPD to a potential

% pot = convert_to_pot(CPD, pot_type, domain, evidence)

%

% pots = CPD evaluated using evidence(domain)

ncases = size(domain,2);

assert(ncases==1); % not yet vectorized

sz = dom_sizes(CPD);

ns = zeros(1, max(domain));

ns(domain) = sz;

CPT1 = CPD_to_CPT(CPD);

spar = issparse(CPT1);

odom = domain(~isemptycell(evidence(domain)));

if spar

   T = convert_to_sparse_table(CPD, domain, evidence);

else 

   T = convert_to_table(CPD, domain, evidence);

end

switch pot_type

 case 'u',

  pot = upot(domain, sz, T, 0*myones(sz));  

 case 'd',

  ns(odom) = 1;

  pot = dpot(domain, ns(domain), T);          

 case {'c','g'},

  % Since we want the output to be a Gaussian, the whole family must be observed.

  % In other words, the potential is really just a constant.

  p = T;

  %p = prob_node(CPD, evidence(domain(end)), evidence(domain(1:end-1)));

  ns(domain) = 0;

  pot = cpot(domain, ns(domain), log(p));       

 case 'cg',

  T = T(:);

  ns(odom) = 1;

  can = cell(1, length(T));

  for i=1:length(T)

    if T(i) == 0 

      can{i} = cpot([], [], -Inf); % bug fix by Bob Welch 20/2/04

    else

      can{i} = cpot([], [], log(T(i)));

    end;

  end

  pot = cgpot(domain, [], ns, can); 

 case 'scg'

  T = T(:);

  ns(odom) = 1;

  pot_array = cell(1, length(T));

  for i=1:length(T)

    pot_array{i} = scgcpot([], [], T(i));

  end

  pot = scgpot(domain, [], [], ns, pot_array);   

 otherwise,

  error(['unrecognized pot type ' pot_type])

end

/* convert_to_sparse_table.c  convert a sparse discrete CPD with evidence into sparse table */

/* convert_to_pot.m located in ../CPDs/discrete_CPD call it */

/* 3 input */

/* CPD      prhs[0] with 1D sparse CPT */

/* domain   prhs[1]                    */

/* evidence prhs[2]                    */

/* 1 output */

/* T        plhs[0] sparse table       */

#include <math.h>

#include "mex.h"

void ind_subv(int index, const int *cumprod, const int n, int *bsubv){

	int i;

	for (i = n-1; i >= 0; i--) {

		bsubv[i] = ((int)floor(index / cumprod[i]));

		index = index % cumprod[i];

	}

}

int subv_ind(const int n, const int *cumprod, const int *subv){

	int i, index=0;

	for(i=0; i<n; i++){

		index += subv[i] * cumprod[i];

	}

	return index;

}

void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){

	double *ptr;

	void   *newptr;

	int    *ir, *jc;

	int    nbytes;

	if(new_nzmax == old_nzmax) return;

	nbytes = new_nzmax * sizeof(*ptr);

	ptr = mxGetPr(spArray);

	newptr = mxRealloc(ptr, nbytes);

	mxSetPr(spArray, newptr);

	nbytes = new_nzmax * sizeof(*ir);

	ir = mxGetIr(spArray);

	newptr = mxRealloc(ir, nbytes);

	mxSetIr(spArray, newptr);

	jc = mxGetJc(spArray);

	jc[0] = 0;

	jc[1] = new_nzmax;

	mxSetNzmax(spArray, new_nzmax);

}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){

	int     i, j, NS, NZB, count, bdim, match, domain, bindex, sindex, nzCounts=0;

	int     *observed, *bsubv, *ssubv, *bir, *sir, *bjc, *sjc, *mask, *ssize, *bcumprod, *scumprod;

	double  *pDomain, *pSize, *bpr, *spr;

	mxArray *pTemp;

	pTemp = mxGetField(prhs[0], 0, "CPT");

	bpr = mxGetPr(pTemp);

	bir = mxGetIr(pTemp);

	bjc = mxGetJc(pTemp);

	NZB = bjc[1];

	pTemp = mxGetField(prhs[0], 0, "sizes");

	pSize = mxGetPr(pTemp);

	pDomain = mxGetPr(prhs[1]);

	bdim = mxGetNumberOfElements(prhs[1]);

	mask = malloc(bdim * sizeof(int));

	ssize = malloc(bdim * sizeof(int));

	observed = malloc(bdim * sizeof(int));

	for(i=0; i<bdim; i++){

		ssize[i] = (int)pSize[i];

	}

	count = 0;

	for(i=0; i<bdim; i++){

		domain = (int)pDomain[i] - 1;

		pTemp = mxGetCell(prhs[2], domain);

		if(pTemp){

			mask[count] = i;

			ssize[i] = 1;

			observed[count] = (int)mxGetScalar(pTemp) - 1;

			count++;

		}

	}

	if(count == 0){

		pTemp = mxGetField(prhs[0], 0, "CPT");

		plhs[0] = mxDuplicateArray(pTemp);

		free(mask);

		free(ssize);

		free(observed);

		return;

	}

	bsubv = malloc(bdim * sizeof(int));

	ssubv = malloc(count * sizeof(int));

	bcumprod = malloc(bdim * sizeof(int));

	scumprod = malloc(bdim * sizeof(int));

	NS = 1;

	for(i=0; i<bdim; i++){

		NS *= ssize[i];

	}

	plhs[0] = mxCreateSparse(NS, 1, NS, mxREAL);

	spr = mxGetPr(plhs[0]);

	sir = mxGetIr(plhs[0]);

	sjc = mxGetJc(plhs[0]);

	sjc[0] = 0;

	sjc[1] = NS;

	bcumprod[0] = 1;

	scumprod[0] = 1;

	for(i=0; i<bdim-1; i++){

		bcumprod[i+1] = bcumprod[i] * (int)pSize[i];

		scumprod[i+1] = scumprod[i] * ssize[i];

	}

	nzCounts = 0;

	for(i=0; i<NZB; i++){

		bindex = bir[i];

		ind_subv(bindex, bcumprod, bdim, bsubv);

		for(j=0; j<count; j++){

			ssubv[j] = bsubv[mask[j]];

		}

		match = 1;

		for(j=0; j<count; j++){

			if((ssubv[j]) != observed[j]){

				match = 0;

				break;

			}

		}

		if(match){

			spr[nzCounts] = bpr[i];

			sindex = subv_ind(bdim, scumprod, bsubv);

			sir[nzCounts] = sindex;

			nzCounts++;

		}

	}

	reset_nzmax(plhs[0], NS, nzCounts);

	free(mask);

	free(ssize);

	free(observed);

	free(bsubv);

	free(ssubv);

	free(bcumprod);

	free(scumprod);

}

function T = convert_to_table(CPD, domain, evidence)

% CONVERT_TO_TABLE Convert a discrete CPD to a table

% T = convert_to_table(CPD, domain, evidence)

%

% We convert the CPD to a CPT, and then lookup the evidence on the discrete parents.

% The resulting table can easily be converted to a potential.

domain = domain(:);

CPT = CPD_to_CPT(CPD);

odom = domain(~isemptycell(evidence(domain)));

vals = cat(1, evidence{odom});

map = find_equiv_posns(odom, domain);

index = mk_multi_index(length(domain), map, vals);

T = CPT(index{:});

T = T(:);

function CPD = discrete_CPD(clamped, dom_sizes)

% DISCRETE_CPD Virtual constructor for generic discrete CPD

% CPD = discrete_CPD(clamped, dom_sizes)

CPD.dom_sizes = dom_sizes;

CPD = class(CPD, 'discrete_CPD', generic_CPD(clamped));

function sz = dom_sizes(CPD)

% DOM_SIZES Return the size of each node in the domain

% sz = dom_sizes(CPD)

sz = CPD.dom_sizes;

function L = log_prob_node(CPD, self_ev, pev)

% LOG_PROB_NODE Compute sum_m log P(x(i,m)| x(pi_i,m), theta_i) for node i (discrete)

% L = log_prob_node(CPD, self_ev, pev)

%

% self_ev(m) is the evidence on this node in case m.

% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).

% (These may also be cell arrays.)

[P, p] = prob_node(CPD, self_ev, pev); % P may underflow, so we use p

tiny = exp(-700);

p = p + (p==0)*tiny; % replace 0s by tiny

L = sum(log(p));

/prod_CPT_and_pi_msgs.m/1.1.1.1/Wed May 29 15:59:52 2002//

D

FullBNT/BNT/CPDs/@discrete_CPD/private

:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt

function T = prod_CPT_and_pi_msgs(CPD, n, ps, msgs, except)

% PROD_CPT_AND_PI_MSGS Multiply the CPD and all the pi messages from parents, perhaps excepting one

% T = prod_CPY_and_pi_msgs(CPD, n, ps, msgs, except)

if nargin < 5, except = -1; end

dom = [ps n];

%ns = sparse(1, max(dom));

ns = zeros(1, max(dom));

CPT = CPD_to_CPT(CPD);

ns(dom) = mysize(CPT);

T = dpot(dom, ns(dom), CPT);

for i=1:length(ps)

  p = ps(i);

  if p ~= except

    T = multiply_by_pot(T, dpot(p, ns(p), msgs{n}.pi_from_parent{i}));

  end

end         

function [P, p] = prob_node(CPD, self_ev, pev)

% PROB_NODE Compute prod_m P(x(i,m)| x(pi_i,m), theta_i) for node i (discrete)

% [P, p] = prob_node(CPD, self_ev, pev)

%

% self_ev(m) is the evidence on this node in case m.

% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).

% (These may also be cell arrays.)

%

% p(m) = P(x(i,m)| x(pi_i,m), theta_i) 

% P = prod p(m)

if iscell(self_ev), usecell = 1; else usecell = 0; end

ncases = length(self_ev);

sz = dom_sizes(CPD);

nparents = length(sz)-1;

if nparents == 0

  assert(isempty(pev));

else

  assert(isequal(size(pev), [nparents ncases]));

end

n = length(sz);

dom = 1:n;

p = zeros(1, ncases);

if isa(CPD, 'tabular_CPD')

  % speed up by looking up CPT using index Zhang Yimin  2001-12-31

  if usecell

    if nparents == 0

      data = [cell2num(self_ev)]; 

    else

      data = [cell2num(pev); cell2num(self_ev)]; 

    end

  else

    if nparents == 0

      data = [self_ev];

    else

      data = [pev; self_ev];

    end

  end

  indices = subv2ind(sz, data'); % each row of data' is a case 

  CPT=CPD_to_CPT(CPD);

  p = CPT(indices);

  %get the prob list

  %cpt_size = prod(sz);

  %prob_list=reshape(CPT, cpt_size, 1);

  %for m=1:ncases  %here we assume we get evidence for node and all its parents

  %  idx=indices(m);

  %  p(m)=prob_list(idx); 

  %end

else % eg. softmax

  for m=1:ncases

    if usecell

      if nparents == 0

	evidence = {self_ev{m}};

      else

	evidence = cell(1,n);

	evidence(1:n-1) = pev(:,m);

	evidence(n) = self_ev(m);

      end

    else

      if nparents == 0

	evidence = num2cell(self_ev(m));

      else

	evidence = num2cell([pev(:,m)', self_ev(m)]);

      end

    end

    T = convert_to_table(CPD, dom, evidence);

    p(m) = T;

  end

end

P = prod(p);

function y = sample_node(CPD, pvals)

% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i)  (discrete)

% y = sample_node(CPD, parent_evidence)

%

% parent_evidence{i} is the value of the i'th parent

if 0

n = length(pvals)+1;

dom = 1:n;

evidence = cell(1,n);

evidence(1:n-1) = pvals;

T = convert_to_table(CPD, dom, evidence);

y = sample_discrete(T);

end

CPT = CPD_to_CPT(CPD);

sz = mysize(CPT);

nparents = length(sz)-1;

switch nparents

 case 0, T = CPT;

 case 1, T = CPT(pvals{1}, :);

 case 2, T = CPT(pvals{1}, pvals{2}, :);

 case 3, T = CPT(pvals{1}, pvals{2}, pvals{3}, :);

 case 4, T = CPT(pvals{1}, pvals{2}, pvals{3}, pvals{4}, :);

 otherwise,

  pvals = cat(1, pvals{:});

  psz = sz(1:end-1);

  ssz = sz(end);

  i = subv2ind(psz, pvals(:)');

  T = reshape(CPT, [prod(psz) ssz]);

  T = T(i,:);

end

y = sample_discrete(T);

function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)

% CPD_TO_LAMBDA_MSG Compute lambda message (gaussian)

% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)

% Pearl p183 eq 4.52

switch msg_type

 case 'd',

  error('gaussian_CPD can''t create discrete msgs')

 case 'g',

  cps = ps(CPD.cps);

  cpsizes = CPD.sizes(CPD.cps);

  self_size = CPD.sizes(end);

  i = find_equiv_posns(p, cps); % p is n's i'th cts parent

  psz = cpsizes(i);

  if all(msg{n}.lambda.precision == 0) % no info to send on

    lam_msg.precision = zeros(psz, psz);

    lam_msg.info_state = zeros(psz, 1);

    return;

  end

  [m, Q, W] = gaussian_CPD_params_given_dps(CPD, [ps n], evidence);

  Bmu = m;

  BSigma = Q;

  for k=1:length(cps) % only get pi msgs from cts parents

    pk = cps(k);

    if pk ~= p

      %bk = block(k, cpsizes);

      bk = CPD.cps_block_ndx{k};

      Bk = W(:, bk);

      m = msg{n}.pi_from_parent{k}; 

      BSigma = BSigma + Bk * m.Sigma * Bk';

      Bmu = Bmu + Bk * m.mu;

    end

  end

  % BSigma = Q + sum_{k \neq i} B_k Sigma_k B_k'

  %bi = block(i, cpsizes);

  bi = CPD.cps_block_ndx{i};

  Bi = W(:,bi);

  P = msg{n}.lambda.precision;

  if (rcond(P) > 1e-3) | isinf(P)

    if isinf(P) % Y is observed

      Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance

      mu_lambda = msg{n}.lambda.mu; % observed_value;

    else

      Sigma_lambda = inv(P);

      mu_lambda = Sigma_lambda * msg{n}.lambda.info_state;

    end

    C = inv(Sigma_lambda + BSigma);

    lam_msg.precision = Bi' * C * Bi;

    lam_msg.info_state = Bi' * C * (mu_lambda - Bmu);

  else

    % method that uses matrix inversion lemma to avoid inverting P

    A = inv(P + inv(BSigma));

    C = P - P*A*P;

    lam_msg.precision = Bi' * C * Bi;

    D = eye(self_size) - P*A;

    z = msg{n}.lambda.info_state;

    lam_msg.info_state = Bi' * (D*z - D*P*Bmu);

  end

end

function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence)

% CPD_TO_PI Compute the pi vector (gaussian)

% function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence)

switch msg_type

 case 'd',

  error('gaussian_CPD can''t create discrete msgs')

 case 'g',

  [m, Q, W] = gaussian_CPD_params_given_dps(CPD, [ps n], evidence);

  cps = ps(CPD.cps);

  cpsizes = CPD.sizes(CPD.cps);

  pi.mu = m;

  pi.Sigma = Q;

  for k=1:length(cps) % only get pi msgs from cts parents

    %bk = block(k, cpsizes);

    bk = CPD.cps_block_ndx{k};

    Bk = W(:, bk);

    m = msg{n}.pi_from_parent{k}; 

    pi.Sigma = pi.Sigma + Bk * m.Sigma * Bk';

    pi.mu = pi.mu + Bk * m.mu; % m.mu = u(k)

  end

end

function pot = CPD_to_scgpot(CPD, domain, ns, cnodes, evidence)

% CPD_TO_CGPOT Convert a Gaussian CPD to a CG potential, incorporating any evidence   

% pot = CPD_to_cgpot(CPD, domain, ns, cnodes, evidence)

self = CPD.self;

dnodes = mysetdiff(1:length(ns), cnodes);

odom = domain(~isemptycell(evidence(domain)));

cdom = myintersect(cnodes, domain);

cheaddom = myintersect(self, domain);

ctaildom = mysetdiff(cdom,cheaddom);

ddom = myintersect(dnodes, domain);

cobs = myintersect(cdom, odom);

dobs = myintersect(ddom, odom);

ens = ns; % effective node size

ens(cobs) = 0;

ens(dobs) = 1;

% Extract the params compatible with the observations (if any) on the discrete parents (if any)

% parents are all but the last domain element

ps = domain(1:end-1);

dps = myintersect(ps, ddom);

dops = myintersect(dps, odom);

map = find_equiv_posns(dops, dps);

dpvals = cat(1, evidence{dops});

index = mk_multi_index(length(dps), map, dpvals);

dpsize = prod(ens(dps));

cpsize = size(CPD.weights(:,:,1), 2); % cts parents size

ss = size(CPD.mean, 1); % self size

% the reshape acts like a squeeze

m = reshape(CPD.mean(:, index{:}), [ss dpsize]);

C = reshape(CPD.cov(:, :, index{:}), [ss ss dpsize]);

W = reshape(CPD.weights(:, :, index{:}), [ss cpsize dpsize]);

% Convert each conditional Gaussian to a canonical potential

pot = cell(1, dpsize);

for i=1:dpsize

  %pot{i} = linear_gaussian_to_scgcpot(m(:,i), C(:,:,i), W(:,:,i), cdom, ns, cnodes, evidence);

  pot{i} = scgcpot(ss, cpsize, 1, m(:,i), W(:,:,i), C(:,:,i));

end

pot = scgpot(ddom, cheaddom, ctaildom, ens, pot);

function pot = linear_gaussian_to_scgcpot(mu, Sigma, W, domain, ns, cnodes, evidence)

% LINEAR_GAUSSIAN_TO_CPOT Convert a linear Gaussian CPD  to a stable conditional potential element.

% pot = linear_gaussian_to_cpot(mu, Sigma, W, domain, ns, cnodes, evidence)

p = 1;

A = mu;

B = W;

C = Sigma;

ns(odom) = 0;

%pot = scgcpot(, ns(domain), p, A, B, C);

/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002//

/CPD_to_pi.m/1.1.1.1/Wed May 29 15:59:52 2002//

/CPD_to_scgpot.m/1.1.1.1/Wed May 29 15:59:52 2002//

/adjustable_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002//

/convert_CPD_to_table_hidden_ps.m/1.1.1.1/Wed May 29 15:59:52 2002//

/convert_to_pot.m/1.1.1.1/Sun Mar  9 23:03:16 2003//

/convert_to_table.m/1.1.1.1/Sun May 11 23:31:54 2003//

/display.m/1.1.1.1/Wed May 29 15:59:52 2002//

/gaussian_CPD.m/1.1.1.1/Wed Jun 15 21:13:06 2005//

/gaussian_CPD_params_given_dps.m/1.1.1.1/Sun May 11 23:13:40 2003//

/get_field.m/1.1.1.1/Wed May 29 15:59:52 2002//

/learn_params.m/1.1.1.1/Thu Jun 10 01:28:10 2004//

/log_prob_node.m/1.1.1.1/Tue Sep 10 17:44:00 2002//

/maximize_params.m/1.1.1.1/Tue May 20 14:10:06 2003//

/maximize_params_debug.m/1.1.1.1/Fri Jan 31 00:13:10 2003//

/reset_ess.m/1.1.1.1/Wed May 29 15:59:52 2002//

/sample_node.m/1.1.1.1/Wed May 29 15:59:52 2002//

/set_fields.m/1.1.1.1/Wed May 29 15:59:52 2002//

/update_ess.m/1.1.1.1/Tue Jul 22 22:55:46 2003//

D

A D/Old////

A D/private////

FullBNT/BNT/CPDs/@gaussian_CPD

:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt

function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p)

% CPD_TO_LAMBDA_MSG Compute lambda message (gaussian)

% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p)

% Pearl p183 eq 4.52

switch msg_type

 case 'd',

  error('gaussian_CPD can''t create discrete msgs')

 case 'g',

  self_size = CPD.sizes(end);

  if all(msg{n}.lambda.precision == 0) % no info to send on

    lam_msg.precision = zeros(self_size);

    lam_msg.info_state = zeros(self_size, 1);

    return;

  end

  cpsizes = CPD.sizes(CPD.cps);

  dpval = 1;

  Q = CPD.cov(:,:,dpval);

  Sigmai = Q;

  wmu = zeros(self_size, 1);

  for k=1:length(ps)

    pk = ps(k);

    if pk ~= p

      bk = block(k, cpsizes);

      Bk = CPD.weights(:, bk, dpval);

      m = msg{n}.pi_from_parent{k};

      Sigmai = Sigmai + Bk * m.Sigma * Bk';

      wmu = wmu + Bk * m.mu; % m.mu = u(k)

    end

  end

  % Sigmai = Q + sum_{k \neq i} B_k Sigma_k B_k'

  i = find_equiv_posns(p, ps);

  bi = block(i, cpsizes);

  Bi = CPD.weights(:,bi, dpval);

  if 0

  P = msg{n}.lambda.precision;

  if isinf(P) % inv(P)=Sigma_lambda=0

    precision_temp = inv(Sigmai);

    lam_msg.precision = Bi' * precision_temp * Bi;

    lam_msg.info_state = precision_temp * (msg{n}.lambda.mu - wmu);

  else

    A = inv(P + inv(Sigmai));

    precision_temp = P + P*A*P;

    lam_msg.precision = Bi' * precision_temp * Bi;

    self_size = length(P);

    C = eye(self_size) + P*A;

    z = msg{n}.lambda.info_state;

    lam_msg.info_state = C*z - C*P*wmu;

  end

  end

  if isinf(msg{n}.lambda.precision)

    Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance

    mu_lambda = msg{n}.lambda.mu; % observed_value;

  else

    Sigma_lambda = inv(msg{n}.lambda.precision);

    mu_lambda = Sigma_lambda * msg{n}.lambda.info_state;

  end

  precision_temp = inv(Sigma_lambda + Sigmai);

  lam_msg.precision = Bi' * precision_temp * Bi;

  lam_msg.info_state = Bi' * precision_temp * (mu_lambda - wmu);

end

/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002//

/gaussian_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002//

/log_prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002//

/maximize_params.m/1.1.1.1/Thu Jan 30 22:38:16 2003//

/update_ess.m/1.1.1.1/Wed May 29 15:59:52 2002//

/update_tied_ess.m/1.1.1.1/Wed May 29 15:59:52 2002//

D

FullBNT/BNT/CPDs/@gaussian_CPD/Old

:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt

function CPD = gaussian_CPD(varargin)

% GAUSSIAN_CPD Make a conditional linear Gaussian distrib.

%

% To define this CPD precisely, call the continuous (cts) parents (if any) X,

% the discrete parents (if any) Q, and this node Y. Then the distribution on Y is:

% - no parents: Y ~ N(mu, Sigma)

% - cts parents : Y|X=x ~ N(mu + W x, Sigma)

% - discrete parents: Y|Q=i ~ N(mu(i), Sigma(i))

% - cts and discrete parents: Y|X=x,Q=i ~ N(mu(i) + W(i) x, Sigma(i))

%

% CPD = gaussian_CPD(bnet, node, ...) will create a CPD with random parameters,

% where node is the number of a node in this equivalence class.

%

% The list below gives optional arguments [default value in brackets].

% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y) and Q = prod(ns(Q)).)

%

% mean       - mu(:,i) is the mean given Q=i [ randn(Y,Q) ]

% cov        - Sigma(:,:,i) is the covariance given Q=i [ repmat(eye(Y,Y), [1 1 Q]) ]

% weights    - W(:,:,i) is the regression matrix given Q=i [ randn(Y,X,Q) ]

% cov_type   - if 'diag', Sigma(:,:,i) is diagonal [ 'full' ]

% tied_cov   - if 1, we constrain Sigma(:,:,i) to be the same for all i [0]

% clamp_mean - if 1, we do not adjust mu(:,i) during learning [0]

% clamp_cov  - if 1, we do not adjust Sigma(:,:,i) during learning [0]

% clamp_weights - if 1, we do not adjust W(:,:,i) during learning [0]

% cov_prior_weight - weight given to I prior for estimating Sigma [0.01]

%

% e.g., CPD = gaussian_CPD(bnet, i, 'mean', [0; 0], 'clamp_mean', 'yes')

%

% For backwards compatibility with BNT2, you can also specify the parameters in the following order

%   CPD = gaussian_CPD(bnet, self, mu, Sigma, W, cov_type, tied_cov, clamp_mean, clamp_cov, clamp_weight)

%

% Sometimes it is useful to create an "isolated" CPD, without needing to pass in a bnet.

% In this case, you must specify the discrete and cts parents (dps, cps) and the family sizes, followed

% by the optional arguments above:

%   CPD = gaussian_CPD('self', i, 'dps', dps, 'cps', cps, 'sz', fam_size, ...)

if nargin==0

  % This occurs if we are trying to load an object from a file.

  CPD = init_fields;

  clamp = 0;

  CPD = class(CPD, 'gaussian_CPD', generic_CPD(clamp));

  return;

elseif isa(varargin{1}, 'gaussian_CPD')

  % This might occur if we are copying an object.

  CPD = varargin{1};

  return;

end

CPD = init_fields;

CPD = class(CPD, 'gaussian_CPD', generic_CPD(0));

% parse mandatory arguments

if ~isstr(varargin{1}) % pass in bnet

  bnet = varargin{1};

  self = varargin{2};

  args = varargin(3:end);

  ns = bnet.node_sizes;

  ps = parents(bnet.dag, self);

  dps = myintersect(ps, bnet.dnodes);

  cps = myintersect(ps, bnet.cnodes);

  fam_sz = ns([ps self]);

else

  disp('parsing new style')

  for i=1:2:length(varargin)

    switch varargin{i},

     case 'self', self = varargin{i+1}; 

     case 'dps',  dps = varargin{i+1};

     case 'cps',  cps = varargin{i+1};

     case 'sz',   fam_sz = varargin{i+1};

    end

  end

  ps = myunion(dps, cps);

  args = varargin;

end

CPD.self = self;

CPD.sizes = fam_sz;

% Figure out which (if any) of the parents are discrete, and which cts, and how big they are