clc
fclose('all');
clear all;
clf;
close all;
N= 1701;

%matfiles = dir(fullfile('C:', 'Users', 'Siavash', 'Downloads', 'RNA_STRAND_data', 'RNA_STRAND_data','all_ct_files2','*.ct'));
matfiles = dir(fullfile(cd,'*.ct'));
[tmp ind] = sort([matfiles.bytes]);
matfiles= matfiles(ind);
%d = zeros(1,N);
%len = zeros(1,N);

%choose the file to be plotted
dist = zeros(N,1);
d1 = zeros(N,1);
d2 = zeros(N,1);
vertices = zeros(N,1);
face = zeros(N,1);

%numVertices = zeros(1,N)+2;

for lm = 1:N
    %extraX = [];
    %extraY = [];
    
    %str = ['C:\Users\Siavash\Downloads\RNA_STRAND_data\RNA_STRAND_data\all_ct_files2\',matfiles(lm).name];
    str = [cd,'/',matfiles(lm).name];
    fid = fopen(str);
    inFile = textscan(fid,'%s');
    idx = find(strcmp(inFile{1,1}(1:min(100,length(inFile{1,1}))), '1'));
    fclose(fid);
    pc = 1;
    while (~(strcmp(inFile{1,1}(idx(pc)+6), '2')) || ~(strcmp(inFile{1,1}(idx(pc)+5), '1')));
        pc = pc + 1;
    end
    idx1 = idx(pc);
    in = reshape(inFile{1,1}(idx1:length(inFile{1,1})),6,(length(inFile{1,1})-idx1+1)/6);
    dist(lm) = -1;
    if (sum(strcmpi(in(2,:),'T')) == 0)
    
    %counting number of possible pairs:
    len(lm) = size(in,2);
    connectNum = in(3:5,:);
    connectType = in(2,:);
    ind_CG = find(strcmpi(connectType,'C'));
    CG_pairs = 0;
    for i = ind_CG
        toBeCompared = str2double(connectNum(:,i));
        cond = find(toBeCompared ~= 0 & toBeCompared ~= i-1 & toBeCompared ~=i+1);
        if ~isempty(cond)
            if strcmpi(connectType(toBeCompared(cond)),'G')
               CG_pairs = CG_pairs + 1;  
            end
        end
    end
    
    ind_AT = find(strcmpi(connectType,'A'));
    AT_pairs = 0;
    for i = ind_AT
        toBeCompared = str2double(connectNum(:,i));
        cond = find(toBeCompared ~= 0 & toBeCompared ~= i-1 & toBeCompared ~=i+1);
        if ~isempty(cond)
            if (strcmpi(connectType(toBeCompared(cond)),'T')) || (strcmpi(connectType(toBeCompared(cond)),'U'))
               AT_pairs = AT_pairs + 1;  
            end
        end
    end
    
    ind_GU = find(strcmpi(connectType,'G'));
    GU_pairs = 0;
    for i = ind_GU
        toBeCompared = str2double(connectNum(:,i));
        cond = find(toBeCompared ~= 0 & toBeCompared ~= i-1 & toBeCompared ~=i+1);
        if ~isempty(cond)
            if (strcmpi(connectType(toBeCompared(cond)),'U')) || (strcmpi(connectType(toBeCompared(cond)),'T'))
               GU_pairs = GU_pairs + 1;  
            end
        end
    end 
    
           

    
        
    %reading the results for analysis
    %str2 = ['C:\Users\Siavash\Downloads\RNA_STRAND_data\Result3D\',matfiles(lm).name,'.txt'];
    str2 = [cd,'/resultGU/',matfiles(lm).name,'.txt'];
    fid2 = fopen(str2,'r');
    inFile2 = textscan(fid2,'%d %d %d');
    fclose(fid2);
    xv = double(inFile2{1,1});
    yv = double(inFile2{1,2});
    zv = double(inFile2{1,3});
    poly3D = [xv, yv, zv];
    poly_sum = sum(poly3D,1);
    dist(lm) = 0;
    vertices(lm) = size(poly3D,1);
    if size(poly3D,1)>2 
        switch (num2str([poly_sum(1) ~= 0, poly_sum(2) ~=0, poly_sum(3) ~=0]))
            case '1  1  1'
                [k, volume] = convhulln(poly3D);
                face(lm) = size(k,1);
                poly = poly3D;
                point = [AT_pairs, CG_pairs, GU_pairs];
                dist(lm) = point2_3Dpoly_dist(k, poly3D, point);
                vec_len = sqrt(AT_pairs^2 + CG_pairs^2 + GU_pairs^2);
                dim = 3;
                d1(lm) = dist(lm)/nthroot(volume,3);
                d2(lm) = dist(lm)/vec_len;
                side_num(lm) = count_3Dpoly_sides(k,poly3D);
            case '1  1  0'
                dim = 2;
                dist(lm) = (-1)*p_poly_dist(AT_pairs, CG_pairs, xv, yv);
               %dist(lm) =0;
                vec_len = sqrt(AT_pairs^2 + CG_pairs^2);
                poly = poly3D(:,1:2);
                area = polyarea(xv,yv);
                d1(lm) = dist(lm)/sqrt(area);
                d2(lm) = dist(lm)/vec_len;
                
            case '1  0  1';
                dim = 2;
                dist(lm) = (-1)*p_poly_dist(AT_pairs, GU_pairs, xv, zv);
                %dist(lm) =0;
                vec_len = sqrt(AT_pairs^2 + GU_pairs^2);
                poly = poly3D(:,[1,3]);
                area = polyarea(xv,zv);
                d1(lm) = dist(lm)/sqrt(area);
                d2(lm) = dist(lm)/vec_len;
                
            case '0  1  1';
                dim = 2;
                dist(lm) = (-1)*p_poly_dist(CG_pairs, GU_pairs, yv, zv);
                 %dist(lm) =0;
                vec_len = sqrt(CG_pairs^2 + GU_pairs^2);
                poly = poly3D(:,2:3);
                area = polyarea(yv,zv);
                d1(lm) = dist(lm)/sqrt(area);
                d2(lm) = dist(lm)/vec_len;

                
        end
        
        
    end
    %if (dist(lm) > 6) 
     %   lm
   %end
   
    %d1(lm) = p_poly_dist(AT_pairs, CG_pairs, xv, yv);
    %d2(lm) = p_poly_dist(AT_pairs, CG_pairs, xv, yv)/sqrt(max(polyarea(xv,yv),1));
    %d3(lm) = p_poly_dist(AT_pairs, CG_pairs, xv, yv)/sqrt(CG_pairs^2+AT_pairs^2);
    
%     XY = [xv,yv];
%     if length(xv)>2
%         dx = diff(xv);
%         dy = diff(yv);
%         for i = 1: length(dx)-1
%             if (dx(i) == dx(i+1)) && (dy(i) == dy(i+1))
%                 extraX = [extraX; xv(i+1)];
%                 extraY = [extraY; yv(i+1)];
%             end
%         end
%         extraXY = [extraX,extraY];
%         
%         if length(extraXY)>0
%             newXY = setxor(XY, extraXY,'rows');
%         else
%             newXY = XY;
%         end
%         
%         numVertices(lm) = size(newXY,1);
%     end
    end
end


%to plot the analysis results for the specified file


%no_inOrOut_points = sum(dist>0)     % number of files for which the point is inside or outside the polygon;
%index_inOrOut = find(dist>0)        %index of the files for which the point is inside or outside the polygon;
% %subplot(2,2,1);
 
%figure
%hist(d1(d1>=0));    %normalized to sguare root of area or 3rd root of volume(whichever applicable)
%title('Normalized Distance');
%xlabel('Distance/3rd Root(Volume)');
%ylabel('Number of Sequences')
 

% figure
% plot(len,vertices,'.');

 % title('Without Normalization');
% %subplot(2,2,2);
 

%figure
%hist(dist(dist>=0)); % not normalized
%title('Distance');
%xlabel('Distance');
%ylabel('Number of Sequences')


%figure
%hist(d2(d2>=0)); % not normalized
%title('Normalized Distance');
%xlabel('Distance/Feature Vector Length');
%ylabel('Number of Sequences')

hist(side_num(side_num>0));
%to plot the Shapes for some files we use plot_analysis_results
%function.(refer to its help). As an Example:
%plot_analysis_results([5,64]);



% title('Normalized to Square Root of Polygon Area');
% %subplot(2,2,3);
% %hist(abs(d3));
% %title('Normalized to Feature Vector Length');
% 
% %subplot(2,2,4);
% %plot(newXY(:,1),newXY(:,2),'r*');
% %title('Polygone')
% length(find(d1>0))
% length(find(d1<0))
% length(find(d1 == 0))
% figure
% plot(len,abs(d1),'.');
% title('Distance vs. Sequence Length');
% figure
% plot(len,numVertices,'.');
% title('Number of Polygon Vertices vs. Strand Length');
