drwnFeatureMaps.h

/******************************************************************************
** DARWIN: A FRAMEWORK FOR MACHINE LEARNING RESEARCH AND DEVELOPMENT
** Distributed under the terms of the BSD license (see the LICENSE file)
** Copyright (c) 2007-2015, Stephen Gould
** All rights reserved.
**
******************************************************************************
** FILENAME:    drwnFeatureMaps.h
** AUTHOR(S):   Stephen Gould <stephen.gould@anu.edu.au>
**
*****************************************************************************/

#pragma once

#include <cstdlib>
#include <vector>
#include <cmath>

using namespace std;

#define DRWN_UNROLL_FEATURE_MAP_LOOPS

//
// classless feature mappings ----------------------------------------------
//

// drwnFeatureMap ----------------------------------------------------------

class drwnFeatureMap : public drwnTypeable, public drwnCloneable {
 protected:
    int _nFeatures; 

 public:
    drwnFeatureMap() : _nFeatures(0) { /* do nothing */ }
    drwnFeatureMap(int nFeatures) : _nFeatures(nFeatures) { /* do nothing */ }
    virtual ~drwnFeatureMap() { /* do nothing */ }

    virtual drwnFeatureMap *clone() const = 0;

    int numFeatures() const { return _nFeatures; }
    virtual int numParameters() const = 0;

    virtual void initialize(int nFeatures) {
        _nFeatures = nFeatures;
    }

    virtual vector<double> operator()(const vector<double>& x) const = 0;

    virtual double dot(const vector<double>& theta, const vector<double>& x) const {
        vector<double> phi((*this)(x));
        double v = 0.0;
        for (unsigned i = 0; i < phi.size(); i++) {
            v += phi[i] * theta[i];
        }
        return v;
    }

    virtual void mac(vector<double>& theta, const vector<double>& x, double alpha) const {
        const vector<double> phi((*this)(x));
        for (unsigned i = 0; i < phi.size(); i++) {
            theta[i] += alpha * phi[i];
        }
    }
};

// drwnIdentityFeatureMap ---------------------------------------------------

class drwnIdentityFeatureMap : public drwnFeatureMap {
 public:
    drwnIdentityFeatureMap() : drwnFeatureMap() { /* do nothing */ }
    drwnIdentityFeatureMap(int nFeatures) : drwnFeatureMap(nFeatures) { /* do nothing */ }
    ~drwnIdentityFeatureMap() { /* do nothing */ }

    // type and cloning
    const char *type () const { return "drwnIdentityFeatureMap"; }
    drwnIdentityFeatureMap *clone() const { return new drwnIdentityFeatureMap(*this); }

    int numParameters() const { return _nFeatures; }

    vector<double> operator()(const vector<double>& x) const {
        return x;
    }

    double dot(const vector<double>& theta, const vector<double>& x) const {
        return Eigen::Map<const VectorXd>(&theta[0], 
            theta.size()).dot(Eigen::Map<const VectorXd>(&x[0], x.size()));
    }

    void mac(vector<double>& theta, const vector<double>& x, double alpha) const {
        Eigen::Map<VectorXd>(&theta[0], theta.size()) += alpha *
            Eigen::Map<const VectorXd>(&x[0], x.size());
    }
};

// drwnBiasFeatureMap ------------------------------------------------------

class drwnBiasFeatureMap : public drwnFeatureMap {
 public:
    drwnBiasFeatureMap() : drwnFeatureMap() { /* do nothing */ }
    drwnBiasFeatureMap(int nFeatures) : drwnFeatureMap(nFeatures) { /* do nothing */ }
    ~drwnBiasFeatureMap() { /* do nothing */ }

    // type and cloning
    const char *type () const { return "drwnBiasFeatureMap"; }
    drwnBiasFeatureMap *clone() const { return new drwnBiasFeatureMap(*this); }

    int numParameters() const { return _nFeatures + 1; }

    vector<double> operator()(const vector<double>& x) const {
        vector<double> phi(x.size() + 1);
        copy(x.begin(), x.end(), phi.begin());
        phi.back() = 1.0;
        return phi;
    }

    double dot(const vector<double>& theta, const vector<double>& x) const {
        return theta.back() +
            Eigen::Map<const VectorXd>(&theta[0], 
                theta.size() - 1).dot(Eigen::Map<const VectorXd>(&x[0], x.size()));
    }

