---

gnn_fir.c

Go to the documentation of this file.

/***************************************************************************
 *  @file  gnn_fir.c
 *  @brief gnn_fir Implementation.
 *
 *  @date   : 27-09-03 01:20
 *  @author : Pedro Ortega C. <peortega@dcc.uchile.cl>
 *  Copyright  2003  Pedro Ortega C.
 ****************************************************************************/
/*
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Library General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */

/**
 * @defgroup gnn_fir_doc gnn_fir : FIR Synaptic Filters
 * @ingroup gnn_atomic_doc
 *
 * The \ref gnn_fir node type implements the function given by
 * \f[ Y = \sum_{k=0}^\tau W(k) x(t-k) \f]
 * where \f$W(k)\f$ with \f$k=0,\ldots,\tau\f$ are matrices of rank
 * \f$m \times n\f$ associated with the delays \f$k\f$.
 * The vectors \f$x(k)\f$ with \f$k=0,\ldots,\tau\f$ area the input
 * vectors to the node. It is important to note that the real input
 * \f$x\f$ is equal to \f$x(t)\f$, and that the vectors
 * \f$x(t-1), x(t-2), \dots, x(t-\tau)\f$ are the input vectors \f$x\f$
 * presented to the node in the previous \f$\tau\f$ iterations, which
 * have been buffered by the node.
 *
 * In other words, each output \f$y_j\f$, \f$j=0,\ldots,m-1\f$ consists of
 *
 * \f[ y_j = \sum_{i=0}^{n-1} \sum_{k=0}^\tau W_{ji}(k) x_i(t-k)
 *         = \sum_{i=0}^{n-1} W_{ji} \ast x_i \f]
 *
 * where the \f$W_{ji}\f$ can be seen as casual discrete filters of order
 * \f$\tau\f$, and so \f$W_{ji} \ast x_i\f$ are the convolutions of the
 * signal \f$x\f$ with the FIR (finite impulse response) filters \f$W_{ji}\f$.
 *
 * The corresponding gradients are:
 *
 * \f[\frac{\partial E}{\partial x_j} =
 *     \sum_j W_{ji}(0) \frac{\partial E}{\partial y_j}\f]
 * or, written in matrix notation:
 * \f[ \frac{\partial E}{\partial x} =
 *     W(k)^t \frac{\partial E}{\partial y} \f]
 *
 * and
 *
 * \f[ \frac{\partial E}{\partial W_{ji}(k)} =
 *     x_i(t-k) \frac{\partial E}{\partial y_j} \f]
 * or:
 * \f[ \frac{\partial E}{\partial W(k)} =
 *     \frac{\partial E}{\partial y}^t x(t-k) \f]
 */

/******************************************/
/* Include Files                          */
/******************************************/

#include <math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_randist.h>
#include "gnn_utilities.h"
#include "gnn_weight.h"



/******************************************/
/* Static Declaration                     */
/******************************************/

/**
 * @brief The structure for a \ref gnn_fir node.
 * @ingroup gnn_fir_doc
 *
 * This datatype holds the information for a \ref gnn_fir node. Basically,
 * it extends the \ref gnn_node with special pointers to get fast accesses
 * to the matrices \f$W(k), k=0,\ldots,\tau\f$ and to the buffered input
 * vector \f$X(t), X(t-1), \ldots, X(t-\tau)\f$.
 *
 * These views exist for the raw parameters, their derivatives and their
 * frozen flags, in order to improve the efficiency.
 */
typedef struct _gnn_fir gnn_fir;

struct _gnn_fir
{
    gnn_node node;   /**< The underlying \ref gnn_node.                   */

    size_t delay;    /**< The amount of delays handled by the node.       */
    size_t offset;   /**< The current offset of \f$x(t)\f$ in the buffer. */

    gsl_matrix_view      *W_view; /**< Views of matrices W(k).            */
    gsl_matrix_view     *dW_view; /**< Views of matrices dW(k).           */
    gsl_matrix_int_view *Wf_view; /**< Views of flags W(k).               */

