forked from omerdr/ensemble-regression
-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathsingle_run.m
256 lines (234 loc) · 13.2 KB
/
single_run.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
% Comparing approaches for ensemble regression
%% Init
clear all; close all;
clc;
addpath 'Ensemble_Regressors'
addpath 'DataGeneration'
[~,hostname] = system('hostname'); hostname = strtrim(hostname);
random_init_iterations = 10;
if ~isempty(strfind(hostname,'weiz')); delete(gcp); parpool(24); random_init_iterations = 1000; end; % if in WIS run worker pool
%%%%%%%% LOAD DATA SET BEFORE RUNNING THIS CODE %%%%%%%%%%%%%%%
%%%
%%% The code assumes existance of Z,y,Ztrain,ytrain, Ey, Ey2
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ROOT = '~/code/github/ensemble-regression/auto/';
%ROOT = './Datasets/auto_mlp5/';
%ROOT = './Datasets/auto/';
ROOT = './Datasets/final/misc';
file_list = dir([ROOT '*.mat']);
datasets = cell(1,length(file_list));
for i=1:length(file_list)
datasets{i}=[ROOT file_list(i).name];
end;
% datasets = {'./Datasets/auto_mlp5/auto_basketball.mat'};
% datasets = {'~/code/github/ensemble-regression/auto/auto_ccpp.mat'};
%% For each dataset
for dataset_name=datasets
%%
fprintf([dataset_name{1} '\n']);
load(dataset_name{1})
if isempty(strfind(dataset_name{1},'MVGaussianDepData'))
y_true = double(y); clear y; % renmae y to y_true and make sure both y and y_true are double
ytrain = double(ytrain); % (in the sweets dataset y is integer)
end;
[m,n] = size(Z);
%if n < 10000; continue; end;
n_training_set = size(Ztrain,2);
var_y = Ey2 - Ey.^2;
%% Predict
% Oracle Prediction
[y_oracle, beta_oracle] = ER_linear_regression_oracle(y_true, Z);
% Perrone & Cooper Supervised GEM
[y_pc, w_pc, C_pc] = ER_PerroneCooperGEM(Ztrain, ytrain, Z);
% Breiman Star
[y_breiman, w_breiman, fval_breiman] = ER_BreimanStar(Ztrain, ytrain, Z);
%% 2nd Moment Estimator (assumes Ey, Ey^2 are given)
Ey = mean(y_true); Ey2 = mean(y_true.^2);
[y_mean, beta_mean] = ER_MeanWithBiasCorrection(Z,Ey);
[y_varw, beta_varw] = ER_VarianceWeightingWithBiasCorrection(Z,Ey);
[y_invc, beta_invc] = ER_InverseCov(Z,Ey);
[y_2, beta_2] = ER_SecondMoment(Z,Ey,Ey2);
[y_2oracle, beta_2oracle] = ER_SecondMomentOracle(Z,Ey,Ey2,y_true);
[y_cvx_m, beta_cvx_m] = ER_ConvexProgram(Z, Ey, Ey2, 2, 2, beta_mean);
[y_cvx_v, beta_cvx_v] = ER_ConvexProgram(Z, Ey, Ey2, 2, 2, beta_varw);
[y_cvx_2, beta_cvx_2] = ER_ConvexProgram(Z, Ey, Ey2, 2, 2, beta_2);
[y_cvx_cv_m, beta_cvx_cv_m, lambda_mean_m, lambda_var_m, fval_m, mean_reg_mse_m] = ...
ER_ConvexProgramWithCV(Z, Ey, Ey2, beta_mean); %set(gcf,'Name',[dataset_name{1} ' (init mean)']);
[y_cvx_cv_v, beta_cvx_cv_v, lambda_mean_v, lambda_var_v, fval_v, mean_reg_mse_v] = ...
ER_ConvexProgramWithCV(Z, Ey, Ey2, beta_varw); %set(gcf,'Name',[dataset_name{1} ' (init varw)']);
[y_cvx_cv_2, beta_cvx_cv_2, lambda_mean_2, lambda_var_2, fval_2, mean_reg_mse_2] = ...
ER_ConvexProgramWithCV(Z, Ey, Ey2, beta_2); %set(gcf,'Name',[dataset_name{1} ' (init 2ME)']);
%[y_cvx2, beta_cvx2] = ER_ConvexProgram2(Z, Ey, Ey2, lambda_mean, lambda_var);
%[y_cvx_cv2, beta_cvx_cv2, lambda_mean2, lambda_var2] = ER_ConvexProgramWithCV2(Z, Ey, Ey2);
%[y_b, beta_b] = ER_Boaz(Z,Ey,Ey2);
%[y_qp, beta_qp] = ER_QuadProg(Z,Ey,Ey2);
%% Test random init
% ER_ConvexProgram INTERNALS
% Init params
b_hat = mean(Z,2) - Ey; % approximate bias
C = cov(Z');
Zc = Z - repmat(b_hat,1,n);
Z1c = [ones(1,n); Zc];
K1 = Zc*Z1c';
K2 = Z1c*Z1c';
% constraints w_i >= 0
A = [zeros(m,1), -eye(m)];
b = zeros(m,1);
options = optimoptions('fmincon', 'Display','None');%'iter-detailed');%,'MaxFunEvals',1e4);
%
fval_cur = zeros(random_init_iterations,1); lambda_rnd_cur = ceil(rand(random_init_iterations,1)*8)/2; %zeros(iters,1); % lambda in .5:.5:4
mean_reg_mse_cur = zeros(random_init_iterations,1);
beta_0_cur = cell(random_init_iterations,1); beta_cur = cell(random_init_iterations,1); y_cur = cell(random_init_iterations,1);
system('rm -f progress.txt; touch progress.txt');
parfor j=1:random_init_iterations % average 50 seconds per 100 iterations
beta_0_cur{j} = rand(m+1,1);
%lambda_rnd_cur(j) = ceil(rand*8)/2; % lambda in .5:.5:4
obj = @(beta) cvx_opt_for_w(beta,m,n,Ey,var_y, C, Z1c, K1, K2, 0, lambda_rnd_cur(j));
%[y_cur{j}, beta_cur{j}, fval_cur(j)] = ER_ConvexProgram(Z, Ey, Ey2, 0, lambda_rnd_cur(j), beta_0_cur{j});
[beta_cur{j},fval_cur(j)] = fmincon(obj, beta_0_cur{j}, A, b,[],[],[],[], [], options); % ER_ConvexProgram INTERNALS
mean_reg_mse_cur(j) = mean(mean((Zc - repmat(beta_cur{j}'*Z1c,m,1)).^2,2)); % the mean regressor error
system('echo "." >> progress.txt');
end;
j = find(fval_cur == min(fval_cur));
beta_cvx_cv_0_rnd = beta_0_cur{j};
beta_cvx_cv_rnd = beta_cur{j};
fval_rnd = fval_cur(j);
lambda_var_rnd = lambda_rnd_cur(j);
%y_cvx_cv_rnd=y_cur{j};
y_cvx_cv_rnd = Z1c'*beta_cvx_cv_rnd; % ER_ConvexProgram INTERNALS
mean_reg_mse_rnd = mean_reg_mse_cur(j);
%% Print Results
fprintf('\nDataset: %s, Var(y) = %.4g, n = %d, m = %d\n',dataset_name{1},var_y,n,m);
fprintf('=====================================================================================\n');
fprintf('Init Mean:\tLambda mean = %.2f, lambda var = %.2f, fval = %.4g, mean epsilon_i^2 = %.4g, var = %.4g\n', lambda_mean_m, lambda_var_m, fval_m, mean_reg_mse_m, var(y_cvx_cv_m));
fprintf('Init Var: \tLambda mean = %.2f, lambda var = %.2f, fval = %.4g, mean epsilon_i^2 = %.4g, var = %.4g\n', lambda_mean_v, lambda_var_v, fval_v, mean_reg_mse_v, var(y_cvx_cv_v));
fprintf('Init 2ME: \tLambda mean = %.2f, lambda var = %.2f, fval = %.4g, mean epsilon_i^2 = %.4g, var = %.4g\n', lambda_mean_2, lambda_var_2, fval_2, mean_reg_mse_2, var(y_cvx_cv_2));
fprintf('Init RND: \tLambda mean = %.2f, lambda var = %.2f, fval = %.4g, mean epsilon_i^2 = %.4g, var = %.4g\n', 0, lambda_var_rnd, fval_rnd, mean_reg_mse_rnd, var(y_cvx_cv_rnd));
results = [y_cvx_cv_m, y_cvx_cv_v, y_cvx_cv_2, y_cvx_cv_rnd];
fvals = [fval_m fval_v fval_2 fval_rnd];
y_our_method = results(:,find(fvals == min(fvals)));
mse_oracle = mean((y_true'-y_oracle).^2);
mse = @(x) (mean((y_true'-x).^2)) / mse_oracle; % show MSE relative to oracle.
calc_g = @(x) (mean((Zc - repmat(x',m,1)).*repmat(x',m,1),2)); % calc g assuming x is the true response (y_true)
calc_g_rel = @(x) (calc_g(x) ./ calc_g(y_oracle));
[y_2g, beta_2g] = ER_SecondMomentGiven_g(Z,Ey,Ey2, calc_g(y_our_method));
results_summary_current_dataset = ...
{'oracle', mse_oracle ; ...
'e2ME Oracle', mse(y_2oracle) ; ...
'PC', mse(y_pc) ; ...
'Breiman Star', mse(y_breiman) ; ...
'Mean', mse(y_mean) ; ...
'Var Weighted', mse(y_varw) ; ...
'Inv Cov', mse(y_invc) ; ...
'e2ME', mse(y_2) ; ...
'g 2ME', mse(y_2g) ; ...
'Our Method', mse(y_our_method) ; ...
'CVX init mean', mse(y_cvx_m) ; ...
'CVX init varw', mse(y_cvx_v) ; ...
'CVX init 2ME', mse(y_cvx_2) ; ...
'CVX with CV init mean', mse(y_cvx_cv_m); ...
'CVX with CV init varw', mse(y_cvx_cv_v); ...
'CVX with CV init 2ME', mse(y_cvx_cv_2); ...
'CVX with CV init RND', mse(y_cvx_cv_rnd); ...
'fval_mean', fval_m; ...
'fval_varw', fval_v; ...
'fval_2ME', fval_2; ...
'fval_rnd', fval_rnd; ...
'n', n; 'm', m; 'Ey',Ey; 'Var_y',var_y; ...
'Lambda init mean', lambda_var_m; ...
'Lambda init varw', lambda_var_v; ...
'Lambda init 2ME', lambda_var_2; ...
'Lambda init RND', lambda_var_rnd; ...
'Mean Regressor MSE init mean', mean_reg_mse_m; ...
'Mean Regressor MSE init varw', mean_reg_mse_v; ...
'Mean Regressor MSE init 2ME', mean_reg_mse_2; ...
'Mean Regressor MSE init RND', mean_reg_mse_rnd; ...
'cond Cov Z transpose', cond(cov(Z')); ...
... % CALC mean(g) over m regressors
'oracle mean g', mean(calc_g(y_oracle)) ; ...
'e2ME Oracle mean g', mean(calc_g(y_2oracle)) ; ...
'PC mean g', mean(calc_g(y_pc)) ; ...
'Breiman Star mean g', mean(calc_g(y_breiman)) ; ...
'Mean mean g', mean(calc_g(y_mean)) ; ...
'Var Weighted mean g', mean(calc_g(y_varw)) ; ...
'Inv Cov mean g', mean(calc_g(y_invc)) ; ...
'e2ME mean g', mean(calc_g(y_2)) ; ...
'g 2ME mean g', mean(calc_g(y_2g)) ; ...
'Our Method mean g', mean(calc_g(y_our_method)) ; ...
'CVX mean mean g', mean(calc_g(y_cvx_m)) ; ...
'CVX varw mean g', mean(calc_g(y_cvx_v)) ; ...
'CVX 2ME mean g', mean(calc_g(y_cvx_2)) ; ...
'CVX CV iMean mean g', mean(calc_g(y_cvx_cv_m)); ...
'CVX CV iVarw mean g', mean(calc_g(y_cvx_cv_v)); ...
'CVX CV i2ME mean g', mean(calc_g(y_cvx_cv_2)); ...
'CVX CV iRND mean g', mean(calc_g(y_cvx_cv_rnd)); ...
... % CALC var(g)
'oracle var g', var(calc_g(y_oracle)) ; ...
'e2ME Oracle var g', var(calc_g(y_2oracle)) ; ...
'PC var g', var(calc_g(y_pc)) ; ...
'Breiman var g', var(calc_g(y_breiman)) ; ...
'Mean var g', var(calc_g(y_mean)) ; ...
'Var Weighted var g', var(calc_g(y_varw)) ; ...
'Inv Cov var g', var(calc_g(y_invc)) ; ...
'e2ME var g', var(calc_g(y_2)) ; ...
'g 2ME var g', var(calc_g(y_2g)) ; ...
'Our Method var g', var(calc_g(y_our_method)) ; ...
'CVX mean var g', var(calc_g(y_cvx_m)) ; ...
'CVX varw var g', var(calc_g(y_cvx_v)) ; ...
'CVX 2ME var g', var(calc_g(y_cvx_2)) ; ...
'CVX CV iMean var g', var(calc_g(y_cvx_cv_m)); ...
'CVX CV iVarw var g', var(calc_g(y_cvx_cv_v)); ...
'CVX CV i2ME var g', var(calc_g(y_cvx_cv_2)); ...
'CVX CV iRND var g', var(calc_g(y_cvx_cv_rnd)); ...
... % PRINT g
'oracle g', num2str(calc_g(y_oracle)','%g\t') ; ...
'e2ME Oracle g', num2str(calc_g_rel(y_2oracle)','%g\t') ; ...
'PC g', num2str(calc_g_rel(y_pc)','%g\t') ; ...
'Breiman Star g', num2str(calc_g_rel(y_breiman)','%g\t') ; ...
'Mean g', num2str(calc_g_rel(y_mean)','%g\t') ; ...
'Var Weighted g', num2str(calc_g_rel(y_varw)','%g\t') ; ...
'Inv Cov g', num2str(calc_g_rel(y_invc)','%g\t') ; ...
'e2ME g', num2str(calc_g_rel(y_2)','%g\t') ; ...
'g 2ME g', num2str(calc_g_rel(y_2g)','%g\t') ; ...
'Our Method g', num2str(calc_g_rel(y_our_method)','%g\t') ; ...
'CVX mean g', num2str(calc_g_rel(y_cvx_m)','%g\t') ; ...
'CVX varw g', num2str(calc_g_rel(y_cvx_v)','%g\t') ; ...
'CVX 2ME g', num2str(calc_g_rel(y_cvx_2)','%g\t') ; ...
'CVX CV iMean g', num2str(calc_g_rel(y_cvx_cv_m)','%g\t') ; ...
'CVX CV iVarw g', num2str(calc_g_rel(y_cvx_cv_v)','%g\t') ; ...
'CVX CV i2ME g', num2str(calc_g_rel(y_cvx_cv_2)','%g\t') ; ...
'CVX CV iRND g', num2str(calc_g_rel(y_cvx_cv_rnd)','%g\t') ; ...
};
%% print g's
% fprintf('### g estimation\n');
% fprintf('### ============\n');
% fprintf('### oracle:\t %s\n', num2str(calc_g(y_oracle)'/var_y));
% fprintf('### e2ME Oracle:\t %s\n', num2str(calc_g(y_2oracle)'/var_y));
% fprintf('### PC:\t %s\n', num2str(calc_g(y_pc)'/var_y));
% fprintf('### Mean:\t %s\n', num2str(calc_g(y_mean)'/var_y));
% fprintf('### VarW:\t %s\n', num2str(calc_g(y_varw)'/var_y));
% fprintf('### e2ME:\t %s\n', num2str(calc_g(y_2)'/var_y));
% fprintf('### g 2ME:\t %s\n', num2str(calc_g(y_2g)'/var_y));
% fprintf('### Our Method:\t %s\n', num2str(calc_g(y_our_method)'/var_y));
% fprintf('### CVX mean:\t %s\n', num2str(calc_g(y_cvx_m)'/var_y));
% fprintf('### CVX varw:\t %s\n', num2str(calc_g(y_cvx_v)'/var_y));
% fprintf('### CVX 2ME:\t %s\n', num2str(calc_g(y_cvx_2)'/var_y));
% fprintf('### CVX CV iMean:\t %s\n', num2str(calc_g(y_cvx_cv_m)'/var_y));
% fprintf('### CVX CV iVarw:\t %s\n', num2str(calc_g(y_cvx_cv_v)'/var_y));
% fprintf('### CVX CV i2ME:\t %s\n', num2str(calc_g(y_cvx_cv_2)'/var_y));
% fprintf('### CVX CV iRND:\t %s\n', num2str(calc_g(y_cvx_cv_rnd)'/var_y));
if ~exist('results_summary','var')
results_summary = cell(length(datasets), size(results_summary_current_dataset,1));
end;
for i=1:length(results_summary_current_dataset)
fprintf('%25s \t%g\n',results_summary_current_dataset{i,1}, results_summary_current_dataset{i,2});
results_summary{find(strcmp(dataset_name,datasets)),i} = results_summary_current_dataset{i,2};
end;
%%
end % for dataset_name in datasets
%% Print Result Summary
cols = {results_summary_current_dataset{:,1}};
for i=1:length(cols); cols{i}=strrep(cols{i},' ','_');end;
out=cell2table(results_summary, 'RowNames', datasets', 'VariableNames', cols)
writetable(out, ['results/last_summary-' hostname '.csv'],'WriteRowNames',true)