demo gaussian process models

src/constitutive/GaussianProcessModel.cpp

@@ -0,0 +1,114 @@
+#include "GaussianProcessModel.h"
+
+GaussianProcessModel::GaussianProcessModel(int n_train, int n_param, TacsScalar Xtrain[], TacsScalar alpha[]) {
+    // constructor for the base GaussianProcessModel class
+    this->n_train = n_train;
+    this->n_param = n_param;
+
+    // create a new copy of the Xtrain, alpha arrays on this process
+    this->Xtrain = new TacsScalar[n_train*n_param];
+    for (int ii = 0; ii < n_train*n_param; ii++) {
+        this->Xtrain[ii] = Xtrain[ii];
+    }
+
+    this->alpha = new TacsScalar[n_train];
+    for (int jj = 0; jj < n_train; jj++) {
+        this->alpha[jj] = alpha[jj];
+    }
+}
+
+GaussianProcessModel::~GaussianProcessModel() {
+    // destructor for the base GaussianProcessModel class
+    delete[] this->Xtrain;
+    this->Xtrain = nullptr;
+
+    delete[] this->alpha;
+    this->alpha = nullptr;
+}
+
+TacsScalar GaussianProcessModel::predictMeanTestData(const TacsScalar* Xtest) {
+    // Xtest is an array of size n_param (for one test data point)
+    // use the equation mean(Ytest) = cov(Xtest,X_train) @ alpha [this is a dot product]
+    // where Ytest is a scalar
+    TacsScalar Ytest = 0.0;
+
+    // iterate over each training data point, get the cross-term covariance and add the coefficient alpha for it
+    for (int itrain = 0; itrain < n_train; itrain++) {
+        const TacsScalar* loc_Xtrain = &Xtrain[n_param*itrain];
+        Ytest += kernel(Xtest, loc_Xtrain) * alpha[itrain];
+    }
+    return Ytest;
+}
+
+TacsScalar GaussianProcessModel::predictMeanTestDataSens(const TacsScalar* Xtest, TacsScalar* Xtestsens) {
+    // Xtest is an array of size n_param (for one test data point)
+    // use the equation mean(Ytest) = cov(Xtest,X_train) @ alpha [this is a dot product]
+    // where Ytest is a scalar
+    TacsScalar Ytest = 0.0;
+    TacsScalar* Xtestsens = new TacsScalar[n_param];
+    memset(Xtestsens, 0, n_param*sizeof(TacsScalar));
+
+    // iterate over each training data point, get the cross-term covariance and add the coefficient alpha for it
+    for (int itrain = 0; itrain < n_train; itrain++) {
+        const TacsScalar* loc_Xtrain = &Xtrain[n_param*itrain];
+        Ytest += kernel(Xtest, loc_Xtrain) * alpha[itrain];
+        // add the kernel sensitivity for this training point
+        kernelSens(Xtest, loc_Xtrain, alpha[itrain], Xtestsens);
+    }
+    return Ytest;
+}
+
+AxialGaussianProcessModel::AxialGaussianProcessModel() {
+    // overwrite the debug constructor
+    // use only 1 training point
+    TacsScalar* Xtrain = new TacsScalar[N_PARAM*1];
+    memset(Xtrain, 0, N_PARAM*sizeof(TacsScalar));
+    TacsScalar* alpha = new TacsScalar[1];
+    alpha[0] = 1.0;
+
+    // call the main constructor
+    AxialGaussianProcessModel(1,Xtrain,alpha);
+}
+
+void AxialGaussianProcessModel::setDefaultHyperParameters() {
+    S1 = 1e-1;
+    S2 = 3e-1;
+    c = -1;
+    L1 = 0.2;
+    S4 = 1.0;
+    S5 = 1.0;
+    L2 = 0.3;
+    alpha1 = 2.0;
+    L3 = 4.0;
+    S6 = 1.0;
+}
+
+TacsScalar AxialGaussianProcessModel::kernel(const TacsScalar* Xtest, const TacsScalar* Xtrain) {
+    // define the kernel function k(*,*) on training and testing points for one training point
+    // the entries are [log(rho_0), log(xi), log(gamma), log(zeta)]
+
+    // log(rho_0) direction 0
+    TacsScalar kernel0 = S1*S1 + S2*S2*soft_relu(Xtest[0]-c, alpha1) * soft_relu(Xtrain[0] - c, alpha1);
+
+    // log(xi) direction 1
+    TacsScalar d1 = Xtest[1] - Xtrain[1];
+    TacsScalar fact1 = soft_relu(1 - soft_relu(Xtest[1], 1.0), 10);
+    TacsScalar fact2 = soft_relu(1 - soft_relu(Xtrain[1], 1.0), 10);
+    TacsScalar rho0term3 = soft_relu(-Xtest[1], 10) * soft_relu(-Xtrain[1], 10);
+    TacsScalar kernel1 = exp(-0.5 * (d1*d1/L1*L1)) * fact1 * fact2 + S4 + S5 * rho0term3;
+
+    // log(gamma) direction 2
+    TacsScalar d2 = Xtest[2] - Xtrain[2];
+    TacsScalar kernel2 = exp(-0.5 * d2 * d2 / L2 / L2);
+
+    // log(zeta) direction 3
+    TacsScalar d3 = Xtest[3] - Xtrain[3];
+    TacsScalar kernel3 = S6 * exp(-0.5 * d3 * d3 / L3 / L3);
+
+    return kernel0 * kernel1 * kernel2 + 2.0 * kernel3
+}
+
+TacsScalar kernelSens(const TacsScalar* Xtest, const TacsScalar* Xtrain, TacsScalar alpha, TacsScalar* Xtestsens) {
+    // add into the Xtestsens (don't reset to zero) for x_test = log[nondim params] vector
+
+}