    void mac(vector<double>& theta, const vector<double>& x, double alpha) const {
        Eigen::Map<VectorXd>(&theta[0], theta.size() - 1) += alpha *
            Eigen::Map<const VectorXd>(&x[0], x.size());
        theta.back() += alpha;
    }
};

// drwnSquareFeatureMap ----------------------------------------------------

class drwnSquareFeatureMap : public drwnFeatureMap {
 public:
    drwnSquareFeatureMap() : drwnFeatureMap() { /* do nothing */ }
    drwnSquareFeatureMap(int nFeatures) : drwnFeatureMap(nFeatures) { /* do nothing */ }
    ~drwnSquareFeatureMap() { /* do nothing */ }

    // type and cloning
    const char *type () const { return "drwnSquareFeatureMap"; }
    drwnSquareFeatureMap *clone() const { return new drwnSquareFeatureMap(*this); }

    int numParameters() const { return 2 * _nFeatures + 1; }

    vector<double> operator()(const vector<double>& x) const {
        vector<double> phi(2 * x.size() + 1);
        copy(x.begin(), x.end(), phi.begin());
        for (unsigned i = 0; i < x.size(); i++) {
            phi[x.size() + i] = M_SQRT1_2 * (x[i] * x[i] - 1.0);
        }
        phi.back() = 1.0;
        return phi;
    }

    double dot(const vector<double>& theta, const vector<double>& x) const {
        double v = theta.back();
        vector<double>::const_iterator it = theta.begin();
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            v += (*it) * (*ix);
        }
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            v += M_SQRT1_2 * (*it) * ((*ix) * (*ix) - 1.0);
        }
        return v;
    }

    void mac(vector<double>& theta, const vector<double>& x, double alpha) const {
        Eigen::Map<ArrayXd>(&theta[0], _nFeatures) += alpha *
            Eigen::Map<const ArrayXd>(&x[0], x.size());
        Eigen::Map<ArrayXd>(&theta[_nFeatures], _nFeatures) += M_SQRT1_2 * alpha *
            (Eigen::Map<const ArrayXd>(&x[0], x.size()).array().square() - 1.0);
        theta.back() += alpha;
    }
};

// drwnQuadraticFeatureMap -------------------------------------------------

class drwnQuadraticFeatureMap : public drwnFeatureMap {
 public:
    drwnQuadraticFeatureMap() : drwnFeatureMap() { /* do nothing */ }
    drwnQuadraticFeatureMap(int nFeatures) : drwnFeatureMap(nFeatures) { /* do nothing */ }
    ~drwnQuadraticFeatureMap() { /* do nothing */ }

    // type and cloning
    const char *type () const { return "drwnQuadraticFeatureMap"; }
    drwnQuadraticFeatureMap *clone() const { return new drwnQuadraticFeatureMap(*this); }

    int numParameters() const { return (_nFeatures + 3) * _nFeatures / 2 + 1; }

    vector<double> operator()(const vector<double>& x) const {
        vector<double> phi(2 * x.size() + 1);
        copy(x.begin(), x.end(), phi.begin());
        int indx = x.size();
        for (unsigned i = 0; i < x.size(); i++) {
            for (unsigned j = 0; j < i; j++) {
                phi[indx++] = x[i] * x[j];
            }
            phi[indx++] = M_SQRT1_2 * (x[i] * x[i] - 1.0);
        }
        phi.back() = 1.0;
        return phi;
    }

    double dot(const vector<double>& theta, const vector<double>& x) const {
        double v = theta.back();
        vector<double>::const_iterator it = theta.begin();
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            v += (*it) * (*ix);
        }
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            for (vector<double>::const_iterator jx = x.begin(); jx != ix; ++jx, ++it) {
                v += (*it) * (*ix) * (*jx);
            }
            v += M_SQRT1_2 * (*it) * ((*ix) * (*ix) - 1.0);
        }
        return v;
    }

    void mac(vector<double>& theta, const vector<double>& x, double alpha) const {
        vector<double>::iterator it = theta.begin();
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            (*it) += alpha * (*ix);
        }
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            for (vector<double>::const_iterator jx = x.begin(); jx != ix; ++jx, ++it) {
                (*it) += alpha * (*ix) * (*jx);
            }
            (*it) += alpha * M_SQRT1_2 * ((*ix) * (*ix) - 1.0);
        }
        theta.back() += alpha;
    }
};

