1+ import numpy as np
2+ import sys
3+ # import warnings
4+ # warnings.filterwarnings("ignore")
5+ from Orange .classification .fasterrisk .utils import normalize_X , compute_logisticLoss_from_ExpyXB
6+
7+ class logRegModel :
8+ def __init__ (self , X , y , lambda2 = 1e-8 , intercept = True , original_lb = - 5 , original_ub = 5 ):
9+ self .X = X
10+ self .X_normalized , self .X_mean , self .X_norm , self .scaled_feature_indices = normalize_X (self .X )
11+ self .n , self .p = self .X_normalized .shape
12+ self .y = y .reshape (- 1 ).astype (float )
13+ self .yX = y .reshape (- 1 , 1 ) * self .X_normalized
14+ self .yXT = np .zeros ((self .p , self .n ))
15+ self .yXT [:] = np .transpose (self .yX )[:]
16+ self .beta0 = 0
17+ self .betas = np .zeros ((self .p , ))
18+ self .ExpyXB = np .exp (self .y * self .beta0 + self .yX .dot (self .betas ))
19+
20+ self .intercept = intercept
21+ self .lambda2 = lambda2
22+ self .twoLambda2 = 2 * self .lambda2
23+
24+ self .Lipschitz = 0.25 + self .twoLambda2
25+ self .lbs = original_lb * np .ones (self .p )
26+ self .lbs [self .scaled_feature_indices ] *= self .X_norm [self .scaled_feature_indices ]
27+ self .ubs = original_ub * np .ones (self .p )
28+ self .ubs [self .scaled_feature_indices ] *= self .X_norm [self .scaled_feature_indices ]
29+
30+ self .total_child_added = 0
31+
32+ def warm_start_from_original_beta0_betas (self , original_beta0 , original_betas ):
33+ # betas_initial has dimension (p+1, 1)
34+ self .original_beta0 = original_beta0
35+ self .original_betas = original_betas
36+ self .beta0 , self .betas = self .transform_coefficients_to_normalized_space (self .original_beta0 , self .original_betas )
37+ print ("warmstart solution in normalized space is {} and {}" .format (self .beta0 , self .betas ))
38+ self .ExpyXB = np .exp (self .y * self .beta0 + self .yX .dot (self .betas ))
39+
40+ def warm_start_from_beta0_betas (self , beta0 , betas ):
41+ self .beta0 , self .betas = beta0 , betas
42+ self .ExpyXB = np .exp (self .y * self .beta0 + self .yX .dot (self .betas ))
43+
44+ def warm_start_from_beta0_betas_ExpyXB (self , beta0 , betas , ExpyXB ):
45+ self .beta0 , self .betas , self .ExpyXB = beta0 , betas , ExpyXB
46+
47+ def get_beta0_betas (self ):
48+ return self .beta0 , self .betas
49+
50+ def get_beta0_betas_ExpyXB (self ):
51+ return self .beta0 , self .betas , self .ExpyXB
52+
53+ def get_original_beta0_betas (self ):
54+ return self .transform_coefficients_to_original_space (self .beta0 , self .betas )
55+
56+ def transform_coefficients_to_original_space (self , beta0 , betas ):
57+ original_betas = betas .copy ()
58+ original_betas [self .scaled_feature_indices ] = original_betas [self .scaled_feature_indices ]/ self .X_norm [self .scaled_feature_indices ]
59+ original_beta0 = beta0 - np .dot (self .X_mean , original_betas )
60+ return original_beta0 , original_betas
61+
62+ def transform_coefficients_to_normalized_space (self , original_beta0 , original_betas ):
63+ betas = original_betas .copy ()
64+ betas [self .scaled_feature_indices ] = betas [self .scaled_feature_indices ] * self .X_norm [self .scaled_feature_indices ]
65+ beta0 = original_beta0 + self .X_mean .dot (original_betas )
66+ return beta0 , betas
67+
68+ def get_grad_at_coord (self , ExpyXB , betas_j , yX_j , j ):
69+ # return -np.dot(1/(1+ExpyXB), self.yX[:, j]) + self.twoLambda2 * betas_j
70+ # return -np.inner(1/(1+ExpyXB), self.yX[:, j]) + self.twoLambda2 * betas_j
71+ # return -np.inner(np.reciprocal(1+ExpyXB), self.yX[:, j]) + self.twoLambda2 * betas_j
72+ return - np .inner (np .reciprocal (1 + ExpyXB ), yX_j ) + self .twoLambda2 * betas_j
73+ # return -yX_j.dot(np.reciprocal(1+ExpyXB)) + self.twoLambda2 * betas_j
74+
75+ def update_ExpyXB (self , ExpyXB , yX_j , diff_betas_j ):
76+ ExpyXB *= np .exp (yX_j * diff_betas_j )
77+
78+ def optimize_1step_at_coord (self , ExpyXB , betas , yX_j , j ):
79+ # in-place modification, heck that ExpyXB and betas are passed by reference
80+ prev_betas_j = betas [j ]
81+ current_betas_j = prev_betas_j
82+ grad_at_j = self .get_grad_at_coord (ExpyXB , current_betas_j , yX_j , j )
83+ step_at_j = grad_at_j / self .Lipschitz
84+ current_betas_j = prev_betas_j - step_at_j
85+ # current_betas_j = np.clip(current_betas_j, self.lbs[j], self.ubs[j])
86+ current_betas_j = max (self .lbs [j ], min (self .ubs [j ], current_betas_j ))
87+ diff_betas_j = current_betas_j - prev_betas_j
88+ betas [j ] = current_betas_j
89+
90+ # ExpyXB *= np.exp(yX_j * diff_betas_j)
91+ self .update_ExpyXB (ExpyXB , yX_j , diff_betas_j )
92+
93+ def finetune_on_current_support (self , ExpyXB , beta0 , betas , total_CD_steps = 100 ):
94+
95+ support = np .where (np .abs (betas ) > 1e-9 )[0 ]
96+ grad_on_support = - self .yXT [support ].dot (np .reciprocal (1 + ExpyXB )) + self .twoLambda2 * betas [support ]
97+ abs_grad_on_support = np .abs (grad_on_support )
98+ support = support [np .argsort (- abs_grad_on_support )]
99+
100+ loss_before = compute_logisticLoss_from_ExpyXB (ExpyXB ) + self .lambda2 * betas [support ].dot (betas [support ])
101+ for steps in range (total_CD_steps ): # number of iterations for coordinate descent
102+
103+ if self .intercept :
104+ grad_intercept = - np .reciprocal (1 + ExpyXB ).dot (self .y )
105+ step_at_intercept = grad_intercept / (self .n * 0.25 ) # lipschitz constant is 0.25 at the intercept
106+ beta0 = beta0 - step_at_intercept
107+ ExpyXB *= np .exp (self .y * (- step_at_intercept ))
108+
109+ for j in support :
110+ self .optimize_1step_at_coord (ExpyXB , betas , self .yXT [j , :], j ) # in-place modification on ExpyXB and betas
111+
112+ if steps % 10 == 0 :
113+ loss_after = compute_logisticLoss_from_ExpyXB (ExpyXB ) + self .lambda2 * betas [support ].dot (betas [support ])
114+ if abs (loss_before - loss_after )/ loss_after < 1e-8 :
115+ # print("break after {} steps; support size is {}".format(steps, len(support)))
116+ break
117+ loss_before = loss_after
118+
119+ return ExpyXB , beta0 , betas
120+
121+ def compute_yXB (self , beta0 , betas ):
122+ return self .y * (beta0 + np .dot (self .X_normalized , betas ))
123+
0 commit comments