STK_LinearAAModel.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2016  Serge Iovleff
 
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser 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 Lesser General Public License for more details.
 
    You should have received a copy of the GNU Lesser 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
 
    Contact : S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)                                   */
 
/*
 * Project:  stkpp::AAModels
 * Purpose:  class for AA Linear models.
 * Author:   Serge Iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 **/
 
#ifndef STK_LINEARAAMODEL_H
#define STK_LINEARAAMODEL_H
 
#include "STK_GaussianAAModel.h"
#include <Reduct/include/STK_ILinearReduct.h>
#include <Regress/include/STK_MultidimRegression.h>
#include <Algebra/include/STK_GramSchmidt.h>
#include <STatistiK/include/STK_Law_Normal.h>
 
#ifdef STK_AAMODELS_VERBOSE
#include <Arrays/include/STK_Display.h>
#endif
 
 
namespace STK
{
 
namespace Law
{
template<class TYPE>
class IUnivLaw;
}
 
template<class Array>
class LinearAAModel: public IRunnerUnsupervised<Array, VectorX>
                    , public GaussianAAModel<Array>
{
  public :
    typedef IRunnerUnsupervised<Array, VectorX> Runner;
    typedef GaussianAAModel<Array> Base;
    using Base::p_workData_;
    using Base::p_regressor_;
    using Base::p_reducer_;
    using Base::p_reduced_;
    using Base::p_predicted_;
    using Base::p_residuals_;
    using Base::dim;
    using Base::regressionStep;
    using Base::reductionStep;
    LinearAAModel( Array const& data);
    LinearAAModel( LinearAAModel const& model);
    virtual ~LinearAAModel();
    inline virtual LinearAAModel* clone() const { return new LinearAAModel(*this);}
    bool run( int const& dim);
    bool run( VectorX const& weights, int const& dim);
    virtual bool run();
    virtual bool run(VectorX const& weights);
    static void simul( const Law::IUnivLaw<Real>& law
                     , VectorX const& mu
                     , Real const& std
                     , Array& proj
                     , Array& data
                     );
    protected:
    Array workData_;
};
 
/* Constructors.                                                             */
/* Default constructor : compute the Linear AA models of the matrix X
 *  using the local variance as criteria.
**/
template<class Array>
LinearAAModel<Array>::LinearAAModel( Array const& data)
                                   : Runner(), Base(0), workData_(data)
{
  Runner::setData(workData_);
  Base::setWorkData(workData_);
  p_regressor_ = new MultidimRegression<Array, VectorX>();
}
 
template<class Array>
LinearAAModel<Array>::LinearAAModel( LinearAAModel const& model)
                                   : Runner(), Base(0), workData_(model.workData_)
{
  Runner::setData(workData_);
  Base::setWorkData(workData_);
  p_regressor_ = new MultidimRegression<Array, VectorX>();
}
 
 
 
/* Virtual destructor */
template<class Array>
LinearAAModel<Array>::~LinearAAModel()
{ if (p_regressor_) delete p_regressor_;}
 
/* run the probabilistic AAModel */
template<class Array>
bool LinearAAModel<Array>::run( int const& dim)
{
  this->setDimension(dim);
  bool res = run();
#ifdef STK_DEBUG
  if (!res)
  { stk_cerr << _T("An error occur.\nWhat: ") << this->error();}
#endif
  return res;
}
 
 
/* run the probabilistic AAModel */
template<class Array>
bool LinearAAModel<Array>::run()
{
  // compute AAM
  try
  {
    if (!p_reducer_)
    {
      this->msg_error_ = STKERROR_NO_ARG(LinearAAModel::run(),reducer is not set);
      return false;
    }
    if (!p_regressor_)
    {
      this->msg_error_ = STKERROR_NO_ARG(LinearAAModel::run(),regressor is not be set);
      return false;
    }
    // set p_workData to the reducer
    p_reducer_->setData(*p_workData_);
    // compute the projected data set
    reductionStep();
    // compute the projected covariance
    this->computeProjectedCovariance();
#ifdef STK_AAMODELS_VERBOSE
    stk_cout << _T("In LinearAAModel::run(), reduction done.\n");
#endif
    // set data
    p_regressor_->setY(*p_workData_);
    p_regressor_->setX(*p_reduced_);
    // compute the regression function
    regressionStep();
#ifdef STK_AAMODELS_VERBOSE
    stk_cout << _T("In LinearAAModel::run(), regression done.\n");
#endif
    this->computeModelParameters();
    // check if data have been standardized or centered
    if (this->isStandardized()) { this->unstandardizeResults();}
    else if (this->isCentered()){ this->uncenterResults();}
  }
  catch (Exception const& error)
  {
    this->msg_error_ = _T("Error in LinearAAModel::run():\nWhat: ");
    this->msg_error_ += error.error();
    return false;
  }
  return true;
}
 
/* run the probabilistic AAModel */
template<class Array>
bool LinearAAModel<Array>::run( VectorX const& weights, int const& dim)
{
  this->setDimension(dim);
  return run(weights);
}
template<class Array>
bool LinearAAModel<Array>::run( VectorX const& weights)
{
  try
  {
    if (!p_reducer_) STKRUNTIME_ERROR_NO_ARG(LinearAAModel::run,reducer is not set);
    if (!p_regressor_) STKRUNTIME_ERROR_NO_ARG(LinearAAModel::run,regressor is not set);
    // set p_workData to the reducer
    p_reducer_->setData(*p_workData_);
    // compute the weighted reduced data set
    reductionStep(weights);
    // compute the projected covariance
    this->computeProjectedCovariance();
#ifdef STK_AAMODELS_VERBOSE
    stk_cout << _T("In LinearAAModel::run(weights), reduction done.\n");
#endif
    // set data
    p_regressor_->setY(*p_workData_);
    p_regressor_->setX(*p_reduced_);
    // compute the weighted regression vectors
    regressionStep(weights);
#ifdef STK_AAMODELS_VERBOSE
    stk_cout << _T("In LinearAAModel::run(weights), regression done.\n");
#endif
    this->computeModelParameters();
    // check if data have been standardized or centered
    if (this->isStandardized()) { this->unstandardizeResults();}
    else if (this->isCentered()){ this->uncenterResults();}
  }
  catch (Exception const& error)
  {
    this->msg_error_ = _T("Error in LinearAAModel::run(weights): ");
    this->msg_error_ += error.error();
    return false;
  }
  return true;
}
 
/* Simulate a centered auto-associative linear model of the form
 * \f[
 *    X = X.P.P' + \epsilon
 * \f]
 * with \f$ P'P = I_d \f$ and d < p.
 * @param law, the law to use in order to simulate the data.
 * @param std the standard deviation of the gaussian noise
 * @param data the data to simulate. The dimension of the container
 * give the number of the samples and variables.
 * @param proj the simulated projection matrix. The dimension of the
 * container give the dimension of the AA model.
 **/
template<class Array>
void LinearAAModel<Array>::simul( const Law::IUnivLaw<Real>& law
                         , VectorX const& mu, Real const& std
                         , Array& proj, Array& data
                         )
{
  if (data.cols() != mu.range())
    STKRUNTIME_ERROR_NO_ARG(LinearAAModel::simul,data.cols()!=mu.range());
  // simul AA model
  data.rand(law);
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("data simulated\n");
#endif
  // choose arbitrary coefficients for proj
  proj.randGauss();
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("proj simulated\n");
#endif
  // orthonormalize proj
  gramSchmidt(proj);
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("proj orthonormalized\n");
#endif
 
#ifdef STK_AAMODELS_DEBUG
  stk_cout << _T("proj =") << proj;
#endif
  // compute data
  data = data * (proj * proj.transpose());
 
  // add noise to the model
  for (int j= data.beginCols(); j< data.endCols(); j++)
  {
    Law::Normal noise(mu[j], std);
    for (int i= data.beginRows(); i< data.endRows(); i++)
      data(i, j) += noise.rand();
  }
}
 
} // namespace STK
 
#endif // STK_LINEARAAMODEL_H