//
// classful feature mappings -----------------------------------------------
//

// drwnJointFeatureMap -----------------------------------------------------

class drwnJointFeatureMap : public drwnTypeable, public drwnCloneable {
 protected:
    int _nFeatures; 
    int _nClasses;  

 public:
    drwnJointFeatureMap() : _nFeatures(0), _nClasses(0) { /* do nothing */ }
    drwnJointFeatureMap(int nFeatures, int nClasses) :
        _nFeatures(nFeatures), _nClasses(nClasses) { /* do nothing */ }
    virtual ~drwnJointFeatureMap() { /* do nothing */ }

    virtual drwnJointFeatureMap *clone() const = 0;

    int numFeatures() const { return _nFeatures; }
    int numClasses() const { return _nClasses; }
    virtual int numParameters() const = 0;

    virtual void initialize(int nFeatures, int nClasses) {
        _nFeatures = nFeatures; _nClasses = nClasses;
    }

    virtual vector<double> operator()(const vector<double>& x, int y) const = 0;

    virtual vector<double> operator()(const vector<double>& x) const {
        vector<double> phi(numParameters());
        for (int y = 0; y < _nClasses; y++) {
            vector<double> phi_y((*this)(x, y));
            for (unsigned i = 0; i < phi.size(); i++) {
                phi[i] += phi_y[i];
            }
        }
        return phi;
    }

    virtual double dot(const vector<double>& theta, const vector<double>& x, int y) const {
        vector<double> phi((*this)(x, y));
        double v = 0.0;
        for (unsigned i = 0; i < phi.size(); i++) {
            v += phi[i] * theta[i];
        }
        return v;
    }

    double dot(const vector<double>& theta, const vector<double>& x) const {
        double v = 0.0;
        for (int y = 0; y < _nClasses; y++) {
            v += dot(theta, x, y);
        }
        return v;
    }

    virtual void mac(vector<double>& theta, const vector<double>& x, double alpha, int y) const {
        const vector<double> phi((*this)(x, y));
        for (unsigned i = 0; i < phi.size(); i++) {
            theta[i] += alpha * phi[i];
        }
    }

    virtual void mac(vector<double>& theta, const vector<double>& x, const vector<double>& alpha) const {
        for (int y = 0; y < _nClasses; y++) {
            const vector<double> phi((*this)(x, y));
            for (unsigned i = 0; i < phi.size(); i++) {
                theta[i] += alpha[y] * phi[i];
            }
        }
    }
};

// drwnIdentityJointFeatureMap ----------------------------------------------

class drwnIdentityJointFeatureMap : public drwnJointFeatureMap {
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
 protected:
    int _nFeaturesDiv4;
    int _nFeaturesMod4;
#endif

 public:
    drwnIdentityJointFeatureMap() : drwnJointFeatureMap() {
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
        _nFeaturesDiv4 = 0;
        _nFeaturesMod4 = 0;
#endif
    }
    drwnIdentityJointFeatureMap(int nFeatures, int nClasses) :
        drwnJointFeatureMap(nFeatures, nClasses) {
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
            _nFeaturesDiv4 = nFeatures / 4;
            _nFeaturesMod4 = nFeatures % 4;
#endif
    }
    ~drwnIdentityJointFeatureMap() { /* do nothing */ }

    // type and cloning
    const char *type () const { return "drwnIdentityJointFeatureMap"; }
    drwnIdentityJointFeatureMap *clone() const { return new drwnIdentityJointFeatureMap(*this); }

    int numParameters() const { return std::max(_nFeatures * (_nClasses - 1), 0); }

#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
    void initialize(int nFeatures, int nClasses) {
        _nFeatures = nFeatures; _nClasses = nClasses;
        _nFeaturesDiv4 = nFeatures / 4;
        _nFeaturesMod4 = nFeatures % 4;
    }
#endif

    vector<double> operator()(const vector<double>& x, int y) const {
        vector<double> phi(numParameters(), 0.0);
        if (y != _nClasses - 1) {
            copy(x.begin(), x.end(), phi.begin() + y * _nFeatures);
        }
        return phi;
    }

