Commit 1902fece authored by Thomas Purcell's avatar Thomas Purcell
Browse files

Adding: Multiple residuals fitting and leave-out CV keys

Did them one after another, adding both here now
parent 0d445757
#include <descriptor_identifier/Model/Model.hpp>
Model::Model(std::vector<double> prop, std::vector<std::shared_ptr<Node>> feats) :
_n_samp(feats[0]->n_samp()),
Model::Model(std::vector<double> prop_train, std::vector<double> prop_test, std::vector<std::shared_ptr<Node>> feats) :
_n_samp_train(feats[0]->n_samp()),
_n_samp_test(feats[0]->n_test_samp()),
_n_dim(feats.size() + 1),
_feats(feats),
_coefs(new double[_n_dim]),
_prop(new double[_n_samp]),
_error(new double[_n_samp]),
_D(new double[_n_samp * _n_dim]),
_prop_est(_n_samp, 0.0)
_coefs(_n_dim),
_prop_train(prop_train),
_prop_test(prop_test),
_train_error(_n_samp_train),
_test_error(_n_samp_test),
_D_train(_n_samp_train * _n_dim),
_D_test(_n_samp_test * _n_dim),
_prop_train_est(_n_samp_train, 0.0),
_prop_test_est(_n_samp_test, 0.0)
{
_prop_est.reserve(_n_samp);
_prop_train_est.reserve(_n_samp_train);
_prop_test_est.reserve(_n_samp_test);
std::copy_n(prop.begin(), _n_samp, _prop.get());
std::vector<double> a(_n_samp * _n_dim, 1.0);
std::vector<double> a(_n_samp_train * _n_dim, 1.0);
for(int ff = 0; ff < feats.size(); ++ff)
{
std::copy_n(feats[ff]->value_ptr(), _n_samp, _D.get() + ff * _n_samp);
std::copy_n(feats[ff]->value_ptr(), _n_samp, a.begin() + ff * _n_samp);
std::copy_n(feats[ff]->value_ptr(), _n_samp_train, _D_train.data() + ff * _n_samp_train);
std::copy_n(feats[ff]->test_value_ptr(), _n_samp_test, _D_test.data() + ff * _n_samp_test);
std::copy_n(feats[ff]->value_ptr(), _n_samp_train, a.data() + ff * _n_samp_train);
}
std::copy_n(a.begin() + feats.size() * _n_samp, _n_samp, _D.get() + feats.size() * _n_samp);
std::copy_n(a.data() + feats.size() * _n_samp_train, _n_samp_train, _D_train.data() + feats.size() * _n_samp_train);
std::copy_n(a.data() + feats.size() * _n_samp_train, _n_samp_test, _D_test.data() + feats.size() * _n_samp_test);
std::vector<double> s(_n_dim, 0.0);
std::vector<double> work(_n_dim * _n_samp, 0.0);
std::vector<double> work(_n_dim * _n_samp_train, 0.0);
int rank = 0;
int info = 0;
dgelss_(_n_samp, _n_dim, 1, a.data(), _n_samp, prop.data(), _n_samp, s.data(), 1e-13, &rank, work.data(), work.size(), &info);
std::copy_n(prop.begin(), _n_dim, _coefs.get());
dgemv_('N', _n_samp, _n_dim, 1.0, _D.get(), _n_samp, _coefs.get(), 1, 0.0, _prop_est.data(), 1);
dgelss_(_n_samp_train, _n_dim, 1, a.data(), _n_samp_train, prop_train.data(), _n_samp_train, s.data(), 1e-13, &rank, work.data(), work.size(), &info);
std::copy_n(prop_train.begin(), _n_dim, _coefs.data());
std::transform(_prop_est.begin(), _prop_est.end(), _prop.get(), _error.get(), std::minus<double>());
}
dgemv_('N', _n_samp_train, _n_dim, 1.0, _D_train.data(), _n_samp_train, _coefs.data(), 1, 0.0, _prop_train_est.data(), 1);
dgemv_('N', _n_samp_test, _n_dim, 1.0, _D_test.data(), _n_samp_test, _coefs.data(), 1, 0.0, _prop_test_est.data(), 1);
Model::Model(Model &o) :
_n_samp(o._n_samp),
_n_dim(o._n_dim),
_feats(o._feats),
_coefs(new double[_n_dim]),
_prop(new double[_n_samp]),
_error(new double[_n_samp]),
_D(new double[_n_samp * _n_dim]),
_prop_est(o._prop_est)
{
std::copy_n(o._coefs.get(), o._n_dim, _coefs.get());
std::copy_n(o._prop.get(), o._n_samp, _prop.get());
std::copy_n(o._error.get(), o._n_samp, _error.get());
std::copy_n(o._D.get(), o._n_samp * o._n_dim, _D.get());
}
Model::Model(Model &&o) :
_n_samp(o._n_samp),
_n_dim(o._n_dim),
_feats(o._feats),
_coefs(new double[_n_dim]),
_prop(new double[_n_samp]),
_error(new double[_n_samp]),
_D(new double[_n_samp * _n_dim]),
_prop_est(o._prop_est)
{
std::copy_n(o._coefs.get(), o._n_dim, _coefs.get());
std::copy_n(o._prop.get(), o._n_samp, _prop.get());
std::copy_n(o._error.get(), o._n_samp, _error.get());
std::copy_n(o._D.get(), o._n_samp * o._n_dim, _D.get());
std::transform(_prop_train_est.begin(), _prop_train_est.end(), _prop_train.data(), _train_error.data(), std::minus<double>());
std::transform(_prop_test_est.begin(), _prop_test_est.end(), _prop_test.data(), _test_error.data(), std::minus<double>());
}
std::string Model::toString() const
......@@ -81,3 +57,59 @@ std::ostream& operator<< (std::ostream& outStream, const Model& model)
outStream << model.toString();
return outStream;
}
void Model::train_to_file(std::string filename)
{
boost::filesystem::path p(filename.c_str());
boost::filesystem::create_directories(p.remove_filename());
std::ofstream out_file_stream = std::ofstream();
out_file_stream.open(filename);
out_file_stream << "# " << toString() << std::endl;
out_file_stream << "# RMSE: " << rmse() << "; Max AE( " << max_ae() << std::endl;
out_file_stream << "# coeffs:";
for(auto& coef: _coefs)
out_file_stream << " " << std::setw(24) << std::setprecision(18) << coef << ";";
out_file_stream << "\n# " << std::setw(23) << "Property Value," << std::setw(24) << "Property Value (EST),";
for(int ff = 0; ff < _feats.size(); ++ff)
out_file_stream << " Feature " << ff << " Value,";
out_file_stream << std::endl;
for(int ss = 0; ss < _n_samp_train; ++ss)
{
out_file_stream << std::setw(24) << std::setprecision(18) << _prop_train[ss] << std::setw(24) << std::setprecision(18) << _prop_train_est[ss];
for(int ff = 0; ff < _n_dim - 1; ++ff)
out_file_stream << std::setw(24) << std::setprecision(18) << _D_train[ss + ff * _n_samp_train];
out_file_stream << std::endl;
}
out_file_stream.close();
}
void Model::test_to_file(std::string filename)
{
boost::filesystem::path p(filename.c_str());
boost::filesystem::create_directories(p.remove_filename());
std::ofstream out_file_stream = std::ofstream();
out_file_stream.open(filename);
out_file_stream << "# " << toString() << std::endl;
out_file_stream << "# RMSE: " << test_rmse() << "; Max AE: " << test_max_ae() << std::endl;
out_file_stream << "# coeffs:";
for(auto& coef: _coefs)
out_file_stream << " " << std::setw(24) << std::setprecision(18) << coef << ";";
out_file_stream << "\n# " << std::setw(23) << "Property Value," << std::setw(24) << "Property Value (EST),";
for(int ff = 0; ff < _feats.size(); ++ff)
out_file_stream << " Feature " << ff << " Value,";
out_file_stream << std::endl;
for(int ss = 0; ss < _n_samp_test; ++ss)
{
out_file_stream << std::setw(24) << std::setprecision(18) << _prop_test[ss] << std::setw(24) << std::setprecision(18) << _prop_test_est[ss];
for(int ff = 0; ff < _n_dim - 1; ++ff)
out_file_stream << std::setw(24) << std::setprecision(18) << _D_test[ss + ff * _n_samp_train];
out_file_stream << std::endl;
}
out_file_stream.close();
}
#ifndef MODEL
#define MODEL
#include <boost/filesystem.hpp>
#include<iomanip>
#include<fstream>
#include<iostream>
#include <feature_creation/node/Node.hpp>
/**
......@@ -9,17 +15,22 @@
*/
class Model
{
int _n_samp; //!< The number of samples per feature
int _n_samp_train; //!< The number of samples per feature
int _n_samp_test; //!< The number of test samples per feature
int _n_dim; //!< Dimension of the model
std::vector<std::shared_ptr<Node>> _feats; //!< List of features in the model
std::unique_ptr<double[]> _coefs; //!< Coefficients for teh features
std::unique_ptr<double[]> _prop; //!< The property to be modeled
std::unique_ptr<double[]> _error; //!< The error of the model
std::unique_ptr<double[]> _D; //!< The Descriptor matrix
std::vector<double> _coefs; //!< Coefficients for teh features
std::vector<double> _prop_train; //!< The property to be modeled
std::vector<double> _prop_test; //!< The property to be modeled
std::vector<double> _train_error; //!< The error of the model
std::vector<double> _test_error; //!< The error of the model
std::vector<double> _D_train; //!< The Descriptor matrix
std::vector<double> _D_test; //!< The Descriptor matrix
std::vector<double> _prop_est; //!< The estimated Property
std::vector<double> _prop_train_est; //!< The estimated Property
std::vector<double> _prop_test_est; //!< The estimated Property
public:
/**
......@@ -28,21 +39,8 @@ public:
* @param prop The property
* @param feats The features for the model
*/
Model(std::vector<double> prop, std::vector<std::shared_ptr<Node>> feats);
/**
* @brief The copy constructor
*
* @param o The model to be copied
*/
Model(Model& o);
Model(std::vector<double> prop_train, std::vector<double> prop_test, std::vector<std::shared_ptr<Node>> feats);
/**
* @brief The Move constructor
*
* @param o The Model to be moved
*/
Model(Model&& o);
/**
* @brief Convert the model to a string
......@@ -54,27 +52,49 @@ public:
/**
* @brief Accessor function to _prop_est
*/
inline std::vector<double>& predict(){return _prop_est;}
inline std::vector<double>& predict(){return _prop_test_est;}
/**
* @brief Accessor function to _prop_est
*/
inline std::vector<double>& predict_train(){return _prop_train_est;}
/**
* @brief Copy the error into a new array
*
* @param res pointer to the beginning of the array
*/
inline void copy_error(double* res){std::copy_n(_error.get(), _n_samp, res);}
inline void copy_error(double* res){std::copy_n(_train_error.data(), _n_samp_train, res);}
/**
* @brief The rmes of the model
*/
inline double rmse(){return util_funcs::norm(_error.get(), _n_samp);}
inline double rmse(){return util_funcs::norm(_train_error.data(), _n_samp_train) / std::sqrt(static_cast<double>(_n_samp_train));}
inline double test_rmse(){return util_funcs::norm(_test_error.data(), _n_samp_test) / std::sqrt(static_cast<double>(_n_samp_test));}
/**
* @brief The max Absolute error of the array
*/
inline double max_ae()
{
return *std::max_element(_error.get(), _error.get() + _n_samp, [](double d1, double d2){return std::abs(d1) < std::abs(d2);});
return *std::max_element(_train_error.data(), _train_error.data() + _n_samp_train, [](double d1, double d2){return std::abs(d1) < std::abs(d2);});
}
inline double test_max_ae()
{
return *std::max_element(_test_error.data(), _test_error.data() + _n_samp_test, [](double d1, double d2){return std::abs(d1) < std::abs(d2);});
}
/**
* @brief Print model to a file
*/
void test_to_file(std::string filename);
/**
* @brief Print model to a file
*/
void train_to_file(std::string filename);
};
/**
......
#include <descriptor_identifier/SISSORegressor.hpp>
SISSORegressor::SISSORegressor(std::shared_ptr<FeatureSpace> feat_space, std::vector<double> prop, int n_dim):
SISSORegressor::SISSORegressor(std::shared_ptr<FeatureSpace> feat_space, std::vector<double> prop, std::vector<double> prop_test, int n_dim, int n_residual):
_feat_space(feat_space),
_mpi_comm(feat_space->mpi_comm()),
_n_samp(prop.size()),
_n_dim(n_dim),
_n_residual(n_residual),
_lwork(-1),
_rank(0),
_a(new double[(_n_dim + 1) * _n_samp]),
......@@ -14,7 +15,8 @@ SISSORegressor::SISSORegressor(std::shared_ptr<FeatureSpace> feat_space, std::ve
_work(nullptr),
_s(new double[n_dim + 1]),
_models(),
_prop(prop)
_prop(prop),
_prop_test(prop_test)
{
// Initialize a, b, ones, s, and _error arrays
......@@ -77,20 +79,28 @@ void SISSORegressor::set_error(std::vector<int>& inds, double* coeffs)
void SISSORegressor::fit()
{
std::vector<double> residual(_n_samp);
std::vector<double> residual(_n_samp * _n_residual);
std::clock_t start;
double duration;
start = std::clock();
_feat_space->sis(_prop);
_mpi_comm->barrier();
duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
if(_mpi_comm->rank() == 0)
std::cout << "Time for SIS: " << duration << std::endl;
start = std::clock();
std::vector<node_ptr> min_node(1, _feat_space->phi_selected()[0]);
_models.push_back(Model(_prop, min_node));
_models[_models.size() - 1].copy_error(residual.data());
std::vector<Model> models;
for(int rr = 0; rr < _n_residual; ++rr)
{
models.push_back(Model(_prop, _prop_test, {_feat_space->phi_selected()[rr]}));
models.back().copy_error(&residual[rr * _n_samp]);
}
_models.push_back(models);
_mpi_comm->barrier();
duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
if(_mpi_comm->rank() == 0)
......@@ -98,20 +108,24 @@ void SISSORegressor::fit()
for(int dd = 2; dd <= _n_dim; ++dd)
{
start = std::clock();
start = std::clock();
_feat_space->sis(residual);
_mpi_comm->barrier();
duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
_mpi_comm->barrier();
duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
if(_mpi_comm->rank() == 0)
std::cout << "Time for SIS: " << duration << std::endl;
start = std::clock();
l0_norm(residual, dd);
_mpi_comm->barrier();
l0_norm(_prop, dd);
_mpi_comm->barrier();
duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
if(_mpi_comm->rank() == 0)
std::cout << "Time for l0-norm: " << duration << std::endl;
_models[_models.size() - 1].copy_error(residual.data());
for(int rr = 0; rr < _n_residual; ++rr)
_models.back()[rr].copy_error(&residual[rr * _n_samp]);
}
}
......@@ -120,9 +134,9 @@ void SISSORegressor::l0_norm(std::vector<double>& prop, int n_dim)
std::vector<double> coefs(n_dim + 1, 0.0);
std::vector<int> inds(n_dim, 0);
std::vector<int> inds_min(n_dim);
std::vector<std::vector<int>> inds_min(_n_residual);
double min_error = util_funcs::norm(prop.data(), _n_samp);
std::vector<double> min_errors(_n_residual, util_funcs::norm(_prop.data(), _n_samp));
int n = _feat_space->phi_selected().size();
......@@ -152,28 +166,40 @@ void SISSORegressor::l0_norm(std::vector<double>& prop, int n_dim)
least_squares(inds, coefs.data());
set_error(inds, coefs.data());
if(util_funcs::norm(_error.get(), _n_samp) < min_error)
double error = util_funcs::norm(_error.get(), _n_samp);
if(error < min_errors.back())
{
std::copy_n(inds.begin(), inds.size(), inds_min.begin());
min_error = util_funcs::norm(_error.get(), _n_samp);
int rr = 0;
while((rr < _n_residual) && (error > min_errors[rr]))
++rr;
inds_min.insert(inds_min.begin() + rr, inds);
min_errors.insert(min_errors.begin() + rr, error);
inds_min.pop_back();
min_errors.pop_back();
}
++rr;
} while(std::next_permutation(v.begin(), v.end()) && (rr < start_end[1]));
std::vector<std::vector<double>> all_min_error;
std::vector<std::vector<std::vector<int>>> all_inds_min;
std::vector<double> all_min_error;
std::vector<std::vector<int>> all_inds_min;
mpi::all_gather(*_mpi_comm, min_error, all_min_error);
mpi::all_gather(*_mpi_comm, min_errors, all_min_error);
mpi::all_gather(*_mpi_comm, inds_min, all_inds_min);
int argmin = std::min_element(all_min_error.begin(), all_min_error.end()) - all_min_error.begin();
min_errors = std::vector<double>(all_min_error.size() * _n_residual);
for(int me = 0; me < all_min_error.size(); ++me)
std::copy_n(all_min_error[me].begin(), all_min_error[me].size(), &min_errors[me * _n_residual]);
inds = util_funcs::argsort(min_errors);
std::vector<node_ptr> min_nodes(n_dim);
for(int ii = 0; ii < n_dim; ++ii)
min_nodes[ii] = _feat_space->phi_selected()[all_inds_min[argmin][ii]];
std::vector<Model> models;
for(int rr = 0; rr < _n_residual; ++rr)
{
for(int ii = 0; ii < n_dim; ++ii)
min_nodes[ii] = _feat_space->phi_selected()[all_inds_min[inds[rr] / _n_residual][inds[rr] % _n_residual][ii]];
models.push_back(Model(_prop, _prop_test, min_nodes));
}
_models.push_back(Model(_prop, min_nodes));
_models.push_back(models);
}
......@@ -17,6 +17,7 @@ protected:
int _n_samp; //!< the number of samples per feature
int _n_dim; //!< Number of dimensions to calculate
int _n_residual; //!< Number of residuals to pass to the next sis model
int _lwork; //!< size of the work array
int _rank; //!< Ranks for the least squares problem
......@@ -27,8 +28,9 @@ protected:
std::unique_ptr<double[]> _work; //!< The work array for least squares problems
std::unique_ptr<double[]> _s; //!< The S array for least squares problems
std::vector<Model> _models; //!< List of models
std::vector<std::vector<Model>> _models; //!< List of models
std::vector<double> _prop; //!< Property array
std::vector<double> _prop_test; //!< Property array
public:
/**
......@@ -37,7 +39,7 @@ public:
* @param prop Property to model
* @param n_dim Maximum dimension of the model
*/
SISSORegressor(std::shared_ptr<FeatureSpace> feat_space, std::vector<double> prop, int n_dim);
SISSORegressor(std::shared_ptr<FeatureSpace> feat_space, std::vector<double> prop, std::vector<double> prop_test, int n_dim, int n_residual);
/**
* @brief Get the optimal size of the working array
......@@ -49,7 +51,7 @@ public:
/**
* @brief Preform Least squares optimization
\ *
*
* @param inds Feature indexes to get the model of
* @param coeffs Coefficients for the model
*/
......@@ -96,7 +98,7 @@ public:
/**
* @brief Acessor function for models
*/
inline std::vector<Model>& models(){return _models;}
inline std::vector<std::vector<Model>>& models(){return _models;}
/**
* @brief Acessor function for n_samp
......@@ -108,6 +110,11 @@ public:
*/
inline int n_dim(){return _n_dim;}
/**
* @brief Acessor function for n_residual
*/
inline int n_residual(){return _n_residual;}
/**
* @brief Acessor function for lwork
*/
......
......@@ -50,8 +50,8 @@ FeatureSpace::FeatureSpace(
_un_operators.push_back(allowed_op_maps::unary_operator_map[op]);
}
generate_feature_space();
_scores.resize(_phi.size());
_scores.reserve(_phi.size());
_scores.resize(_phi.size());
}
void FeatureSpace::generate_feature_space()
......@@ -139,13 +139,15 @@ void FeatureSpace::generate_feature_space()
++feat_ind;
}
}
if(nn <= _n_rung_store)
{
bool use_temp = (nn != _max_phi) || (_max_phi > _n_rung_store);
node_value_arrs::resize_values_arr(_n_rung_store, _phi.size(), use_temp);
for(int ff = _start_gen[0]; ff < _phi.size(); ++ff)
{
_phi[ff]->set_value();
_phi[ff]->set_test_value();
}
}
}
else
......@@ -261,30 +263,35 @@ void FeatureSpace::generate_feature_space()
{
_phi[ff]->reindex(ff + n_feat_below_rank, ff);
_phi[ff]->set_value();
_phi[ff]->set_test_value();
}
}
}
_n_feat = _phi.size();
}
void FeatureSpace::project_r(double* prop)
void FeatureSpace::project_r(double* prop, int size)
{
std::vector<double> scores(_phi.size(), 0.0);
for(int ff = 0; ff < _phi.size(); ++ff)
_scores[ff] = -1.0 * std::abs(util_funcs::r(prop, _phi[ff]->value_ptr(), _n_samp));
_scores[ff] = -1.0 * std::abs(util_funcs::r(&prop[0], _phi[ff]->value_ptr(), _n_samp));
for(int pp = 1; pp < size / _n_samp; ++pp)
{
for(int ff = 0; ff < _phi.size(); ++ff)
scores[ff] = -1.0 * std::abs(util_funcs::r(&prop[_n_samp*pp], _phi[ff]->value_ptr(), _n_samp));
std::transform(scores.begin(), scores.end(), _scores.begin(), _scores.begin(), [](double s1, double s2){return std::min(s1, s2);});
}
}
void FeatureSpace::sis(std::vector<double>& prop)
{
// while(true)
// {}
int cur_feat = node_value_arrs::N_SELECTED;
node_value_arrs::resize_d_matrix_arr(_n_sis_select);
_phi_selected.reserve(_phi_selected.size() + _n_sis_select);
project_r(prop.data());
project_r(prop.data(), prop.size());
std::vector<int> inds = util_funcs::argsort(_scores);
std::vector<double> scores_selected(_n_sis_select, 0.0);
......@@ -305,7 +312,7 @@ void FeatureSpace::sis(std::vector<double>& prop)
// Check the feature against those selected from previous SIS iterations
for(int dd = 0; dd < cur_feat; ++dd)
{
if(1.0 - std::abs(util_funcs::r(node_value_arrs::get_d_matrix_ptr(dd), _phi[inds[ii]]->value_ptr(), prop.size())) < 1e-13)
if(1.0 - std::abs(util_funcs::r(node_value_arrs::get_d_matrix_ptr(dd), _phi[inds[ii]]->value_ptr(), _n_samp)) < 1e-10)
{
is_valid = false;
break;
......@@ -320,7 +327,7 @@ void FeatureSpace::sis(std::vector<double>& prop)
if(*std::min_element(scores_comp.begin(), scores_comp.begin() + cur_feat_local) < 1e-10)
{
int dd = std::min_element(scores_comp.begin(), scores_comp.begin() + cur_feat_local) - scores_comp.begin();
if(1.0 - std::abs(util_funcs::r(node_value_arrs::get_d_matrix_ptr(cur_feat + dd), _phi[inds[ii]]->value_ptr(), prop.size())) < 1e-13)
if(1.0 - std::abs(util_funcs::r(node_value_arrs::get_d_matrix_ptr(cur_feat + dd), _phi[inds[ii]]->value_ptr(), _n_samp)) < 1e-10)
is_valid = false;
}
......@@ -328,7 +335,7 @@ void FeatureSpace::sis(std::vector<double>& prop)
{
inds_selected[cur_feat_local] = _phi[inds[ii]]->feat_ind();
scores_selected[cur_feat_local] = _scores[inds[ii]];
phi_selected.push_back(std::make_shared<FeatureNode>(cur_feat + cur_feat_local, _phi[inds[ii]]->expr(), _phi[inds[ii]]->value(), _phi[inds[ii]]->unit(), true));