clear all;
filename1='/home/sdasgup2/cluster/vs_rfe/data';
filename2='/home/sdasgup2/cluster/vs_rfe';
filename3='/home/sdasgup2/cluster/vs_rfe/data/crime/crime_cluster_new.mat';
cd(filename1);
addpath(genpath(filename2))

[fid err] = fopen('communities_data.txt');
dd=textscan(fid,['%f%f%f%s',repmat('%f',1,124)],'TreatAsEmpty', '?','Delimiter',',');
fclose(fid);
dt=horzcat(dd{5:end});
header={'fold ','population ', 'householdsize ','racepctblack ', 'racePctWhite ','racePctAsian ', 'racePctHisp ', 'agePct12t21 ', 'agePct12t29 ', 'agePct16t24 ', 'agePct65up ', 'numbUrban ', 'pctUrban ', 'medIncome ', 'pctWWage ', 'pctWFarmSelf ', 'pctWInvInc ', 'pctWSocSec ', 'pctWPubAsst ', 'pctWRetire ', 'medFamInc ', 'perCapInc ', 'whitePerCap ', 'blackPerCap ', 'indianPerCap ', 'AsianPerCap ', 'OtherPerCap ', 'HispPerCap ', 'NumUnderPov ', 'PctPopUnderPov ', 'PctLess9thGrade ', 'PctNotHSGrad ', 'PctBSorMore ', 'PctUnemployed ', 'PctEmploy ', 'PctEmplManu ', 'PctEmplProfServ ', 'PctOccupManu ', 'PctOccupMgmtProf ', 'MalePctDivorce ', 'MalePctNevMarr ', 'FemalePctDiv ', 'TotalPctDiv ', 'PersPerFam ', 'PctFam2Par ', 'PctKids2Par ', 'PctYoungKids2Par ', 'PctTeen2Par' , 'PctWorkMomYoungKids', 'PctWorkMom ', 'NumIlleg ', 'PctIlleg ', 'NumImmig ', 'PctImmigRecent ', 'PctImmigRec5 ', 'PctImmigRec8 ', 'PctImmigRec10 ', 'PctRecentImmig ', 'PctRecImmig5 ', 'PctRecImmig8 ', 'PctRecImmig10 ', 'PctSpeakEnglOnly ', 'PctNotSpeakEnglWell ', 'PctLargHouseFam ', 'PctLargHouseOccup ', 'PersPerOccupHous ', 'PersPerOwnOccHous ', 'PersPerRentOccHous ', 'PctPersOwnOccup ', 'PctPersDenseHous ', 'PctHousLess3BR ', 'MedNumBR ', 'HousVacant ','PctHousOccup ', 'PctHousOwnOcc ', 'PctVacantBoarded ', 'PctVacMore6Mos ', 'MedYrHousBuilt ', 'PctHousNoPhone ', 'PctWOFullPlumb ', 'OwnOccLowQuart ', 'OwnOccMedVal ', 'OwnOccHiQuart ', 'RentLowQ ', 'RentMedian ', 'RentHighQ ', 'MedRent ', 'MedRentPctHousInc ', 'MedOwnCostPctInc ', 'MedOwnCostPctIncNoMtg ', 'NumInShelters ', 'NumStreet ', 'PctForeignBorn ', 'PctBornSameState ', 'PctSameHouse85 ', 'PctSameCity85 ', 'PctSameState85 ', 'LemasSwornFT ', 'LemasSwFTPerPop ', 'LemasSwFTFieldOps ', 'LemasSwFTFieldPerPop ', 'LemasTotalReq ', 'LemasTotReqPerPop ', 'PolicReqPerOffic ', 'PolicPerPop ', 'RacialMatchCommPol ', 'PctPolicWhite ', 'PctPolicBlack ', 'PctPolicHisp ', 'PctPolicAsian ', 'PctPolicMinor ', 'OfficAssgnDrugUnits ', 'NumKindsDrugsSeiz ', 'PolicAveOTWorked ', 'LandArea ', 'PopDens ', 'PctUsePubTrans ', 'PolicCars ', 'PolicOperBudg','LemasPctPolicOnPatr ','LemasGangUnitDeploy ','LemasPctOfficDrugUn','PolicBudgPerPop','ViolentCrimesPerPop'};

%% getting rid of missing data
pp=isnan(dt);
qq=sum(pp);
vec=find(qq>0);
ind=find(pp(:,vec(1))>0);
dt(ind,:)=[];

pp=isnan(dt);
qq=sum(pp);
vec=find(qq>0);
dt(:,vec)=[];
header(vec)=[];

crime=dataset({dt,header{:}});


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

