jiri/LA_library
1
0 Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

bugfixes in reg and regsurf

Browse Source
This commit is contained in:
jiri committed 2026-10-29 22:48:21 +01:00
1 parent 8e82597927
commit 4f3ec9aef9
4 files changed
+135 -101
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

reg.cc
+33 -43
View File
@@ -70,50 +70,32 @@ const T* Matern_comb_number)
unsigned int i; unsigned int i;
unsigned int k; unsigned int k;
T part_1 = 0; // sum(x^2)^1/2 T part_1 = 0; // sum(x^2)^1/2
T exponential_part; T z;
T Matern_part = 0; T Matern_part;
T result;
for (i = 0; i < dimensions; i++) for (i = 0; i < dimensions; i++)
part_1 = part_1 + (x[i] * x[i]); part_1 = part_1 + (x[i] * x[i]);
part_1 = std::pow(part_1, 0.5); part_1 = std::sqrt(part_1);
exponential_part = std::exp(-part_1/(sigma)); z = 2 * part_1/sigma;
// Horner scheme for sum_{k=0}^{n} c_k z^(n-k)
Matern_part = Matern_comb_number[0];
for (k = 1; k <= n; k++)
Matern_part = Matern_part * z + Matern_comb_number[k];
for (k = 0; k < n; k++) return(std::exp(-part_1/sigma) * Matern_part);
Matern_part += Matern_comb_number[i] * std::pow(2 * part_1/sigma, n - k);
result = exponential_part * Matern_part;
return(result);
} }
template <typename T> template <typename T>
unsigned int Matern_coefficient(const unsigned int n, T* comb_number) // values of (n + k)!/(2 n)! * n!/(k!(n-k)!) unsigned int Matern_coefficient(const unsigned int n, T* comb_number) // values of (n + k)!/(2 n)! * n!/(k!(n-k)!), k = 0..n
{ {
unsigned int i , j;
unsigned int k; unsigned int k;
unsigned int part_1; // reverse value of (n + k)!/(2 n)!
unsigned int part_2; // n!/(n-k)!
unsigned int part_3; // k!
for (k = 0; k < n; k++) // downward recurrence c_n = 1, c_{k-1} = c_k k/((n+k)(n-k+1)), no factorials (no overflow)
{ comb_number[n] = 1;
part_1 = 1; for (k = n; k > 0; k--)
part_2 = 1; comb_number[k - 1] = comb_number[k] * T(k)/(T(n + k) * T(n - k + 1));
part_3 = 1;
for (i = n + k + 1; i <= 2 * n; i++)
part_1 = part_1 * i;
for (i = n - k + 1; i <= n; i++)
part_2 = part_2 * i;
for (i = 1; i < k; i++)
part_3 = part_3 * i;
comb_number[k] = T(part_2)/T(part_1 * part_3);
}
return(0); return(0);
} }
@@ -260,7 +242,7 @@ const string kernel)
{ {
for (i = 0; i < 4; i++) for (i = 0; i < 4; i++)
REG_params[i] = params[i]; REG_params[i] = params[i];
REG_Matern_comb_number_vec = NRVec<T>(REG_n); REG_Matern_comb_number_vec = NRVec<T>(REG_n + 1);
REG_Matern_comb_number = REG_Matern_comb_number_vec.operator T*(); REG_Matern_comb_number = REG_Matern_comb_number_vec.operator T*();
} }
else else
@@ -273,13 +255,15 @@ const string kernel)
} }
} }
// raw pointers to the data of the vectors below are cached for speed; the vectors must therefore never be
// reallocated: a |= b frees the old buffer of a (pointers would dangle), use a.storesubvector(0,b) (copy in place)
// copy data // copy data
REG_coordinates_vec = NRVec<T>(x, data_size * REG_count_dimensions); REG_coordinates_vec = NRVec<T>(x, data_size * REG_count_dimensions);
REG_coordinates = REG_coordinates_vec.operator T*(); REG_coordinates = REG_coordinates_vec.operator T*();
REG_y_vec = NRVec<T>(y , data_size); REG_y_vec = NRVec<T>(y , data_size);
REG_Y = REG_y_vec.operator T*(); REG_Y = REG_y_vec.operator T*();
// for fitting model // for fitting model
REG_distances_mat = NRMat<T>(data_size, data_size); REG_distances_mat = NRMat<T>(data_size, data_size*REG_count_dimensions); //difference vectors of all pairs
REG_distances = REG_distances_mat.operator T*(); REG_distances = REG_distances_mat.operator T*();
REG_kernels_mat = NRMat<T>(data_size, data_size); REG_kernels_mat = NRMat<T>(data_size, data_size);
REG_kernels = REG_kernels_mat.operator T*(); REG_kernels = REG_kernels_mat.operator T*();
@@ -288,7 +272,7 @@ const string kernel)
REG_weights_vec = NRVec<T>(data_size); REG_weights_vec = NRVec<T>(data_size);
REG_weights = REG_weights_vec.operator T*(); REG_weights = REG_weights_vec.operator T*();
// for predictions // for predictions
REG_coordinates_diff_vec = NRVec<T>(data_size); REG_coordinates_diff_vec = NRVec<T>(data_size*REG_count_dimensions);
REG_coordinates_diff = REG_coordinates_diff_vec.operator T*(); REG_coordinates_diff = REG_coordinates_diff_vec.operator T*();
REG_kernels_diff_vec = NRVec<T>(data_size); REG_kernels_diff_vec = NRVec<T>(data_size);
REG_kernels_diff = REG_kernels_diff_vec.operator T*(); REG_kernels_diff = REG_kernels_diff_vec.operator T*();
@@ -324,7 +308,7 @@ T REG<T>::Kernel(const T* x, const unsigned int dimensions, const unsigned int k
result = Laplacian_kernel(x, dimensions, REG_sigma); result = Laplacian_kernel(x, dimensions, REG_sigma);
break; break;
case 2: case 2:
result = Laplacian_kernel(x, dimensions, REG_sigma); result = Exponential_kernel(x, dimensions, REG_sigma);
break; break;
case 3: case 3:
result = Matern_kernel(x, dimensions, REG_sigma, REG_n, REG_Matern_comb_number); result = Matern_kernel(x, dimensions, REG_sigma, REG_n, REG_Matern_comb_number);
@@ -443,7 +427,7 @@ T REG<T>::REG_RMSE()
T y_max = REG_y_vec.amax(); T y_max = REG_y_vec.amax();
T y_min = REG_y_vec.amin(); T y_min = REG_y_vec.amin();
REG_y_delta_vec |= REG_y_preds_vec; REG_y_delta_vec.storesubvector(0, REG_y_preds_vec); //in place, not |= (cached raw pointer)
REG_y_delta_vec -= REG_y_vec; REG_y_delta_vec -= REG_y_vec;
REG_y_delta_vec *= REG_y_delta_vec; REG_y_delta_vec *= REG_y_delta_vec;
delta_square = REG_y_delta_vec.sum(); delta_square = REG_y_delta_vec.sum();
@@ -509,7 +493,7 @@ unsigned int REG<T>::REG_Fit(const T max_Loss, const unsigned int max_learning_c
REG_mat_init = true; REG_mat_init = true;
} }
// initialitation of weight coefficients vector // initialitation of weight coefficients vector
REG_weights_vec |= REG_y_vec; REG_weights_vec.storesubvector(0, REG_y_vec); //in place, not |= (cached raw pointer)
// Linear learning cycle of weight coefficient for learning // Linear learning cycle of weight coefficient for learning
for (i = 0; i < max_learning_cycles; i++) for (i = 0; i < max_learning_cycles; i++)
{ {
@@ -523,7 +507,7 @@ unsigned int REG<T>::REG_Fit(const T max_Loss, const unsigned int max_learning_c
REG_y_preds_vec += Lasso_part; REG_y_preds_vec += Lasso_part;
// calculating new weights // calculating new weights
REG_y_temp_vec |= REG_y_vec; REG_y_temp_vec.storesubvector(0, REG_y_vec); //in place, not |= (cached raw pointer)
REG_y_temp_vec /= REG_y_preds_vec; REG_y_temp_vec /= REG_y_preds_vec;
REG_weights_vec *= REG_y_temp_vec; REG_weights_vec *= REG_y_temp_vec;
// checking the Loss_Function value // checking the Loss_Function value
@@ -580,8 +564,8 @@ unsigned int REG<T>::REG_Fit_grad(const T max_Loss, const unsigned int max_learn
REG_mat_init = true; REG_mat_init = true;
} }
// initialitation of weight coefficients vector // initialitation of weight coefficients vector
REG_weights_vec |= REG_y_vec; REG_weights_vec.storesubvector(0, REG_y_vec); //in place, not |= (cached raw pointer)
kernel_ridge_part = REG_lambda * (1.00 - REG_alpha); kernel_ridge_part = REG_lambda * (1.00 - REG_alpha) * 0.5;
Lasso_part = REG_lambda * REG_alpha; Lasso_part = REG_lambda * REG_alpha;
sqrt_count_x = std::sqrt(count_x); sqrt_count_x = std::sqrt(count_x);
// gradient descent Learning cycle of weight coefficient for learning // gradient descent Learning cycle of weight coefficient for learning
@@ -760,6 +744,9 @@ unsigned int REG<T>::REG_Get_Score(const T* x, T* y_preds, const unsigned int d
template <typename T> template <typename T>
unsigned int REG<T>::REG_Set_params(const T* params, const string kernel) unsigned int REG<T>::REG_Set_params(const T* params, const string kernel)
{ {
const unsigned int old_kernel_id = REG_kernel_id;
const unsigned int old_n = (REG_kernel_id == 3) ? REG_n : 0; // REG_n is set only for Matern
const T old_sigma = REG_sigma;
// setting type of kernel // setting type of kernel
REG_kernel_type = kernel; REG_kernel_type = kernel;
if (kernel == "Gaussian") if (kernel == "Gaussian")
@@ -784,10 +771,13 @@ unsigned int REG<T>::REG_Set_params(const T* params, const string kernel)
// calculating Matern coefficients, if Matern kernel is selected // calculating Matern coefficients, if Matern kernel is selected
if (REG_kernel_id == 3) if (REG_kernel_id == 3)
{ {
REG_Matern_comb_number_vec.resize(REG_n); REG_Matern_comb_number_vec.resize(REG_n + 1);
REG_Matern_comb_number = REG_Matern_comb_number_vec.operator T*(); REG_Matern_comb_number = REG_Matern_comb_number_vec.operator T*();
Matern_coefficient(REG_n, REG_Matern_comb_number); Matern_coefficient(REG_n, REG_Matern_comb_number);
} }
// kernel matrices depend on kernel, sigma and n (not on lambda, alpha): rebuild them in the next fit
if (REG_kernel_id != old_kernel_id || REG_sigma != old_sigma || (REG_kernel_id == 3 && REG_n != old_n))
REG_mat_init = false;
return(0); return(0);
} }
@@ -800,7 +790,7 @@ unsigned int REG<T>::REG_Get_params(T* params, string* kernel)
params[1] = REG_alpha; params[1] = REG_alpha;
params[2] = REG_sigma; params[2] = REG_sigma;
if (REG_kernel_type == "Mattern") if (REG_kernel_id == 3)
params[3] = REG_n; params[3] = REG_n;
return(0); return(0);
reg.h
+1 -1
View File
@@ -59,7 +59,7 @@ public:
// Init function for constructor // Init function for constructor
void Init(const T* x, const T* y, const T* params, const unsigned int data_size, const unsigned int x_dimensions, const string kernel); void Init(const T* x, const T* y, const T* params, const unsigned int data_size, const unsigned int x_dimensions, const string kernel);
// General kernel function // General kernel function
inline T Kernel(const T* x, const unsigned int dimensions, const unsigned int kernel_id); T Kernel(const T* x, const unsigned int dimensions, const unsigned int kernel_id);
// Function, that build matrix of distances between kernels // Function, that build matrix of distances between kernels
unsigned int Build_REG_distances_mat(); unsigned int Build_REG_distances_mat();
// Function, that build matrix of coefficients for kernel coordinates // Function, that build matrix of coefficients for kernel coordinates
regsurf.cc
+99 -55
View File
@@ -229,14 +229,12 @@ REGSURF<T>::REGSURF(const string geometries_file, const string y_file, const str
if(!s_geom_file) laerror("allocation error of s_geom_file"); if(!s_geom_file) laerror("allocation error of s_geom_file");
s_y_file = new T[count_geometries * count_levels]; s_y_file = new T[count_geometries * count_levels];
if(!s_y_file) laerror("allocation error of s_y_file"); if(!s_y_file) laerror("allocation error of s_y_file");
if (descriptor == "Coulomb_matrix" or Z_file != "") if (Z_file != "")
{ {
s_Z_file = new unsigned int[count_kernels]; s_Z_file = new unsigned int[count_kernels];
if(!s_Z_file) laerror("allocation error of s_Z_file"); if(!s_Z_file) laerror("allocation error of s_Z_file");
} }
else s_Z_file=NULL; else s_Z_file=NULL;
s_y_min = new unsigned int[count_levels];
if(!s_y_min) laerror("allocation error of s_y_min");
} }
// loading geometrical data for regression model from file // loading geometrical data for regression model from file
l = 0; l = 0;
@@ -349,6 +347,15 @@ REGSURF<T>::REGSURF(const string geometries_file, const string y_file, const str
// callling init function with loaded parameters, y and geometries data // callling init function with loaded parameters, y and geometries data
Init(count_dimensions, count_kernels, count_geometries, count_levels, descriptor, kernel, s_param_file, s_geom_file, s_y_file, Init(count_dimensions, count_kernels, count_geometries, count_levels, descriptor, kernel, s_param_file, s_geom_file, s_y_file,
s_Z_file); s_Z_file);
// Init copied the data, the file buffers are not needed any more
delete[] s_param_file;
s_param_file = nullptr;
delete[] s_geom_file;
s_geom_file = nullptr;
delete[] s_y_file;
s_y_file = nullptr;
delete[] s_Z_file;
s_Z_file = nullptr;
} }
@@ -376,6 +383,16 @@ void REGSURF<T>::Init(const unsigned int count_dimensions, const unsigned int co
s_kernel_id = 3; s_kernel_id = 3;
s_descriptor_type = descriptor; s_descriptor_type = descriptor;
if (descriptor != "RE" and descriptor != "Coulomb_matrix" and descriptor != "Inverse_pairs")
{
laerror("invalid type of descriptor");
return;
}
if (descriptor == "Coulomb_matrix" and Z == nullptr)
{
laerror("Coulomb_matrix descriptor needs Z");
return;
}
if (descriptor == "RE") if (descriptor == "RE")
s_descriptor_id = 0; s_descriptor_id = 0;
if (descriptor == "Coulomb_matrix") if (descriptor == "Coulomb_matrix")
@@ -392,7 +409,8 @@ void REGSURF<T>::Init(const unsigned int count_dimensions, const unsigned int co
{ {
s_lambda = new T[count_levels]; s_lambda = new T[count_levels];
s_alpha = new T[count_levels]; s_alpha = new T[count_levels];
if(!s_lambda || !s_alpha) laerror("allocation error in REGSURF::Init"); s_y_min = new unsigned int[count_levels];
if(!s_lambda || !s_alpha || !s_y_min) laerror("allocation error in REGSURF::Init");
} }
for (i = 0; i < count_levels; i++) for (i = 0; i < count_levels; i++)
s_lambda[i] = params[i]; s_lambda[i] = params[i];
@@ -400,10 +418,10 @@ void REGSURF<T>::Init(const unsigned int count_dimensions, const unsigned int co
for (i = count_levels; i < (2 * count_levels); i++) for (i = count_levels; i < (2 * count_levels); i++)
{ {
s_alpha[i - count_levels] = params[i]; s_alpha[i - count_levels] = params[i];
if (s_alpha[i] > 1) if (s_alpha[i - count_levels] > 1)
s_alpha[i] = 1; s_alpha[i - count_levels] = 1;
if (s_alpha[i] < 0) if (s_alpha[i - count_levels] < 0)
s_alpha[i] = 0; s_alpha[i - count_levels] = 0;
} }
s_sigma = params[(2 * count_levels)]; s_sigma = params[(2 * count_levels)];
if (kernel == "Matern") if (kernel == "Matern")
@@ -411,19 +429,18 @@ void REGSURF<T>::Init(const unsigned int count_dimensions, const unsigned int co
s_aloc_init = false; s_aloc_init = false;
{ // for loaded parameters { // for loaded parameters
if (s_kernel_id == 3) // always room for n, so that s_Set_params can switch to Matern
{
s_params = new T[(2 * count_levels) + 2]; s_params = new T[(2 * count_levels) + 2];
if(!s_params) laerror("allocation error in REGSURF::Init"); if(!s_params) laerror("allocation error in REGSURF::Init");
s_Matern_comb_number = new T[s_n]; if (s_kernel_id == 3)
{
s_Matern_comb_number = new T[s_n + 1];
if(!s_Matern_comb_number) laerror("allocation error in REGSURF::Init"); if(!s_Matern_comb_number) laerror("allocation error in REGSURF::Init");
for (i = 0; i < ((2 * count_levels) + 2); i++) for (i = 0; i < ((2 * count_levels) + 2); i++)
s_params[i] = params[i]; s_params[i] = params[i];
} }
else else
{ {
s_params = new T[(2 * count_levels) + 1];
if(!s_params) laerror("allocation error in REGSURF::Init");
for (i = 0; i < ((2 * count_levels) + 1); i++) for (i = 0; i < ((2 * count_levels) + 1); i++)
s_params[i] = params[i]; s_params[i] = params[i];
} }
@@ -435,9 +452,9 @@ void REGSURF<T>::Init(const unsigned int count_dimensions, const unsigned int co
min_value = values[i * count_geometries]; min_value = values[i * count_geometries];
s_y_min[i] = 0; s_y_min[i] = 0;
for (j = 1; j < count_geometries; j++) for (j = 1; j < count_geometries; j++)
if (values[(i * count_kernels) + j] < min_value) if (values[(i * count_geometries) + j] < min_value)
{ {
min_value = values[(i * count_kernels) + j]; min_value = values[(i * count_geometries) + j];
s_y_min[i] = j; s_y_min[i] = j;
} }
} }
@@ -447,7 +464,10 @@ void REGSURF<T>::Init(const unsigned int count_dimensions, const unsigned int co
s_y_mat = NRMat<T>(values, count_levels, count_geometries); s_y_mat = NRMat<T>(values, count_levels, count_geometries);
// NRVec for vector of Z for coulomb matrix descriptor [s_count_kernels] // NRVec for vector of Z for coulomb matrix descriptor [s_count_kernels]
if (Z != nullptr) if (Z != nullptr)
{
s_Z_vec = NRVec<unsigned int>(Z, count_kernels); s_Z_vec = NRVec<unsigned int>(Z, count_kernels);
s_Z = s_Z_vec.operator unsigned int*();
}
// NRMat for coordinate distances for geometries [s_count_kernels, s_count_kernels * s_count_geometries] // NRMat for coordinate distances for geometries [s_count_kernels, s_count_kernels * s_count_geometries]
s_distances_mat = NRMat<T>(count_kernels, count_kernels * count_geometries); s_distances_mat = NRMat<T>(count_kernels, count_kernels * count_geometries);
@@ -494,11 +514,11 @@ void REGSURF<T>::Init(const unsigned int count_dimensions, const unsigned int co
s_kernels_pred_vec_mat = NRVec<NRMat<T>>(1); s_kernels_pred_vec_mat = NRVec<NRMat<T>>(1);
s_kernels_pred_vec_mat[0] = NRMat<T>(s_count_kernels, s_count_kernels); s_kernels_pred_vec_mat[0] = NRMat<T>(s_count_kernels, s_count_kernels);
// raw pointers to the data of the vectors below are cached for speed; the vectors must therefore never be
// reallocated: a |= b frees the old buffer of a (pointers would dangle), use a.storesubvector(0,b) (copy in place)
// for loaded geometries and surfaces // for loaded geometries and surfaces
s_geometries = s_geometries_mat.operator T*(); s_geometries = s_geometries_mat.operator T*();
s_y = s_y_mat.operator T*(); s_y = s_y_mat.operator T*();
// NRVec for vector of Z for coulomb matrix descriptor
s_Z = s_Z_vec.operator unsigned int*();
// for computing weights coefficiensts of kernels and making predictions // for computing weights coefficiensts of kernels and making predictions
s_distances = s_distances_mat.operator T*(); s_distances = s_distances_mat.operator T*();
s_desc = s_desc_mat.operator T*(); s_desc = s_desc_mat.operator T*();
@@ -539,7 +559,7 @@ T REGSURF<T>::Kernel(const T* x, const unsigned int dimensions, const unsigned i
result = Laplacian_kernel(x, dimensions, s_sigma); result = Laplacian_kernel(x, dimensions, s_sigma);
break; break;
case 2: case 2:
result = Laplacian_kernel(x, dimensions, s_sigma); result = Exponential_kernel(x, dimensions, s_sigma);
break; break;
case 3: case 3:
result = Matern_kernel(x, dimensions, s_sigma, s_n, s_Matern_comb_number); result = Matern_kernel(x, dimensions, s_sigma, s_n, s_Matern_comb_number);
@@ -609,7 +629,8 @@ unsigned int REGSURF<T>::Build_s_descriptor_mat(const unsigned int count_geometr
case 3: case 3:
extern_desc(s_count_kernels, count_geometries, distances, desc); extern_desc(s_count_kernels, count_geometries, distances, desc);
break; break;
laerror("unknown decsriptor identification number"); default:
laerror("unknown descriptor identification number");
return(1); return(1);
} }
} }
@@ -644,8 +665,8 @@ unsigned int REGSURF<T>::Build_s_kernels_mat(const unsigned int count_geometries
if (desc[(i * count_kernels * count_kernels) + j * (count_kernels + 1)] != 0) if (desc[(i * count_kernels * count_kernels) + j * (count_kernels + 1)] != 0)
{ {
descriptor_difference = desc[(i * count_kernels * count_kernels) + j * (count_kernels + 1)] - descriptor_difference = desc[(i * count_kernels * count_kernels) + j * (count_kernels + 1)] -
s_desc[(i * count_kernels * count_kernels) + j * (count_kernels + 1)]; s_desc[j * (count_kernels + 1)];
matrix[j * (count_kernels + 1)] = Kernel(&descriptor_difference, 1, s_kernel_id);; matrix[j * (count_kernels + 1)] = Kernel(&descriptor_difference, 1, s_kernel_id);
} }
else else
matrix[j * (count_kernels + 1)] = 0; matrix[j * (count_kernels + 1)] = 0;
@@ -701,7 +722,7 @@ unsigned int REGSURF<T>::Build_s_kernels_mat_sums()
temp[j] = sum; temp[j] = sum;
} }
s_kernels_diff_vec_vec[i] |= s_kernels_rows_sums_vec_vec[0]; s_kernels_diff_vec_vec[i].storesubvector(0, s_kernels_rows_sums_vec_vec[0]); //in place, not |= (cached raw pointer)
s_kernels_diff_vec_vec[i] -= s_y_temp_vec; s_kernels_diff_vec_vec[i] -= s_y_temp_vec;
} }
// create vector of kernel sums differences for each geometry // create vector of kernel sums differences for each geometry
@@ -747,7 +768,7 @@ T REGSURF<T>::s_RMSE(const unsigned int level)
// predictions for all geometries // predictions for all geometries
for (i = 0; i < s_count_geometries; i++) for (i = 0; i < s_count_geometries; i++)
{ {
s_y_temp_part_vec |= s_kernels_rows_sums_vec_vec[i]; s_y_temp_part_vec.storesubvector(0, s_kernels_rows_sums_vec_vec[i]); //in place, not |= (cached raw pointer)
s_y_temp_part_vec *= s_weights_vec_vec[level]; s_y_temp_part_vec *= s_weights_vec_vec[level];
y_pred = s_y_temp_part_vec.sum(); y_pred = s_y_temp_part_vec.sum();
if (kernel_ridge_part != 0) if (kernel_ridge_part != 0)
@@ -758,7 +779,7 @@ T REGSURF<T>::s_RMSE(const unsigned int level)
s_y_preds[i] = y_pred; s_y_preds[i] = y_pred;
} }
s_y_delta_vec |= s_y_preds_vec; s_y_delta_vec.storesubvector(0, s_y_preds_vec); //in place, not |= (cached raw pointer)
// subtracting real values // subtracting real values
s_y_delta_vec -= s_y_mat.row(level); s_y_delta_vec -= s_y_mat.row(level);
// the second power // the second power
@@ -858,7 +879,7 @@ const bool element_invariant, unsigned int* kernel_classes)
Lasso_part = s_lambda[i] * s_alpha[i]; Lasso_part = s_lambda[i] * s_alpha[i];
kernel_ridge_part = s_lambda[i] * (1 - s_alpha[i]) * 0.5 * count_kernels; kernel_ridge_part = s_lambda[i] * (1 - s_alpha[i]) * 0.5 * count_kernels;
// initialitation of weight coefficients vector // initialitation of weight coefficients vector
init_weight = (s_y[i] - Lasso_part)/(s_kernels_sum_sum + kernel_ridge_part); init_weight = (s_y[i * count_geometries] - Lasso_part)/(s_kernels_sum_sum + kernel_ridge_part);
for (j = 0; j < count_kernels; j++) for (j = 0; j < count_kernels; j++)
s_weights_vec_vec[i][j] = init_weight; s_weights_vec_vec[i][j] = init_weight;
@@ -867,7 +888,7 @@ const bool element_invariant, unsigned int* kernel_classes)
for (k = 1; k < count_geometries; k++) for (k = 1; k < count_geometries; k++)
{ {
// get the model prediction for another geometry // get the model prediction for another geometry
s_y_temp_part_vec |= s_kernels_rows_sums_vec_vec[k]; s_y_temp_part_vec.storesubvector(0, s_kernels_rows_sums_vec_vec[k]); //in place, not |= (cached raw pointer)
s_y_temp_part_vec *= s_weights_vec_vec[i]; s_y_temp_part_vec *= s_weights_vec_vec[i];
y_pred = s_y_temp_part_vec.sum(); y_pred = s_y_temp_part_vec.sum();
// calculating kernel ridge part of prediction // calculating kernel ridge part of prediction
@@ -880,12 +901,12 @@ const bool element_invariant, unsigned int* kernel_classes)
// calculating difference // calculating difference
y_diff = y_pred_3 - s_y[k + (i * count_geometries)]; y_diff = y_pred_3 - s_y[k + (i * count_geometries)];
// calculating new weights // calculating new weights
s_y_temp_part_vec |= s_kernels_diff_vec_vec[k]; s_y_temp_part_vec.storesubvector(0, s_kernels_diff_vec_vec[k]); //in place, not |= (cached raw pointer)
s_y_temp_part_vec /= s_kernels_diff_sums_vec[k]; s_y_temp_part_vec /= s_kernels_diff_sums_vec[k];
// optimalization of kernel ridge regression // optimalization of kernel ridge regression
if (kernel_ridge_part != 0) if (kernel_ridge_part != 0)
{ {
s_y_temp_vec |= s_weights_vec_vec[i]; s_y_temp_vec.storesubvector(0, s_weights_vec_vec[i]); //in place, not |= (cached raw pointer)
s_y_temp_vec /= s_weights_vec_vec[i].sum()/y_pred_2; s_y_temp_vec /= s_weights_vec_vec[i].sum()/y_pred_2;
s_y_temp_part_vec -= s_y_temp_vec; s_y_temp_part_vec -= s_y_temp_vec;
} }
@@ -902,7 +923,7 @@ const bool element_invariant, unsigned int* kernel_classes)
break; break;
if (kernel_classes_index[m] == kernel_classes_index[l] and m >= l) if (kernel_classes_index[m] == kernel_classes_index[l] and m >= l)
{ {
element_sum += s_y_temp_part[l]; element_sum += s_y_temp_part[m];
element_count += 1; element_count += 1;
} }
} }
@@ -918,7 +939,7 @@ const bool element_invariant, unsigned int* kernel_classes)
} }
} }
} }
s_y_temp_part_vec /= (T(count_geometries)/(learning_rate * y_diff)); s_y_temp_part_vec *= learning_rate * y_diff/T(count_geometries);
s_weights_vec_vec[i] -= s_y_temp_part_vec; s_weights_vec_vec[i] -= s_y_temp_part_vec;
} }
// convergence criteria check // convergence criteria check
@@ -1066,7 +1087,7 @@ const T* geom, const unsigned int* surfaces, T* predictions)
// getting prediction // getting prediction
if (surfaces[j] < s_count_levels) if (surfaces[j] < s_count_levels)
{ {
s_y_temp_part_vec |= s_y_temp_vec; s_y_temp_part_vec.storesubvector(0, s_y_temp_vec); //in place, not |= (cached raw pointer)
s_y_temp_part_vec *= s_weights_vec_vec[surfaces[j]]; s_y_temp_part_vec *= s_weights_vec_vec[surfaces[j]];
} }
else else
@@ -1093,10 +1114,10 @@ const T* geom, const unsigned int* surfaces, T* predictions)
template <typename T> template <typename T>
unsigned int REGSURF<T>::s_Score(const unsigned int count_geometries, const unsigned int count_levels, unsigned int REGSURF<T>::s_Score(const unsigned int count_geometries, const unsigned int count_levels,
const T* geom, const unsigned int* surfaces, T* predictions, T* y_values, T* score) const T* geom, const unsigned int* surfaces, T* predictions, const T* y_values, T* score)
{ {
unsigned int i; unsigned int i;
T RMSE; T RMSE = 0;
s_Get_predictions(count_geometries, count_levels, geom, surfaces, predictions); s_Get_predictions(count_geometries, count_levels, geom, surfaces, predictions);
for (i = 0; i < count_geometries * count_levels; i++) for (i = 0; i < count_geometries * count_levels; i++)
@@ -1140,6 +1161,16 @@ unsigned int REGSURF<T>::s_Set_params(const string descriptor, const string kern
new_kernel_id = 3; new_kernel_id = 3;
s_descriptor_type = descriptor; s_descriptor_type = descriptor;
if (descriptor != "RE" and descriptor != "Coulomb_matrix" and descriptor != "Inverse_pairs")
{
laerror("invalid type of descriptor");
return(1);
}
if (descriptor == "Coulomb_matrix" and s_Z == nullptr)
{
laerror("Coulomb_matrix descriptor needs Z");
return(1);
}
if (descriptor == "RE") if (descriptor == "RE")
new_descriptor_id = 0; new_descriptor_id = 0;
if (descriptor == "Coulomb_matrix") if (descriptor == "Coulomb_matrix")
@@ -1171,6 +1202,9 @@ unsigned int REGSURF<T>::s_Set_params(const string descriptor, const string kern
if (s_alpha[i - count_levels] < 0) if (s_alpha[i - count_levels] < 0)
s_alpha[i - count_levels] = 0; s_alpha[i - count_levels] = 0;
} }
// kernel matrices depend on sigma and n (not on lambda, alpha): rebuild them in the next fit
if (params[2 * count_levels] != s_sigma)
s_kernels_mat_init = false;
s_sigma = params[2 * count_levels]; s_sigma = params[2 * count_levels];
if (s_kernel_id == 3) if (s_kernel_id == 3)
@@ -1182,6 +1216,9 @@ unsigned int REGSURF<T>::s_Set_params(const string descriptor, const string kern
if (s_kernel_id == 3) if (s_kernel_id == 3)
{ {
// s_n is read only if the kernel was Matern before (otherwise s_kernels_mat_init is already false)
if (s_kernels_mat_init and static_cast<unsigned int>(params[(2 * count_levels) + 1] + 0.5) != s_n)
s_kernels_mat_init = false;
s_n = static_cast<unsigned int>(params[(2 * count_levels) + 1] + 0.5); s_n = static_cast<unsigned int>(params[(2 * count_levels) + 1] + 0.5);
if (s_Matern_comb_number != nullptr) if (s_Matern_comb_number != nullptr)
{ {
@@ -1189,7 +1226,7 @@ unsigned int REGSURF<T>::s_Set_params(const string descriptor, const string kern
s_Matern_comb_number = nullptr; s_Matern_comb_number = nullptr;
} }
{ {
s_Matern_comb_number = new T[s_n]; s_Matern_comb_number = new T[s_n + 1];
if(!s_Matern_comb_number) laerror("allocation error of s_Matern_comb_number"); if(!s_Matern_comb_number) laerror("allocation error of s_Matern_comb_number");
Matern_coefficient(s_n, s_Matern_comb_number); Matern_coefficient(s_n, s_Matern_comb_number);
} }
@@ -1381,25 +1418,25 @@ unsigned int REGSURF<T>::s_Predict_ML_NEA_geometries()
template <typename T> template <typename T>
unsigned int REGSURF<T>::s_Compute_ML_NEA_spectra(const T E_min, const T step, const unsigned int count_steps, const T delta) unsigned int REGSURF<T>::s_Compute_ML_NEA_spectra(const T E_min, const T step, const unsigned int count_steps, const T delta)
{ // sigma(E) = (e^2 *h)/(4*m_e*c*E_0*E) * sum_{n=1}^N\(1/Np * sum_{i=1}^Np * delta_E_0n(x_i) * f_0n(x_i) * f_g(E - delta_E_0n, delta)) { // sigma(E) = (e^2 *h)/(4*m_e*c*E_0*E) * sum_{n=1}^N\(1/Np * sum_{i=1}^Np * delta_E_0n(x_i) * f_0n(x_i) * f_g(E - delta_E_0n, delta))
// delta_E_0n(x_i) = delta_E_min_n/delta_E_0n // f_0n(x_i) = delta_E_min_n/delta_E_0n(x_i) is a placeholder for the true oscillator strength (not predicted by the model)
// f_g = 1/(sqrt(Phi) * delta * delta) * exp(-((E-delta_E_0n)^2)/(0.5 * delta * delta)) // f_g = sqrt(2/pi)/delta * exp(-2*(E-delta_E_0n)^2/(delta * delta)) (normalized Gaussian, standard deviation delta/2)
unsigned int i, j, k; unsigned int i, j;
unsigned int count_kernels, count_levels, count_geometries; unsigned int count_levels, count_geometries;
long k, k_lo, k_hi;
T y_min; T y_min;
T h = 6.62607015E-34; T h = 6.62607015E-34;
T e = 1.602176634E-19; T e = 1.602176634E-19;
T m_e = 9.1093837E-31; T m_e = 9.1093837E-31;
T c = 299796000; T c = 299792458;
T E_0 = 8.8541878128E-12; T E_0 = 8.8541878128E-12;
T Phi = (1.00 + std::sqrt(5))/2.00;
T pre_sum_part; T pre_sum_part;
T f_g_norm, f_g_exp, cutoff;
T E; T E;
T delta_E_0n; T delta_E_0n;
T f_0n; T f_0n;
T f_g; T f_g;
T sigma_part; T sigma_part;
count_kernels = s_count_kernels;
count_levels = s_count_levels; count_levels = s_count_levels;
count_geometries = ML_NEA_count_geometries; count_geometries = ML_NEA_count_geometries;
if (s_base_init == false) if (s_base_init == false)
@@ -1413,16 +1450,16 @@ unsigned int REGSURF<T>::s_Compute_ML_NEA_spectra(const T E_min, const T step, c
return(1); return(1);
} }
// calculating minima of surfaces // calculating minima of surfaces, s_ML_NEA_y is [geometry * count_levels + level] as stored by s_Get_predictions
s_ML_NEA_min_y_vec = NRVec<T>(s_count_levels); s_ML_NEA_min_y_vec = NRVec<T>(s_count_levels);
s_ML_NEA_min_y = s_ML_NEA_min_y_vec.operator T*(); s_ML_NEA_min_y = s_ML_NEA_min_y_vec.operator T*();
for (i = 0; i < count_levels; i++) for (i = 0; i < count_levels; i++)
{ {
y_min = s_ML_NEA_y[i * count_kernels]; y_min = s_ML_NEA_y[i];
for (j = 1; j < count_kernels; j++) for (j = 1; j < count_geometries; j++)
{ {
if (s_ML_NEA_y[(i * count_kernels) + j] < y_min) if (s_ML_NEA_y[(j * count_levels) + i] < y_min)
y_min = s_ML_NEA_y[(i * count_kernels) + j]; y_min = s_ML_NEA_y[(j * count_levels) + i];
} }
s_ML_NEA_min_y[i] = y_min; s_ML_NEA_min_y[i] = y_min;
} }
@@ -1430,7 +1467,10 @@ unsigned int REGSURF<T>::s_Compute_ML_NEA_spectra(const T E_min, const T step, c
s_ML_NEA_spectra_vec = NRVec<T>(count_steps); s_ML_NEA_spectra_vec = NRVec<T>(count_steps);
s_ML_NEA_spectra = s_ML_NEA_spectra_vec.operator T*(); s_ML_NEA_spectra = s_ML_NEA_spectra_vec.operator T*();
pre_sum_part = (e * e * h)/(4 * m_e * c * E_0); pre_sum_part = (e * e * h)/(4 * m_e * c * E_0);
E = E_min; f_g_norm = std::sqrt(2/M_PI)/delta;
f_g_exp = 2/(delta * delta);
// exp(-f_g_exp*x^2) underflows to 0 for |x| > 19.3 delta: grid points outside 20 delta are skipped
cutoff = 20 * delta;
for (i = 0; i < count_steps; i++) for (i = 0; i < count_steps; i++)
s_ML_NEA_spectra[i] = 0; s_ML_NEA_spectra[i] = 0;
@@ -1439,24 +1479,28 @@ unsigned int REGSURF<T>::s_Compute_ML_NEA_spectra(const T E_min, const T step, c
for (j = 0; j < count_levels; j++) for (j = 0; j < count_levels; j++)
{ {
// calculating the oscilator strenght // calculating the oscilator strenght
delta_E_0n = s_ML_NEA_y[(j * count_geometries) + i]; delta_E_0n = s_ML_NEA_y[(i * count_levels) + j];
f_0n = s_ML_NEA_min_y[j]/delta_E_0n; f_0n = s_ML_NEA_min_y[j]/delta_E_0n;
E = E_min;
if (delta_E_0n > 0 and f_0n > 0) if (delta_E_0n > 0 and f_0n > 0)
for (k = 0; k < count_steps; k++) {
// grid points E_k = E_min + k*step within the cutoff
k_lo = static_cast<long>(std::ceil((delta_E_0n - cutoff - E_min)/step));
k_hi = static_cast<long>(std::floor((delta_E_0n + cutoff - E_min)/step));
if (k_lo < 0)
k_lo = 0;
if (k_hi >= static_cast<long>(count_steps))
k_hi = static_cast<long>(count_steps) - 1;
for (k = k_lo; k <= k_hi; k++)
{ {
// calculating gaussian function f_g partial intensity and adding to ML_NEA intensities vector // calculating gaussian function f_g partial intensity and adding to ML_NEA intensities vector
if (abs(E - delta_E_0n)/delta < 100) E = E_min + k * step;
{ f_g = f_g_norm * std::exp(-f_g_exp * (E - delta_E_0n)*(E - delta_E_0n));
f_g = 1/(std::sqrt(Phi) * delta);
f_g *= std::exp(-((E - delta_E_0n)*(E - delta_E_0n))/(0.5 * delta * delta));
if (E > 0 and f_g > 0) if (E > 0 and f_g > 0)
{ {
sigma_part = pre_sum_part/E * delta_E_0n * f_0n * f_g; sigma_part = pre_sum_part/E * delta_E_0n * f_0n * f_g;
s_ML_NEA_spectra[k] += sigma_part; s_ML_NEA_spectra[k] += sigma_part;
} }
} }
E += step;
} }
} }
} }
regsurf.h
+2 -2
View File
@@ -59,7 +59,7 @@ public:
const unsigned int count_levels, const string descriptor, const string kernel, const T* params, const T* geometries, const unsigned int count_levels, const string descriptor, const string kernel, const T* params, const T* geometries,
const T* values, const unsigned int* Z); const T* values, const unsigned int* Z);
// General kernel function // General kernel function
inline T Kernel(const T* x, const unsigned int dimensions, const unsigned int kernel_id); T Kernel(const T* x, const unsigned int dimensions, const unsigned int kernel_id);
// function for building matrix of geometry distances between kernels // function for building matrix of geometry distances between kernels
unsigned int Build_s_distances_mat(const unsigned int count_geometries, const T* geometries, T* distances); unsigned int Build_s_distances_mat(const unsigned int count_geometries, const T* geometries, T* distances);
// function for building descriptor matrix // function for building descriptor matrix
@@ -84,7 +84,7 @@ public:
const T* geom, const unsigned int* surfaces, T* predictions); const T* geom, const unsigned int* surfaces, T* predictions);
// Function, that uses the regression surface to predict y for new data and investigate the model precision // Function, that uses the regression surface to predict y for new data and investigate the model precision
unsigned int s_Score(const unsigned int count_geometries, const unsigned int count_levels, unsigned int s_Score(const unsigned int count_geometries, const unsigned int count_levels,
const T* geom, const unsigned int* surfaces, T* predictions, T* y_values, T* score); const T* geom, const unsigned int* surfaces, T* predictions, const T* y_values, T* score);
// Setting parameters of regression surface models // Setting parameters of regression surface models
unsigned int s_Set_params(const string descriptor, const string kernel, const T* params); unsigned int s_Set_params(const string descriptor, const string kernel, const T* params);
// Getting parameters of regression surface models // Getting parameters of regression surface models