From 073664e092fe75dd5b7feb63c746b1627d813347 Mon Sep 17 00:00:00 2001 From: TheOnePath Date: Wed, 11 Sep 2024 13:10:41 +0100 Subject: [PATCH] Added src files --- src/demo.c | 101 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/gp.c | 62 +++++++++++++++++++++++++++++++ src/kernels.c | 55 +++++++++++++++++++++++++++ src/solve.c | 47 +++++++++++++++++++++++ 4 files changed, 265 insertions(+) create mode 100644 src/demo.c create mode 100644 src/gp.c create mode 100644 src/kernels.c create mode 100644 src/solve.c diff --git a/src/demo.c b/src/demo.c new file mode 100644 index 0000000..54fc9e1 --- /dev/null +++ b/src/demo.c @@ -0,0 +1,101 @@ +#include + +#include +#include + +#include "kernels.h" +#include "solve.h" +#include "gp.h" + +#define RNG_MAX 10 +#define X_MAX 10 + +int linspace ( double start, double end, int length, gsl_matrix *out ) { + if ( out->size2 != length ) + return 1; + + double by = end / length; + for ( int i = 0 ; i < length ; i++ ) + (void)gsl_matrix_set(out, 0, i, i*by); + + return 0; +} + +void print_matrix ( gsl_matrix *matrix ) { + printf("\n"); + for ( int i = 0 ; i < matrix->size1 ; i++ ) { + for ( int j = 0 ; j < matrix->size2 ; j++ ) { + printf(" %.6e ", gsl_matrix_get(matrix, i, j)); + } + printf("\n"); + } + printf("\n"); +} + +void print_flatten ( gsl_matrix *matrix ) { + printf("\n"); + if ( matrix->size1 == 1 ) { + for ( int i = 0 ; i < matrix->size2 ; i++ ) { + if ( i == matrix->size2-1 ) { + printf("%g\n", gsl_matrix_get(matrix, 0, i)); + return; + } + printf("%g, ", gsl_matrix_get(matrix, 0, i)); + } + printf("\n"); + return; + } + + for ( int i = 0 ; i < matrix->size1 ; i++ ) { + if ( i == matrix->size1-1 ) { + printf("%g\n", gsl_matrix_get(matrix, i, 0)); + return; + } + printf("%g, ", gsl_matrix_get(matrix, i, 0)); + } + printf("\n"); +} + + +int main(void) +{ + /* RNG is always seeded with 0 */ + gsl_rng *rng = gsl_rng_alloc(gsl_rng_taus); + + double sigmaf = 5.0160194, ell = 0.8328418, noise = 0.05; + + /* design matrix and training set */ + gsl_matrix *X = gsl_matrix_alloc(1, X_MAX); + gsl_matrix *Y = gsl_matrix_alloc(1, X_MAX); + gsl_matrix *Xs = gsl_matrix_alloc(1, 100); + + /* Initialise our X and Y with some random values */ + for ( int i = 0 ; i < X->size2 ; i++ ) { + (void)gsl_matrix_set(X, 0, i, gsl_rng_uniform_int(rng, RNG_MAX)); + (void)gsl_matrix_set(Y, 0, i, gsl_rng_uniform_int(rng, RNG_MAX)); + } + (void)linspace(0, gsl_matrix_max(X)+1, 100, Xs); + + kernel_info info = { + .kern = &squared_exp, + .length_scale = ell, + .sigmaf = sigmaf, + }; + gp_model model = gpfit(X, Y, Xs, info, noise, &inv_chol); + + print_flatten(X); + print_flatten(Y); + print_flatten(Xs); + print_flatten(model.meanf); + + /* free all allocations made */ + gsl_rng_free(rng); + gsl_matrix_free(X); + gsl_matrix_free(Y); + gsl_matrix_free(Xs); + + gsl_matrix_free(model.meanf); + gsl_matrix_free(model.covf); + + return EXIT_SUCCESS; +} diff --git a/src/gp.c b/src/gp.c new file mode 100644 index 0000000..d374a5d --- /dev/null +++ b/src/gp.c @@ -0,0 +1,62 @@ +#include +#include + +#include "solve.h" +#include "kernels.h" +#include "gp.h" + +gp_model +gpfit +( gsl_matrix *x, gsl_matrix *y, gsl_matrix *xs, kernel_info info, double noise, + invert_t inv_t ) +{ + /* kernel allocations for: K(X,X), K(XS,XS), K(X,XS) and K(XS,X) */ + gsl_matrix *Kxx = gsl_matrix_alloc(x->size2, x->size2); + gsl_matrix *Kxs = gsl_matrix_alloc(x->size2, xs->size2); + gsl_matrix *Ksx = gsl_matrix_alloc(xs->size2, x->size2); + /* the noise to add with the training data */ + gsl_matrix *sigma_noise = gsl_matrix_calloc(Kxx->size1, Kxx->size2); + + gsl_matrix *inv = gsl_matrix_alloc(Kxx->size1, Kxx->size2); + gsl_matrix *tmp = gsl_matrix_alloc(Kxs->size2, inv->size2); + + /* create new GP struct for model */ + gp_model model = { + .meanf = gsl_matrix_alloc(xs->size2, xs->size1), + .covf = gsl_matrix_alloc(xs->size2, xs->size2), + }; + + for ( int i = 0 ; i < sigma_noise->size1 ; i++ ) + gsl_matrix_set(sigma_noise, i, i, noise*noise); + + /* compute the covariance kernels */ + (void)info.kern(info.sigmaf, info.length_scale, x, x, Kxx); + (void)info.kern(info.sigmaf, info.length_scale, xs, xs, model.covf); + (void)info.kern(info.sigmaf, info.length_scale, x, xs, Kxs); + gsl_matrix_transpose_memcpy(Ksx, Kxs); + + /* add the noise to the K(X,X) kernel */ + gsl_matrix_add(Kxx, sigma_noise); + + /* invert the K(X,X) matrix with selected method */ + (void)inv_t(Kxx, inv); + + /* Matrix multiplication of K(XS,X)*[K(X,X)+s^2I]^-1*y */ + (void)gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Ksx, inv, 0.0, tmp); + (void)gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, tmp, y, 0.0, model.meanf); + + /* calculate the cov(f) */ + gsl_matrix *tmp2 = gsl_matrix_alloc(model.covf->size1, model.covf->size2); + (void)gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, tmp, Kxs, 0.0, tmp2); + (void)gsl_matrix_sub(model.covf, tmp2); + + gsl_matrix_free(Kxx); + gsl_matrix_free(Ksx); + gsl_matrix_free(Kxs); + gsl_matrix_free(tmp); + gsl_matrix_free(tmp2); + gsl_matrix_free(inv); + gsl_matrix_free(sigma_noise); + + return model; +} diff --git a/src/kernels.c b/src/kernels.c new file mode 100644 index 0000000..9044fae --- /dev/null +++ b/src/kernels.c @@ -0,0 +1,55 @@ +#include +#include + +#include "kernels.h" + +/* estimation of sqrt(3) */ +#define SQRT3 1.73205080757 + +int squared_exp ( double sigmaf, double ell, gsl_matrix *A, gsl_matrix *B, + gsl_matrix *result ) +{ + int i, j, n = A->size2, m = B->size2; + /* ensure that the matrix result is the same size as the matrixs */ + if ( result->size1 != n || result->size2 != m ) + return 1; + + gsl_matrix_set_zero(result); + double r, rval, x, y; + for ( i = 0 ; i < n ; i++ ) { + for ( j = 0 ; j < m ; j++ ) { + x = gsl_matrix_get(A, 0, i); + y = gsl_matrix_get(B, 0, j); + + r = sqrt((x - y)*(x - y)); + rval = sigmaf*sigmaf*exp(-(r*r)/(2*ell*ell)); + gsl_matrix_set(result, i, j, rval); + } + } + + return 0; +} + +int matern32 ( double sigmaf, double ell, gsl_matrix *A, gsl_matrix *B, + gsl_matrix *result ) +{ + int i, j, n = A->size2, m = B->size2; + /* ensure that the matrix result is the same size as the matrixs */ + if ( result->size1 != n || result->size2 != m ) + return 1; + + double r, rval, x, y; + for ( i = 0 ; i < n ; i++ ) { + for ( j = 0 ; j < m ; j++ ) { + x = gsl_matrix_get(A, 0, i); + y = gsl_matrix_get(B, 0, j); + + r = sqrt((x - y)*(x - y)); + /* sigmaf^2*(1+sqrt(3)*r/ell)*exp(-sqrt(3)*r/ell) */ + rval = sigmaf*sigmaf*(1+SQRT3*r/ell)*exp(-SQRT3*r/ell); + gsl_matrix_set(result, i, j, rval); + } + } + + return 0; +} diff --git a/src/solve.c b/src/solve.c new file mode 100644 index 0000000..ab66048 --- /dev/null +++ b/src/solve.c @@ -0,0 +1,47 @@ +#include +#include +#include +#include +#include + +/* + * https://gist.github.com/yaopu/b6e0151e790ecbdcaa269e633685eae4 + * */ +int inv_lu ( gsl_matrix *A, gsl_matrix *inverse ) { + /* square matrix required */ + if ( A->size1 != A->size2 ) + return 1; + + gsl_permutation *perms = gsl_permutation_alloc(A->size1); + int s; + + // Compute the LU decomposition of this matrix + gsl_linalg_LU_decomp(A, perms, &s); + + // Compute the inverse of the LU decomposition + gsl_linalg_LU_invert(A, perms, inverse); + gsl_permutation_free(perms); + + return 0; +} + +int inv_chol ( gsl_matrix *A, gsl_matrix *inverse ) { + if ( A->size1 != A->size2 || + A->size1 * A->size2 != inverse->size1 * inverse->size2 ) + return 1; + + gsl_permutation *p = gsl_permutation_alloc(A->size1); + gsl_matrix *ldlt = gsl_matrix_alloc(A->size1, A->size2); + + /* perform the pivoted cholesky decomposition to find LDLT */ + gsl_matrix_memcpy(ldlt, A); + gsl_linalg_pcholesky_decomp(ldlt, p); + + /* invert A with LDLT and permutation */ + gsl_linalg_pcholesky_invert(ldlt, p, inverse); + + gsl_permutation_free(p); + gsl_matrix_free(ldlt); + return 0; +} +