root/exactrbn/trunk/src/rbn_solve.h

/* $Id$ $URL$ */
#ifndef RBN_SOLVE_DOT_AITCH
#define RBN_SOLVE_DOT_AITCH

/**\file rbn_solve.h
 * Template classes MyVec, Example and ExactRBFNet for implementing
 * and solving exact-interpolation radial basis function neural networks.
 * Copyright Trent Apted 2008 All rights reserved.
 *
 * \author Trent Apted <tapted@it.usyd.edu.au>
 * \version $Revision$
 * \date $Date$
 *
 * Macros to influence behaviour (if defined before # include):
 * - # define RBFNET_READ_ONLY // make a "read only" network that must be read from disk
 *                              (removes dependencies on lapack)
 * - # define RBFNET_DEBUG     // prints debugging information to stderr
 * - # define RBFNET_VDEBUG    // prints more debugging information to stderr
 * \n\n
 */

#include <vector>
namespace stl = std;

#include <unistd.h>
#include <valarray>
#include <istream>
#include <ostream>
#include <fstream>
#include <typeinfo>

/**
 * Simple N-dimensional vector of type T.
 * \param T the type of elements in the vector
 * \param DIM the dimension (size) of the vector
 */
template <class T, int DIM>
    struct MyVec {
        T v[DIM]; ///< The components
        /** Indexation */
        T& operator[](size_t i) {return v[i];}
        /** Const indexation */
        const T& operator[](size_t i) const {return v[i];}
    };

/**
 * A "training" example -- an instance of a function mapping a value
 * from the domain of T^INDIM to the range T^OUTDIM.
 * \param T the type of element (e.g. float, double, complex)
 * \param INDIM the dimension of the domain
 * \param OUTDIM the dimension of the range
 */
template <class T, int INDIM, int OUTDIM>
struct Example {
    MyVec<T, INDIM> in;   ///< The input vector
    MyVec<T, OUTDIM> out; ///< The resulting output vector
};

/** Gaussian function */
static inline double gauss(const double sigma, const double distsq)
;//__attribute__((const, fastcall, nothrow, unused));

/** gaussian distance for floats */
static inline float gauss(const float sigma, const float distsq)
;//__attribute__((const, fastcall, nothrow, unused));

/** Gaussian "distance" between two points -- float version */
template <class T, int DIM>
static inline T gdist(const T *const v1, const T *const v2, const T sigma)
;//__attribute__((regparm(3), nonnull, nothrow, pure, unused));


/** Gaussian "distance" between two points -- double version */
template <class T, int DIM>
static inline T dgdist(const T *const v1, const T *const v2, const T sigma)
;//__attribute__((regparm(3), nonnull, nothrow, pure, unused));

/** The euler distance between two poitns */
template <class T, int DIM>
static inline T distsq(const T* const v1, const T *const v2) {
    T sum = 0;
    for (unsigned i = 0; i < DIM; ++i) {
        const T d = v1[i] - v2[i];
        sum += d*d;
    }
    return sum;
}

/** Gaussian distribution functor */
template <class T, int DIM>
struct GDist {
    /** Calls gdist() */
    T operator() (const T *const v1, const T *const v2, const T sigma) const {
        return gdist<T, DIM>(v1, v2, sigma);
    }
};

/** Euler/linear distribution functor */
template <class T, int DIM>
struct EulerDistSq {
    /** Calls sqrt() of distsq() */
    T operator() (const T *const v1, const T *const v2, const T /*sigma*/) const {
        return sqrt(distsq<T, DIM>(v1, v2));
    }
};


