clear all
set(0, 'DefaultTextInterpreter', 'none')


% FolderName = ''; 


% % % All togheter
% Domain_ResNos = {1:348};
% Domain_names = {'all'};

% % Community
Domain_ResNos = { [1 35:52  56:73  107:110  117:123  319:322  328:334]; [53:55  74:85  88  111:116  335:348]; [16:34  86:87  89:101  103:106  150:164  185:190] ; [124:138  170:184  314:318  323:327];  [14  102  139:149  298:313]; [165:169  191:228  254  269  271:283]; [229:237  255:260]; [238:253] ; [261:268  270 284:297]};
Domain_names = {'A' 'B' 'C' 'D' 'E' 'F' 'F1' 'G' 'H'};

% DOMAIN - Hunter nomunclature
% Domain_ResNos = {1:42;  43:64; 65:83;  84:98;  99:113; 114:138; 139:162; 163:178; 179:193; 194:211; 212:240; 241:260; 261:297; 298:349};
% Domain_names = {'n-t'   'I'    'II'   'III'    'IV'     'V'     'VIA'    'VIB'    'VII'    'VIII'    'IX'     'X'      'XI'     'XII'};

ReferenceResidue = 311; %%     %     '+ resno_offset' 
ReferenceState	 = 1;   

PC2_Filter  	=  100;   %% use lower number ~ 0.2 to pick only near linear points trajectories
 
%%
Assignment_Column = 1; w1_Column = 2; w2_Column = 3;    
Calculate_ChemicalShiftScaling = 0; 

resno_offset_for_NegativeResidues = 0;

fnames = dir([FolderName '\*.list']);

