%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%   read out cyana restraints   %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  eNOES are set to 8 Angstroms if any proton involved is a aromatic  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sigmaData = cyanaUplScaling(sigmaData);

user.toleranceUplUniue2     = 1;
user.toleranceLolUniue2     = 1;
user.toleranceUplUniue1     = 1; %1.20;
user.toleranceLolUniue1     = 1; %0.80;

%methyl are allways applied also for bi-directional eNOEs!!!
%8.5% communication Beat V.
user.toleranceUplMeth       = 1.085;
user.toleranceLolMeth       = 0.915;

%methylene pseudo are allways applied also for bi-directional eNOEs!!!
%5.0% communicaton Beat V.
user.toleranceUplMethylenePseudo       = 1.05;
user.toleranceLolMethylenePseudo       = 0.95;

cyanaOutput                 = [];
cyanaUpl                    = [];
cyanaLol                    = [];

clearvars cyanaOutput cyanaUpl cyanaLol

user.chiNcutOff             = 0;

user.chiNcutOff             = input('Specify a chiN cut off value: ');
% % % user.bi_eNOE                = input('Output eNOEs derrived from two cross peaks only (2) or also if only one cross peak was computed (1)? Choose (2/1): ');
user.bi_eNOE                = 1;
user.fixedNOES              = input('Remove all fixed distances (intra methelene, rings etc) (1) or keep all eNOEs (2)? Choose (1/2): ');



count.cyana                 = 1;

for j=find(ismember([sigmaData.unique], [1 2]))
        
        if sigmaData(j).chiN <= user.chiNcutOff
            
            if sigmaData(j).unique == 1 && abs(sigmaData(j).buildups.rawIzero) ~= 999
                
                cyanaOutput(count.cyana).numcode            = sigmaData(j).numcode;
                cyanaOutput(count.cyana).ASSi_resName       = sigmaData(j).ASSi_resName;
                cyanaOutput(count.cyana).ATOMi              = sigmaData(j).ATOMi;
                cyanaOutput(count.cyana).ASSj_resName       = sigmaData(j).ASSj_resName;
                cyanaOutput(count.cyana).ATOMj              = sigmaData(j).ATOMj;
                
                cyanaOutput(count.cyana).rCyana             = sigmaData(j).rEffij     * (1/sigmaData(j).numcode.cyanaUplScaling)^(1/6);
                cyanaOutput(count.cyana).unique             = 1;
                
                count.cyana                                 = count.cyana + 1;
                
            
            elseif sigmaData(j).unique == 2
                
                
                if sigmaData(sigmaData(j).indxji).chiN <= user.chiNcutOff && abs(sigmaData(j).buildups.rawIzero) ~= 999 && abs(sigmaData(sigmaData(j).indxji).buildups.rawIzero) ~= 999
                    
                    cyanaOutput(count.cyana).numcode            = sigmaData(j).numcode;
                    cyanaOutput(count.cyana).ASSi_resName       = sigmaData(j).ASSi_resName;
                    cyanaOutput(count.cyana).ATOMi              = sigmaData(j).ATOMi;
                    cyanaOutput(count.cyana).ASSj_resName       = sigmaData(j).ASSj_resName;
                    cyanaOutput(count.cyana).ATOMj              = sigmaData(j).ATOMj;
                    
                    cyanaOutput(count.cyana).rCyana             = sigmaData(j).avgrEff    * (1/sigmaData(j).numcode.cyanaUplScaling)^(1/6);
                    cyanaOutput(count.cyana).unique             = 2;
                    
                    count.cyana                                 = count.cyana + 1;
                    
                end
                
                if abs(sigmaData(j).buildups.rawIzero) ~= 999 && (abs(sigmaData(sigmaData(j).indxji).buildups.rawIzero) == 999 || sigmaData(sigmaData(j).indxji).chiN > user.chiNcutOff)
                
                    cyanaOutput(count.cyana).numcode            = sigmaData(j).numcode;
                    cyanaOutput(count.cyana).ASSi_resName       = sigmaData(j).ASSi_resName;
                    cyanaOutput(count.cyana).ATOMi              = sigmaData(j).ATOMi;
                    cyanaOutput(count.cyana).ASSj_resName       = sigmaData(j).ASSj_resName;
                    cyanaOutput(count.cyana).ATOMj              = sigmaData(j).ATOMj;
                    
                    cyanaOutput(count.cyana).rCyana             = sigmaData(j).rEffij     * (1/sigmaData(j).numcode.cyanaUplScaling)^(1/6);
                    cyanaOutput(count.cyana).unique             = 1;
                    
                    count.cyana                                 = count.cyana + 1;
                
                end
            end
            
        elseif sigmaData(j).chiN > user.chiNcutOff             
            
            if sigmaData(j).unique == 2
            
                if sigmaData(sigmaData(j).indxji).chiN <= user.chiNcutOff && abs(sigmaData(sigmaData(j).indxji).buildups.rawIzero) ~= 999
            
                cyanaOutput(count.cyana).numcode                = sigmaData(sigmaData(j).indxji).numcode;
                cyanaOutput(count.cyana).ASSi_resName           = sigmaData(sigmaData(j).indxji).ASSi_resName;
                cyanaOutput(count.cyana).ATOMi                  = sigmaData(sigmaData(j).indxji).ATOMi; 
                cyanaOutput(count.cyana).ASSj_resName           = sigmaData(sigmaData(j).indxji).ASSj_resName;
                cyanaOutput(count.cyana).ATOMj                  = sigmaData(sigmaData(j).indxji).ATOMj;
                
                
                cyanaOutput(count.cyana).rCyana                 = sigmaData(sigmaData(j).indxji).rEffij * (1/sigmaData(sigmaData(j).indxji).numcode.cyanaUplScaling)^(1/6);
                
                cyanaOutput(count.cyana).unique                 = 1; % since it is not an avgrEff that is used but only one of them
            
                count.cyana                                     = count.cyana + 1;
                
                end
            end
        end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% scrub

