function best_trace = renumberAlignment(trace, step_penalties)
%best_trace = renumberAlignment(trace, step_penalties)
%
%A filter to reorder an alignment to optimize stepping.
%
%The self-aligner only maps together levels that it thinks are actually the 
%same. But it doesn't put them in any particular order, just the order in 
%which they first appeared. This finds the right labeling of levels so they 
%progress in the most reasonable way given a series of step counts.
%
%trace is a series of integers expressing the mapping of a measured
%sequence of levels to the self-alignment filtered sequence of levels.
%
%step_penalties are log(p) penalties for [step, skip, back, hold].


lookout_distance = 5; %how many levels out to check for permutations

best_trace = trace;
best_score = totalPenalty(trace, step_penalties);
changeflag = true;

%loop until nothing changes
while changeflag
    changeflag = false;
    
    %loop over all pairs of labvels that are no more than lookout_distance
    %apart
    for ii = 1:(max(best_trace)-1)
        for jj = ii+1:min(ii+lookout_distance,max(best_trace))
            
            %perform the permutation
            candidate_trace = best_trace;
            candidate_trace(best_trace == ii) = jj;
            candidate_trace(best_trace == jj) = ii;
            
            %calculate the new score
            this_score = totalPenalty(candidate_trace, step_penalties);
            
            %if it's better, update the best one and remember to loop again
            if this_score > best_score
                best_trace = candidate_trace;
                best_score = this_score;
                changeflag = true;
            end
            
        end
    end
end


end

function total_penalty = totalPenalty(trace, step_penalties)
%add up all the penalties for step counts

step_sizes = diff(trace);
total_penalty =   step_penalties(1)*sum(step_sizes == 1) ...
    + sum(step_penalties(2)*(step_sizes(step_sizes > 1) - 1)) ...
    + sum(step_penalties(3)*(-step_sizes(step_sizes < 0))) ...
    + step_penalties(4)*sum(step_sizes == 0);



end