Finish 2002-11-12 00:00:00 UTC

SuperTweak4

by SuperTweak

Status: Passed
Results: 2160.360
CPU Time: 73.556
Score: 2166.01
Submitted at: 2002-11-13 16:19:23 UTC
Scored at: 2002-11-13 16:31:35 UTC

Current Rank: 58th
Based on: FF01 (diff)
Basis for: SuperTweak4b (diff)
Basis for: SuperTweak5 (diff)
Basis for: SuperTweak6b (diff)

Comments
Please login or create a profile.
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);
if M == 36,
  c = [ -i -i -1 -1 i -1 -i -1 -1 -i 1 1 1 -i -i -1 -i 1 -i 1 i i i i 1 -i -i 1 1 1 i -1 -1 i 1]'; 
  return
elseif M == 34
 %a =[ 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1]'
  c = [ i 1 -i -i -i -i -i -i 1 i i i i i i 1 1 i i i -1 -i -i -1 i -1 -1 -1 -1 -i 1 1 1 ]';
  return
elseif M == 10,
 %a = [ 1 0 0 0 0 1 1 1 0 1]'
  c = [ 1 1 -i -1 -1 -1 i i 1 ]';
  return
end

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(hcl,10);
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

SCM=[ -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
      ].';
CCk=CCk+1;CC(:,CCk)= SCM(1:hcl,1);

%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

% new line from Mr. Bond
%CC = [CC(:,1:CCk) CC(hcl:-1:1,1:CCk)];

c = minenergy(CC(:,1:CCk));
imp=1;
al=length(a);
hcl=length(c);
t=[1 -1 1i -1i];
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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
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);

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 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 = 16;
numbranches = 4;

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 z(2) && z(4) && z(6) && z(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