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
global I J % I1 J1
al=length(a);
s=sum(a);
hcl=al-1;
c=ones(hcl,1);
if s<3
if s==2
q=find(a);
brk=floor((q(1)+q(2))/2); c(brk)=-1i; c(brk+1:end)=-1;
end
return
end
% Permutator, written by Yuval Cohen
aOnes=find(a);
aCropL=min(al-aOnes(1)+1,aOnes(end));
if aCropL<13
dl=al-aCropL;
if dl
head = aOnes(1)-1 >= al-aOnes(end);
if head
at=a(aOnes(1):end);
else
at=a(1:aOnes(end));
end
else
at=a;
end
ctl=aCropL-1;
m=1;
t=[ones(ctl,1) ones(ctl,3^(ctl-1)-1)*(1+1i)];
for n=2:ctl;
bendLeft=[ones(n-1,1); 1i+zeros(ctl-n+1,1)]*ones(1,m);
bendRight=[ones(n-1,1); -1i+zeros(ctl-n+1,1)]*ones(1,m);
t(:,m+1:2*m)=t(:,1:m).*bendLeft;
t(:,2*m+1:3*m)=t(:,1:m).*bendRight;
m=3*m;
end
p = zeros(aCropL,m);
p(2:aCropL,:)=cumsum(t,1);
pS=sort(p,1);
pD=diff(pS,1,1);
[n,bads]=find(pD==0);
p(:,bads)=[];
t(:,bads)=[];
pE=p(logical(at),:);
id = logical(ones(s));
[It,Jt] = find(triu(id,1));
e=sum(abs(pE(It,:)-pE(Jt,:)));
if ~dl
[m,n]=min(e);
c=t(:,n);
return
end
[m,emin]=sort(e);
m=1;
while m<length(e);
c=t(:,emin(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))];
pS=sort(p);
pD=pS(2:end)-pS(1:end-1);
if all(pD)
c=ct(:,n);
return;
end
end
m=m+1;
end
end
if s==3
q=find(a);
if (q(2)==q(1)+1)
c(q(1))=1; c(q(2))=1;
elseif (q(2)==q(1)+2)
c(q(1))=1i; c(q(1)+1)=1; c(q(2))=-1i;
else
brk=floor((q(1)+q(2))/2); c(q(1):brk)=1i; c(brk)=1; c(brk+1:q(2))=-1i;
end
beven=(mod(q(2)-q(1),2)==0);
if (q(3)==q(2)+1)
c(q(2):end)=-1i;
elseif (q(3)==q(2)+2)
if (beven)
c(q(2))=1; c(q(2)+1:end)=-1i;
else
c(q(2))=-1i; c(q(2)+1:end)=-1;
end
else
brk=floor((q(2)+q(3))/2); c(q(2):brk)=1; c(brk)=-1i;
c(brk+1:q(3)-1)=-1; c(q(3)-1:end)=-1i;
end
return
elseif s==4
q=find(a);
if (q(2)==q(1)+1)
c(q(1))=1; c(q(2))=1;
elseif (q(2)==q(1)+2)
c(q(1))=1i; c(q(1)+1)=1; c(q(2))=-1i;
else
brk=floor((q(1)+q(2))/2); c(q(1):brk)=1i; c(brk)=1;
c(brk+1:q(2))=-1i;
end
if (q(3)==q(2)+1)
c(q(2))=-1i; c(q(3))=-1;
elseif (q(3)==q(2)+2)
c(q(2))=1; c(q(2)+1)=-1i; c(q(3))=-1;
else
brk=floor((q(2)+q(3))/2); c(q(2):brk)=1; c(brk)=-1i;
c(brk+1:q(3))=-1;
end
if (q(4)==q(3)+1)
c(q(3))=-1; c(q(4)-1:end)=-1;
elseif (q(4)==q(3)+2)
c(q(3))=-1i; c(q(3)+1)=-1; c(q(4)-1:end)=-1;
else
brk=floor((q(3)+q(4))/2); c(q(3):brk)=-1i; c(brk)=-1;
c(brk+1:q(4)-1)=1i; c(q(4)-1:end)=-1;
end
return
end
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(end)-1;
ionesind=iones1:iones2;
a1=a(iones1:iones2+1);
a1l=length(a1);
a1l1=a1l-1;
rev1ind=a1l:-1:1;
rev1ind1=a1l1:-1:1;
%I1=I-iones1+1;
%J1=J-iones1+1;
csal=ceil(sqrt(al));
bc=Inf;
dum=[ 1; 1i; 1;-1i; 1; -1i; -1; -1i; -1; 1i;-1;-1; 1i; 1i; 1; ...
1i; 1; 1; 1; -1i; 1; -1i; -1i; -1i;-1; -1i;-1;-1;-1; ...
1i;-1;-1; 1i; 1i; 1i; 1i; 1; 1i; 1; 1; 1; 1; 1; -1i; ...
1; -1i; -1i; -1i; -1i; -1i;-1; -1i;-1;-1;-1;-1;-1; 1i; ...
-1;-1; 1i; 1i; 1i; 1i; 1i; 1i; 1; 1i; 1; 1; 1; 1; 1; ...
1; 1; -1i; 1; -1i; -1i; -1i; -1i; -1i; -1i; -1i;-1; -1i; ...
-1;-1;-1;-1;-1;-1;-1; 1i;-1;-1; 1i; 1i; 1i; 1i; 1i; ...
1i; 1i; 1i; 1; 1i; 1; 1; 1; 1; 1; 1; 1; 1; 1; -1i; ...
1; -1i; -1i; -1i; -1i; -1i; -1i; -1i; -1i; -1i;-1; -1i;-1; ...
-1;-1;-1;-1;-1;-1;-1;-1; 1i;-1;-1; 1i; 1i; 1i; 1i;...
1i; 1i; 1i; 1i; 1i; 1i; 1; 1i; 1; 1; 1; 1; 1; 1; 1; ...
1; 1; 1; 1; -1i; 1; -1i; -1i; -1i; -1i; -1i; -1i; -1i; ...
-1i; -1i; -1i; -1i;-1; -1i;-1;-1;-1;-1;-1;-1;-1;-1; ...
-1;-1;-1; 1i;-1;-1; 1i; 1i; 1i; 1i; 1i; 1i; 1i; 1i; 1i; ...
1i; 1i; 1i; 1; 1i; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; ...
1; 1; 1; -1i; 1; -1i; -1i; -1i; -1i; -1i; -1i];
C=[dum(ind1) dum(revind1) zeros(hcl,250)];
Ck=3;
% 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(end));
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(end));
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(end));
C(:,Ck)=hc1;
Ck=Ck+1;
end
t=zigzag3(a1(rev1ind),n);
i=find(real(t));
hc1(1:iones1-1)=t(i(end));
hc1(ionesind)=t(rev1ind1);
hc1(iones2+1:hcl)=1;
C(:,Ck)=hc1;
Ck=Ck+1;
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
dum=[-1i;1;1i;1;1i;-1;1i;-1;-1i;-1;-1i;-1;-1i;1;-1i;1;-1i;1;1i;1;...
1i;1;1i;1;1i;-1;1i;-1;1i;-1;1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;...
-1i;1;-1i;1;-1i;1;-1i;1;-1i;1;1i;1;1i;1;1i;1;1i;1;1i;1;...
1i;-1;1i;-1;1i;-1;1i;-1;1i;-1;1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;...
-1i;-1;-1i;-1;-1i;1;-1i;1;-1i;1;-1i;1;-1i;1;-1i;1;-1i;1;1i;1;...
1i;1;1i;1;1i;1;1i;1;1i;1;1i;1;1i;-1;1i;-1;1i;-1;1i;-1;...
1i;-1;1i;-1;1i;-1;1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;...
-1i;-1;-1i;-1;-1i;1;-1i;1;-1i;1;-1i;1;-1i;1;-1i;1;-1i;1;-1i;1;...
-1i;1;1i;1;1i;1;1i;1;1i;1;1i;1;1i;1;1i;1;1i;1;1i;1;...
1i;-1;1i;-1;1i;-1;1i;-1;1i;-1;1i;-1;1i;-1;1i;-1;1i;-1;1i;-1;...
-1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;-1i;-1;...
-1i;1;-1i;1;-1i;1;-1i];
C(:,Ck:Ck+1)=[dum(ind1) dum(revind1)];
Ck=Ck+2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
% rle snake
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);
p=[0; cumsum(c)];
C(:,1)=c;
Ck=2;
x=real(p);
y=imag(p);
jc=c;
for i = 3:length(c)-2
if ~isreal(c(i))
if (max(x(y==y(i)))==x(i)) & (max(x(y==y(i+1)))==x(i+1))
c2=[c(2:i-1);1 ;c(i) ;-1;c(i+1:hcl-1)];
C(:,Ck)=c2;
Ck=Ck+1;
end
if (min(x(y==y(i)))==x(i)) & (min(x(y==y(i+1)))==x(i+1))
c2=[c(2:i-1);-1 ;c(i) ;1;c(i+1:hcl-1)];
C(:,Ck)=c2;
Ck=Ck+1;
end
else
if (max(y(x==x(i)))==y(i)) & (max(y(x==x(i+1)))==y(i+1))
c2=[c(2:i-1);1i ;c(i); -1i;c(i+1:hcl-1)];
C(:,Ck)=c2;
Ck=Ck+1;
end
if (min(y(x==x(i)))==y(i)) & (min(y(x==x(i+1)))==y(i+1))
c2=[c(2:i-1);-1i; c(i); 1i;c(i+1:hcl-1)];
C(:,Ck)=c2;
Ck=Ck+1;
end
end
end
c=minenergy(C(:,1:Ck-1));
t=[1 -1 1i -1i];
p=[0; cumsum(c)];
a1=~(sum(p(:,[1 1 1 1])==t(ones(length(p),1),:)));
c=minenergy([c [-t(find(a1));c(1:hcl-1,ones(1,sum(a1)))]]);
a=randn(1,63);
%_____________________________________________
function [c1,me] = msolver(a,c1,me)
al=length(a);
hcl=al-1;
hc=ones(hcl,1);
r0 = sum(diff(a)>0);
s = ceil(sqrt(r0+sum(a)));
f = find(a);
k = 10;
step = 1;
while k > 0
k = k-step;
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 sign(ones(1,t0+s+t1)) 1i];
t0 = t1;
sign = -sign;
end
if (length(c) < al-1)
c = [c -sign(ones(1,al-1-length(c)))];
else
c = c(1:(al-1));
end
e=energy(c(:));
if e<me
if me-e > 1
k = min(10,k+4);
end
c1 = c(:);
me = e;
end
end
c = ones(al-1,1);
n = cumsum([3 10:8:al-1]);
temp = min(find(n>al-1))-1;
if ~isempty(temp)
n = n(1:temp);
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 = min(find(n>al-1))-1;
if ~isempty(temp)
n = n(1:temp);
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 = conv([1;1;1],a);
a = a(2:end-1)>0;
end;
d = diff([0;a;0]);
i1 = find(d>0)-1;
N = find(d<0) - i1 - 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;
result = energy(c2);
if result<me
c1=c2;
me=result;
end
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]=minenergy(c)
global I J
p = [zeros(1,size(c,2));cumsum(c)];
[e,x]=min(sum(abs(p(I,:)-p(J,:))));
c=c(:,x);
%%% Energy
function ec=energy(c)
global I J
p = [0;cumsum(c)];
ec=sum(abs(p(I)-p(J)));
%%% Energy1
% function ec=energy1(c)
% global I1 J1
% p = [0;cumsum(c)];
% ec=sum(abs(p(I1)-p(J1)));
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=0;
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)<2
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
if godown
p=sum(c);
pr=real(p);
if pr<-1&imag(p)+pr>0
j=al+pr+1;
c(j:al-1)=-1i;
elseif pr>l+2
j=al-pr+l+1;
c(j:al-1)=-1i;
end
end
|
