Test_Sergio_Acop_mvnrnd.asv

clear all
close all
clc

setenv('PATH', [getenv('PATH') ';C:\Program Files\R\R-3.3.2\bin\']);

%%

electrode_selected_number=[1:5:20];
%%

for t=electrode_selected_number
    
    
    Mu=[zeros(1000,t)];
    SIGMA=diag([ones(t,1)],0);
    Rb=mvnrnd(Mu,SIGMA);
    
    Mu=[ones(1000,t)];
    Mu=Mu*5;
    SIGMA=diag([ones(t,1)],0);
%     
%     SIGMA(1,2)=0;
%     SIGMA(2,1)=0;
    
    Rf=mvnrnd(Mu,SIGMA);
    
    
    Bind=Rb;
    Feat=Rf;
        
    Complete=[Bind' Feat'];
    
    fit_size=2*t+sum([1:1:t])+1;
    
    csvwrite(['CSV/Binding.csv'],Bind);
    csvwrite(['CSV/Features.csv'],Feat);
    Complete=[Bind' Feat'];
    csvwrite(['CSV/Complete.csv'],Complete');
    
    csvwrite(['CSV/FitSize.csv'],fit_size);
    
    %%
    %!unset DYLD_LIBRARY_PATH; Rscript  skewNromalFitBind.R
    
    !Rscript skewNromalFitBind.R
    snParamBind = csvread('skewNromalFitedDataBind.csv');
    % snParam(snParam==-999999)=NaN;
    !Rscript skewNromalFitFeat.R
    snParamFeat = csvread('skewNromalFitedDataFeat.csv');
    % snParamBind(snParam==-999999)=NaN;
    !Rscript  skewNromalFit.R
    snParam = csvread('skewNromalFitedData.csv');
    % snParamFeat(snParamGO==-999999)=NaN;
    
    %%

    size(snParam,1)
    
    for n=1:size(snParam,1)
        
        alfainit=t+sum([1:1:t])+1;
        alfaend=alfainit+t-1;
        
        alfa=snParam(n,alfainit:alfaend);
        
        likelihood=snParam(end);
        
        init=electrode_selected_number+1;
        for i=1:electrode_selected_number %for 10 channels
            
            l=electrode_selected_number-i ;% #components
            
            omega(i,i:electrode_selected_number)=snParam(1,init:(init+l));
            omega(i:electrode_selected_number,i)=snParam(1,init:(init+l));
            
            init=(init+l+1);
            
        end
        
        M=100000;
        a=randn(100000,1);
        b=randn(100000,1);
        W(find(sqrt(alfa*alfa')*a>b))=a(find(sqrt(alfa*alfa')*a>b));
        W(find(sqrt(alfa*alfa')*a<=b))=-a(find(sqrt(alfa*alfa')*a<=b));
        H(n) = 1/2*log((det(omega))) + 1 + log(2*pi) - mean(log(2*normcdf(sqrt(alfa*alfa')*W)));
        
        NoNorm(t)=H(n);
        
        sknormBind(t)= mean(2*log(normcdf(sqrt(alfa*alfa')*W)));
        
        %BIND
        
        alfaBind=snParamBind(n,alfainit:alfaend);
        
        
        likelihoodbind=snParamBind(end);
        
        init=electrode_selected_number+1;
        for i=1:electrode_selected_number
            
            l=electrode_selected_number-i ;% #components
            
            omegaBind(i,i:electrode_selected_number)=snParamBind(1,init:(init+l));
            omegaBind(i:electrode_selected_number,i)=snParamBind(1,init:(init+l));
            
            init=(init+l+1);
            
        end
        
        M=100000;
        a=randn(100000,1);
        b=randn(100000,1);
        W(find(sqrt(alfaBind*alfaBind')*a>b))=a(find(sqrt(alfaBind*alfaBind')*a>b));
        W(find(sqrt(alfaBind*alfaBind')*a<=b))=-a(find(sqrt(alfaBind*alfaBind')*a<=b));
        HBind(n) = 1/2*log((det(omegaBind))) + 1 + log(2*pi) - mean(log(2*normcdf(sqrt(alfaBind*alfaBind')*W)));
        
        BindNoNorm(t)=HBind(n);
        
        sknormBind(t)= mean(2*log(normcdf(sqrt(alfaBind*alfaBind')*W)));
        
        %FEAT
        
        alfaFeat=snParamFeat(n,alfainit:alfaend);
        
        likelihoodFeat=snParamFeat(end);
        init=electrode_selected_number+1;
        for i=1:electrode_selected_number %for 10 channels
            
            l=electrode_selected_number-i ;% #components
            
            omegaFeat(i,i:electrode_selected_number)=snParamFeat(1,init:(init+l));
            omegaFeat(i:electrode_selected_number,i)=snParamFeat(1,init:(init+l));
            
            init=(init+l+1);
            
        end
        
        M=100000;
        a=randn(100000,1);
        b=randn(100000,1);
        W(find(sqrt(alfaFeat*alfaFeat')*a>b))=a(find(sqrt(alfaFeat*alfaFeat')*a>b));
        W(find(sqrt(alfaFeat*alfaFeat')*a<=b))=-a(find(sqrt(alfaFeat*alfaFeat')*a<=b));
        HFeat(n) = 1/2*log((det(omegaFeat))) + 1 + log(2*pi) - mean(log(2*normcdf(sqrt(alfaFeat*alfaFeat')*W)));
        
        FeatNoNorm(t)=HFeat(n);
        sknormFeat(t)= mean(log(2*normcdf(sqrt(alfaFeat*alfaFeat')*W)));
    end
    
    MI(t)=H - 1/2*HFeat - 1/2*HBind;
    
    HH(t)=H;
    HHf(t)=HFeat;
    HHb(t)=HBind;



end

%%

Miprom=mean(MI);
HHprom=mean(HH);
HHfprom=mean(HHf);
HHbprom=mean(HHb);