n=size(crime,1);
p=size(crime,2)-1;

n_simu=5;
svr_test_error_nrfe=zeros(n_simu,1);  
svr_train_error_nrfe=zeros(n_simu,1);
svr_test_error_rfe=zeros(n_simu,1);  
svr_train_error_rfe=zeros(n_simu,1);
mes_err=zeros(n_simu,1);
mes_err_train=zeros(n_simu,1);
vars_sel_rfe=cell(n_simu,1);
vars_sel_lasso=cell(n_simu,1);

seeds=zeros(n_simu,1);


for i=1:n_simu

    
seeds(i)=rand(1)*i;    
RandStream.setGlobalStream(RandStream('mcg16807', 'Seed',seeds(i)))    
    
    
%% test and training data

perm=(randperm(size(crime,1)));
train_ind=perm(1:ceil(end/2));
test_ind=perm(ceil(end/2)+1:end);
crime_train=crime(train_ind,:);
crime_test=crime(test_ind,:);


%% non linear regression

loss_name='epsilon_insensitive_loss';
loss_gridsel='linear_loss';
algo=svr(kernel('rbf'));
algo.optimizer='andre';
dcrime=crime_train;
dcrime_test=crime_test;
d=data('d',double(dcrime(:,1:end-1)),dcrime.ViolentCrimesPerPop);
d_test=data('d_test',double(dcrime_test(:,1:end-1)),dcrime_test.ViolentCrimesPerPop);



% %motor score
% d=d_motor;
% d_test=d_motor_test;


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;

%% 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');
svr_test_error_nrfe(i)=svr_test_nrfe.Y;

svr_train_nrfe=test(svr_object,d,'quadratic_loss');
svr_train_error_nrfe(i)=svr_train_nrfe.Y;

%% SVR with 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=rfe_object.val; ranks=rfe_object.rank;
svr_object_rfe=algo; svr_object_rfe.loss=loss_name;
train_X=d.X;

% Test and training error for RFE using change point model
dat=double(train_X);
y=d.Y;
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);
ranks_rest=setdiff(1:size(dat,2),ranks);
vars_sel_rfe{i}=[ranks_rest,ranks(1:end-j_imp)];

train_X_rfe=train_X(:,vars_sel_rfe{i});
train_Y=d.Y;
d_rfe_train=data(train_X_rfe,train_Y);
[d_rfe_trained,svr_object_rfe]=train(svr_object_rfe,d_rfe_train);
% misclassification rates
test_X=d_test.X;
test_X_rfe=test_X(:,vars_sel_rfe{i});
test_Y=d_test.Y;
d_rfe_test=data(test_X_rfe,test_Y);
svr_test=test(svr_object_rfe,d_rfe_test,'quadratic_loss');
svr_test_error_rfe(i)=svr_test.Y;
svr_train_rfe=test(svr_object_rfe,d_rfe_train,'quadratic_loss');
svr_train_error_rfe(i)=svr_train_rfe.Y;
save(filename3);


save(filename3);

Cov=get_x(d); Y=get_y(d);
model_lasso=glmnet(Cov, Y);
cv=cvglmnet(Cov,Y,5,[],'response','gaussian',glmnetSet,0);
lam_opt=cv.lambda_min;
lasso_fit = glmnet(Cov, Y); % regression

beta= glmnetCoef(lasso_fit,lam_opt);
lasso_set=find(beta(2:end)~=0);
vars_sel_lasso{i}=lasso_set;

%ranks_corr_reg(:,i)=lasso_set(1:imp_dim);
Cov_lasso=Cov(:,lasso_set);
beta_fin = regress(Y,[ones(size(Cov_lasso,1),1),Cov_lasso]);
Cov_lasso_test=d_test.X(:,lasso_set);
yfit = [ones(size(Cov_lasso_test,1),1),Cov_lasso_test]*beta_fin;
yresid = d_test.Y - yfit;
mes_err(i) = mean(yresid.^2);
yfit_train = [ones(size(Cov_lasso,1),1),Cov_lasso]*beta_fin;
yresid_train = d.Y - yfit_train;
mes_err_train(i) = mean(yresid_train.^2);


end

save(filename3);


%% total training error
untrained=algo;
sel_data=d;
numEx=size(d.X,1);
lss=loss(loss_name);
lss.param=0.1;
[res,trained]=train(untrained,sel_data);
res=test(trained,sel_data);
r=train(lss,res);
%r=test(trained,sel_data,lss);
alpha = trained.alpha;
kernel_sv=train(trained.child,trained.Xsv);
Ct=0.5*alpha'*kernel_sv.X*alpha/numEx+get_y(r)*trained.C;

save(filename3);