ex-7: implement perceptron

ex-7/main.c

@@ -1,12 +1,14 @@
 #include <stdio.h>
 #include <string.h>
 #include "fisher.h"
+#include "percep.h"
 
 /* Options for the program */
 struct options {
   char *mode;
   size_t nsig;
   size_t nnoise;
+  int iter;
 };
 
 
@@ -16,18 +18,20 @@ int main(int argc, char **argv) {
   opts.mode   = "fisher";
   opts.nsig   =  800;
   opts.nnoise = 1000;
+  opts.iter   = 5;
 
   /* Process CLI arguments */
   for (size_t i = 1; i < argc; i++) {
          if (!strcmp(argv[i], "-m")) opts.mode = argv[++i];
     else if (!strcmp(argv[i], "-s")) opts.nsig = atol(argv[++i]);
     else if (!strcmp(argv[i], "-n")) opts.nnoise = atol(argv[++i]);
+    else if (!strcmp(argv[i], "-i")) opts.nnoise = atoi(argv[++i]);
     else {
       fprintf(stderr, "Usage: %s -[hiIntp]\n", argv[0]);
       fprintf(stderr, "\t-h\tShow this message.\n");
       fprintf(stderr, "\t-m MODE\tThe disciminant to use: 'fisher' for "
                       "Fisher linear discriminant, 'percep' for perceptron.\n");
-      fprintf(stderr, "\t-s N\tThe number of events in signal class.\n");
+      fprintf(stderr, "\t-i N\tThe number of training iterations (for perceptron).\n");
       fprintf(stderr, "\t-n N\tThe number of events in noise class.\n");
       return EXIT_FAILURE;
     }
@@ -47,6 +51,9 @@ int main(int argc, char **argv) {
   sample_t *signal = generate_normal(r, opts.nsig, &par_sig);
   sample_t *noise  = generate_normal(r, opts.nnoise, &par_noise);
 
+  gsl_vector *w;
+  double t_cut;
+
   if (!strcmp(opts.mode, "fisher")) {
     /* Fisher linear discriminant
      *
@@ -55,19 +62,36 @@ int main(int argc, char **argv) {
      * cut which determines the class for each
      * projected point.
      */
-    double ratio = opts.nsig / (double)opts.nnoise;
-    gsl_vector *w = fisher_proj(signal, noise);
-    double t_cut  = fisher_cut(ratio, w, signal, noise);
-
     fputs("# Linear Fisher discriminant\n\n", stderr);
-    fprintf(stderr, "* w: [%.3f, %.3f]\n",
-        gsl_vector_get(w, 0),
-        gsl_vector_get(w, 1));
-    fprintf(stderr, "* t_cut: %.3f\n", t_cut);
-
-    gsl_vector_fprintf(stdout, w, "%g");
-    printf("%f\n", t_cut);
+    double ratio = opts.nsig / (double)opts.nnoise;
+    w = fisher_proj(signal, noise);
+    t_cut = fisher_cut(ratio, w, signal, noise);
   }
+  else if (!strcmp(opts.mode, "percep")) {
+    /* Perceptron
+     *
+     * Train a single perceptron on the
+     * dataset to get an approximate
+     * solution in `iter` iterations.
+     */
+    fputs("# Perceptron \n\n", stderr);
+    w = percep_train(signal, noise, opts.iter, &t_cut);
+  }
+  else {
+    fputs("\n\nerror: invalid mode. select either"
+          " 'fisher' or 'percep'\n", stderr);
+    return EXIT_FAILURE;
+  }
+
+  /* Print the results of the method
+   * selected: weights and threshold.
+   */
+  fprintf(stderr, "* w: [%.3f, %.3f]\n",
+      gsl_vector_get(w, 0),
+      gsl_vector_get(w, 1));
+  fprintf(stderr, "* t_cut: %.3f\n", t_cut);
+  gsl_vector_fprintf(stdout, w, "%g");
+  printf("%f\n", t_cut);
 
   /* Print data to stdout for plotting.
    * Note: we print the sizes to be able

ex-7/percep.c

@@ -0,0 +1,86 @@
+#include "common.h"
+#include "percep.h"
+#include <math.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+
+/* `iterate(data, w, b, d, r)` performs one
+ * iteration of the perceptron training.
+ * 
+ * For each point xi compute:
+ *
+ *   1. yi = θ((w, xi) + b)
+ *
+ *   2. Δi = r⋅(d - yi)
+ *
+ *   3. w = w + Δi⋅xi, b = b + Δi
+ *
+ * The results are computed in-place.
+ */
+void iterate(
+  gsl_matrix *data,
+  gsl_vector *weight,
+  double *bias,
+  int expected,
+  double rate) {
+
+  double proj, delta;
+  gsl_vector *x = gsl_vector_alloc(2);
+
+  for (size_t i = 0; i < data->size1; i++) {
+    /* Get a vector view of the 
+     * current row in the data matrix.
+     */
+    gsl_vector row = gsl_matrix_const_row(data, i).vector;
+    gsl_vector_memcpy(x, &row);
+
+    /* Project x onto the weight vector. */
+    gsl_blas_ddot(weight, x, &proj);
+
+    /* Calculate Δ
+     * Note: the step functions θ(x) is computed
+     * by negating the sign bit of the floating
+     * point.
+     */
+    delta = rate * (expected - !signbit(proj + *bias));
+
+    /* Update weight and bias. */
+    *bias += delta;
+    gsl_vector_scale(x, delta);
+    gsl_vector_add(weight, x);
+  }
+
+  gsl_vector_free(x);
+}
+
+/* `percep_train(sig, noise, iter, &cut)`
+ * train a single perceptron to discriminate `sig`
+ * from `noise`.
+ *
+ * The weigths are adjusted `iter` times and
+ * returned, the bias/threshold is stored in the
+ * `cut` argument.
+ */
+gsl_vector *percep_train(
+  sample_t *signal, sample_t *noise,
+  int iter, double *cut) {
+
+  /* Initially set weights/bias to zero
+   * and a high learning rate.
+   */
+  gsl_vector *w = gsl_vector_calloc(2);
+  double bias = 0;
+  double rate = 0.8;
+
+  /* Go trough the sample `iter` times
+   * and recalculate the weights.
+   */
+  for (int i = 0; i < iter; i++) {
+    iterate( noise->data, w, &bias, 0, rate);
+    iterate(signal->data, w, &bias, 1, rate);
+  }
+
+  *cut = -bias;
+
+  return w;
+}

ex-7/percep.h

@@ -0,0 +1,14 @@
+#include "common.h"
+#include <gsl/gsl_vector.h>
+
+/* `percep_train(sig, noise, iter, &cut)`
+ * train a single perceptron to discriminate `sig`
+ * from `noise`.
+ *
+ * The weigths are adjusted `iter` times and
+ * returned, the bias/threshold is stored in the
+ * `cut` argument.
+ */
+gsl_vector *percep_train(
+  sample_t *signal, sample_t *noise,
+  int iter, double *cut);