function song = segmentation_auto(song, param, varargin)


%% calculating similarity matrix
% exp of minus the cosine distance

if any(strcmp(fieldnames(song.syncchroma),'wide'))
	v = song.syncchroma.wide;
else
	v = song.syncchroma.treble + song.syncchroma.bass(:,1:12) * 0.8;
end
nBeat = song.nBeat;
flatPCP = std(v')<.000001;
v(flatPCP,1) = 1000; % arbitrary
simmat0 = squareform(pdist(v,'correlation'));
% simmat0(flatPCP,:) = 2;
% simmat0(:,flatPCP) = 2;
% simmat = qnormalise(exp(-simmat0),inf,1);
% simmat(flatPCP,:) = 0;
% simmat(:,flatPCP) = 0;

simmat0(flatPCP,:) = 1;
simmat0(:,flatPCP) = 1;
simmat = 1-simmat0/2;
simmat(flatPCP,:) = 0;
simmat(:,flatPCP) = 0;

%%
median_simmat = zeros(size(simmat));
for  i = 1:nBeat
    temp1 = diag(medfilt1(diag(simmat,i-1),param.seg.medfilt_length,[],1),i-1);
    median_simmat = median_simmat + temp1;
end
median_simmat = median_simmat + median_simmat' - diag(diag(median_simmat));
median_simmat(isnan(median_simmat)) = 0;
%%median_simmat = (median_simmat - mean(median_simmat(:)))/std(median_simmat(:));
%%%%%%%% ATTENTION!!!!!!!! using MEDIAN %%%%%%%%%%%%%
% median_simmat = (median_simmat - median(median_simmat(:)))/std(median_simmat(:));
if param.seg.standardise
    med_median_simmat = repmat(median(median_simmat),nBeat,1);
    std_median_simmat = repmat(std(median_simmat),nBeat,1);
    median_simmat = (median_simmat - med_median_simmat) ./ std_median_simmat;
end
%% get bar boundaries!!!

potential_duplicates = triu(median_simmat > param.seg.thresh_beat);

%%
partlengths = param.seg.minlength:4:param.seg.maxlength;
nPartlengths = length(partlengths);

% initialise arrays
simArray = zeros(nBeat,nBeat,nPartlengths);
decisionArray2 = zeros(nBeat,nBeat,nPartlengths);



fprintf(1,'\n');
fprintf(2,'    ');
for iLength = 1:nPartlengths
    len = partlengths(iLength);
    nUsedBeat = nBeat - len + 1; % number of potential rep beginnings: they can't overlap at the end of the song
    fprintf(1,'\rLength: %1.0d',len)
    fprintf(2,'\b\b\b\b%3.0f%%',iLength/nPartlengths*100);
    for iBeat = 1:nUsedBeat % looping over all columns (arbitrarily chosen columns)
        help2 = find(potential_duplicates(1:nUsedBeat,iBeat));
        for kBeat = help2' % check only potential duplicates
            % measure how well two length len segments go together
            simArray(iBeat,kBeat,iLength) = ...
                quantile(diag(median_simmat(iBeat+(0:len-1),kBeat+(0:len-1))),param.seg.quantile);
        end
    end
    simArray(1:nUsedBeat,1:nUsedBeat,iLength) = simArray(1:nUsedBeat,1:nUsedBeat,iLength) ...
        + simArray(1:nUsedBeat,1:nUsedBeat,iLength)' ...
        - eye(nUsedBeat) .* simArray(1:nUsedBeat,1:nUsedBeat,iLength);

    simArray(:,:,iLength) = conv2(simArray(:,:,iLength),[.01 .98 .01], 'same');
    % take only over-average bars that do not overlap
    for iBeat = 1:nUsedBeat
        temp = simArray(:,iBeat,iLength) > param.seg.thresh_seg;
        decisionArray2(temp,iBeat,iLength) = ...
            simArray(temp,iBeat,iLength);
    end
    decisionArray2(:,:,iLength) = decisionArray2(:,:,iLength) .* (decisionArray2(:,:,iLength) >= maxfilt1(decisionArray2(:,:,iLength),len-1));
    decisionArray2(:,:,iLength) = decisionArray2(:,:,iLength) .* decisionArray2(:,:,iLength)';
%     potential_duplicates = triu(simArray(:,:,iLength) > param.seg.thresh_seg * 0.7);
    potential_duplicates = potential_duplicates .* (simArray(:,:,iLength) > param.seg.thresh_seg);
end
%% milk the data

bestvals = [];
for iLength = 1:nPartlengths
    currLogicSum = sum(decisionArray2(:,:,iLength)>0,2);
    for iBeat = 1:nBeat
        if currLogicSum(iBeat) > 1.0
            currSum = (mean(decisionArray2(decisionArray2(:,iBeat,iLength)>0,iBeat,iLength)))/2;
            bestvals = [bestvals; (currLogicSum(iBeat)-1) * partlengths(iLength), currSum, iLength, iBeat, currLogicSum(iBeat)];
        end
    end
end
%%
%%
% make a table of all valid sets of parts
partletters = {'A','B','C','D','E','F','G', 'H','I','J','K','L','M','N','O','P','Q','R','S'};
valid_sets = ones(1,size(bestvals,1));
parts = [];
if ~isempty(bestvals)
    bestvals(:,2) = bestvals(:,2) / max(bestvals(:,2)*2);
    bestvals(:,1) = bestvals(:,1) + bestvals(:,2);
    bestvals(:,2) = [];

    for kSeg = 1:7
        currbestvals = bestvals(logical(valid_sets),:);
        [m,i] = max(currbestvals(:,1));
        if isempty(m)
            break
        end

        bestLength = partlengths(currbestvals(i,2));
        bestIndices = decisionArray2(currbestvals(i,3),:,currbestvals(i,2));
        islands = conv2(1*(bestIndices > 0), [zeros(1,bestLength-1) ones(1,bestLength)], 'same');
        %%
        if param.seg.newcheck
            iBetter = i;
            bestSubval = 0;
            subm1 = find(currbestvals(:,4) > currbestvals(i,4)); % has to have more part instances
            for kSub = 1:length(subm1)
                iSub = subm1(kSub);
                subIndices = decisionArray2(currbestvals(iSub,3),:,currbestvals(iSub,2));
                subLength = partlengths(currbestvals(iSub,2));
                subIslands = conv2(1*(subIndices > 0), [zeros(1,subLength-1) ones(1,subLength)], 'same');
                if sum(islands.*subIslands) == currbestvals(i,4) * subLength ...
                        && currbestvals(iSub,1) > bestSubval
                    iBetter = iSub;
                    bestSubval = currbestvals(iSub,1);
                end
            end
            i = iBetter;
            bestLength = partlengths(currbestvals(i,2));
            bestIndices = decisionArray2(currbestvals(i,3),:,currbestvals(i,2));
            islands = conv2(1*(bestIndices > 0), [zeros(1,bestLength-1) ones(1,bestLength)], 'same');
        end
        %%
        part.n = bestLength;
        part.indices = find(bestIndices);
        part.letter = partletters{kSeg};
        part.level = kSeg;
        parts = [parts part];

        for iSet = find(valid_sets)
            currislands = conv2(2*(decisionArray2(bestvals(iSet,3),:,bestvals(iSet,2)) > 0), ...
                [zeros(1,partlengths(bestvals(iSet,2))-1) ones(1,partlengths(bestvals(iSet,2)))], 'same');
            if any(islands .* currislands > 0)
                valid_sets(iSet) = false;
            else
                %             stem(islands+currislands*.5)
                %             aaa = 1;
            end
        end
    end
else
    part.n = nBeat;
    part.indices = 1;
    part.letter = 'A';
    part.level = 1;
    parts = [parts part];
end

parts = [parts, nullpart(parts,1:nBeat)];
if param.seg.editor
    [pa, ta] = partarray(parts);
    parts = editorssearch(pa, ta, parts);
    parts = [parts, nullpart(parts,1:nBeat)];
end
parts = mergenulls(parts);

song.parts = parts;
% song.parts


