function [fsm, file] = read_htk_slf(filename)
%  fsm = read_htk_slf(filename)
%
%
% This function is broken as of 2006-12-14

NULL = '!NULL';

% Read the M-file into a cell array of strings: 
[fid, message] = fopen(filename, 'rt');
warning(message)
file = textscan(fid, '%s', 'delimiter', '\n', 'whitespace', '', 'bufSize', 16000);
fclose(fid);

file = file{1};
% Remove any empty lines
file = file(cellfun('length', file) > 0);

%wdnet = strvcat(file);
%wdnet = strcat(file{:});
wdnet = file;

tok = regexpi(wdnet, '^VERSION=(\S+)', 'tokens', 'once');
idx = find(~cellfun('isempty', tok));
if ~strcmp(tok{idx}{1}, '1.0')
  error('This function only supports HTK SLF files version 1.0');
end

% get number of states
tok = regexpi(wdnet, '^\s*N=(\d+)\s+L=(\d+)', 'tokens', 'once');
idx = find(~cellfun('isempty', tok));
nstates = str2num(tok{idx}{1});
narcs = str2num(tok{idx}{2});

% get state labels
names = regexpi(wdnet, '^I=(?<index>\d+)\s+W=(?<label>\S+)', 'names');
names = names(~cellfun('isempty', names));
names = [names{:}];
labels = cellstr(strvcat(names.label));
null_states = strmatch(NULL, strvcat(labels));

% get transition matrix states
transmat = regexpi(wdnet,  '^J=(?<arc_num>\d+)\s+S=(?<start>\d+)\s+E=(?<end>\d)', 'names');
transmat = transmat(~cellfun('isempty', transmat));
transmat = [transmat{:}];

fsm.transmat = zeros(nstates, nstates);
fsm.start_prob = zeros(nstates, 1);
fsm.end_prob = zeros(nstates, 1);
for x = 1:narcs
  t = transmat(x);
  t.start = str2num(t.start) + 1;
  t.end = str2num(t.end) + 1;
  t.arc_num = str2num(t.arc_num) + 1;

  fsm.transmat(t.start+1, t.end+1) = 1;
end
fsm.nstates = nstates;

% null states with no incoming arcs (and thus only outgoing arcs)
% are start states
% if strcmp(labels(t.start), NULL)
%   disp(['start state: ' labels(t.start), ' -> ', labels(t.end), ...
%         '  (' t.arc_num ')'])
%   fsm.start_prob(t.end+1) = 1;
%   % null states with no outgoing arcs are end states
% elseif strcmp(labels(t.end), NULL)
%   disp(['end state: ' labels(t.start), ' -> ', labels(t.end), ...
%         '  (' t.arc_num ')'])
%   fsm.end_prob(t.start+1) = 1;
% else
%   disp(['transmat: ' labels(t.start), ' -> ', labels(t.end), ...
%         '  (' t.arc_num ')'])

% need to toss away non-emitting null states
%idx = setdiff([1:nstates], null_states);
%labels = labels(idx);
%fsm.transmat = fsm.transmat(idx, idx);
%fsm.start_prob = fsm.start_prob(idx);
%fsm.end_prob = fsm.end_prob(idx);

%fsm.nstates = length(idx);
fsm.labels = labels;