src/constitutive/GaussianProcessModel.h

@@ -0,0 +1,122 @@
+/*
+========================================================================
+Gaussian Process Model for TACS Buckling Constraints
+========================================================================
+@File   :   GaussianProcessModel.h
+@Date   :   2024/05/10
+@Author :   Sean Phillip Engelstad
+@Description : Use Gaussian Processes for machine learning techniques to interpolate
+learned buckling constraints on a training dataset to new test data points. This
+approach is implemented into the TACSGPBladeStiffenedShellConstitutive class in TACS
+for more physically accurate buckling constraints of stiffened panels.
+*/
+
+#pragma once
+
+// =============================================================================
+// Extension Includes
+// =============================================================================
+#include "TacsUtilities.h"
+
+// =============================================================================
+// Class Declaration
+// =============================================================================
+
+class GaussianProcessModel {
+public:
+  GaussianProcessModel(int n_train, int n_param, TacsScalar Xtrain[], TacsScalar alpha[]);
+  ~GaussianProcessModel();
+
+  // predict the test data at a single point using a matrix-vector product
+  // this is the mean test data prediction. The offline training beforehand trains the mean
+  // surface of the training set to match mu - 3 * sigma, a bounding curve for buckling 
+  // so that way we don't need to invert a covariance matrix each time.
+  // here Xtest is one point prediction with an array of length n_param
+  TacsScalar predictMeanTestData(const TacsScalar* Xtest);
+  TacsScalar predictMeanTestDataSens(const TacsScalar* Xtest, TacsScalar* Xtestsens);
+
+  static inline TacsScalar soft_relu(TacsScalar x, TacsScalar rho) {
+    return 1.0 / rho * log(1 + exp(rho * x));
+  };
+
+protected:
+    // virtual functions for the kernel definition and its sensitivity
+    virtual TacsScalar kernel(const TacsScalar* Xtest, const TacsScalar* Xtrain);
+    virtual void kernelSens(const TacsScalar* Xtest, const TacsScalar* Xtrain, TacsScalar alpha, TacsScalar* Xtestsens);
+
+    int n_train;
+    int n_param;
+    // rank 1-tensor of length [n_param*n_train] [X]
+    // if each point of Xtrain has data [rho0, xi, gamma, delta, zeta] with n_Train=5 then
+    // the entries are basically [rho01, xi1, gamma1, delta1, zeta1, rho02, xi2, gamma2, delta2, zeta2, ..., zetaN]
+    TacsScalar *Xtrain; 
+    TacsScalar *alpha;
+};
+
+class AxialGaussianProcessModel : GaussianProcessModel {
+public:  
+  AxialGaussianProcessModel(int n_train, const TacsScalar Xtrain[], const TacsScalar alpha[]) : GaussianProcessModel(n_train, N_PARAM, Xtrain, alpha) {};
+  AxialGaussianProcessModel(); // empty constructor for debugging purposes
+  ~AxialGaussianProcessModel() : ~GaussianProcessModel() {};
+  void setDefaultHyperParameters();
+protected:
+    // here Xtest, Xtrain are each length 5 arrays (N_PARAM) [just uses one train and one test point here]
+    // these are overwritten from subclass
+    TacsScalar kernel(const TacsScalar* Xtest, const TacsScalar* Xtrain) {};
+    void kernelSens(const TacsScalar* Xtest, const TacsScalar* Xtrain, TacsScalar alpha, TacsScalar* Xtestsens);
+
+    // set the default hyperparameters of the model
+    const TacsScalar S1, S2, c, L1, S4, S5, L2, alpha1, L3, S6;
+
+    // there are 5 parameters [rho_0, xi, gamma, zeta] for the axial model
+    static int N_PARAM = 4;
+};
+
+class ShearGaussianProcessModel : AxialGaussianProcessModel {
+public:  
+  ShearGaussianProcessModel(int n_train, const TacsScalar Xtrain[], const TacsScalar alpha[]) : GaussianProcessModel(n_train, N_PARAM, Xtrain, alpha);
+  ShearGaussianProcessModel(); // empty constructor for debugging purposes
+  ~ShearGaussianProcessModel();
+  void setdefaultHyperParameters();
+protected:
+    // here Xtest, Xtrain are each length 5 arrays (N_PARAM) [just uses one train and one test point here]
+    TacsScalar kernel(const TacsScalar* Xtest, const TacsScalar* Xtrain) {};
+    void kernelSens(const TacsScalar* Xtest, const TacsScalar* Xtrain, TacsScalar alpha, TacsScalar* Xtestsens);
+
+    // set the default hyperparameters of the model
+    const TacsScalar S1, S2, c, L1, S4, S5, L2, alpha1, L3, S6;
+
+    // there are 5 parameters [rho_0, xi, gamma, zeta] for the axial model
+    static int N_PARAM = 4;
+};
+
+class CripplingGaussianProcessModel : GaussianProcessModel {
+public:  
+  CripplingGaussianProcessModel(int n_train, const TacsScalar Xtrain[], const TacsScalar alpha[]) : GaussianProcessModel(n_train, N_PARAM, Xtrain, alpha);
+  CripplingGaussianProcessModel(); // empty constructor for debugging purposes
+  ~CripplingGaussianProcessModel();
+
+  void setdefaultHyperParameters() {
+    S1 = 1e-1;
+    S2 = 3e-1;
+    c = -1;
+    L1 = 0.2;
+    S4 = 1.0;
+    S5 = 1.0;
+    L2 = 0.3;
+    alpha1 = 2.0;
+    L3 = 4.0;
+    S6 = 1.0;
+  };
+
+protected:
+    // here Xtest, Xtrain are each length 5 arrays (N_PARAM) [just uses one train and one test point here]
+    TacsScalar kernel(const TacsScalar* Xtest, const TacsScalar* Xtrain) {};
+    TacsScalar kernelSens(const TacsScalar* Xtest, const TacsScalar* Xtrain, TacsScalar* Xtestsens);
+
+    // set the default hyperparameters of the model
+    const TacsScalar S1, S2, c, L1, S4, S5, L2, alpha1, L3, S6;
+
+    // there are 5 parameters [rho_0, xi, gen_eps, zeta] for the axial model
+    static int N_PARAM = 3;
+};