Give in to mse

function moves = solver( seq, targ, budget )

moves = ones(budget,4);

% 1. Split sequence and target in budget + 1 segments, the last one can be
%    shorter or longer than the others.
% 2. Compare which of the sequence segments (normal and flipped version)
%    corresponds best with the first one of the target (using mse)
% 3. Reorder seq, moving this segment to the beginning.
% 4. budget -= 1;
% 5. Repeat for the rest of seq and targ.

from = 1;
to = length( seq );
iSegTarg = 1;
iMove = 1;

while budget > 0
    [segPoss, segLens] = split( from, to, budget + 1 );
    minMse = inf;
    targFrom = segPoss( iSegTarg );
    for iSeg = 1:length( segPoss )
        targTo = targFrom + segLens( iSeg ) - 1;
        segFrom = segPoss( iSeg );
        segTo = segFrom + segLens( iSeg) - 1;
        curMse = mse( targ( targFrom:targTo ), seq( segFrom:segTo ) );
        if curMse < minMse
            minMse = curMse;
            posBestSeg = segFrom;
            lenBestSeg = segLens( iSeg );
            flip = false;
        end
        curMse = mse( targ( targFrom:targTo ), fliplr( seq( segFrom:segTo ) ) );
        if curMse < minMse
            minMse = curMse;
            posBestSeg = segFrom;
            lenBestSeg = segLens( iSeg );
            flip = true;
        end
    end
    moves( iMove, : ) = [lenBestSeg, posBestSeg, targFrom, flip];
    seq = splice( seq, lenBestSeg, posBestSeg, targFrom, flip );
    iMove = iMove + 1;
    budget = budget - 1;
    from = from + lenBestSeg;
end

function b = splice(a, len, ai, bi, flip)
%SPLICE  Perform the gene splice
%   b = splice(a, len, ai, bi, flip) where a is the input sequence,
%   len is the length of the extracted fragment, ai is the fragment's
%   starting index in a, bi is the fragment's starting index in b, and
%   flip is a boolean flag (true = flip, false = don't flip)

b = zeros(size(a));

aSpan = ai + (1:len) - 1;
bSpan = aSpan - ai + bi;

if flip
    b(bSpan) = fliplr(a(aSpan));
else
    b(bSpan) = a(aSpan);
end

a(aSpan) = 0;
b(~b) = a(~~a);


function [segPoss, segLens] = split( from, to, parts )

len = to - from + 1;
% sequence can not be split exactly
if mod( len, parts ) ~= 0
    partLen = round( len / parts );
    segPoss = from : partLen : min( from + partLen * parts - 1, to - 1 );
    parts = length( segPoss );
    segLens = partLen * ones( 1, parts );
    segLens( parts ) = segLens( parts ) - (parts*partLen - len);
else
    partLen = len / parts;
    segPoss = from: partLen : to;
    segLens = partLen * ones( 1, parts );
end


function m = mse( a, b )

m = sum( ( a - b ).^2 ) / length( a );