    gsl_matrix     **W;           /**< Pointers to matrices W(k).         */
    gsl_matrix     **dW;          /**< Pointers to matrices dW(k).        */
    gsl_matrix_int **Wf;          /**< Pointers to flags for the W(k)'s.  */

    gsl_matrix      *xb;          /**< Input buffer matrix.               */
    gsl_vector_view *X_view;      /**< Views of vectors x(t-k).           */
    gsl_vector     **X;           /**< Pointers to vectors x(t-k).        */

};


static int
gnn_fir_f (gnn_node *node,
           const gsl_vector *x,
           const gsl_vector *w,
           gsl_vector *y);


static int
gnn_fir_dx (gnn_node *node,
            const gsl_vector *x,
            const gsl_vector *w,
            const gsl_vector *dy,
            gsl_vector *dx);

static int
gnn_fir_dw (gnn_node *node,
            const gsl_vector *x,
            const gsl_vector *w,
            const gsl_vector *dy,
            gsl_vector *dw);


static void
gnn_fir_destroy (gnn_node *node);




/******************************************/
/* Static Implementation                  */
/******************************************/

/**
 * @brief Computes the evaluation.
 * @ingroup gnn_fir_doc
 *
 * @param node A pointer to a \ref gnn_fir node.
 * @param x    A pointer to the input vector \f$x\f$.
 * @param w    A pointer to the parameter vector \f$w\f$.
 * @param y    A pointer to the output vector \f$y\f$ where the result should
 *             be placed.
 * @return Returns 0 if suceeded.
 */
static int
gnn_fir_f (gnn_node *node,
           const gsl_vector *x,
           const gsl_vector *w,
           gsl_vector *y)
{
    size_t j;
    size_t k;
    gnn_fir *fnode;

    fnode = (gnn_fir *) node;

    /* compute new offset */
    fnode->offset = (fnode->offset > 0)? fnode->offset - 1 : fnode->delay;

    /* copy input into buffer */
    gsl_vector_memcpy (fnode->X[fnode->offset], x);

    /* compute output */
    for (k=0; k<=fnode->delay; ++k)
    {
        j = (fnode->offset + k) % (fnode->delay + 1);
        gsl_blas_dgemv (CblasNoTrans, 1.0, fnode->W[k], fnode->X[j], 0.0, y);
    }

    return 0;
}

/**
 * @brief Computes the \ref gnn_fir gradient with respect to its inputs.
 * @ingroup gnn_fir_doc
 *
 * @param node A pointer to a \ref gnn_fir node.
 * @param x    A pointer to the input vector \f$x\f$.
 * @param w    A pointer to the parameter vector \f$w\f$.
 * @param dy   A pointer to the vector \f$\frac{\partial E}{\partial y}\f$.
 * @param dx   A pointer to the vector \f$\frac{\partial E}{\partial x}\f$
 *             where the result should be placed.
 * @return Returns 0 if suceeded.
 */
static int
gnn_fir_dx (gnn_node *node,
            const gsl_vector *x,
            const gsl_vector *w,
            const gsl_vector *dy,
            gsl_vector *dx)
{
    size_t j;
    size_t k;
    gnn_fir *fnode;

    fnode = (gnn_fir *) node;

    /* compute output */
    gsl_blas_dgemv (CblasTrans, 1.0, fnode->W[0], dy, 0.0, dx);

    return 0;
}


/**
 * @brief Computes the \ref gnn_fir gradient with respect to its parameters.
 * @ingroup gnn_fir_doc
 *
 * @param node A pointer to a \ref gnn_fir node.
 * @param x    A pointer to the input vector \f$x\f$.
 * @param w    A pointer to the parameter vector \f$w\f$.
 * @param dy   A pointer to the vector \f$\frac{\partial E}{\partial y}\f$.
 * @param dw   A pointer to the vector \f$\frac{\partial E}{\partial w}\f$
 *             where the result should be placed.
 * @return Returns 0 if suceeded.
 */
