la.h

/* MLPACK 0.2
 *
 * Copyright (c) 2008, 2009 Alexander Gray,
 *                          Garry Boyer,
 *                          Ryan Riegel,
 *                          Nikolaos Vasiloglou,
 *                          Dongryeol Lee,
 *                          Chip Mappus, 
 *                          Nishant Mehta,
 *                          Hua Ouyang,
 *                          Parikshit Ram,
 *                          Long Tran,
 *                          Wee Chin Wong
 *
 * Copyright (c) 2008, 2009 Georgia Institute of Technology
 *
 * 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
 * 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., 51 Franklin Street, Fifth Floor, Boston, MA
 * 02110-1301, USA.
 */
#ifndef LINEAR_H
#define LINEAR_H

#include "fastlib/la/matrix.h"
#include "fastlib/la/uselapack.h"
//#include "matrix.h"
//#include "uselapack.h"

#include "fastlib/math/math_lib.h"

#include <math.h>

namespace la {
  inline double DistanceSqEuclidean(
      index_t length, const double *va, const double *vb) {
    double s = 0;
    do {
      double d = *va++ - *vb++;
      s += d * d;
    } while (--length);
    return s;
  }
  inline double DistanceSqEuclidean(const Vector& x, const Vector& y) {
    DEBUG_SAME_SIZE(x.length(), y.length());
    return DistanceSqEuclidean(x.length(), x.ptr(), y.ptr());
  }
  template<int t_pow>
  inline double RawLMetric(
      index_t length, const double *va, const double *vb) {
    double s = 0;
    do {
      double d = *va++ - *vb++;
      s += math::PowAbs<t_pow, 1>(d);
    } while (--length);
    return s;
  }
  template<int t_pow>
  inline double LMetric(
      index_t length, const double *va, const double *vb) {
    return math::Pow<1, t_pow>(RawLMetric<t_pow>(length, va, vb));
  }
  inline double Trace(Matrix &a) {
    // trace has meaning only for square matrices
    DEBUG_SAME_SIZE(a.n_cols(), a.n_rows());
    double trace=0;
    for(index_t i=0; i<a.n_cols(); i++) {
      trace+=a.get(i, i);
    }
    return trace;
  }
  
  inline success_t LeastSquareFit(Vector &y, Matrix &x, Vector *a) {
    DEBUG_SAME_SIZE(y.length(), x.n_rows());
    DEBUG_ASSERT(x.n_rows() >= x.n_cols());
    Vector r_xy_vec;
    Matrix r_xx_mat;
    la::MulTransAInit(x, x, &r_xx_mat);
    la::MulInit(y, x, &r_xy_vec);
    success_t status = la::SolveInit(r_xx_mat, r_xy_vec, a);
    if unlikely(status != SUCCESS_PASS) {
      if (status==SUCCESS_FAIL) {
        FATAL("Least square fit failed \n");
      } else {
        NONFATAL("Least square fit returned a warning \n");
      }
    }
    return status;
  }

  inline success_t LeastSquareFit(Matrix &y, Matrix &x, Matrix *a) {
    DEBUG_SAME_SIZE(y.n_rows(), x.n_rows());
    DEBUG_ASSERT(x.n_rows() >= x.n_cols());
    Matrix r_xy_mat;
    Matrix r_xx_mat;
    la::MulTransAInit(x, x, &r_xx_mat);
    la::MulTransAInit(x, y, &r_xy_mat);
    success_t status = la::SolveInit(r_xx_mat, r_xy_mat, a);
    if unlikely(status != SUCCESS_PASS) {
      if (status==SUCCESS_FAIL) {
        FATAL("Least square fit failed \n");
      } else {
        NONFATAL("Least square fit returned a warning \n");
      }
    }
    return status;
  }
  
  inline success_t LeastSquareFitTrans(Matrix &y, Matrix &x, Matrix *a) {
    DEBUG_SAME_SIZE(y.n_rows(), x.n_cols());
    DEBUG_ASSERT(x.n_cols() >= x.n_rows());
    Matrix r_xy_mat;
    Matrix r_xx_mat;
    la::MulTransBInit(x, x, &r_xx_mat);
    la::MulInit(x, y, &r_xy_mat);
    success_t status = la::SolveInit(r_xx_mat, r_xy_mat, a);
    if unlikely(status != SUCCESS_PASS) {
      if (status==SUCCESS_FAIL) {
        FATAL("Least square fit failed \n");
      } else {
        NONFATAL("Least square fit returned a warning \n");
      }
    }
    return status;
  }
};

#endif