% Generate Data

clear all;
warning('off','MATLAB:dispatcher:InexactCaseMatch')
warning('off','MATLAB:dispatcher:nameConflict')
sample=200; dim=10; imp_dim=2; 
%coef_range=8; cov_range=10;
coef_range=5; cov_range=2;
ranks=zeros(n_sim,dim);
ranks_g=zeros(n_sim,dim);
ranks_l=zeros(n_sim,dim);
coef=unifrnd(coef_range-2,coef_range,imp_dim,1);

%coef=[0.39,0.39];

n_sim=50;
rates_reg=zeros(n_sim,4);

for j=1:n_sim;
    corrset=1:imp_dim;    
%     coef=unifrnd(coef_range-2,coef_range,imp_dim,1);
    
    loss_name='epsilon_insensitive_loss'; 
    loss_gridsel='linear_loss';
    algo=svr(kernel('rbf')); 
    algo.optimizer='andre';
    algo2=svm(kernel('linear'));
    algo2.optimizer='andre';

    example_object=example_lasso_reg; 
    example_object.seed=50+j;
    % example_object.seed=10*j;
    example_object.n=sample;
    example_object.dim=dim;
    example_object.corrdim=imp_dim;
    example_object.cov_r=cov_range;
    example_object.coef=coef;
    [d,example_object]=train(example_object);
    
    % test data
    test_object=example_lasso_reg; 
    test_object.seed=200+j;
    test_object.n=sample;
    test_object.dim=dim;
    test_object.corrdim=imp_dim;
    test_object.coef=coef;
    test_object.cov_r=cov_range;

    [d_test,test_object]=train(test_object);

    cross_validation_names={'rbf','C'};        
    cross_validation_values= {[1,2,3,4,5],[0.01,0.1,1,10,100]};
    g=gridsel(param(algo,cross_validation_names,cross_validation_values),{'score=cv','score.folds=5'});
    g.loss=loss_gridsel;
    [ret_d,ret_algo]=train(g,d);
    algo=ret_algo.best;
%     
cross_validation_names={'C'};        
cross_validation_values= {[0.01,0.1,1,10,100]};
g=gridsel(param(algo2,cross_validation_names,cross_validation_values),{'score=cv','score.folds=5'});
g.loss=loss_gridsel;
[ret_d,ret_algo]=train(g,d);
algo2=ret_algo.best;


    % SVM with risk RFE

    rfe_object=risk_rfe; 
    rfe_object.child=algo;
    rfe_object.loss=loss_name;
    rfe_object.speed=10;
    rfe_object.output_rank=1;
    rfe_object.feat=1; 

    [d_rfe,rfe_object]=train(rfe_object,d);
    val=rfe_object.val; ranks(j,:)=rfe_object.rank;

    svr_train=algo; svr_train.loss=loss_name;
    train_X=d.X;
    train_X_rfe=train_X(:,ranks(j,1:imp_dim));
    train_Y=d.Y;
    d_rfe_train=data(train_X_rfe,train_Y);
    [d_rfe_train,svr_train]=train(svr_train,d_rfe_train);

    % misclassification rates

    

    test_X=d_test.X;
    test_X_rfe=test_X(:,ranks(j,1:imp_dim));
    test_Y=d_test.Y;
    d_rfe_test=data(test_X_rfe,test_Y);
    svr_test=test(svr_train,d_rfe_test,'quadratic_loss');


    % SVM with Guyon's RFE

    rfe_object=rfe; 
    rfe_object.child=algo;
    rfe_object.loss=loss_name;
    rfe_object.speed=10;
    rfe_object.output_rank=1;
    rfe_object.feat=1; 

    [d_rfe,rfe_object]=train(rfe_object,d);
    ranks_g(j,:)=rfe_object.rank;

    svr_train=algo; svr_train.loss=loss_name;
    train_X=d.X;
    train_X_rfe=train_X(:,ranks_g(j,1:imp_dim));
    train_Y=d.Y;
    d_rfe_train=data(train_X_rfe,train_Y);
    [d_rfe_train,svr_train]=train(svr_train,d_rfe_train);

    % error rates

    test_X=d_test.X;
    test_X_rfe=test_X(:,ranks_g(j,1:imp_dim));
    test_Y=d_test.Y;
    d_rfe_test=data(test_X_rfe,test_Y);
    svr_test_grfe=test(svr_train,d_rfe_test,'quadratic_loss');





   % SVM without RFE

    svr_object=algo; svr_object.loss=loss_name;
    [d_svr_test,svr_object]=train(svr_object,d);

    svr_test_nrfe=test(svr_object,d_test,'quadratic_loss');


    % Regression with LASSO

    Cov=get_X(d); Y=get_Y(d);
    model_lasso=glmnet(Cov, Y);
    beta=model_lasso.beta;

    for i=1:size(beta,2)
        l=find(beta(:,i)~=0);
        if length(l)>=imp_dim;
            right_col=i;
            break;
        end
    end
    

    lasso_set=find(beta(:,right_col)~=0);
    ranks_l(j,1:length(lasso_set))=lasso_set;
 
    Cov_lasso=Cov(:,lasso_set);
    p = regress(Y,[ones(size(Cov_lasso,1),1),Cov_lasso]);
    Cov_lasso_test=d_test.X(:,lasso_set);
    yfit = [ones(size(Cov_lasso,1),1),Cov_lasso_test]*p;
    yresid = d_test.Y - yfit;
    MSresid = mean(yresid.^2);
    r_grfe(j)=svr_test_grfe.Y;
   rates_reg(j,:)=[svr_test.Y,svr_test_nrfe.Y, svr_test_grfe.Y, MSresid];
    
