%HMM Hidden Markov Model adapted for PRtools
%
% [W,R,S,M] = QDC(A,R,S,M)
% W = A*QDC([],R,S)
%
% INPUT
% A Dataset
% R,S Regularization parameters, 0 <= R,S <= 1
% (optional; default: no regularization, i.e. R,S = 0)
% M Dimension of subspace structure in covariance matrix (default: K,
% all dimensions)
%
% OUTPUT
% W Quadratic Bayes Normal Classifier mapping
% R Value of regularization parameter R as used
% S Value of regularization parameter S as used
% M Value of regularization parameter M as used
%
% DESCRIPTION
% Computation of the quadratic classifier between the classes of the dataset
% A assuming normal densities. R and S (0 <= R,S <= 1) are regularization
% parameters used for finding the covariance matrix by
%
% G = (1-R-S)*G + R*diag(diag(G)) + S*mean(diag(G))*eye(size(G,1))
%
% This covariance matrix is then decomposed as G = W*W' + sigma^2 * eye(K),
% where W is a K x M matrix containing the M leading principal components
% and sigma^2 is the mean of the K-M smallest eigenvalues.
%
%
%
% The use of soft labels is supported. The classification A*W is computed by
% NORMAL_MAP.
%
% If R, S or M is NaN the regularisation parameter is optimised by REGOPTC.
% The best result are usually obtained by R = 0, S = NaN, M = [], or by
% R = 0, S = 0, M = NaN (which is for problems of moderate or low dimensionality
% faster). If no regularisation is supplied a pseudo-inverse of the
% covariance matrix is used in case it is close to singular.
%
% EXAMPLES
% See PREX_MCPLOT, PREX_PLOTC.
%
% REFERENCES
% 1. R.O. Duda, P.E. Hart, and D.G. Stork, Pattern classification, 2nd
% edition, John Wiley and Sons, New York, 2001.
% 2. A. Webb, Statistical Pattern Recognition, John Wiley & Sons,
% New York, 2002.
%
% SEE ALSO
% MAPPINGS, DATASETS, REGOPTC, NMC, NMSC, LDC, UDC, QUADRC, NORMAL_MAP
% Copyright: R.P.W. Duin, r.p.w.duin@prtools.org
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands
% $Id: qdc.m,v 1.7 2008/03/20 09:25:10 duin Exp $
function result_out = hmm_prtools(varargin)
prtrace(mfilename);
cEmissionCalcMode = { 'normal' 'perPreviousClass' 'custom'};
% 'normal': One classifier per state or hidden node (e.g. one cov marix and mean vector per state)
% 'perPreviousClass': For a C class classification problem, one classifier given that class i happened (e.g. Given one Supraventricular beat, the same cov matrix for all C classes, and one mean vector for each class)
% 'custom': Custom calculation of cov matrices and mean vectors for each state.
%input parsing
p = inputParser; % Create instance of inputParser class.
p.addOptional('dsAny', [], @(x)(isdataset(x) ) );
p.addOptional('wMapping', [], @(x)(ismapping(x) ) );
p.addParamValue('CantPrevBeats', 0, @(x)(isnumeric(x) && x >= 0 ) );
p.addParamValue('TransMat', [], @(x)(isnumeric(x) && all(all(x >= 0)) ) );
p.addParamValue('InitialStateProb', [], @(x)(isnumeric(x) && all(x >= 0) ) );
p.addParamValue('EmmissionMapping', weighted_ldc([],true), @(x)( ismapping(x) ) );
p.addParamValue('TiedStates', [], @(x)(isnumeric(x)) );
p.addParamValue('Details', false, @(x)(islogical(x)) );
p.addParamValue('PooledCovMat', false, @(x)(islogical(x)) );
p.addParamValue('EmissionEstimationMode', 'normal', @(x)any(strcmpi(x,cEmissionCalcMode)) );
p.parse( varargin{:} );
dsAny = p.Results.dsAny;
wMapping = p.Results.wMapping;
CantPrevBeats = p.Results.CantPrevBeats;
TransMat = p.Results.TransMat;
InitialStateProb = p.Results.InitialStateProb;
wEmmissionMapping = p.Results.EmmissionMapping;
TiedStates = p.Results.TiedStates;
EmissionEstimationMode = p.Results.EmissionEstimationMode;
bDetails = p.Results.Details;
bPooledCovMat = p.Results.PooledCovMat;
if( nargin < 1 || isempty(dsAny) )
%% No input arguments:
result_out = mapping(mfilename,varargin); % return an untrained mapping.
result_out = setname(result_out,'HMM');
elseif( ~isempty(wMapping) && istrained(wMapping) && isdataset(dsAny) )
%% mapping
m = getsize(dsAny,1);
[k c] = getsize(wMapping);
result_out = zeros(m,c);
[ obs_labels, ~, dsAny ] = ds_convert_prtools_pmtk(dsAny);
recording_index = getident(dsAny, 'recording_index');
% hare una observacion de cada registro
[recordings, ~, rec_idx] = unique(recording_index);
cant_recordings = length(recordings);
model = +wMapping;
predicted_seq = cell(cant_recordings,1);
lStateGroups = length(model.TiedStatesGroups);
for ii = 1:cant_recordings
each_rec_idx = find( ii == rec_idx );
logB = nan(model.hmm_states, length(each_rec_idx) );
%protect from NaNs
aux = +dsAny(each_rec_idx,:);
aux_mean = nanmean(aux);
aux_nan_idx = find(isnan(aux));
[aux_nan_row aux_nan_col] = ind2sub( size(aux), aux_nan_idx);
for jj = 1:length(aux_nan_row)
aux(aux_nan_row(jj), aux_nan_col(jj)) = aux_mean(aux_nan_col(jj));
end
dsAny(each_rec_idx,:) = aux;
clear aux
%END protect from NaNs
for jj = 1:lStateGroups
this_states = find(model.TiedStatesGroups_idx == jj);
switch(model.EmissionEstimationMode)
case 'normal'
logB = dsAny(each_rec_idx,:) * model.emission;
% logB = normalize( (+logB)', 1);
logB = (+logB)';
case 'perPreviousClass'
aux = dsAny(each_rec_idx,:) * model.emission{jj};
% logB(this_states,:) = normalize( (+aux)', 1);
logB(this_states,:) = (+aux)';
case 'custom'
aux = dsAny(each_rec_idx,:) * model.emission{jj};
% aux = normalize( (+aux)', 1);
aux = (+aux)';
tied_lobB = size(aux, 1);
tied_classes = length(this_states);
if( tied_lobB == tied_classes )
logB(this_states,:) = aux;
else
logB(this_states,:) = repmat(aux, tied_classes/tied_lobB, 1);
end
otherwise
error('Unrecognized case.');
end
end
[predicted_seq{ii}, ~, ~] = hmmViterbiC(log(model.pi+eps), log(model.A+eps), log(logB));
% [predicted_seq{ii}, ~, ~] = hmmViterbiCM(log(model.pi+eps), log(model.A+eps), log(logB));
% % conv_idx = [4 1 2 4 3];
% figure(99)
% % plot(conv_idx(obs_labels{ii}), 'bo-' )
% plot(obs_labels{ii}, 'bo-' )
% hold on
% ps = unchange_labels(predicted_seq{ii}, model.CantPrevBeats+1, c);
% ps2 = unchange_labels(predicted_seq2{ii}, model.CantPrevBeats+1, c);
% plot(ps, 'rx-' )
% plot(ps2, 'gs-' )
% hold off
end
predicted_seq = cell2mat(predicted_seq')';
predicted_changed = predicted_seq;
predicted_seq = unchange_labels(predicted_changed, model.CantPrevBeats+1, c);
if(model.bDetails)
true_labels = getnlab(dsAny);
%take class strings initials letters
learned_lablist = getlabels(wMapping);
test_lablist = getlablist(dsAny);
orig_lablist = (learned_lablist(:,1));
lablist = orig_lablist';
beat_sequences = orig_lablist;
cant_labels = length(lablist);
[~, test_lablist_idx, orig_lablist_idx] = intersect(cellstr(test_lablist), cellstr(learned_lablist) );
test_cant_labels = size(test_lablist,1);
conversion_idx = repmat(c+1,test_cant_labels,1);
if( ~isempty(test_lablist_idx) )
conversion_idx(test_lablist_idx) = orig_lablist_idx;
end
if( test_cant_labels > c )
conv_lablist = [ orig_lablist; '~'];
% warning('More classes than trained in evaluation dataset, "~" means not learned classes.')
else
conv_lablist = orig_lablist;
end
conv_lablist = conv_lablist(:);
% descarto las clases no aprendidas
true_labels = conversion_idx(true_labels);
learned_examples_idx = find(true_labels<=c);
%Generate those labels for each state
dsAux = setlabels(dsAny, true_labels);
[dsAux selected_idx] = seldat(dsAux, 1:3);
test_true_labels = ds_convert_prtools_pmtk(dsAux);
true_states_labels = change_labels(test_true_labels, model.CantPrevBeats+1, c);
true_states_labels = cell2mat(true_states_labels')';
true_labels = true_labels(learned_examples_idx);
predicted_seq_str = conv_lablist(predicted_seq(learned_examples_idx));
true_labels_str = conv_lablist( true_labels );
%Generate a lablist of strings for each hidden state
for ii = 2:(model.CantPrevBeats+1)
aux = repmat(lablist,size(beat_sequences,1),1);
beat_sequences = [ aux(:) repmat(beat_sequences,cant_labels,1) ];
end
lablist = beat_sequences;
true_states_labels_str = lablist(true_states_labels,:);
predicted_changed = lablist(predicted_changed(selected_idx),:);
%Displays confmat of hidden states
confmat(true_states_labels_str, predicted_changed);
%Generate filters to group by previous class sequence history,
%and to measure the relative importance of each group in the
%evaluation dataset.
cant_labels_in = [];
labels_in = true(size(true_labels));
jj = 1;
for ii = 1:c:model.hmm_states
labels_in(:,jj) = true_states_labels >= ii & true_states_labels < (ii+c);
cant_labels_in(jj) = sum(labels_in(:,jj));
jj = jj + 1;
end
[~, importance_idx] = sort(cant_labels_in, 'descend');
c_aux = [];
jj = 1;
for ii = importance_idx
prev_beats(jj,:) = true_states_labels_str( find(labels_in(:,ii), 1,'first') , 1:model.CantPrevBeats );
disp(prev_beats(jj,:))
confmat(true_labels_str(labels_in(:,ii)), predicted_seq_str(labels_in(:,ii)) );
[ConfusionMat,NE,lab_list_en_register_out,lab_list_en_trained_classifier ] = confmat(true_labels_str(labels_in(:,ii)), predicted_seq_str(labels_in(:,ii)) );
% Tengo que compatibilizar la matriz de confusion resultante ya que
% habra distintas clases en register_out y trained_classifier.
[NLAB11,NLAB12] = renumlab(orig_lablist, lab_list_en_register_out );
[NLAB21,NLAB22] = renumlab(orig_lablist, lab_list_en_trained_classifier );
cc_aux = zeros(c);
cc_aux(NLAB12,NLAB22) = ConfusionMat;
c_aux(:,:,jj) = cc_aux;
jj = jj + 1;
end
c_aux = [ c_aux; sum(c_aux) ];
c_aux = [ c_aux sum(c_aux,2) ];
for ii = 1:c
disp(orig_lablist(ii))
importance_i = squeeze( round( c_aux(ii,end,:) * 1/sum(c_aux(ii,end,:)) * 100));
sens_i = squeeze(round( c_aux(ii,ii,:) ./ c_aux(ii,end,:) * 100));
pospred_i = squeeze(round(c_aux(ii,ii,:) ./ c_aux(end,ii,:) * 100));
aux = [importance_i sens_i pospred_i]';
[dummy, importance_idx] = sort(importance_i', 'descend');
aux_str = [prev_beats(importance_idx,:) repmat(' ', size(prev_beats,1) , 1)]';
fprintf(1, ' %s\n', aux_str(:) )
disp([['importance '; ...
'sensitiv '; ...
'pospred '] num2str(aux(:,importance_idx)) ] );
fprintf(1, '\n' )
fprintf(1, 'Total Se: %d\n', round( sum(c_aux(ii,ii,:)) / sum(c_aux(ii,end,:)) * 100) )
fprintf(1, 'Total +P: %d\n\n', round( sum(c_aux(ii,ii,:)) / sum(c_aux(end,ii,:)) * 100) )
end
end
%Simulates a soft output, given that viterbi algorithm does not
%output probabilities.
result_out( sub2ind(size(result_out), 1:size(result_out,1), predicted_seq' ) ) = 1;
%build the result dataset.
result_out = setdata(dsAny, result_out, getlabels(wMapping));
elseif( isdataset(dsAny) )
%% training
islabtype(dsAny,'crisp'); % Assert A has the right labtype.
isvaldfile(dsAny,2,2); % at least 2 objects per class, 2 classes
[m,k,c] = getsize(dsAny);
hmm_states = c ^ (CantPrevBeats+1);
model.hmm_states = hmm_states;
model.CantPrevBeats = CantPrevBeats;
model.bDetails = bDetails;
[obs_labels_train, dummy, dsAny] = ds_convert_prtools_pmtk(dsAny);
obs_labels_train = change_labels(obs_labels_train, CantPrevBeats+1, c);
obs_labels_train_stacked = cell2mat(obs_labels_train')';
%% transition matrix
if( isempty(TransMat) )
TransMat = countTransitions(obs_labels_train, hmm_states );
end
if( bDetails )
disp('Estimated Transition Matrix')
disp(TransMat)
end
TransMat(TransMat == 0) = 1;
% TransMat(:,1:c:hmm_states) = 0.5*TransMat(:,1:c:hmm_states);
% TransMat(:,2:c:hmm_states) = 1.5*TransMat(:,2:c:hmm_states);
% TransMat(:,3:c:hmm_states) = 2*TransMat(:,3:c:hmm_states);
model.A = normalize(TransMat, 2);
%% initial state prior
if( isempty(InitialStateProb) )
seqidx = cumsum([1, cellfun(@(seq)size(seq, 2), obs_labels_train')]);
InitialStateProb = rowvec(histc(obs_labels_train_stacked(seqidx(1:end-1)), 1:hmm_states));
end
if( bDetails )
disp('Estimated Initial State Probability')
disp(InitialStateProb)
end
model.pi = normalize(InitialStateProb);
%% state emission distributions estimation
switch(EmissionEstimationMode)
case 'normal'
TiedStates = ones(hmm_states,1);
[StateGroups, dummy, StateGroups_idx ] = unique(TiedStates);
train_lablist = getlablist(dsAny);
aux_prior = getprior(dsAny);
new_prior = nan(1,hmm_states);
for ii = 1:c
new_prior(ii:c:hmm_states) = aux_prior(ii);
end
%Add new labels, one per state, to estimate one emission
%per state.
dsAny = addlabels(dsAny, obs_labels_train_stacked, 'dummy_lab');
new_prior = new_prior *1/(sum(new_prior));
dsAny = setprior(dsAny, new_prior);
case 'perPreviousClass'
TiedStates = repmat(1:(c^CantPrevBeats),c,1);
TiedStates = (TiedStates(:))';
[StateGroups, ~, StateGroups_idx ] = unique(TiedStates);
case 'custom'
%Check for TiedStates
if( ~isempty(TiedStates) )
if( length(TiedStates) ~= hmm_states )
error(['TiedStates is a numeric vector of length ' num2str(hmm_states) ])
end
[StateGroups, ~, StateGroups_idx ] = unique(TiedStates);
else
warning( 'TiedStates option not defined. Assuming "full" mode.' )
StateGroups = 1:hmm_states;
StateGroups_idx = StateGroups;
end
otherwise
error('Unrecognized case.');
end
model.TiedStatesGroups = StateGroups;
model.TiedStatesGroups_idx = StateGroups_idx;
model.EmissionEstimationMode = EmissionEstimationMode;
lStateGroups = length(StateGroups);
model.emission = cell(lStateGroups,1);
if( bDetails )
disp('Estimated Emisssion Probabilities')
disp('---------------------------------')
end
if( bPooledCovMat )
wAux = weighted_ldc( dsAny, true);
wAux = +wAux;
PooledCovMat = wAux.cov;
PooledCovMatDets = wAux.det;
end
for ii = 1:lStateGroups
if( strcmpi(EmissionEstimationMode, 'normal') )
examples_per_state_filter = true(size(obs_labels_train_stacked));
else
examples_per_state_filter = false(size(obs_labels_train_stacked));
this_states = rowvec(find(StateGroups_idx == ii));
for jj = this_states
examples_per_state_filter = examples_per_state_filter | obs_labels_train_stacked == jj;
end
end
%Actual state emission estimation
wMapping = dsAny(examples_per_state_filter,:) * wEmmissionMapping;
if( bPooledCovMat )
if( strcmpi(EmissionEstimationMode, 'custom') )
emis_aux = +wMapping;
emis_aux.cov = PooledCovMat;
emis_aux.det = PooledCovMatDets;
wMapping = setdata(wMapping, emis_aux);
else
warning('Ignoring PooledCovMat parameter');
end
end
if( bDetails )
disp(['State ' num2str(ii) ])
disp('Data: ' )
dsAny(examples_per_state_filter,:)
emis_aux = +wMapping;
disp('Emission mean: ' )
disp(emis_aux.mean)
disp('Emission cov: ' )
disp(emis_aux.cov)
end
switch(EmissionEstimationMode)
case 'normal'
model.emission = wMapping;
case 'perPreviousClass'
model.emission{ii} = wMapping;
case 'custom'
model.emission{ii} = wMapping;
otherwise
error('Unrecognized case.');
end
end
if( strcmpi(EmissionEstimationMode, 'normal') )
result_out = mapping(mfilename, 'trained', model, train_lablist, k, c);
else
result_out = mapping(mfilename, 'trained', model, getlablist(dsAny), k, c);
end
else
error('Illegal call') % this should not happen
end
return;
function obs_labels_train = change_labels(obs_labels_train, hmm_mem, cant_clases)
%asumo el mismo label a los hmm_mem primeros latidos.
cant_patients = length(obs_labels_train);
for ii = 1:cant_patients
aux = obs_labels_train{ii};
obs_labels_train{ii} = [repmat(aux(1),1,hmm_mem) aux];
end
first_idx = cumsum([1; cellfun( @length, obs_labels_train)]);
last_idx = cumsum([0; cellfun( @length, obs_labels_train)]);
first_idx = first_idx(1:end-1);
last_idx = last_idx(2:end);
obs_labels_train_orig = cell2mat(obs_labels_train')';
obs_labels_train = nan(size(obs_labels_train_orig));
all_labels = (1:cant_clases)';
cant_labels = length(all_labels);
beat_sequences = all_labels;
for ii = 2:hmm_mem
% beat_sequences = [ [repmat(1,3,1);repmat(2,3,1);repmat(3,3,1)] repmat((1:3)',3,1) ];
aux = repmat(all_labels',size(beat_sequences,1),1);
beat_sequences = [ aux(:) repmat(beat_sequences,cant_labels,1) ];
end
seq_length = size(beat_sequences,2);
first_beats = repmat(first_idx(:)',hmm_mem,1) + repmat((0:(hmm_mem-1))', 1, size(first_idx,1));
first_beats = first_beats(:);
ii = 1;
for jj = 1:size(beat_sequences,1)
beat_sequence = beat_sequences(jj,:);
Beats_idx = find(obs_labels_train_orig == beat_sequence(seq_length));
Beats_idx = setdiff(Beats_idx, first_beats);
bAux = true(length(Beats_idx),1);
for kk = 1:seq_length-1
bAux = bAux & obs_labels_train_orig(Beats_idx-kk) == beat_sequence(seq_length-kk);
end
Beats_idx = Beats_idx(bAux);
obs_labels_train(Beats_idx) = ii;
ii = ii + 1;
end
cant_patients = length(first_idx);
aux = cell(cant_patients,1);
for ii = 1:cant_patients
aux{ii} = obs_labels_train(first_idx(ii)+hmm_mem: last_idx(ii))';
end
obs_labels_train = aux;
function new_lab = unchange_labels(old_lab, hmm_mem, cant_clases)
new_lab = nan(size(old_lab));
for ii = 1:cant_clases
for jj = ii:cant_clases:cant_clases^hmm_mem
new_lab(jj == old_lab) = ii;
end
end
function [ obs_labels, observations, dsAny ] = ds_convert_prtools_pmtk(dsAny)
labels = getnlab(dsAny);
data = +dsAny;
% patient_index = getident(dsAny, 'patient_index');
recording_index = getident(dsAny, 'recording_index');
beat_time = getident(dsAny, 'beat_time');
user = getuser(dsAny);
ResultsFeatlab = getfeatlab(dsAny);
% hare una observacion de cada registro
[recordings, dummy, rec_idx] = unique(recording_index);
cant_recordings = length(recordings);
observations = cell(cant_recordings,1);
obs_labels = cell(cant_recordings,1);
for ii = 1:cant_recordings
each_rec_idx = find( ii == rec_idx );
% bNanObservations = any(isnan(data(each_rec_idx,:)),2);
%
% %para no perder la sequencia cambiamos los nan por las medianas en cada feature por cada registro.
% if( sum(bNanObservations) > 0 )
%
% this_median = nanmedian(data(each_rec_idx,:));
%
% for jj = find(any(isnan(data(each_rec_idx,:))))
% thisNanObservations_idx = find(isnan(data(each_rec_idx,jj)));
% data(each_rec_idx(thisNanObservations_idx),jj) = this_median(jj);
% end
%
% end
%la ordeno temporalmente
[dummy, time_idx] = sort(beat_time(each_rec_idx));
observations{ii} = data(each_rec_idx(time_idx),:)';
obs_labels{ii} = labels(each_rec_idx(time_idx),:)';
if( nargout > 2 )
%saco el dataset ordenado temporalmente.
dsAny(each_rec_idx,:) = dsAny(each_rec_idx(time_idx),:);
end
end