static int
gnn_fir_dw (gnn_node *node,
            const gsl_vector *x,
            const gsl_vector *w,
            const gsl_vector *dy,
            gsl_vector *dw)
{
    size_t j;
    size_t k;
    size_t input_size;
    size_t output_size;
    gsl_matrix_view X;
    gsl_matrix_view DY;
    gnn_fir *fnode;

    fnode = (gnn_fir *) node;

    /* get sizes */
    input_size  = gnn_node_input_get_size (node);
    output_size = gnn_node_output_get_size (node);

    /* get dy matrix view */
    DY = gsl_matrix_view_vector ((gsl_vector *) dy, output_size, 1);
    
    /* compute output */
    for (k=0; k<=fnode->delay; ++k)
    {
        j = (fnode->offset + k) % (fnode->delay + 1);
        
        /* get matrix and vector view */
        X  = gsl_matrix_view_vector ((gsl_vector *) fnode->X[j], 1, input_size);
        
        /* multiply */
        gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                        1.0, &(DY.matrix), &(X.matrix), 1.0, fnode->dW[k]);
    }

    return 0;
}

/**
 * @brief The destroy function for a \ref gnn_fir.
 * @ingroup gnn_fir_doc
 *
 * This functions is the \ref gnn_fir destroy function, which
 * destroys the type's specific data.
 *
 * @param node A pointer to a \ref gnn_fir.
 */
static void
gnn_fir_destroy (gnn_node *node)
{
    gnn_fir *fnode;
    assert (node != NULL);
    
    fnode = (gnn_fir *) node;
    
    if (fnode->W_view  != NULL) free (fnode->W_view);
    if (fnode->dW_view != NULL) free (fnode->dW_view);
    if (fnode->Wf_view != NULL) free (fnode->Wf_view);
    if (fnode->W       != NULL) free (fnode->W);
    if (fnode->dW      != NULL) free (fnode->dW);
    if (fnode->Wf      != NULL) free (fnode->Wf);
    
    if (fnode->X_view  != NULL) free (fnode->X_view);
    if (fnode->X       != NULL) free (fnode->X);
    if (fnode->xb      != NULL) gsl_matrix_free (fnode->xb);

    return;
}



/******************************************/
/* Public Implementation                  */
/******************************************/

/**
 * @brief Creates a new \ref gnn_fir node.
 * @ingroup gnn_fir_doc
 *
 * This function creates a new \ref gnn_node of \ref gnn_fir type.
 *
 * @param input_size The node's input size.
 * @param output_size The node's output size.
 * @return A new node.
 */
