function [path,field,trans,decision]=pebls091(target,tokens,c,thresh,p)

[trans,field]=gentransnet(tokens,c,thresh);

target=[ones(13,1)*1e4,target,ones(13,1)*2e4];
[r,M]=size(target);
N=size(field,2);

s=zeros(N,2);
frame1=target(:,1);frame1=frame1(:,ones(N,1));
s(:,1)=exp(-c*sqrt(sum((field(1:r,:)-frame1).^2,1))');
s=(s>thresh).*s;
clear frame1;
S=s;
decision=zeros(N,M);

% Compute cumulative distance matrix
for m=2:M
   tic
   for n=1:N
      % DP step
      s(n,2)=exp(-c*sqrt(sum((field(1:r,n)-target(:,m)).^2)));
      gt=S(:,1).*trans(:,n);
      % instead of selecting max, select the mean of the 'best' cluster
      % where cluster score is ((1+mean)^p)*N/(1+variance).
      [increment,decision(n,m)]=best(gt,p,n,m);
      S(n,2)=s(n,2)+increment;
   end
   s(:,1)=s(:,2);
   S(:,1)=S(:,2);
end

% traceback
path(m)=N;  %start from the stop state
while m>=2
   path(m-1)=decision(path(m),m);
   m=m-1;
end
path=path(2:end-1); % exclude start, stop states

% subroutine 'best'
function [increment,decision]=best(gt,p,n,m)
gts=nonzeros(sort(gt));
if nnz(gts)>1
   C=hclust(gts);
   score=((1+C(:,1)).^p).*C(:,3)./(1+C(:,2)); 
   [o,bestclust]=max(score);
   increment=C(bestclust,1);
   [o,decision]=min(abs(gt-increment));
else
   [increment,decision]=max(gt);
end