%%%%%read in%%%%%%%
clear all;
%fd=fopen('in_100.txt');
%fd=fopen('test2.txt');
fd=fopen('randpts_50_3');
data=fscanf(fd,'%d');
fclose(fd);

%fd=fopen('congestion.txt');
%congestion=fscanf(fd,'%d');
%fclose(fd);

%[a,b]=size(congestion);
%congestionx=congestion(1:x/2,:);
%congestiony=congestion(x/2+1:end,:);

s=1;
delta=0;
N0=10;

for i=2:2:length(data)
    x(i/2)=s*(data(i)+delta);
    y(i/2)=s*(data(i+1)+delta);
end

Output = [0];
h = 1e-11;
Num_Iterations = 1000;
MeshSize = (max(max(x),max(y))+2*delta); % MeshSize by MeshSize mes
MatrixSize = MeshSize * MeshSize;

t=h:h:h*Num_Iterations;
%u=sin(2*pi/5/h*t);
%u=[0 zeros(1,length(t)-1)];
u=ones(Num_Iterations,1);
%u=[0 zeros(1,length(t)-1)];
G = sparse([],[],[],MatrixSize,MatrixSize);
C = sparse([],[],[],MatrixSize,MatrixSize);
N=MatrixSize;

%initialize variables in G matrix due to resistors in mesh (other than
%diagonal of matrix)
for i=1:N 
    if (i+1<=N)
        G(i,i+1)=-1;
        G(i+1,i)=-1;
        G(i,i)=G(i,i)+1;
        G(i+1,i+1)=G(i+1,i+1)+1;
    end
    if (i-1>0)
        G(i,i-1)=-1;
        G(i-1,i)=-1;
        G(i,i)=G(i,i)+1;
        G(i-1,i-1)=G(i-1,i-1)+1;
    end
    if  (i+sqrt(N)<=N)
        G(i,i+sqrt(N))=-1;
        G(i+sqrt(N),i)=-1;
        G(i,i)=G(i,i)+1;
        G(i+sqrt(N),i+sqrt(N))=G(i+sqrt(N),i+sqrt(N))+1;
    end
    if (i-sqrt(N)>0)
        G(i,i-sqrt(N))=-1;
        G(i-sqrt(N),i)=-1;
        G(i,i)=G(i,i)+1;
        G(i-sqrt(N),i-sqrt(N))=G(i-sqrt(N),i-sqrt(N))+1;
    end
end

%initialize diagonal of G and C matrices
 for i=1:MatrixSize
     %G(i,i)=5;
     C(i,i)=1e-10;
 end

%input values of B
B=sparse([]);
i=1;
j=1;
B=sparse([],[],[],MatrixSize,1);
L=sparse([]);
for i=1:length(x)
    B((x(i)-1)*MeshSize+y(i))=1;
end



L=sparse(1:MatrixSize,1:MatrixSize,ones(1,MatrixSize),MatrixSize,MatrixSize);