if user.fixedNOES == 1
    
    fprintf('\n\n%s\r\n','REMOVING FIXED DISTANCES, INCLUDING PROLINE INTRA')
    
    % % % scrub intra methylene etc, i.e. HB2, HB3 in same res
    [cyanaOutput, cyanaOutputScrub] = scrubIntraResMethyleneNew(cyanaOutput);
    
    numScrubed = length(cyanaOutputScrub);
    
    fprintf('\n\n%d%s\r\n',numScrubed,' eNOES were scrubed from the final output')
    
else
    
    fprintf('\n\n%s\r\n','KEEPING ALL eNOEs')
    
end

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


%enter sort function here

count.cyana = 1;

for j = 1:length(cyanaOutput)
    
    if user.bi_eNOE == 2 && cyanaOutput(j).unique == 2
        
        cyanaUpl(count.cyana) = cyanaOutput(j);
        cyanaLol(count.cyana) = cyanaOutput(j);
            
    	cyanaUpl(count.cyana).rCyana  = cyanaUpl(count.cyana).rCyana * user.toleranceUplUniue2;
        cyanaLol(count.cyana).rCyana  = cyanaLol(count.cyana).rCyana * user.toleranceLolUniue2;
        
        count.cyana = count.cyana + 1;
        
    elseif user.bi_eNOE == 1
        
        cyanaUpl(count.cyana) = cyanaOutput(j);
        cyanaLol(count.cyana) = cyanaOutput(j);
        
            if cyanaOutput(j).unique == 2 
        
            cyanaUpl(count.cyana).rCyana  = cyanaUpl(count.cyana).rCyana * user.toleranceUplUniue2;
            cyanaLol(count.cyana).rCyana  = cyanaLol(count.cyana).rCyana * user.toleranceLolUniue2;
                    
            elseif cyanaOutput(j).unique == 1
        
            cyanaUpl(j).rCyana  = cyanaUpl(count.cyana).rCyana * user.toleranceUplUniue1;
            cyanaLol(j).rCyana  = cyanaLol(count.cyana).rCyana * user.toleranceLolUniue1;
            end
        
        count.cyana = count.cyana + 1;
                    
    end
    
end

count.cyana =  count.cyana - 1;


%methyl are allways applied also for bi-directional eNOEs!!!