gnn_node *
gnn_fir_new (size_t input_size, size_t output_size, size_t delay)
{
    size_t k;
    int status;
    size_t pblock_size;
    size_t param_size;
    gnn_node   *node;
    gnn_fir    *fnode;
    gsl_vector *w;
    gsl_vector *dw;
    gsl_vector_int *f;

    /* check if sizes are positive */
    if (input_size < 1 || output_size < 1)
    {
        GSL_ERROR_VAL ("gnn_fir's sizes should be stricly positive",
                       GSL_EINVAL, NULL);
    }
    if (delay < 1)
    {
        GSL_ERROR_VAL ("gnn_fir's dealy should be stricly positive",
                       GSL_EINVAL, NULL);
    }

    /* compute amount of needed parameters */
    pblock_size = input_size * output_size;
    param_size  = (delay + 1) * pblock_size;

    /* allocate the node */
    node = (gnn_node *) malloc (sizeof (gnn_fir));
    if (node == NULL)
    {
        GSL_ERROR_VAL ("couldn't allocate memory for gnn_fir's node",
                       GSL_ENOMEM, NULL);
    }
    
    /* Initialize the node */
    status = gnn_node_init (node,
                            "gnn_fir",
                            gnn_fir_f,
                            gnn_fir_dx,
                            gnn_fir_dw,
                            gnn_fir_destroy);
    if (status)
    {
        GSL_ERROR_VAL ("could not initialize gnn_fir node", GSL_EINVAL, NULL);
    }
    
    status = gnn_node_set_sizes (node, input_size, output_size, param_size);
    if (status)
    {
        GSL_ERROR_VAL ("could not set sizes for gnn_fir node",
                       GSL_EINVAL, NULL);
    }

    /* get references for building views */
    fnode = (gnn_fir *) node;
    w     = gnn_node_local_get_w (node);
    dw    = gnn_node_local_get_dw (node);
    f     = gnn_node_local_get_f (node);

    /* build input buffer */
    fnode->xb = gsl_matrix_calloc (delay + 1, input_size);
    if (fnode->xb == NULL)
    {
        gnn_node_destroy (node);
        GSL_ERROR_VAL ("could not allocate FIR buffers",
                       GSL_EINVAL, NULL);
    }

    /* build view arrays */
    fnode->W_view
        = (gsl_matrix_view *) malloc (sizeof (gsl_matrix_view) * (delay + 1));
    fnode->dW_view
        = (gsl_matrix_view *) malloc (sizeof (gsl_matrix_view) * (delay + 1));
    fnode->Wf_view
        = (gsl_matrix_int_view *)
          malloc (sizeof (gsl_matrix_int_view) * (delay + 1));
    fnode->W
        = (gsl_matrix **) malloc (sizeof (gsl_matrix *) * (delay + 1));
    fnode->dW
        = (gsl_matrix **) malloc (sizeof (gsl_matrix *) * (delay + 1));
    fnode->Wf
        = (gsl_matrix_int **) malloc (sizeof (gsl_matrix_int *) * (delay + 1));

    fnode->X_view
        = (gsl_vector_view *) malloc (sizeof (gsl_vector_view) * (delay + 1));
    fnode->X
        = (gsl_vector **) malloc (sizeof (gsl_vector *) * (delay + 1));

        
    if (    (fnode->W_view  == NULL)
         || (fnode->dW_view == NULL)
         || (fnode->Wf_view == NULL)
         || (fnode->W       == NULL)
         || (fnode->dW      == NULL)
         || (fnode->Wf      == NULL)
         || (fnode->X_view  == NULL)
         || (fnode->X       == NULL) )
    {
        gnn_node_destroy (node);
        GSL_ERROR_VAL ("could not set sizes for gnn_fir node",
                       GSL_EINVAL, NULL);
    }

    /* build views and pointers */
    for (k=0; k<=delay; ++k)
    {
        gsl_vector_view v;
        gsl_vector_view dv;
        gsl_vector_int_view vf;
        
        v  = gsl_vector_subvector (w, k * pblock_size, pblock_size);
        dv = gsl_vector_subvector (dw, k * pblock_size, pblock_size);
        vf = gsl_vector_int_subvector (f, k * pblock_size, pblock_size);

        fnode->W_view[k]
            = gsl_matrix_view_vector (&(v.vector), output_size, input_size);
        fnode->dW_view[k]
            = gsl_matrix_view_vector (&(dv.vector), output_size, input_size);
        fnode->Wf_view[k]
            = gsl_matrix_int_view_vector (&(vf.vector), output_size, input_size);
        fnode->W[k]  = &(fnode->W_view[k].matrix);
        fnode->dW[k] = &(fnode->dW_view[k].matrix);
        fnode->Wf[k] = &(fnode->Wf_view[k].matrix);

        fnode->X_view[k]
            = gsl_matrix_row (fnode->xb, k);
        fnode->X[k]
            = &(fnode->X_view[k].vector);
    }
    
    /* init members */
    fnode->delay  = delay;
    fnode->offset = delay;
    
    return node;
}

/**
 * @brief Initializes the node's parameters.
 * @ingroup gnn_fir_doc
 *
 * This functions initializes the layer's parameters by values drawn from
 * a uniform distribution in \f$[0; 0.5)\f$.
 *
 * @param layer A pointer to a \ref gnn_fir node.
 * @return Returns 0 on success.
 */
int
gnn_fir_init (gnn_node *node)
{
    int    i;
    int    l;
    gsl_rng    *r;
    gsl_vector *w;

    /* get the random number generator */
    r = gnn_get_rng ();

    /* get the parameters */
    w = gnn_node_local_get_w (node);

    /* initialize weights */
    for (i=0; i<w->size; ++i)
    {
        double rnd;
        rnd = gsl_rng_uniform (r) - 0.5;
        gsl_vector_set (w, i, rnd);
    }

    /* update change */
    gnn_node_local_update (node);

    return 0;
}

---