From e4791356c7c6ac63aeddc430c99069696d916ecc Mon Sep 17 00:00:00 2001
From: exapde <exapde@gmail.com>
Date: Fri, 2 Dec 2022 10:29:46 -0500
Subject: [PATCH] Add regularization parameter to make the fitting more robust.

---
 examples/PACKAGES/pod/Ta/Ta_data.pod           |  1 +
 examples/PACKAGES/pod/Ta_Quadratic/Ta_data.pod |  1 +
 src/ML-POD/fitpod_command.cpp                  | 12 ++++++++----
 src/ML-POD/fitpod_command.h                    |  2 +-
 4 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/examples/PACKAGES/pod/Ta/Ta_data.pod b/examples/PACKAGES/pod/Ta/Ta_data.pod
index 6aa256a61d..b345496f71 100644
--- a/examples/PACKAGES/pod/Ta/Ta_data.pod
+++ b/examples/PACKAGES/pod/Ta/Ta_data.pod
@@ -6,6 +6,7 @@ path_to_test_data_set "XYZ"
 
 fitting_weight_energy 100.0
 fitting_weight_force 1.0
+fitting_regularization_parameter 1e-10
 
 error_analysis_for_training_data_set 1
 error_analysis_for_test_data_set 0
diff --git a/examples/PACKAGES/pod/Ta_Quadratic/Ta_data.pod b/examples/PACKAGES/pod/Ta_Quadratic/Ta_data.pod
index 5012baab22..c7436ea26e 100644
--- a/examples/PACKAGES/pod/Ta_Quadratic/Ta_data.pod
+++ b/examples/PACKAGES/pod/Ta_Quadratic/Ta_data.pod
@@ -6,6 +6,7 @@ path_to_test_data_set "../Ta/XYZ"
 
 fitting_weight_energy 100.0
 fitting_weight_force 1.0
+fitting_regularization_parameter 1e-10
 
 error_analysis_for_training_data_set 1
 error_analysis_for_test_data_set 0
diff --git a/src/ML-POD/fitpod_command.cpp b/src/ML-POD/fitpod_command.cpp
index 9b66339c48..7d11da503c 100644
--- a/src/ML-POD/fitpod_command.cpp
+++ b/src/ML-POD/fitpod_command.cpp
@@ -188,7 +188,7 @@ void FitPOD::read_data_file(double *fitting_weights, std::string &file_format,
 
     if (keywd == "fitting_weight_energy") fitting_weights[0] = utils::numeric(FLERR,words[1],false,lmp);
     if (keywd == "fitting_weight_force") fitting_weights[1] = utils::numeric(FLERR,words[1],false,lmp);
-    if (keywd == "fitting_weight_stress") fitting_weights[2] = utils::numeric(FLERR,words[1],false,lmp);
+    if (keywd == "fitting_weight_stress") fitting_weights[2] = utils::numeric(FLERR,words[1],false,lmp);   
     if (keywd == "error_analysis_for_training_data_set") fitting_weights[3] = utils::numeric(FLERR,words[1],false,lmp);
     if (keywd == "error_analysis_for_test_data_set") fitting_weights[4] = utils::numeric(FLERR,words[1],false,lmp);
     if (keywd == "energy_force_calculation_for_training_data_set") fitting_weights[5] = utils::numeric(FLERR,words[1],false,lmp);
@@ -197,7 +197,8 @@ void FitPOD::read_data_file(double *fitting_weights, std::string &file_format,
     if (keywd == "fraction_test_data_set") fitting_weights[8] = utils::numeric(FLERR,words[1],false,lmp);
     if (keywd == "randomize_training_data_set") fitting_weights[9] = utils::numeric(FLERR,words[1],false,lmp);
     if (keywd == "randomize_test_data_set") fitting_weights[10] = utils::numeric(FLERR,words[1],false,lmp);
-
+    if (keywd == "fitting_regularization_parameter") fitting_weights[11] = utils::numeric(FLERR,words[1],false,lmp);
+    
     // other settings
 
     if (keywd == "file_format") file_format = words[1];
@@ -223,6 +224,7 @@ void FitPOD::read_data_file(double *fitting_weights, std::string &file_format,
     utils::logmesg(lmp, "fitting weight for energy: {}\n", fitting_weights[0]);
     utils::logmesg(lmp, "fitting weight for force: {}\n", fitting_weights[1]);
     utils::logmesg(lmp, "fitting weight for stress: {}\n", fitting_weights[2]);
+    utils::logmesg(lmp, "fitting regularization parameter: {}\n", fitting_weights[11]);
     utils::logmesg(lmp, "**************** End of Data File ****************\n");
   }
 }
@@ -1269,10 +1271,12 @@ void FitPOD::least_squares_fit(datastruct data)
     for (int i = 0; i<nd*nd; i++)
       desc.A[i] = desc.A[i]*maxb;
 
+    double regularizing_parameter = data.fitting_weights[11];
+    
     for (int i = 0; i<nd; i++) {
       desc.c[i] = desc.b[i];
-      desc.A[i + nd*i] = desc.A[i + nd*i]*(1.0 + SMALL);
-      if (desc.A[i + nd*i] < SMALL) desc.A[i + nd*i] = SMALL;
+      desc.A[i + nd*i] = desc.A[i + nd*i]*(1.0 + regularizing_parameter);
+      if (desc.A[i + nd*i] < regularizing_parameter) desc.A[i + nd*i] = regularizing_parameter;
     }
 
     // solving the linear system A * c = b
diff --git a/src/ML-POD/fitpod_command.h b/src/ML-POD/fitpod_command.h
index 7401515684..5b836a44fe 100644
--- a/src/ML-POD/fitpod_command.h
+++ b/src/ML-POD/fitpod_command.h
@@ -62,7 +62,7 @@ public:
     int randomize = 1;
     double fraction = 1.0;
 
-    double fitting_weights[12] = {0.0, 0.0, 0.0, 1, 1, 0, 0, 1, 1, 1, 1, 0};
+    double fitting_weights[12] = {100.0, 1.0, 0.0, 1, 1, 0, 0, 1, 1, 1, 1, 1e-10};
 
     void copydatainfo(datastruct &data) {
       data.data_path = data_path;