for j = 1:length(cyanaUpl)
    
        cyanaUpl(j).bothAtoms = [cyanaUpl(j).numcode.ATOMi, cyanaUpl(j).numcode.ATOMj];
                
    if  size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.methylPseudo(1:8))), 2) == 1
        
        cyanaUpl(j).rCyana = cyanaUpl(j).rCyana * user.toleranceUplMeth;
    
    elseif  size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.methylPseudo(9:10))), 2) == 1
        
        cyanaUpl(j).rCyana = cyanaUpl(j).rCyana * (user.toleranceUplMeth *  user.toleranceUplMethylenePseudo);
        
    elseif size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.methylPseudo(1:8))), 2) == 2
        
        cyanaUpl(j).rCyana = cyanaUpl(j).rCyana * user.toleranceUplMeth^2;
        
    elseif size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.methylPseudo(9:10))), 2) == 2
        
        cyanaUpl(j).rCyana = cyanaUpl(j).rCyana * (user.toleranceUplMeth *  user.toleranceUplMethylenePseudo)^2;
        
    elseif size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.methylPseudo)), 2) == 0
        cyanaUpl(j).rCyana = cyanaUpl(j).rCyana;
        
    end
end

for j = 1:length(cyanaLol)
    
        cyanaLol(j).bothAtoms = [cyanaLol(j).numcode.ATOMi, cyanaLol(j).numcode.ATOMj];
                
    if  size(find(ismember(cyanaLol(j).bothAtoms, conditions.all.methylPseudo(1:8))), 2) == 1
        
        cyanaLol(j).rCyana = cyanaLol(j).rCyana * user.toleranceLolMeth;
        
    elseif  size(find(ismember(cyanaLol(j).bothAtoms, conditions.all.methylPseudo(9:10))), 2) == 1
        
        cyanaLol(j).rCyana = cyanaLol(j).rCyana * (user.toleranceLolMeth * user.toleranceLolMethylenePseudo);
        
    elseif size(find(ismember(cyanaLol(j).bothAtoms, conditions.all.methylPseudo(1:8))), 2) == 2
        
        cyanaLol(j).rCyana = cyanaLol(j).rCyana * user.toleranceLolMeth^2;
        
    elseif size(find(ismember(cyanaLol(j).bothAtoms, conditions.all.methylPseudo(9:10))), 2) == 2
        
        cyanaLol(j).rCyana = cyanaLol(j).rCyana * (user.toleranceLolMeth * user.toleranceLolMethylenePseudo)^2;
        
    elseif size(find(ismember(cyanaLol(j).bothAtoms, conditions.all.methylPseudo)), 2) == 0
        
        cyanaLol(j).rCyana = cyanaLol(j).rCyana;
        
    end
end

% methylene Pseudo corrections are allways applied also for bi-directional eNOEs!!!

for j = 1:length(cyanaUpl)
    
        cyanaUpl(j).bothAtoms = [cyanaUpl(j).numcode.ATOMi, cyanaUpl(j).numcode.ATOMj];
                
    if  size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.methylenePseudo)), 2) == 1
        
        cyanaUpl(j).rCyana = cyanaUpl(j).rCyana * user.toleranceUplMethylenePseudo;
        
    elseif size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.methylenePseudo)), 2) == 2
        
        cyanaUpl(j).rCyana = cyanaUpl(j).rCyana * user.toleranceUplMethylenePseudo^2;
        
    elseif size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.methylenePseudo)), 2) == 0
        cyanaUpl(j).rCyana = cyanaUpl(j).rCyana;
        
    end
end

for j = 1:length(cyanaLol)
    
        cyanaLol(j).bothAtoms = [cyanaLol(j).numcode.ATOMi, cyanaLol(j).numcode.ATOMj];
                
    if  size(find(ismember(cyanaLol(j).bothAtoms, conditions.all.methylenePseudo)), 2) == 1
        
        cyanaLol(j).rCyana = cyanaLol(j).rCyana * user.toleranceLolMethylenePseudo;
        
    elseif size(find(ismember(cyanaLol(j).bothAtoms, conditions.all.methylenePseudo)), 2) == 2
        
        cyanaLol(j).rCyana = cyanaLol(j).rCyana * user.toleranceLolMethylenePseudo^2;
        
    elseif size(find(ismember(cyanaLol(j).bothAtoms, conditions.all.methylenePseudo)), 2) == 0
        
        cyanaLol(j).rCyana = cyanaLol(j).rCyana;
        
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  eNOES are set to 8 Angstroms if any proton involved is a aromatic  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% conditions.all.aromatics

