function [means_train means_test std_devs_train std_devs_test]= nonlinear_class_col_last(n_sim,sample,dim)

% Generate Data


filename2='/home/sdasgup2/cluster/vs_rfe';
filename3=strcat('/home/sdasgup2/cluster/vs_rfe/new_nonlinear/batch_n800/nonlinear_class_col_sample',num2str(sample),'_dim',num2str(dim),'_nsim',num2str(n_sim),'_2feat.mat');
addpath(genpath(filename2))


warning('off','MATLAB:dispatcher:InexactCaseMatch')
%n_sim=100; sample=200; dim=10; 
imp_dim=2; 

svm_train_nrfe_vec=zeros(n_sim,1);
svm_test_nrfe_vec=zeros(n_sim,1);

vars_sel_rfe=cell(n_sim,1);
svm_test_rfe_vec=zeros(n_sim,1);
svm_train_rfe_vec=zeros(n_sim,1);

vars_sel_rfe_cp=cell(n_sim,1);
svm_test_rfe_cp_vec=zeros(n_sim,1);
svm_train_rfe_cp_vec=zeros(n_sim,1);

vars_sel_rfe_g=cell(n_sim,1);
svm_test_grfe_vec=zeros(n_sim,1);
svm_train_grfe_vec=zeros(n_sim,1);

vars_sel_lasso=cell(n_sim,1);
log_test_lasso_vec=zeros(n_sim,1);
log_train_lasso_vec=zeros(n_sim,1);

vars_sel_scad=cell(n_sim,1);
svm_test_scad_vec=zeros(n_sim,1);
svm_train_scad_vec=zeros(n_sim,1);

vars_sel_scad_all=cell(n_sim,1);
svm_test_scad_all_vec=zeros(n_sim,1);
svm_train_scad_all_vec=zeros(n_sim,1);

seeds=zeros(n_sim,1);
seeds_test=zeros(n_sim,1);


