Skip to content

Commit

Permalink
fix: bug when support_set is range(k)
Browse files Browse the repository at this point in the history 
del 2 lines code:
        // If A_U not change, U will not change and we can stop.
        if (A_U.size() == 0 || A_U.maxCoeff() == T0 - 1)
            break;
  • Loading branch information
@bbayukari
bbayukari committed Apr 25, 2024
1 parent 4764108 commit 450f023
Show file tree
Hide file tree
Showing 4 changed files with 152 additions and 120 deletions.
23 changes: 15 additions & 8 deletions src/Algorithm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,8 +241,8 @@ void Algorithm::get_A(UniversalData &X, MatrixXd &y, VectorXi &A, VectorXi &I, i
}

// If A_U not change, U will not change and we can stop.
if (A_U.size() == 0 || A_U.maxCoeff() == T0 - 1)
break;
// if (A_U.size() == 0 || A_U.maxCoeff() == T0 - 1)
// break;

// Update & Restore beta, A from U
slice_restore(beta_U, U_ind, beta);
Expand Down Expand Up @@ -320,10 +320,10 @@ bool Algorithm::splicing(UniversalData &X, MatrixXd &y, VectorXi &A, VectorXi &I

for (int k = C_max; k >= 1;)
{
//SPDLOG_INFO("exchange num is {}", k);
// SPDLOG_INFO("exchange num is {}", k);
A_exchange = diff_union(A, s1, s2);
A_ind_exchage = find_ind(A_exchange, g_index, g_size, (this->beta).rows(), N);
X_A_exchage = X.slice_by_para(A_ind_exchage);
X_A_exchage = X.slice_by_para(A_ind_exchage);
beta_A_exchange = beta(A_ind_exchage);
coef0_A_exchange = coef0;
bool success = this->primary_model_fit(X_A_exchage, y, weights, beta_A_exchange, coef0_A_exchange,
Expand All @@ -346,7 +346,7 @@ bool Algorithm::splicing(UniversalData &X, MatrixXd &y, VectorXi &A, VectorXi &I
s2 = s2.head(k).eval();
}

if (train_loss - best_loss <= tau)
if (train_loss <= best_loss)
return false;

train_loss = best_loss;
Expand All @@ -356,7 +356,7 @@ bool Algorithm::splicing(UniversalData &X, MatrixXd &y, VectorXi &A, VectorXi &I
coef0 = best_coef0_A_exchange;
if (this->is_dynamic_exchange_num)
C_max = best_exchange_num;
//SPDLOG_INFO("best exchange num is {}", best_exchange_num);
// SPDLOG_INFO("best exchange num is {}", best_exchange_num);
return true;
};

Expand All @@ -369,6 +369,10 @@ VectorXi Algorithm::inital_screening(UniversalData &X, MatrixXd &y, VectorXd &be
// variable initialization
int beta_size = X.cols();
bd = VectorXd::Zero(N);
if (A.size() == 0)
{
A = max_k(beta.cwiseAbs(), this->sparsity_level);
}

// calculate beta & d & h
VectorXi A_ind = find_ind(A, g_index, g_size, beta_size, N);
Expand Down Expand Up @@ -452,12 +456,15 @@ bool Algorithm::primary_model_fit(UniversalData &active_data, MatrixXd &y, Vecto
void Algorithm::sacrifice(UniversalData &data, UniversalData &XA, MatrixXd &y, VectorXd &para, VectorXd &beta_A, VectorXd &aux_para, VectorXi &A, VectorXi &I, VectorXd &weights, VectorXi &g_index, VectorXi &g_size, int g_num, VectorXi &A_ind, VectorXd &sacrifice, VectorXi &U, VectorXi &U_ind, int num)
{
SPDLOG_DEBUG("sacrifice begin");
SPDLOG_DEBUG("active set is {}", A.transpose());
VectorXd gradient_full;
MatrixXd hessian_full;
if (this->use_hessian){
if (this->use_hessian)
{
data.gradient_and_hessian(para, gradient_full, hessian_full);
}
else{
else
{
data.loss_and_gradient(para, gradient_full);
hessian_full = MatrixXd::Identity(para.size(), para.size());
}
Expand Down
10 changes: 1 addition & 9 deletions src/Algorithm.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,15 +115,7 @@ class Algorithm

void update_tau(int n, int p)
{
if (n == 1)
{
this->tau = 0.0;
}
else
{
this->tau =
0.01 * (double)this->sparsity_level * log((double)p) * log(log((double)n)) / (double)n;
}
this->tau = 0.0;
}

bool get_warm_start() { return this->warm_start; }
Expand Down
158 changes: 97 additions & 61 deletions src/utilities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,8 @@
* Licensed under the MIT License.
*/


#include <Eigen/Eigen>



#include <string.h>

#include <algorithm>
Expand All @@ -22,14 +19,19 @@ using namespace Eigen;

std::mt19937 GLOBAL_RNG(1);

Eigen::VectorXi find_ind(Eigen::VectorXi &L, Eigen::VectorXi &gindex, Eigen::VectorXi &gsize, int beta_size, int N) {
if (L.size() == N) {
Eigen::VectorXi find_ind(Eigen::VectorXi &L, Eigen::VectorXi &gindex, Eigen::VectorXi &gsize, int beta_size, int N)
{
if (L.size() == N)
{
return Eigen::VectorXi::LinSpaced(beta_size, 0, beta_size - 1);
} else {
}
else
{
int mark = 0;
Eigen::VectorXi ind = Eigen::VectorXi::Zero(beta_size);

for (int i = 0; i < L.size(); i++) {
for (int i = 0; i < L.size(); i++)
{
ind.segment(mark, gsize(L(i))) =
Eigen::VectorXi::LinSpaced(gsize(L(i)), gindex(L(i)), gindex(L(i)) + gsize(L(i)) - 1);
mark = mark + gsize(L(i));
Expand All @@ -38,15 +40,17 @@ Eigen::VectorXi find_ind(Eigen::VectorXi &L, Eigen::VectorXi &gindex, Eigen::Vec
}
}

UniversalData X_seg(UniversalData& X, int n, Eigen::VectorXi& ind, int model_type){
UniversalData X_seg(UniversalData &X, int n, Eigen::VectorXi &ind, int model_type)
{
return X.slice_by_para(ind);
}



void slice_assignment(Eigen::VectorXd &nums, Eigen::VectorXi &ind, double value) {
if (ind.size() != 0) {
for (int i = 0; i < ind.size(); i++) {
void slice_assignment(Eigen::VectorXd &nums, Eigen::VectorXi &ind, double value)
{
if (ind.size() != 0)
{
for (int i = 0; i < ind.size(); i++)
{
nums(ind(i)) = value;
}
}
Expand All @@ -55,11 +59,15 @@ void slice_assignment(Eigen::VectorXd &nums, Eigen::VectorXi &ind, double value)

// replace B by C in A
// to do : binary search
Eigen::VectorXi diff_union(Eigen::VectorXi A, Eigen::VectorXi &B, Eigen::VectorXi &C) {
Eigen::VectorXi diff_union(Eigen::VectorXi A, Eigen::VectorXi &B, Eigen::VectorXi &C)
{
unsigned int k;
for (unsigned int i = 0; i < B.size(); i++) {
for (k = 0; k < A.size(); k++) {
if (B(i) == A(k)) {
for (unsigned int i = 0; i < B.size(); i++)
{
for (k = 0; k < A.size(); k++)
{
if (B(i) == A(k))
{
A(k) = C(i);
break;
}
Expand All @@ -69,32 +77,42 @@ Eigen::VectorXi diff_union(Eigen::VectorXi A, Eigen::VectorXi &B, Eigen::VectorX
return A;
}

Eigen::VectorXi min_k(Eigen::VectorXd &vec, int k, bool sort_by_value) {
Eigen::VectorXi ind = Eigen::VectorXi::LinSpaced(vec.size(), 0, vec.size() - 1); // [0 1 2 3 ... N-1]
Eigen::VectorXi min_k(Eigen::VectorXd &vec, int k, bool sort_by_value)
{
Eigen::VectorXi ind = Eigen::VectorXi::LinSpaced(vec.size(), 0, vec.size() - 1); // [0 1 2 3 ... N-1]
// shuffle index to avoid repeat results when there are several equal values in vec
std::shuffle(ind.data(), ind.data() + ind.size(), GLOBAL_RNG);

auto rule = [vec](int i, int j) -> bool { return vec(i) < vec(j); }; // sort rule
auto rule = [vec](int i, int j) -> bool
{ return vec(i) < vec(j); }; // sort rule
std::nth_element(ind.data(), ind.data() + k, ind.data() + ind.size(), rule);
if (sort_by_value) {
if (sort_by_value)
{
std::sort(ind.data(), ind.data() + k, rule);
} else {
}
else
{
std::sort(ind.data(), ind.data() + k);
}

return ind.head(k).eval();
}

Eigen::VectorXi max_k(Eigen::VectorXd &vec, int k, bool sort_by_value) {
Eigen::VectorXi ind = Eigen::VectorXi::LinSpaced(vec.size(), 0, vec.size() - 1); // [0 1 2 3 ... N-1]
Eigen::VectorXi max_k(const Eigen::VectorXd &vec, int k, bool sort_by_value)
{
Eigen::VectorXi ind = Eigen::VectorXi::LinSpaced(vec.size(), 0, vec.size() - 1); // [0 1 2 3 ... N-1]
// shuffle index to avoid repeat results when there are several equal values in vec
std::shuffle(ind.data(), ind.data() + ind.size(), GLOBAL_RNG);

auto rule = [vec](int i, int j) -> bool { return vec(i) > vec(j); }; // sort rule

auto rule = [vec](int i, int j) -> bool
{ return vec(i) > vec(j); }; // sort rule
std::nth_element(ind.data(), ind.data() + k, ind.data() + ind.size(), rule);
if (sort_by_value) {
if (sort_by_value)
{
std::sort(ind.data(), ind.data() + k, rule);
} else {
}
else
{
std::sort(ind.data(), ind.data() + k);
}
return ind.head(k).eval();
Expand All @@ -113,27 +131,38 @@ Eigen::VectorXi max_k(Eigen::VectorXd &vec, int k, bool sort_by_value) {
// }

// complement
Eigen::VectorXi complement(Eigen::VectorXi &A, int N) {
Eigen::VectorXi complement(Eigen::VectorXi &A, int N)
{
int A_size = A.size();
if (A_size == 0) {
if (A_size == 0)
{
return Eigen::VectorXi::LinSpaced(N, 0, N - 1);
} else if (A_size == N) {
}
else if (A_size == N)
{
Eigen::VectorXi I(0);
return I;
} else {
}
else
{
Eigen::VectorXi I(N - A_size);
int cur_index = 0;
int A_index = 0;
for (int i = 0; i < N; i++) {
if (A_index >= A_size) {
for (int i = 0; i < N; i++)
{
if (A_index >= A_size)
{
I(cur_index) = i;
cur_index += 1;
continue;
}
if (i != A(A_index)) {
if (i != A(A_index))
{
I(cur_index) = i;
cur_index += 1;
} else {
}
else
{
A_index += 1;
}
}
Expand Down Expand Up @@ -177,52 +206,63 @@ Eigen::VectorXi complement(Eigen::VectorXi &A, int N) {
// }
// }

void slice(Eigen::VectorXd &nums, Eigen::VectorXi &ind, Eigen::VectorXd &A) {
if (ind.size() == 0) {
void slice(Eigen::VectorXd &nums, Eigen::VectorXi &ind, Eigen::VectorXd &A)
{
if (ind.size() == 0)
{
A = Eigen::VectorXd::Zero(0);
} else {
}
else
{
A = Eigen::VectorXd::Zero(ind.size());
for (int i = 0; i < ind.size(); i++) {
for (int i = 0; i < ind.size(); i++)
{
A(i) = nums(ind(i));
}
}
}

void slice(Eigen::MatrixXd &nums, Eigen::VectorXi &ind, Eigen::MatrixXd &A) {
void slice(Eigen::MatrixXd &nums, Eigen::VectorXi &ind, Eigen::MatrixXd &A)
{
A = Eigen::MatrixXd::Zero(ind.size(), nums.cols());
if (ind.size() != 0) {
for (int i = 0; i < ind.size(); i++) {
if (ind.size() != 0)
{
for (int i = 0; i < ind.size(); i++)
{
A.row(i) = nums.row(ind(i));
}
}
}


void slice(UniversalData& nums, Eigen::VectorXi& ind, UniversalData& A)
void slice(UniversalData &nums, Eigen::VectorXi &ind, UniversalData &A)
{
A = nums.slice_by_sample(ind);
}

void slice_restore(Eigen::VectorXd &A, Eigen::VectorXi &ind, Eigen::VectorXd &nums, int axis) {
if (ind.size() == 0) {
void slice_restore(Eigen::VectorXd &A, Eigen::VectorXi &ind, Eigen::VectorXd &nums, int axis)
{
if (ind.size() == 0)
{
nums = Eigen::VectorXd::Zero(nums.size());
} else {
}
else
{
nums = Eigen::VectorXd::Zero(nums.size());
for (int i = 0; i < ind.size(); i++) {
for (int i = 0; i < ind.size(); i++)
{
nums(ind(i)) = A(i);
}
}
return;
}

void coef_set_zero(int p, int M, Eigen::VectorXd& beta, Eigen::VectorXd& coef0) {
void coef_set_zero(int p, int M, Eigen::VectorXd &beta, Eigen::VectorXd &coef0)
{
beta = Eigen::VectorXd::Zero(p);
coef0 = Eigen::VectorXd::Zero(M);
return;
}



// Eigen::SparseMatrix<double> array_product(Eigen::SparseMatrix<double> &A, Eigen::VectorXd &B, int axis)
// {
// for (int i = 0; i < A.cols(); i++)
Expand All @@ -232,9 +272,6 @@ void coef_set_zero(int p, int M, Eigen::VectorXd& beta, Eigen::VectorXd& coef0)
// return A;
// }




// void matrix_sqrt(Eigen::MatrixXd &A, Eigen::MatrixXd &B)
// {
// A.sqrt().evalTo(B);
Expand Down Expand Up @@ -264,7 +301,6 @@ void coef_set_zero(int p, int M, Eigen::VectorXd& beta, Eigen::VectorXd& coef0)
// }
// }


// void set_nonzeros(Eigen::MatrixXd &X, Eigen::MatrixXd &x)
// {
// return;
Expand Down Expand Up @@ -323,20 +359,20 @@ void coef_set_zero(int p, int M, Eigen::VectorXd& beta, Eigen::VectorXd& coef0)
// return ((l2 == 0 || l1 / l2 > 1e+10) ? true : false);
// }


void init_spdlog(int console_log_level, int file_log_level, const char* log_file_name)
void init_spdlog(int console_log_level, int file_log_level, const char *log_file_name)
{
std::vector<spdlog::sink_ptr> sinks;

if(console_log_level != SPDLOG_LEVEL_OFF){
if (console_log_level != SPDLOG_LEVEL_OFF)
{
auto console_sink = std::make_shared<spdlog::sinks::stdout_color_sink_mt>();
console_sink->set_level(spdlog::level::level_enum(console_log_level));
console_sink->set_pattern("[%T.%e][%s:%#, %!][%^%l%$]: %v");
sinks.push_back(console_sink);
}


if(file_log_level != SPDLOG_LEVEL_OFF){
if (file_log_level != SPDLOG_LEVEL_OFF)
{
auto rotating_sink = std::make_shared<spdlog::sinks::rotating_file_sink_mt>(log_file_name, 1024 * 1024 * 10, 10);
rotating_sink->set_level(spdlog::level::level_enum(file_log_level));
rotating_sink->set_pattern("[%Y/%m/%d][%T.%e][elapsed %o][Process %P Thread %t][%s:%#, %!][%^%l%$]: %v");
Expand Down
Loading

0 comments on commit 450f023

Please sign in to comment.