417 lines
		
	
	
		
			14 KiB
		
	
	
	
		
			C++
		
	
	
			
		
		
	
	
			417 lines
		
	
	
		
			14 KiB
		
	
	
	
		
			C++
		
	
	
/* ----------------------------------------------------------------------------
 | 
						|
 | 
						|
 * 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    ConcurrentBatchSmoother.cpp
 | 
						|
 * @brief   A Levenberg-Marquardt Batch Smoother that implements the
 | 
						|
 *          Concurrent Filtering and Smoothing interface.
 | 
						|
 * @author  Stephen Williams
 | 
						|
 */
 | 
						|
 | 
						|
#include <gtsam_unstable/nonlinear/ConcurrentBatchSmoother.h>
 | 
						|
#include <gtsam/nonlinear/LinearContainerFactor.h>
 | 
						|
#include <gtsam/linear/GaussianJunctionTree.h>
 | 
						|
#include <gtsam/base/timing.h>
 | 
						|
#include <gtsam/base/debug.h>
 | 
						|
 | 
						|
namespace gtsam {
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
void ConcurrentBatchSmoother::print(const std::string& s, const KeyFormatter& keyFormatter) const {
 | 
						|
  std::cout << s;
 | 
						|
  std::cout << "  Factors:" << std::endl;
 | 
						|
  BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, factors_) {
 | 
						|
    PrintNonlinearFactor(factor, "    ", keyFormatter);
 | 
						|
  }
 | 
						|
  theta_.print("Values:\n");
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
bool ConcurrentBatchSmoother::equals(const ConcurrentSmoother& rhs, double tol) const {
 | 
						|
  const ConcurrentBatchSmoother* smoother = dynamic_cast<const ConcurrentBatchSmoother*>(&rhs);
 | 
						|
  return smoother
 | 
						|
      && factors_.equals(smoother->factors_)
 | 
						|
      && theta_.equals(smoother->theta_)
 | 
						|
      && ordering_.equals(smoother->ordering_)
 | 
						|
      && delta_.equals(smoother->delta_)
 | 
						|
      && separatorValues_.equals(smoother->separatorValues_);
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
ConcurrentBatchSmoother::Result ConcurrentBatchSmoother::update(const NonlinearFactorGraph& newFactors, const Values& newTheta,
 | 
						|
    const boost::optional< std::vector<size_t> >& removeFactorIndices) {
 | 
						|
 | 
						|
  gttic(update);
 | 
						|
 | 
						|
  // Create the return result meta-data
 | 
						|
  Result result;
 | 
						|
 | 
						|
  // Update all of the internal variables with the new information
 | 
						|
  gttic(augment_system);
 | 
						|
  {
 | 
						|
    // Add the new variables to theta
 | 
						|
    theta_.insert(newTheta);
 | 
						|
 | 
						|
    // Add new variables to the end of the ordering
 | 
						|
    BOOST_FOREACH(const Values::ConstKeyValuePair& key_value, newTheta) {
 | 
						|
      ordering_.push_back(key_value.key);
 | 
						|
    }
 | 
						|
 | 
						|
    // Augment Delta
 | 
						|
    delta_.insert(newTheta.zeroVectors());
 | 
						|
 | 
						|
    // Add the new factors to the graph, updating the variable index
 | 
						|
    insertFactors(newFactors);
 | 
						|
 | 
						|
    if(removeFactorIndices)
 | 
						|
      removeFactors(*removeFactorIndices);
 | 
						|
  }
 | 
						|
  gttoc(augment_system);
 | 
						|
 | 
						|
  if(factors_.size() > 0) {
 | 
						|
    // Reorder the system to ensure efficient optimization (and marginalization) performance
 | 
						|
    gttic(reorder);
 | 
						|
    reorder();
 | 
						|
    gttoc(reorder);
 | 
						|
 | 
						|
    // Optimize the factors using a modified version of L-M
 | 
						|
    gttic(optimize);
 | 
						|
    result = optimize();
 | 
						|
    gttoc(optimize);
 | 
						|
  }
 | 
						|
 | 
						|
  // TODO: The following code does considerable work, much of which could be redundant given the previous optimization step
 | 
						|
  // Refactor this code to reduce computational burden
 | 
						|
 | 
						|
  // Calculate the marginal on the separator from the smoother factors
 | 
						|
  if(separatorValues_.size() > 0) {
 | 
						|
    gttic(presync);
 | 
						|
    updateSmootherSummarization();
 | 
						|
    gttoc(presync);
 | 
						|
  }
 | 
						|
 | 
						|
  gttoc(update);
 | 
						|
 | 
						|
  return result;
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
void ConcurrentBatchSmoother::presync() {
 | 
						|
 | 
						|
  gttic(presync);
 | 
						|
 | 
						|
  gttoc(presync);
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
void ConcurrentBatchSmoother::getSummarizedFactors(NonlinearFactorGraph& summarizedFactors, Values& separatorValues) {
 | 
						|
 | 
						|
  gttic(get_summarized_factors);
 | 
						|
 | 
						|
  // Copy the previous calculated smoother summarization factors into the output
 | 
						|
  summarizedFactors.push_back(smootherSummarization_);
 | 
						|
 | 
						|
  // Copy the separator values into the output
 | 
						|
  separatorValues.insert(separatorValues_);
 | 
						|
 | 
						|
  gttoc(get_summarized_factors);
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
void ConcurrentBatchSmoother::synchronize(const NonlinearFactorGraph& smootherFactors, const Values& smootherValues,
 | 
						|
    const NonlinearFactorGraph& summarizedFactors, const Values& separatorValues) {
 | 
						|
 | 
						|
  gttic(synchronize);
 | 
						|
 | 
						|
  // Remove the previous filter summarization from the graph
 | 
						|
  removeFactors(filterSummarizationSlots_);
 | 
						|
 | 
						|
  // Insert new linpoints into the values, augment the ordering, and store new dims to augment delta
 | 
						|
  BOOST_FOREACH(const Values::ConstKeyValuePair& key_value, smootherValues) {
 | 
						|
    std::pair<Values::iterator, bool> iter_inserted = theta_.tryInsert(key_value.key, key_value.value);
 | 
						|
    if(iter_inserted.second) {
 | 
						|
      // If the insert succeeded
 | 
						|
      ordering_.push_back(key_value.key);
 | 
						|
      delta_.insert(key_value.key, Vector::Zero(key_value.value.dim()));
 | 
						|
    } else {
 | 
						|
      // If the element already existed in theta_
 | 
						|
      iter_inserted.first->value = key_value.value;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  BOOST_FOREACH(const Values::ConstKeyValuePair& key_value, separatorValues) {
 | 
						|
    std::pair<Values::iterator, bool> iter_inserted = theta_.tryInsert(key_value.key, key_value.value);
 | 
						|
    if(iter_inserted.second) {
 | 
						|
      // If the insert succeeded
 | 
						|
      ordering_.push_back(key_value.key);
 | 
						|
      delta_.insert(key_value.key, Vector::Zero(key_value.value.dim()));
 | 
						|
    } else {
 | 
						|
      // If the element already existed in theta_
 | 
						|
      iter_inserted.first->value = key_value.value;
 | 
						|
    }
 | 
						|
  }
 | 
						|
 | 
						|
  // Insert the new smoother factors
 | 
						|
  insertFactors(smootherFactors);
 | 
						|
 | 
						|
  // Insert the new filter summarized factors
 | 
						|
  filterSummarizationSlots_ = insertFactors(summarizedFactors);
 | 
						|
 | 
						|
  // Update the list of root keys
 | 
						|
  separatorValues_ = separatorValues;
 | 
						|
 | 
						|
  gttoc(synchronize);
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
void ConcurrentBatchSmoother::postsync() {
 | 
						|
 | 
						|
  gttic(postsync);
 | 
						|
 | 
						|
  gttoc(postsync);
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
std::vector<size_t> ConcurrentBatchSmoother::insertFactors(const NonlinearFactorGraph& factors) {
 | 
						|
 | 
						|
  gttic(insert_factors);
 | 
						|
 | 
						|
  // create the output vector
 | 
						|
  std::vector<size_t> slots;
 | 
						|
  slots.reserve(factors.size());
 | 
						|
 | 
						|
  // Insert the factor into an existing hole in the factor graph, if possible
 | 
						|
  BOOST_FOREACH(const NonlinearFactor::shared_ptr& factor, factors) {
 | 
						|
    size_t slot;
 | 
						|
    if(availableSlots_.size() > 0) {
 | 
						|
      slot = availableSlots_.front();
 | 
						|
      availableSlots_.pop();
 | 
						|
      factors_.replace(slot, factor);
 | 
						|
    } else {
 | 
						|
      slot = factors_.size();
 | 
						|
      factors_.push_back(factor);
 | 
						|
    }
 | 
						|
    slots.push_back(slot);
 | 
						|
  }
 | 
						|
 | 
						|
  gttoc(insert_factors);
 | 
						|
 | 
						|
  return slots;
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
void ConcurrentBatchSmoother::removeFactors(const std::vector<size_t>& slots) {
 | 
						|
 | 
						|
  gttic(remove_factors);
 | 
						|
 | 
						|
  // For each factor slot to delete...
 | 
						|
  BOOST_FOREACH(size_t slot, slots) {
 | 
						|
 | 
						|
    // Remove the factor from the graph
 | 
						|
    factors_.remove(slot);
 | 
						|
 | 
						|
    // Mark the factor slot as available
 | 
						|
    availableSlots_.push(slot);
 | 
						|
  }
 | 
						|
 | 
						|
  gttoc(remove_factors);
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
void ConcurrentBatchSmoother::reorder() {
 | 
						|
 | 
						|
  // Recalculate the variable index
 | 
						|
  variableIndex_ = VariableIndex(factors_);
 | 
						|
 | 
						|
  FastList<Key> separatorKeys = separatorValues_.keys();
 | 
						|
  ordering_ = Ordering::COLAMDConstrainedLast(variableIndex_, std::vector<Key>(separatorKeys.begin(), separatorKeys.end()));
 | 
						|
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
ConcurrentBatchSmoother::Result ConcurrentBatchSmoother::optimize() {
 | 
						|
 | 
						|
  // Create output result structure
 | 
						|
  Result result;
 | 
						|
  result.nonlinearVariables = theta_.size() - separatorValues_.size();
 | 
						|
  result.linearVariables = separatorValues_.size();
 | 
						|
 | 
						|
  // Pull out parameters we'll use
 | 
						|
  const LevenbergMarquardtParams::VerbosityLM lmVerbosity = parameters_.verbosityLM;
 | 
						|
  double lambda = parameters_.lambdaInitial;
 | 
						|
 | 
						|
  // Create a Values that holds the current evaluation point
 | 
						|
  Values evalpoint = theta_.retract(delta_);
 | 
						|
  result.error = factors_.error(evalpoint);
 | 
						|
  if(result.error < parameters_.errorTol) {
 | 
						|
    return result;
 | 
						|
  }
 | 
						|
 | 
						|
  // Use a custom optimization loop so the linearization points can be controlled
 | 
						|
  double previousError;
 | 
						|
  VectorValues newDelta;
 | 
						|
  do {
 | 
						|
    previousError = result.error;
 | 
						|
 | 
						|
    // Do next iteration
 | 
						|
    gttic(optimizer_iteration);
 | 
						|
    {
 | 
						|
      // Linearize graph around the linearization point
 | 
						|
      GaussianFactorGraph linearFactorGraph = *factors_.linearize(theta_);
 | 
						|
 | 
						|
      // Keep increasing lambda until we make make progress
 | 
						|
      while(true) {
 | 
						|
        if (lmVerbosity >= LevenbergMarquardtParams::TRYLAMBDA)
 | 
						|
          std::cout << "trying lambda = " << lambda << std::endl;
 | 
						|
        
 | 
						|
        // Add prior factors at the current solution
 | 
						|
        gttic(damp);
 | 
						|
        GaussianFactorGraph dampedFactorGraph(linearFactorGraph);
 | 
						|
        dampedFactorGraph.reserve(linearFactorGraph.size() + delta_.size());
 | 
						|
        {
 | 
						|
          // for each of the variables, add a prior at the current solution
 | 
						|
          BOOST_FOREACH(const VectorValues::KeyValuePair& key_value, delta_) {
 | 
						|
            size_t dim = key_value.second.size();
 | 
						|
            Matrix A = Matrix::Identity(dim,dim);
 | 
						|
            Vector b = key_value.second;
 | 
						|
            SharedDiagonal model = noiseModel::Isotropic::Sigma(dim, 1.0 / std::sqrt(lambda));
 | 
						|
            GaussianFactor::shared_ptr prior(new JacobianFactor(key_value.first, A, b, model));
 | 
						|
            dampedFactorGraph.push_back(prior);
 | 
						|
          }
 | 
						|
        }
 | 
						|
        gttoc(damp);
 | 
						|
        if (lmVerbosity >= LevenbergMarquardtParams::DAMPED) 
 | 
						|
          dampedFactorGraph.print("damped");
 | 
						|
        result.lambdas++;
 | 
						|
 | 
						|
        gttic(solve);
 | 
						|
        // Solve Damped Gaussian Factor Graph
 | 
						|
        newDelta = dampedFactorGraph.optimize(ordering_, parameters_.getEliminationFunction());
 | 
						|
        // update the evalpoint with the new delta
 | 
						|
        evalpoint = theta_.retract(newDelta);
 | 
						|
        gttoc(solve);
 | 
						|
 | 
						|
        if (lmVerbosity >= LevenbergMarquardtParams::TRYLAMBDA) 
 | 
						|
          std::cout << "linear delta norm = " << newDelta.norm() << std::endl;
 | 
						|
        if (lmVerbosity >= LevenbergMarquardtParams::TRYDELTA) 
 | 
						|
          newDelta.print("delta");
 | 
						|
 | 
						|
        // Evaluate the new error
 | 
						|
        gttic(compute_error);
 | 
						|
        double error = factors_.error(evalpoint);
 | 
						|
        gttoc(compute_error);
 | 
						|
 | 
						|
        if (lmVerbosity >= LevenbergMarquardtParams::TRYLAMBDA) 
 | 
						|
          std::cout << "next error = " << error << std::endl;
 | 
						|
        
 | 
						|
        if(error < result.error) {
 | 
						|
          // Keep this change
 | 
						|
          // Update the error value
 | 
						|
          result.error = error;
 | 
						|
          // Update the linearization point
 | 
						|
          theta_ = evalpoint;
 | 
						|
          // Reset the deltas to zeros
 | 
						|
          delta_.setZero();
 | 
						|
          // Put the linearization points and deltas back for specific variables
 | 
						|
          if(separatorValues_.size() > 0) {
 | 
						|
            theta_.update(separatorValues_);
 | 
						|
            BOOST_FOREACH(const Values::ConstKeyValuePair& key_value, separatorValues_) {
 | 
						|
              delta_.at(key_value.key) = newDelta.at(key_value.key);
 | 
						|
            }
 | 
						|
          }
 | 
						|
 | 
						|
          // Decrease lambda for next time
 | 
						|
          lambda /= parameters_.lambdaFactor;
 | 
						|
          // End this lambda search iteration
 | 
						|
          break;
 | 
						|
        } else {
 | 
						|
          // Reject this change
 | 
						|
          if(lambda >= parameters_.lambdaUpperBound) {
 | 
						|
            // The maximum lambda has been used. Print a warning and end the search.
 | 
						|
            std::cout << "Warning:  Levenberg-Marquardt giving up because cannot decrease error with maximum lambda" << std::endl;
 | 
						|
            break;
 | 
						|
          } else {
 | 
						|
            // Increase lambda and continue searching
 | 
						|
            lambda *= parameters_.lambdaFactor;
 | 
						|
          }
 | 
						|
        }
 | 
						|
      } // end while
 | 
						|
    }
 | 
						|
    gttoc(optimizer_iteration);
 | 
						|
 | 
						|
    if (lmVerbosity >= LevenbergMarquardtParams::LAMBDA)
 | 
						|
      std::cout << "using lambda = " << lambda << std::endl;
 | 
						|
 | 
						|
    result.iterations++;
 | 
						|
  } while(result.iterations < (size_t)parameters_.maxIterations &&
 | 
						|
      !checkConvergence(parameters_.relativeErrorTol, parameters_.absoluteErrorTol, parameters_.errorTol, previousError, result.error, NonlinearOptimizerParams::SILENT));
 | 
						|
 | 
						|
  return result;
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
void ConcurrentBatchSmoother::updateSmootherSummarization() {
 | 
						|
 | 
						|
  // The smoother summarization factors are the resulting marginal factors on the separator
 | 
						|
  // variables that result from marginalizing out all of the other variables
 | 
						|
  // These marginal factors will be cached for later transmission to the filter using
 | 
						|
  // linear container factors
 | 
						|
 | 
						|
  // Create a nonlinear factor graph without the filter summarization factors
 | 
						|
  NonlinearFactorGraph graph(factors_);
 | 
						|
  BOOST_FOREACH(size_t slot, filterSummarizationSlots_) {
 | 
						|
    graph.remove(slot);
 | 
						|
  }
 | 
						|
 | 
						|
  // Get the set of separator keys
 | 
						|
  gtsam::FastSet<Key> separatorKeys;
 | 
						|
  BOOST_FOREACH(const Values::ConstKeyValuePair& key_value, separatorValues_) {
 | 
						|
    separatorKeys.insert(key_value.key);
 | 
						|
  }
 | 
						|
 | 
						|
  // Calculate the marginal factors on the separator
 | 
						|
  smootherSummarization_ = internal::calculateMarginalFactors(graph, theta_, separatorKeys, parameters_.getEliminationFunction());
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
void ConcurrentBatchSmoother::PrintNonlinearFactor(const NonlinearFactor::shared_ptr& factor, const std::string& indent, const KeyFormatter& keyFormatter) {
 | 
						|
  std::cout << indent;
 | 
						|
  if(factor) {
 | 
						|
    if(boost::dynamic_pointer_cast<LinearContainerFactor>(factor)) {
 | 
						|
      std::cout << "l( ";
 | 
						|
    } else {
 | 
						|
      std::cout << "f( ";
 | 
						|
    }
 | 
						|
    BOOST_FOREACH(Key key, *factor) {
 | 
						|
      std::cout << keyFormatter(key) << " ";
 | 
						|
    }
 | 
						|
    std::cout << ")" << std::endl;
 | 
						|
  } else {
 | 
						|
    std::cout << "{ NULL }" << std::endl;
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
void ConcurrentBatchSmoother::PrintLinearFactor(const GaussianFactor::shared_ptr& factor, const std::string& indent, const KeyFormatter& keyFormatter) {
 | 
						|
  std::cout << indent;
 | 
						|
  if(factor) {
 | 
						|
    std::cout << "g( ";
 | 
						|
    BOOST_FOREACH(Key key, *factor) {
 | 
						|
      std::cout << keyFormatter(key) << " ";
 | 
						|
    }
 | 
						|
    std::cout << ")" << std::endl;
 | 
						|
  } else {
 | 
						|
    std::cout << "{ NULL }" << std::endl;
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
/* ************************************************************************* */
 | 
						|
}/// namespace gtsam
 |