view toolboxes/FullBNT-1.0.7/bnt/examples/dynamic/HHMM/Motif/mk_motif_hhmm.m @ 0:e9a9cd732c1e tip

first hg version after svn

function [bnet, Qnodes, Fnodes, Onode] = mk_motif_hhmm(varargin)
% [bnet, Qnodes, Fnodes, Onode] = mk_motif_hhmm(...)
%
% Make the following HHMM
%
%    S2 <----------------------> S1
%    |                           |
%    |                           |
%   M1 -> M2 -> M3 -> end        B1 -> end
%
% where Mi represents the i'th letter in the motif
% and B is the background state.
% Si chooses between running the motif or the background.
% The Si and B states have self loops (not shown).
%
% The transition params are defined to respect the above topology.
% The background is uniform; each motif state has a random obs. distribution.
%
% Optional params:
% motif_length  - required, unless we specify motif_pattern
% motif_pattern - if specified, we make the motif submodel deterministically
%                  emit this pattern
% background    - if specified, we make the background submodel
%                  deterministically emit this (makes the motif easier to see!)


args = varargin;
nargs = length(args);

% extract pattern, if any
motif_pattern = [];
for i=1:2:nargs
  switch args{i},
   case 'motif_pattern', motif_pattern = args{i+1}; 
  end
end

% set defaults
motif_length = length(motif_pattern);
background_char = [];

% get params
for i=1:2:nargs
  switch args{i},
   case 'motif_length', motif_length = args{i+1}; 
   case 'background', background_char = args{i+1};
  end
end


chars = ['a', 'c', 'g', 't'];
Osize = length(chars);

Qsize = [2 motif_length];
Qnodes = 1:2;
D = 2;
transprob = cell(1,D);
termprob = cell(1,D);
startprob = cell(1,D);

% startprob{d}(k,j), startprob{1}(1,j)
% transprob{d}(i,k,j), transprob{1}(i,j)
% termprob{d}(k,j)


% LEVEL 1

startprob{1} = zeros(1, 2);
startprob{1} = [1 0]; % always start in the background model

% When in the background state, we stay there with high prob
% When in the motif state, we immediately return to the background state.
transprob{1} = [0.8 0.2;
		1.0 0.0];


% LEVEL 2
startprob{2} = 'leftstart'; % both submodels start in substate 1
transprob{2} = zeros(motif_length, 2, motif_length);
termprob{2} = zeros(2, motif_length);

% In the background model, we only use state 1.
transprob{2}(1,1,1) = 1; % self loop
termprob{2}(1,1) = 0.2; % prob transition to end state

% Motif model
transprob{2}(:,2,:) = mk_leftright_transmat(motif_length, 0); % no self loops
termprob{2}(2,end) = 1.0; % last state immediately terminates


% OBS LEVEl

obsprob = zeros([Qsize Osize]);
if isempty(background_char)
  % uniform background model
  %obsprob(1,1,:) = normalise(ones(Osize,1));
  obsprob(1,1,:) = normalise(rand(Osize,1));
else
  % deterministic background model (easy to see!)
  m = find(chars==background_char);
  obsprob(1,1,m) = 1.0;
end

if ~isempty(motif_pattern)
  % initialise with true motif (cheating)
  for i=1:motif_length
    m = find(chars == motif_pattern(i));
    obsprob(2,i,m) = 1.0;
  end
else
  obsprob(2,:,:) = mk_stochastic(rand(motif_length, Osize));
end

if 0
  Oargs = {'CPT', obsprob};
else
  % We use a minent prior for the emission distribution for the states in the motif model
  % (but not the background model). This encourages nearly deterministic distributions.
  % We create an index matrix  (where M = motif length)
  %  [2 1
  %   2 2
  %   ...
  %   2 M]
  % and then convert this to a list of integers, which
  % specifies when to use the minent prior (Q1=2 specifies motif model).
  M = motif_length;
  ndx = [2*ones(M,1) (1:M)'];
  pcases = subv2ind([2 motif_length], ndx);
  Oargs = {'CPT', obsprob, 'prior_type', 'entropic', 'entropic_pcases', pcases};
end



[bnet, Qnodes, Fnodes, Onode] = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ...
	       'Oargs', Oargs, 'Ops', Qnodes(1:2), ...
	       'startprob', startprob, 'transprob', transprob, 'termprob', termprob);