% Generate Data

clear all;
n_sim=100;
rates=zeros(n_sim,5);
rrr=zeros(n_sim,2);
warning('off','MATLAB:dispatcher:InexactCaseMatch')
sample=200; dim=50; imp_dim=3; x_range=2^(1/imp_dim); y_range=1;
ranks=zeros(n_sim,dim);
ranks_l=zeros(n_sim,dim);
ranks_s=zeros(n_sim,dim);
ranks_g=zeros(n_sim,dim);
val=zeros(dim-1,n_sim);
a={[4,5],[5,10],[1,4,10],[6,8,10],[5,9,10],[2,4,6],[3,9],[4,6],[5,6,8],[2,3]};
indices=[1,3,5,6,7,8,9,10,11,12];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Seeds to consider 21-70 %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for j=1:n_sim;
corrset=1:imp_dim;    
loss_name='hinge_loss';
loss_gridsel='class_loss';
algo=svm(kernel('rbf',1));
algo.optimizer='andre';
algo2=svm(kernel('linear'));
algo2.optimizer='andre';

example_object=example_lasso_cor; 
example_object.seed=j;
example_object.n=sample;
example_object.dim=dim;
example_object.corrdim=imp_dim;
example_object.coef_r=x_range;
example_object.output_r=y_range;


test_object=example_lasso_cor; 
test_object.seed=100+j;
test_object.n=sample;
test_object.dim=dim;
test_object.corrdim=imp_dim;
test_object.coef_r=x_range;
test_object.output_r=y_range;

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



cross_validation_names={'rbf','C'};        
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;

% 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=100;
rfe_object.output_rank=1;
rfe_object.feat=1; 

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

svm_train=algo; svm_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,svm_train]=train(svm_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);
svm_test=test(svm_train,d_rfe_test,'class_loss');


% SVM with Guyon's 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);
ranks_g(j,:)=rfe_object.rank;

svm_train=algo; svm_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,svm_train]=train(svm_train,d_rfe_train);

% misclassification 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);
svm_test_grfe=test(svm_train,d_rfe_test,'class_loss');


% % % SVM without RFE
%  
%  svm_object=algo; svm_object.loss=loss_name;
%  [d_svm,svm_object]=train(svm_object,d);
% svm_test_nrfe=test(svm_object,d_test,'class_loss');

% % logistic regression with LASSO
% 
% Cov=get_X(d); Y=double(get_Y(d)==1)+1;
% lasso_fit = glmnet(Cov, Y,'binomial');
% 
% beta=lasso_fit.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;
%  
% % cv=cvglmnet(Cov,Y,5,[],'response','gaussian',glmnetSet,1);
% % lam_opt=cv.lambda_1se;
% % lasso_fit = glmnet(Cov, Y); % regression
% % 
% % beta= glmnetCoef(lasso_fit,lam_opt);
% % lasso_set=find(beta(2:end)~=0);
% % vars_sel(i)=length(lasso_set);
% 
% Cov_lasso=Cov(:,lasso_set);
% p = glmfit(Cov_lasso,(get_Y(d)==1),'binomial');
% Cov_lasso_test=d_test.X(:,lasso_set);
% y_hat = glmval(p, Cov_lasso_test, 'logit');
% class=(y_hat > 0.5)-(y_hat <= 0.5);
% misclass = mean(d_test.Y ~= class);
% 
% 
% 
% 
% 
% 