/**
 * \brief An Exact-interpolation radial basis function (RBF) neural network (NN).
 *
 * The interpolation task:
 * <ul>
 *   <li>Given</li>
 *   <ol>
 *     <li>A black box with \f$n\f$ inputs \f$x_1, x_2, \dots, x_n\f$ and
 *        an \f$m\f$ outputs \f$y_1, y_2, \dots, y_m\f$</li>
 *     <li>A database of \f$k\f$ sample input-output \f$(\boldsymbol{x}, \boldsymbol{y}(\boldsymbol{x}))\f$
 *        combinations (training examples)</li>
 *   </ol>
 *   <li>Goal: Create a NN that predicts the output given an input vector \f$\boldsymbol{x}\f$</li>
 *   <ul>
 *     <li>if \f$\boldsymbol{x}\f$ is a vector from the sample data, it must
 *       produce the exact output of the black box</li>
 *     <li>otherwise it must produce output clost to the output of the black
 *        box (i.e. a result from the [most likely non-linear] interpolation of
 *        the training examples.)</li>
 *   </ul>
 * </ul>
 *
 * \author Trent Apted <tapted@it.usyd.edu.au>
 *
 * \param T The type of the floating point number to use. Currently \a double and
 *          \a float are supported. \a complex is also a possibility.
 * \param INDIM The dimension of the domain (inputs)
 * \param OUTDIM The dimension of the range (outputs)
 */
template <class T, int INDIM, int OUTDIM, class DISTFUNC = GDist<T, INDIM> >
class ExactRBFNet {
public:
    typedef MyVec<T, INDIM> INVEC;   ///< An input vector -- \f$\boldsymbol{x}\f$
    typedef MyVec<T, OUTDIM> OUTVEC; ///< An output vector -- \f$\boldsymbol{y}(\boldsymbol{x})\f$
    typedef Example<T, INDIM, OUTDIM> EXAMPLE; ///< An example \f$(\boldsymbol{x}, \boldsymbol{y}(\boldsymbol{x}))\f$
    typedef stl::vector<EXAMPLE> EXAMPLEVEC;   ///< A collection of examples/centres
#ifndef RBFNET_DEBUG
protected:
#endif
    EXAMPLEVEC examples;              ///< The examples given AND the "centres" of the gaussians
    std::valarray<T> weights[OUTDIM]; ///< The weights vector -- \f$\boldsymbol{w}\f$
    T sigma;                          ///< The \f$\sigma\f$ value to use for the gaussian function
    DISTFUNC distfunc;                ///< The distance function
    bool badbit;                      ///< Set if/when there are any read errors
public:
    /** Construct a (previously output) ExactRBFNet from an input stream. */
    ExactRBFNet(std::istream& in, const DISTFUNC &df = DISTFUNC()) : sigma(0), distfunc(df), badbit(false) { input(in); }

    /** Construct a (previously output) ExactRBFNet from a file. */
    ExactRBFNet(const char* file, const DISTFUNC &df = DISTFUNC()) : sigma(0), distfunc(df), badbit(false) { std::ifstream f(file); input(f); }

#ifndef RBFNET_READ_ONLY
    /**
     * Construct a new ExactRBFNet with a given value for \f$\sigma\f$. \f$\sigma\f$ is the
     * width / influence of the Gaussian function. A larger \f$\sigma\f$ means that each
     * example will have a greater influence over the domain -- each sub() will be influenced
     * by a greater number of training examples.
     *
     * \param sig \a sigma, \f$\sigma\f$
     * \param df The distance function to use
     */
    ExactRBFNet(const T& sig = 4, const DISTFUNC &df = DISTFUNC()) : sigma(sig), distfunc(df), badbit(false) {}

    /**
     * Add a training example / centre.
     */
    void addExample(const EXAMPLE &ex);