    vector<double> operator()(const vector<double>& x) const {
        vector<double> phi(numParameters());
        for (int y = 0; y < _nClasses - 1; y++) {
            copy(x.begin(), x.end(), phi.begin() + y * _nFeatures);
        }
        return phi;
    }

    double dot(const vector<double>& theta, const vector<double>& x, int y) const {
        if (y == _nClasses - 1) return 0.0;
        double v = 0.0;
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
        const double *it = &theta[y * _nFeatures];
        const double *ix = &x[0];
        for (int i = _nFeaturesDiv4; i != 0; i--) {
            v += it[0] * ix[0] + it[1] * ix[1] + it[2] * ix[2] + it[3] * ix[3];
            it += 4; ix += 4;
        }
        for (int i = 0; i < _nFeaturesMod4; i++, ++it, ++ix) {
            v += (*it) * (*ix);
        }
#else
        vector<double>::const_iterator it = theta.begin() + y * _nFeatures;
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            v += (*it) * (*ix);
        }
#endif
        return v;
    }

    double dot(const vector<double>& theta, const vector<double>& x) const {
        double v = 0.0;
        vector<double>::const_iterator it = theta.begin();
        for (int y = 0; y < _nClasses - 1; y++) {
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                v += (*it) * (*ix);
            }
        }
        return v;
    }

    void mac(vector<double>& theta, const vector<double>& x, double alpha, int y) const {
        if (y == _nClasses - 1) return;
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
        double *it = &theta[y * _nFeatures];
        const double *ix = &x[0];
        for (int i = _nFeaturesDiv4; i != 0; i--) {
            it[0] += alpha * ix[0];
            it[1] += alpha * ix[1];
            it[2] += alpha * ix[2];
            it[3] += alpha * ix[3];
            it += 4; ix += 4;
        }
        for (int i = 0; i < _nFeaturesMod4; i++, ++it, ++ix) {
            (*it) += alpha * (*ix);
        }
#else
        vector<double>::iterator it = theta.begin() + y * _nFeatures;
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            (*it) += alpha * (*ix);
        }
#endif
    }

    void mac(vector<double>& theta, const vector<double>& x, const vector<double>& alpha) const {
        vector<double>::iterator it = theta.begin();
        for (int y = 0; y < _nClasses - 1; y++) {
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                (*it) += alpha[y] * (*ix);
            }
        }
    }
};

// drwnBiasJointFeatureMap -------------------------------------------------

class drwnBiasJointFeatureMap : public drwnJointFeatureMap {
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
 protected:
    int _nFeaturesDiv4;
    int _nFeaturesMod4;
#endif

 public:
    drwnBiasJointFeatureMap() : drwnJointFeatureMap() {
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
        _nFeaturesDiv4 = 0;
        _nFeaturesMod4 = 0;
#endif
    }
    drwnBiasJointFeatureMap(int nFeatures, int nClasses) :
        drwnJointFeatureMap(nFeatures, nClasses) {
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
            _nFeaturesDiv4 = nFeatures / 4;
            _nFeaturesMod4 = nFeatures % 4;
#endif
    }
    ~drwnBiasJointFeatureMap() { /* do nothing */ }

    // type and cloning
    const char *type () const { return "drwnBiasJointFeatureMap"; }
    drwnBiasJointFeatureMap *clone() const { return new drwnBiasJointFeatureMap(*this); }

    int numParameters() const { return std::max((_nFeatures + 1) * (_nClasses - 1), 0); }

#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
    void initialize(int nFeatures, int nClasses) {
        _nFeatures = nFeatures; _nClasses = nClasses;
        _nFeaturesDiv4 = nFeatures / 4;
        _nFeaturesMod4 = nFeatures % 4;
    }
#endif

    vector<double> operator()(const vector<double>& x, int y) const {
        vector<double> phi(numParameters(), 0.0);
        if (y != _nClasses - 1) {
            copy(x.begin(), x.end(), phi.begin() + y * (_nFeatures + 1));
            phi[y * (_nFeatures + 1) + _nFeatures] = 1.0;
        }
        return phi;
    }

    vector<double> operator()(const vector<double>& x) const {
        vector<double> phi(numParameters());
        vector<double>::iterator it = phi.begin();
        for (int y = 0; y < _nClasses - 1; y++) {
            copy(x.begin(), x.end(), it);
            *(it += _nFeatures)++ = 1.0;
        }
        return phi;
    }

