clear all, close all
%load /home/camporea/geospace_work/feat_selection.mat
load /home/camporea/geospace_work/iqa/iqa_all_features.mat
T=Table_features;clear Table_features

load /home/camporea/geospace_work/frn/frn_all_features.mat
T=[T; Table_features];clear Table_features

load /home/camporea/geospace_work/ott/ott_all_features.mat
T=[T; Table_features];clear Table_features

%T=very_large_T; clear very_large_T;
%T.Properties.VariableNames(2)={'log10(maxdBdt)'};
%T.Properties.VariableNames(6:22)={'log10(maxBgeo-minBgeo)' 
%  'log10(maxBxgeo-minBxgeo)' 
%  'log10(maxBygeo-minBygeo)' 
%  'log10(maxBmhd-minBmhd)' 
%  'log10(maxBfac-minBfac)'
%  'log10(maxBhal-minBhal)'
%  'log10(maxBped-minBped)'
%  'maxBgeo' 
%  'minBgeo' 
%  'maxBmhd' 
%  'minBmhd' 
%  'maxBfac' 
%  'minBfac' 
%  'maxBhal' 
%  'minBhal' 
%  'maxBped' 
%  'minBped'}';
%  T(:,13:22)=array2table(log10(abs(table2array(T(:,13:22)))));
%  %T(:,13:end)=[];
%  
%  T=T(:,[3:end 2]);

Time_interval = T.Time;

f=find(isnan(T.("log10(maxB-minB)_geo")));T(f,:)=[];
f=find(T.Time>='2019-01-01');
T(f,:)=[];

T.Time=[];
T.log10_maxB_minB_ =[];
T.log10_maxBx_minBx_ =[];
T.log10_maxByminBy_ =[];
T.maxdBdt=log10(T.maxdBdt);
T.maxdBdt_before=log10(T.maxdBdt_before);
T.n=log10(abs(T.n));
%T.T=log10(abs(T.T));
T.Vx=log10(abs(T.Vx));

T.Properties.VariableNames(1:11)={'range(B)','range(Bn)','range(Be)','dB/dt_{geo}','V_x','B_z','n','E','Sym-H','dB/dt_{past}','range(B)_{past}'};


[i j]=find(isnan(table2array(T)));
T(i,:)=[];
[i j]=find(isinf(table2array(T)));
T(i,:)=[];
f=find(T.maxdBdt>4.9);
T(f,:)=[];

s=size(T,1);


N=normalize(T);
%y=table2array(T(:,end));
%[ysort yperc]=sort(y);
%yperc=yperc/max(yperc);
%N(:,end)=table(yperc);

list_features={'all'};


%return
% backward elimination
for i=1:11
c=cvpartition(height(T),'Holdout',0.5);
Mdl=fitlm(N(c.training,:),'RobustOpts','on');
p=Mdl.Coefficients.tStat(2:end);
[dum m]=min(abs(p));
disp([num2str(size(N,2)-1) ' features. R_sq = ' num2str(Mdl.Rsquared.Ordinary) '; RMSE= ' num2str(Mdl.RMSE) ' Eliminate tstat = ' num2str(abs(p(m))) '; ' Mdl.CoefficientNames{m+1} ])
list_features=[list_features;Mdl.CoefficientNames{m+1}];
ind(i) = find(strcmp((T.Properties.VariableNames),Mdl.CoefficientNames{m+1}));
N(:,m)=[];
Rsq(i)=Mdl.Rsquared.Ordinary;
Rmse(i)=Mdl.RMSE;
end

list_features=list_features(end:-1:1)
ind=ind(end:-1:1);Rsq=Rsq(end:-1:1);
l=1:length(Rsq);

plot(l,Rsq,'ko-','linewidth',2), hold on
%text(l+0.1,Rsq+0.1,char(list_features(1:12)),'Fontsize',16)
set(gca,'fontsize',14)
set(gca,'xtick',[0:11])
ticks=[{''};list_features];
set(gca,'XTickLabel',ticks,'XTickLabelRotation',90)
ylabel('R^2','Rotation',0)
axis([0 12.5 0.8 0.83])

return
% Fisher score

station_list={'ott','iqa','frn'};
station_number=0;

for station=station_list

    station_number=station_number+1;
    
    clear T
    filename=(['/home/camporea/geospace_work/' char(station) '/' char(station) '_all_features_new.mat'])
    load(filename)
    T=Table_features;clear Table_features
    f=find(isnan(T.("log10(maxB-minB)_geo")));T(f,:)=[];
    T.Time=[];
    T.log10_maxB_minB_ =[];
    T.log10_maxBx_minBx_ =[];
    T.log10_maxByminBy_ =[];
    T.maxdBdt=log10(T.maxdBdt);
    T.maxdBdt_before=log10(T.maxdBdt_before);
    T.n=log10(abs(T.n));
    T.T=log10(abs(T.T));
    T.Vx=log10(abs(T.Vx));

    y=normalize(T.maxdBdt); thre=quantile(y,0.95);

    t=(y>=thre);

    for i=1:length(ind)
        X=table2array(T(:,ind(i)));
       [~,featureScore(i,station_number)] = fisher_rank(X',t); 
    end
    
end

%N=[];
%for i=1:length(index)

%N=[N T(:,index(i))];
%N=normalize(N);
%Mdl=fitlm(table2array(N),y,'RobustOpts','on');
%disp([num2str(i) ' features (added ' char(T.Properties.VariableNames(index(i))) ') Rsq= ' num2str(Mdl.Rsquared.Ordinary) ' RMSE = ' num2str(Mdl.RMSE)])
%end




