eigen/test/sparseLM.cpp

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include <iostream>
#include <fstream>
#include <iomanip>

#include "main.h"
#include <Eigen/LevenbergMarquardt>

using namespace std;
using namespace Eigen;

template <typename Scalar>
struct sparseGaussianTest : SparseFunctor<Scalar, int>
{
  typedef Matrix<Scalar,Dynamic,1> VectorType;
  typedef SparseFunctor<Scalar,int> Base;
  typedef typename Base::JacobianType JacobianType;
  sparseGaussianTest(int inputs, int values) : SparseFunctor<Scalar,int>(inputs,values)
  { }
  
  VectorType model(const VectorType& uv, VectorType& x)
  {
    VectorType y; //Change this to use expression template
    int m = Base::values(); 
    int n = Base::inputs();
    eigen_assert(uv.size()%2 == 0);
    eigen_assert(uv.size() == n);
    eigen_assert(x.size() == m);
    y.setZero(m);
    int half = n/2;
    VectorBlock<const VectorType> u(uv, 0, half);
    VectorBlock<const VectorType> v(uv, half, half);
    Scalar coeff;
    for (int j = 0; j < m; j++)
    {
      for (int i = 0; i < half; i++) 
      {
        coeff = (x(j)-i)/v(i);
        coeff *= coeff;
        if (coeff < 1. && coeff > 0.)
          y(j) += u(i)*std::pow((1-coeff), 2);
      }
    }
    return y;
  }
  void initPoints(VectorType& uv_ref, VectorType& x)
  {
    m_x = x;
    m_y = this->model(uv_ref,x);
  }
  int operator()(const VectorType& uv, VectorType& fvec)
  {
    int m = Base::values(); 
    int n = Base::inputs();
    eigen_assert(uv.size()%2 == 0);
    eigen_assert(uv.size() == n);
    int half = n/2;
    VectorBlock<const VectorType> u(uv, 0, half);
    VectorBlock<const VectorType> v(uv, half, half);
    fvec = m_y;
    Scalar coeff;
    for (int j = 0; j < m; j++)
    {
      for (int i = 0; i < half; i++)
      {
        coeff = (m_x(j)-i)/v(i);
        coeff *= coeff;
        if (coeff < 1. && coeff > 0.)
          fvec(j) -= u(i)*std::pow((1-coeff), 2);
      }
    }
    return 0;
  }
  
  int df(const VectorType& uv, JacobianType& fjac)
  {
    int m = Base::values(); 
    int n = Base::inputs();
    eigen_assert(n == uv.size());
    eigen_assert(fjac.rows() == m);
    eigen_assert(fjac.cols() == n);
    int half = n/2;
    VectorBlock<const VectorType> u(uv, 0, half);
    VectorBlock<const VectorType> v(uv, half, half);
    Scalar coeff;
    
    //Derivatives with respect to u
    for (int col = 0; col < half; col++)
    {
      for (int row = 0; row < m; row++)
      {
        coeff = (m_x(row)-col)/v(col);
          coeff = coeff*coeff;
        if(coeff < 1. && coeff > 0.)
        {
          fjac.coeffRef(row,col) = -(1-coeff)*(1-coeff);
        }
      }
    }
    //Derivatives with respect to v
    for (int col = 0; col < half; col++)
    {
      for (int row = 0; row < m; row++)
      {
        coeff = (m_x(row)-col)/v(col);
        coeff = coeff*coeff;
        if(coeff < 1. && coeff > 0.)
        {
          fjac.coeffRef(row,col+half) = -4 * (u(col)/v(col))*coeff*(1-coeff);
        }
      }
    }
    return 0;
  }
  
  VectorType m_x, m_y; //Data points
};


template<typename T>
void test_sparseLM_T()
{
  typedef Matrix<T,Dynamic,1> VectorType;
  
  int inputs = 10;
  int values = 2000;
  sparseGaussianTest<T> sparse_gaussian(inputs, values);
  VectorType uv(inputs),uv_ref(inputs);
  VectorType x(values);
  // Generate the reference solution 
  uv_ref << -2, 1, 4 ,8, 6, 1.8, 1.2, 1.1, 1.9 , 3;
  //Generate the reference data points
  x.setRandom();
  x = 10*x;
  x.array() += 10;
  sparse_gaussian.initPoints(uv_ref, x);
  
  
  // Generate the initial parameters 
  VectorBlock<VectorType> u(uv, 0, inputs/2); 
  VectorBlock<VectorType> v(uv, inputs/2, inputs/2);
  v.setOnes();
  //Generate u or Solve for u from v
  u.setOnes();
  
  // Solve the optimization problem
  LevenbergMarquardt<sparseGaussianTest<T> > lm(sparse_gaussian);
  int info;
//   info = lm.minimize(uv);
  
  VERIFY_IS_EQUAL(info,1);
    // Do a step by step solution and save the residual 
  int maxiter = 200;
  int iter = 0;
  MatrixXd Err(values, maxiter);
  MatrixXd Mod(values, maxiter);
  LevenbergMarquardtSpace::Status status; 
  status = lm.minimizeInit(uv);
  if (status==LevenbergMarquardtSpace::ImproperInputParameters)
      return ;

}
void test_sparseLM()
{
  CALL_SUBTEST_1(test_sparseLM_T<double>());
  
  // CALL_SUBTEST_2(test_sparseLM_T<std::complex<double>());
}