analistica/ex-7/percep.c

#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 function θ(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);
  }

  /* Normalise the weights and force
   * positiveness for easier
   * comparison with other methods.
   */
  double norm = gsl_blas_dnrm2(w);
  if (gsl_vector_isneg(w))
    norm = -norm;
  gsl_vector_scale(w, 1/norm);
  *cut = -bias/norm;

  return w;
}