- function moves = solver(sequence,target,budget)
+ function moves = solver(ss,tt,bb)
last = 2;
- if budget > last
- moves = oldconvolutedsolver(sequence,target,budget-last);
- moves2 = new2solver(doMoves(sequence,moves,length(sequence)),target,budget-size(moves,1));
- score2 = findscore(doMoves(sequence,moves,length(sequence)),target,budget-size(moves,1),moves2);
- moves3 = anothersolver(doMoves(sequence,moves,length(sequence)),target,budget-size(moves,1));
- score3 = findscore(doMoves(sequence,moves,length(sequence)),target,budget-size(moves,1),moves3);
- if score3<score2,
- moves2 = moves3;
- end
- moves = [moves; moves2];
+ if bb > last
+ moves = oldconvolutedsolver(ss, tt, bb-last);
+ sequence2 = doMoves(ss, moves, length(ss));
+ bb = bb - size(moves,1);
+ moves2 = new2solver(sequence2, tt, bb);
+ score2 = findscore(sequence2, tt, bb, moves2);
+ moves3 = anothersolver(sequence2, tt, bb);
+ score3 = findscore(sequence2, tt, bb, moves3);
+ if score3 < score2
+ moves = [moves; moves3];
else
- moves = markussolver(sequence,target,budget);
+ moves = [moves; moves2];
end
+ else
+ moves = markussolver(ss, tt, bb);
+ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function [moves,score] = swapper(sequence,target,budget,rfac)
+ function [moves,score] = swapper(ss,tt,bb,rfac)
- moves = zeros(budget,4);
- n=numel(sequence);
- runningScore=zeros(budget,1);
+ moves = zeros(bb,4);
+ n=numel(ss);
+ runningScore=zeros(bb,1);
- targrows = target(ones(n,1),:);
+ targrows = repmat(tt,n,1);
targcols=targrows';
- seqrows = sequence(ones(n,1),:);
+ seqrows = repmat(ss,n,1);
altered=false(n,1);
i=1;
- while i<budget
- seqrows(:,altered) = repmat(sequence(altered),n,1);
- seqcols = seqrows';
- % Delta(i,j) = score increase resulting from swap of i with j
- Delta = abs(targcols-seqrows) + abs(seqcols-targrows) - abs(targrows-seqrows) - abs(targcols-seqcols) ;
- if rfac>0
- Delta = Delta.*(1 + rfac*rand(n));
- end
- [dmin,i1] = min(Delta);
- [dmin,j1] = min(dmin);
- i1=i1(j1);
- ij =[i1 j1]; i1 = min(ij); j1=max(ij);
+ while i<bb
+ seqrows(:,altered) = repmat(ss(altered),n,1);
+ seqcols = seqrows';
+ % Delta(i,j) = score increase resulting from swap of i with j
+ Delta = abs(targcols-seqrows) + abs(seqcols-targrows) - abs(targrows-seqrows) - abs(targcols-seqcols) ;
+ if rfac>0
+ Delta = Delta.*(1 + rfac*rand(n));
+ end
+ [dmin,i1] = min(Delta);
+ [dmin,j1] = min(dmin);
+ i1=i1(j1);
+ ij =[i1 j1]; i1 = min(ij); j1=max(ij);
- % now see if it pays to include neighbours so that we swap a band i1 to i2 with a band j1 to j2
- dcum1 = Delta(i1,j1);
- i2=i1;j2=j1;
- % go R (i2>i1) first:
- keepgoing = (j2<n) * (i2<j1-1);
- while keepgoing
- d = Delta(i2+1,j2+1);
- if d<0
- i2=i2+1;
- j2=j2+1;
- dcum1 = dcum1+d; % increment cumulative delta
- end
- keepgoing = ~( (d>=0) + (j2==n) + (i2>=j1-1) );
- end
- % go L (decreasing i1) :
- keepgoing = (i1>1 )*(i2<j1-1);
- while keepgoing
- d = Delta(i1-1,j1-1);
- if d<0
- i1=i1-1;
- j1=j1-1;
- dcum1 = dcum1+d; % increment cumulative delta
- end
- keepgoing= ~((d>=0) + (i1==1) + (i2>=j1-1) );
- end
+ % now see if it pays to include neighbours so that we swap a band i1 to i2 with a band j1 to j2
+ dcum1 = Delta(i1,j1);
+ i2=i1;j2=j1;
+ % go R (i2>i1) first:
+ keepgoing = (j2<n) * (i2<j1-1);
+ while keepgoing
+ d = Delta(i2+1,j2+1);
+ if d<0
+ i2=i2+1;
+ j2=j2+1;
+ dcum1 = dcum1+d; % increment cumulative delta
+ end
+ keepgoing = ~( (d>=0) + (j2==n) + (i2>=j1-1) );
+ end
+ % go L (decreasing i1) :
+ keepgoing = (i1>1 )*(i2<j1-1);
+ while keepgoing
+ d = Delta(i1-1,j1-1);
+ if d<0
+ i1=i1-1;
+ j1=j1-1;
+ dcum1 = dcum1+d; % increment cumulative delta
+ end
+ keepgoing= ~((d>=0) + (i1==1) + (i2>=j1-1) );
+ end
- L = i2-i1+1;
- moves1 = [L j1 i1 0; L i1+L j1 0];
+ L = i2-i1+1;
+ moves1 = [L j1 i1 0; L i1+L j1 0];
- % now repeat the exercise for a flipped sequence
- dcum2 = Delta(i1,j1);
- i2=i1;j2=j1;
- % go R (increasing i2) first:
- keepgoing = i2<j1-1;
- while keepgoing
- % d = swapdelta(sequence,target,i2+1,j1-1);
- d = Delta(i2+1,j1-1);
- if d<0
- i2=i2+1;
- j1=j1-1;
- dcum2 = dcum2+d; % increment cumulative delta
- end
- keepgoing = ~( (d>=0) + (i2>=j1-1) );
- end
- % go L (decreasing i1) :
- keepgoing = (i1>1)*(i2<j1-1)*( j2<n);
- while keepgoing
- % d = swapdelta(sequence,target,i1-1,j2+1);
- d = Delta(i1-1,j2+1);
- if d<0
- i1=i1-1;
- j2=j2+1;
- dcum2 = dcum2+d; % increment cumulative delta
- end
- keepgoing = ~((d>=0) + (i1==1) + (j2==n));
- end
- L = i2-i1+1;
- moves2 = [L j1 i1 1; L i1+L j1 1];
+ % now repeat the exercise for a flipped ss
+ dcum2 = Delta(i1,j1);
+ i2=i1;j2=j1;
+ % go R (increasing i2) first:
+ keepgoing = i2<j1-1;
+ while keepgoing
+ % d = swapdelta(ss,tt,i2+1,j1-1);
+ d = Delta(i2+1,j1-1);
+ if d<0
+ i2=i2+1;
+ j1=j1-1;
+ dcum2 = dcum2+d; % increment cumulative delta
+ end
+ keepgoing = ~( (d>=0) + (i2>=j1-1) );
+ end
+ % go L (decreasing i1) :
+ keepgoing = (i1>1)*(i2<j1-1)*( j2<n);
+ while keepgoing
+ % d = swapdelta(ss,tt,i1-1,j2+1);
+ d = Delta(i1-1,j2+1);
+ if d<0
+ i1=i1-1;
+ j2=j2+1;
+ dcum2 = dcum2+d; % increment cumulative delta
+ end
+ keepgoing = ~((d>=0) + (i1==1) + (j2==n));
+ end
+ L = i2-i1+1;
+ moves2 = [L j1 i1 1; L i1+L j1 1];
- if dcum1<=dcum2
- moves(i:i+1,:) = moves1;
- else
- moves(i:i+1,:) = moves2;
- end
- L= moves(i,1); j1 = moves(i,2); i1=moves(i,3); % note the affected genes
+ if dcum1<=dcum2
+ moves(i:i+1,:) = moves1;
+ else
+ moves(i:i+1,:) = moves2;
+ end
+ L= moves(i,1); j1 = moves(i,2); i1=moves(i,3); % note the affected genes
- sequence1(1,:) = doMoves(sequence,moves(i,:),n);
- score1(1) = sum(abs(sequence1(1,:)-target));
- sequence1(2,:) = doMoves(sequence1(1,:),moves(i+1,:),n);
- score1(2) = sum(abs(sequence1(2,:)-target));
- sequence1(3,:) = doMoves(sequence,moves(i+1,:),n);
- score1(3) = sum(abs(sequence1(3,:)-target));
+ sequence1(1,:) = doMoves(ss,moves(i,:),n);
+ score1(1) = sum(abs(sequence1(1,:)-tt));
+ sequence1(2,:) = doMoves(sequence1(1,:),moves(i+1,:),n);
+ score1(2) = sum(abs(sequence1(2,:)-tt));
+ sequence1(3,:) = doMoves(ss,moves(i+1,:),n);
+ score1(3) = sum(abs(sequence1(3,:)-tt));
- [minscore1,ims] = min(score1);
- sequence=sequence1(ims,:);
- runningScore(i:i+1) = score1(1:2);
+ [minscore1,ims] = min(score1);
+ ss=sequence1(ims,:);
+ runningScore(i:i+1) = score1(1:2);
- if ims==3
- moves(i,:) = moves(i+1,:);
- runningScore(i) = score1(3);
- i=i+1;
- else
- i=i+ims;
- end
+ if ims==3
+ moves(i,:) = moves(i+1,:);
+ runningScore(i) = score1(3);
+ i=i+1;
+ else
+ i=i+ims;
+ end
- % flag the altered bands in the sequence:
- altered = false(1,n);
- altered(i1:i1+L-1) = true;
- altered(j1:j1+L-1) = true;
+ % flag the altered bands in the ss:
+ altered = false(1,n);
+ altered(i1:i1+L-1) = true;
+ altered(j1:j1+L-1) = true;
end
- if i==budget
- s = abs(sequence-target);
- [maxs,imaxs]=max(s); % worst current gene
- [mins,imins1] = min(abs(target-sequence(imaxs))) ; % best place to put the worst gene
- move1 = [1 imaxs imins1 false];
- seq1 = doMoves(sequence,move1,n);
- s1= sum(abs(seq1-target));
- [mins,imins2] = min(abs(target(imaxs)-sequence)) ; % best gene to put in the worst place
- move2 = [1 imins2 imaxs false];
- seq2 = doMoves(sequence,move2,n);
- s2= sum(abs(seq2-target));
- if s1<s2
- moves(i,:)=move1;
- sequence = seq1;
- else
- moves(i,:)=move2;
- sequence = seq2;
- end
- runningScore(i)= sum(abs(sequence-target));
+ if i==bb
+ s = abs(ss-tt);
+ [maxs,imaxs]=max(s); % worst current gene
+ [mins,imins1] = min(abs(tt-ss(imaxs))) ; % best place to put the worst gene
+ move1 = [1 imaxs imins1 false];
+ seq1 = doMoves(ss,move1,n);
+ s1= sum(abs(seq1-tt));
+ [mins,imins2] = min(abs(tt(imaxs)-ss)) ; % best gene to put in the worst place
+ move2 = [1 imins2 imaxs false];
+ seq2 = doMoves(ss,move2,n);
+ s2= sum(abs(seq2-tt));
+ if s1<s2
+ moves(i,:)=move1;
+ ss = seq1;
+ else
+ moves(i,:)=move2;
+ ss = seq2;
end
+ runningScore(i)= sum(abs(ss-tt));
+ end
- score= runningScore(budget) ;
+ score= runningScore(bb) ;
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function moves=oldconvolutedsolver(sequence,target,budget)
- seq1=sequence; score1 = 67108864;
+ function moves=oldconvolutedsolver(ss,tt,bb)
+ seq1=ss; score1 = 67108864;
- L = length(target);
+ L = length(tt);
if L < 97
- moves1 = markussolver(sequence, target, budget);
- score1 = findscore(sequence, target, budget, moves1);
- if score1<=50
- moves=moves1;
- return
- end
+ moves1 = markussolver(ss, tt, bb);
+ score1 = findscore(ss, tt, bb, moves1);
+ if score1<=50
+ moves=moves1;
+ return
end
+ end
- moves = new2solver(sequence, target, budget);
- score = findscore(sequence,target,budget,moves);
+ moves = new2solver(ss, tt, bb);
+ score = findscore(ss,tt,bb,moves);
if score<=50
- return
+ return
end
if score1<score,
- moves = moves1;
- score=score1;
+ moves = moves1;
+ score=score1;
end
seq3 = doMoves(seq1,moves,L);
- score3 = sum(abs(seq3-target));
- mybudget = budget-size(moves,1);
+ score3 = sum(abs(seq3-tt));
+ mybudget = bb-size(moves,1);
if mybudget > 0
- tempmoves=ukilaccuratesolver(seq3,target,mybudget);
- seq4=doMoves(seq3,tempmoves,L);
- score4=sum(abs(seq4-target));
- if score4<score3
- moves=[moves;tempmoves];
- end
+ tempmoves=ukilaccuratesolver(seq3,tt,mybudget);
+ seq4=doMoves(seq3,tempmoves,L);
+ score4=sum(abs(seq4-tt));
+ if score4<score3
+ moves=[moves;tempmoves];
end
+ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function moves=ukilaccuratesolver(s,t,budget)
- moves = zeros(budget,4);
+ function moves=ukilaccuratesolver(s,t,bb)
+ moves = zeros(bb,4);
L = length(s);
- for iter=1:budget
- diff=abs(s-t);
- [val,index]=sort(diff,'descend');
+ for iter=1:bb
+ diff=abs(s-t);
+ [val,index]=sort(diff,'descend');
- % first do simple
- value = inf(4,5);
- for i=1:2
- diff2 = abs(s(index(i))-t);
- [val2, index2] = sort(diff2);
- for j=1:3
- [localmoves, score] = removefromto(s,t,index(i),index2(j));
- if ~isempty(localmoves)
- value(i,:) = [localmoves score];
- end
- end
- end
- [mv,mi]=min(value(:,5));
+ % first do simple
+ value = inf(4,5);
+ for i=1:2
+ diff2 = abs(s(index(i))-t);
+ [val2, index2] = sort(diff2);
+ for j=1:3
+ [localmoves, score] = removefromto(s,t,index(i),index2(j));
+ if ~isempty(localmoves)
+ value(i,:) = [localmoves score];
+ end
+ end
+ end
+ [mv,mi]=min(value(:,5));
- if ~isempty(mi)
- localmoves=value(mi,1:4);
- moves(iter,:) = localmoves;
- s = doMove(s,localmoves(1), localmoves(2), localmoves(3), localmoves(4), L);
- else
- moves = moves(1:iter-1,:);
- return
- end
+ if ~isempty(mi)
+ localmoves=value(mi,1:4);
+ moves(iter,:) = localmoves;
+ s = doMove(s,localmoves(1), localmoves(2), localmoves(3), localmoves(4), L);
+ else
+ moves = moves(1:iter-1,:);
+ return
end
+ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function moves=anothersolver(seq2,MC2,budget)
- moves=zeros(budget,4);
+ function moves=anothersolver(seq2,MC2,bb)
+ moves=zeros(bb,4);
L = length(seq2);
- for i=1:budget
- diff=abs(seq2-MC2);
- [Dim,index]=max(diff);
- diff2=abs(seq2(index)-MC2);[N0T,J_m]=min(diff2);
- newmoves=removefromto(seq2,MC2,index,J_m);
- if ~isempty(newmoves)
- moves(i,:) = newmoves;
- seq2=doMoves(seq2,newmoves,L);
- end
+ for i=1:bb
+ diff=abs(seq2-MC2);
+ [Dim,index]=max(diff);
+ diff2=abs(seq2(index)-MC2);[N0T,J_m]=min(diff2);
+ newmoves=removefromto(seq2,MC2,index,J_m);
+ if ~isempty(newmoves)
+ moves(i,:) = newmoves;
+ seq2=doMoves(seq2,newmoves,L);
end
+ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [moves,INITIAL]=removefromto(seq2,tar2,index,index2)
L = length(seq2);
INITIAL=sum(abs(seq2-tar2));
tempmoves=removefromtosub1(seq2,index,index2);
if isempty(tempmoves)
- moves=ones(0,4);
+ moves=ones(0,4);
else
- n1 = size(tempmoves,1);
- sn1=zeros(n1,L);
- tscore1=zeros(n1,L);
- for i=1:n1
- sn1(i,:)=oldsplice(seq2,tempmoves(i,1),tempmoves(i,2),tempmoves(i,3),tempmoves(i,4));
- tscore1(i,:)=sn1(i,:)-tar2;
- end
- tscore1=sum(abs(tscore1),2);
- [minval,minindex1]=min(tscore1);
- moves=tempmoves(minindex1,:);
- INITIAL=minval;
+ n1 = size(tempmoves,1);
+ sn1=zeros(n1,L);
+ tscore1=zeros(n1,L);
+ for i=1:n1
+ sn1(i,:)=oldsplice(seq2,tempmoves(i,1),tempmoves(i,2),tempmoves(i,3),tempmoves(i,4));
+ tscore1(i,:)=sn1(i,:)-tar2;
end
+ tscore1=sum(abs(tscore1),2);
+ [minval,minindex1]=min(tscore1);
+ moves=tempmoves(minindex1,:);
+ INITIAL=minval;
+ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tempmoves=removefromtosub1(seq2,index,index2)
if index<index2
- N=abs(index2-index)+1;
- c = (1:N)';
- w = ones(N,1);
- tempmoves1 = [c index*w (index2+1)-c w];
+ N=abs(index2-index)+1;
+ c = (1:N)';
+ w = ones(N,1);
+ tempmoves1 = [c index*w (index2+1)-c w];
- N=min([index-1,numel(seq2)-index2])+1;
- c = (2:N)';
- w = ones(N-1,1);
- tempmoves2 = [c (index+1)-c index2*w w];
+ N=min([index-1,numel(seq2)-index2])+1;
+ c = (2:N)';
+ w = ones(N-1,1);
+ tempmoves2 = [c (index+1)-c index2*w w];
- N=index;
- c = (2:N)';
- z = zeros(N-1,1);
- tempmoves3 = [c (index+1)-c (index2+1)-c z];
+ N=index;
+ c = (2:N)';
+ z = zeros(N-1,1);
+ tempmoves3 = [c (index+1)-c (index2+1)-c z];
- N=numel(seq2)-index2+1;
- c = (2:N)';
- z = zeros(N-1,1);
- tempmoves=[tempmoves1;tempmoves2;tempmoves3;[c index+z index2+z z]];
+ N=numel(seq2)-index2+1;
+ c = (2:N)';
+ z = zeros(N-1,1);
+ tempmoves=[tempmoves1;tempmoves2;tempmoves3;[c index+z index2+z z]];
elseif index>=index2
- N=abs(index2-index)+1;
- c = (1:N)';
- w = ones(N,1);
- tempmoves1 = [c (index+1)-c index2*w w];
+ N=abs(index2-index)+1;
+ N=N+min([numel(seq2)-index,index2-1]);
+ c = (1:N)';
+ w = ones(N,1);
+ tempmoves1 = [c (index+1)-c index2*w w];
- N=min([numel(seq2)-index,index2-1])+1;
- c = (2:N)';
- w = ones(N-1,1);
- tempmoves2 = [c index*w (index2+1)-c w];
+ N=min([numel(seq2)-index,index2-1])+1;
+ c = (2:N)';
+ w = ones(N-1,1);
+ tempmoves2 = [c index*w (index2+1)-c w];
- N=numel(seq2)-index+1;
- c = (2:N)';
- z = zeros(N-1,1);
- tempmoves3 = [c index+z index2+z z];
+ N=numel(seq2)-index+1;
+ c = (2:N)';
+ z = zeros(N-1,1);
+ tempmoves3 = [c index+z index2+z z];
- N=index2;
- c = (2:N)';
- z = zeros(N-1,1);
- tempmoves=[tempmoves1;tempmoves2;tempmoves3;[c (index+1)-c (index2+1)-c z]];
+ N=index2;
+ c = (2:N)';
+ z = zeros(N-1,1);
+ tempmoves=[tempmoves1;tempmoves2;tempmoves3;[c (index+1)-c (index2+1)-c z]];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function moves = new2solver(sequence, target, budget)
+ function moves = new2solver(ss, tt, bb)
swap_counter=2;
- sequence_ini=sequence;
- L = length(sequence);cumSumShiftLeftMatrix = zeros(L, L+1);
- cumSumShiftRightMatrix = cumSumShiftLeftMatrix;cumSumShiftLeftRevMatrix = cumSumShiftLeftMatrix;cumSumShiftRightRevMatrix = cumSumShiftLeftMatrix;moves = zeros(budget,4);moveNr = 1;
- while moveNr <= budget
- for k = 1:L
- index = 2:L-k+2;
- target1 = target(index-1);
- target2 = target(k:L);
- cumSumShiftLeftMatrix(k,index) = sequence(k:L) - target1;
- cumSumShiftRightMatrix(k,index) = sequence(1:L-k+1) - target2;
- cumSumShiftLeftRevMatrix(k,index) = sequence(L-k+1:-1:1) - target1;
- cumSumShiftRightRevMatrix(k,index) = sequence(L:-1:k) - target2;
- end
- cumSumShiftLeftMatrix=cumsum(abs(cumSumShiftLeftMatrix),2);
- cumSumShiftRightMatrix=cumsum(abs(cumSumShiftRightMatrix),2);
- cumSumShiftLeftRevMatrix=cumsum(abs(cumSumShiftLeftRevMatrix),2);
- cumSumShiftRightRevMatrix=cumsum(abs(cumSumShiftRightRevMatrix),2);
- cumSumFromleft = cumSumShiftLeftMatrix(1,:);cumSumFromRight = cumSumFromleft(L+1) - cumSumFromleft;
- bestscore = cumSumFromleft(L+1);
- relMove = moveNr/budget;
- tmpscore = bestscore;
- curMoveLengthVector = [(33:4:L-1), 32, 29:-1:1];
- curMoveLengthVector(curMoveLengthVector >= L) = [];
- zeroIndex = (abs(sequence - target) < 1e-3)*floor(relMove);
- blen=1;bkFrom=1;bkTo=1;bFlip=false;
+ sequence_ini=ss;
+ L = length(ss);cumSumShiftLeftMatrix = zeros(L, L+1);
+ cumSumShiftRightMatrix = cumSumShiftLeftMatrix;cumSumShiftLeftRevMatrix = cumSumShiftLeftMatrix;cumSumShiftRightRevMatrix = cumSumShiftLeftMatrix;moves = zeros(bb,4);moveNr = 1;
+ while moveNr <= bb
+ for k = 1:L
+ index = 2:L-k+2;
+ target1 = tt(index-1);
+ target2 = tt(k:L);
+ cumSumShiftLeftMatrix(k,index) = ss(k:L) - target1;
+ cumSumShiftRightMatrix(k,index) = ss(1:L-k+1) - target2;
+ cumSumShiftLeftRevMatrix(k,index) = ss(L-k+1:-1:1) - target1;
+ cumSumShiftRightRevMatrix(k,index) = ss(L:-1:k) - target2;
+ end
+ cumSumShiftLeftMatrix=cumsum(abs(cumSumShiftLeftMatrix),2);
+ cumSumShiftRightMatrix=cumsum(abs(cumSumShiftRightMatrix),2);
+ cumSumShiftLeftRevMatrix=cumsum(abs(cumSumShiftLeftRevMatrix),2);
+ cumSumShiftRightRevMatrix=cumsum(abs(cumSumShiftRightRevMatrix),2);
+ cumSumFromleft = cumSumShiftLeftMatrix(1,:);cumSumFromRight = cumSumFromleft(L+1) - cumSumFromleft;
+ bestscore = cumSumFromleft(L+1);
+ relMove = moveNr/bb;
+ tmpscore = bestscore;
+ curMoveLengthVector = [(33:4:L-1), 32, 29:-1:1];
+ curMoveLengthVector(curMoveLengthVector >= L) = [];
+ zeroIndex = (abs(ss - tt) < 1e-3)*floor(relMove);
+ blen=1;bkFrom=1;bkTo=1;bFlip=false;
- for len = curMoveLengthVector
- cssLm = cumSumShiftLeftMatrix(len+1,:);
- cssRm = cumSumShiftRightMatrix(len+1,:);
- kMax = L - len + 1;
- for kFrom = 1:kMax
- for kTo = 1:kMax*~zeroIndex(kFrom)
+ for len = curMoveLengthVector
+ cssLm = cumSumShiftLeftMatrix(len+1,:);
+ cssRm = cumSumShiftRightMatrix(len+1,:);
+ kMax = L - len + 1;
+ for kFrom = 1:kMax
+ for kTo = 1:kMax*~zeroIndex(kFrom)
- if kFrom <= kTo
- score1 =cumSumFromleft(kFrom) + cssLm(kTo) +cumSumFromRight(kTo+len) - cssLm(kFrom);
- else
- score1 =cumSumFromleft(kTo) + cssRm(kFrom) + cumSumFromRight(kFrom+len) - cssRm(kTo);
- end
+ if kFrom <= kTo
+ score1 =cumSumFromleft(kFrom) + cssLm(kTo) +cumSumFromRight(kTo+len) - cssLm(kFrom);
+ else
+ score1 =cumSumFromleft(kTo) + cssRm(kFrom) + cumSumFromRight(kFrom+len) - cssRm(kTo);
+ end
- if len ~= 1
- seqRevIndex = L-kFrom-len+2;
- if kTo >= seqRevIndex
- curScore = score1 +cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex+len) - ...
- cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex);
- else
- curScore = score1 + ...
- cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo+len) -cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo);
- end
- if curScore < bestscore
- bestscore = curScore;
- blen=len;bkFrom=kFrom;bkTo=kTo; bFlip=true;
- end
- end
+ if len ~= 1
+ seqRevIndex = L-kFrom-len+2;
+ if kTo >= seqRevIndex
+ curScore = score1 +cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex+len) - ...
+ cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex);
+ else
+ curScore = score1 + ...
+ cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo+len) -cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo);
+ end
+ if curScore < bestscore
+ bestscore = curScore;
+ blen=len;bkFrom=kFrom;bkTo=kTo; bFlip=true;
+ end
+ end
- if len<=L/2
- if kFrom <= kTo
- curScore = score1 +cumSumShiftRightMatrix(kTo-kFrom+1, kFrom+len) -cumSumShiftRightMatrix(kTo-kFrom+1, kFrom);
- else
- curScore = score1 + cumSumShiftLeftMatrix(kFrom-kTo+1, kTo+len) -cumSumShiftLeftMatrix(kFrom-kTo+1, kTo);
- end
- if curScore < bestscore
- bestscore = curScore;blen=len;bkFrom=kFrom;bkTo=kTo;
- bFlip=false;
- end
- end
- end
- end
- end
- moves(moveNr,:)=[blen bkFrom bkTo bFlip];sequence = doMove(sequence, blen, bkFrom, bkTo, bFlip, L);
- swap_counter=swap_counter-1;
+ if len<=L/2
+ if kFrom <= kTo
+ curScore = score1 +cumSumShiftRightMatrix(kTo-kFrom+1, kFrom+len) -cumSumShiftRightMatrix(kTo-kFrom+1, kFrom);
+ else
+ curScore = score1 + cumSumShiftLeftMatrix(kFrom-kTo+1, kTo+len) -cumSumShiftLeftMatrix(kFrom-kTo+1, kTo);
+ end
+ if curScore < bestscore
+ bestscore = curScore;blen=len;bkFrom=kFrom;bkTo=kTo;
+ bFlip=false;
+ end
+ end
+ end
+ end
+ end
+ moves(moveNr,:)=[blen bkFrom bkTo bFlip];ss = doMove(ss, blen, bkFrom, bkTo, bFlip, L);
+ swap_counter=swap_counter-1;
- if swap_counter<=0
- % if moveNr>1
+ if swap_counter<=0
+ % if moveNr>1
- for jj = 1:2;
- rfac = (jj>1)*0.8;
- seq1 = splice(sequence_ini,moves,moveNr-2);
- [movesS,score]=swapper(seq1,target,2,rfac);
- if score < bestscore
- moves(moveNr-1:moveNr,:) = movesS;
- bestscore = score;
- sequence = splice(seq1,movesS,2);
- swap_counter=2;
- end
- end
- end
+ for jj = 1:2;
+ rfac = (jj>1)*0.8;
+ seq1 = splice(sequence_ini,moves,moveNr-2);
+ [movesS,score]=swapper(seq1,tt,2,rfac);
+ if score < bestscore
+ moves(moveNr-1:moveNr,:) = movesS;
+ bestscore = score;
+ ss = splice(seq1,movesS,2);
+ swap_counter=2;
+ end
+ end
+ end
- moveNr = moveNr + 1;
- if tmpscore==bestscore
- moves = moves(1:moveNr-1,:);
- return
- end
+ moveNr = moveNr + 1;
+ if tmpscore==bestscore
+ moves = moves(1:moveNr-1,:);
+ return
end
+ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function scorecheck = findscore(sequence,target,budget,moves)
- seq = splice(sequence,moves,budget);
- scorecheck = sum(abs(seq-target));
+ function scorecheck = findscore(ss,tt,bb,moves)
+ seq = splice(ss,moves,bb);
+ scorecheck = sum(abs(seq-tt));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function seq = splice(seq,moves,budget)
- if size(moves,1)>budget
- moves = moves(1:budget,:);
+ function seq = splice(seq,moves,bb)
+ if size(moves,1)>bb
+ moves = moves(1:bb,:);
end
Z65 = numel(seq);moves(:,1) = max(1,min(Z65,moves(:,1)));
moves(:,2) = max(1,min(Z65-moves(:,1)+1,moves(:,2)));
moves(:,3) = max(1,min(Z65-moves(:,1)+1,moves(:,3)));
moves(:,4) = moves(:,4) ~= 0;
for i =1:size(moves,1)
- le = moves(i,1);ai = moves(i,2);bi = moves(i,3);
- flip = moves(i,4);
- if flip
- if ai>bi
- seq = seq([1:bi-1 ai+le-1:-1:ai bi:ai-1 ai+le:end]);
- else
- seq = seq([1:ai-1 ai+le:bi+le-1 ai+le-1:-1:ai bi+le:end]);
- end
- else
- if ai>bi
- seq = seq([1:bi-1 ai:ai+le-1 bi:ai-1 ai+le:end]);
- else
- seq = seq([1:ai-1 ai+le:bi+le-1 ai:ai+le-1 bi+le:end]);
- end
- end
+ le = moves(i,1);ai = moves(i,2);bi = moves(i,3);
+ flip = moves(i,4);
+ if flip
+ if ai>bi
+ seq = seq([1:bi-1 ai+le-1:-1:ai bi:ai-1 ai+le:end]);
+ else
+ seq = seq([1:ai-1 ai+le:bi+le-1 ai+le-1:-1:ai bi+le:end]);
end
+ else
+ if ai>bi
+ seq = seq([1:bi-1 ai:ai+le-1 bi:ai-1 ai+le:end]);
+ else
+ seq = seq([1:ai-1 ai+le:bi+le-1 ai:ai+le-1 bi+le:end]);
+ end
+ end
+ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function moves = markussolver(sequence, target, budget)
- sequence_ini=sequence;
- moveLengthVector = [19 12 9 8 7];
+ function moves = markussolver(ss, tt, bb)
+ sequence_ini=ss;
+ moveLengthVector = [17 12 9 8 7];
moveLengthThresh = [0.9 0.6 0.3 0.1 -1];
- L = length(sequence);
+ L = length(ss);
index = moveLengthVector > L;
moveLengthVector(index) = [];
moveLengthThresh(index) = [];
cumSumShiftLeftMatrix = zeros(L, L+1);
cumSumShiftRightMatrix = cumSumShiftLeftMatrix;
cumSumShiftLeftRevMatrix = cumSumShiftLeftMatrix;
cumSumShiftRightRevMatrix = cumSumShiftLeftMatrix;
- moves = zeros(budget,4);
+ moves = zeros(bb,4);
moveNr = 1;
- while moveNr <= budget
- for k = 1:L
- index = 2:L-k+2;
- target1 = target(index-1);
- target2 = target(k:L);
- cumSumShiftLeftMatrix(k,index) = sequence(k:L) - target1;
- cumSumShiftRightMatrix(k,index) = sequence(1:L-k+1) - target2;
- cumSumShiftLeftRevMatrix(k,index) = sequence(L-k+1:-1:1) - target1;
- cumSumShiftRightRevMatrix(k,index) = sequence(L:-1:k) - target2;
- end
- cumSumShiftLeftMatrix=cumsum(abs(cumSumShiftLeftMatrix),2);
- cumSumShiftRightMatrix=cumsum(abs(cumSumShiftRightMatrix),2);
- cumSumShiftLeftRevMatrix=cumsum(abs(cumSumShiftLeftRevMatrix),2);
- cumSumShiftRightRevMatrix=cumsum(abs(cumSumShiftRightRevMatrix),2);
- cumSumFromleft = cumSumShiftLeftMatrix(1,:);
- cumSumFromRight = cumSumFromleft(L+1) - cumSumFromleft;
- bestscore = cumSumFromleft(L+1);
- tmpscore=bestscore;
- relMove = 1 - moveNr / budget;
- moveLengthIndex = find(relMove >= moveLengthThresh, 1, 'first');
- curLen = moveLengthVector(moveLengthIndex);
- curMoveLengthVector = [31:-2:21 curLen+4 curLen+3 curLen+2 curLen 6 5 4 3 2 1];
- curMoveLengthVector(curMoveLengthVector >= L) = [];
- curMoveLengthVector = fliplr(unique(curMoveLengthVector));
- zeroIndex = (abs(sequence - target) < 1e-3)*round(relMove-0.3);
- blen=1;bkFrom=1;bkTo=1;bFlip=false;
- for len = curMoveLengthVector
- cssLm = cumSumShiftLeftMatrix(len+1,:);
- cssRm = cumSumShiftRightMatrix(len+1,:);
- kMax = L - len + 1;
- for kFrom = 1:kMax
- for kTo = 1:(kFrom-1)*~zeroIndex(kFrom)
- score1 =cumSumFromleft(kTo) + cssRm(kFrom) + ...
- cumSumFromRight(kFrom+len) - cssRm(kTo);
- if len == 1
- curScore=score1+cumSumShiftLeftMatrix(kFrom-kTo+1, kTo+len) -cumSumShiftLeftMatrix(kFrom-kTo+1, kTo);
- else
- seqRevIndex = L-kFrom-len+2;
- if kTo >= seqRevIndex
- curScore = score1 +cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex+len) - ...
- cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex);
- else
- curScore = score1 + ...
- cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo+len) -cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo);
- end
- end
- if curScore < bestscore
- bestscore = curScore;
- blen=len;
- bkFrom=kFrom;
- bkTo=kTo;
- bFlip=true;
- end
- end
- for kTo = kFrom:kMax
- score1 =cumSumFromleft(kFrom) + cssLm(kTo) +cumSumFromRight(kTo+len) - cssLm(kFrom);
- if len == 1
- curScore = score1 +cumSumShiftRightMatrix(kTo-kFrom+1, kFrom+len) -cumSumShiftRightMatrix(kTo-kFrom+1, kFrom);
- else
- seqRevIndex = L-kFrom-len+2;
- if kTo >= seqRevIndex
- curScore = score1 +cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex+len) - ...
- cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex);
- else
- curScore = score1 + ...
- cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo+len) -cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo);
- end
- end
- if curScore < bestscore
- bestscore = curScore;blen=len;bkFrom=kFrom;bkTo=kTo;
- bFlip=true;
- end
- end
- end
- end
- moves(moveNr,:)=[blen bkFrom bkTo bFlip];
- sequence = doMove(sequence, blen, bkFrom, bkTo, bFlip, L);
+ while moveNr <= bb
+ for k = 1:L
+ index = 2:L-k+2;
+ target1 = tt(index-1);
+ target2 = tt(k:L);
+ cumSumShiftLeftMatrix(k,index) = ss(k:L) - target1;
+ cumSumShiftRightMatrix(k,index) = ss(1:L-k+1) - target2;
+ cumSumShiftLeftRevMatrix(k,index) = ss(L-k+1:-1:1) - target1;
+ cumSumShiftRightRevMatrix(k,index) = ss(L:-1:k) - target2;
+ end
+ cumSumShiftLeftMatrix=cumsum(abs(cumSumShiftLeftMatrix),2);
+ cumSumShiftRightMatrix=cumsum(abs(cumSumShiftRightMatrix),2);
+ cumSumShiftLeftRevMatrix=cumsum(abs(cumSumShiftLeftRevMatrix),2);
+ cumSumShiftRightRevMatrix=cumsum(abs(cumSumShiftRightRevMatrix),2);
+ cumSumFromleft = cumSumShiftLeftMatrix(1,:);
+ cumSumFromRight = cumSumFromleft(L+1) - cumSumFromleft;
+ bestscore = cumSumFromleft(L+1);
+ tmpscore=bestscore;
+ relMove = 1 - moveNr / bb;
+ moveLengthIndex = find(relMove >= moveLengthThresh, 1, 'first');
+ curLen = moveLengthVector(moveLengthIndex);
+ curMoveLengthVector = [31:-2:21 curLen+4 curLen+3 curLen+2 curLen 6 5 4 3 2 1];
+ curMoveLengthVector(curMoveLengthVector >= L) = [];
+ curMoveLengthVector = fliplr(unique(curMoveLengthVector));
+ zeroIndex = (abs(ss - tt) < 1e-3)*round(relMove-0.3);
+ blen=1;bkFrom=1;bkTo=1;bFlip=false;
+ for len = curMoveLengthVector
+ cssLm = cumSumShiftLeftMatrix(len+1,:);
+ cssRm = cumSumShiftRightMatrix(len+1,:);
+ kMax = L - len + 1;
+ for kFrom = 1:kMax
+ for kTo = 1:(kFrom-1)*~zeroIndex(kFrom)
+ score1 =cumSumFromleft(kTo) + cssRm(kFrom) + ...
+ cumSumFromRight(kFrom+len) - cssRm(kTo);
+ if len == 1
+ curScore=score1+cumSumShiftLeftMatrix(kFrom-kTo+1, kTo+len) -cumSumShiftLeftMatrix(kFrom-kTo+1, kTo);
+ else
+ seqRevIndex = L-kFrom-len+2;
+ if kTo >= seqRevIndex
+ curScore = score1 +cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex+len) - ...
+ cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex);
+ else
+ curScore = score1 + ...
+ cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo+len) -cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo);
+ end
+ end
+ if curScore < bestscore
+ bestscore = curScore;
+ blen=len;
+ bkFrom=kFrom;
+ bkTo=kTo;
+ bFlip=true;
+ end
+ end
+ for kTo = kFrom:kMax
+ score1 =cumSumFromleft(kFrom) + cssLm(kTo) +cumSumFromRight(kTo+len) - cssLm(kFrom);
+ if len == 1
+ curScore = score1 +cumSumShiftRightMatrix(kTo-kFrom+1, kFrom+len) -cumSumShiftRightMatrix(kTo-kFrom+1, kFrom);
+ else
+ seqRevIndex = L-kFrom-len+2;
+ if kTo >= seqRevIndex
+ curScore = score1 +cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex+len) - ...
+ cumSumShiftRightRevMatrix(kTo - seqRevIndex + 1, seqRevIndex);
+ else
+ curScore = score1 + ...
+ cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo+len) -cumSumShiftLeftRevMatrix(seqRevIndex - kTo + 1, kTo);
+ end
+ end
+ if curScore < bestscore
+ bestscore = curScore;blen=len;bkFrom=kFrom;bkTo=kTo;
+ bFlip=true;
+ end
+ end
+ end
+ end
+ moves(moveNr,:)=[blen bkFrom bkTo bFlip];
+ ss = doMove(ss, blen, bkFrom, bkTo, bFlip, L);
- % if swap_counter<=0
- if moveNr>1
- for jj = 1:2
- rfac = (jj>1)*0.15;
- seq1 = splice(sequence_ini,moves,moveNr-2);
- [movesS,score]=swapper(seq1,target,2,rfac);
- if score < bestscore
- moves(moveNr-1:moveNr,:) = movesS;
- bestscore = score;
- sequence = splice(seq1,movesS,2);
- end
- end
- end
+ % if swap_counter<=0
+ if moveNr>1
+ for jj = 1:2
+ rfac = (jj>1)*0.15;
+ seq1 = splice(sequence_ini,moves,moveNr-2);
+ [movesS,score]=swapper(seq1,tt,2,rfac);
+ if score < bestscore
+ moves(moveNr-1:moveNr,:) = movesS;
+ bestscore = score;
+ ss = splice(seq1,movesS,2);
+ end
+ end
+ end
- moveNr = moveNr + 1;
- if tmpscore==bestscore
- moves = moves(1:moveNr-1,:);
- return
- end
+ moveNr = moveNr + 1;
+ if tmpscore==bestscore
+ moves = moves(1:moveNr-1,:);
+ return
end
+ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function seq = oldsplice(seq, le, ai, bi, flip)
if flip
- if ai>bi
- seq = seq([1:bi-1 ai+le-1:-1:ai bi:ai-1 ai+le:end]);
- else
- seq = seq([1:ai-1 ai+le:bi+le-1 ai+le-1:-1:ai bi+le:end]);
- end
+ if ai>bi
+ seq = seq([1:bi-1 ai+le-1:-1:ai bi:ai-1 ai+le:end]);
else
- if ai>bi
- seq = seq([1:bi-1 ai:ai+le-1 bi:ai-1 ai+le:end]);
- else
- seq = seq([1:ai-1 ai+le:bi+le-1 ai:ai+le-1 bi+le:end]);
- end
+ seq = seq([1:ai-1 ai+le:bi+le-1 ai+le-1:-1:ai bi+le:end]);
end
+ else
+ if ai>bi
+ seq = seq([1:bi-1 ai:ai+le-1 bi:ai-1 ai+le:end]);
+ else
+ seq = seq([1:ai-1 ai+le:bi+le-1 ai:ai+le-1 bi+le:end]);
+ end
+ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function seq = doMove(seq, le, ai, bi, flip, L)
if flip
- if ai>bi
- seq = seq([1:bi-1 ai+le-1:-1:ai bi:ai-1 ai+le:L]);
- else
- seq = seq([1:ai-1 ai+le:bi+le-1 ai+le-1:-1:ai bi+le:L]);
- end
+ if ai>bi
+ seq = seq([1:bi-1 ai+le-1:-1:ai bi:ai-1 ai+le:L]);
else
- if ai>bi
- seq = seq([1:bi-1 ai:ai+le-1 bi:ai-1 ai+le:L]);
- else
- seq = seq([1:ai-1 ai+le:bi+le-1 ai:ai+le-1 bi+le:L]);
- end
+ seq = seq([1:ai-1 ai+le:bi+le-1 ai+le-1:-1:ai bi+le:L]);
end
+ else
+ if ai>bi
+ seq = seq([1:bi-1 ai:ai+le-1 bi:ai-1 ai+le:L]);
+ else
+ seq = seq([1:ai-1 ai+le:bi+le-1 ai:ai+le-1 bi+le:L]);
+ end
+ end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function seq = doMoves(seq,moves,L)
for i =1:size(moves,1)
- le = moves(i,1);
- ai = moves(i,2);
- bi = moves(i,3);
- if moves(i,4)
- if ai>bi
- seq = seq([1:bi-1 ai+le-1:-1:ai bi:ai-1 ai+le:L]);
- else
- seq = seq([1:ai-1 ai+le:bi+le-1 ai+le-1:-1:ai bi+le:L]);
- end
- else
- if ai>bi
- seq = seq([1:bi-1 ai:ai+le-1 bi:ai-1 ai+le:L]);
- else
- seq = seq([1:ai-1 ai+le:bi+le-1 ai:ai+le-1 bi+le:L]);
- end
- end
+ le = moves(i,1);
+ ai = moves(i,2);
+ bi = moves(i,3);
+ if moves(i,4)
+ if ai>bi
+ seq = seq([1:bi-1 ai+le-1:-1:ai bi:ai-1 ai+le:L]);
+ else
+ seq = seq([1:ai-1 ai+le:bi+le-1 ai+le-1:-1:ai bi+le:L]);
end
-
+ else
-
+ if ai>bi
+ seq = seq([1:bi-1 ai:ai+le-1 bi:ai-1 ai+le:L]);
+ else
+ seq = seq([1:ai-1 ai+le:bi+le-1 ai:ai+le-1 bi+le:L]);
+ end
+ end
+ end
|