function [means_train means_test std_devs_train std_devs_test]= nonlinear_reg_col(n_sim,sample,dim,imp_dim,cov_r,coef,script)

% Generate Data


filename1='/home/sdasgup2/cluster/vs_rfe/new_nonlinear';
filename2='/home/sdasgup2/cluster/vs_rfe';
%filename3='/home/sdasgup2/cluster/vs_rfe/new_nonlinear/regression/nonlinear_reg.mat';
filename3=strcat('/home/sdasgup2/cluster/vs_rfe/new_nonlinear/regression/nonlinear_reg_sample',num2str(sample),'_dim',num2str(dim),'_impdim',num2str(imp_dim),'_nsim',num2str(n_sim),'_covr',num2str(cov_r),'_coef',num2str(coef),'_script',num2str(script),'.mat');
cd(filename1);
addpath(genpath(filename2))


warning('off','MATLAB:dispatcher:InexactCaseMatch')
%n_sim=100; sample=200; dim=50; imp_dim=3; coef_range=5; cov_r=2;

coef=coef*ones(imp_dim,1);

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);
lin_test_lasso_vec=zeros(n_sim,1);
lin_train_lasso_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='epsilon_insensitive_loss'; 
    loss_gridsel='linear_loss';
    algo=svr(kernel('rbf',1)); 
    algo.optimizer='andre';
    
    example_object=nonlinear_reg_col; 
    example_object.seed=seeds(j);
    example_object.n=sample;
    example_object.dim=dim;
    example_object.corrdim=imp_dim;
    example_object.cov_r=cov_r;
    example_object.coef=coef;
    
    [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 data
    test_object=nonlinear_reg_col; 
    test_object.seed=seeds_test(j);
    test_object.n=sample;
    test_object.dim=dim;
    test_object.corrdim=imp_dim;
    test_object.coef=coef;
    test_object.cov_r=cov_r;

    [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,'quadratic_loss');
    svm_test_nrfe_vec(j)=svm_test_nrfe.Y;
    svm_train_nrfe=test(svm_object_trained,d,'quadratic_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);
    % measurement error 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,'quadratic_loss');
    svm_test_rfe_vec(j)=svm_test_rfe.Y;
    svm_train_rfe=test(svm_object_rfe_trained,d_rfe_train,'quadratic_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);

    % measurement error 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,'quadratic_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,'quadratic_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);

    % measurement error 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,'quadratic_loss');
    svm_test_grfe_vec(j)=svm_test_grfe.Y;
    svm_train_grfe=test(svm_object_rfe_g_trained,d_rfe_train_g,'quadratic_loss');
    svm_train_grfe_vec(j)=svm_train_grfe.Y;

    %% linear regression with LASSO

    Cov=get_x(d); Y=get_y(d);
    lasso_cv = cvglmnet(Cov, Y, 5,[],'response','gaussian',glmnetSet,0);
    options=struct('lambda',lasso_cv.lambda_min);
    lasso_fit = glmnet(Cov, Y,'gaussian',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,'normal');
    Cov_lasso_test=d_test.X(:,lasso_set);
    y_hat = glmval(beta_fin, Cov_lasso_test, 'identity');
    y_train_hat = glmval(beta_fin, Cov_lasso, 'identity');
    y_resid=d_test.Y - y_hat;
    mes_err = mean(y_resid.^2);
    y_train_resid=d.Y - y_train_hat;
    mes_err_train = mean(y_train_resid.^2);

    lin_test_lasso_vec(j)=mes_err;
    lin_train_lasso_vec(j)=mes_err_train;



    save(filename3);
        
end

rates_train=[svm_train_nrfe_vec svm_train_rfe_vec svm_train_rfe_cp_vec svm_train_grfe_vec lin_train_lasso_vec];
rates_test=[svm_test_nrfe_vec svm_test_rfe_vec svm_test_rfe_cp_vec svm_test_grfe_vec lin_test_lasso_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);


end

