function [Incipits,Starts,Ends,Beats] = make_incipits(ENAstruct,NBeat,NIncip)
% [Incipits,Starts,Ends,Beats] = make_incipits(ENAstruct,NBeat,NIncip)
%    Return a set of "Incipits" - beat-chroma patches that
%    represent a full piece.  Algorithm is to take the Bar closest 
%    to each section marker, then take the next NBeat beats
%    (default 32) as the incipit.  Do this for all Sections; if
%    there are more than NIncip (default 4) of them, repeatedly
%    find the two most similar (by Euclidean distance), and remove
%    the chronologically later of the two.
%    Return NIncips x 12 x NBeat array of Incipits, and vectors of
%    the start and end times (in secs) of each selected incipit, as
%    well as indices of starting beats.
% 2010-04-13 Dan Ellis dpwe@ee.columbia.edu

% Make sure args 2 and 3 are not empty
if nargin < 2;  NBeat = 32; end
if nargin < 3;  NIncip = 4; end

% key-normalized beat-chroma
BChroma = chromrot(en_beatchroma(ENAstruct),ENAstruct.key);

nchr = size(BChroma,1);

% Convert Bars into Beat indices, and Sections into Bar indices
BarIndex = find_nearest(ENAstruct.bar, ENAstruct.beat);
SectBarIndex = find_nearest(ENAstruct.section, ENAstruct.bar);
% Map SectBarIndex to beats
SectBeatIndex = BarIndex(SectBarIndex);

nsect = length(SectBeatIndex);
nbeats = length(ENAstruct.beat);

% Output always has NIncip entries, some may be zero
Incipits = zeros(max(NIncip,nsect),nchr,NBeat);
Starts = zeros(1,max(NIncip,nsect));
Ends = zeros(1,max(NIncip,nsect));
Beats = zeros(1,max(NIncip,nsect));
for i = 1:nsect
  nbts = min(NBeat, nbeats-SectBeatIndex(i)+1);
  Incipits(i,:,1:nbts) = BChroma(:,SectBeatIndex(i)-1+[1:nbts]);
  Starts(i) = ENAstruct.beat(SectBeatIndex(i));
  Ends(i) = ENAstruct.beat(SectBeatIndex(i)+nbts-1);
  Beats(i) = SectBeatIndex(i);
end

% Now, maybe prune incipits
while size(Incipits,1) > NIncip
  % Make distance matrix
  currnincip = size(Incipits,1);
  dists = zeros(currnincip,currnincip);
  mindist = Inf;
  for i = 1:currnincip-1
    for j = i+1:currnincip
      dist = norm(squeeze(Incipits(i,:,:)-Incipits(j,:,:)));
%      disp(sprintf('%d - %d = %f',i,j,dist));
      if dist < mindist
        mindist = dist;
        minpair = [i j];
      end
    end
  end
%  disp(sprintf('min: %d - %d = %f',minpair(1),minpair(2),mindist));
  % Now remove later of min pair
  inkeep = find([1:currnincip] ~= minpair(2));
  Incipits = Incipits(inkeep,:,:);
  Starts = Starts(inkeep);
  Ends = Ends(inkeep);
  Beats = Beats(inkeep);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function I = find_nearest(times, base)
%  Return I such that base(I) = times, or as close as possible
nt = length(times);
nb = length(base);

diffs = repmat(times,nb,1) - repmat(base',1,nt);

[vv,I] = min(abs(diffs));