for FileNo = 1:length(fnames) 
    FileName        =  fnames((FileNo)).name;
    TotalColumns    = GetNum_Clmns([FolderName '\' FileName]);
    fid             = fopen([FolderName '\' FileName],'rt');

    ListData          = textscan(fid, ['           %s' repmat('    %f', 1, TotalColumns-1)] , 'HeaderLines',2);
    fclose(fid);  
    s(FileNo).FileName = FileName;
    try
        [OrderNo, Field_MHz]    = Get_OrderNo_Field_FromName_1(FileName);
        s(FileNo).OrderNo       = OrderNo; s(FileNo).Field_MHz     = Field_MHz; 
        s(FileNo).ListData      = ListData;
    end
    
    temp_resNos_Array = [];
    for k = 1:length(ListData{1}) 
        temp_resNos_Array  = [temp_resNos_Array ExtractNumFromString_1(ListData{Assignment_Column}{k})];
    end
    resNos_Array{FileNo} = temp_resNos_Array;
end

 
    
    
%%
if min([resNos_Array{:}]) > 1
    ResNo2Idx = @(resno) resno   ;
    ResIdx2No = @(ResIdx)  ResIdx ;
else
    ResNo2Idx = @(resno) resno - min([resNos_Array{:}]) + 1;
    ResIdx2No = @(ResIdx)  ResIdx  + min([resNos_Array{:}]) - 1;
end


a = zeros(length(s) ,  ResNo2Idx(max([resNos_Array{:}])) -  ResNo2Idx(min([resNos_Array{:}])) );
for FileNo = 1:length(s) 
    a(FileNo, ResNo2Idx(   resNos_Array{FileNo} ))     = 1;
end
Idx2Consider = find(prod(a, 1)==1);


%%

for FileNo = 1:length(s) 
    ResIdx_Array = [];
    for k = 1:length(s(FileNo).ListData{1}) 
        resno = ExtractNumFromString_1(s(FileNo).ListData{Assignment_Column}{k});  
        resno = resno(1) ;
        ResIdx = ResNo2Idx(resno);
        
        
        if ~isempty( find(Idx2Consider == ResIdx, 1))
            s(FileNo).IsPresent(ResIdx) = 1;

    %         *******************************
            if isempty(find( ResIdx_Array == ResIdx, 1 ))
                s(FileNo).resno(ResIdx) = resno;
                s(FileNo).resname{ResIdx} = s(FileNo).ListData{Assignment_Column}{k}; 
                s(FileNo).x(ResIdx) = s(FileNo).ListData{w2_Column}(k);   s(FileNo).y(ResIdx) = s(FileNo).ListData{w1_Column}(k);      

                ResIdx_Array = [ResIdx_Array ResIdx];
            else
                fprintf('repeated res no : %s :::   %d\n', FileName, resno);
                s_(FileNo).x(ResIdx) = s(FileNo).ListData{w2_Column}(k);   s_(FileNo).y(ResIdx) = s(FileNo).ListData{w1_Column}(k);      
                s_(FileNo).resname{ResIdx} = s(FileNo).ListData{Assignment_Column}{k};
            end        
        end
    end
end

%% REFERENCING for residue
ReferenceResIdx = ResNo2Idx(ReferenceResidue);
for k = 1:length(s)  
    x_off = (mean(s(ReferenceState).x(ReferenceResIdx)) - mean(s(k).x(ReferenceResIdx)));
    y_off = (mean(s(ReferenceState).y(ReferenceResIdx)) - mean(s(k).y(ReferenceResIdx)));
    for ki = 1:size(s(k).x, 2) 
        if s(k).x(ki)
            s(k).x(ki) =  x_off  + s(k).x(ki) ; 
            s(k).y(ki) =  y_off  + s(k).y(ki); 
        end
    end    
end 

%%


%%
s1 = s;
clear X1 Y1
for k = 1:length(s)
    X1(k,:) =  s1(k).x;
    Y1(k,:) =  s1(k).y;
end
% if assignmes are not present in all states remove it 
temp_idx = find(~isnan(sum(X1)));
X1 =  X1(:, temp_idx);
Y1 =  Y1(:, temp_idx);


for k_d = 1:length(Domain_ResNos)
    resNumbers = Domain_ResNos{k_d};
    Idx =  ResNo2Idx(resNumbers); 
     
    X1_{k_d} = X1(:,Idx(  s(1).IsPresent(Idx) == 1));
    Y1_{k_d} = Y1(:,Idx(  s(1).IsPresent(Idx) == 1));
end
%%  
figure
YShift4Plot_Array = 0;
YShift4Plot  = 0;
counter = 1;
for k_d = 1:length(Domain_ResNos)
    

    X1 = X1_{k_d}; 
    Y1 = Y1_{k_d}; 

    if size(X1, 2) >= 1


        clear PC1_all  PC2_all
%         figure;    
        hold on
        for kr  =  1:size(X1,2) 
            XY = [X1(:,kr) Y1(:,kr)];
            XY = XY - mean(XY);

            W = pca(XY);
            PC = XY*W; 
            PC = XY*W/std(PC(:,1)); 

            if (PC(1,1) - PC(end,1)) > 0
                PC = (-PC);
            end
            PC1 = PC(:,1); PC2 = PC(:,2);
            if (max(PC2) - min(PC2)) < PC2_Filter
                PC1_all(:,kr) = PC1;PC2_all(:,kr) = PC2;
%                     plot(PC1,PC2,'*');
            end        
        end

        %
        clr1 =  getColor1_21(length(s1));
        xdtbn = linspace(-3,3,200);
%         figure(1)
        for k =  1:size(PC1_all,1)
            st      = PC1_all(k,:); st = st(st~=0);   
%             figure(1)
%             plot(xdtbn, normpdf(xdtbn,mean(st),std(st)), 'linewidth', 2, 'color', clr1(k,:)); hold on 
%             xlabel('PC1');ylabel('Probability Density');
%             ConciseToPrint(k,:) = normpdf(xdtbn,mean(st),std(st));
%             figure(2)

    % comment this line if you dont want the 'dots'
%             plot3(PC1_all(k,:),  YShift4Plot + PC2_all(k,:), 8*ones(size(PC1_all(k,:))), '.', 'MarkerSize', 10, 'color', clr1(k,:)) ; hold on
        end
%         clc
%         disp([xdtbn' ConciseToPrint'])



%         figure(2)
        N1 = 201; N2 = 201;
        Yarray = ones(size(PC1_all,1), N1, N2);
        clr1 =  getColor1_21(length(s1));
        xdtbn = linspace(-3,3,200);
        x0_4_LinearReg = []; y0_4_LinearReg = [];
        for k =  1:size(PC1_all,1)

            x = PC1_all(k,:);%/std( PC1_all(k,:));
            y = PC2_all(k,:);%/std( PC2_all(k,:));

            a = std(x); % horizontal radius
            b = std(y); % vertical radius
            ex0 = mean(x); % x0,y0 ellipse centre coordinates
            ey0 = mean(y);
            t=-pi:0.01:pi;

            mu = [ex0 ey0];
            sigma = [a 0; 0 b];
            x1 = linspace(-2, 2, N1);
            x2 = linspace(-2, 2, N2);
            [X1,X2] = meshgrid(x1,x2);
            X = [X1(:) X2(:)];

            Y = mvnpdf(X,mu,sigma);
            Y = reshape(Y,length(x2),length(x1));


            Adata = (Y/max(Y(:))).^8;
            % Adata(Adata < 0.1) = 0;
 %         mesh(x1,  YShift4Plot + x2, 0.5*k+Y.^2, 'AlphaData',  Adata,    'FaceAlpha','flat',   'FaceColor',clr1(k,:), 'linestyle', 'none'); hold on 

%             plot3(ex0, YShift4Plot + ey0, 10, '.', 'markersize', 35, 'color',clr1(k,:) )
            plot3(ex0, YShift4Plot + ey0, 10, 'o', 'markeredgecolor', 'k', 'linewidth', 1, 'markersize', 5, 'MarkerFaceColor',  clr1(k,:)  )


            ex=ex0+a*cos(t);
            ey=ey0+b*sin(t);
            plot3(ex, YShift4Plot + ey, 9*ones(size(ey)), 'LINEWIDTH', 1, 'color',clr1(k,:)); hold on 

%              ylim(0.5*[-1 1]);
            Yarray(k, :, :) =   Y;
            x0_4_LinearReg = [x0_4_LinearReg  ex0];
            y0_4_LinearReg = [y0_4_LinearReg  ey0]; 
        end
        
        plot3([min(x0_4_LinearReg) max(x0_4_LinearReg) ]+0.1*[-1 1], ... 
            YShift4Plot *[1 1], 11*[1 1], 'k');
        
%         b1 = x0_4_LinearReg'\y0_4_LinearReg';
%         yCalc1 = b1*x0_4_LinearReg;  
%         plot3(x0_4_LinearReg, YShift4Plot + yCalc1, 7*ones(size(x0_4_LinearReg)), 'k');
        
        YShift4Plot = YShift4Plot + 0.75;
        YShift4Plot_Array = [YShift4Plot_Array YShift4Plot];
        DomainName_Plot{counter} = Domain_names{k_d};
        counter = counter +1;
    end
end
 
 
ylim([min(YShift4Plot_Array)-0.75 max(YShift4Plot_Array)])
view(2)
xlim([-1.5 1.5]);

yticks(YShift4Plot_Array(1:(end-1)));
yticklabels(DomainName_Plot);
xlabel('PC1')
ylabel('Domains')
[~, Title_name, ~]=fileparts(FolderName);
title(Title_name)
set(gcf,'units', 'inches')
set(gcf, 'position', [0 0 4 0.6*length(YShift4Plot_Array)]);
box on
grid on

%% 



function clr1 = getColor1_21(n)
    if n < 21
        Clr = ...
        [        0.00   0.00    1.00
                 0.9843137255	0.8117647059	0.01568627451
                                0.8745098039	0.1764705882	0.568627451];  
        clr1 = Clr([1:(n-1) end], :);


    else 
        disp('getColor_21  has only 21 colors '); 
    end
end
function numCols = GetNum_Clmns(FileName)
    fid = fopen(FileName,'rt');
    fgets(fid); fgets(fid); tLines = fgets(fid); 
    idx_temp = find(tLines==' ');numCols = length(find(diff(idx_temp) > 1)) ; 
    fclose(fid);
end

    
function [OrderNo, Field_MHz] = Get_OrderNo_Field_FromName_1(Name_string)

%     a ='12_850_APO_WT.list';   output is 12 and 850
    pos_        = find(Name_string =='_');
    pos_        = pos_-1;
    OrderNo     = str2num( Name_string(1:pos_(1)) );
    Field_MHz   = str2num( Name_string((pos_(1)+2):pos_(2)) );
end


function n = ExtractNumFromString_1(str)
 
    tn=  regexp(str,'\s+\.?\d*|-\d+\.?\d*|\.?\d+|-\.?\d+','match'); 
 
    for k = 1:length(tn)
        tn_(k)  = str2double(tn{k});
    end

    n = tn_;
    
end