n.sim=100;
sample=200;
dim=30;
imp_dim=7;
val=zeros(dim-1,n.sim);
val_reg=zeros(dim-1,n.sim);
ranks=zeros(dim,n.sim);
ranks_reg=zeros(dim,n.sim);
ranks_g=zeros(dim,n.sim);
ranks_s=zeros(dim,n.sim);
ranks_reg_g=zeros(dim,n.sim);
svm_test=zeros(n.sim,1);
svm_test_reg=zeros(n.sim,1);
svm_test_scadsvc=zeros(n.sim,1);
svm_test_g=zeros(n.sim,1);
svm_test_reg_g=zeros(n.sim,1);



% population=[-1,-0.5,0.5,1];
% coeff=randsample(population,imp_dim,true);
% coef=[1,-1,-0.5,0.5]; % SVM with Gauss d=15, d0=4;
% coeff=[-0.5,1,0.5,-1]; % SVR with linear d=15, d0=4;
% coeff=[-0.5,-0.5,-1,0.5,-0.5,0.5,-1]; 
 coef=[-1,-1.5,-1,1.5,-1.5,1.5,-1]; % SVM with Gauss d=30, d0=7;
 coeff=[1,1,-1,0.5,0.5,0.5,-0.5]; % SVR with linear d=30, d0=7;
%  coef=[-0.5,-1,1]; % SVM with Gauss d=50, d0=3;
%  coeff=[1,0.5,-1]; % SVR with Linear d=50, d0=3;


for i=1:n.sim;
    
