1c1f9011ceba380967cde2019b352f963f90f9589reed@google.com// This file is part of Eigen, a lightweight C++ template library 2c1f9011ceba380967cde2019b352f963f90f9589reed@google.com// for linear algebra. 3c1f9011ceba380967cde2019b352f963f90f9589reed@google.com// 4c1f9011ceba380967cde2019b352f963f90f9589reed@google.com// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> 5c1f9011ceba380967cde2019b352f963f90f9589reed@google.com// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> 6c1f9011ceba380967cde2019b352f963f90f9589reed@google.com// 7c1f9011ceba380967cde2019b352f963f90f9589reed@google.com// This Source Code Form is subject to the terms of the Mozilla 8f168b86d7fafc5c20c87bebc6fd393cb17e120catfarina// Public License v. 2.0. If a copy of the MPL was not distributed 9c1f9011ceba380967cde2019b352f963f90f9589reed@google.com// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10c1f9011ceba380967cde2019b352f963f90f9589reed@google.com#include <iostream> 11c1f9011ceba380967cde2019b352f963f90f9589reed@google.com#include <fstream> 12c1f9011ceba380967cde2019b352f963f90f9589reed@google.com#include <iomanip> 13c1f9011ceba380967cde2019b352f963f90f9589reed@google.com 14c1f9011ceba380967cde2019b352f963f90f9589reed@google.com#include "main.h" 15c1f9011ceba380967cde2019b352f963f90f9589reed@google.com#include <Eigen/LevenbergMarquardt> 16c1f9011ceba380967cde2019b352f963f90f9589reed@google.com 17c1f9011ceba380967cde2019b352f963f90f9589reed@google.comusing namespace std; 18c1f9011ceba380967cde2019b352f963f90f9589reed@google.comusing namespace Eigen; 19c1f9011ceba380967cde2019b352f963f90f9589reed@google.com 20c1f9011ceba380967cde2019b352f963f90f9589reed@google.comtemplate <typename Scalar> 21c1f9011ceba380967cde2019b352f963f90f9589reed@google.comstruct sparseGaussianTest : SparseFunctor<Scalar, int> 22c1f9011ceba380967cde2019b352f963f90f9589reed@google.com{ 23c1f9011ceba380967cde2019b352f963f90f9589reed@google.com typedef Matrix<Scalar,Dynamic,1> VectorType; 24c1f9011ceba380967cde2019b352f963f90f9589reed@google.com typedef SparseFunctor<Scalar,int> Base; 25c1f9011ceba380967cde2019b352f963f90f9589reed@google.com typedef typename Base::JacobianType JacobianType; 26c1f9011ceba380967cde2019b352f963f90f9589reed@google.com sparseGaussianTest(int inputs, int values) : SparseFunctor<Scalar,int>(inputs,values) 27c1f9011ceba380967cde2019b352f963f90f9589reed@google.com { } 28c1f9011ceba380967cde2019b352f963f90f9589reed@google.com 29c1f9011ceba380967cde2019b352f963f90f9589reed@google.com VectorType model(const VectorType& uv, VectorType& x) 30c1f9011ceba380967cde2019b352f963f90f9589reed@google.com { 31c1f9011ceba380967cde2019b352f963f90f9589reed@google.com VectorType y; //Change this to use expression template 32c1f9011ceba380967cde2019b352f963f90f9589reed@google.com int m = Base::values(); 33c1f9011ceba380967cde2019b352f963f90f9589reed@google.com int n = Base::inputs(); 34c1f9011ceba380967cde2019b352f963f90f9589reed@google.com eigen_assert(uv.size()%2 == 0); 35c1f9011ceba380967cde2019b352f963f90f9589reed@google.com eigen_assert(uv.size() == n); 36c1f9011ceba380967cde2019b352f963f90f9589reed@google.com eigen_assert(x.size() == m); 37c1f9011ceba380967cde2019b352f963f90f9589reed@google.com y.setZero(m); 38c1f9011ceba380967cde2019b352f963f90f9589reed@google.com int half = n/2; 39c1f9011ceba380967cde2019b352f963f90f9589reed@google.com VectorBlock<const VectorType> u(uv, 0, half); 40c1f9011ceba380967cde2019b352f963f90f9589reed@google.com VectorBlock<const VectorType> v(uv, half, half); 41c1f9011ceba380967cde2019b352f963f90f9589reed@google.com Scalar coeff; 42c1f9011ceba380967cde2019b352f963f90f9589reed@google.com for (int j = 0; j < m; j++) 43c1f9011ceba380967cde2019b352f963f90f9589reed@google.com { 44c1f9011ceba380967cde2019b352f963f90f9589reed@google.com for (int i = 0; i < half; i++) 45c1f9011ceba380967cde2019b352f963f90f9589reed@google.com { 46c1f9011ceba380967cde2019b352f963f90f9589reed@google.com coeff = (x(j)-i)/v(i); 47c1f9011ceba380967cde2019b352f963f90f9589reed@google.com coeff *= coeff; 48c1f9011ceba380967cde2019b352f963f90f9589reed@google.com if (coeff < 1. && coeff > 0.) 49c1f9011ceba380967cde2019b352f963f90f9589reed@google.com y(j) += u(i)*std::pow((1-coeff), 2); 50c1f9011ceba380967cde2019b352f963f90f9589reed@google.com } 51c1f9011ceba380967cde2019b352f963f90f9589reed@google.com } 52c1f9011ceba380967cde2019b352f963f90f9589reed@google.com return y; 53c1f9011ceba380967cde2019b352f963f90f9589reed@google.com } 54c1f9011ceba380967cde2019b352f963f90f9589reed@google.com void initPoints(VectorType& uv_ref, VectorType& x) 55c1f9011ceba380967cde2019b352f963f90f9589reed@google.com { 56f168b86d7fafc5c20c87bebc6fd393cb17e120catfarina m_x = x; 57c1f9011ceba380967cde2019b352f963f90f9589reed@google.com m_y = this->model(uv_ref,x); 58c1f9011ceba380967cde2019b352f963f90f9589reed@google.com } 59c1f9011ceba380967cde2019b352f963f90f9589reed@google.com int operator()(const VectorType& uv, VectorType& fvec) 60c1f9011ceba380967cde2019b352f963f90f9589reed@google.com { 61c1f9011ceba380967cde2019b352f963f90f9589reed@google.com int m = Base::values(); 62410e6e80f00a6c660675c80904807a041c7b7d2amtklein@google.com int n = Base::inputs(); 63c1f9011ceba380967cde2019b352f963f90f9589reed@google.com eigen_assert(uv.size()%2 == 0); 64c1f9011ceba380967cde2019b352f963f90f9589reed@google.com eigen_assert(uv.size() == n); 65c1f9011ceba380967cde2019b352f963f90f9589reed@google.com int half = n/2; 66c1f9011ceba380967cde2019b352f963f90f9589reed@google.com VectorBlock<const VectorType> u(uv, 0, half); 67c1f9011ceba380967cde2019b352f963f90f9589reed@google.com VectorBlock<const VectorType> v(uv, half, half); 68c1f9011ceba380967cde2019b352f963f90f9589reed@google.com fvec = m_y; 69c1f9011ceba380967cde2019b352f963f90f9589reed@google.com Scalar coeff; 700c38ed3b1d704a0ed6147299046f51fd52e841a2skia.committer@gmail.com for (int j = 0; j < m; j++) 71c1f9011ceba380967cde2019b352f963f90f9589reed@google.com { 72c1f9011ceba380967cde2019b352f963f90f9589reed@google.com for (int i = 0; i < half; i++) 73c1f9011ceba380967cde2019b352f963f90f9589reed@google.com { 74c1f9011ceba380967cde2019b352f963f90f9589reed@google.com coeff = (m_x(j)-i)/v(i); 75c1f9011ceba380967cde2019b352f963f90f9589reed@google.com coeff *= coeff; 76c1f9011ceba380967cde2019b352f963f90f9589reed@google.com if (coeff < 1. && coeff > 0.) 77c1f9011ceba380967cde2019b352f963f90f9589reed@google.com fvec(j) -= u(i)*std::pow((1-coeff), 2); 78c1f9011ceba380967cde2019b352f963f90f9589reed@google.com } 79c1f9011ceba380967cde2019b352f963f90f9589reed@google.com } 80c1f9011ceba380967cde2019b352f963f90f9589reed@google.com return 0; 813361471a3504ecd0351ff70f4c42d8d6fee963d4commit-bot@chromium.org } 82c1f9011ceba380967cde2019b352f963f90f9589reed@google.com 83c1f9011ceba380967cde2019b352f963f90f9589reed@google.com int df(const VectorType& uv, JacobianType& fjac) 84c1f9011ceba380967cde2019b352f963f90f9589reed@google.com { 85c1f9011ceba380967cde2019b352f963f90f9589reed@google.com int m = Base::values(); 86c1f9011ceba380967cde2019b352f963f90f9589reed@google.com int n = Base::inputs(); 873361471a3504ecd0351ff70f4c42d8d6fee963d4commit-bot@chromium.org eigen_assert(n == uv.size()); 88c1f9011ceba380967cde2019b352f963f90f9589reed@google.com eigen_assert(fjac.rows() == m); 89c1f9011ceba380967cde2019b352f963f90f9589reed@google.com eigen_assert(fjac.cols() == n); 90c1f9011ceba380967cde2019b352f963f90f9589reed@google.com int half = n/2; 91f168b86d7fafc5c20c87bebc6fd393cb17e120catfarina VectorBlock<const VectorType> u(uv, 0, half); 92c1f9011ceba380967cde2019b352f963f90f9589reed@google.com VectorBlock<const VectorType> v(uv, half, half); 93c1f9011ceba380967cde2019b352f963f90f9589reed@google.com Scalar coeff; 94410e6e80f00a6c660675c80904807a041c7b7d2amtklein@google.com 95410e6e80f00a6c660675c80904807a041c7b7d2amtklein@google.com //Derivatives with respect to u 96410e6e80f00a6c660675c80904807a041c7b7d2amtklein@google.com for (int col = 0; col < half; col++) 97c1f9011ceba380967cde2019b352f963f90f9589reed@google.com { 98410e6e80f00a6c660675c80904807a041c7b7d2amtklein@google.com for (int row = 0; row < m; row++) 99410e6e80f00a6c660675c80904807a041c7b7d2amtklein@google.com { 100410e6e80f00a6c660675c80904807a041c7b7d2amtklein@google.com coeff = (m_x(row)-col)/v(col); 101c1f9011ceba380967cde2019b352f963f90f9589reed@google.com coeff = coeff*coeff; 102410e6e80f00a6c660675c80904807a041c7b7d2amtklein@google.com if(coeff < 1. && coeff > 0.) 103410e6e80f00a6c660675c80904807a041c7b7d2amtklein@google.com { 104410e6e80f00a6c660675c80904807a041c7b7d2amtklein@google.com fjac.coeffRef(row,col) = -(1-coeff)*(1-coeff); 105 } 106 } 107 } 108 //Derivatives with respect to v 109 for (int col = 0; col < half; col++) 110 { 111 for (int row = 0; row < m; row++) 112 { 113 coeff = (m_x(row)-col)/v(col); 114 coeff = coeff*coeff; 115 if(coeff < 1. && coeff > 0.) 116 { 117 fjac.coeffRef(row,col+half) = -4 * (u(col)/v(col))*coeff*(1-coeff); 118 } 119 } 120 } 121 return 0; 122 } 123 124 VectorType m_x, m_y; //Data points 125}; 126 127 128template<typename T> 129void test_sparseLM_T() 130{ 131 typedef Matrix<T,Dynamic,1> VectorType; 132 133 int inputs = 10; 134 int values = 2000; 135 sparseGaussianTest<T> sparse_gaussian(inputs, values); 136 VectorType uv(inputs),uv_ref(inputs); 137 VectorType x(values); 138 // Generate the reference solution 139 uv_ref << -2, 1, 4 ,8, 6, 1.8, 1.2, 1.1, 1.9 , 3; 140 //Generate the reference data points 141 x.setRandom(); 142 x = 10*x; 143 x.array() += 10; 144 sparse_gaussian.initPoints(uv_ref, x); 145 146 147 // Generate the initial parameters 148 VectorBlock<VectorType> u(uv, 0, inputs/2); 149 VectorBlock<VectorType> v(uv, inputs/2, inputs/2); 150 v.setOnes(); 151 //Generate u or Solve for u from v 152 u.setOnes(); 153 154 // Solve the optimization problem 155 LevenbergMarquardt<sparseGaussianTest<T> > lm(sparse_gaussian); 156 int info; 157// info = lm.minimize(uv); 158 159 VERIFY_IS_EQUAL(info,1); 160 // Do a step by step solution and save the residual 161 int maxiter = 200; 162 int iter = 0; 163 MatrixXd Err(values, maxiter); 164 MatrixXd Mod(values, maxiter); 165 LevenbergMarquardtSpace::Status status; 166 status = lm.minimizeInit(uv); 167 if (status==LevenbergMarquardtSpace::ImproperInputParameters) 168 return ; 169 170} 171void test_sparseLM() 172{ 173 CALL_SUBTEST_1(test_sparseLM_T<double>()); 174 175 // CALL_SUBTEST_2(test_sparseLM_T<std::complex<double>()); 176} 177