blob: 8e148f9bc71bff99133e586ceca1416ae8a88c74 [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2012 Desire Nuentsa // Copyright (C) 2012 Gael Guennebaud // // 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 #include #include #include "main.h" #include using namespace std; using namespace Eigen; template struct sparseGaussianTest : SparseFunctor { typedef Matrix VectorType; typedef SparseFunctor Base; typedef typename Base::JacobianType JacobianType; sparseGaussianTest(int inputs, int values) : SparseFunctor(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 u(uv, 0, half); VectorBlock 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 u(uv, 0, half); VectorBlock 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 u(uv, 0, half); VectorBlock 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 void test_sparseLM_T() { typedef Matrix VectorType; int inputs = 10; int values = 2000; sparseGaussianTest 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 u(uv, 0, inputs/2); VectorBlock v(uv, inputs/2, inputs/2); v.setOnes(); //Generate u or Solve for u from v u.setOnes(); // Solve the optimization problem LevenbergMarquardt > 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()); // CALL_SUBTEST_2(test_sparseLM_T()); }