From b4c5237925744ab0639425142cfbc65d93118ff6 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Mon, 15 Dec 2014 15:08:41 -0500 Subject: [PATCH 001/130] Added structure of test and class. Left with implementation. --- ...tLinearlyConstrainedNonlinearOptimizer.cpp | 137 ++++++++++++++++++ 1 file changed, 137 insertions(+) create mode 100644 gtsam_unstable/nonlinear/tests/testLinearlyConstrainedNonlinearOptimizer.cpp diff --git a/gtsam_unstable/nonlinear/tests/testLinearlyConstrainedNonlinearOptimizer.cpp b/gtsam_unstable/nonlinear/tests/testLinearlyConstrainedNonlinearOptimizer.cpp new file mode 100644 index 000000000..78b511e2c --- /dev/null +++ b/gtsam_unstable/nonlinear/tests/testLinearlyConstrainedNonlinearOptimizer.cpp @@ -0,0 +1,137 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ + +/** + * @file testLinearlyConstrainedNonlinearOptimizer.cpp + * @brief Unit tests for LinearlyConstrainedNonlinearOptimizer + * @author Krunal Chande + * @author Duy-Nguyen Ta + * @author Luca Carlone + * @date Dec 15, 2014 + */ + +#include +#include +#include +#include +#include +#include + + +namespace gtsam { +struct LinearlyConstrainedNLP { + NonlinearFactorGraph cost; + LinearEqualityFactorGraph equalities; + LinearInequalityFactorGraph inequalities; +}; + +struct LinearlyConstrainedNLPState { + Values values; + VectorValues duals; + bool converged; + LinearlyConstrainedNLPState(const Values& initialValues) : + values(initialValues), duals(VectorValues()), converged(false) { + } +}; +class LinearlyConstrainedNonLinearOptimizer { + LinearlyConstrainedNLP lcNLP_; +public: + LinearlyConstrainedNonLinearOptimizer(const LinearlyConstrainedNLP& lcNLP): lcNLP_(lcNLP) {} + + LinearlyConstrainedNLPState iterate(const LinearlyConstrainedNLPState& state) const { + QP qp; + qp.cost = lcNLP_.cost.linearize(state.values); + qp.equalities = lcNLP_.equalities; + qp.inequalities = lcNLP_.inequalities; + QPSolver qpSolver(qp); + VectorValues delta, duals; + boost::tie(delta, duals) = qpSolver.optimize(); + LinearlyConstrainedNLPState newState; + newState.values = state.values.retract(delta); + newState.duals = duals; + newState.converged = checkConvergence(newState.values, newState.duals); + return newState; + } + + /** + * Main optimization function. + */ + std::pair optimize(const Values& initialValues) const { + LinearlyConstrainedNLPState state(initialValues); + while(!state.converged){ + state = iterate(state); + } + + return std::make_pair(initialValues, VectorValues()); + } +}; +} + +using namespace std; +using namespace gtsam::symbol_shorthand; +using namespace gtsam; +const double tol = 1e-10; +//****************************************************************************** +TEST(LinearlyConstrainedNonlinearOptimizer, Problem1 ) { + + // build a quadratic Objective function x1^2 - x1*x2 + x2^2 - 3*x1 + 5 + // Note the Hessian encodes: + // 0.5*x1'*G11*x1 + x1'*G12*x2 + 0.5*x2'*G22*x2 - x1'*g1 - x2'*g2 + 0.5*f + // Hence, we have G11=2, G12 = -1, g1 = +3, G22 = 2, g2 = 0, f = 10 + HessianFactor lf(X(1), X(2), 2.0 * ones(1, 1), -ones(1, 1), 3.0 * ones(1), + 2.0 * ones(1, 1), zero(1), 10.0); + + // build linear inequalities + LinearInequalityFactorGraph inequalities; + inequalities.push_back(LinearInequality(X(1), ones(1,1), X(2), ones(1,1), 2, 0)); // x1 + x2 <= 2 --> x1 + x2 -2 <= 0, --> b=2 + inequalities.push_back(LinearInequality(X(1), -ones(1,1), 0, 1)); // -x1 <= 0 + inequalities.push_back(LinearInequality(X(2), -ones(1,1), 0, 2)); // -x2 <= 0 + inequalities.push_back(LinearInequality(X(1), ones(1,1), 1.5, 3)); // x1 <= 3/2 + + // Instantiate LinearlyConstrainedNLP, pretending that the cost is not quadratic + // (LinearContainerFactor makes a linear factor behave like a nonlinear one) + LinearlyConstrainedNLP lcNLP; + lcNLP.cost.add(LinearContainerFactor(lf)); + lcNLP.inequalities = inequalities; + + // Compare against a QP + QP qp; + qp.cost.add(lf); + qp.inequalities = inequalities; + + // instantiate QPsolver + QPSolver qpSolver(qp); + // create initial values for optimization + VectorValues initialVectorValues; + initialVectorValues.insert(X(1), zero(1)); + initialVectorValues.insert(X(2), zero(1)); + VectorValues expectedSolution = qpSolver.optimize(initialVectorValues).first; + + // instantiate LinearlyConstrainedNonLinearOptimizer + LinearlyConstrainedNonLinearOptimizer lcNLPSolver(lcNLP); + // create initial values for optimization + Values initialValues; + initialValues.insert(X(1), 0.0); + initialValues.insert(X(2), 0.0); + Values actualSolution = lcNLPSolver.optimize(initialValues).first; + + + DOUBLES_EQUAL(expectedSolution.at(X(1))[0], actualSolution.at(X(1)), tol); + DOUBLES_EQUAL(expectedSolution.at(X(2))[0], actualSolution.at(X(2)), tol); +} + +//****************************************************************************** +int main() { + TestResult tr; + return TestRegistry::runAllTests(tr); +} +//****************************************************************************** + From 7d2f69335d366850e5a1955881914bb995b70fa0 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Thu, 18 Dec 2014 20:47:14 -0500 Subject: [PATCH 002/130] Added constructor --- gtsam_unstable/linear/LinearEquality.h | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/gtsam_unstable/linear/LinearEquality.h b/gtsam_unstable/linear/LinearEquality.h index bc1b2bc12..5d96d8d38 100644 --- a/gtsam_unstable/linear/LinearEquality.h +++ b/gtsam_unstable/linear/LinearEquality.h @@ -41,6 +41,16 @@ public: Base() { } + /** + * Construct from a constrained noisemodel JacobianFactor with a dual key. + */ + explicit LinearEquality(const JacobianFactor& jf, Key dualKey) : Base(jf), dualKey_(dualKey){ + if (!jf.isConstrained()) { + throw std::runtime_error("Cannot convert an unconstrained JacobianFactor to LinearEquality"); + } + } + + /** Conversion from HessianFactor (does Cholesky to obtain Jacobian matrix) */ explicit LinearEquality(const HessianFactor& hf) { throw std::runtime_error("Cannot convert HessianFactor to LinearEquality"); From 3af06ef029fb1655a80d1758ee8944dd5e6af1d5 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Thu, 18 Dec 2014 20:48:10 -0500 Subject: [PATCH 003/130] Added nonlinearconstraint and tests --- .../nonlinear/NonlinearConstraint.h | 544 ++++++++++++++++++ .../tests/testNonlinearConstraints.cpp | 176 ++++++ 2 files changed, 720 insertions(+) create mode 100644 gtsam_unstable/nonlinear/NonlinearConstraint.h create mode 100644 gtsam_unstable/nonlinear/tests/testNonlinearConstraints.cpp diff --git a/gtsam_unstable/nonlinear/NonlinearConstraint.h b/gtsam_unstable/nonlinear/NonlinearConstraint.h new file mode 100644 index 000000000..3d1bf7eb0 --- /dev/null +++ b/gtsam_unstable/nonlinear/NonlinearConstraint.h @@ -0,0 +1,544 @@ +/** + * @file NonlinearConstraint.h + * @brief + * @author Duy-Nguyen Ta + * @date Sep 30, 2013 + */ + +#pragma once + +#include +#include +#include +#include +#include +//#include "DualKeyGenerator.h" + +namespace gtsam { + +class NonlinearConstraint { + Key dualKey_; + +public: + typedef boost::shared_ptr shared_ptr; +public: + /// Construct with dual key + NonlinearConstraint(Key dualKey) : dualKey_(dualKey) {} + + /** + * compute the HessianFactor of the (-dual * constraintHessian) for the qp subproblem's objective function + */ + virtual GaussianFactor::shared_ptr multipliedHessian(const Values& x, + const VectorValues& duals) const = 0; + + /// return the dual key + Key dualKey() const { return dualKey_; } +}; + +/* ************************************************************************* */ +/** A convenient base class for creating a nonlinear equality constraint with 1 + * variables. To derive from this class, implement evaluateError(). */ +template +class NonlinearConstraint1: public NoiseModelFactor1, public NonlinearConstraint { + +public: + + // typedefs for value types pulled from keys + typedef VALUE X; + +protected: + + typedef NoiseModelFactor1 Base; + typedef NonlinearConstraint1 This; + +public: + + /** + * Default Constructor for I/O + */ + NonlinearConstraint1() { + } + + /** + * Constructor + * @param j key of the variable + * @param constraintDim number of dimensions of the constraint error function + */ + NonlinearConstraint1(Key key, Key dualKey, size_t constraintDim = 1) : + Base(noiseModel::Constrained::All(constraintDim), key), NonlinearConstraint(dualKey) { + } + + virtual ~NonlinearConstraint1() { + } + + /** + * Override this method to finish implementing a binary factor. + * If any of the optional Matrix reference arguments are specified, it should compute + * both the function evaluation and its derivative(s) in X1 (and/or X2). + */ + virtual Vector + evaluateError(const X&, boost::optional H1 = boost::none) const = 0; + + /** produce a Gaussian factor graph containing all Hessian Factors, scaled by lambda, + * corresponding to this constraint */ + virtual GaussianFactor::shared_ptr multipliedHessian(const Values& x, + const VectorValues& duals) const { + if (!this->active(x)) { + return GaussianFactor::shared_ptr(); + } + const X& x1 = x.at < X > (Base::key()); + const Vector& lambda = duals.at(this->dualKey()); + + std::vector G11; + evaluateHessians(x1, G11); + + if (dim(lambda) != G11.size()) { + throw std::runtime_error( + "Error in evaluateHessians: the number of returned Gij matrices must be the same as the constraint dimension!"); + } + + // Combine the Lagrange-multiplier part into this constraint factor + Matrix lG11sum = zeros(G11[0].rows(), G11[0].cols()); + for (size_t i = 0; i < lambda.rows(); ++i) { + lG11sum += -lambda[i] * G11[i]; + } + + return boost::make_shared(Base::key(), lG11sum, + zero(x1.dim()), 100.0); + } + + /** evaluate Hessians for lambda factors */ + virtual void evaluateHessians(const X& x1, std::vector& G11) const { + + static const bool debug = false; + + boost::function vecH1( + boost::bind(&This::vectorizeH1t, this, _1)); + + Matrix G11all = numericalDerivative11(vecH1, x1, 1e-5); + + if (debug) { + gtsam::print(G11all, "G11all: "); + std::cout << "x1dim: " << x1.dim() << std::endl; + } + + for (size_t i = 0; i < Base::get_noiseModel()->dim(); ++i) { + G11.push_back(G11all.block(i * x1.dim(), 0, x1.dim(), x1.dim())); + if (debug) + gtsam::print(G11[i], "G11: "); + } + } + +private: + /** Vectorize the transpose of Jacobian H1 to compute the Hessian numerically */ + Vector vectorizeH1t(const X& x1) const { + Matrix H1; + evaluateError(x1, H1); + Matrix H1t = H1.transpose(); + H1t.resize(H1t.size(), 1); + return Vector(H1t); + } + + /** Serialization function */ + friend class boost::serialization::access; + template + void serialize(ARCHIVE & ar, const unsigned int version) { + ar + & boost::serialization::make_nvp("NoiseModelFactor", + boost::serialization::base_object(*this)); + } +}; +// \class NonlinearConstraint1 + +/* ************************************************************************* */ +/** A convenient base class for creating your own NonlinearConstraint with 2 + * variables. To derive from this class, implement evaluateError(). */ +template +class NonlinearConstraint2: public NoiseModelFactor2, public NonlinearConstraint { + +public: + + // typedefs for value types pulled from keys + typedef VALUE1 X1; + typedef VALUE2 X2; + +protected: + + typedef NoiseModelFactor2 Base; + typedef NonlinearConstraint2 This; + +public: + + /** + * Default Constructor for I/O + */ + NonlinearConstraint2() { + } + + /** + * Constructor + * @param j1 key of the first variable + * @param j2 key of the second variable + * @param constraintDim number of dimensions of the constraint error function + */ + NonlinearConstraint2(Key j1, Key j2, Key dualKey, size_t constraintDim = 1) : + Base(noiseModel::Constrained::All(constraintDim), j1, j2), NonlinearConstraint(dualKey) { + } + + virtual ~NonlinearConstraint2() { + } + + /** + * Override this method to finish implementing a binary factor. + * If any of the optional Matrix reference arguments are specified, it should compute + * both the function evaluation and its derivative(s) in X1 (and/or X2). + */ + virtual Vector + evaluateError(const X1&, const X2&, boost::optional H1 = boost::none, + boost::optional H2 = boost::none) const = 0; + + /** produce a Gaussian factor graph containing all Hessian Factors, scaled by lambda, + * corresponding to this constraint */ + virtual GaussianFactor::shared_ptr multipliedHessian(const Values& x, + const VectorValues& duals) const { + if (!this->active(x)) { + return GaussianFactor::shared_ptr(); + } + const X1& x1 = x.at < X1 > (Base::keys_[0]); + const X2& x2 = x.at < X2 > (Base::keys_[1]); + const Vector& lambda = duals.at(this->dualKey()); + + std::vector G11, G12, G22; + evaluateHessians(x1, x2, G11, G12, G22); + + if (dim(lambda) != G11.size() || dim(lambda) != G12.size() + || dim(lambda) != G22.size()) { + throw std::runtime_error( + "Error in evaluateHessians: the number of returned Gij matrices must be the same as the constraint dimension!"); + } + + // Combine the Lagrange-multiplier part into this constraint factor + Matrix lG11sum = zeros(G11[0].rows(), G11[0].cols()), lG12sum = zeros( + G12[0].rows(), G12[0].cols()), lG22sum = zeros(G22[0].rows(), + G22[0].cols()); + for (size_t i = 0; i < lambda.rows(); ++i) { + lG11sum += -lambda[i] * G11[i]; + lG12sum += -lambda[i] * G12[i]; + lG22sum += -lambda[i] * G22[i]; + } + + return boost::make_shared(Base::keys_[0], Base::keys_[1], + lG11sum, lG12sum, zero(x1.dim()), lG22sum, zero(x2.dim()), 0.0); + } + + /** evaluate Hessians for lambda factors */ + virtual void evaluateHessians(const X1& x1, const X2& x2, + std::vector& G11, std::vector& G12, + std::vector& G22) const { + + static const bool debug = false; + + boost::function vecH1( + boost::bind(&This::vectorizeH1t, this, _1, _2)), vecH2( + boost::bind(&This::vectorizeH2t, this, _1, _2)); + + Matrix G11all = numericalDerivative21(vecH1, x1, x2, 1e-5); + Matrix G12all = numericalDerivative22(vecH1, x1, x2, 1e-5); + Matrix G22all = numericalDerivative22(vecH2, x1, x2, 1e-5); + + if (debug) { + gtsam::print(G11all, "G11all: "); + gtsam::print(G12all, "G12all: "); + gtsam::print(G22all, "G22all: "); + std::cout << "x1dim: " << x1.dim() << std::endl; + std::cout << "x2dim: " << x2.dim() << std::endl; + } + + for (size_t i = 0; i < Base::get_noiseModel()->dim(); ++i) { + G11.push_back(G11all.block(i * x1.dim(), 0, x1.dim(), x1.dim())); + if (debug) + gtsam::print(G11[i], "G11: "); + + G12.push_back(G12all.block(i * x1.dim(), 0, x1.dim(), x2.dim())); + if (debug) + gtsam::print(G12[i], "G12: "); + + G22.push_back(G22all.block(i * x2.dim(), 0, x2.dim(), x2.dim())); + if (debug) + gtsam::print(G22[i], "G22: "); + } + } + +private: + /** Vectorize the transpose of Jacobian H1 to compute the Hessian numerically */ + Vector vectorizeH1t(const X1& x1, const X2& x2) const { + Matrix H1; + Vector err = evaluateError(x1, x2, H1); + Matrix H1t = H1.transpose(); + H1t.resize(H1t.size(), 1); + return Vector(H1t); + } + + /** Vectorize the transpose of Jacobian H2 to compute the Hessian numerically */ + Vector vectorizeH2t(const X1& x1, const X2& x2) const { + Matrix H2; + Vector err = evaluateError(x1, x2, boost::none, H2); + Matrix H2t = H2.transpose(); + H2t.resize(H2t.size(), 1); + return Vector(H2t); + } + + /** Serialization function */ + friend class boost::serialization::access; + template + void serialize(ARCHIVE & ar, const unsigned int version) { + ar + & boost::serialization::make_nvp("NoiseModelFactor", + boost::serialization::base_object(*this)); + } +}; +// \class NonlinearConstraint2 + +/* ************************************************************************* */ +/** A convenient base class for creating your own NonlinearConstraint with 3 + * variables. To derive from this class, implement evaluateError(). */ +template +class NonlinearConstraint3: public NoiseModelFactor3, public NonlinearConstraint { + +public: + + // typedefs for value types pulled from keys + typedef VALUE1 X1; + typedef VALUE2 X2; + typedef VALUE3 X3; + +protected: + + typedef NoiseModelFactor3 Base; + typedef NonlinearConstraint3 This; + +public: + + /** + * Default Constructor for I/O + */ + NonlinearConstraint3() { + } + + /** + * Constructor + * @param j1 key of the first variable + * @param j2 key of the second variable + * @param constraintDim number of dimensions of the constraint error function + */ + NonlinearConstraint3(Key j1, Key j2, Key j3, Key dualKey, size_t constraintDim = 1) : + Base(noiseModel::Constrained::All(constraintDim), j1, j2, j3), NonlinearConstraint(dualKey) { + } + + virtual ~NonlinearConstraint3() { + } + + /** + * Override this method to finish implementing a binary factor. + * If any of the optional Matrix reference arguments are specified, it should compute + * both the function evaluation and its derivative(s) in X1 (and/or X2). + */ + virtual Vector + evaluateError(const X1&, const X2&, const X3&, boost::optional H1 = + boost::none, boost::optional H2 = boost::none, + boost::optional H3 = boost::none) const = 0; + + /** produce a Gaussian factor graph containing all Hessian Factors, scaled by lambda, + * corresponding to this constraint */ + virtual GaussianFactor::shared_ptr multipliedHessian( + const Values& x, const VectorValues& duals) const { + if (!this->active(x)) { + return GaussianFactor::shared_ptr(); + } + const X1& x1 = x.at < X1 > (Base::keys_[0]); + const X2& x2 = x.at < X2 > (Base::keys_[1]); + const X3& x3 = x.at < X3 > (Base::keys_[2]); + const Vector& lambda = duals.at(this->dualKey()); + + std::vector G11, G12, G13, G22, G23, G33; + evaluateHessians(x1, x2, x3, G11, G12, G13, G22, G23, G33); + + if (dim(lambda) != G11.size() || dim(lambda) != G12.size() + || dim(lambda) != G13.size() || dim(lambda) != G22.size() + || dim(lambda) != G23.size() || dim(lambda) != G33.size()) { + throw std::runtime_error( + "Error in evaluateHessians: the number of returned Gij matrices must be the same as the constraint dimension!"); + } + + // Combine the Lagrange-multiplier part into this constraint factor + Matrix lG11sum = zeros(G11[0].rows(), G11[0].cols()), lG12sum = zeros( + G12[0].rows(), G12[0].cols()), lG13sum = zeros(G13[0].rows(), + G13[0].cols()), lG22sum = zeros(G22[0].rows(), G22[0].cols()), lG23sum = + zeros(G23[0].rows(), G23[0].cols()), lG33sum = zeros(G33[0].rows(), + G33[0].cols()); + for (size_t i = 0; i < lambda.rows(); ++i) { + lG11sum += -lambda[i] * G11[i]; + lG12sum += -lambda[i] * G12[i]; + lG13sum += -lambda[i] * G13[i]; + lG22sum += -lambda[i] * G22[i]; + lG23sum += -lambda[i] * G23[i]; + lG33sum += -lambda[i] * G33[i]; + } + + return boost::shared_ptr( + new HessianFactor(Base::keys_[0], Base::keys_[1], Base::keys_[2], + lG11sum, lG12sum, lG13sum, zero(x1.dim()), lG22sum, lG23sum, + zero(x2.dim()), lG33sum, zero(x3.dim()), 0.0)); + } + + /** + * Default Hessian computation using numerical derivatives + * + * As an example, assuming we have f(x1,x2,x3) where dim(f) = 2, dim(x1) = 3, dim(x2) = 2, dim(x3) = 1 + * + * The Jacobian is: + * f1x1 f1x1 f1x1 | f1x2 f1x2 | f1x3 + * f2x1 f2x1 f2x1 | f2x2 f2x2 | f2x3 + * + * We transpose it to have the gradients: + * f1x1 f2x1 + * f1x1 f2x1 + * f1x1 f2x1 + * f1x2 f2x2 + * f1x2 f2x2 + * f1x3 f2x3 + * Then we vectorize this gradient to have: + * [f1x1 + * f1x1 + * f1x1 + * f1x2 + * f1x2 + * f1x3 + * f2x1 + * f2x1 + * f2x1 + * f2x2 + * f2x2 + * f2x3] + * + * The Derivative of this gradient is then: + * [f1x1x1 f1x1x1 f1x1x1 | f1x1x2 f1x1x2 | f1x1x3 + * f1x1x1 f1x1x1 f1x1x1 | f1x1x2 f1x1x2 | f1x1x3 + * f1x1x1 f1x1x1 f1x1x1 | f1x1x2 f1x1x2 | f1x1x3 + * ---------------------|---------------|------- + * f1x2x1 f1x2x1 f1x2x1 | f1x2x2 f1x2x2 | f1x2x3 + * f1x2x1 f1x2x1 f1x2x1 | f1x2x2 f1x2x2 | f1x2x3 + * ---------------------|---------------|------- + * f1x3x1 f1x3x1 f1x3x1 | f1x3x2 f1x3x2 | f1x3x3 + * ============================================= + * f2x1x1 f2x1x1 f2x1x1 | f2x1x2 f2x1x2 | f2x1x3 + * f2x1x1 f2x1x1 f2x1x1 | f2x1x2 f2x1x2 | f2x1x3 + * f2x1x1 f2x1x1 f2x1x1 | f2x1x2 f2x1x2 | f2x1x3 + * ---------------------|---------------|------- + * f2x2x1 f2x2x1 f2x2x1 | f2x2x2 f2x2x2 | f2x2x3 + * f2x2x1 f2x2x1 f2x2x1 | f2x2x2 f2x2x2 | f2x2x3 + * ---------------------|---------------|------- + * f2x3x1 f2x3x1 f2x3x1 | f2x3x2 f2x3x2 | f2x3x3 ] + * + * It is the combination of the Hessian of each component of f + * stacking on top of each other. + * + * */ + virtual void evaluateHessians(const X1& x1, const X2& x2, const X3& x3, + std::vector& G11, std::vector& G12, + std::vector& G13, std::vector& G22, + std::vector& G23, std::vector& G33) const { + + static const bool debug = false; + + boost::function vecH1( + boost::bind(&This::vectorizeH1t, this, _1, _2, _3)), vecH2( + boost::bind(&This::vectorizeH2t, this, _1, _2, _3)), vecH3( + boost::bind(&This::vectorizeH3t, this, _1, _2, _3)); + + Matrix G11all = numericalDerivative31(vecH1, x1, x2, x3, 1e-5); + Matrix G12all = numericalDerivative32(vecH1, x1, x2, x3, 1e-5); + Matrix G13all = numericalDerivative33(vecH1, x1, x2, x3, 1e-5); + Matrix G22all = numericalDerivative32(vecH2, x1, x2, x3, 1e-5); + Matrix G23all = numericalDerivative33(vecH2, x1, x2, x3, 1e-5); + Matrix G33all = numericalDerivative33(vecH3, x1, x2, x3, 1e-5); + + if (debug) { + gtsam::print(G11all, "G11all: "); + gtsam::print(G12all, "G12all: "); + gtsam::print(G13all, "G13all: "); + gtsam::print(G22all, "G22all: "); + gtsam::print(G23all, "G23all: "); + gtsam::print(G33all, "G33all: "); + std::cout << "x1dim: " << x1.dim() << std::endl; + std::cout << "x2dim: " << x2.dim() << std::endl; + std::cout << "x3dim: " << x3.dim() << std::endl; + } + + for (size_t i = 0; i < Base::get_noiseModel()->dim(); ++i) { + G11.push_back(G11all.block(i * x1.dim(), 0, x1.dim(), x1.dim())); + if (debug) + gtsam::print(G11[i], "G11: "); + + G12.push_back(G12all.block(i * x1.dim(), 0, x1.dim(), x2.dim())); + if (debug) + gtsam::print(G12[i], "G12: "); + + G13.push_back(G13all.block(i * x1.dim(), 0, x1.dim(), x3.dim())); + if (debug) + gtsam::print(G13[i], "G13: "); + + G22.push_back(G22all.block(i * x2.dim(), 0, x2.dim(), x2.dim())); + if (debug) + gtsam::print(G22[i], "G22: "); + + G23.push_back(G23all.block(i * x2.dim(), 0, x2.dim(), x3.dim())); + if (debug) + gtsam::print(G23[i], "G23: "); + + G33.push_back(G33all.block(i * x3.dim(), 0, x3.dim(), x3.dim())); + if (debug) + gtsam::print(G33[i], "G33: "); + } + } + +private: + /** Vectorize the transpose of Jacobian H1 to compute the Hessian numerically */ + Vector vectorizeH1t(const X1& x1, const X2& x2, const X3& x3) const { + Matrix H1; + Vector err = evaluateError(x1, x2, x3, H1); + Matrix H1t = H1.transpose(); + H1t.resize(H1t.size(), 1); + return Vector(H1t); + } + + /** Vectorize the transpose of Jacobian H2 to compute the Hessian numerically */ + Vector vectorizeH2t(const X1& x1, const X2& x2, const X3& x3) const { + Matrix H2; + Vector err = evaluateError(x1, x2, x3, boost::none, H2); + Matrix H2t = H2.transpose(); + H2t.resize(H2t.size(), 1); + return Vector(H2t); + } + + /** Vectorize the transpose of Jacobian H3 to compute the Hessian numerically */ + Vector vectorizeH3t(const X1& x1, const X2& x2, const X3& x3) const { + Matrix H3; + Vector err = evaluateError(x1, x2, x3, boost::none, boost::none, H3); + Matrix H3t = H3.transpose(); + H3t.resize(H3t.size(), 1); + return Vector(H3t); + } + + /** Serialization function */ + friend class boost::serialization::access; + template + void serialize(ARCHIVE & ar, const unsigned int version) { + ar + & boost::serialization::make_nvp("NoiseModelFactor", + boost::serialization::base_object(*this)); + } +}; +// \class NonlinearConstraint3 + +} /* namespace gtsam */ diff --git a/gtsam_unstable/nonlinear/tests/testNonlinearConstraints.cpp b/gtsam_unstable/nonlinear/tests/testNonlinearConstraints.cpp new file mode 100644 index 000000000..522f4096d --- /dev/null +++ b/gtsam_unstable/nonlinear/tests/testNonlinearConstraints.cpp @@ -0,0 +1,176 @@ +/** + * @file testNonlinearConstraints.cpp + * @brief + * @author Duy-Nguyen Ta + * @date Oct 26, 2013 + */ + +#include +#include +#include +#include +#include +#include "../NonlinearConstraint.h" + +using namespace gtsam; +using namespace gtsam::symbol_shorthand; +using namespace std; + +TEST(Vector, vectorize) { + Matrix m = (Matrix(4,3) << 1,2,3,4,5,6,7,8,9,10,11,12).finished(); + Matrix m2 = m.transpose(); + m2.resize(12,1); + Vector v = Vector(m2); + Vector expectedV = (Vector(12) << 1,2,3,4,5,6,7,8,9,10,11,12).finished(); + EXPECT(assert_equal(expectedV, v)); +} + +/* ************************************************************************* */ +/** + * Test 3-way equality nonlinear constraint + * x(1) + x(2)^2 + x(3)^3 = 3 + */ +class Constraint: public NonlinearConstraint3 { + typedef NonlinearConstraint3 Base; + +public: + Constraint(Key key1, Key key2, Key key3, Key dualKey) : + Base(key1, key2, key3, dualKey, 1) { + } + + Vector evaluateError(const LieScalar& x1, const LieScalar& x2, + const LieScalar& x3, boost::optional H1 = boost::none, + boost::optional H2 = boost::none, boost::optional H3 = + boost::none) const { + if (H1) { + *H1 = (Matrix(1, 1) << 1.0).finished(); + } + if (H2) { + *H2 = (Matrix(1, 1) << 2.0 * x2.value()).finished(); + } + if (H3) { + *H3 = (Matrix(1, 1) << 3.0 * x3.value() * x3.value()).finished(); + } + + return (Vector(1) << + x1.value() + x2.value() * x2.value() + + x3.value() * x3.value() * x3.value() - 3.0).finished(); + } + + void expectedHessians(const LieScalar& x1, const LieScalar& x2, + const LieScalar& x3, std::vector& G11, std::vector& G12, + std::vector& G13, std::vector& G22, + std::vector& G23, std::vector& G33) const { + G11.push_back(zeros(1, 1)); + G12.push_back(zeros(1, 1)); + G13.push_back(zeros(1, 1)); + G22.push_back((Matrix(1, 1) << 2.0).finished()); + G23.push_back(zeros(1, 1)); + G33.push_back((Matrix(1, 1) << 6.0 * x3.value()).finished()); + } +}; + +/* ************************************************************************* */ + +TEST(NonlinearConstraint3, Hessian) { + LieScalar x1(2.0), x2(sqrt(2) - 1), x3(sqrt(2) - 1), x4(2.0), x5(0.5); + Constraint constraint(X(1), X(2), X(3), 0); + + std::vector expectedG11, expectedG12, expectedG13, expectedG22, + expectedG23, expectedG33; + constraint.expectedHessians(x1, x2, x3, expectedG11, expectedG12, expectedG13, + expectedG22, expectedG23, expectedG33); + std::vector G11, G12, G13, G22, G23, G33; + constraint.evaluateHessians(x1, x2, x3, G11, G12, G13, G22, G23, G33); + EXPECT(assert_equal(expectedG11[0], G11[0], 1e-5)); + EXPECT(assert_equal(expectedG12[0], G12[0], 1e-5)); + EXPECT(assert_equal(expectedG13[0], G13[0], 1e-5)); + EXPECT(assert_equal(expectedG22[0], G22[0], 1e-5)); + EXPECT(assert_equal(expectedG23[0], G23[0], 1e-5)); + EXPECT(assert_equal(expectedG33[0], G33[0], 1e-5)); +} + +/* ************************************************************************* */ +/** + * Test 3-way nonlinear equality multi-values constraint + * p(1).x + p(2).x^2 + p(3).y^3 = 3 + * p(1).x^2 + p(2).y + p(3).x^3 = 5 + */ +class Constraint2d: public NonlinearConstraint3 { + typedef NonlinearConstraint3 Base; + +public: + Constraint2d(Key key1, Key key2, Key key3, Key dualKey) : + Base(key1, key2, key3, dualKey, 1) { + } + + Vector evaluateError(const Point2& p1, const Point2& p2, const Point2& p3, + boost::optional H1 = boost::none, boost::optional H2 = + boost::none, boost::optional H3 = boost::none) const { + if (H1) { + *H1 = (Matrix(2, 2) << 1.0, 0.0, 2 * p1.x(), 0.0).finished(); + } + if (H2) { + *H2 = (Matrix(2, 2) << 2.0 * p2.x(), 0.0, 0.0, 1.0).finished(); + } + if (H3) { + *H3 = (Matrix(2, 2) << 0.0, 3.0 * p3.y() * p3.y(), 3.0 * p3.x() * p3.x(), 0.0).finished(); + } + + return (Vector(2) << p1.x() + p2.x() * p2.x() + p3.y() * p3.y() * p3.y() - 3.0, + p1.x() * p1.x() + p2.y() + p3.x() * p3.x() * p3.x() - 5.0).finished(); + } + + void expectedHessians(const Point2& x1, const Point2& x2, const Point2& x3, + std::vector& G11, std::vector& G12, + std::vector& G13, std::vector& G22, + std::vector& G23, std::vector& G33) const { + G11.push_back(zeros(2, 2)); + G11.push_back((Matrix(2,2) << 2.0, 0.0, 0.0, 0.0).finished()); + + G12.push_back(zeros(2,2)); + G12.push_back(zeros(2,2)); + + G13.push_back(zeros(2,2)); + G13.push_back(zeros(2,2)); + + G22.push_back((Matrix(2,2) << 2.0, 0.0, 0.0, 0.0).finished()); + G22.push_back(zeros(2,2)); + + G23.push_back(zeros(2, 2)); + G23.push_back(zeros(2, 2)); + + G33.push_back((Matrix(2, 2) << 0.0, 0.0, 0.0, 6.0 * x3.y() ).finished()); + G33.push_back((Matrix(2, 2) << 6.0 * x3.x(), 0.0, 0.0, 0.0 ).finished()); + } +}; + +/* ************************************************************************* */ + +TEST(NonlinearConstraint3, Hessian2) { + Point2 x1(2.0, 2.0), x2(sqrt(2) - 1, 3.0), x3(4.0, sqrt(2) - 2); + Constraint2d constraint(X(1), X(2), X(3), 0); + + std::vector expectedG11, expectedG12, expectedG13, expectedG22, + expectedG23, expectedG33; + constraint.expectedHessians(x1, x2, x3, expectedG11, expectedG12, expectedG13, + expectedG22, expectedG23, expectedG33); + std::vector G11, G12, G13, G22, G23, G33; + constraint.evaluateHessians(x1, x2, x3, G11, G12, G13, G22, G23, G33); + for (size_t i = 0; i Date: Thu, 18 Dec 2014 20:49:08 -0500 Subject: [PATCH 004/130] Added Simple QP solver and test. Unit test doesn't work yet --- .../nonlinear/tests/testSQPSimple.cpp | 273 ++++++++++++++++++ 1 file changed, 273 insertions(+) create mode 100644 gtsam_unstable/nonlinear/tests/testSQPSimple.cpp diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp new file mode 100644 index 000000000..c649cf295 --- /dev/null +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -0,0 +1,273 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ + +/** + * @file testQPSimple.cpp + * @brief Unit tests for testQPSimple + * @author Krunal Chande + * @author Duy-Nguyen Ta + * @author Luca Carlone + * @date Dec 15, 2014 + */ + +#include +#include +#include +#include +#include +#include +#include + +namespace gtsam { + +class LinearEqualityManifoldFactorGraph: public FactorGraph { +public: + /// default constructor + LinearEqualityManifoldFactorGraph() { + } + + /// linearize to a LinearEqualityFactorGraph + LinearEqualityFactorGraph::shared_ptr linearize( + const Values& linearizationPoint) const { + LinearEqualityFactorGraph::shared_ptr linearGraph( + new LinearEqualityFactorGraph()); + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + JacobianFactor::shared_ptr jacobian = boost::dynamic_pointer_cast( + factor->linearize(linearizationPoint)); + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + linearGraph->add(LinearEquality(*jacobian, constraint->dualKey())); + } + return linearGraph; + } + + /** + * Return true if the max absolute error all factors is less than a tolerance + */ + bool checkFeasibility(const Values& values, double tol) const { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( + factor); + Vector error = noiseModelFactor->unwhitenedError(values); + if (error.lpNorm() > tol) { + return false; + } + } + return true; + } +}; + +class NonlinearEqualityFactorGraph: public LinearEqualityManifoldFactorGraph { +public: + /// default constructor + NonlinearEqualityFactorGraph() { + } + + GaussianFactorGraph::shared_ptr multipliedHessians(const Values& values, const VectorValues& duals) const { + GaussianFactorGraph::shared_ptr constrainedHessians(new GaussianFactorGraph()); + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this) { + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + constrainedHessians->push_back(constraint->multipliedHessian(values, duals)); + } + return constrainedHessians; + } +}; + + +struct NLP { + NonlinearFactorGraph cost; + NonlinearEqualityFactorGraph linearEqualities; + NonlinearEqualityFactorGraph nonlinearEqualities; +}; + +struct SQPSimpleState { + Values values; + VectorValues duals; + bool converged; + + /// Default constructor + SQPSimpleState() : values(), duals(), converged(false) {} + + /// Constructor with an initialValues + SQPSimpleState(const Values& initialValues) : + values(initialValues), duals(VectorValues()), converged(false) { + } +}; + +/** + * Simple SQP optimizer to solve nonlinear constrained problems. + * This simple version won't care about nonconvexity, which needs + * more advanced techniques to solve, e.g., merit function, line search, second-order correction etc. + */ +class SQPSimple { + NLP nlp_; + static const double errorTol = 1e-5; +public: + SQPSimple(const NLP& nlp) : + nlp_(nlp) { + } + + /// Check if \nabla f(x) - \lambda * \nabla c(x) == 0 + bool isDualFeasible(const VectorValues& delta) const { + return delta.vector().lpNorm() < errorTol; + } + + /// Check if c(x) == 0 + bool isPrimalFeasible(const SQPSimpleState& state) const { + return nlp_.linearEqualities.checkFeasibility(state.values, errorTol) + && nlp_.nonlinearEqualities.checkFeasibility(state.values, errorTol); + } + + /// Check convergence + bool checkConvergence(const SQPSimpleState& state, const VectorValues& delta) const { + return isPrimalFeasible(state) & isDualFeasible(delta); + } + + /** + * Single iteration of SQP + */ + SQPSimpleState iterate(const SQPSimpleState& state) const { + // construct the qp subproblem + QP qp; + qp.cost = *nlp_.cost.linearize(state.values); + GaussianFactorGraph::shared_ptr multipliedHessians = nlp_.nonlinearEqualities.multipliedHessians(state.values, state.duals); + qp.cost.push_back(*multipliedHessians); + + qp.equalities.add(*nlp_.linearEqualities.linearize(state.values)); + qp.equalities.add(*nlp_.nonlinearEqualities.linearize(state.values)); + + // solve the QP subproblem + VectorValues delta, duals; + QPSolver qpSolver(qp); + boost::tie(delta, duals) = qpSolver.optimize(); + + // update new state + SQPSimpleState newState; + newState.values = state.values.retract(delta); + newState.duals = duals; + newState.converged = checkConvergence(newState, delta); + + return newState; + } + + /** + * Main optimization function. + */ + std::pair optimize(const Values& initialValues) const { + SQPSimpleState state(initialValues); + while (!state.converged) { + state = iterate(state); + } + + return std::make_pair(initialValues, VectorValues()); + } + + +}; +} + +using namespace std; +using namespace gtsam::symbol_shorthand; +using namespace gtsam; +const double tol = 1e-10; +//****************************************************************************** +TEST(testSQPSimple, Problem2) { + // Simple quadratic cost: x1^2 + x2^2 + // Note the Hessian encodes: + // 0.5*x1'*G11*x1 + x1'*G12*x2 + 0.5*x2'*G22*x2 - x1'*g1 - x2'*g2 + 0.5*f + // Hence here we have G11 = 2, G12 = 0, G22 = 2, g1 = 0, g2 = 0, f = 0 + HessianFactor hf(X(1), X(2), 2.0 * ones(1,1), zero(1), zero(1), + 2*ones(1,1), zero(1) , 0); + + LinearEqualityFactorGraph equalities; + equalities.push_back(LinearEquality(X(1), ones(1), X(2), ones(1), -1*ones(1), 0)); // x + y - 1 = 0 + + // Compare against QP + QP qp; + qp.cost.add(hf); + qp.equalities = equalities; + + // instantiate QPsolver + QPSolver qpSolver(qp); + // create initial values for optimization + VectorValues initialVectorValues; + initialVectorValues.insert(X(1), zero(1)); + initialVectorValues.insert(X(2), zero(1)); + VectorValues expectedSolution = qpSolver.optimize(initialVectorValues).first; + cout<<"expectedSolution.at(X(1))[0]: "< x1 + x2 -2 <= 0, --> b=2 +// inequalities.push_back(LinearInequality(X(1), -ones(1, 1), 0, 1)); // -x1 <= 0 +// inequalities.push_back(LinearInequality(X(2), -ones(1, 1), 0, 2)); // -x2 <= 0 +// inequalities.push_back(LinearInequality(X(1), ones(1, 1), 1.5, 3)); // x1 <= 3/2 +// +// // Compare against a QP +// QP qp; +// qp.cost.add(lf); +// qp.inequalities = inequalities; +// +// // instantiate QPsolver +// QPSolver qpSolver(qp); +// // create initial values for optimization +// VectorValues initialVectorValues; +// initialVectorValues.insert(X(1), zero(1)); +// initialVectorValues.insert(X(2), zero(1)); +// VectorValues expectedSolution = qpSolver.optimize(initialVectorValues).first; +// +// +// NonlinearEqualityFactorGraph linearEqualities; +// NonlinearEqualityFactorGraph nonlinearEqualities; +// nonlinearEqualities.push_back(NonlinearEquality(X(1), ones(1, 1), X(2), ones(1, 1), 2, 0))); +// +// NLP nlp; +// nlp.cost.add(lf); +// nlp.linearEqualities.push_back(NonlinearEqualityFactorGraph()); +// // instantiate QPsolver +// SQPSimple sqpSolver(nlp); +// // create initial values for optimization +// Values initialValues; +// initialValues.insert(X(1), zero(1)); +// initialValues.insert(X(2), zero(1)); +// +// std::pair actualSolution = sqpSolver.optimize(initialValues); +// +// DOUBLES_EQUAL(expectedSolution.at(X(1))[0], actualSolution.at(X(1)), +// tol); +// DOUBLES_EQUAL(expectedSolution.at(X(2))[0], actualSolution.at(X(2)), +// tol); +//} From 6d76b5910c7a0a5738a10daf8c34d1af718a60e1 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Fri, 19 Dec 2014 12:47:52 -0500 Subject: [PATCH 005/130] added comments --- gtsam_unstable/nonlinear/tests/testSQPSimple.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index c649cf295..923d7be4c 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -206,7 +206,9 @@ TEST(testSQPSimple, Problem2) { //Instantiate NLP NLP nlp; - nlp.cost.add(); + nlp.cost.add(); // wrap it using linearcontainerfactor + nlp.linearEqualities // for constraint it has to inherit from + // write an evaluate error and return jacobian // Instantiate SQP SQPSimple sqpSimple(nlp); From 147d174a663a9ac3e5d900ac0828099de34cf5f3 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Fri, 19 Dec 2014 12:48:22 -0500 Subject: [PATCH 006/130] test files deprecated, might be removed later --- .../tests/testLinearInequalityManifolds.cpp | 65 +++++++ .../slam/tests/testPositionConstraints.cpp | 176 ++++++++++++++++++ 2 files changed, 241 insertions(+) create mode 100644 gtsam_unstable/nonlinear/tests/testLinearInequalityManifolds.cpp create mode 100644 gtsam_unstable/slam/tests/testPositionConstraints.cpp diff --git a/gtsam_unstable/nonlinear/tests/testLinearInequalityManifolds.cpp b/gtsam_unstable/nonlinear/tests/testLinearInequalityManifolds.cpp new file mode 100644 index 000000000..8b1077ebf --- /dev/null +++ b/gtsam_unstable/nonlinear/tests/testLinearInequalityManifolds.cpp @@ -0,0 +1,65 @@ +///* ---------------------------------------------------------------------------- +// +// * GTSAM Copyright 2010, Georgia Tech Research Corporation, +// * Atlanta, Georgia 30332-0415 +// * All Rights Reserved +// * Authors: Frank Dellaert, et al. (see THANKS for the full author list) +// +// * See LICENSE for the license information +// +// * -------------------------------------------------------------------------- */ +// +///** +// * @file testLinearInequalityManifolds.cpp +// * @brief Unit tests for testLinearInequalityManifolds +// * @author Krunal Chande +// * @author Duy-Nguyen Ta +// * @author Luca Carlone +// * @date Dec 15, 2014 +// */ +// +//#include +//#include +//#include +// +//using namespace std; +//using namespace gtsam::symbol_shorthand; +//using namespace gtsam; +//const double tol = 1e-10; +// +//namespace gtsam { +//class PositionConstraint : public LinearInequalityManifold1 { +//public: +// PositionConstraint() +//}; +//} +////****************************************************************************** +//TEST(testLinearEqualityManifolds, equals ) { +// // Instantiate a class LinearEqualityManifolds +// LinearEqualityManifolds le1(); +// +// // Instantiate another class LinearEqualityManifolds +// LinearEqualityManifolds le2(); +// +// // check equals +// CHECK(assert_equal(le1,le1)); +// CHECK(le2.equals(le2)); +// CHECK(!le2.equals(le1)); +// CHECK(!le1.equals(le2)); +//} +// +//////****************************************************************************** +////TEST(testLinearEqualityManifolds, linearize ) { +////} +//// +//////****************************************************************************** +////TEST(testLinearEqualityManifolds, size ) { +////} +// +////****************************************************************************** +//int main() { +// TestResult tr; +// return TestRegistry::runAllTests(tr); +//} +////****************************************************************************** +// diff --git a/gtsam_unstable/slam/tests/testPositionConstraints.cpp b/gtsam_unstable/slam/tests/testPositionConstraints.cpp new file mode 100644 index 000000000..04eb4a260 --- /dev/null +++ b/gtsam_unstable/slam/tests/testPositionConstraints.cpp @@ -0,0 +1,176 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ + +/** + * @file testPositionUpperBoundX.cpp + * @brief Unit tests for testPositionUpperBoundX + * @author Krunal Chande + * @author Duy-Nguyen Ta + * @author Luca Carlone + * @date Dec 16, 2014 + */ + +#include +#include +#include +#include +#include +#include + +namespace gtsam { + +template +class LinearInequalityManifold1: public NoiseModelFactor1 { + typedef NoiseModelFactor1 Base; + typedef LinearInequalityManifold1 This; +protected: + bool active_; +public: + + /// Constructor + LinearInequalityManifold1(Key key) : + Base(noiseModel::Constrained::All(1), key), active_(true) { + } + + /** implement functions needed for Testable */ + + /** print */ + virtual void print(const std::string& s, const KeyFormatter& keyFormatter = + DefaultKeyFormatter) const { + std::cout << s << " key = { " << keyFormatter(this->key()) << "}" + << std::endl; + if (active_) + std::cout << " Active" << std::endl; + else + std::cout << " Inactive" << std::endl; + } + + /** equals */ + virtual bool equals(const NonlinearFactor& expected, + double tol = 1e-9) const { + const This* e = dynamic_cast(&expected); + return e != NULL && Base::equals(*e, tol) && (active_ == e->active_); + } + + /** Evaluate error MUST return a *one dimensional* vector, + * because we don't support multi-valued inequality factors + */ + virtual Vector evaluateError(const T& x, boost::optional H = + boost::none) const = 0; + +}; + +/** + * This class defines an inequality upper bound on x for a Pose3 + * x <= upperBound + */ +class PositionUpperBoundX: public LinearInequalityManifold1 { + double upperBound_; + typedef LinearInequalityManifold1 Base; + typedef PositionUpperBoundX This; +public: + PositionUpperBoundX(Key key, const double upperBound) : + Base(key), upperBound_(upperBound) { + } + + /// @return a deep copy of this factor + virtual gtsam::NonlinearFactor::shared_ptr clone() const { + return boost::static_pointer_cast( + gtsam::NonlinearFactor::shared_ptr(new This(*this))); + } + + /** implement functions needed for Testable */ + + /** print */ + virtual void print(const std::string& s, const KeyFormatter& keyFormatter = + DefaultKeyFormatter) const { + Base::print(s); + std::cout << "PositionUpperBoundX" << std::endl; + std::cout << "x <= " << upperBound_ << std::endl; + } + + /** equals */ + virtual bool equals(const NonlinearFactor& expected, + double tol = 1e-9) const { + const This* e = dynamic_cast(&expected); + return e != NULL && Base::equals(*e, tol) && (upperBound_ == e->upperBound_); + } + + /** + * error = x - upperBound_ + */ + double computeError(const Pose3& pose) const { + return pose.x() - upperBound_; + } + + /** + * error = x - upperBound_ + */ + Vector evaluateError(const Pose3& pose, boost::optional H = + boost::none) const { + if (H) + *H = (Matrix(1, 6) << 0, 0, 0, 1, 0, 0).finished(); + return (Vector(1) << computeError(pose)).finished(); + } +}; +} + +//****************************************************************************** + +using namespace std; +using namespace gtsam::symbol_shorthand; +using namespace gtsam; + +//****************************************************************************** +TEST(PositionUpperBoundX, equals ) { + // Instantiate a class PositionUpperBoundX + PositionUpperBoundX ineq1(X(1), 1); + + // Instantiate another class PositionUpperBoundX + PositionUpperBoundX ineq2(X(1), 1); + + // check equals + CHECK(ineq1.equals(ineq2, 1e-5)); +} + +//****************************************************************************** +TEST(PositionUpperBoundX, evaluateError ) { + // Instantiate a class PositionUpperBoundX + PositionUpperBoundX ineq(X(1), 45.6); + + Pose3 pose(Rot3::ypr(0.1, 0.3, 0.2), Point3(43.0, 27.8, 91.1)); + Matrix H; + Vector error = ineq.evaluateError(pose, H); + + Matrix expectedH = numericalDerivative11( + boost::bind(&PositionUpperBoundX::evaluateError, ineq, _1, boost::none), + pose, 1e-5); + + cout << "expectedH: " << expectedH << endl; + cout << "H: " << H << endl; + + Vector expectedError = (Vector(1) << 0.0).finished(); +// CHECK(assert_equal(expectedError, error, 1e-100)); +// CHECK(assert_equal(expectedH, H, 1e-100)); + + Matrix hessian = numericalHessian( + boost::bind(&PositionUpperBoundX::computeError, ineq, _1), + pose, 1e-5); + cout << "Hessian: \n" << hessian << endl; +} + +//****************************************************************************** +int main() { + TestResult tr; + return TestRegistry::runAllTests(tr); +} +//****************************************************************************** + From ccc243d37ae1d604e5c5956db9d913fa130777f9 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Fri, 19 Dec 2014 18:05:21 -0500 Subject: [PATCH 007/130] Obtain dim by using traits. --- .../nonlinear/NonlinearConstraint.h | 53 ++++++++++++------- 1 file changed, 33 insertions(+), 20 deletions(-) diff --git a/gtsam_unstable/nonlinear/NonlinearConstraint.h b/gtsam_unstable/nonlinear/NonlinearConstraint.h index 3d1bf7eb0..ce01f89e8 100644 --- a/gtsam_unstable/nonlinear/NonlinearConstraint.h +++ b/gtsam_unstable/nonlinear/NonlinearConstraint.h @@ -7,6 +7,7 @@ #pragma once +#include #include #include #include @@ -51,6 +52,9 @@ protected: typedef NoiseModelFactor1 Base; typedef NonlinearConstraint1 This; +private: + static const int X1Dim = traits::dimension::value; + public: /** @@ -104,7 +108,7 @@ public: } return boost::make_shared(Base::key(), lG11sum, - zero(x1.dim()), 100.0); + zero(X1Dim), 100.0); } /** evaluate Hessians for lambda factors */ @@ -119,11 +123,11 @@ public: if (debug) { gtsam::print(G11all, "G11all: "); - std::cout << "x1dim: " << x1.dim() << std::endl; + std::cout << "x1dim: " << X1Dim << std::endl; } for (size_t i = 0; i < Base::get_noiseModel()->dim(); ++i) { - G11.push_back(G11all.block(i * x1.dim(), 0, x1.dim(), x1.dim())); + G11.push_back(G11all.block(i * X1Dim, 0, X1Dim, X1Dim)); if (debug) gtsam::print(G11[i], "G11: "); } @@ -167,6 +171,10 @@ protected: typedef NoiseModelFactor2 Base; typedef NonlinearConstraint2 This; +private: + static const int X1Dim = traits::dimension::value; + static const int X2Dim = traits::dimension::value; + public: /** @@ -228,7 +236,7 @@ public: } return boost::make_shared(Base::keys_[0], Base::keys_[1], - lG11sum, lG12sum, zero(x1.dim()), lG22sum, zero(x2.dim()), 0.0); + lG11sum, lG12sum, zero(X1Dim), lG22sum, zero(X2Dim), 0.0); } /** evaluate Hessians for lambda factors */ @@ -250,20 +258,20 @@ public: gtsam::print(G11all, "G11all: "); gtsam::print(G12all, "G12all: "); gtsam::print(G22all, "G22all: "); - std::cout << "x1dim: " << x1.dim() << std::endl; - std::cout << "x2dim: " << x2.dim() << std::endl; + std::cout << "x1dim: " << traits::dimension::value << std::endl; + std::cout << "x2dim: " << traits::dimension::value << std::endl; } for (size_t i = 0; i < Base::get_noiseModel()->dim(); ++i) { - G11.push_back(G11all.block(i * x1.dim(), 0, x1.dim(), x1.dim())); + G11.push_back(G11all.block(i * X1Dim, 0, X1Dim, X1Dim)); if (debug) gtsam::print(G11[i], "G11: "); - G12.push_back(G12all.block(i * x1.dim(), 0, x1.dim(), x2.dim())); + G12.push_back(G12all.block(i * X1Dim, 0, X1Dim, X2Dim)); if (debug) gtsam::print(G12[i], "G12: "); - G22.push_back(G22all.block(i * x2.dim(), 0, x2.dim(), x2.dim())); + G22.push_back(G22all.block(i * X2Dim, 0, X2Dim, X2Dim)); if (debug) gtsam::print(G22[i], "G22: "); } @@ -317,6 +325,11 @@ protected: typedef NoiseModelFactor3 Base; typedef NonlinearConstraint3 This; +private: + static const int X1Dim = traits::dimension::value; + static const int X2Dim = traits::dimension::value; + static const int X3Dim = traits::dimension::value; + public: /** @@ -387,8 +400,8 @@ public: return boost::shared_ptr( new HessianFactor(Base::keys_[0], Base::keys_[1], Base::keys_[2], - lG11sum, lG12sum, lG13sum, zero(x1.dim()), lG22sum, lG23sum, - zero(x2.dim()), lG33sum, zero(x3.dim()), 0.0)); + lG11sum, lG12sum, lG13sum, zero(X1Dim), lG22sum, lG23sum, + zero(X2Dim), lG33sum, zero(X3Dim), 0.0)); } /** @@ -470,33 +483,33 @@ public: gtsam::print(G22all, "G22all: "); gtsam::print(G23all, "G23all: "); gtsam::print(G33all, "G33all: "); - std::cout << "x1dim: " << x1.dim() << std::endl; - std::cout << "x2dim: " << x2.dim() << std::endl; - std::cout << "x3dim: " << x3.dim() << std::endl; + std::cout << "x1dim: " << X1Dim << std::endl; + std::cout << "x2dim: " << X2Dim << std::endl; + std::cout << "x3dim: " << X3Dim << std::endl; } for (size_t i = 0; i < Base::get_noiseModel()->dim(); ++i) { - G11.push_back(G11all.block(i * x1.dim(), 0, x1.dim(), x1.dim())); + G11.push_back(G11all.block(i * X1Dim, 0, X1Dim, X1Dim)); if (debug) gtsam::print(G11[i], "G11: "); - G12.push_back(G12all.block(i * x1.dim(), 0, x1.dim(), x2.dim())); + G12.push_back(G12all.block(i * X1Dim, 0, X1Dim, X2Dim)); if (debug) gtsam::print(G12[i], "G12: "); - G13.push_back(G13all.block(i * x1.dim(), 0, x1.dim(), x3.dim())); + G13.push_back(G13all.block(i * X1Dim, 0, X1Dim, X3Dim)); if (debug) gtsam::print(G13[i], "G13: "); - G22.push_back(G22all.block(i * x2.dim(), 0, x2.dim(), x2.dim())); + G22.push_back(G22all.block(i * X2Dim, 0, X2Dim, X2Dim)); if (debug) gtsam::print(G22[i], "G22: "); - G23.push_back(G23all.block(i * x2.dim(), 0, x2.dim(), x3.dim())); + G23.push_back(G23all.block(i * X2Dim, 0, X2Dim, X3Dim)); if (debug) gtsam::print(G23[i], "G23: "); - G33.push_back(G33all.block(i * x3.dim(), 0, x3.dim(), x3.dim())); + G33.push_back(G33all.block(i * X3Dim, 0, X3Dim, X3Dim)); if (debug) gtsam::print(G33[i], "G33: "); } From 29e6e67de774740fb9240afe6a7dfadc2b96fcd8 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Fri, 19 Dec 2014 18:07:11 -0500 Subject: [PATCH 008/130] added debug info, fixed unit test, added nonlinear constraint (circle) test. Doesn't work because of negative definite hessian obtained from multiplying the dual with the constraint hessian. --- .../nonlinear/tests/testSQPSimple.cpp | 220 ++++++++++++------ 1 file changed, 154 insertions(+), 66 deletions(-) diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index 923d7be4c..16bc17b03 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -91,13 +91,14 @@ struct SQPSimpleState { Values values; VectorValues duals; bool converged; + size_t iterations; /// Default constructor - SQPSimpleState() : values(), duals(), converged(false) {} + SQPSimpleState() : values(), duals(), converged(false), iterations(0) {} /// Constructor with an initialValues SQPSimpleState(const Values& initialValues) : - values(initialValues), duals(VectorValues()), converged(false) { + values(initialValues), duals(VectorValues()), converged(false), iterations(0) { } }; @@ -117,6 +118,7 @@ public: /// Check if \nabla f(x) - \lambda * \nabla c(x) == 0 bool isDualFeasible(const VectorValues& delta) const { return delta.vector().lpNorm() < errorTol; +// return false; } /// Check if c(x) == 0 @@ -134,6 +136,8 @@ public: * Single iteration of SQP */ SQPSimpleState iterate(const SQPSimpleState& state) const { + static const bool debug = true; + // construct the qp subproblem QP qp; qp.cost = *nlp_.cost.linearize(state.values); @@ -143,30 +147,56 @@ public: qp.equalities.add(*nlp_.linearEqualities.linearize(state.values)); qp.equalities.add(*nlp_.nonlinearEqualities.linearize(state.values)); + if (debug) + qp.print("QP subproblem:"); + // solve the QP subproblem VectorValues delta, duals; QPSolver qpSolver(qp); boost::tie(delta, duals) = qpSolver.optimize(); + if (debug) + delta.print("delta = "); + + if (debug) + duals.print("duals = "); + // update new state SQPSimpleState newState; newState.values = state.values.retract(delta); newState.duals = duals; newState.converged = checkConvergence(newState, delta); + newState.iterations = state.iterations + 1; return newState; } + VectorValues initializeDuals() const { + VectorValues duals; + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearEqualities) { + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + duals.insert(constraint->dualKey(), zero(factor->dim())); + } + + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.nonlinearEqualities) { + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + duals.insert(constraint->dualKey(), zero(factor->dim())); + } + return duals; + } + /** * Main optimization function. */ std::pair optimize(const Values& initialValues) const { SQPSimpleState state(initialValues); - while (!state.converged) { + state.duals = initializeDuals(); + + while (!state.converged && state.iterations < 100) { state = iterate(state); } - return std::make_pair(initialValues, VectorValues()); + return std::make_pair(state.values, state.duals); } @@ -177,17 +207,38 @@ using namespace std; using namespace gtsam::symbol_shorthand; using namespace gtsam; const double tol = 1e-10; + //****************************************************************************** -TEST(testSQPSimple, Problem2) { +// x + y - 1 = 0 +class ConstraintProblem1 : public NonlinearConstraint2 { + typedef NonlinearConstraint2 Base; + +public: + ConstraintProblem1(Key xK, Key yK, Key dualKey) : Base(xK, yK, dualKey, 1) {} + + // x + y - 1 + Vector evaluateError(const double& x, const double& y, + boost::optional H1 = boost::none, boost::optional H2 = + boost::none) const { + if (H1) *H1 = eye(1); + if (H2) *H2 = eye(1); + return (Vector(1) << x + y - 1.0).finished(); + } +}; + +TEST_DISABLED(testSQPSimple, QPProblem) { + const Key dualKey = 0; + // Simple quadratic cost: x1^2 + x2^2 // Note the Hessian encodes: // 0.5*x1'*G11*x1 + x1'*G12*x2 + 0.5*x2'*G22*x2 - x1'*g1 - x2'*g2 + 0.5*f // Hence here we have G11 = 2, G12 = 0, G22 = 2, g1 = 0, g2 = 0, f = 0 - HessianFactor hf(X(1), X(2), 2.0 * ones(1,1), zero(1), zero(1), + HessianFactor hf(X(1), Y(1), 2.0 * ones(1,1), zero(1), zero(1), 2*ones(1,1), zero(1) , 0); LinearEqualityFactorGraph equalities; - equalities.push_back(LinearEquality(X(1), ones(1), X(2), ones(1), -1*ones(1), 0)); // x + y - 1 = 0 + LinearEquality linearConstraint(X(1), ones(1), Y(1), ones(1), 1*ones(1), dualKey); // x + y - 1 = 0 + equalities.push_back(linearConstraint); // Compare against QP QP qp; @@ -199,77 +250,114 @@ TEST(testSQPSimple, Problem2) { // create initial values for optimization VectorValues initialVectorValues; initialVectorValues.insert(X(1), zero(1)); - initialVectorValues.insert(X(2), zero(1)); + initialVectorValues.insert(Y(1), ones(1)); VectorValues expectedSolution = qpSolver.optimize(initialVectorValues).first; - cout<<"expectedSolution.at(X(1))[0]: "<(X(1)), 1e-100); + DOUBLES_EQUAL(expectedSolution.at(Y(1))[0], actualValues.at(Y(1)), 1e-100); + } +//****************************************************************************** +class CircleConstraint : public NonlinearConstraint2 { + typedef NonlinearConstraint2 Base; +public: + CircleConstraint(Key xK, Key yK, Key dualKey) : Base(xK, yK, dualKey, 1) {} + + Vector evaluateError(const double& x, const double& y, + boost::optional H1 = boost::none, boost::optional H2 = + boost::none) const { + if (H1) *H1 = (Matrix(1,1) << 2*(x-1)).finished(); + if (H2) *H2 = (Matrix(1,1) << 2*y).finished(); + return (Vector(1) << (x-1)*(x-1) + y*y - 0.25).finished(); + } + + void evaluateHessians(const double& x, const double& y, + std::vector& G11, std::vector& G12, + std::vector& G22) const { + G11.push_back((Matrix(1,1) << 2).finished()); + G12.push_back((Matrix(1,1) << 0).finished()); + G22.push_back((Matrix(1,1) << 2).finished()); + } + +}; + +TEST_UNSAFE(testSQPSimple, quadraticCostNonlinearConstraint) { + const Key dualKey = 0; + + //Instantiate NLP + NLP nlp; + + // Simple quadratic cost: x1^2 + x2^2 + // Note the Hessian encodes: + // 0.5*x1'*G11*x1 + x1'*G12*x2 + 0.5*x2'*G22*x2 - x1'*g1 - x2'*g2 + 0.5*f + // Hence here we have G11 = 2, G12 = 0, G22 = 2, g1 = 0, g2 = 0, f = 0 + Values linPoint; + linPoint.insert(X(1), zero(1)); + linPoint.insert(Y(1), zero(1)); + HessianFactor hf(X(1), Y(1), 2.0 * ones(1,1), zero(1), zero(1), + 2*ones(1,1), zero(1) , 0); + nlp.cost.add(LinearContainerFactor(hf, linPoint)); // wrap it using linearcontainerfactor + nlp.nonlinearEqualities.add(CircleConstraint(X(1), Y(1), dualKey)); + + Values initialValues; + initialValues.insert(X(1), 0.0); + initialValues.insert(Y(1), 2.0); + + Values expectedSolution; + expectedSolution.insert(X(1), 0.5); + expectedSolution.insert(Y(1), 0.0); + + // Instantiate SQP + SQPSimple sqpSimple(nlp); + Values actualSolution = sqpSimple.optimize(initialValues).first; + + CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); + actualSolution.print("actualSolution: "); +} + +//****************************************************************************** +// +//TEST_UNSAFE(testSQPSimple, poseOnALine) { +// const Key dualKey = 0; +// +// //Instantiate NLP +// NLP nlp; +// nlp.cost.add(PriorFactor(X(1), Pose3(), noiseModel::Unit::Create(6))); // wrap it using linearcontainerfactor +// nlp.nonlinearEqualities.add(CircleConstraint(X(1), Y(1), dualKey)); +// +// Values initialValues; +// initialValues.insert(X(1), 0.0); +// initialValues.insert(Y(1), 0.0); +// +// Values expectedSolution; +// expectedSolution.insert(X(1), 0.5); +// expectedSolution.insert(Y(1), 0.0); +// +// // Instantiate SQP +// SQPSimple sqpSimple(nlp); +// Values actualSolution = sqpSimple.optimize(initialValues).first; +// +// CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); +// actualSolution.print("actualSolution: "); +//} + //****************************************************************************** int main() { TestResult tr; return TestRegistry::runAllTests(tr); } //****************************************************************************** - -//TEST(testSQPSimple, Problem1 ) { -// -// // build a quadratic Objective function x1^2 - x1*x2 + x2^2 - 3*x1 + 5 -// // Note the Hessian encodes: -// // 0.5*x1'*G11*x1 + x1'*G12*x2 + 0.5*x2'*G22*x2 - x1'*g1 - x2'*g2 + 0.5*f -// // Hence, we have G11=2, G12 = -1, g1 = +3, G22 = 2, g2 = 0, f = 10 -// HessianFactor lf(X(1), X(2), 2.0 * ones(1, 1), -ones(1, 1), 3.0 * ones(1), -// 2.0 * ones(1, 1), zero(1), 10.0); -// -// // build linear inequalities -// LinearInequalityFactorGraph inequalities; -// inequalities.push_back( -// LinearInequality(X(1), ones(1, 1), X(2), ones(1, 1), 2, 0)); // x1 + x2 <= 2 --> x1 + x2 -2 <= 0, --> b=2 -// inequalities.push_back(LinearInequality(X(1), -ones(1, 1), 0, 1)); // -x1 <= 0 -// inequalities.push_back(LinearInequality(X(2), -ones(1, 1), 0, 2)); // -x2 <= 0 -// inequalities.push_back(LinearInequality(X(1), ones(1, 1), 1.5, 3)); // x1 <= 3/2 -// -// // Compare against a QP -// QP qp; -// qp.cost.add(lf); -// qp.inequalities = inequalities; -// -// // instantiate QPsolver -// QPSolver qpSolver(qp); -// // create initial values for optimization -// VectorValues initialVectorValues; -// initialVectorValues.insert(X(1), zero(1)); -// initialVectorValues.insert(X(2), zero(1)); -// VectorValues expectedSolution = qpSolver.optimize(initialVectorValues).first; -// -// -// NonlinearEqualityFactorGraph linearEqualities; -// NonlinearEqualityFactorGraph nonlinearEqualities; -// nonlinearEqualities.push_back(NonlinearEquality(X(1), ones(1, 1), X(2), ones(1, 1), 2, 0))); -// -// NLP nlp; -// nlp.cost.add(lf); -// nlp.linearEqualities.push_back(NonlinearEqualityFactorGraph()); -// // instantiate QPsolver -// SQPSimple sqpSolver(nlp); -// // create initial values for optimization -// Values initialValues; -// initialValues.insert(X(1), zero(1)); -// initialValues.insert(X(2), zero(1)); -// -// std::pair actualSolution = sqpSolver.optimize(initialValues); -// -// DOUBLES_EQUAL(expectedSolution.at(X(1))[0], actualSolution.at(X(1)), -// tol); -// DOUBLES_EQUAL(expectedSolution.at(X(2))[0], actualSolution.at(X(2)), -// tol); -//} From ecc87bdb2b7f27e8ebaa928c367a1de3c2bd8bd9 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Fri, 19 Dec 2014 19:38:10 -0500 Subject: [PATCH 009/130] Added test with pose and a line equality constraint. Works but hessian is incorrect. So basically using nonlinearequality vs linearequality makes no difference. --- .../nonlinear/tests/testSQPSimple.cpp | 85 ++++++++++++------- 1 file changed, 55 insertions(+), 30 deletions(-) diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index 16bc17b03..aa8de00e8 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -203,6 +203,10 @@ public: }; } +#include +#include +#include + using namespace std; using namespace gtsam::symbol_shorthand; using namespace gtsam; @@ -226,7 +230,7 @@ public: } }; -TEST_DISABLED(testSQPSimple, QPProblem) { +TEST(testSQPSimple, QPProblem) { const Key dualKey = 0; // Simple quadratic cost: x1^2 + x2^2 @@ -255,7 +259,10 @@ TEST_DISABLED(testSQPSimple, QPProblem) { //Instantiate NLP NLP nlp; - nlp.cost.add(LinearContainerFactor(hf)); // wrap it using linearcontainerfactor + Values linPoint; + linPoint.insert(X(1), zero(1)); + linPoint.insert(Y(1), zero(1)); + nlp.cost.add(LinearContainerFactor(hf, linPoint)); // wrap it using linearcontainerfactor nlp.linearEqualities.add(ConstraintProblem1(X(1), Y(1), dualKey)); Values initialValues; @@ -301,7 +308,7 @@ TEST_UNSAFE(testSQPSimple, quadraticCostNonlinearConstraint) { //Instantiate NLP NLP nlp; - // Simple quadratic cost: x1^2 + x2^2 + // Simple quadratic cost: x1^2 + x2^2 +1000 // Note the Hessian encodes: // 0.5*x1'*G11*x1 + x1'*G12*x2 + 0.5*x2'*G22*x2 - x1'*g1 - x2'*g2 + 0.5*f // Hence here we have G11 = 2, G12 = 0, G22 = 2, g1 = 0, g2 = 0, f = 0 @@ -309,13 +316,13 @@ TEST_UNSAFE(testSQPSimple, quadraticCostNonlinearConstraint) { linPoint.insert(X(1), zero(1)); linPoint.insert(Y(1), zero(1)); HessianFactor hf(X(1), Y(1), 2.0 * ones(1,1), zero(1), zero(1), - 2*ones(1,1), zero(1) , 0); + 2*ones(1,1), zero(1) , 1000); nlp.cost.add(LinearContainerFactor(hf, linPoint)); // wrap it using linearcontainerfactor nlp.nonlinearEqualities.add(CircleConstraint(X(1), Y(1), dualKey)); Values initialValues; - initialValues.insert(X(1), 0.0); - initialValues.insert(Y(1), 2.0); + initialValues.insert(X(1), 4.0); + initialValues.insert(Y(1), 10.0); Values expectedSolution; expectedSolution.insert(X(1), 0.5); @@ -330,30 +337,48 @@ TEST_UNSAFE(testSQPSimple, quadraticCostNonlinearConstraint) { } //****************************************************************************** -// -//TEST_UNSAFE(testSQPSimple, poseOnALine) { -// const Key dualKey = 0; -// -// //Instantiate NLP -// NLP nlp; -// nlp.cost.add(PriorFactor(X(1), Pose3(), noiseModel::Unit::Create(6))); // wrap it using linearcontainerfactor -// nlp.nonlinearEqualities.add(CircleConstraint(X(1), Y(1), dualKey)); -// -// Values initialValues; -// initialValues.insert(X(1), 0.0); -// initialValues.insert(Y(1), 0.0); -// -// Values expectedSolution; -// expectedSolution.insert(X(1), 0.5); -// expectedSolution.insert(Y(1), 0.0); -// -// // Instantiate SQP -// SQPSimple sqpSimple(nlp); -// Values actualSolution = sqpSimple.optimize(initialValues).first; -// -// CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); -// actualSolution.print("actualSolution: "); -//} +class LineConstraintX : public NonlinearConstraint1 { + typedef NonlinearConstraint1 Base; +public: + LineConstraintX(Key key, Key dualKey) : Base(key, dualKey, 1) { + } + + double computeError(const Pose3& pose) const { + return pose.x(); + } + + Vector evaluateError(const Pose3& pose, boost::optional H = boost::none) const { + if (H) + *H = (Matrix(1,6) << zeros(1,3), pose.rotation().matrix().row(0)).finished(); + return (Vector(1) << pose.x()).finished(); + } +}; +TEST_UNSAFE(testSQPSimple, poseOnALine) { + const Key dualKey = 0; + + + //Instantiate NLP + NLP nlp; + nlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(1, 0, 0)), noiseModel::Unit::Create(6))); + LineConstraintX constraint(X(1), dualKey); + nlp.nonlinearEqualities.add(constraint); + + Values initialValues; + initialValues.insert(X(1), Pose3(Rot3::ypr(0.3, 0.2, 0.3), Point3(1,0,0))); + + Values expectedSolution; + expectedSolution.insert(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3())); + + // Instantiate SQP + SQPSimple sqpSimple(nlp); + Values actualSolution = sqpSimple.optimize(initialValues).first; + + CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); + actualSolution.print("actualSolution: "); + Pose3 pose(Rot3::ypr(0.1, 0.2, 0.3), Point3()); + Matrix hessian = numericalHessian(boost::bind(&LineConstraintX::computeError, constraint, _1), pose, 1e-2); + cout << "hessian: \n" << hessian << endl; +} //****************************************************************************** int main() { From fd461a1c15552eec52273b9baeeb78078db8770a Mon Sep 17 00:00:00 2001 From: krunalchande Date: Mon, 22 Dec 2014 14:41:22 -0500 Subject: [PATCH 010/130] [unfinished] prototyping inequality SQP with Luca. --- gtsam_unstable/linear/LinearInequality.h | 13 +- .../nonlinear/NonlinearConstraint.h | 3 +- .../nonlinear/NonlinearInequality.h | 93 ++++++++++++++ .../nonlinear/tests/testSQPSimple.cpp | 121 +++++++++++++++++- 4 files changed, 225 insertions(+), 5 deletions(-) create mode 100644 gtsam_unstable/nonlinear/NonlinearInequality.h diff --git a/gtsam_unstable/linear/LinearInequality.h b/gtsam_unstable/linear/LinearInequality.h index 7c62c3d54..6aa03eaef 100644 --- a/gtsam_unstable/linear/LinearInequality.h +++ b/gtsam_unstable/linear/LinearInequality.h @@ -44,12 +44,23 @@ public: Base(), active_(true) { } - /** Conversion from HessianFactor (does Cholesky to obtain Jacobian matrix) */ + /** Conversion from HessianFactor */ explicit LinearInequality(const HessianFactor& hf) { throw std::runtime_error( "Cannot convert HessianFactor to LinearInequality"); } + /** Conversion from JacobianFactor */ + explicit LinearInequality(const JacobianFactor& jf) : Base(jf), dualKey_(dualKey), active_(true) { + if (!jf.isConstrained()) { + throw std::runtime_error("Cannot convert an unconstrained JacobianFactor to LinearEquality"); + } + + if (jf.get_model()->dim() != 1) { + throw std::runtime_error("Only support single-valued inequality factor!"); + } + } + /** Construct unary factor */ LinearInequality(Key i1, const RowVector& A1, double b, Key dualKey) : Base(i1, A1, (Vector(1) << b).finished(), noiseModel::Constrained::All(1)), dualKey_( diff --git a/gtsam_unstable/nonlinear/NonlinearConstraint.h b/gtsam_unstable/nonlinear/NonlinearConstraint.h index ce01f89e8..f24517983 100644 --- a/gtsam_unstable/nonlinear/NonlinearConstraint.h +++ b/gtsam_unstable/nonlinear/NonlinearConstraint.h @@ -13,11 +13,12 @@ #include #include #include -//#include "DualKeyGenerator.h" namespace gtsam { class NonlinearConstraint { + +protected: Key dualKey_; public: diff --git a/gtsam_unstable/nonlinear/NonlinearInequality.h b/gtsam_unstable/nonlinear/NonlinearInequality.h new file mode 100644 index 000000000..b9a9e418a --- /dev/null +++ b/gtsam_unstable/nonlinear/NonlinearInequality.h @@ -0,0 +1,93 @@ +/** + * @file NonlinearConstraint.h + * @brief + * @author Duy-Nguyen Ta + * @date Sep 30, 2013 + */ + +#pragma once + +#include + +namespace gtsam { + +class NonlinearInequality : public NonlinearConstraint { + bool active_; + + typedef NonlinearConstraint Base; + +public: + typedef boost::shared_ptr shared_ptr; + +public: + /// Construct with dual key + NonlinearInequality(Key dualKey) : Base(dualKey), active_(true) {} + + /** + * compute the HessianFactor of the (-dual * constraintHessian) for the qp subproblem's objective function + */ + virtual GaussianFactor::shared_ptr multipliedHessian(const Values& x, + const VectorValues& duals) const = 0; +}; + +/* ************************************************************************* */ +/** A convenient base class for creating a nonlinear equality constraint with 1 + * variables. To derive from this class, implement evaluateError(). */ +template +class NonlinearInequality1: public NonlinearConstraint1, public NonlinearInequality { + +public: + + // typedefs for value types pulled from keys + typedef VALUE X; + +protected: + + typedef NonlinearConstraint1 Base; + typedef NonlinearInequality1 This; + +private: + static const int X1Dim = traits::dimension::value; + +public: + + /** + * Default Constructor for I/O + */ + NonlinearInequality1() { + } + + /** + * Constructor + * @param j key of the variable + * @param constraintDim number of dimensions of the constraint error function + */ + NonlinearInequality1(Key key, Key dualKey, size_t constraintDim = 1) : + Base(noiseModel::Constrained::All(constraintDim), key), NonlinearConstraint(dualKey) { + } + + virtual ~NonlinearInequality1() { + } + + /** + * Override this method to finish implementing a binary factor. + * If any of the optional Matrix reference arguments are specified, it should compute + * both the function evaluation and its derivative(s) in X1 (and/or X2). + */ + virtual Vector + evaluateError(const X&, boost::optional H1 = boost::none) const { + return (Vector(1) << computeError(X, H1)); + } + + /** + * Override this method to finish implementing a binary factor. + * If any of the optional Matrix reference arguments are specified, it should compute + * both the function evaluation and its derivative(s) in X1 (and/or X2). + */ + virtual double + computeError(const X&, boost::optional H1 = boost::none) const = 0; + +}; +// \class NonlinearConstraint1 + +} /* namespace gtsam */ diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index aa8de00e8..08bcfb614 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -80,11 +80,76 @@ public: } }; +class NonlinearInequalityFactorGraph : public FactorGraph { + +public: + /// default constructor + NonlinearInequalityFactorGraph() { + } + + /// linearize to a LinearEqualityFactorGraph + LinearInequalityFactorGraph::shared_ptr linearize( + const Values& linearizationPoint) const { + LinearInequalityFactorGraph::shared_ptr linearGraph( + new LinearInequalityFactorGraph()); + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + JacobianFactor::shared_ptr jacobian = boost::dynamic_pointer_cast( + factor->linearize(linearizationPoint)); + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + linearGraph->add(LinearInequality(*jacobian, constraint->dualKey())); + } + return linearGraph; + } + + /** + * Return true if the error is <= 0.0 + */ + bool checkFeasibility(const Values& values, double tol) const { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( + factor); + Vector error = noiseModelFactor->unwhitenedError(values); + // TODO: Do we need to check if it's active or not? + if (error[0] > tol) { + return false; + } + } + return true; + } + + /** + * Return true if the max absolute error all factors is less than a tolerance + */ + bool checkDualFeasibility(const VectorValues& duals, double tol) const { + BOOST_FOREACH(const Vector& dual, duals){ + if (dual[0] < 0.0) { + return false; + } + } + return true; + } + + /** + * Return true if the max absolute error all factors is less than a tolerance + */ + bool checkComplimentaryCondition(const Values& values, const VectorValues& duals, double tol) const { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( + factor); + Vector error = noiseModelFactor->unwhitenedError(values); + if (error[0] > 0.0) { + return false; + } + } + return true; + } +}; struct NLP { NonlinearFactorGraph cost; NonlinearEqualityFactorGraph linearEqualities; NonlinearEqualityFactorGraph nonlinearEqualities; + NonlinearInequalityFactorGraph linearInequalities; }; struct SQPSimpleState { @@ -117,14 +182,16 @@ public: /// Check if \nabla f(x) - \lambda * \nabla c(x) == 0 bool isDualFeasible(const VectorValues& delta) const { - return delta.vector().lpNorm() < errorTol; + return delta.vector().lpNorm() < errorTol + && nlp_.linearInequalities.checkDualFeasibility(errorTol); // return false; } /// Check if c(x) == 0 bool isPrimalFeasible(const SQPSimpleState& state) const { return nlp_.linearEqualities.checkFeasibility(state.values, errorTol) - && nlp_.nonlinearEqualities.checkFeasibility(state.values, errorTol); + && nlp_.nonlinearEqualities.checkFeasibility(state.values, errorTol) + && nlp_.linearInequalities.checkFeasibility(state.values, errorTol); } /// Check convergence @@ -147,6 +214,8 @@ public: qp.equalities.add(*nlp_.linearEqualities.linearize(state.values)); qp.equalities.add(*nlp_.nonlinearEqualities.linearize(state.values)); + qp.inequalities.add(*nlp_.linearInequalities.linearize(state.values)); + if (debug) qp.print("QP subproblem:"); @@ -206,6 +275,7 @@ public: #include #include #include +#include using namespace std; using namespace gtsam::symbol_shorthand; @@ -353,7 +423,8 @@ public: return (Vector(1) << pose.x()).finished(); } }; -TEST_UNSAFE(testSQPSimple, poseOnALine) { + +TEST(testSQPSimple, poseOnALine) { const Key dualKey = 0; @@ -380,6 +451,50 @@ TEST_UNSAFE(testSQPSimple, poseOnALine) { cout << "hessian: \n" << hessian << endl; } + +//****************************************************************************** +/** + * Inequality boundary constraint + * x <= bound + */ +class UpperBoundX : public NonlinearInequality1 { + typedef NonlinearInequality1 Base; + double bound_; +public: + UpperBoundX(Key key, double bound, Key dualKey) : Base(key, dualKey, 1), bound_(bound) { + } + + double computeError(const Pose3& pose, boost::optional H = boost::none) const { + if (H) + *H = (Matrix(1,6) << zeros(1,3), pose.rotation().matrix().row(0)).finished(); + return pose.x() - bound_; + } +}; + +TEST(testSQPSimple, poseOnALine) { + const Key dualKey = 0; + + + //Instantiate NLP + NLP nlp; + nlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(-1, 0, 0)), noiseModel::Unit::Create(6))); + UpperBoundX constraint(X(1), 0, dualKey); + nlp.nonlinearInequalities.add(constraint); + + Values initialValues; + initialValues.insert(X(1), Pose3(Rot3::ypr(0.3, 0.2, 0.3), Point3(-1,0,0))); + + Values expectedSolution; + expectedSolution.insert(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3())); + + // Instantiate SQP + SQPSimple sqpSimple(nlp); + Values actualSolution = sqpSimple.optimize(initialValues).first; + + CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); + actualSolution.print("actualSolution: "); +} + //****************************************************************************** int main() { TestResult tr; From 4f92492e34e4eb60eca3fabbc28e11226f816839 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Mon, 22 Dec 2014 15:47:28 -0500 Subject: [PATCH 011/130] Reapply hacks in EliminatePreferCholesky to deal with negative definite hessians obtained from multiplying dual variables with the hessians of nonlinear constraints needed for SQP. --- gtsam/linear/GaussianFactorGraph.cpp | 23 ++++++ gtsam/linear/GaussianFactorGraph.h | 6 ++ gtsam/linear/HessianFactor.cpp | 106 ++++++++++++++++++++++++++- 3 files changed, 132 insertions(+), 3 deletions(-) diff --git a/gtsam/linear/GaussianFactorGraph.cpp b/gtsam/linear/GaussianFactorGraph.cpp index 6953d2969..1d6066c1c 100644 --- a/gtsam/linear/GaussianFactorGraph.cpp +++ b/gtsam/linear/GaussianFactorGraph.cpp @@ -384,6 +384,29 @@ namespace gtsam { return false; } + /* ************************************************************************* */ + boost::tuple GaussianFactorGraph::splitConstraints() const { + typedef HessianFactor H; + typedef JacobianFactor J; + + GaussianFactorGraph hessians, jacobians, constraints; + BOOST_FOREACH(const GaussianFactor::shared_ptr& factor, *this) { + H::shared_ptr hessian(boost::dynamic_pointer_cast(factor)); + if (hessian) + hessians.push_back(factor); + else { + J::shared_ptr jacobian(boost::dynamic_pointer_cast(factor)); + if (jacobian && jacobian->get_model() && jacobian->get_model()->isConstrained()) { + constraints.push_back(jacobian); + } + else { + jacobians.push_back(factor); + } + } + } + return boost::make_tuple(hessians, jacobians, constraints); + } + /* ************************************************************************* */ // x += alpha*A'*e void GaussianFactorGraph::transposeMultiplyAdd(double alpha, const Errors& e, diff --git a/gtsam/linear/GaussianFactorGraph.h b/gtsam/linear/GaussianFactorGraph.h index 811c0f8b0..067a42d03 100644 --- a/gtsam/linear/GaussianFactorGraph.h +++ b/gtsam/linear/GaussianFactorGraph.h @@ -316,6 +316,12 @@ namespace gtsam { /** In-place version e <- A*x that takes an iterator. */ void multiplyInPlace(const VectorValues& x, const Errors::iterator& e) const; + /** + * Split constraints and unconstrained factors into two different graphs + * @return a pair of graphs + */ + boost::tuple splitConstraints() const; + /// @} private: diff --git a/gtsam/linear/HessianFactor.cpp b/gtsam/linear/HessianFactor.cpp index f2bebcab9..755f9c313 100644 --- a/gtsam/linear/HessianFactor.cpp +++ b/gtsam/linear/HessianFactor.cpp @@ -616,10 +616,110 @@ EliminatePreferCholesky(const GaussianFactorGraph& factors, const Ordering& keys // all factors to JacobianFactors. Otherwise, we can convert all factors // to HessianFactors. This is because QR can handle constrained noise // models but Cholesky cannot. - if (hasConstraints(factors)) - return EliminateQR(factors, keys); - else + + /* Currently, when eliminating a constrained variable, EliminatePreferCholesky + * converts every other factors to JacobianFactor before doing the special QR + * factorization for constrained variables. Unfortunately, after a constrained + * nonlinear graph is linearized, new hessian factors from constraints, multiplied + * with the dual variable (-lambda*\hessian{c} terms in the Lagrangian objective + * function), might become negative definite, thus cannot be converted to JacobianFactors. + * + * Following EliminateCholesky, this version of EliminatePreferCholesky for + * constrained var gathers all unconstrained factors into a big joint HessianFactor + * before converting it into a JacobianFactor to be eliminiated by QR together with + * the other constrained factors. + * + * Of course, this might not solve the non-positive-definite problem entirely, + * because (1) the original hessian factors might be non-positive definite + * and (2) large strange value of lambdas might cause the joint factor non-positive + * definite [is this true?]. But at least, this will help in typical cases. + */ + GaussianFactorGraph hessians, jacobians, constraints; +// factors.print("factors: "); + boost::tie(hessians, jacobians, constraints) = factors.splitConstraints(); +// keys.print("Frontal variables to eliminate: "); +// hessians.print("Hessians: "); +// jacobians.print("Jacobians: "); +// constraints.print("Constraints: "); + + bool hasHessians = hessians.size() > 0; + + // Add all jacobians to gather as much info as we can + hessians.push_back(jacobians); + + if (constraints.size()>0) { +// // Build joint factor +// HessianFactor::shared_ptr jointFactor; +// try { +// jointFactor = boost::make_shared(hessians, Scatter(factors, keys)); +// } catch(std::invalid_argument&) { +// throw InvalidDenseElimination( +// "EliminateCholesky was called with a request to eliminate variables that are not\n" +// "involved in the provided factors."); +// } +// constraints.push_back(jointFactor); +// return EliminateQR(constraints, keys); + + // If there are hessian factors, turn them into conditional + // by doing partial elimination, then use those jacobians when eliminating the constraints + GaussianFactor::shared_ptr unconstrainedNewFactor; + if (hessians.size() > 0) { + bool hasSeparator = false; + GaussianFactorGraph::Keys unconstrainedKeys = hessians.keys(); + BOOST_FOREACH(Key key, unconstrainedKeys) { + if (find(keys.begin(), keys.end(), key) == keys.end()) { + hasSeparator = true; + break; + } + } + + if (hasSeparator) { + // find frontal keys in the unconstrained factor to eliminate + Ordering subkeys; + BOOST_FOREACH(Key key, keys) { + if (unconstrainedKeys.exists(key)) + subkeys.push_back(key); + } + GaussianConditional::shared_ptr unconstrainedConditional; + boost::tie(unconstrainedConditional, unconstrainedNewFactor) + = EliminateCholesky(hessians, subkeys); + constraints.push_back(unconstrainedConditional); + } + else { + if (hasHessians) { + HessianFactor::shared_ptr jointFactor = boost::make_shared< + HessianFactor>(hessians, Scatter(factors, keys)); + constraints.push_back(jointFactor); + } + else { + constraints.push_back(hessians); + } + } + } + + // Now eliminate the constraints + GaussianConditional::shared_ptr constrainedConditional; + GaussianFactor::shared_ptr constrainedNewFactor; + boost::tie(constrainedConditional, constrainedNewFactor) = EliminateQR( + constraints, keys); +// constraints.print("constraints: "); +// constrainedConditional->print("constrainedConditional: "); +// constrainedNewFactor->print("constrainedNewFactor: "); + + if (unconstrainedNewFactor) { + GaussianFactorGraph newFactors; + newFactors.push_back(unconstrainedNewFactor); + newFactors.push_back(constrainedNewFactor); +// newFactors.print("newFactors: "); + HessianFactor::shared_ptr newFactor(new HessianFactor(newFactors)); + return make_pair(constrainedConditional, newFactor); + } else { + return make_pair(constrainedConditional, constrainedNewFactor); + } + } + else { return EliminateCholesky(factors, keys); + } } } // gtsam From cc0e5cd3ca0f29b47398b8bdb3e65556275066f8 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Mon, 22 Dec 2014 18:20:44 -0500 Subject: [PATCH 012/130] Working nonlinear inequality constraints with unit tests. --- gtsam_unstable/linear/LinearInequality.h | 2 +- .../nonlinear/NonlinearConstraint.h | 40 +-- .../nonlinear/NonlinearInequality.h | 95 +++++-- .../nonlinear/tests/testSQPSimple.cpp | 233 ++++++++++++++---- 4 files changed, 272 insertions(+), 98 deletions(-) diff --git a/gtsam_unstable/linear/LinearInequality.h b/gtsam_unstable/linear/LinearInequality.h index 6aa03eaef..7eac21903 100644 --- a/gtsam_unstable/linear/LinearInequality.h +++ b/gtsam_unstable/linear/LinearInequality.h @@ -51,7 +51,7 @@ public: } /** Conversion from JacobianFactor */ - explicit LinearInequality(const JacobianFactor& jf) : Base(jf), dualKey_(dualKey), active_(true) { + explicit LinearInequality(const JacobianFactor& jf, Key dualKey) : Base(jf), dualKey_(dualKey), active_(true) { if (!jf.isConstrained()) { throw std::runtime_error("Cannot convert an unconstrained JacobianFactor to LinearEquality"); } diff --git a/gtsam_unstable/nonlinear/NonlinearConstraint.h b/gtsam_unstable/nonlinear/NonlinearConstraint.h index f24517983..af52ea178 100644 --- a/gtsam_unstable/nonlinear/NonlinearConstraint.h +++ b/gtsam_unstable/nonlinear/NonlinearConstraint.h @@ -70,7 +70,7 @@ public: * @param constraintDim number of dimensions of the constraint error function */ NonlinearConstraint1(Key key, Key dualKey, size_t constraintDim = 1) : - Base(noiseModel::Constrained::All(constraintDim), key), NonlinearConstraint(dualKey) { + Base(noiseModel::Constrained::All(constraintDim), key), NonlinearConstraint(dualKey) { } virtual ~NonlinearConstraint1() { @@ -108,8 +108,12 @@ public: lG11sum += -lambda[i] * G11[i]; } - return boost::make_shared(Base::key(), lG11sum, - zero(X1Dim), 100.0); + + std::cout << "lG11sum: " << lG11sum << std::endl; + HessianFactor::shared_ptr hf(new HessianFactor(Base::key(), lG11sum, + zero(X1Dim), 100.0)); + hf->print("HessianFactor: "); + return hf; } /** evaluate Hessians for lambda factors */ @@ -149,8 +153,8 @@ private: template void serialize(ARCHIVE & ar, const unsigned int version) { ar - & boost::serialization::make_nvp("NoiseModelFactor", - boost::serialization::base_object(*this)); + & boost::serialization::make_nvp("NoiseModelFactor", + boost::serialization::base_object(*this)); } }; // \class NonlinearConstraint1 @@ -191,7 +195,7 @@ public: * @param constraintDim number of dimensions of the constraint error function */ NonlinearConstraint2(Key j1, Key j2, Key dualKey, size_t constraintDim = 1) : - Base(noiseModel::Constrained::All(constraintDim), j1, j2), NonlinearConstraint(dualKey) { + Base(noiseModel::Constrained::All(constraintDim), j1, j2), NonlinearConstraint(dualKey) { } virtual ~NonlinearConstraint2() { @@ -229,7 +233,7 @@ public: // Combine the Lagrange-multiplier part into this constraint factor Matrix lG11sum = zeros(G11[0].rows(), G11[0].cols()), lG12sum = zeros( G12[0].rows(), G12[0].cols()), lG22sum = zeros(G22[0].rows(), - G22[0].cols()); + G22[0].cols()); for (size_t i = 0; i < lambda.rows(); ++i) { lG11sum += -lambda[i] * G11[i]; lG12sum += -lambda[i] * G12[i]; @@ -249,7 +253,7 @@ public: boost::function vecH1( boost::bind(&This::vectorizeH1t, this, _1, _2)), vecH2( - boost::bind(&This::vectorizeH2t, this, _1, _2)); + boost::bind(&This::vectorizeH2t, this, _1, _2)); Matrix G11all = numericalDerivative21(vecH1, x1, x2, 1e-5); Matrix G12all = numericalDerivative22(vecH1, x1, x2, 1e-5); @@ -302,8 +306,8 @@ private: template void serialize(ARCHIVE & ar, const unsigned int version) { ar - & boost::serialization::make_nvp("NoiseModelFactor", - boost::serialization::base_object(*this)); + & boost::serialization::make_nvp("NoiseModelFactor", + boost::serialization::base_object(*this)); } }; // \class NonlinearConstraint2 @@ -346,7 +350,7 @@ public: * @param constraintDim number of dimensions of the constraint error function */ NonlinearConstraint3(Key j1, Key j2, Key j3, Key dualKey, size_t constraintDim = 1) : - Base(noiseModel::Constrained::All(constraintDim), j1, j2, j3), NonlinearConstraint(dualKey) { + Base(noiseModel::Constrained::All(constraintDim), j1, j2, j3), NonlinearConstraint(dualKey) { } virtual ~NonlinearConstraint3() { @@ -387,9 +391,9 @@ public: // Combine the Lagrange-multiplier part into this constraint factor Matrix lG11sum = zeros(G11[0].rows(), G11[0].cols()), lG12sum = zeros( G12[0].rows(), G12[0].cols()), lG13sum = zeros(G13[0].rows(), - G13[0].cols()), lG22sum = zeros(G22[0].rows(), G22[0].cols()), lG23sum = - zeros(G23[0].rows(), G23[0].cols()), lG33sum = zeros(G33[0].rows(), - G33[0].cols()); + G13[0].cols()), lG22sum = zeros(G22[0].rows(), G22[0].cols()), lG23sum = + zeros(G23[0].rows(), G23[0].cols()), lG33sum = zeros(G33[0].rows(), + G33[0].cols()); for (size_t i = 0; i < lambda.rows(); ++i) { lG11sum += -lambda[i] * G11[i]; lG12sum += -lambda[i] * G12[i]; @@ -467,8 +471,8 @@ public: boost::function vecH1( boost::bind(&This::vectorizeH1t, this, _1, _2, _3)), vecH2( - boost::bind(&This::vectorizeH2t, this, _1, _2, _3)), vecH3( - boost::bind(&This::vectorizeH3t, this, _1, _2, _3)); + boost::bind(&This::vectorizeH2t, this, _1, _2, _3)), vecH3( + boost::bind(&This::vectorizeH3t, this, _1, _2, _3)); Matrix G11all = numericalDerivative31(vecH1, x1, x2, x3, 1e-5); Matrix G12all = numericalDerivative32(vecH1, x1, x2, x3, 1e-5); @@ -549,8 +553,8 @@ private: template void serialize(ARCHIVE & ar, const unsigned int version) { ar - & boost::serialization::make_nvp("NoiseModelFactor", - boost::serialization::base_object(*this)); + & boost::serialization::make_nvp("NoiseModelFactor", + boost::serialization::base_object(*this)); } }; // \class NonlinearConstraint3 diff --git a/gtsam_unstable/nonlinear/NonlinearInequality.h b/gtsam_unstable/nonlinear/NonlinearInequality.h index b9a9e418a..8d41b3ec0 100644 --- a/gtsam_unstable/nonlinear/NonlinearInequality.h +++ b/gtsam_unstable/nonlinear/NonlinearInequality.h @@ -11,30 +11,11 @@ namespace gtsam { -class NonlinearInequality : public NonlinearConstraint { - bool active_; - - typedef NonlinearConstraint Base; - -public: - typedef boost::shared_ptr shared_ptr; - -public: - /// Construct with dual key - NonlinearInequality(Key dualKey) : Base(dualKey), active_(true) {} - - /** - * compute the HessianFactor of the (-dual * constraintHessian) for the qp subproblem's objective function - */ - virtual GaussianFactor::shared_ptr multipliedHessian(const Values& x, - const VectorValues& duals) const = 0; -}; - /* ************************************************************************* */ /** A convenient base class for creating a nonlinear equality constraint with 1 * variables. To derive from this class, implement evaluateError(). */ template -class NonlinearInequality1: public NonlinearConstraint1, public NonlinearInequality { +class NonlinearInequality1: public NonlinearConstraint1 { public: @@ -62,8 +43,8 @@ public: * @param j key of the variable * @param constraintDim number of dimensions of the constraint error function */ - NonlinearInequality1(Key key, Key dualKey, size_t constraintDim = 1) : - Base(noiseModel::Constrained::All(constraintDim), key), NonlinearConstraint(dualKey) { + NonlinearInequality1(Key key, Key dualKey) : + Base(key, dualKey, 1) { } virtual ~NonlinearInequality1() { @@ -74,9 +55,63 @@ public: * If any of the optional Matrix reference arguments are specified, it should compute * both the function evaluation and its derivative(s) in X1 (and/or X2). */ + virtual double + computeError(const X&, boost::optional H1 = boost::none) const = 0; + + /** predefine evaluateError to return a 1-dimension vector */ virtual Vector - evaluateError(const X&, boost::optional H1 = boost::none) const { - return (Vector(1) << computeError(X, H1)); + evaluateError(const X& x, boost::optional H1 = boost::none) const { + return (Vector(1) << computeError(x, H1)).finished(); + } +// +// virtual GaussianFactor::shared_ptr multipliedHessian(const Values& x, +// const VectorValues& duals) const { +// return Base::multipliedHessian(x, duals); +// } + +}; +// \class NonlinearConstraint1 + +/* ************************************************************************* */ +/** A convenient base class for creating your own NonlinearConstraint with 2 + * variables. To derive from this class, implement evaluateError(). */ +template +class NonlinearInequality2: public NonlinearConstraint2 { + +public: + + // typedefs for value types pulled from keys + typedef VALUE1 X1; + typedef VALUE2 X2; + +protected: + + typedef NonlinearConstraint2 Base; + typedef NonlinearInequality2 This; + +private: + static const int X1Dim = traits::dimension::value; + static const int X2Dim = traits::dimension::value; + +public: + + /** + * Default Constructor for I/O + */ + NonlinearInequality2() { + } + + /** + * Constructor + * @param j1 key of the first variable + * @param j2 key of the second variable + * @param constraintDim number of dimensions of the constraint error function + */ + NonlinearInequality2(Key j1, Key j2, Key dualKey) : + Base(j1, j2, 1) { + } + + virtual ~NonlinearInequality2() { } /** @@ -85,9 +120,17 @@ public: * both the function evaluation and its derivative(s) in X1 (and/or X2). */ virtual double - computeError(const X&, boost::optional H1 = boost::none) const = 0; + computeError(const X1&, const X2&, boost::optional H1 = boost::none, + boost::optional H2 = boost::none) const = 0; + /** predefine evaluateError to return a 1-dimension vector */ + virtual Vector + evaluateError(const X1& x1, const X2& x2, boost::optional H1 = boost::none, + boost::optional H2 = boost::none) const { + return (Vector(1) << computeError(x1, x2, H1, H2)).finished(); + } }; -// \class NonlinearConstraint1 +// \class NonlinearConstraint2 + } /* namespace gtsam */ diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index 08bcfb614..d83312ee1 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -102,47 +103,31 @@ public: } /** - * Return true if the error is <= 0.0 + * Return true if the all errors are <= 0.0 */ - bool checkFeasibility(const Values& values, double tol) const { + bool checkFeasibilityAndComplimentary(const Values& values, const VectorValues& duals, double tol) const { BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( factor); Vector error = noiseModelFactor->unwhitenedError(values); - // TODO: Do we need to check if it's active or not? + + // Primal feasibility condition: all constraints need to be <= 0.0 if (error[0] > tol) { return false; } - } - return true; - } - /** - * Return true if the max absolute error all factors is less than a tolerance - */ - bool checkDualFeasibility(const VectorValues& duals, double tol) const { - BOOST_FOREACH(const Vector& dual, duals){ - if (dual[0] < 0.0) { - return false; - } - } - return true; - } - - /** - * Return true if the max absolute error all factors is less than a tolerance - */ - bool checkComplimentaryCondition(const Values& values, const VectorValues& duals, double tol) const { - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ - NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( + // Complimentary condition: errors of active constraints need to be 0.0 + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast( factor); - Vector error = noiseModelFactor->unwhitenedError(values); - if (error[0] > 0.0) { + Key dualKey = constraint->dualKey(); + if (!duals.exists(dualKey)) continue; // if dualKey doesn't exist, it is an inactive constraint! + if (fabs(error[0]) > tol) // for active constraint, the error should be 0.0 return false; - } + } return true; } + }; struct NLP { @@ -181,22 +166,49 @@ public: } /// Check if \nabla f(x) - \lambda * \nabla c(x) == 0 - bool isDualFeasible(const VectorValues& delta) const { - return delta.vector().lpNorm() < errorTol - && nlp_.linearInequalities.checkDualFeasibility(errorTol); -// return false; + bool isStationary(const VectorValues& delta) const { + return delta.vector().lpNorm() < errorTol; } - /// Check if c(x) == 0 + /// Check if c_E(x) == 0 bool isPrimalFeasible(const SQPSimpleState& state) const { return nlp_.linearEqualities.checkFeasibility(state.values, errorTol) - && nlp_.nonlinearEqualities.checkFeasibility(state.values, errorTol) - && nlp_.linearInequalities.checkFeasibility(state.values, errorTol); + && nlp_.nonlinearEqualities.checkFeasibility(state.values, errorTol); + } + + /** + * Dual variables of inequality constraints need to be >=0 + * For active inequalities, the dual needs to be > 0 + * For inactive inequalities, they need to be == 0. However, we don't compute + * dual variables for inactive constraints in the qp subproblem, so we don't care. + */ + bool isDualFeasible(const VectorValues& duals) const { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearInequalities) { + NonlinearConstraint::shared_ptr inequality = boost::dynamic_pointer_cast(factor); + Key dualKey = inequality->dualKey(); + if (!duals.exists(dualKey)) continue; // should be inactive constraint! + double dual = duals.at(dualKey)[0]; // because we only support single-valued inequalities + if (dual < 0.0) + return false; + } + return true; + } + + /** + * Check complimentary slackness condition: + * For all inequality constraints, + * dual * constraintError(primals) == 0. + * If the constraint is active, we need to check constraintError(primals) == 0, and ignore the dual + * If it is inactive, the dual should be 0, regardless of the error. However, we don't compute + * dual variables for inactive constraints in the QP subproblem, so we don't care. + */ + bool isComplementary(const SQPSimpleState& state) const { + return nlp_.linearInequalities.checkFeasibilityAndComplimentary(state.values, state.duals, errorTol); } /// Check convergence bool checkConvergence(const SQPSimpleState& state, const VectorValues& delta) const { - return isPrimalFeasible(state) & isDualFeasible(delta); + return isStationary(delta) && isPrimalFeasible(state) && isDualFeasible(state.duals) && isComplementary(state); } /** @@ -275,7 +287,6 @@ public: #include #include #include -#include using namespace std; using namespace gtsam::symbol_shorthand; @@ -403,7 +414,6 @@ TEST_UNSAFE(testSQPSimple, quadraticCostNonlinearConstraint) { Values actualSolution = sqpSimple.optimize(initialValues).first; CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); - actualSolution.print("actualSolution: "); } //****************************************************************************** @@ -417,6 +427,13 @@ public: return pose.x(); } + void evaluateHessians(const Pose3& pose, std::vector& G11) const { + Matrix G11all = Z_6x6; + Vector rT1 = pose.rotation().matrix().row(0); + G11all.block<3,3>(3,0) = skewSymmetric(rT1); + G11.push_back(G11all); + } + Vector evaluateError(const Pose3& pose, boost::optional H = boost::none) const { if (H) *H = (Matrix(1,6) << zeros(1,3), pose.rotation().matrix().row(0)).finished(); @@ -445,54 +462,164 @@ TEST(testSQPSimple, poseOnALine) { Values actualSolution = sqpSimple.optimize(initialValues).first; CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); - actualSolution.print("actualSolution: "); Pose3 pose(Rot3::ypr(0.1, 0.2, 0.3), Point3()); Matrix hessian = numericalHessian(boost::bind(&LineConstraintX::computeError, constraint, _1), pose, 1e-2); - cout << "hessian: \n" << hessian << endl; +} + +//****************************************************************************** +/// x + y - 1 <= 0 +class InequalityProblem1 : public NonlinearInequality2 { + typedef NonlinearInequality2 Base; +public: + InequalityProblem1(Key xK, Key yK, Key dualKey) : Base(xK, yK, dualKey) {} + + double computeError(const double& x, const double& y, + boost::optional H1 = boost::none, boost::optional H2 = + boost::none) const { + if (H1) *H1 = eye(1); + if (H2) *H2 = eye(1); + return x + y - 1.0; + } +}; + +TEST(testSQPSimple, inequalityConstraint) { + const Key dualKey = 0; + + // Simple quadratic cost: x^2 + y^2 + // Note the Hessian encodes: + // 0.5*x1'*G11*x1 + x1'*G12*x2 + 0.5*x2'*G22*x2 - x1'*g1 - x2'*g2 + 0.5*f + // Hence here we have G11 = 2, G12 = 0, G22 = 2, g1 = 0, g2 = 0, f = 0 + HessianFactor hf(X(1), Y(1), 2.0 * ones(1,1), zero(1), zero(1), + 2*ones(1,1), zero(1) , 0); + + LinearInequalityFactorGraph inequalities; + LinearInequality linearConstraint(X(1), ones(1), Y(1), ones(1), 1.0, dualKey); // x + y - 1 <= 0 + inequalities.push_back(linearConstraint); + + // Compare against QP + QP qp; + qp.cost.add(hf); + qp.inequalities = inequalities; + + // instantiate QPsolver + QPSolver qpSolver(qp); + // create initial values for optimization + VectorValues initialVectorValues; + initialVectorValues.insert(X(1), zero(1)); + initialVectorValues.insert(Y(1), zero(1)); + VectorValues expectedSolution = qpSolver.optimize(initialVectorValues).first; + + //Instantiate NLP + NLP nlp; + Values linPoint; + linPoint.insert(X(1), zero(1)); + linPoint.insert(Y(1), zero(1)); + nlp.cost.add(LinearContainerFactor(hf, linPoint)); // wrap it using linearcontainerfactor + nlp.linearInequalities.add(InequalityProblem1(X(1), Y(1), dualKey)); + + Values initialValues; + initialValues.insert(X(1), 1.0); + initialValues.insert(Y(1), -10.0); + + // Instantiate SQP + SQPSimple sqpSimple(nlp); + Values actualValues = sqpSimple.optimize(initialValues).first; + + DOUBLES_EQUAL(expectedSolution.at(X(1))[0], actualValues.at(X(1)), 1e-10); + DOUBLES_EQUAL(expectedSolution.at(Y(1))[0], actualValues.at(Y(1)), 1e-10); } //****************************************************************************** +const size_t X_AXIS = 0; +const size_t Y_AXIS = 1; +const size_t Z_AXIS = 2; + /** - * Inequality boundary constraint - * x <= bound + * Inequality boundary constraint on one axis (x, y or z) + * axis <= bound */ -class UpperBoundX : public NonlinearInequality1 { +class AxisUpperBound : public NonlinearInequality1 { typedef NonlinearInequality1 Base; + size_t axis_; double bound_; + public: - UpperBoundX(Key key, double bound, Key dualKey) : Base(key, dualKey, 1), bound_(bound) { + AxisUpperBound(Key key, size_t axis, double bound, Key dualKey) : Base(key, dualKey), axis_(axis), bound_(bound) { } double computeError(const Pose3& pose, boost::optional H = boost::none) const { if (H) - *H = (Matrix(1,6) << zeros(1,3), pose.rotation().matrix().row(0)).finished(); - return pose.x() - bound_; + *H = (Matrix(1,6) << zeros(1,3), pose.rotation().matrix().row(axis_)).finished(); + return pose.translation().vector()[axis_] - bound_; } }; -TEST(testSQPSimple, poseOnALine) { - const Key dualKey = 0; +/** + * Inequality boundary constraint on one axis (x, y or z) + * bound <= axis + */ +class AxisLowerBound : public NonlinearInequality1 { + typedef NonlinearInequality1 Base; + size_t axis_; + double bound_; +public: + AxisLowerBound(Key key, size_t axis, double bound, Key dualKey) : Base(key, dualKey), axis_(axis), bound_(bound) { + } + + double computeError(const Pose3& pose, boost::optional H = boost::none) const { + if (H) + *H = (Matrix(1,6) << zeros(1,3), -pose.rotation().matrix().row(axis_)).finished(); + return -pose.translation().vector()[axis_] + bound_; + } +}; + +TEST(testSQPSimple, poseWithABoundary) { + const Key dualKey = 0; //Instantiate NLP NLP nlp; - nlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(-1, 0, 0)), noiseModel::Unit::Create(6))); - UpperBoundX constraint(X(1), 0, dualKey); - nlp.nonlinearInequalities.add(constraint); + nlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(1, 0, 0)), noiseModel::Unit::Create(6))); + AxisUpperBound constraint(X(1), X_AXIS, 0, dualKey); + nlp.linearInequalities.add(constraint); Values initialValues; - initialValues.insert(X(1), Pose3(Rot3::ypr(0.3, 0.2, 0.3), Point3(-1,0,0))); + initialValues.insert(X(1), Pose3(Rot3::ypr(0.3, 0.2, 0.3), Point3(1, 0, 0))); Values expectedSolution; - expectedSolution.insert(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3())); + expectedSolution.insert(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(0, 0, 0))); // Instantiate SQP SQPSimple sqpSimple(nlp); Values actualSolution = sqpSimple.optimize(initialValues).first; CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); - actualSolution.print("actualSolution: "); +} + +TEST(testSQPSimple, poseWithinA2DBox) { + const Key dualKey = 0; + + //Instantiate NLP + NLP nlp; + nlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(10, 0.5, 0)), noiseModel::Unit::Create(6))); + nlp.linearInequalities.add(AxisLowerBound(X(1), X_AXIS, -1, dualKey)); + nlp.linearInequalities.add(AxisUpperBound(X(1), X_AXIS, 1, dualKey)); + nlp.linearInequalities.add(AxisLowerBound(X(1), Y_AXIS, -1, dualKey)); + nlp.linearInequalities.add(AxisUpperBound(X(1), Y_AXIS, 1, dualKey)); + + Values initialValues; + initialValues.insert(X(1), Pose3(Rot3::ypr(0.3, 0.2, 0.3), Point3(1, 0, 0))); + + Values expectedSolution; + expectedSolution.insert(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(1, 0.5, 0))); + + // Instantiate SQP + SQPSimple sqpSimple(nlp); + Values actualSolution = sqpSimple.optimize(initialValues).first; + + CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); + //TODO: remove printing, refactoring, } //****************************************************************************** From d435fa72dfb666fc0612d4ba03827e2073806cb1 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Mon, 22 Dec 2014 18:22:28 -0500 Subject: [PATCH 013/130] removed all printing --- gtsam_unstable/nonlinear/NonlinearConstraint.h | 3 --- gtsam_unstable/nonlinear/tests/testSQPSimple.cpp | 2 +- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/gtsam_unstable/nonlinear/NonlinearConstraint.h b/gtsam_unstable/nonlinear/NonlinearConstraint.h index af52ea178..d4c7f86b5 100644 --- a/gtsam_unstable/nonlinear/NonlinearConstraint.h +++ b/gtsam_unstable/nonlinear/NonlinearConstraint.h @@ -108,11 +108,8 @@ public: lG11sum += -lambda[i] * G11[i]; } - - std::cout << "lG11sum: " << lG11sum << std::endl; HessianFactor::shared_ptr hf(new HessianFactor(Base::key(), lG11sum, zero(X1Dim), 100.0)); - hf->print("HessianFactor: "); return hf; } diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index d83312ee1..b2be55289 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -215,7 +215,7 @@ public: * Single iteration of SQP */ SQPSimpleState iterate(const SQPSimpleState& state) const { - static const bool debug = true; + static const bool debug = false; // construct the qp subproblem QP qp; From c99a8481482034327531d9869bd70d82fae734ec Mon Sep 17 00:00:00 2001 From: krunalchande Date: Mon, 22 Dec 2014 18:35:22 -0500 Subject: [PATCH 014/130] Refactoring. --- .../LinearEqualityManifoldFactorGraph.h | 50 +++++++++ .../nonlinear/NonlinearEqualityFactorGraph.h | 27 +++++ .../NonlinearInequalityFactorGraph.h | 58 ++++++++++ .../nonlinear/tests/testSQPSimple.cpp | 104 +----------------- 4 files changed, 138 insertions(+), 101 deletions(-) create mode 100644 gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h create mode 100644 gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h create mode 100644 gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h diff --git a/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h b/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h new file mode 100644 index 000000000..50f9c3245 --- /dev/null +++ b/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h @@ -0,0 +1,50 @@ +/** + * @file LinearEqualityManifoldFactorGraph.h + * @author Duy-Nguyen Ta + * @author Krunal Chande + * @author Luca Carlone + * @date Dec 15, 2014 + */ + +#pragma once +#include + + +namespace gtsam { + +class LinearEqualityManifoldFactorGraph: public FactorGraph { +public: + /// default constructor + LinearEqualityManifoldFactorGraph() { + } + + /// linearize to a LinearEqualityFactorGraph + LinearEqualityFactorGraph::shared_ptr linearize( + const Values& linearizationPoint) const { + LinearEqualityFactorGraph::shared_ptr linearGraph( + new LinearEqualityFactorGraph()); + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + JacobianFactor::shared_ptr jacobian = boost::dynamic_pointer_cast( + factor->linearize(linearizationPoint)); + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + linearGraph->add(LinearEquality(*jacobian, constraint->dualKey())); + } + return linearGraph; + } + + /** + * Return true if the max absolute error all factors is less than a tolerance + */ + bool checkFeasibility(const Values& values, double tol) const { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( + factor); + Vector error = noiseModelFactor->unwhitenedError(values); + if (error.lpNorm() > tol) { + return false; + } + } + return true; + } +}; +} diff --git a/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h b/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h new file mode 100644 index 000000000..116b273fb --- /dev/null +++ b/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h @@ -0,0 +1,27 @@ +/** + * @file NonlinearEqualityFactorGraph.h + * @author Krunal Chande + * @date Dec 22, 2014 + */ + +#pragma once +#include + + +namespace gtsam { +class NonlinearEqualityFactorGraph: public LinearEqualityManifoldFactorGraph { +public: + /// default constructor + NonlinearEqualityFactorGraph() { + } + + GaussianFactorGraph::shared_ptr multipliedHessians(const Values& values, const VectorValues& duals) const { + GaussianFactorGraph::shared_ptr constrainedHessians(new GaussianFactorGraph()); + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this) { + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + constrainedHessians->push_back(constraint->multipliedHessian(values, duals)); + } + return constrainedHessians; + } +}; +} diff --git a/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h b/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h new file mode 100644 index 000000000..a76caeeb6 --- /dev/null +++ b/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h @@ -0,0 +1,58 @@ +/** + * @file NonlinearInequalityFactorGraph.h + * @author Krunal Chande + * @date Dec 22, 2014 + */ + +#pragma once +#include + +namespace gtsam { +class NonlinearInequalityFactorGraph : public FactorGraph { + +public: + /// default constructor + NonlinearInequalityFactorGraph() { + } + + /// linearize to a LinearEqualityFactorGraph + LinearInequalityFactorGraph::shared_ptr linearize( + const Values& linearizationPoint) const { + LinearInequalityFactorGraph::shared_ptr linearGraph( + new LinearInequalityFactorGraph()); + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + JacobianFactor::shared_ptr jacobian = boost::dynamic_pointer_cast( + factor->linearize(linearizationPoint)); + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + linearGraph->add(LinearInequality(*jacobian, constraint->dualKey())); + } + return linearGraph; + } + + /** + * Return true if the all errors are <= 0.0 + */ + bool checkFeasibilityAndComplimentary(const Values& values, const VectorValues& duals, double tol) const { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( + factor); + Vector error = noiseModelFactor->unwhitenedError(values); + + // Primal feasibility condition: all constraints need to be <= 0.0 + if (error[0] > tol) { + return false; + } + + // Complimentary condition: errors of active constraints need to be 0.0 + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast( + factor); + Key dualKey = constraint->dualKey(); + if (!duals.exists(dualKey)) continue; // if dualKey doesn't exist, it is an inactive constraint! + if (fabs(error[0]) > tol) // for active constraint, the error should be 0.0 + return false; + + } + return true; + } +}; +} diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index b2be55289..e8a14203b 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -24,112 +24,14 @@ #include #include #include +#include +#include +#include #include #include namespace gtsam { -class LinearEqualityManifoldFactorGraph: public FactorGraph { -public: - /// default constructor - LinearEqualityManifoldFactorGraph() { - } - - /// linearize to a LinearEqualityFactorGraph - LinearEqualityFactorGraph::shared_ptr linearize( - const Values& linearizationPoint) const { - LinearEqualityFactorGraph::shared_ptr linearGraph( - new LinearEqualityFactorGraph()); - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ - JacobianFactor::shared_ptr jacobian = boost::dynamic_pointer_cast( - factor->linearize(linearizationPoint)); - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); - linearGraph->add(LinearEquality(*jacobian, constraint->dualKey())); - } - return linearGraph; - } - - /** - * Return true if the max absolute error all factors is less than a tolerance - */ - bool checkFeasibility(const Values& values, double tol) const { - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ - NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( - factor); - Vector error = noiseModelFactor->unwhitenedError(values); - if (error.lpNorm() > tol) { - return false; - } - } - return true; - } -}; - -class NonlinearEqualityFactorGraph: public LinearEqualityManifoldFactorGraph { -public: - /// default constructor - NonlinearEqualityFactorGraph() { - } - - GaussianFactorGraph::shared_ptr multipliedHessians(const Values& values, const VectorValues& duals) const { - GaussianFactorGraph::shared_ptr constrainedHessians(new GaussianFactorGraph()); - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this) { - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); - constrainedHessians->push_back(constraint->multipliedHessian(values, duals)); - } - return constrainedHessians; - } -}; - -class NonlinearInequalityFactorGraph : public FactorGraph { - -public: - /// default constructor - NonlinearInequalityFactorGraph() { - } - - /// linearize to a LinearEqualityFactorGraph - LinearInequalityFactorGraph::shared_ptr linearize( - const Values& linearizationPoint) const { - LinearInequalityFactorGraph::shared_ptr linearGraph( - new LinearInequalityFactorGraph()); - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ - JacobianFactor::shared_ptr jacobian = boost::dynamic_pointer_cast( - factor->linearize(linearizationPoint)); - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); - linearGraph->add(LinearInequality(*jacobian, constraint->dualKey())); - } - return linearGraph; - } - - /** - * Return true if the all errors are <= 0.0 - */ - bool checkFeasibilityAndComplimentary(const Values& values, const VectorValues& duals, double tol) const { - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ - NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( - factor); - Vector error = noiseModelFactor->unwhitenedError(values); - - // Primal feasibility condition: all constraints need to be <= 0.0 - if (error[0] > tol) { - return false; - } - - // Complimentary condition: errors of active constraints need to be 0.0 - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast( - factor); - Key dualKey = constraint->dualKey(); - if (!duals.exists(dualKey)) continue; // if dualKey doesn't exist, it is an inactive constraint! - if (fabs(error[0]) > tol) // for active constraint, the error should be 0.0 - return false; - - } - return true; - } - -}; - struct NLP { NonlinearFactorGraph cost; NonlinearEqualityFactorGraph linearEqualities; From b6d85e83aed239a76225f84d9034fab93e536033 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Mon, 22 Dec 2014 18:56:19 -0500 Subject: [PATCH 015/130] Fixed includes. --- .../LinearEqualityManifoldFactorGraph.h | 2 + .../NonlinearInequalityFactorGraph.h | 1 + gtsam_unstable/nonlinear/SQPSimple.h | 172 ++++++++++++++++++ .../nonlinear/tests/testSQPSimple.cpp | 169 +---------------- 4 files changed, 179 insertions(+), 165 deletions(-) create mode 100644 gtsam_unstable/nonlinear/SQPSimple.h diff --git a/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h b/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h index 50f9c3245..03a01fd3c 100644 --- a/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h +++ b/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h @@ -8,6 +8,8 @@ #pragma once #include +#include +#include namespace gtsam { diff --git a/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h b/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h index a76caeeb6..301c66505 100644 --- a/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h +++ b/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h @@ -6,6 +6,7 @@ #pragma once #include +#include namespace gtsam { class NonlinearInequalityFactorGraph : public FactorGraph { diff --git a/gtsam_unstable/nonlinear/SQPSimple.h b/gtsam_unstable/nonlinear/SQPSimple.h new file mode 100644 index 000000000..df44f1cda --- /dev/null +++ b/gtsam_unstable/nonlinear/SQPSimple.h @@ -0,0 +1,172 @@ +/** + * @file SQPSimple.h + * @author Krunal Chande + * @date Dec 22, 2014 + */ + +#pragma once +#include +#include +#include +#include +#include +#include +#include + + + +namespace gtsam { + +struct NLP { + NonlinearFactorGraph cost; + NonlinearEqualityFactorGraph linearEqualities; + NonlinearEqualityFactorGraph nonlinearEqualities; + NonlinearInequalityFactorGraph linearInequalities; +}; + +struct SQPSimpleState { + Values values; + VectorValues duals; + bool converged; + size_t iterations; + + /// Default constructor + SQPSimpleState() : values(), duals(), converged(false), iterations(0) {} + + /// Constructor with an initialValues + SQPSimpleState(const Values& initialValues) : + values(initialValues), duals(VectorValues()), converged(false), iterations(0) { + } +}; + +/** + * Simple SQP optimizer to solve nonlinear constrained problems. + * This simple version won't care about nonconvexity, which needs + * more advanced techniques to solve, e.g., merit function, line search, second-order correction etc. + */ +class SQPSimple { + NLP nlp_; + static const double errorTol = 1e-5; +public: + SQPSimple(const NLP& nlp) : + nlp_(nlp) { + } + + /// Check if \nabla f(x) - \lambda * \nabla c(x) == 0 + bool isStationary(const VectorValues& delta) const { + return delta.vector().lpNorm() < errorTol; + } + + /// Check if c_E(x) == 0 + bool isPrimalFeasible(const SQPSimpleState& state) const { + return nlp_.linearEqualities.checkFeasibility(state.values, errorTol) + && nlp_.nonlinearEqualities.checkFeasibility(state.values, errorTol); + } + + /** + * Dual variables of inequality constraints need to be >=0 + * For active inequalities, the dual needs to be > 0 + * For inactive inequalities, they need to be == 0. However, we don't compute + * dual variables for inactive constraints in the qp subproblem, so we don't care. + */ + bool isDualFeasible(const VectorValues& duals) const { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearInequalities) { + NonlinearConstraint::shared_ptr inequality = boost::dynamic_pointer_cast(factor); + Key dualKey = inequality->dualKey(); + if (!duals.exists(dualKey)) continue; // should be inactive constraint! + double dual = duals.at(dualKey)[0]; // because we only support single-valued inequalities + if (dual < 0.0) + return false; + } + return true; + } + + /** + * Check complimentary slackness condition: + * For all inequality constraints, + * dual * constraintError(primals) == 0. + * If the constraint is active, we need to check constraintError(primals) == 0, and ignore the dual + * If it is inactive, the dual should be 0, regardless of the error. However, we don't compute + * dual variables for inactive constraints in the QP subproblem, so we don't care. + */ + bool isComplementary(const SQPSimpleState& state) const { + return nlp_.linearInequalities.checkFeasibilityAndComplimentary(state.values, state.duals, errorTol); + } + + /// Check convergence + bool checkConvergence(const SQPSimpleState& state, const VectorValues& delta) const { + return isStationary(delta) && isPrimalFeasible(state) && isDualFeasible(state.duals) && isComplementary(state); + } + + /** + * Single iteration of SQP + */ + SQPSimpleState iterate(const SQPSimpleState& state) const { + static const bool debug = false; + + // construct the qp subproblem + QP qp; + qp.cost = *nlp_.cost.linearize(state.values); + GaussianFactorGraph::shared_ptr multipliedHessians = nlp_.nonlinearEqualities.multipliedHessians(state.values, state.duals); + qp.cost.push_back(*multipliedHessians); + + qp.equalities.add(*nlp_.linearEqualities.linearize(state.values)); + qp.equalities.add(*nlp_.nonlinearEqualities.linearize(state.values)); + + qp.inequalities.add(*nlp_.linearInequalities.linearize(state.values)); + + if (debug) + qp.print("QP subproblem:"); + + // solve the QP subproblem + VectorValues delta, duals; + QPSolver qpSolver(qp); + boost::tie(delta, duals) = qpSolver.optimize(); + + if (debug) + delta.print("delta = "); + + if (debug) + duals.print("duals = "); + + // update new state + SQPSimpleState newState; + newState.values = state.values.retract(delta); + newState.duals = duals; + newState.converged = checkConvergence(newState, delta); + newState.iterations = state.iterations + 1; + + return newState; + } + + VectorValues initializeDuals() const { + VectorValues duals; + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearEqualities) { + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + duals.insert(constraint->dualKey(), zero(factor->dim())); + } + + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.nonlinearEqualities) { + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + duals.insert(constraint->dualKey(), zero(factor->dim())); + } + return duals; + } + + /** + * Main optimization function. + */ + std::pair optimize(const Values& initialValues) const { + SQPSimpleState state(initialValues); + state.duals = initializeDuals(); + + while (!state.converged && state.iterations < 100) { + state = iterate(state); + } + + return std::make_pair(state.values, state.duals); + } + + +}; +} diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index e8a14203b..19e9205cc 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -18,178 +18,17 @@ * @date Dec 15, 2014 */ +#include +#include #include -#include #include +#include #include -#include +#include #include -#include -#include -#include #include #include -namespace gtsam { - -struct NLP { - NonlinearFactorGraph cost; - NonlinearEqualityFactorGraph linearEqualities; - NonlinearEqualityFactorGraph nonlinearEqualities; - NonlinearInequalityFactorGraph linearInequalities; -}; - -struct SQPSimpleState { - Values values; - VectorValues duals; - bool converged; - size_t iterations; - - /// Default constructor - SQPSimpleState() : values(), duals(), converged(false), iterations(0) {} - - /// Constructor with an initialValues - SQPSimpleState(const Values& initialValues) : - values(initialValues), duals(VectorValues()), converged(false), iterations(0) { - } -}; - -/** - * Simple SQP optimizer to solve nonlinear constrained problems. - * This simple version won't care about nonconvexity, which needs - * more advanced techniques to solve, e.g., merit function, line search, second-order correction etc. - */ -class SQPSimple { - NLP nlp_; - static const double errorTol = 1e-5; -public: - SQPSimple(const NLP& nlp) : - nlp_(nlp) { - } - - /// Check if \nabla f(x) - \lambda * \nabla c(x) == 0 - bool isStationary(const VectorValues& delta) const { - return delta.vector().lpNorm() < errorTol; - } - - /// Check if c_E(x) == 0 - bool isPrimalFeasible(const SQPSimpleState& state) const { - return nlp_.linearEqualities.checkFeasibility(state.values, errorTol) - && nlp_.nonlinearEqualities.checkFeasibility(state.values, errorTol); - } - - /** - * Dual variables of inequality constraints need to be >=0 - * For active inequalities, the dual needs to be > 0 - * For inactive inequalities, they need to be == 0. However, we don't compute - * dual variables for inactive constraints in the qp subproblem, so we don't care. - */ - bool isDualFeasible(const VectorValues& duals) const { - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearInequalities) { - NonlinearConstraint::shared_ptr inequality = boost::dynamic_pointer_cast(factor); - Key dualKey = inequality->dualKey(); - if (!duals.exists(dualKey)) continue; // should be inactive constraint! - double dual = duals.at(dualKey)[0]; // because we only support single-valued inequalities - if (dual < 0.0) - return false; - } - return true; - } - - /** - * Check complimentary slackness condition: - * For all inequality constraints, - * dual * constraintError(primals) == 0. - * If the constraint is active, we need to check constraintError(primals) == 0, and ignore the dual - * If it is inactive, the dual should be 0, regardless of the error. However, we don't compute - * dual variables for inactive constraints in the QP subproblem, so we don't care. - */ - bool isComplementary(const SQPSimpleState& state) const { - return nlp_.linearInequalities.checkFeasibilityAndComplimentary(state.values, state.duals, errorTol); - } - - /// Check convergence - bool checkConvergence(const SQPSimpleState& state, const VectorValues& delta) const { - return isStationary(delta) && isPrimalFeasible(state) && isDualFeasible(state.duals) && isComplementary(state); - } - - /** - * Single iteration of SQP - */ - SQPSimpleState iterate(const SQPSimpleState& state) const { - static const bool debug = false; - - // construct the qp subproblem - QP qp; - qp.cost = *nlp_.cost.linearize(state.values); - GaussianFactorGraph::shared_ptr multipliedHessians = nlp_.nonlinearEqualities.multipliedHessians(state.values, state.duals); - qp.cost.push_back(*multipliedHessians); - - qp.equalities.add(*nlp_.linearEqualities.linearize(state.values)); - qp.equalities.add(*nlp_.nonlinearEqualities.linearize(state.values)); - - qp.inequalities.add(*nlp_.linearInequalities.linearize(state.values)); - - if (debug) - qp.print("QP subproblem:"); - - // solve the QP subproblem - VectorValues delta, duals; - QPSolver qpSolver(qp); - boost::tie(delta, duals) = qpSolver.optimize(); - - if (debug) - delta.print("delta = "); - - if (debug) - duals.print("duals = "); - - // update new state - SQPSimpleState newState; - newState.values = state.values.retract(delta); - newState.duals = duals; - newState.converged = checkConvergence(newState, delta); - newState.iterations = state.iterations + 1; - - return newState; - } - - VectorValues initializeDuals() const { - VectorValues duals; - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearEqualities) { - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); - duals.insert(constraint->dualKey(), zero(factor->dim())); - } - - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.nonlinearEqualities) { - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); - duals.insert(constraint->dualKey(), zero(factor->dim())); - } - return duals; - } - - /** - * Main optimization function. - */ - std::pair optimize(const Values& initialValues) const { - SQPSimpleState state(initialValues); - state.duals = initializeDuals(); - - while (!state.converged && state.iterations < 100) { - state = iterate(state); - } - - return std::make_pair(state.values, state.duals); - } - - -}; -} - -#include -#include -#include - using namespace std; using namespace gtsam::symbol_shorthand; using namespace gtsam; From 2523fa2fe5d209340d0070cccc186331b31a73ff Mon Sep 17 00:00:00 2001 From: krunalchande Date: Mon, 22 Dec 2014 19:14:49 -0500 Subject: [PATCH 016/130] removed comment --- gtsam_unstable/nonlinear/tests/testSQPSimple.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index 19e9205cc..94e7b08e4 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -360,7 +360,6 @@ TEST(testSQPSimple, poseWithinA2DBox) { Values actualSolution = sqpSimple.optimize(initialValues).first; CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); - //TODO: remove printing, refactoring, } //****************************************************************************** From de7149af63fc9521701dda30875c570ff658db01 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Mon, 22 Dec 2014 19:49:32 -0500 Subject: [PATCH 017/130] Added licensing information. --- gtsam_unstable/linear/QPSolver.h | 10 ++ .../LinearEqualityManifoldFactorGraph.h | 10 ++ .../nonlinear/NonlinearConstraint.h | 10 ++ .../nonlinear/NonlinearEqualityFactorGraph.h | 19 ++- .../nonlinear/NonlinearInequality.h | 10 ++ .../NonlinearInequalityFactorGraph.h | 19 ++- gtsam_unstable/nonlinear/SQPSimple.cpp | 124 ++++++++++++++++++ gtsam_unstable/nonlinear/SQPSimple.h | 115 ++++------------ .../nonlinear/tests/testSQPSimple.cpp | 2 +- 9 files changed, 221 insertions(+), 98 deletions(-) create mode 100644 gtsam_unstable/nonlinear/SQPSimple.cpp diff --git a/gtsam_unstable/linear/QPSolver.h b/gtsam_unstable/linear/QPSolver.h index f7a575f8c..53cb6ca6d 100644 --- a/gtsam_unstable/linear/QPSolver.h +++ b/gtsam_unstable/linear/QPSolver.h @@ -1,3 +1,13 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ /* * QPSolver.h * @brief: A quadratic programming solver implements the active set method diff --git a/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h b/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h index 03a01fd3c..a20f2e5a3 100644 --- a/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h +++ b/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h @@ -1,3 +1,13 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ /** * @file LinearEqualityManifoldFactorGraph.h * @author Duy-Nguyen Ta diff --git a/gtsam_unstable/nonlinear/NonlinearConstraint.h b/gtsam_unstable/nonlinear/NonlinearConstraint.h index d4c7f86b5..81208872f 100644 --- a/gtsam_unstable/nonlinear/NonlinearConstraint.h +++ b/gtsam_unstable/nonlinear/NonlinearConstraint.h @@ -1,3 +1,13 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ /** * @file NonlinearConstraint.h * @brief diff --git a/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h b/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h index 116b273fb..1d82c6fbf 100644 --- a/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h +++ b/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h @@ -1,7 +1,20 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ + /** - * @file NonlinearEqualityFactorGraph.h - * @author Krunal Chande - * @date Dec 22, 2014 + * @file NonlinearEqualityFactorGraph.h + * @author Duy-Nguyen Ta + * @author Krunal Chande + * @author Luca Carlone + * @date Dec 15, 2014 */ #pragma once diff --git a/gtsam_unstable/nonlinear/NonlinearInequality.h b/gtsam_unstable/nonlinear/NonlinearInequality.h index 8d41b3ec0..439d5e070 100644 --- a/gtsam_unstable/nonlinear/NonlinearInequality.h +++ b/gtsam_unstable/nonlinear/NonlinearInequality.h @@ -1,3 +1,13 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ /** * @file NonlinearConstraint.h * @brief diff --git a/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h b/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h index 301c66505..0eccd1fb0 100644 --- a/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h +++ b/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h @@ -1,7 +1,20 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ + /** - * @file NonlinearInequalityFactorGraph.h - * @author Krunal Chande - * @date Dec 22, 2014 + * @file NonlinearInequalityFactorGraph.h + * @author Duy-Nguyen Ta + * @author Krunal Chande + * @author Luca Carlone + * @date Dec 15, 2014 */ #pragma once diff --git a/gtsam_unstable/nonlinear/SQPSimple.cpp b/gtsam_unstable/nonlinear/SQPSimple.cpp new file mode 100644 index 000000000..35cfb0a1b --- /dev/null +++ b/gtsam_unstable/nonlinear/SQPSimple.cpp @@ -0,0 +1,124 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ + +/** + * @file SQPSimple.cpp + * @author Duy-Nguyen Ta + * @author Krunal Chande + * @author Luca Carlone + * @date Dec 15, 2014 + */ + +#include +namespace gtsam { + + +/* ************************************************************************* */ +bool SQPSimple::isStationary(const VectorValues& delta) const { + return delta.vector().lpNorm() < errorTol; +} + +/* ************************************************************************* */ +bool SQPSimple::isPrimalFeasible(const SQPSimpleState& state) const { + return nlp_.linearEqualities.checkFeasibility(state.values, errorTol) + && nlp_.nonlinearEqualities.checkFeasibility(state.values, errorTol); +} + +/* ************************************************************************* */ +bool SQPSimple::isDualFeasible(const VectorValues& duals) const { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearInequalities) { + NonlinearConstraint::shared_ptr inequality = boost::dynamic_pointer_cast(factor); + Key dualKey = inequality->dualKey(); + if (!duals.exists(dualKey)) continue; // should be inactive constraint! + double dual = duals.at(dualKey)[0]; // because we only support single-valued inequalities + if (dual < 0.0) + return false; + } + return true; +} + +/* ************************************************************************* */ +bool SQPSimple::isComplementary(const SQPSimpleState& state) const { + return nlp_.linearInequalities.checkFeasibilityAndComplimentary(state.values, state.duals, errorTol); +} + +/* ************************************************************************* */ +bool SQPSimple::checkConvergence(const SQPSimpleState& state, const VectorValues& delta) const { + return isStationary(delta) && isPrimalFeasible(state) && isDualFeasible(state.duals) && isComplementary(state); +} + +/* ************************************************************************* */ +VectorValues SQPSimple::initializeDuals() const { + VectorValues duals; + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearEqualities) { + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + duals.insert(constraint->dualKey(), zero(factor->dim())); + } + + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.nonlinearEqualities) { + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + duals.insert(constraint->dualKey(), zero(factor->dim())); + } + return duals; +} + +/* ************************************************************************* */ +std::pair SQPSimple::optimize(const Values& initialValues) const { + SQPSimpleState state(initialValues); + state.duals = initializeDuals(); + while (!state.converged && state.iterations < 100) { + state = iterate(state); + } + return std::make_pair(state.values, state.duals); +} + +/* ************************************************************************* */ +SQPSimpleState SQPSimple::iterate(const SQPSimpleState& state) const { + static const bool debug = false; + + // construct the qp subproblem + QP qp; + qp.cost = *nlp_.cost.linearize(state.values); + GaussianFactorGraph::shared_ptr multipliedHessians = nlp_.nonlinearEqualities.multipliedHessians(state.values, state.duals); + qp.cost.push_back(*multipliedHessians); + + qp.equalities.add(*nlp_.linearEqualities.linearize(state.values)); + qp.equalities.add(*nlp_.nonlinearEqualities.linearize(state.values)); + + qp.inequalities.add(*nlp_.linearInequalities.linearize(state.values)); + + if (debug) + qp.print("QP subproblem:"); + + // solve the QP subproblem + VectorValues delta, duals; + QPSolver qpSolver(qp); + boost::tie(delta, duals) = qpSolver.optimize(); + + if (debug) + delta.print("delta = "); + + if (debug) + duals.print("duals = "); + + // update new state + SQPSimpleState newState; + newState.values = state.values.retract(delta); + newState.duals = duals; + newState.converged = checkConvergence(newState, delta); + newState.iterations = state.iterations + 1; + + return newState; +} +} + + + diff --git a/gtsam_unstable/nonlinear/SQPSimple.h b/gtsam_unstable/nonlinear/SQPSimple.h index df44f1cda..ceda522ea 100644 --- a/gtsam_unstable/nonlinear/SQPSimple.h +++ b/gtsam_unstable/nonlinear/SQPSimple.h @@ -1,7 +1,20 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ + /** - * @file SQPSimple.h - * @author Krunal Chande - * @date Dec 22, 2014 + * @file SQPSimple.h + * @author Duy-Nguyen Ta + * @author Krunal Chande + * @author Luca Carlone + * @date Dec 15, 2014 */ #pragma once @@ -13,8 +26,6 @@ #include #include - - namespace gtsam { struct NLP { @@ -53,15 +64,10 @@ public: } /// Check if \nabla f(x) - \lambda * \nabla c(x) == 0 - bool isStationary(const VectorValues& delta) const { - return delta.vector().lpNorm() < errorTol; - } + bool isStationary(const VectorValues& delta) const; /// Check if c_E(x) == 0 - bool isPrimalFeasible(const SQPSimpleState& state) const { - return nlp_.linearEqualities.checkFeasibility(state.values, errorTol) - && nlp_.nonlinearEqualities.checkFeasibility(state.values, errorTol); - } + bool isPrimalFeasible(const SQPSimpleState& state) const; /** * Dual variables of inequality constraints need to be >=0 @@ -69,17 +75,7 @@ public: * For inactive inequalities, they need to be == 0. However, we don't compute * dual variables for inactive constraints in the qp subproblem, so we don't care. */ - bool isDualFeasible(const VectorValues& duals) const { - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearInequalities) { - NonlinearConstraint::shared_ptr inequality = boost::dynamic_pointer_cast(factor); - Key dualKey = inequality->dualKey(); - if (!duals.exists(dualKey)) continue; // should be inactive constraint! - double dual = duals.at(dualKey)[0]; // because we only support single-valued inequalities - if (dual < 0.0) - return false; - } - return true; - } + bool isDualFeasible(const VectorValues& duals) const; /** * Check complimentary slackness condition: @@ -89,84 +85,21 @@ public: * If it is inactive, the dual should be 0, regardless of the error. However, we don't compute * dual variables for inactive constraints in the QP subproblem, so we don't care. */ - bool isComplementary(const SQPSimpleState& state) const { - return nlp_.linearInequalities.checkFeasibilityAndComplimentary(state.values, state.duals, errorTol); - } + bool isComplementary(const SQPSimpleState& state) const; /// Check convergence - bool checkConvergence(const SQPSimpleState& state, const VectorValues& delta) const { - return isStationary(delta) && isPrimalFeasible(state) && isDualFeasible(state.duals) && isComplementary(state); - } + bool checkConvergence(const SQPSimpleState& state, const VectorValues& delta) const; /** * Single iteration of SQP */ - SQPSimpleState iterate(const SQPSimpleState& state) const { - static const bool debug = false; - - // construct the qp subproblem - QP qp; - qp.cost = *nlp_.cost.linearize(state.values); - GaussianFactorGraph::shared_ptr multipliedHessians = nlp_.nonlinearEqualities.multipliedHessians(state.values, state.duals); - qp.cost.push_back(*multipliedHessians); - - qp.equalities.add(*nlp_.linearEqualities.linearize(state.values)); - qp.equalities.add(*nlp_.nonlinearEqualities.linearize(state.values)); - - qp.inequalities.add(*nlp_.linearInequalities.linearize(state.values)); - - if (debug) - qp.print("QP subproblem:"); - - // solve the QP subproblem - VectorValues delta, duals; - QPSolver qpSolver(qp); - boost::tie(delta, duals) = qpSolver.optimize(); - - if (debug) - delta.print("delta = "); - - if (debug) - duals.print("duals = "); - - // update new state - SQPSimpleState newState; - newState.values = state.values.retract(delta); - newState.duals = duals; - newState.converged = checkConvergence(newState, delta); - newState.iterations = state.iterations + 1; - - return newState; - } - - VectorValues initializeDuals() const { - VectorValues duals; - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearEqualities) { - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); - duals.insert(constraint->dualKey(), zero(factor->dim())); - } - - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.nonlinearEqualities) { - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); - duals.insert(constraint->dualKey(), zero(factor->dim())); - } - return duals; - } + SQPSimpleState iterate(const SQPSimpleState& state) const; + VectorValues initializeDuals() const; /** * Main optimization function. */ - std::pair optimize(const Values& initialValues) const { - SQPSimpleState state(initialValues); - state.duals = initializeDuals(); - - while (!state.converged && state.iterations < 100) { - state = iterate(state); - } - - return std::make_pair(state.values, state.duals); - } - + std::pair optimize(const Values& initialValues) const; }; } diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index 94e7b08e4..a14e52d71 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -12,8 +12,8 @@ /** * @file testQPSimple.cpp * @brief Unit tests for testQPSimple - * @author Krunal Chande * @author Duy-Nguyen Ta + * @author Krunal Chande * @author Luca Carlone * @date Dec 15, 2014 */ From 1dd23ced0211583698a815a9bd81caebcabad2ab Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 14:43:30 -0500 Subject: [PATCH 018/130] Simple unit test for NonlinearInequalityFactorGraph. --- .../testNonlinearInequalityFactorGraph.cpp | 43 +++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 gtsam_unstable/nonlinear/tests/testNonlinearInequalityFactorGraph.cpp diff --git a/gtsam_unstable/nonlinear/tests/testNonlinearInequalityFactorGraph.cpp b/gtsam_unstable/nonlinear/tests/testNonlinearInequalityFactorGraph.cpp new file mode 100644 index 000000000..503528c85 --- /dev/null +++ b/gtsam_unstable/nonlinear/tests/testNonlinearInequalityFactorGraph.cpp @@ -0,0 +1,43 @@ +/* ---------------------------------------------------------------------------- + + * GTSAM Copyright 2010, Georgia Tech Research Corporation, + * Atlanta, Georgia 30332-0415 + * All Rights Reserved + * Authors: Frank Dellaert, et al. (see THANKS for the full author list) + + * See LICENSE for the license information + + * -------------------------------------------------------------------------- */ + +/** + * @file testQPSimple.cpp + * @brief Unit tests for testQPSimple + * @author Duy-Nguyen Ta + * @author Krunal Chande + * @author Luca Carlone + * @date Dec 15, 2014 + */ + +#include +#include +#include +#include + +using namespace std; +using namespace gtsam::symbol_shorthand; +using namespace gtsam; +const double tol = 1e-10; + +//****************************************************************************** +TEST(NonlinearInequalityFactorGraph, constructor) { + NonlinearInequalityFactorGraph nonlinearInequalities; + CHECK(nonlinearInequalities.empty()); +} + + +//****************************************************************************** +int main() { + TestResult tr; + return TestRegistry::runAllTests(tr); +} +//****************************************************************************** From 79b69d44890924df564bb479f5b882e9a38212f3 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 14:44:11 -0500 Subject: [PATCH 019/130] [unfinished] unit test with multiple betweenFactors with box constraint. --- .../nonlinear/tests/testSQPSimple.cpp | 71 +++++++++++++++++-- 1 file changed, 65 insertions(+), 6 deletions(-) diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp index a14e52d71..dc406acd8 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp @@ -23,6 +23,8 @@ #include #include #include +#include + #include #include #include @@ -124,7 +126,8 @@ public: }; -TEST_UNSAFE(testSQPSimple, quadraticCostNonlinearConstraint) { +// Nonlinear constraint not supported currently +TEST_DISABLED(testSQPSimple, quadraticCostNonlinearConstraint) { const Key dualKey = 0; //Instantiate NLP @@ -344,13 +347,13 @@ TEST(testSQPSimple, poseWithinA2DBox) { //Instantiate NLP NLP nlp; nlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(10, 0.5, 0)), noiseModel::Unit::Create(6))); - nlp.linearInequalities.add(AxisLowerBound(X(1), X_AXIS, -1, dualKey)); - nlp.linearInequalities.add(AxisUpperBound(X(1), X_AXIS, 1, dualKey)); - nlp.linearInequalities.add(AxisLowerBound(X(1), Y_AXIS, -1, dualKey)); - nlp.linearInequalities.add(AxisUpperBound(X(1), Y_AXIS, 1, dualKey)); + nlp.linearInequalities.add(AxisLowerBound(X(1), X_AXIS, -1, dualKey)); // -1 <= x + nlp.linearInequalities.add(AxisUpperBound(X(1), X_AXIS, 1, dualKey+1)); // x <= 1 + nlp.linearInequalities.add(AxisLowerBound(X(1), Y_AXIS, -1, dualKey+2)); // -1 <= y + nlp.linearInequalities.add(AxisUpperBound(X(1), Y_AXIS, 1, dualKey+3));// y <= 1 Values initialValues; - initialValues.insert(X(1), Pose3(Rot3::ypr(0.3, 0.2, 0.3), Point3(1, 0, 0))); + initialValues.insert(X(1), Pose3(Rot3::ypr(1, -1, 2), Point3(3, -5, 0))); Values expectedSolution; expectedSolution.insert(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(1, 0.5, 0))); @@ -362,6 +365,62 @@ TEST(testSQPSimple, poseWithinA2DBox) { CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); } +TEST_DISABLED(testSQPSimple, posesInA2DBox) { + const double xLowerBound = -3.0, + xUpperBound = 5.0, + yLowerBound = -1.0, + yUpperBound = 2.0, + zLowerBound = 0.0, + zUpperBound = 2.0; + + //Instantiate NLP + NLP nlp; + + // prior on the first pose + SharedDiagonal priorNoise = noiseModel::Diagonal::Sigmas( + (Vector(6) << 0.001, 0.001, 0.001, 0.001, 0.001, 0.001).finished()); + nlp.cost.add(PriorFactor(X(1), Pose3(), priorNoise)); + + // odometry between factor for subsequent poses + SharedDiagonal odoNoise = noiseModel::Diagonal::Sigmas( + (Vector(6) << 0.01, 0.01, 0.01, 0.1, 0.1, 0.1).finished()); + Pose3 odo12(Rot3::ypr(M_PI/2.0, 0, 0), Point3(10, 0, 0)); + nlp.cost.add(BetweenFactor(X(1), X(2), odo12, odoNoise)); + + Pose3 odo23(Rot3::ypr(M_PI/2.0, 0, 0), Point3(2, 0, 2)); + nlp.cost.add(BetweenFactor(X(2), X(3), odo23, odoNoise)); + + // Box constraints + Key dualKey = 0; + for (size_t i=1; i<=3; ++i) { + nlp.linearInequalities.add(AxisLowerBound(X(i), X_AXIS, xLowerBound, dualKey++)); + nlp.linearInequalities.add(AxisUpperBound(X(i), X_AXIS, xUpperBound, dualKey++)); + nlp.linearInequalities.add(AxisLowerBound(X(i), Y_AXIS, yLowerBound, dualKey++)); + nlp.linearInequalities.add(AxisUpperBound(X(i), Y_AXIS, yUpperBound, dualKey++)); + nlp.linearInequalities.add(AxisLowerBound(X(i), Z_AXIS, zLowerBound, dualKey++)); + nlp.linearInequalities.add(AxisUpperBound(X(i), Z_AXIS, zUpperBound, dualKey++)); + } + + Values initialValues; + initialValues.insert(X(1), Pose3(Rot3(), Point3(0, 0, 0))); + initialValues.insert(X(2), Pose3(Rot3(), Point3(0, 0, 0))); + initialValues.insert(X(3), Pose3(Rot3(), Point3(0, 0, 0))); + + + Values expectedSolution; + expectedSolution.insert(X(1), Pose3()); + expectedSolution.insert(X(2), Pose3()); + expectedSolution.insert(X(3), Pose3()); + + // Instantiate SQP + SQPSimple sqpSimple(nlp); + Values actualSolution = sqpSimple.optimize(initialValues).first; + + actualSolution.print("actual solution: "); + + CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); +} + //****************************************************************************** int main() { TestResult tr; From e0e5e72020c3b5eceaab1aa9e81758c440e6f918 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 14:44:43 -0500 Subject: [PATCH 020/130] Fixed includes --- gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h | 7 ++++--- gtsam_unstable/nonlinear/SQPSimple.cpp | 2 ++ gtsam_unstable/nonlinear/SQPSimple.h | 4 ---- 3 files changed, 6 insertions(+), 7 deletions(-) diff --git a/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h b/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h index 0eccd1fb0..20196c91b 100644 --- a/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h +++ b/gtsam_unstable/nonlinear/NonlinearInequalityFactorGraph.h @@ -18,18 +18,19 @@ */ #pragma once -#include +#include #include +#include namespace gtsam { class NonlinearInequalityFactorGraph : public FactorGraph { public: - /// default constructor + /// Default constructor NonlinearInequalityFactorGraph() { } - /// linearize to a LinearEqualityFactorGraph + /// Linearize to a LinearEqualityFactorGraph LinearInequalityFactorGraph::shared_ptr linearize( const Values& linearizationPoint) const { LinearInequalityFactorGraph::shared_ptr linearGraph( diff --git a/gtsam_unstable/nonlinear/SQPSimple.cpp b/gtsam_unstable/nonlinear/SQPSimple.cpp index 35cfb0a1b..b0e3e8e5a 100644 --- a/gtsam_unstable/nonlinear/SQPSimple.cpp +++ b/gtsam_unstable/nonlinear/SQPSimple.cpp @@ -18,6 +18,8 @@ */ #include +#include + namespace gtsam { diff --git a/gtsam_unstable/nonlinear/SQPSimple.h b/gtsam_unstable/nonlinear/SQPSimple.h index ceda522ea..42b5c451a 100644 --- a/gtsam_unstable/nonlinear/SQPSimple.h +++ b/gtsam_unstable/nonlinear/SQPSimple.h @@ -18,13 +18,9 @@ */ #pragma once -#include #include #include #include -#include -#include -#include namespace gtsam { From 0fdd49ca4e2bacf1196e49aabc6b3eca1c39e7f0 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 14:45:23 -0500 Subject: [PATCH 021/130] Removed LinearEqualityManifoldFactorGraph. --- .../LinearEqualityManifoldFactorGraph.h | 62 ------------------- .../nonlinear/NonlinearEqualityFactorGraph.h | 43 +++++++++++-- 2 files changed, 39 insertions(+), 66 deletions(-) delete mode 100644 gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h diff --git a/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h b/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h deleted file mode 100644 index a20f2e5a3..000000000 --- a/gtsam_unstable/nonlinear/LinearEqualityManifoldFactorGraph.h +++ /dev/null @@ -1,62 +0,0 @@ -/* ---------------------------------------------------------------------------- - - * GTSAM Copyright 2010, Georgia Tech Research Corporation, - * Atlanta, Georgia 30332-0415 - * All Rights Reserved - * Authors: Frank Dellaert, et al. (see THANKS for the full author list) - - * See LICENSE for the license information - - * -------------------------------------------------------------------------- */ -/** - * @file LinearEqualityManifoldFactorGraph.h - * @author Duy-Nguyen Ta - * @author Krunal Chande - * @author Luca Carlone - * @date Dec 15, 2014 - */ - -#pragma once -#include -#include -#include - - -namespace gtsam { - -class LinearEqualityManifoldFactorGraph: public FactorGraph { -public: - /// default constructor - LinearEqualityManifoldFactorGraph() { - } - - /// linearize to a LinearEqualityFactorGraph - LinearEqualityFactorGraph::shared_ptr linearize( - const Values& linearizationPoint) const { - LinearEqualityFactorGraph::shared_ptr linearGraph( - new LinearEqualityFactorGraph()); - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ - JacobianFactor::shared_ptr jacobian = boost::dynamic_pointer_cast( - factor->linearize(linearizationPoint)); - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); - linearGraph->add(LinearEquality(*jacobian, constraint->dualKey())); - } - return linearGraph; - } - - /** - * Return true if the max absolute error all factors is less than a tolerance - */ - bool checkFeasibility(const Values& values, double tol) const { - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ - NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( - factor); - Vector error = noiseModelFactor->unwhitenedError(values); - if (error.lpNorm() > tol) { - return false; - } - } - return true; - } -}; -} diff --git a/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h b/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h index 1d82c6fbf..fbc52f376 100644 --- a/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h +++ b/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h @@ -18,16 +18,49 @@ */ #pragma once -#include - +#include +#include namespace gtsam { -class NonlinearEqualityFactorGraph: public LinearEqualityManifoldFactorGraph { + +class NonlinearEqualityFactorGraph: public FactorGraph { public: - /// default constructor + /// Default constructor NonlinearEqualityFactorGraph() { } + /// Linearize to a LinearEqualityFactorGraph + LinearEqualityFactorGraph::shared_ptr linearize( + const Values& linearizationPoint) const { + LinearEqualityFactorGraph::shared_ptr linearGraph( + new LinearEqualityFactorGraph()); + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + JacobianFactor::shared_ptr jacobian = boost::dynamic_pointer_cast( + factor->linearize(linearizationPoint)); + NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); + linearGraph->add(LinearEquality(*jacobian, constraint->dualKey())); + } + return linearGraph; + } + + /** + * Return true if the max absolute error all factors is less than a tolerance + */ + bool checkFeasibility(const Values& values, double tol) const { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this){ + NoiseModelFactor::shared_ptr noiseModelFactor = boost::dynamic_pointer_cast( + factor); + Vector error = noiseModelFactor->unwhitenedError(values); + if (error.lpNorm() > tol) { + return false; + } + } + return true; + } + + /** + * Additional cost for -lambda*ConstraintHessian for SQP + */ GaussianFactorGraph::shared_ptr multipliedHessians(const Values& values, const VectorValues& duals) const { GaussianFactorGraph::shared_ptr constrainedHessians(new GaussianFactorGraph()); BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this) { @@ -36,5 +69,7 @@ public: } return constrainedHessians; } + }; + } From 276959e39a65c24d8a650d42762e22ecbe3483ce Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 15:12:53 -0500 Subject: [PATCH 022/130] remove support for nonlinear constraints. Refactor SQPSimple to LCNLPSolver. --- gtsam/linear/GaussianFactorGraph.cpp | 23 --- gtsam/linear/GaussianFactorGraph.h | 6 - gtsam/linear/HessianFactor.cpp | 106 +--------- .../{SQPSimple.cpp => LCNLPSolver.cpp} | 48 ++--- .../nonlinear/{SQPSimple.h => LCNLPSolver.h} | 39 ++-- .../nonlinear/NonlinearEqualityFactorGraph.h | 22 +-- ...{testSQPSimple.cpp => testLCNLPSolver.cpp} | 184 ++++++------------ 7 files changed, 114 insertions(+), 314 deletions(-) rename gtsam_unstable/nonlinear/{SQPSimple.cpp => LCNLPSolver.cpp} (60%) rename gtsam_unstable/nonlinear/{SQPSimple.h => LCNLPSolver.h} (73%) rename gtsam_unstable/nonlinear/tests/{testSQPSimple.cpp => testLCNLPSolver.cpp} (63%) diff --git a/gtsam/linear/GaussianFactorGraph.cpp b/gtsam/linear/GaussianFactorGraph.cpp index 1d6066c1c..6953d2969 100644 --- a/gtsam/linear/GaussianFactorGraph.cpp +++ b/gtsam/linear/GaussianFactorGraph.cpp @@ -384,29 +384,6 @@ namespace gtsam { return false; } - /* ************************************************************************* */ - boost::tuple GaussianFactorGraph::splitConstraints() const { - typedef HessianFactor H; - typedef JacobianFactor J; - - GaussianFactorGraph hessians, jacobians, constraints; - BOOST_FOREACH(const GaussianFactor::shared_ptr& factor, *this) { - H::shared_ptr hessian(boost::dynamic_pointer_cast(factor)); - if (hessian) - hessians.push_back(factor); - else { - J::shared_ptr jacobian(boost::dynamic_pointer_cast(factor)); - if (jacobian && jacobian->get_model() && jacobian->get_model()->isConstrained()) { - constraints.push_back(jacobian); - } - else { - jacobians.push_back(factor); - } - } - } - return boost::make_tuple(hessians, jacobians, constraints); - } - /* ************************************************************************* */ // x += alpha*A'*e void GaussianFactorGraph::transposeMultiplyAdd(double alpha, const Errors& e, diff --git a/gtsam/linear/GaussianFactorGraph.h b/gtsam/linear/GaussianFactorGraph.h index 067a42d03..811c0f8b0 100644 --- a/gtsam/linear/GaussianFactorGraph.h +++ b/gtsam/linear/GaussianFactorGraph.h @@ -316,12 +316,6 @@ namespace gtsam { /** In-place version e <- A*x that takes an iterator. */ void multiplyInPlace(const VectorValues& x, const Errors::iterator& e) const; - /** - * Split constraints and unconstrained factors into two different graphs - * @return a pair of graphs - */ - boost::tuple splitConstraints() const; - /// @} private: diff --git a/gtsam/linear/HessianFactor.cpp b/gtsam/linear/HessianFactor.cpp index 755f9c313..f2bebcab9 100644 --- a/gtsam/linear/HessianFactor.cpp +++ b/gtsam/linear/HessianFactor.cpp @@ -616,110 +616,10 @@ EliminatePreferCholesky(const GaussianFactorGraph& factors, const Ordering& keys // all factors to JacobianFactors. Otherwise, we can convert all factors // to HessianFactors. This is because QR can handle constrained noise // models but Cholesky cannot. - - /* Currently, when eliminating a constrained variable, EliminatePreferCholesky - * converts every other factors to JacobianFactor before doing the special QR - * factorization for constrained variables. Unfortunately, after a constrained - * nonlinear graph is linearized, new hessian factors from constraints, multiplied - * with the dual variable (-lambda*\hessian{c} terms in the Lagrangian objective - * function), might become negative definite, thus cannot be converted to JacobianFactors. - * - * Following EliminateCholesky, this version of EliminatePreferCholesky for - * constrained var gathers all unconstrained factors into a big joint HessianFactor - * before converting it into a JacobianFactor to be eliminiated by QR together with - * the other constrained factors. - * - * Of course, this might not solve the non-positive-definite problem entirely, - * because (1) the original hessian factors might be non-positive definite - * and (2) large strange value of lambdas might cause the joint factor non-positive - * definite [is this true?]. But at least, this will help in typical cases. - */ - GaussianFactorGraph hessians, jacobians, constraints; -// factors.print("factors: "); - boost::tie(hessians, jacobians, constraints) = factors.splitConstraints(); -// keys.print("Frontal variables to eliminate: "); -// hessians.print("Hessians: "); -// jacobians.print("Jacobians: "); -// constraints.print("Constraints: "); - - bool hasHessians = hessians.size() > 0; - - // Add all jacobians to gather as much info as we can - hessians.push_back(jacobians); - - if (constraints.size()>0) { -// // Build joint factor -// HessianFactor::shared_ptr jointFactor; -// try { -// jointFactor = boost::make_shared(hessians, Scatter(factors, keys)); -// } catch(std::invalid_argument&) { -// throw InvalidDenseElimination( -// "EliminateCholesky was called with a request to eliminate variables that are not\n" -// "involved in the provided factors."); -// } -// constraints.push_back(jointFactor); -// return EliminateQR(constraints, keys); - - // If there are hessian factors, turn them into conditional - // by doing partial elimination, then use those jacobians when eliminating the constraints - GaussianFactor::shared_ptr unconstrainedNewFactor; - if (hessians.size() > 0) { - bool hasSeparator = false; - GaussianFactorGraph::Keys unconstrainedKeys = hessians.keys(); - BOOST_FOREACH(Key key, unconstrainedKeys) { - if (find(keys.begin(), keys.end(), key) == keys.end()) { - hasSeparator = true; - break; - } - } - - if (hasSeparator) { - // find frontal keys in the unconstrained factor to eliminate - Ordering subkeys; - BOOST_FOREACH(Key key, keys) { - if (unconstrainedKeys.exists(key)) - subkeys.push_back(key); - } - GaussianConditional::shared_ptr unconstrainedConditional; - boost::tie(unconstrainedConditional, unconstrainedNewFactor) - = EliminateCholesky(hessians, subkeys); - constraints.push_back(unconstrainedConditional); - } - else { - if (hasHessians) { - HessianFactor::shared_ptr jointFactor = boost::make_shared< - HessianFactor>(hessians, Scatter(factors, keys)); - constraints.push_back(jointFactor); - } - else { - constraints.push_back(hessians); - } - } - } - - // Now eliminate the constraints - GaussianConditional::shared_ptr constrainedConditional; - GaussianFactor::shared_ptr constrainedNewFactor; - boost::tie(constrainedConditional, constrainedNewFactor) = EliminateQR( - constraints, keys); -// constraints.print("constraints: "); -// constrainedConditional->print("constrainedConditional: "); -// constrainedNewFactor->print("constrainedNewFactor: "); - - if (unconstrainedNewFactor) { - GaussianFactorGraph newFactors; - newFactors.push_back(unconstrainedNewFactor); - newFactors.push_back(constrainedNewFactor); -// newFactors.print("newFactors: "); - HessianFactor::shared_ptr newFactor(new HessianFactor(newFactors)); - return make_pair(constrainedConditional, newFactor); - } else { - return make_pair(constrainedConditional, constrainedNewFactor); - } - } - else { + if (hasConstraints(factors)) + return EliminateQR(factors, keys); + else return EliminateCholesky(factors, keys); - } } } // gtsam diff --git a/gtsam_unstable/nonlinear/SQPSimple.cpp b/gtsam_unstable/nonlinear/LCNLPSolver.cpp similarity index 60% rename from gtsam_unstable/nonlinear/SQPSimple.cpp rename to gtsam_unstable/nonlinear/LCNLPSolver.cpp index b0e3e8e5a..3a2f017f7 100644 --- a/gtsam_unstable/nonlinear/SQPSimple.cpp +++ b/gtsam_unstable/nonlinear/LCNLPSolver.cpp @@ -10,33 +10,32 @@ * -------------------------------------------------------------------------- */ /** - * @file SQPSimple.cpp + * @file LCNLPSolver.cpp * @author Duy-Nguyen Ta * @author Krunal Chande * @author Luca Carlone * @date Dec 15, 2014 */ -#include +#include #include namespace gtsam { /* ************************************************************************* */ -bool SQPSimple::isStationary(const VectorValues& delta) const { +bool LCNLPSolver::isStationary(const VectorValues& delta) const { return delta.vector().lpNorm() < errorTol; } /* ************************************************************************* */ -bool SQPSimple::isPrimalFeasible(const SQPSimpleState& state) const { - return nlp_.linearEqualities.checkFeasibility(state.values, errorTol) - && nlp_.nonlinearEqualities.checkFeasibility(state.values, errorTol); +bool LCNLPSolver::isPrimalFeasible(const LCNLPState& state) const { + return lcnlp_.linearEqualities.checkFeasibility(state.values, errorTol); } /* ************************************************************************* */ -bool SQPSimple::isDualFeasible(const VectorValues& duals) const { - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearInequalities) { +bool LCNLPSolver::isDualFeasible(const VectorValues& duals) const { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, lcnlp_.linearInequalities) { NonlinearConstraint::shared_ptr inequality = boost::dynamic_pointer_cast(factor); Key dualKey = inequality->dualKey(); if (!duals.exists(dualKey)) continue; // should be inactive constraint! @@ -48,33 +47,29 @@ bool SQPSimple::isDualFeasible(const VectorValues& duals) const { } /* ************************************************************************* */ -bool SQPSimple::isComplementary(const SQPSimpleState& state) const { - return nlp_.linearInequalities.checkFeasibilityAndComplimentary(state.values, state.duals, errorTol); +bool LCNLPSolver::isComplementary(const LCNLPState& state) const { + return lcnlp_.linearInequalities.checkFeasibilityAndComplimentary(state.values, state.duals, errorTol); } /* ************************************************************************* */ -bool SQPSimple::checkConvergence(const SQPSimpleState& state, const VectorValues& delta) const { +bool LCNLPSolver::checkConvergence(const LCNLPState& state, const VectorValues& delta) const { return isStationary(delta) && isPrimalFeasible(state) && isDualFeasible(state.duals) && isComplementary(state); } /* ************************************************************************* */ -VectorValues SQPSimple::initializeDuals() const { +VectorValues LCNLPSolver::initializeDuals() const { VectorValues duals; - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.linearEqualities) { + BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, lcnlp_.linearEqualities) { NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); duals.insert(constraint->dualKey(), zero(factor->dim())); } - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, nlp_.nonlinearEqualities) { - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); - duals.insert(constraint->dualKey(), zero(factor->dim())); - } return duals; } /* ************************************************************************* */ -std::pair SQPSimple::optimize(const Values& initialValues) const { - SQPSimpleState state(initialValues); +std::pair LCNLPSolver::optimize(const Values& initialValues) const { + LCNLPState state(initialValues); state.duals = initializeDuals(); while (!state.converged && state.iterations < 100) { state = iterate(state); @@ -83,19 +78,14 @@ std::pair SQPSimple::optimize(const Values& initialValues) } /* ************************************************************************* */ -SQPSimpleState SQPSimple::iterate(const SQPSimpleState& state) const { +LCNLPState LCNLPSolver::iterate(const LCNLPState& state) const { static const bool debug = false; // construct the qp subproblem QP qp; - qp.cost = *nlp_.cost.linearize(state.values); - GaussianFactorGraph::shared_ptr multipliedHessians = nlp_.nonlinearEqualities.multipliedHessians(state.values, state.duals); - qp.cost.push_back(*multipliedHessians); - - qp.equalities.add(*nlp_.linearEqualities.linearize(state.values)); - qp.equalities.add(*nlp_.nonlinearEqualities.linearize(state.values)); - - qp.inequalities.add(*nlp_.linearInequalities.linearize(state.values)); + qp.cost = *lcnlp_.cost.linearize(state.values); + qp.equalities.add(*lcnlp_.linearEqualities.linearize(state.values)); + qp.inequalities.add(*lcnlp_.linearInequalities.linearize(state.values)); if (debug) qp.print("QP subproblem:"); @@ -112,7 +102,7 @@ SQPSimpleState SQPSimple::iterate(const SQPSimpleState& state) const { duals.print("duals = "); // update new state - SQPSimpleState newState; + LCNLPState newState; newState.values = state.values.retract(delta); newState.duals = duals; newState.converged = checkConvergence(newState, delta); diff --git a/gtsam_unstable/nonlinear/SQPSimple.h b/gtsam_unstable/nonlinear/LCNLPSolver.h similarity index 73% rename from gtsam_unstable/nonlinear/SQPSimple.h rename to gtsam_unstable/nonlinear/LCNLPSolver.h index 42b5c451a..6fde771b6 100644 --- a/gtsam_unstable/nonlinear/SQPSimple.h +++ b/gtsam_unstable/nonlinear/LCNLPSolver.h @@ -10,7 +10,7 @@ * -------------------------------------------------------------------------- */ /** - * @file SQPSimple.h + * @file LCNLPSolver.h * @author Duy-Nguyen Ta * @author Krunal Chande * @author Luca Carlone @@ -24,46 +24,51 @@ namespace gtsam { -struct NLP { +/** + * Nonlinear Programming problem with + * only linear constraints, encoded in factor graphs + */ +struct LCNLP { NonlinearFactorGraph cost; NonlinearEqualityFactorGraph linearEqualities; - NonlinearEqualityFactorGraph nonlinearEqualities; NonlinearInequalityFactorGraph linearInequalities; }; -struct SQPSimpleState { +/** + * State of LCNLPSolver before/after each iteration + */ +struct LCNLPState { Values values; VectorValues duals; bool converged; size_t iterations; /// Default constructor - SQPSimpleState() : values(), duals(), converged(false), iterations(0) {} + LCNLPState() : values(), duals(), converged(false), iterations(0) {} /// Constructor with an initialValues - SQPSimpleState(const Values& initialValues) : + LCNLPState(const Values& initialValues) : values(initialValues), duals(VectorValues()), converged(false), iterations(0) { } }; /** - * Simple SQP optimizer to solve nonlinear constrained problems. - * This simple version won't care about nonconvexity, which needs - * more advanced techniques to solve, e.g., merit function, line search, second-order correction etc. + * Simple SQP optimizer to solve nonlinear constrained problems + * ONLY linear constraints are supported. */ -class SQPSimple { - NLP nlp_; +class LCNLPSolver { + LCNLP lcnlp_; static const double errorTol = 1e-5; public: - SQPSimple(const NLP& nlp) : - nlp_(nlp) { + LCNLPSolver(const LCNLP& lcnlp) : + lcnlp_(lcnlp) { } /// Check if \nabla f(x) - \lambda * \nabla c(x) == 0 bool isStationary(const VectorValues& delta) const; /// Check if c_E(x) == 0 - bool isPrimalFeasible(const SQPSimpleState& state) const; + bool isPrimalFeasible(const LCNLPState& state) const; /** * Dual variables of inequality constraints need to be >=0 @@ -81,15 +86,15 @@ public: * If it is inactive, the dual should be 0, regardless of the error. However, we don't compute * dual variables for inactive constraints in the QP subproblem, so we don't care. */ - bool isComplementary(const SQPSimpleState& state) const; + bool isComplementary(const LCNLPState& state) const; /// Check convergence - bool checkConvergence(const SQPSimpleState& state, const VectorValues& delta) const; + bool checkConvergence(const LCNLPState& state, const VectorValues& delta) const; /** * Single iteration of SQP */ - SQPSimpleState iterate(const SQPSimpleState& state) const; + LCNLPState iterate(const LCNLPState& state) const; VectorValues initializeDuals() const; /** diff --git a/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h b/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h index fbc52f376..b4dab147b 100644 --- a/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h +++ b/gtsam_unstable/nonlinear/NonlinearEqualityFactorGraph.h @@ -58,17 +58,17 @@ public: return true; } - /** - * Additional cost for -lambda*ConstraintHessian for SQP - */ - GaussianFactorGraph::shared_ptr multipliedHessians(const Values& values, const VectorValues& duals) const { - GaussianFactorGraph::shared_ptr constrainedHessians(new GaussianFactorGraph()); - BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this) { - NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); - constrainedHessians->push_back(constraint->multipliedHessian(values, duals)); - } - return constrainedHessians; - } +// /** +// * Additional cost for -lambda*ConstraintHessian for SQP +// */ +// GaussianFactorGraph::shared_ptr multipliedHessians(const Values& values, const VectorValues& duals) const { +// GaussianFactorGraph::shared_ptr constrainedHessians(new GaussianFactorGraph()); +// BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, *this) { +// NonlinearConstraint::shared_ptr constraint = boost::dynamic_pointer_cast(factor); +// constrainedHessians->push_back(constraint->multipliedHessian(values, duals)); +// } +// return constrainedHessians; +// } }; diff --git a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp b/gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp similarity index 63% rename from gtsam_unstable/nonlinear/tests/testSQPSimple.cpp rename to gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp index dc406acd8..b8b1fde25 100644 --- a/gtsam_unstable/nonlinear/tests/testSQPSimple.cpp +++ b/gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp @@ -26,7 +26,7 @@ #include #include -#include +#include #include #include #include @@ -54,7 +54,7 @@ public: } }; -TEST(testSQPSimple, QPProblem) { +TEST(testlcnlpSolver, QPProblem) { const Key dualKey = 0; // Simple quadratic cost: x1^2 + x2^2 @@ -81,85 +81,27 @@ TEST(testSQPSimple, QPProblem) { initialVectorValues.insert(Y(1), ones(1)); VectorValues expectedSolution = qpSolver.optimize(initialVectorValues).first; - //Instantiate NLP - NLP nlp; + //Instantiate LCNLP + LCNLP lcnlp; Values linPoint; linPoint.insert(X(1), zero(1)); linPoint.insert(Y(1), zero(1)); - nlp.cost.add(LinearContainerFactor(hf, linPoint)); // wrap it using linearcontainerfactor - nlp.linearEqualities.add(ConstraintProblem1(X(1), Y(1), dualKey)); + lcnlp.cost.add(LinearContainerFactor(hf, linPoint)); // wrap it using linearcontainerfactor + lcnlp.linearEqualities.add(ConstraintProblem1(X(1), Y(1), dualKey)); Values initialValues; initialValues.insert(X(1), 0.0); initialValues.insert(Y(1), 0.0); - // Instantiate SQP - SQPSimple sqpSimple(nlp); - Values actualValues = sqpSimple.optimize(initialValues).first; + // Instantiate LCNLPSolver + LCNLPSolver lcnlpSolver(lcnlp); + Values actualValues = lcnlpSolver.optimize(initialValues).first; DOUBLES_EQUAL(expectedSolution.at(X(1))[0], actualValues.at(X(1)), 1e-100); DOUBLES_EQUAL(expectedSolution.at(Y(1))[0], actualValues.at(Y(1)), 1e-100); } -//****************************************************************************** -class CircleConstraint : public NonlinearConstraint2 { - typedef NonlinearConstraint2 Base; -public: - CircleConstraint(Key xK, Key yK, Key dualKey) : Base(xK, yK, dualKey, 1) {} - - Vector evaluateError(const double& x, const double& y, - boost::optional H1 = boost::none, boost::optional H2 = - boost::none) const { - if (H1) *H1 = (Matrix(1,1) << 2*(x-1)).finished(); - if (H2) *H2 = (Matrix(1,1) << 2*y).finished(); - return (Vector(1) << (x-1)*(x-1) + y*y - 0.25).finished(); - } - - void evaluateHessians(const double& x, const double& y, - std::vector& G11, std::vector& G12, - std::vector& G22) const { - G11.push_back((Matrix(1,1) << 2).finished()); - G12.push_back((Matrix(1,1) << 0).finished()); - G22.push_back((Matrix(1,1) << 2).finished()); - } - -}; - -// Nonlinear constraint not supported currently -TEST_DISABLED(testSQPSimple, quadraticCostNonlinearConstraint) { - const Key dualKey = 0; - - //Instantiate NLP - NLP nlp; - - // Simple quadratic cost: x1^2 + x2^2 +1000 - // Note the Hessian encodes: - // 0.5*x1'*G11*x1 + x1'*G12*x2 + 0.5*x2'*G22*x2 - x1'*g1 - x2'*g2 + 0.5*f - // Hence here we have G11 = 2, G12 = 0, G22 = 2, g1 = 0, g2 = 0, f = 0 - Values linPoint; - linPoint.insert(X(1), zero(1)); - linPoint.insert(Y(1), zero(1)); - HessianFactor hf(X(1), Y(1), 2.0 * ones(1,1), zero(1), zero(1), - 2*ones(1,1), zero(1) , 1000); - nlp.cost.add(LinearContainerFactor(hf, linPoint)); // wrap it using linearcontainerfactor - nlp.nonlinearEqualities.add(CircleConstraint(X(1), Y(1), dualKey)); - - Values initialValues; - initialValues.insert(X(1), 4.0); - initialValues.insert(Y(1), 10.0); - - Values expectedSolution; - expectedSolution.insert(X(1), 0.5); - expectedSolution.insert(Y(1), 0.0); - - // Instantiate SQP - SQPSimple sqpSimple(nlp); - Values actualSolution = sqpSimple.optimize(initialValues).first; - - CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); -} - //****************************************************************************** class LineConstraintX : public NonlinearConstraint1 { typedef NonlinearConstraint1 Base; @@ -171,13 +113,6 @@ public: return pose.x(); } - void evaluateHessians(const Pose3& pose, std::vector& G11) const { - Matrix G11all = Z_6x6; - Vector rT1 = pose.rotation().matrix().row(0); - G11all.block<3,3>(3,0) = skewSymmetric(rT1); - G11.push_back(G11all); - } - Vector evaluateError(const Pose3& pose, boost::optional H = boost::none) const { if (H) *H = (Matrix(1,6) << zeros(1,3), pose.rotation().matrix().row(0)).finished(); @@ -185,15 +120,15 @@ public: } }; -TEST(testSQPSimple, poseOnALine) { +TEST(testlcnlpSolver, poseOnALine) { const Key dualKey = 0; - //Instantiate NLP - NLP nlp; - nlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(1, 0, 0)), noiseModel::Unit::Create(6))); + //Instantiate LCNLP + LCNLP lcnlp; + lcnlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(1, 0, 0)), noiseModel::Unit::Create(6))); LineConstraintX constraint(X(1), dualKey); - nlp.nonlinearEqualities.add(constraint); + lcnlp.linearEqualities.add(constraint); Values initialValues; initialValues.insert(X(1), Pose3(Rot3::ypr(0.3, 0.2, 0.3), Point3(1,0,0))); @@ -201,9 +136,9 @@ TEST(testSQPSimple, poseOnALine) { Values expectedSolution; expectedSolution.insert(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3())); - // Instantiate SQP - SQPSimple sqpSimple(nlp); - Values actualSolution = sqpSimple.optimize(initialValues).first; + // Instantiate LCNLPSolver + LCNLPSolver lcnlpSolver(lcnlp); + Values actualSolution = lcnlpSolver.optimize(initialValues).first; CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); Pose3 pose(Rot3::ypr(0.1, 0.2, 0.3), Point3()); @@ -226,7 +161,7 @@ public: } }; -TEST(testSQPSimple, inequalityConstraint) { +TEST(testlcnlpSolver, inequalityConstraint) { const Key dualKey = 0; // Simple quadratic cost: x^2 + y^2 @@ -253,27 +188,26 @@ TEST(testSQPSimple, inequalityConstraint) { initialVectorValues.insert(Y(1), zero(1)); VectorValues expectedSolution = qpSolver.optimize(initialVectorValues).first; - //Instantiate NLP - NLP nlp; + //Instantiate LCNLP + LCNLP lcnlp; Values linPoint; linPoint.insert(X(1), zero(1)); linPoint.insert(Y(1), zero(1)); - nlp.cost.add(LinearContainerFactor(hf, linPoint)); // wrap it using linearcontainerfactor - nlp.linearInequalities.add(InequalityProblem1(X(1), Y(1), dualKey)); + lcnlp.cost.add(LinearContainerFactor(hf, linPoint)); // wrap it using linearcontainerfactor + lcnlp.linearInequalities.add(InequalityProblem1(X(1), Y(1), dualKey)); Values initialValues; initialValues.insert(X(1), 1.0); initialValues.insert(Y(1), -10.0); - // Instantiate SQP - SQPSimple sqpSimple(nlp); - Values actualValues = sqpSimple.optimize(initialValues).first; + // Instantiate LCNLPSolver + LCNLPSolver lcnlpSolver(lcnlp); + Values actualValues = lcnlpSolver.optimize(initialValues).first; DOUBLES_EQUAL(expectedSolution.at(X(1))[0], actualValues.at(X(1)), 1e-10); DOUBLES_EQUAL(expectedSolution.at(Y(1))[0], actualValues.at(Y(1)), 1e-10); } - //****************************************************************************** const size_t X_AXIS = 0; const size_t Y_AXIS = 1; @@ -319,14 +253,14 @@ public: } }; -TEST(testSQPSimple, poseWithABoundary) { +TEST(testlcnlpSolver, poseWithABoundary) { const Key dualKey = 0; - //Instantiate NLP - NLP nlp; - nlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(1, 0, 0)), noiseModel::Unit::Create(6))); + //Instantiate LCNLP + LCNLP lcnlp; + lcnlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(1, 0, 0)), noiseModel::Unit::Create(6))); AxisUpperBound constraint(X(1), X_AXIS, 0, dualKey); - nlp.linearInequalities.add(constraint); + lcnlp.linearInequalities.add(constraint); Values initialValues; initialValues.insert(X(1), Pose3(Rot3::ypr(0.3, 0.2, 0.3), Point3(1, 0, 0))); @@ -334,23 +268,23 @@ TEST(testSQPSimple, poseWithABoundary) { Values expectedSolution; expectedSolution.insert(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(0, 0, 0))); - // Instantiate SQP - SQPSimple sqpSimple(nlp); - Values actualSolution = sqpSimple.optimize(initialValues).first; + // Instantiate LCNLPSolver + LCNLPSolver lcnlpSolver(lcnlp); + Values actualSolution = lcnlpSolver.optimize(initialValues).first; CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); } -TEST(testSQPSimple, poseWithinA2DBox) { +TEST(testlcnlpSolver, poseWithinA2DBox) { const Key dualKey = 0; - //Instantiate NLP - NLP nlp; - nlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(10, 0.5, 0)), noiseModel::Unit::Create(6))); - nlp.linearInequalities.add(AxisLowerBound(X(1), X_AXIS, -1, dualKey)); // -1 <= x - nlp.linearInequalities.add(AxisUpperBound(X(1), X_AXIS, 1, dualKey+1)); // x <= 1 - nlp.linearInequalities.add(AxisLowerBound(X(1), Y_AXIS, -1, dualKey+2)); // -1 <= y - nlp.linearInequalities.add(AxisUpperBound(X(1), Y_AXIS, 1, dualKey+3));// y <= 1 + //Instantiate LCNLP + LCNLP lcnlp; + lcnlp.cost.add(PriorFactor(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(10, 0.5, 0)), noiseModel::Unit::Create(6))); + lcnlp.linearInequalities.add(AxisLowerBound(X(1), X_AXIS, -1, dualKey)); // -1 <= x + lcnlp.linearInequalities.add(AxisUpperBound(X(1), X_AXIS, 1, dualKey+1)); // x <= 1 + lcnlp.linearInequalities.add(AxisLowerBound(X(1), Y_AXIS, -1, dualKey+2)); // -1 <= y + lcnlp.linearInequalities.add(AxisUpperBound(X(1), Y_AXIS, 1, dualKey+3));// y <= 1 Values initialValues; initialValues.insert(X(1), Pose3(Rot3::ypr(1, -1, 2), Point3(3, -5, 0))); @@ -358,14 +292,14 @@ TEST(testSQPSimple, poseWithinA2DBox) { Values expectedSolution; expectedSolution.insert(X(1), Pose3(Rot3::ypr(0.1, 0.2, 0.3), Point3(1, 0.5, 0))); - // Instantiate SQP - SQPSimple sqpSimple(nlp); - Values actualSolution = sqpSimple.optimize(initialValues).first; + // Instantiate LCNLPSolver + LCNLPSolver lcnlpSolver(lcnlp); + Values actualSolution = lcnlpSolver.optimize(initialValues).first; CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); } -TEST_DISABLED(testSQPSimple, posesInA2DBox) { +TEST_DISABLED(testlcnlpSolver, posesInA2DBox) { const double xLowerBound = -3.0, xUpperBound = 5.0, yLowerBound = -1.0, @@ -373,32 +307,32 @@ TEST_DISABLED(testSQPSimple, posesInA2DBox) { zLowerBound = 0.0, zUpperBound = 2.0; - //Instantiate NLP - NLP nlp; + //Instantiate LCNLP + LCNLP lcnlp; // prior on the first pose SharedDiagonal priorNoise = noiseModel::Diagonal::Sigmas( (Vector(6) << 0.001, 0.001, 0.001, 0.001, 0.001, 0.001).finished()); - nlp.cost.add(PriorFactor(X(1), Pose3(), priorNoise)); + lcnlp.cost.add(PriorFactor(X(1), Pose3(), priorNoise)); // odometry between factor for subsequent poses SharedDiagonal odoNoise = noiseModel::Diagonal::Sigmas( (Vector(6) << 0.01, 0.01, 0.01, 0.1, 0.1, 0.1).finished()); Pose3 odo12(Rot3::ypr(M_PI/2.0, 0, 0), Point3(10, 0, 0)); - nlp.cost.add(BetweenFactor(X(1), X(2), odo12, odoNoise)); + lcnlp.cost.add(BetweenFactor(X(1), X(2), odo12, odoNoise)); Pose3 odo23(Rot3::ypr(M_PI/2.0, 0, 0), Point3(2, 0, 2)); - nlp.cost.add(BetweenFactor(X(2), X(3), odo23, odoNoise)); + lcnlp.cost.add(BetweenFactor(X(2), X(3), odo23, odoNoise)); // Box constraints Key dualKey = 0; for (size_t i=1; i<=3; ++i) { - nlp.linearInequalities.add(AxisLowerBound(X(i), X_AXIS, xLowerBound, dualKey++)); - nlp.linearInequalities.add(AxisUpperBound(X(i), X_AXIS, xUpperBound, dualKey++)); - nlp.linearInequalities.add(AxisLowerBound(X(i), Y_AXIS, yLowerBound, dualKey++)); - nlp.linearInequalities.add(AxisUpperBound(X(i), Y_AXIS, yUpperBound, dualKey++)); - nlp.linearInequalities.add(AxisLowerBound(X(i), Z_AXIS, zLowerBound, dualKey++)); - nlp.linearInequalities.add(AxisUpperBound(X(i), Z_AXIS, zUpperBound, dualKey++)); + lcnlp.linearInequalities.add(AxisLowerBound(X(i), X_AXIS, xLowerBound, dualKey++)); + lcnlp.linearInequalities.add(AxisUpperBound(X(i), X_AXIS, xUpperBound, dualKey++)); + lcnlp.linearInequalities.add(AxisLowerBound(X(i), Y_AXIS, yLowerBound, dualKey++)); + lcnlp.linearInequalities.add(AxisUpperBound(X(i), Y_AXIS, yUpperBound, dualKey++)); + lcnlp.linearInequalities.add(AxisLowerBound(X(i), Z_AXIS, zLowerBound, dualKey++)); + lcnlp.linearInequalities.add(AxisUpperBound(X(i), Z_AXIS, zUpperBound, dualKey++)); } Values initialValues; @@ -412,9 +346,9 @@ TEST_DISABLED(testSQPSimple, posesInA2DBox) { expectedSolution.insert(X(2), Pose3()); expectedSolution.insert(X(3), Pose3()); - // Instantiate SQP - SQPSimple sqpSimple(nlp); - Values actualSolution = sqpSimple.optimize(initialValues).first; + // Instantiate LCNLPSolver + LCNLPSolver lcnlpSolver(lcnlp); + Values actualSolution = lcnlpSolver.optimize(initialValues).first; actualSolution.print("actual solution: "); From bcdeddbda177a983f621c9aa94e1c7837513eb08 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 16:09:18 -0500 Subject: [PATCH 023/130] Compute the least-square values of dual variables instead of forcing them to satisfy the stationarity condition exactly. This fixes an infinite-loop bug in QPSolver when a constraint was continuously added to and removed from the active set. --- gtsam_unstable/linear/QPSolver.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/gtsam_unstable/linear/QPSolver.cpp b/gtsam_unstable/linear/QPSolver.cpp index f0eb7d7fb..e5b4cd801 100644 --- a/gtsam_unstable/linear/QPSolver.cpp +++ b/gtsam_unstable/linear/QPSolver.cpp @@ -58,7 +58,7 @@ JacobianFactor::shared_ptr QPSolver::createDualFactor(Key key, b += factor->gradient(key, delta); } } - return boost::make_shared(Aterms, b, noiseModel::Constrained::All(b.rows())); + return boost::make_shared(Aterms, b); // compute the least-square approximation of dual variables } else { return boost::make_shared(); @@ -156,7 +156,7 @@ boost::tuple QPSolver::computeStepSize( //****************************************************************************** QPState QPSolver::iterate(const QPState& state) const { - static bool debug = false; + static bool debug = true; // Solve with the current working set VectorValues newValues = solveWithCurrentWorkingSet(state.workingSet); @@ -164,7 +164,7 @@ QPState QPSolver::iterate(const QPState& state) const { newValues.print("New solution:"); // If we CAN'T move further - if (newValues.equals(state.values, 1e-5)) { + if (newValues.equals(state.values, 1e-7)) { // Compute lambda from the dual graph if (debug) cout << "Building dual graph..." << endl; From 121feb2281080a31b23966c8396b34eec5b77ef1 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 16:14:07 -0500 Subject: [PATCH 024/130] disable printing --- gtsam_unstable/linear/QPSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gtsam_unstable/linear/QPSolver.cpp b/gtsam_unstable/linear/QPSolver.cpp index e5b4cd801..53d8d5c4c 100644 --- a/gtsam_unstable/linear/QPSolver.cpp +++ b/gtsam_unstable/linear/QPSolver.cpp @@ -156,7 +156,7 @@ boost::tuple QPSolver::computeStepSize( //****************************************************************************** QPState QPSolver::iterate(const QPState& state) const { - static bool debug = true; + static bool debug = false; // Solve with the current working set VectorValues newValues = solveWithCurrentWorkingSet(state.workingSet); From 3653e9333819af75fdbac093624fca7232577adf Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 16:14:50 -0500 Subject: [PATCH 025/130] fix bug in dual feasibility check. We want dual variables < 0. --- gtsam_unstable/nonlinear/LCNLPSolver.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gtsam_unstable/nonlinear/LCNLPSolver.cpp b/gtsam_unstable/nonlinear/LCNLPSolver.cpp index 3a2f017f7..0fa840ab8 100644 --- a/gtsam_unstable/nonlinear/LCNLPSolver.cpp +++ b/gtsam_unstable/nonlinear/LCNLPSolver.cpp @@ -40,8 +40,9 @@ bool LCNLPSolver::isDualFeasible(const VectorValues& duals) const { Key dualKey = inequality->dualKey(); if (!duals.exists(dualKey)) continue; // should be inactive constraint! double dual = duals.at(dualKey)[0]; // because we only support single-valued inequalities - if (dual < 0.0) + if (dual > 0.0) { // See the explanation in QPSolver::identifyLeavingConstraint, we want dual < 0 ? return false; + } } return true; } From e74b737a66987d254e840baa431483d35189cbe1 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 16:17:14 -0500 Subject: [PATCH 026/130] box-constraints unit test passed with prior and between factors --- .../nonlinear/tests/testLCNLPSolver.cpp | 21 ++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp b/gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp index b8b1fde25..1a1053147 100644 --- a/gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp +++ b/gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp @@ -299,7 +299,7 @@ TEST(testlcnlpSolver, poseWithinA2DBox) { CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); } -TEST_DISABLED(testlcnlpSolver, posesInA2DBox) { +TEST(testlcnlpSolver, posesInA2DBox) { const double xLowerBound = -3.0, xUpperBound = 5.0, yLowerBound = -1.0, @@ -312,12 +312,12 @@ TEST_DISABLED(testlcnlpSolver, posesInA2DBox) { // prior on the first pose SharedDiagonal priorNoise = noiseModel::Diagonal::Sigmas( - (Vector(6) << 0.001, 0.001, 0.001, 0.001, 0.001, 0.001).finished()); + (Vector(6) << 0.001, 0.001, 0.001, 0.0001, 0.0001, 0.0001).finished()); lcnlp.cost.add(PriorFactor(X(1), Pose3(), priorNoise)); // odometry between factor for subsequent poses SharedDiagonal odoNoise = noiseModel::Diagonal::Sigmas( - (Vector(6) << 0.01, 0.01, 0.01, 0.1, 0.1, 0.1).finished()); + (Vector(6) << 0.001, 0.001, 0.001, 0.1, 0.1, 0.1).finished()); Pose3 odo12(Rot3::ypr(M_PI/2.0, 0, 0), Point3(10, 0, 0)); lcnlp.cost.add(BetweenFactor(X(1), X(2), odo12, odoNoise)); @@ -343,16 +343,23 @@ TEST_DISABLED(testlcnlpSolver, posesInA2DBox) { Values expectedSolution; expectedSolution.insert(X(1), Pose3()); - expectedSolution.insert(X(2), Pose3()); - expectedSolution.insert(X(3), Pose3()); + expectedSolution.insert(X(2), Pose3(Rot3::ypr(M_PI/2.0, 0, 0), Point3(5, 0, 0))); + expectedSolution.insert(X(3), Pose3(Rot3::ypr(M_PI, 0, 0), Point3(5, 2, 2))); // Instantiate LCNLPSolver LCNLPSolver lcnlpSolver(lcnlp); Values actualSolution = lcnlpSolver.optimize(initialValues).first; - actualSolution.print("actual solution: "); +// cout << "Rotation angles: " << endl; +// for (size_t i = 1; i<=3; i++) { +// cout << actualSolution.at(X(i)).rotation().ypr().transpose()*180/M_PI << endl; +// } - CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); +// cout << "Actual Error: " << lcnlp.cost.error(actualSolution) << endl; +// cout << "Expected Error: " << lcnlp.cost.error(expectedSolution) << endl; +// actualSolution.print("actualSolution: "); + + CHECK(assert_equal(expectedSolution, actualSolution, 1e-5)); } //****************************************************************************** From 2791f2930bf9c3b1adc113e7788d79254ecf1a84 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 18:12:24 -0500 Subject: [PATCH 027/130] fix return value from key to size_t. --- gtsam/linear/VectorValues.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gtsam/linear/VectorValues.h b/gtsam/linear/VectorValues.h index ce33116ab..a599e7269 100644 --- a/gtsam/linear/VectorValues.h +++ b/gtsam/linear/VectorValues.h @@ -135,7 +135,7 @@ namespace gtsam { /// @{ /** Number of variables stored. */ - Key size() const { return values_.size(); } + size_t size() const { return values_.size(); } /** Return the dimension of variable \c j. */ size_t dim(Key j) const { return at(j).rows(); } From f4a4ce43258d7ec773f7d129a7394cf818c5137e Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 18:15:54 -0500 Subject: [PATCH 028/130] added debug print statements --- gtsam_unstable/nonlinear/LCNLPSolver.h | 10 ++++++++++ gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp | 13 ++++++++----- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/gtsam_unstable/nonlinear/LCNLPSolver.h b/gtsam_unstable/nonlinear/LCNLPSolver.h index 6fde771b6..a31a04928 100644 --- a/gtsam_unstable/nonlinear/LCNLPSolver.h +++ b/gtsam_unstable/nonlinear/LCNLPSolver.h @@ -50,6 +50,16 @@ struct LCNLPState { LCNLPState(const Values& initialValues) : values(initialValues), duals(VectorValues()), converged(false), iterations(0) { } + + /// print + void print(const std::string& s = "") const { + std::cout << s << std::endl; + values.print("Values: "); + duals.print("Duals: "); + if (converged) std::cout << "Converged!" << std::endl; + else std::cout << "Not converged" << std::endl; + std::cout << "Iterations: " << iterations << std::endl; + } }; /** diff --git a/gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp b/gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp index 1a1053147..edfc3bee6 100644 --- a/gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp +++ b/gtsam_unstable/nonlinear/tests/testLCNLPSolver.cpp @@ -54,7 +54,7 @@ public: } }; -TEST(testlcnlpSolver, QPProblem) { +TEST_DISABLED(testlcnlpSolver, QPProblem) { const Key dualKey = 0; // Simple quadratic cost: x1^2 + x2^2 @@ -120,7 +120,7 @@ public: } }; -TEST(testlcnlpSolver, poseOnALine) { +TEST_DISABLED(testlcnlpSolver, poseOnALine) { const Key dualKey = 0; @@ -161,7 +161,7 @@ public: } }; -TEST(testlcnlpSolver, inequalityConstraint) { +TEST_DISABLED(testlcnlpSolver, inequalityConstraint) { const Key dualKey = 0; // Simple quadratic cost: x^2 + y^2 @@ -253,7 +253,7 @@ public: } }; -TEST(testlcnlpSolver, poseWithABoundary) { +TEST_DISABLED(testlcnlpSolver, poseWithABoundary) { const Key dualKey = 0; //Instantiate LCNLP @@ -275,7 +275,7 @@ TEST(testlcnlpSolver, poseWithABoundary) { CHECK(assert_equal(expectedSolution, actualSolution, 1e-10)); } -TEST(testlcnlpSolver, poseWithinA2DBox) { +TEST_DISABLED(testlcnlpSolver, poseWithinA2DBox) { const Key dualKey = 0; //Instantiate LCNLP @@ -359,6 +359,9 @@ TEST(testlcnlpSolver, posesInA2DBox) { // cout << "Expected Error: " << lcnlp.cost.error(expectedSolution) << endl; // actualSolution.print("actualSolution: "); + AxisLowerBound factor(X(1), X_AXIS, xLowerBound, dualKey++); + Matrix hessian = numericalHessian(boost::bind(&AxisLowerBound::computeError, factor, _1, boost::none), Pose3(), 1e-3); + cout << "Hessian: \n" << hessian << endl; CHECK(assert_equal(expectedSolution, actualSolution, 1e-5)); } From 3e352f109e096152ff046c3444f123d17eb9f081 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Tue, 23 Dec 2014 18:16:22 -0500 Subject: [PATCH 029/130] Added warm start for initializing active set. --- gtsam_unstable/linear/QPSolver.cpp | 36 +++++++++++++++--------- gtsam_unstable/linear/QPSolver.h | 12 ++++---- gtsam_unstable/nonlinear/LCNLPSolver.cpp | 27 ++++++++++++++---- 3 files changed, 51 insertions(+), 24 deletions(-) diff --git a/gtsam_unstable/linear/QPSolver.cpp b/gtsam_unstable/linear/QPSolver.cpp index 53d8d5c4c..10d248d29 100644 --- a/gtsam_unstable/linear/QPSolver.cpp +++ b/gtsam_unstable/linear/QPSolver.cpp @@ -181,13 +181,13 @@ QPState QPSolver::iterate(const QPState& state) const { // If all inequality constraints are satisfied: We have the solution!! if (leavingFactor < 0) { - return QPState(newValues, duals, state.workingSet, true); + return QPState(newValues, duals, state.workingSet, true, state.iterations+1); } else { // Inactivate the leaving constraint LinearInequalityFactorGraph newWorkingSet = state.workingSet; newWorkingSet.at(leavingFactor)->inactivate(); - return QPState(newValues, duals, newWorkingSet, false); + return QPState(newValues, duals, newWorkingSet, false, state.iterations+1); } } else { @@ -210,23 +210,33 @@ QPState QPSolver::iterate(const QPState& state) const { // step! newValues = state.values + alpha * p; - return QPState(newValues, state.duals, newWorkingSet, false); + return QPState(newValues, state.duals, newWorkingSet, false, state.iterations+1); } } //****************************************************************************** LinearInequalityFactorGraph QPSolver::identifyActiveConstraints( const LinearInequalityFactorGraph& inequalities, - const VectorValues& initialValues) const { + const VectorValues& initialValues, const VectorValues& duals) const { LinearInequalityFactorGraph workingSet; - BOOST_FOREACH(const LinearInequality::shared_ptr& factor, inequalities){ + BOOST_FOREACH(const LinearInequality::shared_ptr& factor, inequalities) { LinearInequality::shared_ptr workingFactor(new LinearInequality(*factor)); - double error = workingFactor->error(initialValues); - if (fabs(error)>1e-7){ - workingFactor->inactivate(); - } else { + if (duals.exists(workingFactor->dualKey())) { workingFactor->activate(); } + else { + if (duals.size() > 0) { + workingFactor->inactivate(); + } else { + double error = workingFactor->error(initialValues); + if (fabs(error)<1e-7) { + workingFactor->activate(); + } + else { + workingFactor->inactivate(); + } + } + } workingSet.push_back(workingFactor); } return workingSet; @@ -234,18 +244,18 @@ LinearInequalityFactorGraph QPSolver::identifyActiveConstraints( //****************************************************************************** pair QPSolver::optimize( - const VectorValues& initialValues) const { + const VectorValues& initialValues, const VectorValues& duals) const { // Initialize workingSet from the feasible initialValues LinearInequalityFactorGraph workingSet = - identifyActiveConstraints(qp_.inequalities, initialValues); - QPState state(initialValues, VectorValues(), workingSet, false); + identifyActiveConstraints(qp_.inequalities, initialValues, duals); + QPState state(initialValues, duals, workingSet, false, 0); /// main loop of the solver while (!state.converged) { state = iterate(state); } - + std::cout << "Number of inner iterations: " << state.iterations << std::endl; return make_pair(state.values, state.duals); } diff --git a/gtsam_unstable/linear/QPSolver.h b/gtsam_unstable/linear/QPSolver.h index 53cb6ca6d..52a6f2f37 100644 --- a/gtsam_unstable/linear/QPSolver.h +++ b/gtsam_unstable/linear/QPSolver.h @@ -31,17 +31,18 @@ struct QPState { VectorValues duals; LinearInequalityFactorGraph workingSet; bool converged; + size_t iterations; /// default constructor QPState() : - values(), duals(), workingSet(), converged(false) { + values(), duals(), workingSet(), converged(false), iterations(0) { } /// constructor with initial values QPState(const VectorValues& initialValues, const VectorValues& initialDuals, - const LinearInequalityFactorGraph& initialWorkingSet, bool _converged) : + const LinearInequalityFactorGraph& initialWorkingSet, bool _converged, size_t _iterations) : values(initialValues), duals(initialDuals), workingSet(initialWorkingSet), converged( - _converged) { + _converged), iterations(_iterations) { } }; @@ -183,7 +184,8 @@ public: */ LinearInequalityFactorGraph identifyActiveConstraints( const LinearInequalityFactorGraph& inequalities, - const VectorValues& initialValues) const; + const VectorValues& initialValues, + const VectorValues& duals = VectorValues()) const; /** Optimize with a provided initial values * For this version, it is the responsibility of the caller to provide @@ -191,7 +193,7 @@ public: * @return a pair of solutions */ std::pair optimize( - const VectorValues& initialValues) const; + const VectorValues& initialValues, const VectorValues& duals = VectorValues()) const; }; diff --git a/gtsam_unstable/nonlinear/LCNLPSolver.cpp b/gtsam_unstable/nonlinear/LCNLPSolver.cpp index 0fa840ab8..720ee569d 100644 --- a/gtsam_unstable/nonlinear/LCNLPSolver.cpp +++ b/gtsam_unstable/nonlinear/LCNLPSolver.cpp @@ -19,6 +19,7 @@ #include #include +#include namespace gtsam { @@ -75,12 +76,13 @@ std::pair LCNLPSolver::optimize(const Values& initialValue while (!state.converged && state.iterations < 100) { state = iterate(state); } + std::cout << "Number of iterations: " << state.iterations << std::endl; return std::make_pair(state.values, state.duals); } /* ************************************************************************* */ LCNLPState LCNLPSolver::iterate(const LCNLPState& state) const { - static const bool debug = false; + static const bool debug = true; // construct the qp subproblem QP qp; @@ -88,19 +90,29 @@ LCNLPState LCNLPSolver::iterate(const LCNLPState& state) const { qp.equalities.add(*lcnlp_.linearEqualities.linearize(state.values)); qp.inequalities.add(*lcnlp_.linearInequalities.linearize(state.values)); - if (debug) - qp.print("QP subproblem:"); +// if (debug) +// qp.print("QP subproblem:"); // solve the QP subproblem VectorValues delta, duals; QPSolver qpSolver(qp); - boost::tie(delta, duals) = qpSolver.optimize(); + if (state.iterations == 0) + boost::tie(delta, duals) = qpSolver.optimize(); + else { + VectorValues zeroInitialValues; + BOOST_FOREACH(const Values::ConstKeyValuePair& key_value, state.values) { + zeroInitialValues.insert(key_value.key, zero(key_value.value.dim())); + } + boost::tie(delta, duals) = qpSolver.optimize(zeroInitialValues, state.duals); + } + if (debug) + state.values.print("state.values: "); if (debug) delta.print("delta = "); - if (debug) - duals.print("duals = "); +// if (debug) +// duals.print("duals = "); // update new state LCNLPState newState; @@ -109,6 +121,9 @@ LCNLPState LCNLPSolver::iterate(const LCNLPState& state) const { newState.converged = checkConvergence(newState, delta); newState.iterations = state.iterations + 1; + if (debug) + newState.print("newState: "); + return newState; } } From 2476bbe8d7e62fc22f1c964978538b1d0094493e Mon Sep 17 00:00:00 2001 From: krunalchande Date: Thu, 25 Dec 2014 00:00:22 -0500 Subject: [PATCH 030/130] can pass debug flag as parameter. --- gtsam_unstable/nonlinear/LCNLPSolver.cpp | 10 ++++++---- gtsam_unstable/nonlinear/LCNLPSolver.h | 6 +++--- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/gtsam_unstable/nonlinear/LCNLPSolver.cpp b/gtsam_unstable/nonlinear/LCNLPSolver.cpp index 720ee569d..99022ff5d 100644 --- a/gtsam_unstable/nonlinear/LCNLPSolver.cpp +++ b/gtsam_unstable/nonlinear/LCNLPSolver.cpp @@ -70,19 +70,21 @@ VectorValues LCNLPSolver::initializeDuals() const { } /* ************************************************************************* */ -std::pair LCNLPSolver::optimize(const Values& initialValues) const { +std::pair LCNLPSolver::optimize(const Values& initialValues, bool debug) const { LCNLPState state(initialValues); state.duals = initializeDuals(); while (!state.converged && state.iterations < 100) { - state = iterate(state); + if (debug) + std::cout << "state: iteration " << state.iterations << std::endl; + state = iterate(state, debug); } + if (debug) std::cout << "Number of iterations: " << state.iterations << std::endl; return std::make_pair(state.values, state.duals); } /* ************************************************************************* */ -LCNLPState LCNLPSolver::iterate(const LCNLPState& state) const { - static const bool debug = true; +LCNLPState LCNLPSolver::iterate(const LCNLPState& state, bool debug) const { // construct the qp subproblem QP qp; diff --git a/gtsam_unstable/nonlinear/LCNLPSolver.h b/gtsam_unstable/nonlinear/LCNLPSolver.h index a31a04928..047fca38b 100644 --- a/gtsam_unstable/nonlinear/LCNLPSolver.h +++ b/gtsam_unstable/nonlinear/LCNLPSolver.h @@ -104,13 +104,13 @@ public: /** * Single iteration of SQP */ - LCNLPState iterate(const LCNLPState& state) const; + LCNLPState iterate(const LCNLPState& state, bool debug = false) const; VectorValues initializeDuals() const; /** - * Main optimization function. + * Main optimization function. new */ - std::pair optimize(const Values& initialValues) const; + std::pair optimize(const Values& initialValues, bool debug = false) const; }; } From 3142f0a9a7fa47055c8a786a3e4b8be860a74d5e Mon Sep 17 00:00:00 2001 From: krunalchande Date: Thu, 25 Dec 2014 01:30:02 -0500 Subject: [PATCH 031/130] disabled test --- .../tests/testLinearInequalityManifolds.cpp | 21 +- ...tLinearlyConstrainedNonlinearOptimizer.cpp | 211 +++++++++--------- 2 files changed, 118 insertions(+), 114 deletions(-) diff --git a/gtsam_unstable/nonlinear/tests/testLinearInequalityManifolds.cpp b/gtsam_unstable/nonlinear/tests/testLinearInequalityManifolds.cpp index 8b1077ebf..3d6c401b0 100644 --- a/gtsam_unstable/nonlinear/tests/testLinearInequalityManifolds.cpp +++ b/gtsam_unstable/nonlinear/tests/testLinearInequalityManifolds.cpp @@ -48,18 +48,21 @@ // CHECK(!le1.equals(le2)); //} // -//////****************************************************************************** -////TEST(testLinearEqualityManifolds, linearize ) { -////} -//// -//////****************************************************************************** -////TEST(testLinearEqualityManifolds, size ) { -////} +////****************************************************************************** +//TEST(testLinearEqualityManifolds, linearize ) { +//} +// +////****************************************************************************** +//TEST(testLinearEqualityManifolds, size ) { +//} // ////****************************************************************************** //int main() { // TestResult tr; // return TestRegistry::runAllTests(tr); //} -////****************************************************************************** -// +//****************************************************************************** + +int main() { + +} diff --git a/gtsam_unstable/nonlinear/tests/testLinearlyConstrainedNonlinearOptimizer.cpp b/gtsam_unstable/nonlinear/tests/testLinearlyConstrainedNonlinearOptimizer.cpp index 78b511e2c..dd716868e 100644 --- a/gtsam_unstable/nonlinear/tests/testLinearlyConstrainedNonlinearOptimizer.cpp +++ b/gtsam_unstable/nonlinear/tests/testLinearlyConstrainedNonlinearOptimizer.cpp @@ -26,112 +26,113 @@ #include -namespace gtsam { -struct LinearlyConstrainedNLP { - NonlinearFactorGraph cost; - LinearEqualityFactorGraph equalities; - LinearInequalityFactorGraph inequalities; -}; - -struct LinearlyConstrainedNLPState { - Values values; - VectorValues duals; - bool converged; - LinearlyConstrainedNLPState(const Values& initialValues) : - values(initialValues), duals(VectorValues()), converged(false) { - } -}; -class LinearlyConstrainedNonLinearOptimizer { - LinearlyConstrainedNLP lcNLP_; -public: - LinearlyConstrainedNonLinearOptimizer(const LinearlyConstrainedNLP& lcNLP): lcNLP_(lcNLP) {} - - LinearlyConstrainedNLPState iterate(const LinearlyConstrainedNLPState& state) const { - QP qp; - qp.cost = lcNLP_.cost.linearize(state.values); - qp.equalities = lcNLP_.equalities; - qp.inequalities = lcNLP_.inequalities; - QPSolver qpSolver(qp); - VectorValues delta, duals; - boost::tie(delta, duals) = qpSolver.optimize(); - LinearlyConstrainedNLPState newState; - newState.values = state.values.retract(delta); - newState.duals = duals; - newState.converged = checkConvergence(newState.values, newState.duals); - return newState; - } - - /** - * Main optimization function. - */ - std::pair optimize(const Values& initialValues) const { - LinearlyConstrainedNLPState state(initialValues); - while(!state.converged){ - state = iterate(state); - } - - return std::make_pair(initialValues, VectorValues()); - } -}; -} - -using namespace std; -using namespace gtsam::symbol_shorthand; -using namespace gtsam; -const double tol = 1e-10; -//****************************************************************************** -TEST(LinearlyConstrainedNonlinearOptimizer, Problem1 ) { - - // build a quadratic Objective function x1^2 - x1*x2 + x2^2 - 3*x1 + 5 - // Note the Hessian encodes: - // 0.5*x1'*G11*x1 + x1'*G12*x2 + 0.5*x2'*G22*x2 - x1'*g1 - x2'*g2 + 0.5*f - // Hence, we have G11=2, G12 = -1, g1 = +3, G22 = 2, g2 = 0, f = 10 - HessianFactor lf(X(1), X(2), 2.0 * ones(1, 1), -ones(1, 1), 3.0 * ones(1), - 2.0 * ones(1, 1), zero(1), 10.0); - - // build linear inequalities - LinearInequalityFactorGraph inequalities; - inequalities.push_back(LinearInequality(X(1), ones(1,1), X(2), ones(1,1), 2, 0)); // x1 + x2 <= 2 --> x1 + x2 -2 <= 0, --> b=2 - inequalities.push_back(LinearInequality(X(1), -ones(1,1), 0, 1)); // -x1 <= 0 - inequalities.push_back(LinearInequality(X(2), -ones(1,1), 0, 2)); // -x2 <= 0 - inequalities.push_back(LinearInequality(X(1), ones(1,1), 1.5, 3)); // x1 <= 3/2 - - // Instantiate LinearlyConstrainedNLP, pretending that the cost is not quadratic - // (LinearContainerFactor makes a linear factor behave like a nonlinear one) - LinearlyConstrainedNLP lcNLP; - lcNLP.cost.add(LinearContainerFactor(lf)); - lcNLP.inequalities = inequalities; - - // Compare against a QP - QP qp; - qp.cost.add(lf); - qp.inequalities = inequalities; - - // instantiate QPsolver - QPSolver qpSolver(qp); - // create initial values for optimization - VectorValues initialVectorValues; - initialVectorValues.insert(X(1), zero(1)); - initialVectorValues.insert(X(2), zero(1)); - VectorValues expectedSolution = qpSolver.optimize(initialVectorValues).first; - - // instantiate LinearlyConstrainedNonLinearOptimizer - LinearlyConstrainedNonLinearOptimizer lcNLPSolver(lcNLP); - // create initial values for optimization - Values initialValues; - initialValues.insert(X(1), 0.0); - initialValues.insert(X(2), 0.0); - Values actualSolution = lcNLPSolver.optimize(initialValues).first; - - - DOUBLES_EQUAL(expectedSolution.at(X(1))[0], actualSolution.at(X(1)), tol); - DOUBLES_EQUAL(expectedSolution.at(X(2))[0], actualSolution.at(X(2)), tol); -} - +//namespace gtsam { +//struct LinearlyConstrainedNLP { +// NonlinearFactorGraph cost; +// LinearEqualityFactorGraph equalities; +// LinearInequalityFactorGraph inequalities; +//}; +// +//struct LinearlyConstrainedNLPState { +// Values values; +// VectorValues duals; +// bool converged; +// LinearlyConstrainedNLPState(const Values& initialValues) : +// values(initialValues), duals(VectorValues()), converged(false) { +// } +//}; +//class LinearlyConstrainedNonLinearOptimizer { +// LinearlyConstrainedNLP lcNLP_; +//public: +// LinearlyConstrainedNonLinearOptimizer(const LinearlyConstrainedNLP& lcNLP): lcNLP_(lcNLP) {} +// +// LinearlyConstrainedNLPState iterate(const LinearlyConstrainedNLPState& state) const { +// QP qp; +// qp.cost = lcNLP_.cost.linearize(state.values); +// qp.equalities = lcNLP_.equalities; +// qp.inequalities = lcNLP_.inequalities; +// QPSolver qpSolver(qp); +// VectorValues delta, duals; +// boost::tie(delta, duals) = qpSolver.optimize(); +// LinearlyConstrainedNLPState newState; +// newState.values = state.values.retract(delta); +// newState.duals = duals; +// newState.converged = checkConvergence(newState.values, newState.duals); +// return newState; +// } +// +// /** +// * Main optimization function. +// */ +// std::pair optimize(const Values& initialValues) const { +// LinearlyConstrainedNLPState state(initialValues); +// while(!state.converged){ +// state = iterate(state); +// } +// +// return std::make_pair(initialValues, VectorValues()); +// } +//}; +//} +// +//using namespace std; +//using namespace gtsam::symbol_shorthand; +//using namespace gtsam; +//const double tol = 1e-10; +////****************************************************************************** +//TEST(LinearlyConstrainedNonlinearOptimizer, Problem1 ) { +// +// // build a quadratic Objective function x1^2 - x1*x2 + x2^2 - 3*x1 + 5 +// // Note the Hessian encodes: +// // 0.5*x1'*G11*x1 + x1'*G12*x2 + 0.5*x2'*G22*x2 - x1'*g1 - x2'*g2 + 0.5*f +// // Hence, we have G11=2, G12 = -1, g1 = +3, G22 = 2, g2 = 0, f = 10 +// HessianFactor lf(X(1), X(2), 2.0 * ones(1, 1), -ones(1, 1), 3.0 * ones(1), +// 2.0 * ones(1, 1), zero(1), 10.0); +// +// // build linear inequalities +// LinearInequalityFactorGraph inequalities; +// inequalities.push_back(LinearInequality(X(1), ones(1,1), X(2), ones(1,1), 2, 0)); // x1 + x2 <= 2 --> x1 + x2 -2 <= 0, --> b=2 +// inequalities.push_back(LinearInequality(X(1), -ones(1,1), 0, 1)); // -x1 <= 0 +// inequalities.push_back(LinearInequality(X(2), -ones(1,1), 0, 2)); // -x2 <= 0 +// inequalities.push_back(LinearInequality(X(1), ones(1,1), 1.5, 3)); // x1 <= 3/2 +// +// // Instantiate LinearlyConstrainedNLP, pretending that the cost is not quadratic +// // (LinearContainerFactor makes a linear factor behave like a nonlinear one) +// LinearlyConstrainedNLP lcNLP; +// lcNLP.cost.add(LinearContainerFactor(lf)); +// lcNLP.inequalities = inequalities; +// +// // Compare against a QP +// QP qp; +// qp.cost.add(lf); +// qp.inequalities = inequalities; +// +// // instantiate QPsolver +// QPSolver qpSolver(qp); +// // create initial values for optimization +// VectorValues initialVectorValues; +// initialVectorValues.insert(X(1), zero(1)); +// initialVectorValues.insert(X(2), zero(1)); +// VectorValues expectedSolution = qpSolver.optimize(initialVectorValues).first; +// +// // instantiate LinearlyConstrainedNonLinearOptimizer +// LinearlyConstrainedNonLinearOptimizer lcNLPSolver(lcNLP); +// // create initial values for optimization +// Values initialValues; +// initialValues.insert(X(1), 0.0); +// initialValues.insert(X(2), 0.0); +// Values actualSolution = lcNLPSolver.optimize(initialValues).first; +// +// +// DOUBLES_EQUAL(expectedSolution.at(X(1))[0], actualSolution.at(X(1)), tol); +// DOUBLES_EQUAL(expectedSolution.at(X(2))[0], actualSolution.at(X(2)), tol); +//} +// //****************************************************************************** int main() { - TestResult tr; - return TestRegistry::runAllTests(tr); + std::cout<<"here"< Date: Thu, 25 Dec 2014 01:32:31 -0500 Subject: [PATCH 032/130] updated with current constructor whch needs number of iterations --- gtsam_unstable/linear/tests/testQPSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gtsam_unstable/linear/tests/testQPSolver.cpp b/gtsam_unstable/linear/tests/testQPSolver.cpp index 8fca61ca4..7898b4391 100644 --- a/gtsam_unstable/linear/tests/testQPSolver.cpp +++ b/gtsam_unstable/linear/tests/testQPSolver.cpp @@ -176,7 +176,7 @@ TEST(QPSolver, iterate) { LinearInequalityFactorGraph workingSet = solver.identifyActiveConstraints(qp.inequalities, currentSolution); - QPState state(currentSolution, VectorValues(), workingSet, false); + QPState state(currentSolution, VectorValues(), workingSet, false, 100); int it = 0; while (!state.converged) { From 37fe405872ba56fe2c1609f4a42362a6e63b8b16 Mon Sep 17 00:00:00 2001 From: krunalchande Date: Sat, 3 Jan 2015 12:32:49 -0500 Subject: [PATCH 033/130] Added warmStart flag. --- .cproject | 3284 ++++++++--------- gtsam_unstable/linear/QPSolver.cpp | 10 +- gtsam_unstable/linear/QPSolver.h | 4 +- gtsam_unstable/nonlinear/LCNLPSolver.cpp | 10 +- gtsam_unstable/nonlinear/LCNLPSolver.h | 4 +- .../nonlinear/tests/testLCNLPSolver.cpp | 4 +- 6 files changed, 1639 insertions(+), 1677 deletions(-) diff --git a/.cproject b/.cproject index ab1d0cdfc..77ccd8a14 100644 --- a/.cproject +++ b/.cproject @@ -5,13 +5,13 @@ - - + + @@ -60,13 +60,13 @@ - - + + @@ -116,13 +116,13 @@ - - + + @@ -130,7 +130,7 @@ - +