% % SCAD SVM 
% 
%  svm_object=algo2; svm_object.loss=loss_name;
%  [d_svm,svm_object]=train(svm_object,d);
% 
% Cov=get_X(d); Y=get_Y(d);
% dim=size(Cov,2);
% w_init= svm_object.Xsv.X'*svm_object.alpha;
% b_init=svm_object.b0;
% 
% lam=[0.1:-0.005:0];
% 
% for k=1:length(lam);
%     [w,b,xind]=scadsvc_spider(Cov,Y,w_init,b_init,3.7,lam(k));
%     if length(xind)>=imp_dim;
%         break;
%     end;
% end;
% 
% 
% if xind==0;
%     ranks_s(j,:)=zeros(dim,1);
%     xind=1:dim;
% else ranks_s(j,:)=ismember(1:dim,xind);
% end
%     svm_train=algo2; svm_train.loss=loss_name;
%     train_X=d.X;
%     train_Y=d.Y;
%     test_X=d_test.X;
%     test_Y=d_test.Y;
% 
%     train_X_scadsvc2=train_X(:,xind);
%     d_scadsvc_train2=data(train_X_scadsvc2,train_Y);
%     [d_scadsvc_train2,svm_train2]=train(svm_train,d_scadsvc_train2);
%     test_X_scadsvc2=test_X(:,xind);
%     d_scadsvc_test2=data(test_X_scadsvc2,test_Y);
%     svm_test_scadsvc2=test(svm_train2,d_scadsvc_test2,'class_loss');
%     scad.err(j)=svm_test_scadsvc2.Y;
%     scad_err=svm_test_scadsvc2.Y;
% 
%     
% 
% 
% % SVM with L1 penalty using LPSVM by Mangasarian
% 
% % A=get_X(d); y=get_Y(d); k=size(A,2);
% % [w gamma trainset testset time nu] = lpsvm(A,y,k);
% % delta=10^-3;
% 
% % in case we need R
% %  data_to_write=[A,y]; 
% %  dlmwrite(strcat('data_2_run_', num2str(j), '.txt'),data_to_write);
% 
% % w_seq = textread('data.txt', '', 'delimiter', ',', 'emptyvalue', NaN);
% % 
% % for i=size(w_seq,2):-1:1
% %     l=find(w_seq(:,i)~=0);
% %     if length(l)>=2;
% %         right_col=i;
% %         break;
% %     end
% % end
% % 
% % lasso_set=find(w_seq(:,right_col)~=0);
% 
% % svm_lasso=algo; svm_lasso.loss=loss_name;
% % lasso_set=a{find(indices==j)};
% % Cov_lasso=Cov(:,lasso_set);
% % d_lasso_train=data(Cov_lasso,train_Y);
% % [d_lasso_trained,svm_lasso_new]=train(svm_lasso,d_lasso_train);
% % 
% % Cov_lasso_test=d_test.X(:,lasso_set);
% % d_lasso_test=data(Cov_lasso_test,test_Y);
% % svm_lasso_test=test(svm_lasso_new,d_lasso_test,'class_loss');
% 
% 
% %rates(j,:)=[svm_test.Y, svm_test_nrfe.Y, svm_test_grfe.Y misclass, svm_lasso_test.Y ];
%rates(j,:)=[svm_test.Y, svm_test_nrfe.Y, svm_test_grfe.Y, misclass, scad_err];
rrr(j,:)=[svm_test.Y, svm_test_grfe.Y];
   
 

end;


means=mean(rates(1:50,:));
std_devs=std(rates(1:50,:));



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);
ranks_corr_s=ranks_s(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);
ind10=zeros(50,1);
ind11=zeros(50,1);
ind12=zeros(50,1);

ind13=zeros(n_sim,1);
ind14=zeros(n_sim,1);
ind15=zeros(n_sim,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)


aa=(mean(ranks(1:imp_dim,:))==2);
bb=(mean(ranks_g(1:imp_dim,:))==2);
rates_n=rates(1:50,:);

means_n=mean(rates);
std_devs_n=std(rates);


fid = fopen('test_cor2.txt');
line = fgetl(fid);

k=0;
b=cell(n_sim,1);
rr=zeros(n_sim,1);
while ischar(line)
  k=k+1;
  yy = str2num(line);
  loss_name='hinge_loss';
loss_gridsel='class_loss';
algo=svm(kernel('rbf',1));
algo.optimizer='andre';
example_object=example_lasso; 
example_object.seed=j;
example_object.n=sample;
example_object.dim=dim;
example_object.corrdim=imp_dim;
example_object.coef=coef;
example_object.output_r=y_range;
  
[d,example_object]=train(example_object);

test_object=example_lasso; 
test_object.seed=100+j;
test_object.n=sample;
test_object.dim=dim;
test_object.corrdim=imp_dim;
test_object.coef_r=x_range;
test_object.output_r=y_range;

[d_test,test_object]=train(test_object);

Cov=get_X(d); Y=get_Y(d);

svm_lasso=algo; svm_lasso.loss=loss_name;

lasso_set=yy;
Cov_lasso=Cov(:,lasso_set);
d_lasso_train=data(Cov_lasso,train_Y);
[d_lasso_trained,svm_lasso_new]=train(svm_lasso,d_lasso_train);

Cov_lasso_test=d_test.X(:,lasso_set);
d_lasso_test=data(Cov_lasso_test,test_Y);
svm_lasso_test=test(svm_lasso_new,d_lasso_test,'class_loss');

rr(k)=svm_lasso_test.Y;

line = fgetl(fid);
end

fclose(fid);

rates=[rates, rr];


% no prediction

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

% 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');


warning('on','MATLAB:dispatcher:InexactCaseMatch')
