diff --git a/src/classification/SVMWrapper.cpp b/src/classification/SVMWrapper.cpp
index 369fcf99338958c79cb4b97bee97bdc9f0cea7a5..097caf2a7d949abef144787c6aa8a694349de8d6 100644
--- a/src/classification/SVMWrapper.cpp
+++ b/src/classification/SVMWrapper.cpp
@@ -178,7 +178,7 @@ void SVMWrapper::copy_data(std::vector<double*> val_ptrs)
}
}
-void SVMWrapper::train()
+void SVMWrapper::train(bool remap_coefs)
{
_model = svm_train(&_prob, &_param);
std::copy_n(_model->rho, _intercept.size(), _intercept.data());
@@ -201,8 +201,11 @@ void SVMWrapper::train()
for(int sv = start_c2; sv < start_c2 + n_sv_c2; ++sv)
_coefs[class_combo][dd] -= _model->sv_coef[c1][sv] * _model->SV[sv][dd].value;
- _coefs[class_combo][dd] *= _w_remap[dd];
- _intercept[class_combo] -= _coefs[class_combo][dd] * _b_remap[dd];
+ if(remap_coefs)
+ {
+ _coefs[class_combo][dd] *= _w_remap[dd];
+ _intercept[class_combo] -= _coefs[class_combo][dd] * _b_remap[dd];
+ }
}
++class_combo;
}
@@ -214,16 +217,16 @@ void SVMWrapper::train()
_n_misclassified += _y[ss] != _y_est[ss];
}
-void SVMWrapper::train(std::vector<int> inds, int task)
+void SVMWrapper::train(std::vector<int> inds, int task, bool remap_coefs)
{
copy_data(inds, task);
- train();
+ train(remap_coefs);
}
-void SVMWrapper::train(std::vector<double*> val_ptrs)
+void SVMWrapper::train(std::vector<double*> val_ptrs, bool remap_coefs)
{
copy_data(val_ptrs);
- train();
+ train(remap_coefs);
}
std::vector<double> SVMWrapper::predict(int n_test_samp, std::vector<double*> val_ptrs)
diff --git a/src/classification/SVMWrapper.hpp b/src/classification/SVMWrapper.hpp
index 3c2e34770ecc7dad072dd4cda5ecbcf24fb6099f..96ec8bd1f108d503733f48d87088f27ef6e2f288 100644
--- a/src/classification/SVMWrapper.hpp
+++ b/src/classification/SVMWrapper.hpp
@@ -132,23 +132,27 @@ public:
/**
* @brief Train the SVM model
+ *
+ * @param remap_coefs If true remap the final coefficients back to the unscaled feature space
*/
- void train();
+ void train(bool remap_coefs=true);
/**
* @brief Copy the data from a set of feature indexes (sorted_dmatrix) into the x_space and train the SVM model
*
* @param inds list of indexes to pull data for
* @param task the task number for the calculation
+ * @param remap_coefs If true remap the final coefficients back to the unscaled feature space
*/
- void train(std::vector<int> inds, int task);
+ void train(std::vector<int> inds, int task, bool remap_coefs=true);
/**
* @brief Copy tthe data from a set data pointers and train the SVM model
*
* @param val_ptrs The pointers to the feature's data
+ * @param remap_coefs If true remap the final coefficients back to the unscaled feature space
*/
- void train(std::vector<double*> val_ptrs);
+ void train(std::vector<double*> val_ptrs, bool remap_coefs=true);
/**
* @brief Predict the class of a set of data from the SVM model
diff --git a/src/descriptor_identifier/SISSO_DI/SISSOClassifier.cpp b/src/descriptor_identifier/SISSO_DI/SISSOClassifier.cpp
index 1b916648831a7efd5d89ef540da474f2e191236a..df62b7d88628160efe7031b0b92ee18747fcc61d 100644
--- a/src/descriptor_identifier/SISSO_DI/SISSOClassifier.cpp
+++ b/src/descriptor_identifier/SISSO_DI/SISSOClassifier.cpp
@@ -125,6 +125,29 @@ void SISSOClassifier::setup_lp()
prop_sorted_d_mat::initialize_sroted_d_matrix_arr(0, _n_task, _n_class, n_samp_per_class);
}
+std::array<double, 2> SISSOClassifier::svm_error(std::vector<int>& feat_inds)
+{
+ double error = 0.0;
+ double dist_error = 0.0;
+ for(int tt = 0; tt < _n_task; ++tt)
+ {
+ _svm[tt].train(feat_inds, tt, false);
+ error += _svm[tt].n_misclassified();
+ for(auto& coefs : _svm[tt].coefs())
+ dist_error += 1.0 / util_funcs::norm(coefs.data(), feat_inds.size());
+ }
+ dist_error /= static_cast<double>(_n_task * _svm[0].n_class());
+ return {error, dist_error};
+}
+
+int SISSOClassifier::get_max_error_ind(int n_models, int* n_convex_overlap, double* svm_score, double* svm_margin, double* scores)
+{
+ std::transform(n_convex_overlap, n_convex_overlap, svm_score, scores, [this](int n_overlap, double score){return static_cast<double>(n_overlap * _n_samp) + score;});
+ double max_dist = *std::max_element(svm_margin, svm_margin) + 0.01;
+ std::transform(svm_margin, svm_margin, scores, scores, [&max_dist](double margin, double score){return score + (1.0 - margin / max_dist);});
+ return std::max_element(scores, scores) - scores;
+}
+
std::vector<std::vector<int>> SISSOClassifier::find_min_svm(std::vector<std::vector<int>>& test_inds, int n_keep)
{
int n_dim = test_inds[0].size();
@@ -183,22 +206,22 @@ void SISSOClassifier::transfer_d_mat_to_sorted()
}
}
-void SISSOClassifier::set_lp_data(int task, int cls)
+void SISSOClassifier::set_lp_data(int task, int cls, std::vector<int>& feat_inds)
{
- int n_dim = _cur_inds.size();
+ int n_dim = feat_inds.size();
for(int ii = 0; ii < n_dim; ++ii)
- dcopy_(_lp_n_col, prop_sorted_d_mat::access_sorted_d_matrix(_cur_inds[ii], task, cls), 1, &_lp_elements[ii], _lp_n_row);
+ dcopy_(_lp_n_col, prop_sorted_d_mat::access_sorted_d_matrix(feat_inds[ii], task, cls), 1, &_lp_elements[ii], _lp_n_row);
int n_copied = 0;
for(auto& cls_task_ind : _class_task_check[task * _n_class + cls])
{
- for(int ii = 0; ii < _cur_inds.size(); ++ii)
+ for(int ii = 0; ii < feat_inds.size(); ++ii)
{
- dcopy_(prop_sorted_d_mat::N_SAMPLES_PER_CLASS[cls_task_ind], prop_sorted_d_mat::access_sorted_d_matrix(_cur_inds[ii], task, cls_task_ind % _n_class), 1, &_lp_row_lower[n_copied + ii], _lp_n_row);
- dcopy_(prop_sorted_d_mat::N_SAMPLES_PER_CLASS[cls_task_ind], prop_sorted_d_mat::access_sorted_d_matrix(_cur_inds[ii], task, cls_task_ind % _n_class), 1, &_lp_row_upper[n_copied + ii], _lp_n_row);
+ dcopy_(prop_sorted_d_mat::N_SAMPLES_PER_CLASS[cls_task_ind], prop_sorted_d_mat::access_sorted_d_matrix(feat_inds[ii], task, cls_task_ind % _n_class), 1, &_lp_row_lower[n_copied + ii], _lp_n_row);
+ dcopy_(prop_sorted_d_mat::N_SAMPLES_PER_CLASS[cls_task_ind], prop_sorted_d_mat::access_sorted_d_matrix(feat_inds[ii], task, cls_task_ind % _n_class), 1, &_lp_row_upper[n_copied + ii], _lp_n_row);
}
- n_copied += prop_sorted_d_mat::N_SAMPLES_PER_CLASS[cls_task_ind] * _cur_inds.size();
+ n_copied += prop_sorted_d_mat::N_SAMPLES_PER_CLASS[cls_task_ind] * feat_inds.size();
}
- _sz_cls_task_check = n_copied / _cur_inds.size();
+ _sz_cls_task_check = n_copied / feat_inds.size();
_lp_simplex.loadProblem(
_lp_n_col,
@@ -229,7 +252,7 @@ int SISSOClassifier::pt_in_convex_hull_lp(int cls_task_ind)
return n_overlap;
}
-int SISSOClassifier::get_n_overlap_lp()
+int SISSOClassifier::get_n_overlap_lp(std::vector<int>& feat_inds)
{
int n_overlap = 0;
for(int tt = 0; tt < _n_task; ++tt)
@@ -237,7 +260,7 @@ int SISSOClassifier::get_n_overlap_lp()
for(int cc = 0; cc < _n_class; ++cc)
{
_lp_n_col = prop_sorted_d_mat::get_class_size(tt, cc);
- set_lp_data(tt, cc);
+ set_lp_data(tt, cc, feat_inds);
for(auto& cls_task_ind : _class_task_check[tt * _n_class + cc])
n_overlap += pt_in_convex_hull_lp(cls_task_ind);
}
@@ -263,10 +286,11 @@ void SISSOClassifier::l0_norm(std::vector<double>& prop, int n_dim)
std::vector<double> coefs(n_dim + 1, 0.0);
_cur_inds = std::vector<int>(n_dim, 0);
- std::vector<std::vector<std::vector<int>>> min_inds(n_get_models);
-
+ std::vector<int> min_inds(n_get_models * n_dim, -1);
std::vector<int> min_n_convex_overlap(n_get_models, ((n_dim > 1) ? models().back().back().n_convex_overlap_train() : _n_samp));
+ std::vector<double> min_svm_score(n_get_models, ((n_dim > 1) ? models().back().back().n_convex_overlap_train() : _n_samp));
+ std::vector<double> min_svm_margin(n_get_models, -1.0);
for(int rr = 0; rr < n_dim; ++rr)
_cur_inds[rr] = _feat_space->phi_selected().size() - 1 - rr;
@@ -274,69 +298,80 @@ void SISSOClassifier::l0_norm(std::vector<double>& prop, int n_dim)
util_funcs::iterate(_cur_inds, _cur_inds.size(), _mpi_comm->rank());
int max_error_ind = 0;
- std::vector<double> gamma(_prop.size(), 0.0);
- std::vector<double> gamma_feat(_prop.size(), 0.0);
+ #pragma omp parallel firstprivate(_svm, max_error_ind)
+ {
+ std::vector<int> min_inds_private(min_inds);
+ std::vector<int> min_n_convex_overlap_private(min_n_convex_overlap);
+ std::vector<double> min_svm_score_private(min_svm_score);
+ std::vector<double> min_svm_margin_private(min_svm_margin);
+ std::array<double, 2> temp_svm_error;
- do {
- int n_convex_overlap = get_n_overlap_lp();
- int el_num = std::find(min_n_convex_overlap.begin(), min_n_convex_overlap.end(), n_convex_overlap) - min_n_convex_overlap.begin();
+ std::vector<double> scores(n_get_models);
- if(el_num < n_get_models)
- {
- min_inds[el_num].push_back(_cur_inds);
- }
- else if(n_convex_overlap < min_n_convex_overlap[max_error_ind])
+ #pragma omp single
{
- min_n_convex_overlap[max_error_ind] = n_convex_overlap;
- min_inds[max_error_ind] = {_cur_inds};
- max_error_ind = std::max_element(min_n_convex_overlap.begin(), min_n_convex_overlap.end()) - min_n_convex_overlap.begin();
+ do {
+ #pragma omp task firstprivate(_cur_inds, _svm)
+ {
+ int n_convex_overlap = get_n_overlap_lp(_cur_inds);
+ if(n_convex_overlap <= min_n_convex_overlap_private[max_error_ind])
+ {
+ temp_svm_error = svm_error(_cur_inds);
+
+ if((n_convex_overlap < min_n_convex_overlap_private[max_error_ind]) || (temp_svm_error[0] < min_svm_score_private[max_error_ind]) || ((temp_svm_error[0] == min_svm_score_private[max_error_ind]) && (temp_svm_error[1] > min_svm_margin_private[max_error_ind])))
+ {
+ min_n_convex_overlap_private[max_error_ind] = n_convex_overlap;
+ min_svm_score_private[max_error_ind] = temp_svm_error[0];
+ min_svm_margin_private[max_error_ind] = temp_svm_error[1];
+ std::copy_n(_cur_inds.begin(), n_dim, min_inds_private.begin() + max_error_ind * n_dim);
+
+ max_error_ind = get_max_error_ind(n_get_models, min_n_convex_overlap_private.data(), min_svm_score_private.data(), min_svm_margin_private.data(), scores.data());
+ }
+ }
+ }
+ } while(util_funcs::iterate(_cur_inds, _cur_inds.size(), _mpi_comm->size()));
}
- } while(util_funcs::iterate(_cur_inds, _cur_inds.size(), _mpi_comm->size()));
-
- std::vector<int> inds_min(n_get_models * n_dim, -1);
- std::vector<int> sorted_overlap_inds = util_funcs::argsort(min_n_convex_overlap);
- std::vector<int> min_n_convex_overlap_copy(min_n_convex_overlap);
-
- int ii = 0;
- int nn = 0;
- while((ii < n_get_models) && (nn < min_n_convex_overlap.size()))
- {
- if(((ii + min_inds[sorted_overlap_inds[nn]].size()) * n_dim <= inds_min.size()))
+ #pragma omp barrier
+ std::cout << min_n_convex_overlap_private[0] << std::endl;
+ std::cout << min_svm_score_private[0] << std::endl;
+ std::cout << min_svm_margin_private[0] << std::endl;
+ std::cout << "check inds" << std::endl;
+ #pragma omp critical
{
- for(auto& mi : min_inds[sorted_overlap_inds[nn]])
+ max_error_ind = get_max_error_ind(n_get_models, min_n_convex_overlap.data(), min_svm_score.data(), min_svm_margin.data(), scores.data());
+ for(int ee = 0; ee < min_n_convex_overlap.size(); ++ee)
{
- std::copy_n(mi.begin(), mi.size(), &inds_min[ii * n_dim]);
- min_n_convex_overlap[ii] = min_n_convex_overlap_copy[sorted_overlap_inds[nn]];
- ++ii;
+ std::cout << min_n_convex_overlap_private[ee] << '\t' << min_n_convex_overlap[max_error_ind] << std::endl;
+ if((min_n_convex_overlap_private[ee] < min_n_convex_overlap[max_error_ind]) || (min_svm_score_private[ee] < min_svm_score[max_error_ind]) || ((min_svm_score_private[ee] == min_svm_score[max_error_ind]) && (min_svm_margin_private[ee] > min_svm_margin[max_error_ind])))
+ {
+ min_n_convex_overlap[max_error_ind] = min_n_convex_overlap_private[ee];
+ min_svm_score[max_error_ind] = min_svm_score_private[ee];
+ min_svm_margin[max_error_ind] = min_svm_margin_private[ee];
+ std::copy_n(min_inds_private.begin() + ee * n_dim, n_dim, min_inds.begin() + max_error_ind * n_dim);
+
+ max_error_ind = get_max_error_ind(n_get_models, min_n_convex_overlap.data(), min_svm_score.data(), min_svm_margin.data(), scores.data());
+ }
}
}
- else
- {
- std::vector<std::vector<int>> test_inds = find_min_svm(min_inds[sorted_overlap_inds[nn]], n_get_models - ii);
- int kk = 0;
- while((ii < inds_min.size()) && (kk < test_inds.size()))
- {
- std::copy_n(test_inds[kk].begin(), test_inds[kk].size(), &inds_min[ii * n_dim]);
- min_n_convex_overlap[ii] = min_n_convex_overlap_copy[sorted_overlap_inds[nn]];
- ++ii;
- ++kk;
- }
- }
- ++nn;
}
-
+ std::cout << "mpi" << std::endl;
std::vector<int> all_min_n_convex_overlap(_mpi_comm->size() * n_get_models);
- std::vector<int> all_inds_min(_mpi_comm->size() * n_get_models * n_dim);
+ std::vector<double> all_min_svm_score(_mpi_comm->size() * n_get_models);
+ std::vector<double> all_min_svm_margin(_mpi_comm->size() * n_get_models);
+ std::vector<int> all_min_inds(_mpi_comm->size() * n_get_models * n_dim);
mpi::all_gather(*_mpi_comm, min_n_convex_overlap.data(), n_get_models, all_min_n_convex_overlap);
- mpi::all_gather(*_mpi_comm, inds_min.data(), n_get_models * n_dim, all_inds_min);
+ mpi::all_gather(*_mpi_comm, min_svm_score.data(), n_get_models, all_min_svm_score);
+ mpi::all_gather(*_mpi_comm, min_svm_margin.data(), n_get_models, all_min_svm_margin);
+ mpi::all_gather(*_mpi_comm, min_inds.data(), n_get_models * n_dim, all_min_inds);
std::vector<int> inds = util_funcs::argsort(all_min_n_convex_overlap);
- std::vector<int> all_inds_min_original(all_inds_min);
+ std::vector<int> all_min_inds_original(all_min_inds);
- ii = 1;
- nn = 0;
+ int ii = 1;
+ int nn = 0;
+ std::cout << "find" << std::endl;
while(nn < n_get_models)
{
while((ii < all_min_n_convex_overlap.size()) && (all_min_n_convex_overlap[inds[ii]] == all_min_n_convex_overlap[inds[nn]]))
@@ -349,26 +384,32 @@ void SISSOClassifier::l0_norm(std::vector<double>& prop, int n_dim)
}
std::vector<std::vector<int>> test_inds(ii - nn, std::vector<int>(n_dim, 0));
for(int cc = nn; cc < ii; ++cc)
- std::copy_n(&all_inds_min_original[inds[cc] * n_dim], n_dim, test_inds[cc - nn].begin());
+ std::copy_n(&all_min_inds_original[inds[cc] * n_dim], n_dim, test_inds[cc - nn].begin());
test_inds = find_min_svm(test_inds, std::min(ii - nn, n_get_models - nn));
for(int cc = 0; cc < test_inds.size(); ++cc)
- std::copy_n(test_inds[cc].begin(), n_dim, &all_inds_min[(cc + nn) * n_dim]);
+ std::copy_n(test_inds[cc].begin(), n_dim, &all_min_inds[(cc + nn) * n_dim]);
nn = ii;
++ii;
}
std::vector<model_node_ptr> min_nodes(n_dim);
std::vector<ModelClassifier> models;
+
+ std::cout << "models" << std::endl;
+ #pragma omp parallel for
for(int rr = 0; rr < n_get_models; ++rr)
{
node_value_arrs::clear_temp_test_reg();
for(int ii = 0; ii < n_dim; ++ii)
{
- int index = all_inds_min[rr * n_dim + ii];
+ int index = all_min_inds[rr * n_dim + ii];
+ std::cout << index << std::endl;
min_nodes[ii] = std::make_shared<ModelNode>(_feat_space->phi_selected()[index]->arr_ind(), _feat_space->phi_selected()[index]->rung(), _feat_space->phi_selected()[index]->expr(), _feat_space->phi_selected()[index]->postfix_expr(), _feat_space->phi_selected()[index]->value(), _feat_space->phi_selected()[index]->test_value(), _feat_space->phi_selected()[index]->unit());
}
+
+ #pragma omp critical
models.push_back(ModelClassifier(_prop_unit, _prop, _prop_test, min_nodes, _task_sizes_train, _task_sizes_test, _fix_intercept));
}
diff --git a/src/descriptor_identifier/SISSO_DI/SISSOClassifier.hpp b/src/descriptor_identifier/SISSO_DI/SISSOClassifier.hpp
index 9ce400a5e4e4dc9d6cdf4450edb99d91134d3d71..cccc8fa9ddd587e55d2cc9fb79d9e7f307dd0f48 100644
--- a/src/descriptor_identifier/SISSO_DI/SISSOClassifier.hpp
+++ b/src/descriptor_identifier/SISSO_DI/SISSOClassifier.hpp
@@ -119,8 +119,9 @@ public:
*
* @param task The task number to set up the data objects for
* @param cls The class number to set up the data objects for
+ * @param feat_inds The indexes of the features to create the model for
*/
- void set_lp_data(int task, int cls);
+ void set_lp_data(int task, int cls, std::vector<int>& feat_inds);
/**
* @brief Check how many points in a given set are inside the convex hull
@@ -134,9 +135,31 @@ public:
/**
* @brief Get the number of points in overlapping convex hull regions
* @details Loop over all task/class combinations and finds the total number of points in overlapping
+ *
+ * @param feat_inds The indexes to develop the SVM model for
* @return The total number of overlapping points in all convex hull regions.
*/
- int get_n_overlap_lp();
+ int get_n_overlap_lp(std::vector<int>& feat_inds);
+
+ /**
+ * @brief Find the SVM error and margin for a given set of indexes
+ *
+ * @param feat_inds The indexes to develop the SVM model for
+ * @return {The number of misclassified points, the SVM margin}
+ */
+ std::array<double, 2> svm_error(std::vector<int>& feat_inds);
+
+ /**
+ * @brief Gets the max error index for the classification problem
+ * @details [long description]
+ *
+ * @param n_models number of models to be stored
+ * @param n_convex_overlap number of points in the overlapping convex hull regions (in all models)
+ * @param svm_score the number of points misclassified by SVM (in all models)
+ * @param svm_margin The margin of the SVM model (in all models)
+ * @return [description]
+ */
+ int get_max_error_ind(int n_models, int* n_convex_overlap, double* svm_score, double* svm_margin, double* scores);
/**
* @brief Preform the l0 normalization for a property or the residual
diff --git a/src/descriptor_identifier/SISSO_DI/SISSORegressor.cpp b/src/descriptor_identifier/SISSO_DI/SISSORegressor.cpp
index 92d1c9500f129a1872898bcd2b6f21406e909011..d0eb107f495fe31b25b5fb394455ed5b4997aa02 100644
--- a/src/descriptor_identifier/SISSO_DI/SISSORegressor.cpp
+++ b/src/descriptor_identifier/SISSO_DI/SISSORegressor.cpp
@@ -27,27 +27,28 @@ SISSORegressor::SISSORegressor(
_work = std::vector<double>(_lwork, 0.0);
}
-void SISSORegressor::set_a(std::vector<int>& inds, int start, int n_samp)
+void SISSORegressor::set_a(std::vector<int>& inds, int start, int n_samp, double* a)
{
for(int ii = 0; ii < inds.size(); ++ii)
- std::copy_n(node_value_arrs::get_d_matrix_ptr(inds[ii]) + start, n_samp, _a.data() + ii * n_samp);
+ std::copy_n(node_value_arrs::get_d_matrix_ptr(inds[ii]) + start, n_samp, a + ii * n_samp);
+
if(!_fix_intercept)
- std::fill_n(_a.data() + inds.size() * n_samp, n_samp, 1.0);
+ std::fill_n(a + inds.size() * n_samp, n_samp, 1.0);
}
-void SISSORegressor::least_squares(std::vector<int>& inds, double* coeffs, int start, int n_samp)
+void SISSORegressor::least_squares(std::vector<int>& inds, double* coeffs, int start, int n_samp, double* a, double*b, double* work)
{
int info;
int n_dim = inds.size() + (1 - _fix_intercept);
- set_a(inds, start, n_samp);
- std::copy_n(_prop.data() + start, n_samp, _b.data());
+ set_a(inds, start, n_samp, a);
+ std::copy_n(_prop.data() + start, n_samp, b);
- dgels_('N', n_samp, n_dim, 1, _a.data(), n_samp, _b.data(), n_samp, _work.data(), _lwork, &info);
+ dgels_('N', n_samp, n_dim, 1, a, n_samp, b, n_samp, work, _lwork, &info);
if(info == 0)
{
- std::copy_n(_b.data(), n_dim, coeffs);
+ std::copy_n(b, n_dim, coeffs);
}
else
{
@@ -71,13 +72,13 @@ int SISSORegressor::get_opt_lwork(int n_dim)
throw std::logic_error("Failed to get lwork.");
}
-void SISSORegressor::set_error(std::vector<int>& inds, double* coeffs, int start, int n_samp)
+void SISSORegressor::set_error(std::vector<int>& inds, double* coefs, int start, int n_samp, double* error)
{
- std::fill_n(_error.data() + start, n_samp, -1 * (!_fix_intercept) * coeffs[inds.size()]);
- daxpy_(n_samp, 1.0, _prop.data() + start, 1, _error.data() + start, 1);
-
+ std::fill_n(error + start, n_samp, -1 * (!_fix_intercept) * coefs[inds.size()]);
for(int ii = 0; ii < inds.size(); ++ii)
- daxpy_(n_samp, -1.0 * coeffs[ii], node_value_arrs::get_d_matrix_ptr(inds[ii]) + start, 1, _error.data() + start, 1);
+ daxpy_(n_samp, -1.0 * coefs[ii], node_value_arrs::get_d_matrix_ptr(inds[ii]) + start, 1, error + start, 1);
+ daxpy_(n_samp, 1.0, _prop.data() + start, 1, error + start, 1);
+
}
void SISSORegressor::l0_norm(std::vector<double>& prop, int n_dim)
@@ -90,36 +91,67 @@ void SISSORegressor::l0_norm(std::vector<double>& prop, int n_dim)
std::vector<double> min_errors(n_get_models, util_funcs::norm(_prop.data(), _n_samp));
+ int n_interactions = 1;
+ int n_dim_fact = 1;
for(int rr = 0; rr < n_dim; ++rr)
+ {
inds[rr] = _feat_space->phi_selected().size() - 1 - rr;
+ n_interactions *= inds[rr];
+ n_dim_fact *= (rr + 1);
+ }
+ n_interactions /= n_dim_fact;
util_funcs::iterate(inds, inds.size(), _mpi_comm->rank());
- int max_error_ind = 0;
if(inds.back() >= 0)
{
- do {
- int start = 0;
- double error = 0.0;
- for(auto& sz : _task_sizes_train)
+ #pragma omp parallel shared(inds_min, min_errors, n_interactions, n_get_models, n_dim) firstprivate(inds, coefs, _error, _a, _b, _work)
+ {
+ int max_error_ind = 0;
+ std::vector<int> inds_min_private(inds_min);
+ std::vector<double> min_errors_private(min_errors);
+
+ int ii_prev = 0;
+ #pragma omp for schedule(dynamic)
+ for(int ii = 0; ii < n_interactions; ii += _mpi_comm->size())
{
- least_squares(inds, coefs.data(), start, sz);
- set_error(inds, coefs.data(), start, sz);
- error += ddot_(sz, _error.data() + start, 1, _error.data() + start, 1) / sz;
- // error += std::pow(util_funcs::norm(_error.data() + start, sz), 2.0) / sz;
- start += sz;
+ util_funcs::iterate(inds, inds.size(), _mpi_comm->size() * (ii - ii_prev));
+
+ int start = 0;
+ for(auto& sz : _task_sizes_train)
+ {
+ least_squares(inds, coefs.data(), start, sz, _a.data(), _b.data(), _work.data());
+ set_error(inds, coefs.data(), start, sz, _error.data());
+ start += sz;
+ }
+ double error = std::sqrt(ddot_(_n_samp, _error.data(), 1, _error.data(), 1) / _n_samp);
+ #pragma omp critical
+ {
+ if(error < min_errors_private[max_error_ind])
+ {
+ min_errors_private[max_error_ind] = error;
+ std::copy_n(inds.begin(), inds.size(), inds_min_private.begin() + max_error_ind * n_dim);
+ max_error_ind = std::max_element(min_errors_private.begin(), min_errors_private.end()) - min_errors_private.begin();
+ }
+ }
+ ii_prev = ii;
}
- error = std::sqrt(error / _n_task);
- if(error < min_errors[max_error_ind])
- {
- min_errors[max_error_ind] = error;
- std::copy_n(inds.begin(), inds.size(), inds_min.begin() + max_error_ind * n_dim);
+ #pragma omp critical
+ {
max_error_ind = std::max_element(min_errors.begin(), min_errors.end()) - min_errors.begin();
+ for(int ee = 0; ee < n_get_models; ++ee)
+ {
+ if(min_errors_private[ee] < min_errors[max_error_ind])
+ {
+ min_errors[max_error_ind] = min_errors_private[ee];
+ std::copy_n(inds_min_private.begin() + ee * n_dim, n_dim, inds_min.begin() + max_error_ind * n_dim);
+ max_error_ind = std::max_element(min_errors.begin(), min_errors.end()) - min_errors.begin();
+ }
+ }
}
- } while(util_funcs::iterate(inds, inds.size(), _mpi_comm->size()));
+ }
}
-
std::vector<double> all_min_error(_mpi_comm->size() * n_get_models);
std::vector<int> all_inds_min(_mpi_comm->size() * n_get_models * n_dim);
@@ -131,6 +163,7 @@ void SISSORegressor::l0_norm(std::vector<double>& prop, int n_dim)
std::vector<model_node_ptr> min_nodes(n_dim);
std::vector<ModelRegressor> models;
+ // #pragma omp parallel for firstprivate(min_nodes)
for(int rr = 0; rr < n_get_models; ++rr)
{
node_value_arrs::clear_temp_test_reg();
@@ -139,6 +172,8 @@ void SISSORegressor::l0_norm(std::vector<double>& prop, int n_dim)
int index = all_inds_min[inds[rr] * n_dim + ii];
min_nodes[ii] = std::make_shared<ModelNode>(_feat_space->phi_selected()[index]->arr_ind(), _feat_space->phi_selected()[index]->rung(), _feat_space->phi_selected()[index]->expr(), _feat_space->phi_selected()[index]->postfix_expr(), _feat_space->phi_selected()[index]->value(), _feat_space->phi_selected()[index]->test_value(), _feat_space->phi_selected()[index]->unit());
}
+
+ // #pragma omp critical
models.push_back(ModelRegressor(_prop_unit, _prop, _prop_test, min_nodes, _task_sizes_train, _task_sizes_test, _fix_intercept));
}
diff --git a/src/descriptor_identifier/SISSO_DI/SISSORegressor.hpp b/src/descriptor_identifier/SISSO_DI/SISSORegressor.hpp
index 939994dfb8e3cbb9163268b4986974086b57f506..83b5d191561af4fd79ce65f7a83a66d20e82a2b1 100644
--- a/src/descriptor_identifier/SISSO_DI/SISSORegressor.hpp
+++ b/src/descriptor_identifier/SISSO_DI/SISSORegressor.hpp
@@ -59,7 +59,7 @@ public:
* @param start The index in the property and feature vectors start copying into _b and _a
* @param n_samp number of samples to perform least squares optimization on
*/
- void set_error(std::vector<int>& inds, double* coeffs, int start, int n_samp);
+ void set_error(std::vector<int>& inds, double* coefs, int start, int n_samp, double* error);
/**
* @brief Get the optimal size of the working array
@@ -77,7 +77,7 @@ public:
* @param start The index in the property and feature vectors start copying into _b and _a
* @param n_samp number of samples to perform least squares optimization on
*/
- void least_squares(std::vector<int>& inds, double* coeffs, int start, int n_samp);
+ void least_squares(std::vector<int>& inds, double* coefs, int start, int n_samp, double* a, double*b, double* work);
/**
* @brief Set the A matrix for the least squares problem
@@ -86,7 +86,7 @@ public:
* @param start The index in the property and feature vectors start copying into _b and _a
* @param n_samp number of samples to perform least squares optimization on
*/
- void set_a(std::vector<int>& inds, int start, int n_samp);
+ void set_a(std::vector<int>& inds, int start, int n_samp, double* a);
/**
* @brief Preform the l0 normalization for a property or the residual
diff --git a/src/feature_creation/feature_space/FeatureSpace.cpp b/src/feature_creation/feature_space/FeatureSpace.cpp
index 7b33ce2c5bc8fc198690ffb0876ca6405f6a3678..f0ab6be199d924df17693cab1165313a70a6db39 100644
--- a/src/feature_creation/feature_space/FeatureSpace.cpp
+++ b/src/feature_creation/feature_space/FeatureSpace.cpp
@@ -471,51 +471,53 @@ void FeatureSpace::project_generated(double* prop, int size, std::vector<node_pt
std::vector<node_ptr> generated_phi;
generate_new_feats(feat, generated_phi, feat_ind, _l_bound, _u_bound);
- std::vector<double> scores(generated_phi.size(), 0.0);
- _project(prop, scores.data(), generated_phi, _task_sizes, size / _n_samp);
-
- std::vector<int> inds = util_funcs::argsort(scores);
+ #pragma omp critical
+ {
+ std::vector<double> scores(generated_phi.size(), 0.0);
+ _project(prop, scores.data(), generated_phi, _task_sizes, size / _n_samp);
- int ii = 0;
+ std::vector<int> inds = util_funcs::argsort(scores);
- while(_scores[inds[ii]] < -1.0)
- ++ii;
+ int ii = 0;
- while((ii < inds.size()) && (scores[inds[ii]] < worst_score))
- {
- double cur_score = scores[inds[ii]];
+ while(scores[inds[ii]] < -1.0)
+ ++ii;
- if(valid_feature_against_selected(generated_phi[inds[ii]]->value_ptr(), node_value_arrs::N_SELECTED - _n_sis_select + scores_sel_private.size()))
+ while((ii < inds.size()) && (scores[inds[ii]] < worst_score))
{
- if(scores_sel_private.size() == _n_sis_select)
+ double cur_score = scores[inds[ii]];
+
+ if(valid_feature_against_selected(generated_phi[inds[ii]]->value_ptr(), node_value_arrs::N_SELECTED - _n_sis_select + scores_sel_private.size()))
{
- phi_sel_private[worst_score_ind]->set_selected(false);
- phi_sel_private[worst_score_ind]->set_d_mat_ind(-1);
+ if(scores_sel_private.size() == _n_sis_select)
+ {
+ phi_sel_private[worst_score_ind]->set_selected(false);
+ phi_sel_private[worst_score_ind]->set_d_mat_ind(-1);
- phi_sel_private[worst_score_ind] = generated_phi[inds[ii]];
- scores_sel_private[worst_score_ind] = cur_score;
+ phi_sel_private[worst_score_ind] = generated_phi[inds[ii]];
+ scores_sel_private[worst_score_ind] = cur_score;
- phi_sel_private[worst_score_ind]->set_selected(true);
- phi_sel_private[worst_score_ind]->set_d_mat_ind(node_value_arrs::N_SELECTED - _n_sis_select + worst_score_ind);
- phi_sel_private[worst_score_ind]->set_value();
- }
- else
- {
+ phi_sel_private[worst_score_ind]->set_selected(true);
+ phi_sel_private[worst_score_ind]->set_d_mat_ind(node_value_arrs::N_SELECTED - _n_sis_select + worst_score_ind);
+ phi_sel_private[worst_score_ind]->set_value();
+ }
+ else
+ {
- phi_sel_private.push_back(generated_phi[inds[ii]]);
- scores_sel_private.push_back(cur_score);
+ phi_sel_private.push_back(generated_phi[inds[ii]]);
+ scores_sel_private.push_back(cur_score);
- phi_sel_private.back()->set_selected(true);
- phi_sel_private.back()->set_d_mat_ind(node_value_arrs::N_SELECTED - _n_sis_select + scores_sel_private.size());
- phi_sel_private.back()->set_value();
+ phi_sel_private.back()->set_selected(true);
+ phi_sel_private.back()->set_d_mat_ind(node_value_arrs::N_SELECTED - _n_sis_select + scores_sel_private.size());
+ phi_sel_private.back()->set_value();
+ }
+ worst_score_ind = std::max_element(scores_sel_private.begin(), scores_sel_private.end()) - scores_sel_private.begin();
+ worst_score = scores_sel_private[worst_score_ind];
}
- worst_score_ind = std::max_element(scores_sel_private.begin(), scores_sel_private.end()) - scores_sel_private.begin();
- worst_score = scores_sel_private[worst_score_ind];
+ ++ii;
}
- ++ii;
}
}
-
#pragma omp critical
{
worst_score_ind = std::max_element(scores_sel.begin(), scores_sel.end()) - scores_sel.begin();
diff --git a/src/utils/project.cpp b/src/utils/project.cpp
index 52e96203fe24fa6e352b8e165549ced849f9e1ca..1fd04c6071c0ef9697e3e23e7f24a36e74e83214 100644
--- a/src/utils/project.cpp
+++ b/src/utils/project.cpp
@@ -18,7 +18,7 @@ void project_funcs::project_r2(double* prop, double* scores, std::vector<node_pt
{
int n_samp = std::accumulate(sizes.begin(), sizes.end(), 0);
- #pragma omp parallel
+ #pragma omp parallel firstprivate(prop)
{
int start = omp_get_thread_num() * (phi.size() / omp_get_num_threads()) + std::min(omp_get_thread_num(), static_cast<int>(phi.size()) % omp_get_num_threads());
int end = (omp_get_thread_num() + 1) * (phi.size() / omp_get_num_threads()) + std::min(omp_get_thread_num() + 1, static_cast<int>(phi.size()) % omp_get_num_threads());
@@ -31,6 +31,8 @@ void project_funcs::project_r2(double* prop, double* scores, std::vector<node_pt
std::transform(phi.begin() + start, phi.begin() + end, scores + start, scores + start, [&prop, &sizes](node_ptr feat, double sc){return std::min(-1 * util_funcs::r2(prop, feat->value_ptr(), sizes.data(), sizes.size()), sc);});
}
std::transform(scores + start, scores + end, scores + start, [](double score){return std::isnan(score) ? 0.0 : score;});
+
+ #pragma omp barrier
}
}
@@ -40,7 +42,7 @@ void project_funcs::project_classify(double* prop, double* scores, std::vector<n
convex_hull::initialize_projection(sizes, prop);
- #pragma omp parallel
+ #pragma omp parallel firstprivate(prop)
{
int start = omp_get_thread_num() * (phi.size() / omp_get_num_threads()) + std::min(omp_get_thread_num(), static_cast<int>(phi.size()) % omp_get_num_threads());
int end = (omp_get_thread_num() + 1) * (phi.size() / omp_get_num_threads()) + std::min(omp_get_thread_num() + 1, static_cast<int>(phi.size()) % omp_get_num_threads());
@@ -54,5 +56,6 @@ void project_funcs::project_classify(double* prop, double* scores, std::vector<n
std::transform(phi.begin() + start, phi.begin() + end, scores + start, scores + start, [](node_ptr feat, double score){return std::min(convex_hull::overlap_1d(feat->value_ptr()), score);});
}
std::transform(scores + start, scores + end, scores + start, [](double score){return std::isnan(score) ? std::numeric_limits<double>::infinity() : score;});
+ #pragma omp barrier
}
}