Added src files

src/demo.c

@@ -0,0 +1,101 @@
+#include <stdio.h>
+
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_rng.h>
+
+#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;
+}

src/gp.c

@@ -0,0 +1,62 @@
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+
+#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;
+}

src/kernels.c

@@ -0,0 +1,55 @@
+#include <math.h>
+#include <gsl/gsl_matrix.h>
+
+#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;
+}

src/solve.c

@@ -0,0 +1,47 @@
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_matrix_double.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_linalg.h>
+
+/*
+ * 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;
+}
+