%% Classification




 corrset=[1:imp_dim];
 loss_name='hinge_loss';
 loss_gridsel='class_loss';
 algo=svm(kernel('rbf',1));
 algo.optimizer='andre';

 example_object=toy_n;
 example_object.l=sample;
 example_object.n=dim;
 example_object.relevant_n=imp_dim;
 example_object.coef=[coef';zeros(example_object.n-example_object.relevant_n,1)];
 example_object.seed=i;
 test_object=toy_n;
 test_object.l=sample;
 test_object.n=dim;
 test_object.relevant_n=imp_dim;
 test_object.coef=[coef';zeros(example_object.n-example_object.relevant_n,1)];
 test_object.seed=100+i;  

[d,example_object]=train(example_object);
[d_test,test_object]=train(test_object);

cross_validation_names={'rbf','C'};              %     cross_validation_names,...
 cross_validation_values= {[1,2,3,4],[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;
     
% Risk rfe

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

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

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

% misclassification rates

test_X=d_test.X;
test_X_rfe=test_X(:,ranks(1:imp_dim,i));
test_Y=d_test.Y;
d_rfe_test=data(test_X_rfe,test_Y);
t_svm=test(svm_train,d_rfe_test,loss_gridsel);
svm_test(i)=t_svm.Y;

% Guyon rfe

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

rfe_object.feat=1;
[d_rfe,rfe_object]=train(rfe_object,d);
%val(:,i)=rfe_object.val;
ranks_g(:,i)=rfe_object.rank;

svm_train=algo; svm_train.loss=loss_name;
train_X=d.X;
train_X_grfe=train_X(:,ranks_g(1:imp_dim,i));
train_Y=d.Y;
d_grfe_train=data(train_X_grfe,train_Y);
[d_grfe_train,svm_train]=train(svm_train,d_grfe_train);

% misclassification rates

test_X=d_test.X;
test_X_grfe=test_X(:,ranks_g(1:imp_dim,i));
test_Y=d_test.Y;
d_grfe_test=data(test_X_grfe,test_Y);
t_svm_g=test(svm_train,d_grfe_test,'class_loss');
svm_test_g(i)=t_svm_g.Y;


% SCAD SVM 

Cov=get_X(d); Y=get_Y(d);
svm_object=algo; svm_object.loss=loss_name;
[d_svm,svm_object]=train(svm_object,d);

dim=size(Cov,2);
w_init= svm_object.Xsv.X'*svm_object.alpha;
b_init=svm_object.b0;

lam=[3:-0.05:0];

for j=1:length(lam);
    [w,b,xind]=scadsvc_spider(Cov,Y,w_init,b_init,3.7,lam(j));
    if length(xind)>=imp_dim;
        break;
    end;
end;

ranks_s(:,i)=ismember(1:dim,xind);


svm_train=algo; svm_train.loss=loss_name;
train_X=d.X;
train_X_scadsvc=train_X(:,xind);
train_Y=d.Y;
d_scadsvc_train=data(train_X_scadsvc,train_Y);
[d_scadsvc_train,svm_train]=train(svm_train,d_scadsvc_train);

test_X=d_test.X;
test_X_scadsvc=test_X(:,xind);
test_Y=d_test.Y;
d_scadsvc_test=data(test_X_scadsvc,test_Y);
t_svm_c=test(svm_train,d_scadsvc_test,loss_gridsel);
svm_test_scadsvc(i)=t_svm_c.Y;

%% Regression

corrset=[1:imp_dim];    
loss_name='epsilon_insensitive_loss'; 
loss_gridsel='linear_loss';
loss_test='quadratic_loss';
algo=svr(kernel('linear')); 
algo.optimizer='andre';

example_object=example_reg; 
example_object.seed=i;
example_object.n=sample;
example_object.dim=dim;
example_object.corrdim=imp_dim;
example_object.coef=[coeff';zeros(example_object.dim-example_object.corrdim,1)];

test_object=example_reg; 
test_object.seed=100+i;
test_object.n=sample;
test_object.dim=dim;
test_object.corrdim=imp_dim;
test_object.coef=[coeff';zeros(example_object.dim-example_object.corrdim,1)];

[d,example_object]=train(example_object);
[d_test,test_object]=train(test_object);

cross_validation_names={'rbf','C'};              %     cross_validation_names,...
 cross_validation_values= {[1,2,3,4],[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;
     

% Risk rfe

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

rfe_object.feat=1;
[d_rfe,rfe_object]=train(rfe_object,d);
val_reg(:,i)=rfe_object.val;
ranks_reg(:,i)=rfe_object.rank;

svm_train=algo; svm_train.loss=loss_name;
train_X=d.X;
train_X_rfe=train_X(:,ranks_reg(1:imp_dim,i));
train_Y=d.Y;
d_rfe_train=data(train_X_rfe,train_Y);
[d_rfe_train,svm_train]=train(svm_train,d_rfe_train);

% measurement error

test_X=d_test.X;
test_X_rfe=test_X(:,ranks_reg(1:imp_dim,i));
test_Y=d_test.Y;
d_rfe_test=data(test_X_rfe,test_Y);
t_svm_reg=test(svm_train,d_rfe_test,loss_test);
svm_test_reg(i)=t_svm_reg.Y;

% Guyon rfe

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

rfe_object.feat=1;
[d_rfe,rfe_object]=train(rfe_object,d);
%val_reg(:,i)=rfe_object.val;
ranks_reg_g(:,i)=rfe_object.rank;

svm_train=algo; svm_train.loss=loss_name;
train_X=d.X;
train_X_rfe=train_X(:,ranks_reg_g(1:imp_dim,i));
train_Y=d.Y;
d_rfe_train=data(train_X_rfe,train_Y);
[d_rfe_train,svm_train]=train(svm_train,d_rfe_train);

% measurement error

test_X=d_test.X;
test_X_rfe=test_X(:,ranks_reg_g(1:imp_dim,i));
test_Y=d_test.Y;
d_rfe_test=data(test_X_rfe,test_Y);
t_svm_reg_g=test(svm_train,d_rfe_test,loss_test);
svm_test_reg_g(i)=t_svm_reg_g.Y;



i
end;

%% Evaluation

mean_val=mean(val,2)
mean_ranks=mean(ranks,2)

ranks_corr=ranks(1:imp_dim,:);
ranks_corr_reg=ranks_reg(1:imp_dim,:);
ranks_corr_reg_g=ranks_reg_g(1:imp_dim,:);
ranks_corr_g=ranks_g(1:imp_dim,:);
ranks_corr_s=ranks_s(1:imp_dim,:);


ind1=zeros(n.sim,1);
ind2=zeros(n.sim,1);
ind3=zeros(n.sim,1);
ind4=zeros(n.sim,1);
ind5=zeros(n.sim,1);
ind6=zeros(n.sim,1);
ind7=zeros(n.sim,1);
ind8=zeros(n.sim,1);
ind9=zeros(n.sim,1);
ind10=zeros(n.sim,1);
ind11=zeros(n.sim,1);
ind12=zeros(n.sim,1);

ind13=zeros(n.sim,1);
ind14=zeros(n.sim,1);
ind15=zeros(n.sim,1);
for j=1:n.sim
    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;
    ind10(j)=isempty(setdiff(corrset,ranks_corr_reg(:,j)));
    ind11(j)=length(setdiff(corrset,ranks_corr_reg(:,j)))==1;
    ind12(j)=length(setdiff(corrset,ranks_corr_reg(:,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;
    ind13(j)=isempty(setdiff(corrset,ranks_corr_reg_g(:,j)));
    ind14(j)=length(setdiff(corrset,ranks_corr_reg_g(:,j)))==1;
    ind15(j)=length(setdiff(corrset,ranks_corr_reg_g(:,j)))>1;
    ind7(j)=sum(ranks_corr_s(:,j))==imp_dim;
    ind8(j)=sum(ranks_corr_s(:,j))==imp_dim-1;
    ind9(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)
mean(ind13)
mean(ind14)
mean(ind15)

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

mean(svm_test)
mean(svm_test_g)
mean(svm_test_scadsvc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mean(svm_test_reg)
mean(svm_test_reg_g)