%%%%%%model order reduction%%%%%%%
if (MatrixSize>=1e5)
    G=sparse(G);
    C=sparse(C);
    B=sparse(B);
    L=sparse(L);
    q=40; %reduced order
    %A=inv(G);
    %V=A*B;
    %V=V/norm(V);
    %for i=1:q-1
    %    V0=-A*C*V(:,end);
    %    for j=1:i-size(V,2)
    %        V0=V0-(V0'*V(:,j))*V(:,j);
    %    end
   %     if(norm(V0)<1e-10)
   %         continue;
   %     end
  %      V=[V,V0/norm(V0)];
  %  end
  
    Vi(1).V= G\B;
    [Vi(1).V,R] = qr(Vi(1).V,0);
    for j=1:q-1
        W = -G\(C*Vi(j).V);
        for i=1:j
            Hi(i,j).H = Vi(i).V.'*W;
            W = W - Vi(i).V*Hi(i,j).H;
        end
        for i=1:j
            Rho = Vi(i).V.'*W;
            W = W - Vi(i).V*Rho;
            Hi(i,j).H = Hi(i,j).H + Rho;
        end;
        [Vi(j+1).V,Hi(i+1,j).H]=qr(W,0);
    end; 
    Vt = []; 
    for i=1:q
        Vt=[Vt,Vi(i).V];
    end
    Gq = V.'*G*V; 
    Cq = V.'*C*V;
    Bq = V.'*B; 
    Lq = V.'*L;  
else
    Gq=G;
    Cq=C;
    Bq=B;
    Lq=L;
end
%Lq=L;
   
%Xnew=sparse([],[],[],MatrixSize,1);
Xnew=B;

%u(1)=0;
%Xnew=B;

%sequence=rand(40,1);
%sort(sequence);


%Xnew = (Gq+ ((1/h) * Cq)) \ ((((1/h) * Cq) * Xold) + Bq);
Xold=zeros(MatrixSize,1);
sequence=[];
flag=0;
for i=1:Num_Iterations
    if (flag==1)
        break;
    end
    Xnew = (Gq+ ((1/h) * Cq)) \ ((((1/h) * Cq) * Xold(:,end))+Bq*u(i));
    Xold=[Xold, Xnew];
    Yold=Lq'*Xold;
    %flag=0;
    for j=1:length(Yold)
        if (length(find(sequence==j))==0)
            if(abs(Yold(j,end)-Yold(j,end-1))<=1e-3)
                sequence=[sequence,j];
%                 if (flag==0)
%                     flag=1;
%                 else
%                     error('Time step too large');
             end
             if (length(sequence)==MatrixSize)
                 flag=1;
             end
         end
     end
end
%     %Create MeshSize x MeshSize Matrix for easier viewing
%     for i=1:MatrixSize
%         OutputXnew(floor((i-1)/MeshSize)+1,rem(i-1,MeshSize)+1) = Xnew(i,1);
%         OutputXdiff(floor((i-1)/MeshSize)+1,rem(i-1,MeshSize)+1) = Xdiff(i,1);
%     end
%end
%Xold=0;
%Xdiff = abs(Lq'*(Xnew - Xold));

%Xnew=L'*Xnew;

% for i=1:sqrt(length(Xnew))
%     for j=1:sqrt(length(Xnew))
%     U(i,j)=Xnew((i-1)*MeshSize+j);
%     end
% end

% [Rx,Ry]=gradient(U);
% Y=abs(Rx)+abs(Ry);
% 
% for i=1:MeshSize
%     for j=1:MeshSize
%         M((i-1)*MeshSize+j)=Y(i,j);
%     end
% end
% 
%[A,I]=sort(M);

I=[];
for i=1:length(sequence)
    b=mod(sequence(i)-1,MeshSize);
    a=(sequence(i)-b-1)/MeshSize;
    I(i,:)=[a+1,b+1];
end
% for i=1:length(x)
%     a=x(i);
%     for j=1:length(y)
%         b=y(j);
%         if ((a>1) && (U(a,b)<U(a-1,b)))
%             continue;
%         end
%         if ((b>1) && (U(a,b)<U(a,b-1)))
%             continue;
%         end
%         if ((a<MeshSize) && (U(a,b)<U(a+1,b)))
%             continue;
%         end
%         if ((b<MeshSize) && (U(a,b)<U(a,b+1)))
%             continue;
%         end
%         I=[I;[a,b]];
%     end
% end

%if (length(sequence)<s^2*N)
%    error('Iteration number too small');
%end
%[Y,I]=sort(abs(Xnew));


% sequence0=[];
% for i=1:length(sequence)
%     a=mod(sequence(i),s*N0);
%     if(mod(a,s)==0)
%         v=a/s;
%         if (a<=s*N0)
%             b=(sequence(i)-a)/s/N0;
%             if (mod(b,s)==s-1)
%             if(a==0)
%                 v=10;
%                 b=b-1;
%             end
%                 sequence0=[sequence0;b,v];
%             end
%         end
%     end
% end
%             
       
%sequence=sequence0;
    
sequence=I;

%toc,

fid=fopen('input.data','w');
fprintf(fid,'%d\n',data(1));
for i=2:2:length(data)
    fprintf(fid, '%d %d\n',data(i),data(i+1));
end
for i=size(sequence,1):-1:1
    %fprintf(fid,'%d %d\n',x0(sequence(i))/s,y0(sequence(i))/s);
    fprintf(fid,'%d %d\n',sequence(i,1),sequence(i,2));
end
fclose(fid);

system('cdctree.exe')
%system('CktSteiner.exe'); 
% for i=1:length(x)
%     for j=1:length(y)
%         U(i,j)=abs(Xnew((i-1)*MeshSize+j));
%     end
% end


%contour(sort(x),sort(y),U,100);

%hold on;

%plot(x,y,'rd');