    double dot(const vector<double>& theta, const vector<double>& x, int y) const {
        if (y == _nClasses - 1) return 0.0;
        double v = 0.0;
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
        const double *it = &theta[y * (_nFeatures + 1)];
        const double *ix = &x[0];
        for (int i = _nFeaturesDiv4; i != 0; i--) {
            v += it[0] * ix[0] + it[1] * ix[1] + it[2] * ix[2] + it[3] * ix[3];
            it += 4; ix += 4;
        }
        for (int i = 0; i < _nFeaturesMod4; i++, ++it, ++ix) {
            v += (*it) * (*ix);
        }
#else
        vector<double>::const_iterator it = theta.begin() + y * (_nFeatures + 1);
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            v += (*it) * (*ix);
        }
#endif
        v += (*it);
        return v;
    }

    double dot(const vector<double>& theta, const vector<double>& x) const {
        double v = 0.0;
        vector<double>::const_iterator it = theta.begin();
        for (int y = 0; y < _nClasses - 1; y++) {
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                v += (*it) * (*ix);
            }
            v += (*it++);
        }
        return v;
    }

    void mac(vector<double>& theta, const vector<double>& x, double alpha, int y) const {
        if (y == _nClasses - 1) return;
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
        double *it = &theta[y * (_nFeatures + 1)];
        const double *ix = &x[0];
        for (int i = _nFeaturesDiv4; i != 0; i--) {
            it[0] += alpha * ix[0];
            it[1] += alpha * ix[1];
            it[2] += alpha * ix[2];
            it[3] += alpha * ix[3];
            it += 4; ix += 4;
        }
        for (int i = 0; i < _nFeaturesMod4; i++, ++it, ++ix) {
            (*it) += alpha * (*ix);
        }
#else
        vector<double>::iterator it = theta.begin() + y * (_nFeatures + 1);
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            (*it) += alpha * (*ix);
        }
#endif
        (*it) += alpha;
    }

    void mac(vector<double>& theta, const vector<double>& x, const vector<double>& alpha) const {
        vector<double>::iterator it = theta.begin();
        for (int y = 0; y < _nClasses - 1; y++) {
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                (*it) += alpha[y] * (*ix);
            }
            (*it++) += alpha[y];
        }
    }
};

// drwnSquareJointFeatureMap -----------------------------------------------

class drwnSquareJointFeatureMap : public drwnJointFeatureMap {
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
 protected:
    int _nFeaturesDiv4;
    int _nFeaturesMod4;
#endif

 public:
    drwnSquareJointFeatureMap() : drwnJointFeatureMap() {
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
        _nFeaturesDiv4 = 0;
        _nFeaturesMod4 = 0;
#endif
    }
    drwnSquareJointFeatureMap(int nFeatures, int nClasses) :
        drwnJointFeatureMap(nFeatures, nClasses) {
#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
            _nFeaturesDiv4 = nFeatures / 4;
            _nFeaturesMod4 = nFeatures % 4;
#endif
    }
    ~drwnSquareJointFeatureMap() { /* do nothing */ }

    // type and cloning
    const char *type () const { return "drwnSquareJointFeatureMap"; }
    drwnSquareJointFeatureMap *clone() const { return new drwnSquareJointFeatureMap(*this); }

    int numParameters() const { return std::max((2 * _nFeatures + 1) * (_nClasses - 1), 0); }

#ifdef DRWN_UNROLL_FEATURE_MAP_LOOPS
    void initialize(int nFeatures, int nClasses) {
        _nFeatures = nFeatures; _nClasses = nClasses;
        _nFeaturesDiv4 = nFeatures / 4;
        _nFeaturesMod4 = nFeatures % 4;
    }
#endif

    vector<double> operator()(const vector<double>& x, int y) const {
        vector<double> phi(numParameters(), 0.0);
        if (y != _nClasses - 1) {
            vector<double>::iterator jt = phi.begin() + y * (2 * _nFeatures + 1);
            for (vector<double>::const_iterator it = x.begin(); it != x.end(); ++it, ++jt) {
                (*jt) = (*it);
            }
            for (vector<double>::const_iterator it = x.begin(); it != x.end(); ++it, ++jt) {
                (*jt) = M_SQRT1_2 * ((*it) * (*it) - 1.0);
            }
            (*jt) = 1.0;
        }
        return phi;
    }

