Daniel@0: function score = score_dags(data, ns, dags, varargin)
Daniel@0: % SCORE_DAGS Compute the score of one or more DAGs
Daniel@0: % score = score_dags(data, ns, dags, varargin)
Daniel@0: %
Daniel@0: % data{i,m} = value of node i in case m (can be a cell array).
Daniel@0: % node_sizes(i) is the number of size of node i.
Daniel@0: % dags{g} is the g'th dag
Daniel@0: % score(g) is the score of the i'th dag
Daniel@0: %
Daniel@0: % The following optional arguments can be specified in the form of name/value pairs:
Daniel@0: % [default value in brackets]
Daniel@0: %
Daniel@0: % scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ]
Daniel@0: %              Currently, only networks with all tabular nodes support Bayesian scoring.
Daniel@0: % type       - type{i} is the type of CPD to use for node i, where the type is a string
Daniel@0: %              of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ]
Daniel@0: % params     - params{i} contains optional arguments passed to the CPD constructor for node i,
Daniel@0: %              or [] if none.  [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ]
Daniel@0: % discrete   - the list of discrete nodes [ 1:N ]
Daniel@0: % clamped    - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ]
Daniel@0: %
Daniel@0: % e.g., score = score_dags(data, ns, mk_all_dags(n), 'scoring_fn', 'bic', 'params', []);
Daniel@0: %
Daniel@0: % If the DAGs have a lot of families in common, we can cache the sufficient statistics,
Daniel@0: % making this potentially more efficient than scoring the DAGs one at a time.
Daniel@0: % (Caching is not currently implemented, however.)
Daniel@0: 
Daniel@0: [n ncases] = size(data);
Daniel@0: 
Daniel@0: % set default params
Daniel@0: type = cell(1,n);
Daniel@0: params = cell(1,n);
Daniel@0: for i=1:n
Daniel@0:   type{i} = 'tabular';
Daniel@0:   params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 };
Daniel@0: end
Daniel@0: scoring_fn = 'bayesian';
Daniel@0: discrete = 1:n;
Daniel@0: 
Daniel@0: u = [1:ncases]'; % DWH
Daniel@0: isclamped = 0; %DWH
Daniel@0: clamped = zeros(n, ncases);
Daniel@0: 
Daniel@0: args = varargin;
Daniel@0: nargs = length(args);
Daniel@0: for i=1:2:nargs
Daniel@0:   switch args{i},
Daniel@0:    case 'scoring_fn', scoring_fn = args{i+1};
Daniel@0:    case 'type',       type = args{i+1}; 
Daniel@0:    case 'discrete',   discrete = args{i+1}; 
Daniel@0:    case 'clamped',    clamped = args{i+1}, isclamped = 1; %DWH
Daniel@0:    case 'params',     if isempty(args{i+1}), params = cell(1,n); else params = args{i+1};  end
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: NG = length(dags);
Daniel@0: score = zeros(1, NG);
Daniel@0: for g=1:NG
Daniel@0:   dag = dags{g};
Daniel@0:   for j=1:n
Daniel@0:     if isclamped %DWH
Daniel@0:         u = find(clamped(j,:)==0);    
Daniel@0:     end
Daniel@0:     ps = parents(dag, j);
Daniel@0:     score(g) = score(g) + score_family(j, ps, type{j}, scoring_fn, ns, discrete, data(:,u), params{j});
Daniel@0:   end
Daniel@0: end