# Diffing "Santa Claus" and "MM03"

 Title: Santa Claus MM03 Author: MR Keenan MadMark Submitted: 2002-11-13 16:50:34 UTC 2002-11-13 16:53:44 UTC Status: Passed Passed Score: 2160.44 2160.61 Result: 2156.474 2156.474 CPU Time: 66.456 67.315 Code: ```function c = solver(a) % SOLVER: Takes the amino acid chain, a, as input and returns the % corresponding configuration sequence, c. The configuration sequence, c, % consists of the directions to the next successive amino acid: % Direction: 1 = east, i = north, -1 = west, -i = south %%%%%%%%%%%%%% Early termination cases %%%%%%%%%%%%%%%%%%%%%%%% M=length(a); hcl=M-1; c=ones(hcl,1); q=find(a); s=length(q); if s<4 if s==2 brk=floor((q(1)+q(2))/2); c(brk)=-1i; c(brk+1:hcl)=-1; elseif s==3 if q(2)>q(1)+1 brk=floor((q(1)+q(2))/2); c(q(1):brk-1)=1i; c(brk+1:q(2)-1)=-1i; end beven=rem(q(2)-q(1),2); if (q(3)==q(2)+1) c(q(2))=-1i; else if beven brk=floor((q(2)+q(3)-1)/2); else brk=floor((q(2)+q(3))/2); end c(brk)=-1i; c(brk+1:q(3)-1)=-1; c(q(3):hcl)=-1i; end end return end if M > 90 & s >= hcl % sliploop for M=100, s = 99 or 100. if (M >= 99) loop = [1; -1i; 1; 1i; 1; -1i; -1i; 1; 1i; 1i; 1; -1i; -1i; 1; 1i; 1i; 1; -1i; 1; 1i; 1i; 1; 1i; -1; 1i; 1; 1i; -1; 1i; -1; 1i; 1; 1i; -1; 1i; -1; -1; -1; -1i; -1i; 1; 1i; 1; -1i; -1i; -1; -1; -1; 1i; 1i; 1i; -1; -1i; -1; -1; -1i; 1; 1; -1i; -1; -1; -1; -1i; 1; 1; 1; 1; 1; 1; 1; 1; -1i; -1; -1; -1; -1; -1; -1; -1; -1; -1i; 1; 1; 1; 1; 1; 1; 1; 1; -1i; -1; -1; -1; -1; -1; -1; -1; -1; -1i; 1]; if s < length(loop) f = find(a==0); f = M-f(1); c = [loop(f+1:end); loop(1:f)]; else c = loop; end else % M=97, s=97 c = [1; 1i; -1; 1i; -1; -1; -1; -1i; -1; -1; -1i; -1i; -1; -1; -1i; 1; -1i; -1; -1i; 1; -1i; -1i; 1; -1i; 1; -1i; 1; 1; 1; 1; 1; 1i; 1; 1i; -1; -1; -1i; -1; 1i; -1; -1i; -1; 1i; -1; 1i; -1; 1i; 1; 1; -1i; 1; 1; 1; 1; 1; 1; 1i; -1; -1; -1; -1; -1; 1i; 1; 1; 1; 1; 1; 1i; -1; -1; -1; -1; -1; -1; -1i; -1; -1; 1i; 1; 1i; 1i; 1; -1i; 1; 1i; 1i; 1; -1i; -1i; 1; 1; 1; 1; 1i; -1; 1i]; end c = c(1:hcl); return end % Permutator, written by Yuval Cohen aCropL=min(M-q(1)+1,q(s)); if aCropL<13 dl=M-aCropL; if dl head = q(1)-1 >= M-q(s); if head at=a(q(1):M); else at=a(1:q(s)); end else at=a; end ctl=aCropL-1; m=1; permN=3^(ctl-1); t=ones(ctl,1); for n=2:ctl; bendLeft=[ones(n-1,1); 1i+zeros(ctl-n+1,1)]; bendRight=[ones(n-1,1); -1i+zeros(ctl-n+1,1)]; t=[t t.*bendLeft(:,ones(1,m)), t.*bendRight(:,ones(1,m))]; goods=logical(ones(1,3*m)); if n>3 p=sum(t(n-3:n,:),1); goods(find(p==0))=0; end nn=5; while nn=10 CCk=CCk+1;CC(:,CCk)= magicball(aori); end setspecialcases %Various authors kk=size(SCM,2); CCk=CCk+kk;CC(:,CCk-kk+1:CCk)= SCM(1:M-1,:); CCk=CCk+1;CC(:,CCk)= SCM(M-1:-1:1,[4]); %if M < 80 %Other guys ... :) (various) CCk=CCk+1;CC(:,CCk)= solverpira(aori); %else %if M > 10 %searchnfold by Rog CCk = CCk+1; CC(:,CCk) = searchnfold(aori); %end CCk=CCk+1;CC(:,CCk)= calccsquare(a); % new line from Mr. Bond CC = [CC(:,1:CCk) CC(hcl:-1:1,1:CCk)]; c = minenergy(CC(:,1:CCk)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %++++++++++++++++++++++++++++++++++++++++++++++++++++++ %++++++++++++++++++++++++++++++++++++++++++++++++++++++ %++++++++++++++++++++++++++++++++++++++++++++++++++++++ function c=fillrows(a,sw) % by Lucio Andrade a=a(:);M=length(a); v=diff(find(diff([0;~a;0]))); vz=v(1:2:end); lvz=length(vz); v=diff(find(diff([0;a;0]))); vo=v(1:2:end); lvo=length(vo); mu=ones(lvo,1);vodmu=ceil(vo./mu); while max(vodmu)>sum(mu)+sw can=find(vodmu==max(vodmu)); [dum0,h]=max(abs(can-(lvo+1)/2)); mu(can(h))=mu(can(h))+1; vodmu=ceil(vo./mu); end nro=sum(mu); voe=zeros(nro,1);vze=zeros(nro,1); vodmu=floor(vo./mu);j=1;k=1; while j<=nro tem=zeros(mu(k),1)+vodmu(k); tem(1:rem(vo(k),mu(k)))=vodmu(k)+1; if (j+j+mu(k)-1)2&xI12&xI22&yI12&yI20); s = ceil(sqrt(r0+sum(a))); f = find(a); k = 10; e = Inf; while k > 0 k = k-1; sign = 1; f1 = f; if (f1(1) > 1) c = ones(1,f1(1)-2); else c = []; end t0 = 0; while ~isempty(f1) t = length(c)+1+t0+s; f1 = f1(f1 > t); if ~isempty(f1) & f1(1) > t+1 t1 = floor((f1(1)-t)/2+randn/1.5); else t1 = 0; end c = [c ones(1,1+t0+s+t1)]; t0 = t1; %sign = -sign; end if (length(c) < al-1) c = [c ones(1,al-1-length(c))]; else c = c(1:(al-1)); end end c = ones(al-1,1); n = cumsum([3 10:8:al-1]); temp = find(n>al-1); if ~isempty(temp) n = n(1:temp(1)-1); end add = 0; for kk = 1:length(n) c(n(kk):n(kk)+kk+add) = -1; c(n(kk)+kk+add+1:n(kk)+2*kk+2*add+1) = 1i; add = add + 1; end n = cumsum([2 8:8:al-1]); temp = find(n>al-1); if ~isempty(temp) n = n(1:temp(1)-1); end add = -1; for kk = 1:length(n) c(n(kk):n(kk)+kk+add) = -1i; add = add + 1; end while sum(diff([0;a])>0) > sqrt(sum(a)); a = (a + [0; a(1:end-1)] + [a(2:al); 0]) > 0; end; d = diff([0;a;0]); i1 = find(d>0)-1; N = find(d<0) - i1 - 1; C = zeros(al-1,11); Ck=1; for k = sqrt(sum(a)) + (-1:.2:1); m = []; for i=1:length(i1); n = max(0, ceil ((N(i)-2.5*k) / k)); r = N(i) - n*k; if r>1.5*k; m = [m, i1(i)+[r/4, r/2-k/2+k*(1:n), 3*r/4+k*n]]; else m = [m, i1(i) + r/2]; end; end; cb = zeros(length(a)-1,1); for i=1:length(m)-1; d = (m(i)+m(i+1))/2; b = round(d); cb(b) = 1; m(i+1) = m(i+1)-2*(d-b); end; c2 = 1-2*rem(cumsum(cb(:)),2); c2(logical(cb))=1i; C(:,Ck)=c2; Ck=Ck+1; end [c2,result]=minenergy(C(:,1:Ck-1)); if result3) secondmoment=myvar(real(p(iones,:)))+myvar(imag(p(iones,:))); [scratch,ixlist] = sort(secondmoment); ixlist=ixlist(1:3); if ~any(ixlist==1) ixlist=[1 ixlist(1:3)]; end else ixlist=1:size(c,2); end [e,x]=min(sum(abs(p(I,ixlist)-p(J,ixlist)))); c=c(:,ixlist(x)); if ixlist(x)==1 valt=0; else valt=1; end function y = myvar(x) m=size(x,1); avg = sum(x)/m; centerx = x - ones(m,1)*avg; y = sum((conj(centerx) .* centerx))/m; function y = uvvar(x) m=size(x,1); avg = sum(x)/m; centerx = x - avg(ones(m,1),:); y = sum((conj(centerx) .* centerx))/m; function y = uvvar(x) m=size(x,1); avg = sum(x)/m; centerx = x - avg(ones(m,1),:); y = sum(centerx.*centerx); function ch = halfchord(r,h) ch = sqrt(r*r-h*h); function c = snakesquare(a) c1 = sum(a); r0 = sum(diff(a)>0); s = ceil(sqrt(r0+c1)); f1 = find(a); al = length(a); sign = 1; if (f1(1) > 1) c = ones(f1(1)-2,1); else c = []; end t0 = 0; while ~isempty(f1) t = length(c)+1+t0+s; f1 = f1(f1 > t); if ~isempty(f1) & f1(1) > t+1 t1 = floor((f1(1)-t)/2+randn/1.5); else t1 = 0; end c = [c;sign(ones(t0+s+t1,1));1i]; t0 = t1; sign = -sign; end if (length(c) < al-1) c = [c;-sign(ones(al-1-length(c),1))]; else c = c(1:(al-1)); end %%% RLE Snake. function c = rlesnake(a,N1) % run length encode a N = length(a); S = round(sqrt(N1)*1.5); S1=S/2; %pminrev=ceil(S/3); % better estimation?(depending on length(edges),N1,....) pminrev=1; edges = [0;find(diff(a));N]; n=1; while nS l=ceil(k/2)-1; c(l+1)=1i; c(l+2:k)=-1; dd=-1; p=rem(k,2)-l/2; % position (around middle) else p=(k-1)/2; dd = 1; end for nn = 2:2:length(edges) n=edges(nn); N=edges(nn+1); Nlong=0; if N>S Nlong=floor(N/2); % should be calculated better (depending on p and n) N=N-Nlong; end pd=N/2; % desired position for next row ap=abs(p); dpd=pd-ap; % offset to current position adpd=abs(dpd); nN=n+N; if dd*p<0 % next step would be farther from middle if n<=ap c(k:k+nN-1)=dd; p=p+dd*(nN); k=k+nN; N=0; else n1=ceil(ap); c(k:k+n1-1)=dd; k=k+n1; p=p+dd*n1; n=n-n1; nN=n+N; % recalculate position values ap=abs(p); dpd=pd-ap; % offset to current position adpd=abs(dpd); end end if N % not yet processed above if adpd>n % higher desired step then possible if ap>pd % step towards middle if ap0 p=p+dd*(n1-rem(n-n1,2)-N+1); c(k:k+n1+n2-1)=dd; dd=-dd; c(k+n1+n2)=1i; c(k+n1+n2+1:k+nN-1)=dd; else p=p-dd*(rem(n-n1,2)-1+n1+N); c(k:k+n2-1)=dd; dd=-dd; c(k+n2)=1i; c(k+n2+1:k+nN-1)=dd; end end % if else k=k+nN; if Nlong c(k)=1i; dd=-dd; c(k+1:k+Nlong-1)=dd; k=k+Nlong; p=p+dd*(Nlong-1); end end end function c=zigzag3(a,l) % zigzag3 al=length(a); k=1; % index in c c=zeros(al-1,1); d=1; % direction while 1 if k==1 % first line % this is put in this loop because the end of the loop has to be used % this hopefully works, but nice code is something completely different.... % sorry about that m=find(a==0); pos=0; if ~isempty(m)&m(1)<=l m=m(1)-1; n=find(a(m+2:al)); if ~isempty(n)&n>1 c(1:m-1)=d; n=floor(n(1)/2); c(m:m+n-1)=-1i; c(m+n)=d; c(m+n+1:m+2*n)=1i; pos=m; k=m+1+2*n; if a(k)==0|pos=l break; end n=n(1)-1; c(k:k+n-1)=d; pos=pos+n; k=k+n; n=find(a(k+1:al)); if isempty(n) break; elseif n(1)==1 c(k)=d; pos=pos+1; k=k+1; else n=floor(n(1)/2); c(k:k+n-1)=-1i; c(k+n)=1; c(k+n+1:k+2*n)=1i; k=k+2*n+1; pos=pos+1; if pos=al-1 break; end l1=find(a(k+1:al)); % can't be empty l1=l1(1); if l1>=3 l1=floor((l1-1)/2); c(k:k+l1-1)=d; k=k+l1; c(k)=1i; c(k+1:k+l1)=-d; k=k+l1+1; else c(k)=1i; k=k+1; end d=-d; end c=c(1:al-1); function c = rlesnakeold(a) N = length(a); % run length encode a edges = [1 1+find(a(2:N) ~= a(1:N-1))']; val = a(edges)'; run = [ edges(2:end)-edges(1:end-1) N-edges(end)+1 ]; N1 = sum(run(val>0)); S = round(sqrt(N1)); c = []; dd = 1; for nn = 1:length(val) V = val(nn); R = run(nn); if V if R < 1.5*S c = [c;dd(ones(1,R),1)]; else if R == 1 c = [c;dd]; elseif R == 2 c = [c;dd;1i]; dd = -dd; else k = round((R-1)/2); c = [c;dd(ones(1,k),1);1i;-dd(ones(1,R-k-1),1)]; dd = -dd; end end else if R == 1 c = [c;1i]; dd = -dd; elseif R == 2 c = [c;dd;1i]; dd = -dd; else k = round((R-1)/2); c = [c;dd(ones(1,k),1);1i;-dd(ones(1,R-k-1),1)]; dd = -dd; end end end if length(c) > N-1; c=c(1:N-1); end; function c=improvezigzagend(c,a,l,godown) i=find(a==0); if isempty(i) c=[]; % no trial return end al=length(a); cl=al-1; i=i(end); % last 0 j=find(a(1:i)); j=j(end); % last 1 before last 0 p=sum(c(1:j)); c1=c(j); n=i-j; if n<2|imag(p)<1 c=[]; return end n=floor(n/2); if c1==1i % up c(j:j+n-1)=c(j-1); c(j+n)=1i; c(j+n+1:cl)=-c(j-1); else c(j:j+n-1)=1i; c(j+n)=c1; c(j+n+1:j+2*n)=-1i; c(j+2*n+1:cl)=c1; end wentdown=0; p=sum(c); pr=real(p); if pr<-1&imag(p)+pr>0 j=al+pr+1; wentdown=1; elseif pr>l+2 j=al-pr+l+1; wentdown=1; end if wentdown c0=c; c(j:al-1)=-1i; if ~godown p=cumsum(c); if any(ismember(p(j:al-1),p(1:j-2))) c=c0; end end end %++++++++++++++++++++++++++++++++++++++++++++++++++++++ %++++++++++++++++++++++++++++++++++++++++++++++++++++++ %++++++++++++++++++++++++++++++++++++++++++++++++++++++ %++++++++++++++++++++++++++++++++++++++++++++++++++++++ function c=magicball(a) %by Lucio Andrade M=length(a); ao=find(a); aol=ao(end)-ao(1)+1; bll=[0 2 2 2 2 2 3 3 3 4 4 4 4 4 4 ... 4 5 5 5 6 5 5 6 6 5 6 6 6 6 6 ... 7 7 7 6 6 6 8 8 7 7 7 7 7 8 8 ... 8 8 8 8 8 8 8 8 8 9 8 8 8 9 9 ... 9 9 9 9 9 9 9 10 10 10 9 9 10 10 10 ... 10 10 10 11 11 11 10 10 10 10 11 10 11 11 11 ... 11 11 11 11 11 11 11 11 12 12]; BL=[0 0 0 0 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 0 0 3 3 0 0 0 0 0 0 0 0 0 0 3 3 0 0 0 0 0 0 0 0 0 0 3 3 3 0 0 0 0 0 0 0 0 0 3 3 3 0 0 0 0 0 0 0 0 0 3 3 3 0 0 0 0 0 0 0 0 0 2 3 3 2 0 0 0 0 0 0 0 0 3 3 3 3 0 0 0 0 0 0 0 0 3 3 3 3 0 0 0 0 0 0 0 0 3 4 4 3 0 0 0 0 0 0 0 0 3 4 4 3 0 0 0 0 0 0 0 0 4 4 4 4 0 0 0 0 0 0 0 0 4 4 4 4 0 0 0 0 0 0 0 0 3 4 5 4 3 0 0 0 0 0 0 0 3 4 5 4 3 0 0 0 0 0 0 0 3 4 5 4 3 0 0 0 0 0 0 0 2 3 5 5 3 2 0 0 0 0 0 0 4 5 5 4 3 0 0 0 0 0 0 0 4 4 5 5 4 0 0 0 0 0 0 0 3 4 5 5 4 3 0 0 0 0 0 0 3 4 5 5 4 3 0 0 0 0 0 0 5 5 5 5 5 0 0 0 0 0 0 0 4 5 6 6 5 4 0 0 0 0 0 0 4 5 6 6 5 4 0 0 0 0 0 0 4 5 6 6 5 4 0 0 0 0 0 0 4 5 6 6 5 4 0 0 0 0 0 0 4 5 6 6 5 4 0 0 0 0 0 0 3 4 6 7 6 4 3 0 0 0 0 0 3 4 6 7 6 4 3 0 0 0 0 0 3 4 6 7 6 4 3 0 0 0 0 0 5 6 7 7 5 4 0 0 0 0 0 0 5 6 7 7 6 5 0 0 0 0 0 0 5 6 7 7 6 5 0 0 0 0 0 0 2 4 6 7 7 6 4 2 0 0 0 0 2 4 6 7 7 6 4 2 0 0 0 0 5 6 7 7 7 6 5 0 0 0 0 0 5 6 7 7 7 6 5 0 0 0 0 0 5 6 7 7 7 6 5 0 0 0 0 0 5 6 7 7 7 6 5 0 0 0 0 0 5 6 7 7 7 6 5 0 0 0 0 0 4 5 7 8 8 7 5 4 0 0 0 0 4 5 7 8 8 7 5 4 0 0 0 0 4 5 7 8 8 7 5 4 0 0 0 0 4 5 7 8 8 7 5 4 0 0 0 0 4 5 7 8 8 7 5 4 0 0 0 0 5 6 7 8 8 7 5 3 0 0 0 0 5 6 7 8 8 7 5 4 0 0 0 0 5 6 7 8 8 7 6 5 0 0 0 0 5 6 7 8 8 7 6 5 0 0 0 0 4 5 7 8 8 8 7 6 0 0 0 0 5 6 7 8 9 8 6 5 0 0 0 0 3 5 7 8 9 8 7 5 3 0 0 0 6 7 8 8 8 8 7 6 0 0 0 0 6 7 8 8 8 8 7 6 0 0 0 0 6 7 8 8 8 8 7 6 0 0 0 0 5 6 8 9 9 8 7 5 3 0 0 0 5 6 8 9 9 8 7 5 3 0 0 0 5 6 8 9 9 9 8 6 5 0 0 0 5 6 8 9 9 9 8 6 5 0 0 0 5 6 8 9 9 9 8 6 5 0 0 0 5 6 8 9 9 9 8 6 5 0 0 0 5 6 8 9 9 9 8 6 5 0 0 0 5 7 8 9 10 9 8 6 5 0 0 0 5 7 8 9 10 9 8 6 5 0 0 0 4 5 7 9 10 10 9 7 5 4 0 0 4 5 7 9 10 10 9 7 5 4 0 0 4 5 7 9 10 10 9 7 5 4 0 0 6 7 9 10 10 9 8 7 6 0 0 0 6 7 9 10 10 9 8 7 6 0 0 0 4 6 8 9 10 10 9 8 6 4 0 0 4 6 8 9 10 10 9 8 6 4 0 0 5 6 8 9 10 10 9 8 6 5 0 0 5 6 8 9 10 10 9 8 6 5 0 0 5 7 9 10 10 10 10 8 5 4 0 0 5 7 9 10 10 10 10 8 5 4 0 0 3 5 8 9 10 11 10 9 8 5 3 0 3 5 8 9 10 11 10 9 8 5 3 0 3 5 8 9 10 11 10 9 8 5 3 0 6 7 9 10 10 10 10 9 7 6 0 0 6 7 9 10 10 10 10 9 7 6 0 0 6 7 9 10 10 10 10 9 7 6 0 0 6 7 9 10 11 11 10 9 7 5 0 0 5 6 8 10 11 11 10 9 8 5 3 0 5 7 9 10 11 11 10 9 8 7 0 0 5 6 8 10 11 11 11 10 8 6 5 0 5 6 8 10 11 11 11 10 8 6 5 0 5 6 8 10 11 11 11 10 8 6 5 0 5 6 8 10 11 11 11 10 8 6 5 0 3 5 8 10 11 11 11 11 9 7 6 0 7 8 9 10 11 11 11 10 8 5 3 0 3 5 8 10 11 11 11 11 10 8 6 0 6 7 9 10 11 11 10 10 9 7 6 0 6 7 9 10 11 11 10 10 9 7 6 0 6 7 9 10 11 12 11 10 9 7 6 0 6 7 9 10 11 12 11 10 9 7 6 0 4 6 9 10 11 12 12 11 10 9 6 4 4 6 9 10 11 12 12 11 10 9 6 4]; bl=cumsum(BL(aol,1:bll(aol))); numben=bll(aol)-1; bl(end+1:end+2)=bl(end)+5; ax=[a(ao(1):ao(end));zeros(6*bll(aol),1)]; P=[ 0 0 0; 0 0 -1; 0 0 1;-1 -1 0; 0 -1 -1; 0 1 1; 1 1 0;-1 -1 -1;-1 -1 1; 0 -1 -2; 0 1 2; 1 1 -1; 1 1 1;-1 -2 -1; 0 -1 -3; 0 1 3; 1 2 1;-1 -2 -2; 1 2 2;-1 -2 -3; 1 2 3;-1 -3 -3; 1 3 3;-1 -3 -4; 1 3 4;-1 -3 -5; 1 3 5]; hbl=[0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2; 0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2; 0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2;0 1 2]; if aol<6 P(end,:)=[];hbl(end,:)=[]; end if aol<5 P(end-1,:)=[];hbl(end-1,:)=[]; end if aol<4 P(end-2,:)=[];hbl(end-2,:)=[]; end j=1; if (bl(1)==1 | bl(2)==3) pr=find(P(:,1)>=0); [du,h]=min(sum(ax(bl(hbl(pr,:)+j)+P(pr,:))+ax(bl(hbl(pr,:)+j)+P(pr,:)+1),2)); else [du,h]=min(sum(ax(bl(hbl+j)+P)+ax(bl(hbl+j)+P+1),2)); end bl(j)=bl(j)+P(h,1); bl(j+1:end)=bl(j+1:end)+P(h,1)+P(h,1); for j=2:numben [du,h]=min(sum(ax(bl(hbl+j)+P)+ax(bl(hbl+j)+P+1),2)); bl(j)=bl(j)+P(h,1); bl(j+1:end)=bl(j+1:end)+P(h,1)+P(h,1); end c=ones(M-1,1); bl=bl+ao(1)-1; bl(bl>=ao(end))=[]; c(bl)=1i; sign=0; for i=1:M-1 if c(i)==1i sign=~sign; elseif sign c(i)=-1; end end if length(bl)>2 if sum(a((1:bl(1))))>sum(a([bl(2)-1:-4:1 bl(2)-2:-4:1])) c(bl(2)-2:-2:1)=-1; c(bl(2)-1:-4:1)=1i; c(bl(2)-3:-4:1)=-1i; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function c = pleats(aarg,type) % by MR Keenan ? al=length(aarg); s=sum(aarg); hcl=al-1; c=ones(hcl,1); if s<3 if s==2 q=find(aarg); brk=floor((q(1)+q(2))/2); c(brk)=-1i; c(brk+1:end)=-1; end return end %remove the a's that have less than one zeros in a row. a = aarg; for j=2:length(aarg)-1 if( aarg(j) == 0 ) if( aarg(j-1) == 1 & aarg(j+1) == 1 ) a(j) = 1; end end end %remove the a's that have less than two zeros in a row. if( type == 2 ) for j=2:length(aarg)-2 if( aarg(j) == 0 & aarg(j+1) == 0 ) if( aarg(j-1) == 1 & aarg(j+2) == 1 ) a(j) = 1; a(j+1) = 1; end end end end c = [ones(length(a)-1,1)]; a = a(:); while sum(diff([0;a;0])==1) > sqrt(sum(a)); a = conv([1;1;1],a); a = a(2:end-1)>0; end; k = sqrt(sum(a)); d = diff([0;a;0]); i1 = find(d==1)-1; N = find(d==-1) - i1 - 1; m = []; for i=1:length(i1); n = max(0, ceil ((N(i)-2.5*k) / k)); r = N(i) - n*k; if r>1.5*k; m = [m, i1(i)+[r/4, r/2-k/2+k*cumsum(ones(1,n)), 3*r/4+k*n]]; else m = [m, i1(i) + r/2]; end; end; cb = zeros(length(a)-1,1); for i=1:length(m)-1; d = (m(i)+m(i+1))/2; b = round(d); cb(b) = 1; m(i+1) = m(i+1)-2*(d-b); end; c = (-1).^cumsum(cb); c(logical(cb))=1i; function setspecialcases global SCM %%%% add strings at your convenience, reverse the strings if you %%%% want to test both directions SCM=[i 1 -i -i -1 -1 i i i 1 1 1 -i -i -i -i -1 -i -1 i ... -1 -1 i i -1 i 1 i i 1 1 i 1 -i 1 1 -i -i 1 -i -1 ... -i -i -i -1 -i -1 -1 -1 i -1 -1 i i -1 i i i 1 i i ... 1 1 i 1 1 1 -i 1 1 -i -i 1 -i -i -i -1 -i -i -i -1 ... -i -1 -1 -1 -1 -1 i -1 -1 i i -1 i i i i i 1; %%% dumbunny1 1 -i -1 -1 i i 1 1 1 -i -i -i -1 -1 -1 -1 i i i i 1 ... 1 1 1 1 -i -i -i -i 1 i i i i i i -1 -i -1 i i -1 ... -i -i -1 i i -1 -i -i -1 i -1 -i -1 -i 1 -i -i -1 ... i -1 -i -i 1 1 -i -1 -i 1 -i 1 i 1 -i -i 1 i i 1 ... -i -i 1 i i 1 -i 1 i 1 i i 1 i -1 i 1 i -1; %%% dumbunny2 1 1 1 1 1 1 1 i -1 -1 -1 -1 -1 -1 -1 -1 i 1 1 1 1 1 ... 1 1 1 1 i -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 i 1 1 1 1 1 ... 1 1 1 1 1 1 i -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 i 1 ... 1 1 1 1 1 1 1 1 1 i -1 -1 -1 -1 -1 -1 -1 -1 -1 i 1 i ... 1 -i 1 1 1 1 1 1 i -1 -1 -1 -1 -1; %%% dumbunny3 % 1 1 1 i -1 -1 -1 -1 -i -i 1 1 1 1 1 i i i -1 -1 -1 -1 ... % -1 -1 -i -i -i -i 1 1 1 1 1 1 1 i i i i i -1 -1 -1 -1 ... % -1 -1 -1 -1 -i -i -i -i -i -i 1 1 1 1 1 1 1 1 1 i i i ... % i i i i -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -i -i -i -i -i ... % -i -i -i 1 1 1 1 1 1 1 1 1 1 1; %%% this spiral works -i 1 -i -1 -i -1 i -1 i 1 i i 1 1 -i 1 -i -i -i -1 -i ... -1 -1 -1 i -1 i i i 1 i i 1 1 1 1 -i 1 -i -i -i -i ... -i -1 -i -1 -1 -1 -1 -1 i -1 i i i i i 1 i i 1 1 1 ... 1 1 1 -i 1 -i -i -i -i -i -i -i -1 -i -1 -1 -1 -1 -1 ... -1 -1 i -1 i i i i i i i 1 i i 1 1 1; %%% cross pattern ].'; function hc = searchnfold(a) global I J I1 J1 global xc yc al al = length(a); s = sum(a); hcl = al-1; hc = zeros(hcl,1); id = ones(s); [I,J] = find(triu(id,1)); Ia = logical(a); iones = find(Ia); I = iones(I); J = iones(J); fsal = floor(al/sqrt(2)); Z = int8(zeros(2*(fsal+2)+1)); xc = fsal+3; yc = fsal+3; k = 0; x = 0; y = 0; Hcc = foldbranch(a,hc,k,x,y,1,Z); hcc = foldbranch(a(al:-1:1),hc,k,x,y,1,Z); Hcc = [Hcc,hcc(hcl:-1:1,:)]; hcf = logical(~any(isnan(Hcc))); if ~any(hcf), hc = ones(hcl,1); return, end Hcc = Hcc(:,hcf); hc = minenergy(Hcc); function Hcc = foldbranch(ha,hc,k,x,y,mine,Z) global xc yc al maxbranches = 8; numbranches = 8; Hcc = []; if k == 0 s = 0; kf = find(ha); kf = kf(1); if kf < 3 k = 1; Z(yc,xc) = 1; if ha(2) hc(k) = 1; hZ = Z; hZ(yc,xc+1) = 1; Hcc = foldbranch(ha,hc,k+1,x+1,y,1,hZ); else Hcc = foldbranch(ha,hc,k,x,y,1,Z); end else hc(1:kf-2) = 1; Z(yc,[max(1,xc-kf+2):xc]) = 1; hc(kf-1) = 1; hZ = Z; hZ(yc,xc+1) = 1; Hcc = foldbranch(ha,hc,kf,x+1,y,1,hZ); end return end if ha(k+1) while (k < al) & ha(k+1) z = Z([y-1:y+1]+yc,[x-1:x+1]+xc); z(1,1) = 1; z(1,3) = 1; z(2,2) = 1; z(3,1) = 1; z(3,3) = 1; if mine bc = inf; else bc = -inf; end if all(z([2,4,6,8])) hc(k) = NaN; % stuck! Hcc = hc; return end hcc = hc; p = [0;cumsum(hcc(1:k-1))]; p = p(logical(ha(1:k))); if ~z(2,3) hcc(k) = 1; phc = p(end)+hcc(k); ec = sum(abs(phc-p)); if (mine & (ec < bc)) | (~mine & (ec > bc)) bc = ec; hc = hcc; xz= 1; yz = 0; end end if ~z(3,2) hcc(k) = -1i; phc = p(end)+hcc(k); ec = sum(abs(phc-p)); if (mine & (ec < bc)) | (~mine & (ec > bc)) bc = ec; hc = hcc; xz= 0; yz = 1; end end if ~z(2,1) hcc(k) = -1; phc = p(end)+hcc(k); ec = sum(abs(phc-p)); if (mine & (ec < bc)) | (~mine & (ec > bc)) bc = ec; hc = hcc; xz= -1; yz = 0; end end if ~z(1,2) hcc(k) = 1i; phc = p(end)+hcc(k); ec = sum(abs(phc-p)); if (mine & (ec < bc)) | (~mine & (ec > bc)) bc = ec; hc = hcc; xz= 0; yz = -1; end end k = k+1; x = x+xz; y = y+yz; mine = 1; Z(y+yc,x+xc) = 1; end if k < al Hcc = foldbranch(ha,hc,k,x,y,1,Z); else Hcc = hc; end else kf = find(ha(k+1:end)); yyc = y+yc; xxc = x+xc; if isempty(kf) zl = al-k; if ~any(Z(yyc,xxc+[1:zl])) hc(k:al-1) = ones(1,zl); Hcc = hc; elseif ~any(Z(yyc+[1:zl],xxc)) hc(k:al-1) = -1i*ones(1,zl); Hcc = hc; elseif ~any(Z(yyc,xxc+[-zl:-1])) hc(k:al-1) = -ones(1,zl); Hcc = hc; elseif ~any(Z(yyc+[-zl:-1],xxc)) hc(k:al-1) = 1i*ones(1,zl); Hcc = hc; else ha(k+1) = 1; Hcc = foldbranch(ha,hc,k,x,y,0,Z); end return end kf = k+kf(1); zl = kf-k-1; zl1 = floor(zl/2); zl2 = ceil(zl/2); zlr = rem(zl,2); if zl == 1 s = 0; if (~Z(yyc,xxc+1)) hc(k) = 1; s = s+1; hZ = Z; hZ(yyc,xxc+1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+1,y,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~Z(yyc+1,xxc)) hc(k) = -1i; s = s+1; hZ = Z; hZ(yyc+1,xxc) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x,y+1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~Z(yyc,xxc-1)) hc(k) = -1; s = s+1; hZ = Z; hZ(yyc,xxc-1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-1,y,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~Z(yyc-1,xxc)) hc(k) = 1i; s = s+1; hZ = Z; hZ(yyc-1,xxc) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x,y-1,1,hZ)]; end if s < 1 hc(k) = NaN; % stuck! Hcc = hc; end return end zl1i = [1:zl1]; zl2i = [1:zl2]; if (y >= 0) & (x >= 0) s = 0; if (~(any(Z(yyc,xxc+zl2i))|any(Z(yyc-1,xxc+zlr+zl1i)))) hc(k:kf-2) = [ones(1,zl2),1i,-ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc,xxc+zl2i) = 1; hZ(yyc-1,xxc+zlr+zl1i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+zlr+1,y-1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc+zl2i,xxc))|any(Z(yyc+zlr+zl1i,xxc-1)))) hc(k:kf-2) = [-1i*ones(1,zl2),-1,1i*ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc+zl2i,xxc) = 1; hZ(yyc+zlr+zl1i,xxc-1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-1,y+zlr+1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc+zl2i))|any(Z(yyc+1,xxc+zlr+zl1i)))) hc(k:kf-2) = [ones(1,zl2),-1i,-ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc,xxc+zl2i) = 1; hZ(yyc+1,xxc+zlr+zl1i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+zlr+1,y+1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc+zl2i,xxc))|any(Z(yyc+zlr+zl1i,xxc+1)))) hc(k:kf-2) = [-1i*ones(1,zl2),1,1i*ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc+zl2i,xxc) = 1; hZ(yyc+zlr+zl1i,xxc+1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+1,y+zlr+1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc+zl1i))|any(Z(yyc-1,xxc-zlr+zl2i)))) hc(k:kf-2) = [ones(1,zl1),1i,-ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc,xxc+zl1i) = 1; hZ(yyc-1,xxc-zlr+zl2i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-zlr+1,y-1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc+zl1i,xxc))|any(Z(yyc-zlr+zl2i,xxc-1)))) hc(k:kf-2) = [-1i*ones(1,zl1),-1,1i*ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc+zl1i,xxc) = 1; hZ(yyc-zlr+zl2i,xxc-1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-1,y-zlr+1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc+zl1i))|any(Z(yyc+1,xxc-zlr+zl2i)))) hc(k:kf-2) = [ones(1,zl1),-1i,-ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc,xxc+zl1i) = 1; hZ(yyc+1,xxc-zlr+zl2i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-zlr+1,y+1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc+zl1i,xxc))|any(Z(yyc-zlr+zl2i,xxc+1)))) hc(k:kf-2) = [-1i*ones(1,zl1),1,1i*ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc+zl1i,xxc) = 1; hZ(yyc-zlr+zl2i,xxc+1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+1,y-zlr+1,1,hZ)]; end if s < 1 ha(k+1) = 1; Hcc = foldbranch(ha,hc,k,x,y,0,Z); end elseif (y >= 0) & (x < 0) s = 0; if (~(any(Z(yyc+zl2i,xxc))|any(Z(yyc+zlr+zl1i,xxc+1)))) hc(k:kf-2) = [-1i*ones(1,zl2),1,1i*ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc+zl2i,xxc) = 1; hZ(yyc+zlr+zl1i,xxc+1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+1,y+zlr+1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc-zl2i))|any(Z(yyc-1,xxc-zlr-zl1i)))) hc(k:kf-2) = [-ones(1,zl2),1i,ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc,xxc-zl2i) = 1; hZ(yyc-1,xxc-zlr-zl1i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-zlr-1,y-1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc+zl2i,xxc))|any(Z(yyc+zlr+zl1i,xxc-1)))) hc(k:kf-2) = [-1i*ones(1,zl2),-1,1i*ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc+zl2i,xxc) = 1; hZ(yyc+zlr+zl1i,xxc-1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-1,y+zlr+1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc-zl2i))|any(Z(yyc+1,xxc-zlr-zl1i)))) hc(k:kf-2) = [-ones(1,zl2),-1i,ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc,xxc-zl2i) = 1; hZ(yyc+1,xxc-zlr-zl1i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-zlr-1,y+1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc+zl1i,xxc))|any(Z(yyc-zlr+zl2i,xxc+1)))) hc(k:kf-2) = [-1i*ones(1,zl1),1,1i*ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc+zl1i,xxc) = 1; hZ(yyc-zlr+zl2i,xxc+1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+1,y-zlr+1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc-zl1i))|any(Z(yyc-1,xxc+zlr-zl2i)))) hc(k:kf-2) = [-ones(1,zl1),1i,ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc,xxc-zl1i) = 1; hZ(yyc-1,xxc+zlr-zl2i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+zlr-1,y-1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc+zl1i,xxc))|any(Z(yyc-zlr+zl2i,xxc-1)))) hc(k:kf-2) = [-1i*ones(1,zl1),-1,1i*ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc+zl1i,xxc) = 1; hZ(yyc-zlr+zl2i,xxc-1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-1,y-zlr+1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc-zl1i))|any(Z(yyc+1,xxc+zlr-zl2i)))) hc(k:kf-2) = [-ones(1,zl1),-1i,ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc,xxc-zl1i) = 1; hZ(yyc+1,xxc+zlr-zl2i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+zlr-1,y+1,1,hZ)]; end if s < 1 ha(k+1) = 1; Hcc = foldbranch(ha,hc,k,x,y,0,Z); end elseif (y < 0) & (x < 0) s = 0; if (~(any(Z(yyc,xxc-zl2i))|any(Z(yyc+1,xxc-zlr-zl1i)))) hc(k:kf-2) = [-ones(1,zl2),-1i,ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc,xxc-zl2i) = 1; hZ(yyc+1,xxc-zlr-zl1i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-zlr-1,y+1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc-zl2i,xxc))|any(Z(yyc-zlr-zl1i,xxc+1)))) hc(k:kf-2) = [1i*ones(1,zl2),1,-1i*ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc-zl2i,xxc) = 1; hZ(yyc-zlr-zl1i,xxc+1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+1,y-zlr-1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc-zl2i))|any(Z(yyc-1,xxc-zlr-zl1i)))) hc(k:kf-2) = [-ones(1,zl2),1i,ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc,xxc-zl2i) = 1; hZ(yyc-1,xxc-zlr-zl1i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-zlr-1,y-1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc-zl2i,xxc))|any(Z(yyc-zlr-zl1i,xxc-1)))) hc(k:kf-2) = [1i*ones(1,zl2),-1,-1i*ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc-zl2i,xxc) = 1; hZ(yyc-zlr-zl1i,xxc-1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-1,y-zlr-1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc-zl1i))|any(Z(yyc+1,xxc+zlr-zl2i)))) hc(k:kf-2) = [-ones(1,zl1),-1i,ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc,xxc-zl1i) = 1; hZ(yyc+1,xxc+zlr-zl2i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+zlr-1,y+1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc-zl1i,xxc))|any(Z(yyc+zlr-zl2i,xxc+1)))) hc(k:kf-2) = [1i*ones(1,zl1),1,-1i*ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc-zl1i,xxc) = 1; hZ(yyc+zlr-zl2i,xxc+1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+1,y+zlr-1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc-zl1i))|any(Z(yyc-1,xxc+zlr-zl2i)))) hc(k:kf-2) = [-ones(1,zl1),1i,ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc,xxc-zl1i) = 1; hZ(yyc-1,xxc+zlr-zl2i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+zlr-1,y-1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc-zl1i,xxc))|any(Z(yyc+zlr-zl2i,xxc-1)))) hc(k:kf-2) = [1i*ones(1,zl1),-1,-1i*ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc-zl1i,xxc) = 1; hZ(yyc+zlr-zl2i,xxc-1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-1,y+zlr-1,1,hZ)]; end if s < 1 ha(k+1) = 1; Hcc = foldbranch(ha,hc,k,x,y,0,Z); end elseif (y < 0) & (x >= 0) s = 0; if (~(any(Z(yyc-zl2i,xxc))|any(Z(yyc-zlr-zl1i,xxc-1)))) hc(k:kf-2) = [1i*ones(1,zl2),-1,-1i*ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc-zl2i,xxc) = 1; hZ(yyc-zlr-zl1i,xxc-1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-1,y-zlr-1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc+zl2i))|any(Z(yyc+1,xxc+zlr+zl1i)))) hc(k:kf-2) = [ones(1,zl2),-1i,-ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc,xxc+zl2i) = 1; hZ(yyc+1,xxc+zlr+zl1i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+zlr+1,y+1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc-zl2i,xxc))|any(Z(yyc-zlr-zl1i,xxc+1)))) hc(k:kf-2) = [1i*ones(1,zl2),1,-1i*ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc-zl2i,xxc) = 1; hZ(yyc-zlr-zl1i,xxc+1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+1,y-zlr-1,1,hZ)]; end if (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc+zl2i))|any(Z(yyc-1,xxc+zlr+zl1i)))) hc(k:kf-2) = [ones(1,zl2),1i,-ones(1,zl1-1)]; s = s+1; hZ = Z; hZ(yyc,xxc+zl2i) = 1; hZ(yyc-1,xxc+zlr+zl1i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+zlr+1,y-1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc-zl1i,xxc))|any(Z(yyc+zlr-zl2i,xxc-1)))) hc(k:kf-2) = [1i*ones(1,zl1),-1,-1i*ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc-zl1i,xxc) = 1; hZ(yyc+zlr-zl2i,xxc-1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-1,y+zlr-1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc+zl1i))|any(Z(yyc+1,xxc-zlr+zl2i)))) hc(k:kf-2) = [ones(1,zl1),-1i,-ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc,xxc+zl1i) = 1; hZ(yyc+1,xxc-zlr+zl2i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-zlr+1,y+1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc-zl1i,xxc))|any(Z(yyc+zlr-zl2i,xxc+1)))) hc(k:kf-2) = [1i*ones(1,zl1),1,-1i*ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc-zl1i,xxc) = 1; hZ(yyc+zlr-zl2i,xxc+1) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x+1,y+zlr-1,1,hZ)]; end if zlr & (s < numbranches) & ((s < 1) | (size(Hcc,2) < maxbranches)) & (~(any(Z(yyc,xxc+zl1i))|any(Z(yyc-1,xxc-zlr+zl2i)))) hc(k:kf-2) = [ones(1,zl1),1i,-ones(1,zl2-1)]; s = s+1; hZ = Z; hZ(yyc,xxc+zl1i) = 1; hZ(yyc-1,xxc-zlr+zl2i) = 1; Hcc = [Hcc,foldbranch(ha,hc,kf-1,x-zlr+1,y-1,1,hZ)]; end end end % Calculate distribution of zeros and ones for % the binary string function da=distribution(a) el=length(a); curr=0; da=zeros(el+2,1); dal=1; for l=1:el if ( a(l)==curr ) da(dal)=da(dal)+1; else curr=~curr; dal=dal+1; da(dal)=da(dal)+1; end end if ( ~mod(dal,2) ), dal=dal+1;,end da=da(1:dal); % Calculate adjusted distribution of zeros and ones % for the binary string such that there is always % an even number of ones and zeros for the sequences % in between the end zeros in the sequence function da=adistribution(da) das=length(da); if ( mod(sum(da(2:end-1)),2) ) if ( da(1)>0 ) da(1:2)=da(1:2)+[-1;1]; elseif ( da(end)>0 ) da(end-1:end)=da(end-1:end)+[1;-1]; else da(end-1:end)=da(end-1:end)+[-1;1]; end end l=2; while (l2 ) % Skip next even l=l+2; if (l>=das-2), l=l+2;break, end else da(l)=da(l)+1; da(l+1)=da(l+1)-1; end end if ( mod(da(l+1),2) ) da(l+1)=da(l+1)-1; da(l+2)=da(l+2)+1; end l=l+2; end if ( mod(da(das-1),2) & l=4) & csqs*(csqs-1)>=num ) CCk=CCk+1;CC(:,CCk)=fillcsquare(ada,csqs,csqs-1); CCk=CCk+1;CC(:,CCk)=fillcsquare(ada,csqs-1,csqs); end ada=adistribution(flipud(da)); CCk=CCk+1;CC(:,CCk)=flipud(fillcsquare(ada,csqs,csqs)); if ( (csqs>=4) & csqs*(csqs-1)>=num ) CCk=CCk+1;CC(:,CCk)=flipud(fillcsquare(ada,csqs,csqs-1)); CCk=CCk+1;CC(:,CCk)=flipud(fillcsquare(ada,csqs-1,csqs)); end CC = [CC(:,1:CCk) CC(hcl:-1:1,1:CCk)]; c = minenergy(CC); % Calculate filling of a compact square where we try to fit the % ones within the compact square and the zeros outside of it function c=fillcsquare(da,cscols,csrows) ci=1; % Index of current position for c stepsize=0; % Step size for each filling % Help counters north(1:cscols,1)=csrows; south(1:cscols,1)=csrows-1; west(1:csrows,1)=cscols-1; % Start by adding first zeros to c if ( da(1)>0 ) c=1i*ones(da(1)-1,1); da(1)=0; istep=1i; else c=[]; da(1)=-2; istep=0; end % Add the bottom row to the square if ( mod(cscols,2) ) [c,da,steps,lsteps,rsteps]=fillcsquareside(1i,c,da,south,cscols,csrows,istep,1); else [c,da,steps,lsteps,rsteps]=fillcsquareside(1i,c,da,south,cscols,csrows,istep,2); end west=min(west,lsteps); east=flipud(rsteps); north=north-flipud(steps); % Then fill the west side of the square if ( west(1) ) istep=-1; else istep=0; end if ( mod(csrows+west(1),2) ) west(end)=1; west(end-1)=1; end [c,da,steps,lsteps,rsteps]=fillcsquareside(1,c,da,west,csrows,cscols,istep,1); east=min(east,flipud(cscols-steps)); north=min(north,lsteps); % Fill top row of the square if ( north(1) ) istep=1i; else istep=0; end if ( mod(cscols+north(1),2) ) north(end)=1; north(end-1)=1; end [c,da,steps,lsteps,rsteps]=fillcsquareside(-1i,c,da,north,cscols,csrows,istep,1); east=min(east,lsteps); % Fill east side of the square if ( east(1) ) istep=1; else istep=0; end [c,da]=fillcsquareside(-1,c,da,east,csrows,cscols,istep,1); % Add the rest of the atoms if ( c(end)==-1 ) c=[c;-1i*ones(sum(da),1)]; else c=[c;ones(sum(da),1)]; end % Fill one direction with zeros and ones. Begin with % the zeros and the continue with ones until the entire % length of the square has been filled. Return the actual % depths to the caller. Fdir indicates in which direction % we are filling the square and istep indicate whether we % should initiate the filling with an initial step. If istep % has been defined, the step is taken directly if there are no % zeros, otherwise the zeros are extended in a direction % 90 degrees with istep (istep*(-i)) function [c,da,steps,lsteps,rsteps]=fillcsquareside(fdir,c,da,depth,size,lrsize,istep,di) adir=1i*fdir; % Direction to move when advancing if ( depth(di)==0 ), di=di+1;, end steps=zeros(size,1); steps(di)=1; lsteps=ones(lrsize,1)*size; rsteps=lsteps; lsteps(1)=size-di; rsteps(1)=di-1; minstep=1; if (length(da)<3), return, end % Check if we should quit if ( depth(di)==0 ) % Quit here. Cannot fill any more of the rectangle c=[c;-fdir*ones(sum(da),1)]; da=[0]; return; end % If istep defined, start with taking the istep if ( abs(istep)>0 ) stepsize=da(1)/2; if ( stepsize > 0 ) if ( c(end)~=(-1i)*istep ) c=[c;1i*istep*ones(stepsize,1);istep;(-1i)*istep*ones(stepsize,1)]; da(1)=-2; else c=[c;istep]; if ( da(end)>0 ) da(end)=da(end)-1; else da(end-1)=da(end-1)-1; end end else c=[c;istep]; da(1)=-2; end end % Fill out ones until we have filled the square; while ( di0 ) if ( (di>=size-1) | (depth(di+2)>0 & depth(di+1)>0 ) ) c=[c;-fdir*ones(stepsize,1);adir;fdir*ones(stepsize,1)]; else c=[c;-fdir*ones(sum(da),1)]; da=[0;0;0]; break; end elseif ( stepsize==0 ) if ( depth(di+1) > 0 ) c=[c;adir]; else c=[c;-fdir*ones(sum(da),1)]; da=[0;0;0]; break; end end da(1)=0; if ( da(2)==0 ), da(3)=da(3)-1;, end while ( da(2)>0 ) stepsize=min(depth(di),depth(di+1)); if ( da(2)/2<=stepsize ) stepsize=da(2)/2; end if ( stepsize==0 ) c=[c;-fdir*ones(sum(da)-1,1)]; da=[0;0;0]; break; end stepsize=floor(stepsize); da(2)=da(2)-2*stepsize; if ( stepsize==0 ) if ( (isempty(c) | c(end)~=fdir) & ((di>=size-1) | depth(di+2)>0) ) c=[c;-fdir*ones(da(3)/2,1);adir;fdir*ones(da(3)/2,1)]; steps(di)=1; steps(di+1)=1; lsteps(1)=size-di-1; if ( minstep==0 ) rsteps(1)=di-1; minstep=1; end di=di+2; da=da(3:end); da(1)=0; da(2)=da(2)-1; else c=[c;adir]; steps(di)=1; steps(di+1)=1; lsteps(1)=size-di-1; if ( minstep==0 ) rsteps(1)=di-1; minstep=1; end da(2)=0; da(3)=da(3)-2; da(4)=da(4)+1; di=di+2; end else steps(di)=stepsize; steps(di+1)=stepsize; c=[c;fdir*ones(stepsize-1,1);adir;-fdir*ones(stepsize-1,1)]; lsteps(1:stepsize)=size-di-1; if ( stepsize>minstep ) rsteps(minstep+1:stepsize)=di-1; minstep=stepsize; end di=di+2; end if ( di>=size ), break, end if ( da(2)>0 ) if ( depth(di)>0 ) c=[c;adir]; else c=[c;-fdir]; end end end if ( da(2)==0 ) da=da(3:end); if ( length(da)==1 ) while ( da(1)>0 & di<=size & depth(di)>0 ) c=[c;adir]; da(1)=da(1)-1; di=di+1; end if (da(1)>0) c=[c;-fdir*ones(da(1),1)]; end da(1)=0; return; end end end```