    vector<double> operator()(const vector<double>& x) const {
        vector<double> phi(numParameters());
        vector<double>::iterator jt = phi.begin();
        for (vector<double>::const_iterator it = x.begin(); it != x.end(); ++it, ++jt) {
            (*jt) = (*it);
        }
        for (vector<double>::const_iterator it = x.begin(); it != x.end(); ++it, ++jt) {
            (*jt) = M_SQRT1_2 * ((*it) * (*it) - 1.0);
        }
        (*jt) = 1.0;

        for (int y = 1; y < _nClasses - 1; y++) {
            copy(jt - (2 * _nFeatures + 1), jt, jt + 1);
            jt += (2 * _nFeatures + 1);
        }
        return phi;
    }

    double dot(const vector<double>& theta, const vector<double>& x, int y) const {
        if (y == _nClasses - 1) return 0.0;
#if 1
        double v = 0.0;
        vector<double>::const_iterator it = theta.begin() + y * (2 * _nFeatures + 1);
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            v += (*it) * (*ix);
        }
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            v += M_SQRT1_2 * (*it) * ((*ix) * (*ix) - 1.0);
        }
        v += (*it);

        return v;
#else
        return Eigen::Map<const VectorXd>(&theta[y * (2 * _nFeatures + 1)],
            _nFeatures).dot(Eigen::Map<const VectorXd>(&x[0], _nFeatures)) +
            Eigen::Map<const VectorXd>(&theta[y * (2 * _nFeatures + 1) + _nFeatures],
                _nFeatures).dot(M_SQRT1_2 *
                    (Eigen::Map<const VectorXd>(&x[0], x.size()).array().square().array() - 1.0)) +
            theta[y * (2 * _nFeatures + 1) + 2 * _nFeatures];
#endif
    }

    double dot(const vector<double>& theta, const vector<double>& x) const {
        double v = 0.0;
        vector<double>::const_iterator it = theta.begin();
        for (int y = 0; y < _nClasses - 1; y++) {
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                v += (*it) * (*ix);
            }
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                v += M_SQRT1_2 * (*it) * ((*ix) * (*ix) - 1.0);
            }
            v += (*it++);
        }
        return v;
    }

    void mac(vector<double>& theta, const vector<double>& x, double alpha, int y) const {
        if (y == _nClasses - 1) return;
        vector<double>::iterator it = theta.begin() + y * (2 * _nFeatures + 1);
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            (*it) += alpha * (*ix);
        }
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            (*it) += M_SQRT1_2 * alpha * ((*ix) * (*ix) - 1.0);
        }
        (*it) += alpha;
    }

    void mac(vector<double>& theta, const vector<double>& x, const vector<double>& alpha) const {
        vector<double>::iterator it = theta.begin();
        for (int y = 0; y < _nClasses - 1; y++) {
#if 1
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                (*it) += alpha[y] * (*ix);
            }
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                (*it) += M_SQRT1_2 * alpha[y] * ((*ix) * (*ix) - 1.0);
            }
            (*it++) += alpha[y];
#else
            Eigen::Map<VectorXd>(&theta[y * (2 * _nFeatures + 1)], _nFeatures) += alpha[y] *
                Eigen::Map<const VectorXd>(&x[0], x.size());
            Eigen::Map<VectorXd>(&theta[y * (2 * _nFeatures + 1) + _nFeatures], _nFeatures) += M_SQRT1_2 * alpha[y] *
                (Eigen::Map<const VectorXd>(&x[0], x.size()).array().square().array() - 1.0);
            theta[y * (2 * _nFeatures + 1) + 2 * _nFeatures] += alpha[y];
#endif
        }
    }
};

// drwnQuadraticJointFeatureMap --------------------------------------------

class drwnQuadraticJointFeatureMap : public drwnJointFeatureMap {
 public:
    drwnQuadraticJointFeatureMap() : drwnJointFeatureMap() { /* do nothing */ }
    drwnQuadraticJointFeatureMap(int nFeatures, int nClasses) :
        drwnJointFeatureMap(nFeatures, nClasses) { /* do nothing */ }
    ~drwnQuadraticJointFeatureMap() { /* do nothing */ }