    /**
     * \brief Solve this RBF Network for the weights.
     *
     * That is, solve the set of linear equations:
     *
     * \f{eqnarray*}
     *   y(\boldsymbol{x}_1) & = & w_1\mathrm{e}^\frac{\scriptstyle \|\boldsymbol{x}_1-\boldsymbol{x}_1\|^2}{\scriptstyle 2\sigma} +
     *                             w_2\mathrm{e}^\frac{\scriptstyle \|\boldsymbol{x}_1-\boldsymbol{x}_2\|^2}{\scriptstyle 2\sigma} + \cdots +
     *                             w_k\mathrm{e}^\frac{\scriptstyle \|\boldsymbol{x}_1-\boldsymbol{x}_k\|^2}{\scriptstyle 2\sigma} \\
     *   y(\boldsymbol{x}_2) & = & w_1\mathrm{e}^\frac{\scriptstyle \|\boldsymbol{x}_2-\boldsymbol{x}_1\|^2}{\scriptstyle 2\sigma} +
     *                             w_2\mathrm{e}^\frac{\scriptstyle \|\boldsymbol{x}_2-\boldsymbol{x}_2\|^2}{\scriptstyle 2\sigma} + \cdots +
     *                             w_k\mathrm{e}^\frac{\scriptstyle \|\boldsymbol{x}_2-\boldsymbol{x}_k\|^2}{\scriptstyle 2\sigma} \\
     *    & \vdots & \\
     *   y(\boldsymbol{x}_k) & = & w_1\mathrm{e}^\frac{\scriptstyle \|\boldsymbol{x}_k-\boldsymbol{x}_1\|^2}{\scriptstyle 2\sigma} +
     *                             w_2\mathrm{e}^\frac{\scriptstyle \|\boldsymbol{x}_k-\boldsymbol{x}_2\|^2}{\scriptstyle 2\sigma} + \cdots +
     *                             w_k\mathrm{e}^\frac{\scriptstyle \|\boldsymbol{x}_k-\boldsymbol{x}_k\|^2}{\scriptstyle 2\sigma}
     * \f}
     * for \f$\boldsymbol{w}\f$ for each of the \f$\boldsymbol{m}\f$ outputs (i.e. for each set of \f$y\f$ s).
     *
     * This needs to be called before sub() if there have been any calls to addExample().
     *
     * If RBFNET_DEBUG is defined, the matrix passed to lapack is printed on stderr.
     * If RBFNET_VDEBUG is defined, the formula substitutions are also printed.
     */
    int solve();

    /**
     * Set a new value of \f$\sigma\f$. Afterwards, solve() needs to be called again.
     */
    void setSigma(const T &newsig) {sigma = newsig;}
#endif
    /**
     * Return the number of examples == the dimension of each weight vector, \f$\boldsymbol{w}\f$.
     */
    size_t dim() {return examples.size();}

    /**
     * Substitute an input into the NN, for the interpolated output.
     * \param out The output that is generated (assigned by reference)
     * \param in The input to present
     */
    void sub(OUTVEC &out, const INVEC &in) const ;
    /**
     * A debugging version of sub(), that shows the formula substitution on
     * sterr if RBFNET_VDEBUG is defined.
     */
    void dsub(OUTVEC &out, const INVEC &in) const ;

    /**
     * Output a human and machine readable representation of the network to \a o.
     */
    std::ostream &output(std::ostream& o) const;

    /**
     * Read in a previously output() representation of the network from \a in.
     */
    std::istream &input(std::istream& i);

    /**
     * Are we ready for use (i.e. no read errors from the fstream constructor)
     */
    bool ready() { return !badbit; }
};

/** Returns rhs.output(o) */
template <class T, int INDIM, int OUTDIM, class DISTFUNC>
std::ostream& operator<< (std::ostream &o, const ExactRBFNet<T, INDIM, OUTDIM, DISTFUNC> &rhs) {
    return rhs.output(o);
}

/** Returns rhs.input(i) */
template <class T, int INDIM, int OUTDIM, class DISTFUNC>
std::istream& operator>> (std::istream &i, const ExactRBFNet<T, INDIM, OUTDIM, DISTFUNC> &rhs) {
    return rhs.input(i);
}

/**
 * A template instantiation for 2D to 2D interpolation (e.g. calibrating a
 * projected touch display).
 */
typedef ExactRBFNet<double, 2, 2> XYInterp;