end;

%rates_new=rates_reg([1,3,5,6,7,8,9,10,11,12],:);
means=mean(rates_reg);
std_devs=std(rates_reg);


ranks_corr=ranks(1:50,1:imp_dim);
ranks_corr_g=ranks_g(1:50,1:imp_dim);
ranks_corr_l=ranks_l(1:50,1:imp_dim);

ind1=zeros(50,1);
ind2=zeros(50,1);
ind3=zeros(50,1);
ind4=zeros(50,1);
ind5=zeros(50,1);
ind6=zeros(50,1);
ind7=zeros(50,1);
ind8=zeros(50,1);
ind9=zeros(50,1);
for j=1:50
    ind1(j)=isempty(setdiff(corrset,ranks_corr(j,:)));
    ind2(j)=length(setdiff(corrset,ranks_corr(j,:)))==1;
    ind3(j)=length(setdiff(corrset,ranks_corr(j,:)))>1;
    ind4(j)=isempty(setdiff(corrset,ranks_corr_g(j,:)));
    ind5(j)=length(setdiff(corrset,ranks_corr_g(j,:)))==1;
    ind6(j)=length(setdiff(corrset,ranks_corr_g(j,:)))>1;
    ind7(j)=isempty(setdiff(corrset,ranks_l(j,find(ranks_l(j,:)~=0))));
    ind8(j)=length(setdiff(corrset,ranks_l(j,find(ranks_l(j,:)~=0))))==1;
    ind9(j)=length(setdiff(corrset,ranks_l(j,find(ranks_l(j,:)~=0))))>1;
%     ind13(j)=isempty(setdiff(corrset,ranks_corr_reg_(:,j)));
%     ind14(j)=length(setdiff(corrset,ranks_corr_reg_g(:,j)))==1;
%     ind15(j)=length(setdiff(corrset,ranks_corr_reg_g(:,j)))>1;
%      ind10(j)=sum(ranks_corr_s(j,:))==imp_dim;
%      ind11(j)=sum(ranks_corr_s(j,:))==imp_dim-1;
%      ind12(j)=sum(ranks_corr_s(j,:))<imp_dim-1;
%     

end


mean(ind1)
mean(ind2)
mean(ind3)
mean(ind4)
mean(ind5)
mean(ind6)



mean(ind7)
mean(ind8)
mean(ind9)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mean(ind10)
mean(ind11)
mean(ind12)



% set(0,'DefaultAxesFontSize',12);
% plot(val(end:-1:1),'.')
% hold on;
% plot(val(end:-1:1))
% title('Reverse Scree Graph');
% xlabel('Cycle number');
% ylabel('Value of the Objective Function');

for j=1:n_sim;
    corrset=1:imp_dim; 
    loss_name='epsilon_insensitive_loss'; 
    loss_gridsel='linear_loss';
    algo=svr(kernel('rbf')); 
    algo.optimizer='andre';
    example_object=example_lasso_reg; 
    example_object.seed=j;
    % example_object.seed=10*j;
    example_object.n=sample;
    example_object.dim=dim;
    example_object.corrdim=imp_dim;
    example_object.coef_r=coef_range;
    example_object.cov_r=cov_range;
    [d,example_object]=train(example_object);

    Cov=get_X(d); Y=get_Y(d);
    % in case we need R
    data_to_write=[Cov,Y]; 
    dlmwrite(strcat('data_2_run_reg', num2str(j), '.txt'),data_to_write);
end




% SVR with L1 penalty 
% 
% X=get_X(d); y=get_Y(d);
% nVars=size(X,2);
% X = [ones(size(X,1),1),X];
% lambda = [0;10*ones(nVars,1)];
% epsilon=0.1;
% % funObj=@(w)LogisticLoss(w,X,y); 
% funObj = @(w)eps_insen_loss(w,X,y,epsilon);
% wSVR_l1 = L1General2_PSSgb(funObj,zeros(nVars+1,1),lambda);