for j = 1:length(cyanaUpl)
    
        cyanaUpl(j).bothAtoms = [cyanaUpl(j).numcode.ATOMi, cyanaUpl(j).numcode.ATOMj];
                
    if  size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.aromatics)), 2) == 1 ...
        || size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.aromatics)), 2) == 2    
        
        cyanaUpl(j).rCyana = 8;
        
    elseif size(find(ismember(cyanaUpl(j).bothAtoms, conditions.all.aromatics)), 2) == 0
        
        cyanaUpl(j).rCyana = cyanaUpl(j).rCyana;
        
    end
end

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

% % % % % % % % scrub risk groups before writing to file!!
% % % % % % % 
% % % % % % % if ~isempty(user.riskGROUP) 
% % % % % % % 
% % % % % % %     fprintf('\n\n%s\r\n','SCRUBING THE RISK GROUPS')
% % % % % % %     user.riskGROUP
% % % % % % %     
% % % % % % %     
% % % % % % %     for j = 1:length(cyanaUpl)
% % % % % % %             
% % % % % % %         cyanaUpl(j).bothASS = [cyanaUpl(j).numcode.ASSi, cyanaUpl(j).numcode.ASSj];
% % % % % % %  
% % % % % % %         
% % % % % % %     end
% % % % % % %     
% % % % % % %     count = 1;
% % % % % % %     numRiskGroupScrubed = 0;
% % % % % % %         
% % % % % % %     for j = 1:length(cyanaUpl)
% % % % % % %     
% % % % % % %         if size(find(ismember(cyanaUpl(j).bothASS, user.riskGROUP)), 2) == 0
% % % % % % %         
% % % % % % %             cyanaUpl(count) = cyanaUpl(j);
% % % % % % %             count =  count + 1;
% % % % % % %         
% % % % % % %         else
% % % % % % %             
% % % % % % %                     numRiskGroupScrubed = numRiskGroupScrubed + 1;
% % % % % % %         end
% % % % % % %     
% % % % % % %     end
% % % % % % %  	fprintf('\n\n%d%s\r\n',numRiskGroupScrubed,' NOES were in risk group and scrubed')
% % % % % % % end

%%%%%%%%%%%%%%%%%%%%%% upl and lol write out %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if user.bi_eNOE == 1
    fid = fopen(sprintf('./%s/%s_all_genNOEs.upl', work.finalCollected, user.projectname),'wt');
elseif user.bi_eNOE == 2
    fid = fopen(sprintf('./%s/%s_bi_eNOE_genNOEs.upl', work.finalCollected, user.projectname),'wt');
end
    
for j = 1:length(cyanaUpl)
    
            fprintf(fid, '%4d%5s%5s%8d%5s%5s%16.2f\r\n', cyanaUpl(j).numcode.ASSi...
                , cyanaUpl(j).ASSi_resName,  cyanaUpl(j).ATOMi...
                , cyanaUpl(j).numcode.ASSj,  cyanaUpl(j).ASSj_resName...
                , cyanaUpl(j).ATOMj,         cyanaUpl(j).rCyana);
end

fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%% chimera eNOE network %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if user.bi_eNOE == 1
    fid = fopen(sprintf('./%s/%s_chimera_all_genNOEs_eNOEnetwork', work.finalCollected, user.projectname),'wt');
elseif user.bi_eNOE == 2
    fid = fopen(sprintf('./%s/%s_chimera_bi_genNOEs_eNOE_eNOEnetwork', work.finalCollected, user.projectname),'wt');
end

fprintf('\n\n%s\r\n','TO DISPLAY THE eNOE NETOWORK in CHIMERA:')
fprintf('%s\r\n','Edit a copy of your pdb to contain only one confomer')
fprintf('%s\r\n\n','Start chimera; load structure; Actions > Atoms/Bonds > Show then Tools > Depiction > pseudoBond Reader')

for j = 1:length(cyanaUpl)
            
    if user.bi_eNOE == 1

        fprintf(fid, ':%d@%s :%d@%s yellow \n', cyanaUpl(j).numcode.ASSi...
                , cyanaUpl(j).ATOMi, cyanaUpl(j).numcode.ASSj, cyanaUpl(j).ATOMj);
        
    elseif user.bi_eNOE == 2

        fprintf(fid, ':%d@%s :%d@%s green \n', cyanaUpl(j).numcode.ASSi...
                , cyanaUpl(j).ATOMi, cyanaUpl(j).numcode.ASSj, cyanaUpl(j).ATOMj);
        
    end
    
end

fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%