/**
 * A template instantiation of a linear interpolation network (works poorly)
 */
typedef ExactRBFNet<double, 2, 2, EulerDistSq<double, 2> > XYLinInterp;


//IMPLEMENTATION
#ifndef RBFNET_READ_ONLY
#include "lapack_leqn.h"
#endif

#include <stdio.h>
#include <math.h>

static inline double gauss(const double sigma, const double distsq) {
    return exp((distsq)/(2.0*sigma));
}

static inline float gauss(const float sigma, const float distsq) {
    return expf((distsq)/(2.0*sigma));
}

template <class T, int DIM>
static inline T gdist(const T *const v1, const T *const v2, const T sigma) {
    return gauss(sigma, distsq<T, DIM>(v1, v2));
}

#ifndef RBFNET_READ_ONLY
template <class T, int INDIM, int OUTDIM, class DISTFUNC>
void ExactRBFNet<T, INDIM, OUTDIM, DISTFUNC>::addExample(const EXAMPLE &ex) {
    examples.push_back(ex);
}

template <class T, int INDIM, int OUTDIM, class DISTFUNC>
int ExactRBFNet<T, INDIM, OUTDIM, DISTFUNC>::solve() {
    int sz = examples.size();
    int oks = 0;
    T *A = new T[sz*sz];
    T *R = new T[sz];
    for (int r = 0; r < sz; ++r) {
#ifdef RBFNET_VDEBUG
        fprintf(stderr, "\ny = ");
#endif
        int c = 0;
#ifdef RBFNET_EXPONE_OPTIMIZATION
        for (; c < r; ++c) {
            A[r*sz + c] = distfunc(&examples[r].in[0], &examples[c].in[0], sigma);
        }
        A[r*sz + c++] = 1.0; //e^0 == 1
#endif
        for (; c < sz; ++c) {
            A[r*sz + c] = distfunc(&examples[r].in[0], &examples[c].in[0], sigma);
        }
    }
    for (int o = 0; o < OUTDIM; ++o) {
        for (int r = 0; r < sz; ++r)
            R[r] = examples[r].out[o];
#ifdef RBFNET_DEBUG
        fprintf(stderr, "\nSolving :\n");
        for (int _r = 0; _r < sz; ++_r) {
            fprintf(stderr, "[ %8f", A[_r*sz]);
            for (int _c = 1; _c < sz; ++ _c) {
                fprintf(stderr, ", %8f", A[_r*sz + _c]);
            }
            fprintf(stderr, "][%c] = [%8f]\n", 'a'+_r, R[_r]);
        }
#endif
        oks += lapack_solve(R, A, sz);
#ifdef RBFNET_DEBUG
        if (oks) {
            fprintf(stderr, "\nNo Solution\n");
            oks = 0;
        } else {
            fputc('\n', stderr);
            for (int _r = 0; _r < sz; ++_r) {
                fprintf(stderr, "[%c] = [%8f]\n", 'a'+_r, R[_r]);
            }
        }
#endif
        weights[o].resize(sz);
        //weights[o] = std::valarray<T>(R, sz);
        for (int r = 0; r < sz; ++r)
            weights[o][r] = R[r];
    }
    delete[] R;
    delete[] A;
    return oks;
}
#endif //#ifndef RBFNET_READ_ONLY

template <class T, int INDIM, int OUTDIM, class DISTFUNC>
void ExactRBFNet<T, INDIM, OUTDIM, DISTFUNC>::sub(OUTVEC &out, const INVEC &in) const {
    int sz = examples.size();
    for (int o = 0; o < OUTDIM; ++o) {
        out.v[o] = weights[o][0] * distfunc(&in[0], &examples[0].in[0], sigma);
        for (int i = 1; i < sz; ++i) {
            out.v[o] += weights[o][i] * distfunc(&in[0], &examples[i].in[0], sigma);
        }
    }
}