for j=1:n_sim;
    seeds(j)=j*rand(1);
    seeds_test(j)=(100+j)*rand(1);
    corrset=1:imp_dim;    
    loss_name='hinge_loss';
    loss_gridsel='class_loss';
    algo=svm(kernel('rbf',1));
    algo.optimizer='andre';
    example_object=nonlinear_class_2feat_int_col; 
    example_object.seed=seeds(j);
    example_object.n=sample;
    example_object.dim=dim;
    
    [d,example_object_trained]=train(example_object);
    cross_validation_names={'rbf','C'};        
    cross_validation_values= {[1,2,3,4,5],[0.001,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;


    test_object=nonlinear_class_2feat_int_col; 
    test_object.seed=seeds_test(j);
    test_object.n=sample;
    test_object.dim=dim;
    
    [d_test,test_object_trained]=train(test_object);

    % SVM without RFE

    svm_object=algo; svm_object.loss=loss_name;
    [d_svm,svm_object_trained]=train(svm_object,d);
    svm_test_nrfe=test(svm_object_trained,d_test,'class_loss');
    svm_test_nrfe_vec(j)=svm_test_nrfe.Y;
    svm_train_nrfe=test(svm_object_trained,d,'class_loss');
    svm_train_nrfe_vec(j)=svm_train_nrfe.Y;

    %SVM with risk RFE

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

    [d_rfe,rfe_object_trained]=train(rfe_object,d);
    val=rfe_object_trained.val; ranks=rfe_object_trained.rank;
    vars_sel_rfe{j}=ranks(1:imp_dim);

    % using naive ranks
    svm_object_rfe=algo; svm_object_rfe.loss=loss_name;
    train_X=d.X;
    train_X_rfe=train_X(:,ranks(1:imp_dim));
    train_Y=d.Y;
    d_rfe_train=data(train_X_rfe,train_Y);
    [d_rfe_trained,svm_object_rfe_trained]=train(svm_object_rfe,d_rfe_train);

    % misclassification rates

    test_X=d_test.X;
    test_X_rfe=test_X(:,ranks(1:imp_dim));
    test_Y=d_test.Y;
    d_rfe_test=data(test_X_rfe,test_Y);
    svm_test_rfe=test(svm_object_rfe_trained,d_rfe_test,'class_loss');
    svm_test_rfe_vec(j)=svm_test_rfe.Y;
    svm_train_rfe=test(svm_object_rfe_trained,d_rfe_train,'class_loss');
    svm_train_rfe_vec(j)=svm_train_rfe.Y;


    % Test and training error for RFE using change point model
    val2=val(end:-1:1);
    x=1:length(val);
    err_sos=zeros(length(val),1);
    for k=1:length(val)-1;
        x0=k;
        X_x0=[ones(length(val),1),x',(x'-x0).*(x'>x0),(x'.^2-x0^2).*(x'>x0)];
        beta_hat=(X_x0'*X_x0)\(X_x0'*val2');
        err_sos(k)=sum((val2'-X_x0*beta_hat).^2);
    end
    X_full=[ones(length(val),1),x'];
    beta_full=(X_full'*X_full)\(X_full'*val2');
    err_sos(length(val))=sum((val2'-X_full*beta_full).^2);
    [c,j_imp]=min(err_sos);
    vars_sel_rfe_cp{j}=ranks(1:end-j_imp);

    train_X_rfe_cp=train_X(:,vars_sel_rfe_cp{j});
    train_Y=d.Y;
    d_rfe_train_cp=data(train_X_rfe_cp,train_Y);
    [d_rfe_trained_cp,svm_object_rfe_trained_cp]=train(svm_object_rfe,d_rfe_train_cp);

    % misclassification rates
    test_X=d_test.X;
    test_X_rfe_cp=test_X(:,vars_sel_rfe_cp{j});
    test_Y=d_test.Y;
    d_rfe_test_cp=data(test_X_rfe_cp,test_Y);
    svm_test_rfe_cp=test(svm_object_rfe_trained_cp,d_rfe_test_cp,'class_loss');
    svm_test_rfe_cp_vec(j)=svm_test_rfe_cp.Y;
    svm_train_rfe_cp=test(svm_object_rfe_trained_cp,d_rfe_train_cp,'class_loss');
    svm_train_rfe_cp_vec(j)=svm_train_rfe_cp.Y;


    % SVM with Guyon's RFE

    rfe_object_g=rfe; 
    rfe_object_g.child=algo;
    rfe_object_g.loss=loss_name;
    rfe_object_g.speed=dim;
    rfe_object_g.output_rank=1;
    rfe_object_g.feat=1; 

    [d_rfe_g,rfe_object_trained_g]=train(rfe_object_g,d);
    ranks_g=rfe_object_trained_g.rank;
    vars_sel_rfe_g{j}=ranks_g(1:imp_dim);


    svm_object_rfe_g=algo; svm_object_rfe_g.loss=loss_name;
    train_X=d.X;
    train_X_rfe_g=train_X(:,ranks_g(1:imp_dim));
    train_Y=d.Y;
    d_rfe_train_g=data(train_X_rfe_g,train_Y);
    [d_rfe_trained_g,svm_object_rfe_g_trained]=train(svm_object_rfe_g,d_rfe_train_g);

    % misclassification rates

    test_X=d_test.X;
    test_X_rfe_g=test_X(:,ranks_g(1:imp_dim));
    test_Y=d_test.Y;
    d_rfe_test_g=data(test_X_rfe_g,test_Y);
    svm_test_grfe=test(svm_object_rfe_g_trained,d_rfe_test_g,'class_loss');
    svm_test_grfe_vec(j)=svm_test_grfe.Y;
    svm_train_grfe=test(svm_object_rfe_g_trained,d_rfe_train_g,'class_loss');
    svm_train_grfe_vec(j)=svm_train_grfe.Y;

    % logistic regression with LASSO

    Cov=get_x(d); Y=double(get_y(d)==1); Y_lass=Y+1; 
    Cov_test=get_x(d_test); Y_test=double(get_y(d_test)==1);
    lasso_cv = cvglmnet(Cov, Y_lass, 5,[],'response','binomial',glmnetSet,0);
    options=struct('lambda',lasso_cv.lambda_min);
    lasso_fit = glmnet(Cov, Y_lass,'binomial',glmnetSet(options));

    beta=lasso_fit.beta;

    lasso_set=find(beta~=0);
    vars_sel_lasso{j}=lasso_set;

    Cov_lasso=Cov(:,lasso_set);
    beta_fin = glmfit(Cov_lasso,Y,'binomial');
    Cov_lasso_test=Cov_test(:,lasso_set);
    y_hat = glmval(beta_fin, Cov_lasso_test, 'logit');
    y_train_hat = glmval(beta_fin, Cov_lasso, 'logit');
    class=(y_hat > 0.5);
    class_train=(y_train_hat > 0.5);
    misclass = mean(Y_test ~= class);
    misclass_train = mean(Y ~= class_train);
    log_test_lasso_vec(j)=misclass;
    log_train_lasso_vec(j)=misclass_train;



    % SCAD SVM 

    Cov=get_x(d); Y=get_y(d);
    dim=size(Cov,2);
    w_init= svm_object_trained.Xsv.X'*svm_object_trained.alpha;
    b_init=svm_object_trained.b0;
    [w,b,iter]=scadsvc_spider_all(Cov,Y,w_init,b_init);
    [w1,b1,xind]=scadsvc_spider(Cov,Y,w_init,b_init,3.7,0.12);
    tol2=0.001;
    ind=abs(w/norm(w))>tol2;

    vars_sel_scad{j}=xind;
    vars_sel_scad_all{j}=ind;

    svm_object_scad=algo; svm_object_scad.loss=loss_name;
    train_X=d.X;
    train_X_scadsvc_all=train_X(:,ind);
    train_Y=d.Y;
    d_scadsvc_train_all=data(train_X_scadsvc_all,train_Y);
    [d_scadsvc_trained_all,svm_object_scad_trained_all]=train(svm_object_scad,d_scadsvc_train_all);

    test_X=d_test.X;
    test_X_scadsvc_all=test_X(:,ind);
    test_Y=d_test.Y;
    d_scadsvc_test_all=data(test_X_scadsvc_all,test_Y);
    svm_test_scadsvc_all=test(svm_object_scad_trained_all,d_scadsvc_test_all,'class_loss');
    svm_train_scadsvc_all=test(svm_object_scad_trained_all,d_scadsvc_train_all,'class_loss');

    svm_test_scad_all_vec(j)=svm_test_scadsvc_all.Y;
    svm_train_scad_all_vec(j)=svm_train_scadsvc_all.Y;


    if xind~=0
        train_X_scadsvc=train_X(:,xind);
        d_scadsvc_train=data(train_X_scadsvc,train_Y);
        [d_scadsvc_trained,svm_object_scad_trained]=train(svm_object_scad,d_scadsvc_train);
        test_X_scadsvc=test_X(:,xind);
        d_scadsvc_test=data(test_X_scadsvc,test_Y);
        svm_test_scadsvc=test(svm_object_scad_trained,d_scadsvc_test,'class_loss');
        svm_train_scadsvc=test(svm_object_scad_trained,d_scadsvc_train,'class_loss');
        svm_test_scad_vec(j)=svm_test_scadsvc.Y;
        svm_train_scad_vec(j)=svm_train_scadsvc.Y;

    end
    save(filename3);
        
end

rates_train=[svm_train_nrfe_vec svm_train_rfe_vec svm_train_rfe_cp_vec svm_train_grfe_vec log_train_lasso_vec svm_train_scad_vec];
rates_test=[svm_test_nrfe_vec svm_test_rfe_vec svm_test_rfe_cp_vec svm_test_grfe_vec log_test_lasso_vec svm_test_scad_vec];

means_train=mean(rates_train);
std_devs_train=std(rates_train);
means_test=mean(rates_test);
std_devs_test=std(rates_test);

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

save(filename3);

 % CLASSIFICATION


 corr=1:imp_dim;
 
 count_l=zeros(1,n_sim);
 count_r=zeros(1,n_sim);
 count_rcp=zeros(1,n_sim);
 count_g=zeros(1,n_sim);
 count_s=zeros(1,n_sim);
 count_sall=zeros(1,n_sim);
  for i=1:n_sim;
     vars_sel_l=vars_sel_lasso{i};
     vars_sel_r=vars_sel_rfe{i};
     vars_sel_rcp=vars_sel_rfe_cp{i};
     vars_sel_g=vars_sel_rfe_g{i};
     vars_sel_s=vars_sel_scad{i};
     vars_sel_sall=vars_sel_scad_all{i};
     
     count_l(i)=sum(ismember(corr,vars_sel_l));
     count_r(i)=sum(ismember(corr,vars_sel_r));
     count_rcp(i)=sum(ismember(corr,vars_sel_rcp));
     count_g(i)=sum(ismember(corr,vars_sel_g));
     count_s(i)=sum(ismember(corr,vars_sel_s));
     count_sall(i)=sum(ismember(corr,vars_sel_sall));
    
  end
 
count_l=count_l(rates_test(:,5)~=0);
count_r=count_r(rates_test(:,2)~=0);
count_rcp=count_rcp(rates_test(:,3)~=0);
count_g=count_g(rates_test(:,4)~=0);
count_s=count_s(rates_test(:,6)~=0);  
  
 prop_none_l=mean(count_l==imp_dim);
 prop_1_l=mean(count_l==imp_dim-1);
 prop_more_l=mean(count_l<imp_dim-1);
 
 prop_none_r=mean(count_r==imp_dim);
 prop_1_r=mean(count_r==imp_dim-1);
 prop_more_r=mean(count_r<imp_dim-1);
 
 prop_none_rcp=mean(count_rcp==imp_dim);
 prop_1_rcp=mean(count_rcp==imp_dim-1);
 prop_more_rcp=mean(count_rcp<imp_dim-1);
 
 prop_none_g=mean(count_g==imp_dim);
 prop_1_g=mean(count_g==imp_dim-1);
 prop_more_g=mean(count_g<imp_dim-1);
 
 prop_none_s=mean(count_s==imp_dim);
 prop_1_s=mean(count_s==imp_dim-1);
 prop_more_s=mean(count_s<imp_dim-1);
 
 prop_none_sall=mean(count_sall==imp_dim);
 prop_1_sall=mean(count_sall==imp_dim-1);
 prop_more_sall=mean(count_sall<imp_dim-1);
 
 prop_l=[prop_none_l prop_1_l prop_more_l];
 prop_r=[prop_none_r prop_1_r prop_more_r];
 prop_rcp=[prop_none_rcp prop_1_rcp prop_more_rcp];
 prop_g=[prop_none_g prop_1_g prop_more_g];
 prop_s=[prop_none_s prop_1_s prop_more_s];
 prop_sall=[prop_none_sall prop_1_sall prop_more_sall];
 
 prop_l
 prop_r
 prop_rcp
 prop_g
 prop_s
 prop_sall
 
 
[  mean(rates_test(rates_test(:,2)~=0,2))
    mean(rates_test(rates_test(:,3)~=0,3))
    mean(rates_test(rates_test(:,4)~=0,4))
    mean(rates_test(rates_test(:,5)~=0,5))
    mean(rates_test(rates_test(:,6)~=0,6))
    mean(rates_test(rates_test(:,1)~=0,1))]

[  std(rates_test(rates_test(:,2)~=0,2))
    std(rates_test(rates_test(:,3)~=0,3))
    std(rates_test(rates_test(:,4)~=0,4))
    std(rates_test(rates_test(:,5)~=0,5))
    std(rates_test(rates_test(:,6)~=0,6))
    std(rates_test(rates_test(:,1)~=0,1))]

 save(filename3);


end

