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<n
p=p+t(n-nn,:)+t(n-nn+1,:);
goods(find(p==0))=0;
nn=nn+2;
end
t=t(:,find(goods));
m=size(t,2);
end
p = [zeros(1,m); cumsum(t,1)];
% Sum of square distances costs only O(iones) operations. /PR
atones=find(at);
secondmoment=uvvar(real(p(atones,:)))+uvvar(imag(p(atones,:)));
[scratch,ixlist] = sort(secondmoment);
ixlist=ixlist(1:5);
if ~dl
c=t(:,ixlist(1));
return
end
for m=1:length(ixlist);
c=t(:,ixlist(m));
if head
ct=[ones(dl,1)*[1 1i -1i -1]; c*ones(1,4)];
else
ct=[c*ones(1,4); ones(dl,1)*[1 1i -1i -1]];
end
for n=1:4
p=[0; cumsum(ct(:,n))];
pD=diff(sort(p));
if all(pD)
c=ct(:,n);
return;
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% End early termination cases %%%%%%%%%%%%%%%%%%%
c=zeros(M-1,1);
global I J iones
global C hv jv K Kh Kj SCM
global aori aoriM aoriN1 aoriP aorif;
aori=logical(a);aoriM=M;aoriN1=sum(a);aoriP=sum(1:aoriN1-1);
aorif=find(aori);iones=aorif;
[I,J] = find(triu(ones(aoriN1),1));Ia=logical(a);
I=aorif(I);J=aorif(J);
%"Hydrophilic go away" this is the basic solver didn't work
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
% OVERALL COMPARISON
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
CC=zeros(M-1,12);
CCk=0;
%fillrows by LAC
CCk=CCk+2;
CC(:,[CCk-1:CCk])=[fillrows(aori,-1) fillrows(aori,0)];
%magicball by LAC
if M>=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)<nro tem=tem(end:-1:1,:); end
voe(j:j+mu(k)-1)=tem;
j=j+mu(k);k=k+1;
end
switch sum(a([1 1 M]))
case 0, vze(cumsum(mu))=vz(2:end);
case 1, vze(cumsum(mu(1:end-1)))=vz(2:end);
case 2, vze(cumsum(mu))=vz;
case 3, vze(cumsum(mu(1:end-1)))=vz;
end
voe2=voe-1;
voe2(2:2:end)=-voe2(2:2:end);
tra=ceil((vze+1)/2);
va=ones(sum(mu),1);
nva=testtra(tra,voe2,vze,1,0);
ava=inf;
while nva<ava
ava=nva;va(:)=inf;
whe1=find(vze+1-tra);
for i=1:length(whe1)
va(i)=testtra(tra,voe2,vze,whe1(i),1);
end
[nva,h]=min(va);
if nva<=ava
tra(whe1(h))=tra(whe1(h))+1;
end
end
nva=testtra(tra,voe2,vze,1,0);
ava=inf;
while nva<ava
ava=nva;va(:)=inf;
whe1=find(tra-1);
for i=1:length(whe1)
va(i)=testtra(tra,voe2,vze,whe1(i),-1);
end
[nva,h]=min(va);
if nva<=ava
tra(whe1(h))=tra(whe1(h))-1;
end
end
ltra=length(tra);
lv=zeros(ltra,1);rv=lv;
for j=1:ltra-1
lv(j)=rv(j)+voe2(j);
rv(j+1)=lv(j)+tra(j)+tra(j)-vze(j)-2;
end
lv(ltra)=rv(ltra)+voe2(ltra);
vze2=vze;
vze2(2:2:end)=-vze2(2:2:end);
pui=[(lv(1:nro-1)+rv(2:nro)+vze2(1:nro-1))/2;lv(nro)];
locs=abs(diff([0;pui]))+1;
if ~a(1) locs(1)=locs(1)+vz(1); end
if ~a(M) locs(end)=locs(end)+vz(end); end
locs=cumsum(locs);
c=ones(M-1,1);
c(locs(1:end-1))=1i;
lp=locs(1:2:end);
rp=locs(2:2:end);
for i=1:length(rp)
c(lp(i)+1:rp(i)-1)=-1;
end
function va=testtra(tra,voe2,vze,w,s)
tra(w)=tra(w)+s;ltra=length(tra);
lv=zeros(ltra,1);rv=lv;
for j=1:ltra-1
lv(j)=rv(j)+voe2(j);
rv(j+1)=lv(j)+tra(j)+tra(j)-vze(j)-2;
end
lv(ltra)=rv(ltra)+voe2(ltra);
va=myvar(lv+rv);
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
% From this point ideas from others the above written
% all by Lucio Andrade
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
function c = solverpira(a)
global I J iones
al=length(a);
s=sum(a);
c=ones((al-1),1);
hcl=al-1;
ind1=1:hcl;
revind=al:-1:1;
revind1=hcl:-1:1;
id = ones(s);
[I,J] = find(triu(id,1));
Ia=logical(a);
iones=find(Ia);
I=iones(I);
J=iones(J);
iones1=iones(1);
iones2=iones(length(iones))-1;
ionesind=iones1:iones2;
a1=a(iones1:iones2+1);
a1l=length(a1);
a1l1=a1l-1;
rev1ind=a1l:-1:1;
rev1ind1=a1l1:-1:1;
csal=ceil(sqrt(al));
bc=Inf;
C=[zeros(hcl,2*al)];
Ck=1;
if 1
% zigzag
hc=1i+zeros(hcl,1);
hc1=hc;
for n=max(2,min(a1l1,3)):min([13,csal+2,floor(a1l/3)]);
t=[ones(n-1,1);1i; -ones(n-1,1);1i];
t=t(:,ones(floor(a1l/n),1));
t=t(:);
t=t(1:a1l1);
hc(ionesind)=t;
C(:,Ck)=hc;
Ck=Ck+1;
t=improvezigzagend(t,a1,n-1,1);
if ~isempty(t)
hc1=hc;
hc1(ionesind)=t;
if t(a1l1)==-1i
i=find(real(t));
hc1(iones2+1:hcl)=t(i(length(i)));
end
C(:,Ck)=hc1;
Ck=Ck+1;
end
t=zigzag3(a1,n);
hc1(1:iones1-1)=1;
hc1(ionesind)=t;
i=find(real(t));
hc1(iones2+1:hcl)=t(i(length(i)));
C(:,Ck)=hc1;
Ck=Ck+1;
t=improvezigzagend(t,a1,n,0);
if ~isempty(t)
hc1(ionesind)=t;
i=find(real(t));
hc1(iones2+1:hcl)=t(i(length(i)));
C(:,Ck)=hc1;
Ck=Ck+1;
end
t=zigzag3(a1(rev1ind),n);
i=find(real(t));
hc1(1:iones1-1)=t(i(length(i)));
hc1(ionesind)=t(rev1ind1);
hc1(iones2+1:hcl)=1;
C(:,Ck)=hc1;
Ck=Ck+1;
end
end
% random walk
for jgh=1:ceil(al/8),
for k1=1:ceil(ceil(al/8)/4),
hc=ones(hcl,1);
hc(k1+1)=1i;
hc(k1+2:2*k1+2)=-1;
hc(2*k1+3)=-1i;
k=2*k1+4; d=-1i; pos=-1; m=(k1+1)+1i;
while k<=hcl,
if randn<=0,
hc(k)=d;
m=m+abs(real(d))+1i*abs(imag(d));
pos=pos+d;
q=sign(real(pos))*sign(imag(pos));
if real(d)==0;
d=1i*d*q;
hc(k+1:hcl)=d;
k=k+real(m)+1;
pos=pos+d*real(m);
else
d=-1i*d*q;
hc(k+1:hcl)=d;
k=k+imag(m)+1;
pos=pos+d*imag(m);
end
else
m=m+1i*abs(real(d))+abs(imag(d));
q=sign(real(pos))*sign(imag(pos));
if q==0;
q=-sign(real(pos));
end
if real(d)==0;
d=-1i*d*q;
pos=pos+d;
hc(k)=d;
d=-1i*d*q;
hc(k+1:hcl)=d;
k=k+imag(m)+1;
pos=pos+d*imag(m);
else
d=1i*d*q;
pos=pos+d;
hc(k)=d;
d=1i*d*q;
hc(k+1:hcl)=d;
k=k+real(m)+1;
pos=pos+d*real(m);
end
end
end
hc=hc(ind1);
C(:,Ck:Ck+1)=[hc hc(revind1)];
Ck=Ck+2;
end
end % random walk
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% snake disc
hc = snakedisc(al);
hc=hc(revind1);
C(:,Ck)=hc;
Ck=Ck+1;
hc = rlesnakeold(a(revind));
hc=hc(revind1);
C(:,Ck)=hc;
Ck=Ck+1;
cc=c;
if s<a1l
hc = rlesnake(a1,s);
cc(:)=1i;
cc(ionesind)=hc;
C(:,Ck)=cc;
Ck=Ck+1;
hc = rlesnake(a1(rev1ind),s);
hc=hc(rev1ind1);
cc(:)=1i;
cc(ionesind)=hc;
C(:,Ck)=cc;
Ck=Ck+1;
end
for kk = 1:12
hc = snakesquare(a);
C(:,Ck)=hc;
Ck=Ck+1;
end
for kk = 1:12
hc = snakesquare(a(revind));
hc = hc(revind1);
C(:,Ck)=hc;
Ck=Ck+1;
end
[c,bc]=minenergy(C(:,1:Ck-1));
[c,bc] = msolver(a,c,bc);
%Sragner Laszlo sragner@mit.bme.hu
t=[1 -1 1i -1i];
imp=1;
hossz=length(c);
while imp
p=[0; cumsum(c)];
t2=t+p(end);
C(:,1)=c;
Ck=2;
eleje=-t(find(~(sum(p(:,[1 1 1 1])==t(ones(length(p),1),:)))));
vege=t(find(~(sum(p(:,[1 1 1 1])==t2(ones(length(p),1),:)))));
if ((~isempty(eleje))&(~isempty(vege)))
nyom=1;
else
nyom=0;
end
if nyom
for i=2:hossz-1
if (c(i-1)==-c(i+1))
for el=eleje
for ve=vege
C(:,Ck)=[el;c(1:i-2);c(i);c(i+2:hossz);ve];
Ck=Ck+1;
end
end
end
end
[c,e,imp]=minenergy(C(:,1:Ck-1));
else
imp=0;
end
end
t=[1 -1 1i -1i];
imp=1;
while imp
p=[0; cumsum(c)];
t2=t+p(end);
C(:,1)=c;
Ck=2;
eleje=-t(find(~(sum(p(:,[1 1 1 1])==t(ones(length(p),1),:)))));
vege=t(find(~(sum(p(:,[1 1 1 1])==t2(ones(length(p),1),:)))));
if ((~isempty(eleje))&(~isempty(vege)))
nyom=1;
else
nyom=0;
end
x=real(p);
y=imag(p);
X1=x(:,ones(1,al));
X2=X1';
Y1=y(:,ones(1,al));
Y2=Y1';
dY=(Y1==Y2);
dX=(X1==X2);
maxX=find(max(X1.*dY)==x');
xI1=maxX(find(diff(maxX)==1));xI1=xI1(~isreal(c(xI1))&xI1>2&xI1<hcl-1);
minX=find(min(X1+(~dY*1000))==x');
xI2=minX(find(diff(minX)==1));xI2=xI2(~isreal(c(xI2))&xI2>2&xI2<hcl-1);
maxY=find(max(Y1.*dX)==y');
yI1=maxY(find(diff(maxY)==1));yI1=yI1(isreal(c(yI1))&yI1>2&yI1<hcl-1);
minY=find(min(Y1+(~dX*1000))==y');
yI2=minY(find(diff(minY)==1));yI2=yI2(isreal(c(yI2))&yI2>2&yI2<hcl-1);
for i=xI1,
C(:,Ck)=[c(2:i-1);1 ;c(i) ;-1;c(i+1:hcl-1)];
Ck=Ck+1;
if (nyom)&(c(i-1)==-c(i+1))
for el=eleje
for ve=vege
C(:,Ck)=[el;c(1:i-2);c(i);c(i+2:hcl);ve];
Ck=Ck+1;
end
end
end
end
for i=xI2,
C(:,Ck)=[c(2:i-1);-1 ;c(i) ;1;c(i+1:hcl-1)];
Ck=Ck+1;
if (nyom)&(c(i-1)==-c(i+1))
for el=eleje
for ve=vege
C(:,Ck)=[el;c(1:i-2);c(i);c(i+2:hcl);ve];
Ck=Ck+1;
end
end
end
end
for i=yI1,
C(:,Ck)=[c(2:i-1);1i ;c(i) ;-1i;c(i+1:hcl-1)];
Ck=Ck+1;
if (nyom)&(c(i-1)==-c(i+1))
for el=eleje
for ve=vege
C(:,Ck)=[el;c(1:i-2);c(i);c(i+2:hcl);ve];
Ck=Ck+1;
end
end
end
end
for i=yI2,
C(:,Ck)=[c(2:i-1);-1i ;c(i) ;1i;c(i+1:hcl-1)];
Ck=Ck+1;
if (nyom)&(c(i-1)==-c(i+1))
for el=eleje
for ve=vege
C(:,Ck)=[el;c(1:i-2);c(i);c(i+2:hcl);ve];
Ck=Ck+1;
end
end
end
end
[c,e,imp]=minenergy(C(:,1:Ck-1));
end
imp=1;
t=[1 -1 1i -1i];
while imp
p=[0; cumsum(c)];
t2=t+p(end);
a1=~(sum(p(:,[1 1 1 1])==t(ones(length(p),1),:)));
a2=~(sum(p(:,[1 1 1 1])==t2(ones(length(p),1),:)));
[c,e,imp]=minenergy([c [-t(find(a1));c(1:hcl-1,ones(1,sum(a1)))] [c(2:hcl,ones(1,sum(a2)));t(find(a2))]]);
end
a=randn(1,63);
%_____________________________________________
function [c1,me] = msolver(a,c1,me)
global I J
al=length(a);
hc=ones((al-1),1);
r0 = sum(diff(a)>0);
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 result<me
c1=c2;
me=result;
end
%%% Snake Disc.
function c = snakedisc(N)
R = sqrt(N/pi);
c = [];
dd = 1;
ll = 0;
for hh = -R:R
ll = max(round(abs(halfchord(R,hh+0.5))),1);
c = [ c;dd(ones(1,ll-1),1);1i;-dd(ones(1,ll),1) ];
dd = -dd;
end
if length(c) < N-1; c(end+1:N-1)=c(end); end;
c = c(1:N-1);
function [c,e,valt]=minenergy(c)
global I J iones
if(size(c,2)<2)
return
end
p = [zeros(1,size(c,2));cumsum(c)];
% Sum of square distances costs only O(iones) operations. /PR
if(size(c,2)>3)
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 - 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 n<size(edges,1)-3
if edges(n+3)-edges(n)<S1
% combine rows to one
edges(n+1:n+2)=[];
else
n=n+2;
end
end
edges = diff(edges);
k=edges(1);
c = ones(N-1,1);
if k>S
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 ap<pminrev % not yet reversing
if ap+nN<pminrev % don't reverse at all
c(k+1:k+nN-1)=dd;
p=p+dd*nN;
else % reverse after zeros
c(k+1:k+n-1)=dd;
c(k+n)=1i;
dd=-dd;
c(k+n+1:k+nN-1)=dd;
p=p+dd*(n-N-1);
end
else % reverse
c(k)=1i;
dd=-dd;
c(k+1:k+n-1+N)=dd;
p=p+dd*(nN-1);
end % reversed
else % step out
p=p+dd*(n-N+1);
c(k:k+n-1)=dd;
c(k+n)=1i;
dd=-dd;
c(k+n+1:k+nN-1)=dd;
end
else
n1=floor(adpd);
n2=floor((n-n1)/2);
if dpd>0
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
c(k)=d;
k=k+1;
pos=pos+1;
end
while pos<l
n=find(a(k+1:al)==0);
if isempty(n)|pos+n(1)>=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<l
c(k)=d;
pos=pos+1;
k=k+1;
end
end
end % while
end % if
end
if pos<l
c(k:k+l-pos-1)=d;
k=k+l-pos;
end
else % not the first line
c(k:k+l-1)=d;
k=k+l;
end
if k==al-1
c(k)=1i;
break;
elseif k>=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 (l<das-2)
if ( mod(da(l),2) )
if ( mod(da(l+2),2) & ~mod(da(l+1),2) & da(l+2)>2 )
% 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<das )
da(das-1)=da(das-1)+1;
da(das-2)=da(das-2)-1;
end
for l=das-2:-2:3
if ( da(l)==0 )
da(l-1)=da(l-1)+da(l+1);
da=[da(1:l-1);da(l+2:end)];
end
end
if ( da(end-1)==0 )
da(end-2)=da(end-2)+da(end);
da=da(1:end-2);
end
% Calculate compact square representation for the genetic string
function c=calccsquare(a)
hcl=length(a)-1;
CC=zeros(hcl,12);
CCk=0;
da=distribution(a);
ada=adistribution(da);
num=sum(ada(2:2:end));
csqs=ceil(sqrt(num));
csqs=max(csqs,3);
CCk=CCk+1;CC(:,CCk)=fillcsquare(ada,csqs,csqs);
if ( (csqs>=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 ( di<size )
% Start by extending zeros
stepsize=da(1)/2;
if ( stepsize>0 )
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
|