template <class T, int INDIM, int OUTDIM, class DISTFUNC>
void ExactRBFNet<T, INDIM, OUTDIM, DISTFUNC>::dsub(OUTVEC &out, const INVEC &in) const {
    int sz = examples.size();
    for (int o = 0; o < OUTDIM; ++o) {
#ifdef RBFNET_VDEBUG
        fprintf(stderr, "\ny = (%f)", weights[o][0]);
#endif
        out.v[o] = weights[o][0] * distfunc(&in[0], &examples[0].in[0], sigma);
        for (int i = 1; i < sz; ++i) {
#ifdef RBFNET_VDEBUG
            fprintf(stderr, "(%f)", weights[o][i]);
#endif
            out.v[o] += weights[o][i] * distfunc(&in[0], &examples[i].in[0], sigma);
        }
    }
}

template <class T, int DIM>
static inline T dgdist(const T *const v1, const T *const v2, const T sigma) {
    T sum = 0;
    for (unsigned i = 0; i < DIM; ++i) {
        const T d = v1[i] - v2[i];
        sum += d*d;
    }
#ifdef RBFNET_VDEBUG
    fprintf(stderr, "we^(%f/2*%f) ", sum, sigma);
#endif
    return gauss(sigma, sum);
}

/** Output a MyVec */
template <class T, int DIM>
std::ostream &operator<< (std::ostream &o, const MyVec<T, DIM> &v) {
    o << v[0];
    for (int i = 1; i < DIM; ++i)
        o << ' ' << v[i];
    return o;
}

/** Input a MyVec */
template <class T, int DIM>
std::istream &operator>> (std::istream &i, MyVec<T, DIM> &v) {
    i >> v[0];
    for (int j = 1; j < DIM; ++j)
        i >> v[j];
    return i;
}

/** Output an Example */
template <class T, int INDIM, int OUTDIM>
std::ostream &operator<< (std::ostream &o, const Example<T, INDIM, OUTDIM> &e) {
    return o << e.in << "\n" << e.out;
}

/** Input an Example */
template <class T, int INDIM, int OUTDIM>
std::istream &operator>> (std::istream &i, Example<T, INDIM, OUTDIM> &e) {
    return i >> e.in >> e.out;
}

template <class T, int INDIM, int OUTDIM, class DISTFUNC>
std::ostream &ExactRBFNet<T, INDIM, OUTDIM, DISTFUNC>::output(std::ostream& o) const {
    o.flags(std::ios_base::scientific);
    o.precision(15);
    o << "ExactRBFNet< " << typeid(T).name() << " , " << INDIM << " , " << OUTDIM << " >\n"
        << "sigma= " << sigma << "\n";
    o << "weights==examples.size()= " << examples.size();
    for (int p = 0; p < OUTDIM; ++p) {
        o << "\nweights[" << p << "]= ";
        for (size_t i = 0; i < weights[p].size(); ++i)
            o << weights[p][i] << " ";
    }
    for (typename EXAMPLEVEC::const_iterator it = examples.begin(); it != examples.end(); ++it) {
        o << "\n\n" << *it;
    }
    return o;
}

template <class T, int INDIM, int OUTDIM, class DISTFUNC>
std::istream &ExactRBFNet<T, INDIM, OUTDIM, DISTFUNC>::input(std::istream& i) {
    std::string d, type;
    int indim, outdim, size;
    i >> d >> type >> d >> indim >> d >> outdim >> d >> d >> sigma;
    i >> d >> size;
    if (i) {
        for (int o = 0; i && o < OUTDIM; ++o) {
            i >> d;
            weights[o].resize(size);
            for (int j = 0; i && j < size; ++j)
                i >> weights[o][j];
        }
        examples.reserve(size);
        for (int j = 0; i && j < size; ++j) {
            EXAMPLE eg;
            i >> eg;
            examples.push_back(eg);
        }
    }
    if (!i)
        badbit = true;
    return i;
}

#endif //#ifndef RBN_SOLVE_DOT_AITCH