    // type and cloning
    const char *type () const { return "drwnQuadraticJointFeatureMap"; }
    drwnQuadraticJointFeatureMap *clone() const { return new drwnQuadraticJointFeatureMap(*this); }

    int numParameters() const { return std::max(((_nFeatures + 3) * _nFeatures / 2 + 1) *
            (_nClasses - 1), 0); }

    vector<double> operator()(const vector<double>& x, int y) const {
        vector<double> phi(numParameters(), 0.0);
        if (y != _nClasses - 1) {
            vector<double>::iterator it = phi.begin() + y * ((_nFeatures + 3) * _nFeatures / 2 + 1);
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                (*it) = (*ix);
            }
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                for (vector<double>::const_iterator jx = x.begin(); jx != ix; ++jx, ++it) {
                    (*it) = (*ix) * (*jx);
                }
                (*it) = M_SQRT1_2 * ((*ix) * (*ix) - 1.0);
            }
            (*it) = 1.0;
        }
        return phi;
    }

    vector<double> operator()(const vector<double>& x) const {
        vector<double> phi(numParameters());
        vector<double>::iterator it = phi.begin();
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            (*it) = (*ix);
        }
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            for (vector<double>::const_iterator jx = x.begin(); jx != ix; ++jx, ++it) {
                (*it) = (*ix) * (*jx);
            }
            (*it) = M_SQRT1_2 * ((*ix) * (*ix) - 1.0);
        }
        (*it) = 1.0;

        vector<double>::iterator jt = ++it;
        while (jt != phi.end()) {
            copy(phi.begin(), it, jt);
            jt += it - phi.begin();
        }
        return phi;
    }

    double dot(const vector<double>& theta, const vector<double>& x, int y) const {
        if (y == _nClasses - 1) return 0.0;
        double v = 0.0;
        vector<double>::const_iterator it = theta.begin() +
            y * ((_nFeatures + 3) * _nFeatures / 2 + 1);
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            v += (*it) * (*ix);
        }
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            for (vector<double>::const_iterator jx = x.begin(); jx != ix; ++jx, ++it) {
                v += (*it) * (*ix) * (*jx);
            }
            v += M_SQRT1_2 * (*it) * ((*ix) * (*ix) - 1.0);
        }
        v += (*it);

        return v;
    }

    double dot(const vector<double>& theta, const vector<double>& x) const {
        double v = 0.0;
        vector<double>::const_iterator it = theta.begin();
        for (int y = 0; y < _nClasses - 1; y++) {
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                v += (*it) * (*ix);
            }
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                for (vector<double>::const_iterator jx = x.begin(); jx != ix; ++jx, ++it) {
                    v += (*it) * (*ix) * (*jx);
                }
                v += M_SQRT1_2 * (*it) * ((*ix) * (*ix) - 1.0);
            }
            v += (*it++);
        }
        return v;
    }

    void mac(vector<double>& theta, const vector<double>& x, double alpha, int y) const {
        if (y == _nClasses - 1) return;
        vector<double>::iterator it = theta.begin() + y * ((_nFeatures + 3) * _nFeatures / 2 + 1);
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            (*it) += alpha * (*ix);
        }
        for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
            for (vector<double>::const_iterator jx = x.begin(); jx != ix; ++jx, ++it) {
                (*it) += alpha * (*ix) * (*jx);
            }
            (*it) += M_SQRT1_2 * alpha * ((*ix) * (*ix) - 1.0);
        }
        (*it) += alpha;
    }

    void mac(vector<double>& theta, const vector<double>& x, const vector<double>& alpha) const {
        vector<double>::iterator it = theta.begin();
        for (int y = 0; y < _nClasses - 1; y++) {
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                (*it) += alpha[y] * (*ix);
            }
            for (vector<double>::const_iterator ix = x.begin(); ix != x.end(); ++ix, ++it) {
                for (vector<double>::const_iterator jx = x.begin(); jx != ix; ++jx, ++it) {
                    (*it) += alpha[y] * (*ix) * (*jx);
                }
                (*it) += M_SQRT1_2 * alpha[y] * ((*ix) * (*ix) - 1.0);
            }
            (*it++) += alpha[y];
        }
    }
};