parameter_block.h 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388
  1. // Ceres Solver - A fast non-linear least squares minimizer
  2. // Copyright 2023 Google Inc. All rights reserved.
  3. // http://ceres-solver.org/
  4. //
  5. // Redistribution and use in source and binary forms, with or without
  6. // modification, are permitted provided that the following conditions are met:
  7. //
  8. // * Redistributions of source code must retain the above copyright notice,
  9. // this list of conditions and the following disclaimer.
  10. // * Redistributions in binary form must reproduce the above copyright notice,
  11. // this list of conditions and the following disclaimer in the documentation
  12. // and/or other materials provided with the distribution.
  13. // * Neither the name of Google Inc. nor the names of its contributors may be
  14. // used to endorse or promote products derived from this software without
  15. // specific prior written permission.
  16. //
  17. // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18. // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19. // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20. // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  21. // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22. // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23. // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24. // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25. // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26. // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27. // POSSIBILITY OF SUCH DAMAGE.
  28. //
  29. // Author: keir@google.com (Keir Mierle)
  30. #ifndef CERES_INTERNAL_PARAMETER_BLOCK_H_
  31. #define CERES_INTERNAL_PARAMETER_BLOCK_H_
  32. #include <algorithm>
  33. #include <cstdint>
  34. #include <cstdlib>
  35. #include <limits>
  36. #include <memory>
  37. #include <string>
  38. #include <unordered_set>
  39. #include "ceres/array_utils.h"
  40. #include "ceres/internal/disable_warnings.h"
  41. #include "ceres/internal/eigen.h"
  42. #include "ceres/internal/export.h"
  43. #include "ceres/manifold.h"
  44. #include "ceres/stringprintf.h"
  45. #include "glog/logging.h"
  46. namespace ceres::internal {
  47. class ProblemImpl;
  48. class ResidualBlock;
  49. // The parameter block encodes the location of the user's original value, and
  50. // also the "current state" of the parameter. The evaluator uses whatever is in
  51. // the current state of the parameter when evaluating. This is inlined since the
  52. // methods are performance sensitive.
  53. //
  54. // The class is not thread-safe, unless only const methods are called. The
  55. // parameter block may also hold a pointer to a manifold; the parameter block
  56. // does not take ownership of this pointer, so the user is responsible for the
  57. // proper disposal of the manifold.
  58. class CERES_NO_EXPORT ParameterBlock {
  59. public:
  60. using ResidualBlockSet = std::unordered_set<ResidualBlock*>;
  61. // Create a parameter block with the user state, size, and index specified.
  62. // The size is the size of the parameter block and the index is the position
  63. // of the parameter block inside a Program (if any).
  64. ParameterBlock(double* user_state, int size, int index)
  65. : user_state_(user_state),
  66. size_(size),
  67. state_(user_state),
  68. index_(index) {}
  69. ParameterBlock(double* user_state, int size, int index, Manifold* manifold)
  70. : user_state_(user_state),
  71. size_(size),
  72. state_(user_state),
  73. index_(index) {
  74. if (manifold != nullptr) {
  75. SetManifold(manifold);
  76. }
  77. }
  78. // The size of the parameter block.
  79. int Size() const { return size_; }
  80. // Manipulate the parameter state.
  81. bool SetState(const double* x) {
  82. CHECK(x != nullptr) << "Tried to set the state of constant parameter "
  83. << "with user location " << user_state_;
  84. CHECK(!IsConstant()) << "Tried to set the state of constant parameter "
  85. << "with user location " << user_state_;
  86. state_ = x;
  87. return UpdatePlusJacobian();
  88. }
  89. // Copy the current parameter state out to x. This is "GetState()" rather than
  90. // simply "state()" since it is actively copying the data into the passed
  91. // pointer.
  92. void GetState(double* x) const {
  93. if (x != state_) {
  94. std::copy(state_, state_ + size_, x);
  95. }
  96. }
  97. // Direct pointers to the current state.
  98. const double* state() const { return state_; }
  99. const double* user_state() const { return user_state_; }
  100. double* mutable_user_state() { return user_state_; }
  101. const Manifold* manifold() const { return manifold_; }
  102. Manifold* mutable_manifold() { return manifold_; }
  103. // Set this parameter block to vary or not.
  104. void SetConstant() { is_set_constant_ = true; }
  105. void SetVarying() { is_set_constant_ = false; }
  106. bool IsConstant() const { return (is_set_constant_ || TangentSize() == 0); }
  107. double UpperBound(int index) const {
  108. return (upper_bounds_ ? upper_bounds_[index]
  109. : std::numeric_limits<double>::max());
  110. }
  111. double LowerBound(int index) const {
  112. return (lower_bounds_ ? lower_bounds_[index]
  113. : -std::numeric_limits<double>::max());
  114. }
  115. bool IsUpperBounded() const { return (upper_bounds_ == nullptr); }
  116. bool IsLowerBounded() const { return (lower_bounds_ == nullptr); }
  117. // This parameter block's index in an array.
  118. int index() const { return index_; }
  119. void set_index(int index) { index_ = index; }
  120. // This parameter offset inside a larger state vector.
  121. int state_offset() const { return state_offset_; }
  122. void set_state_offset(int state_offset) { state_offset_ = state_offset; }
  123. // This parameter offset inside a larger delta vector.
  124. int delta_offset() const { return delta_offset_; }
  125. void set_delta_offset(int delta_offset) { delta_offset_ = delta_offset; }
  126. // Methods relating to the parameter block's manifold.
  127. // The local to global jacobian. Returns nullptr if there is no manifold for
  128. // this parameter block. The returned matrix is row-major and has Size() rows
  129. // and TangentSize() columns.
  130. const double* PlusJacobian() const { return plus_jacobian_.get(); }
  131. int TangentSize() const {
  132. return (manifold_ == nullptr) ? size_ : manifold_->TangentSize();
  133. }
  134. // Set the manifold. The parameter block does not take ownership of
  135. // the manifold.
  136. void SetManifold(Manifold* new_manifold) {
  137. // Nothing to do if the new manifold is the same as the old
  138. // manifold.
  139. if (new_manifold == manifold_) {
  140. return;
  141. }
  142. if (new_manifold == nullptr) {
  143. manifold_ = nullptr;
  144. plus_jacobian_ = nullptr;
  145. return;
  146. }
  147. CHECK_EQ(new_manifold->AmbientSize(), size_)
  148. << "The parameter block has size = " << size_
  149. << " while the manifold has ambient size = "
  150. << new_manifold->AmbientSize();
  151. CHECK_GE(new_manifold->TangentSize(), 0)
  152. << "Invalid Manifold. Manifolds must have a "
  153. << "non-negative dimensional tangent space.";
  154. manifold_ = new_manifold;
  155. plus_jacobian_ = std::make_unique<double[]>(manifold_->AmbientSize() *
  156. manifold_->TangentSize());
  157. CHECK(UpdatePlusJacobian())
  158. << "Manifold::PlusJacobian computation failed for x: "
  159. << ConstVectorRef(state_, Size()).transpose();
  160. }
  161. void SetUpperBound(int index, double upper_bound) {
  162. CHECK_LT(index, size_);
  163. if (upper_bound >= std::numeric_limits<double>::max() && !upper_bounds_) {
  164. return;
  165. }
  166. if (!upper_bounds_) {
  167. upper_bounds_ = std::make_unique<double[]>(size_);
  168. std::fill(upper_bounds_.get(),
  169. upper_bounds_.get() + size_,
  170. std::numeric_limits<double>::max());
  171. }
  172. upper_bounds_[index] = upper_bound;
  173. }
  174. void SetLowerBound(int index, double lower_bound) {
  175. CHECK_LT(index, size_);
  176. if (lower_bound <= -std::numeric_limits<double>::max() && !lower_bounds_) {
  177. return;
  178. }
  179. if (!lower_bounds_) {
  180. lower_bounds_ = std::make_unique<double[]>(size_);
  181. std::fill(lower_bounds_.get(),
  182. lower_bounds_.get() + size_,
  183. -std::numeric_limits<double>::max());
  184. }
  185. lower_bounds_[index] = lower_bound;
  186. }
  187. // Generalization of the addition operation. This is the same as
  188. // Manifold::Plus() followed by projection onto the
  189. // hyper cube implied by the bounds constraints.
  190. bool Plus(const double* x, const double* delta, double* x_plus_delta) {
  191. if (manifold_ != nullptr) {
  192. if (!manifold_->Plus(x, delta, x_plus_delta)) {
  193. return false;
  194. }
  195. } else {
  196. VectorRef(x_plus_delta, size_) =
  197. ConstVectorRef(x, size_) + ConstVectorRef(delta, size_);
  198. }
  199. // Project onto the box constraints.
  200. if (lower_bounds_.get() != nullptr) {
  201. for (int i = 0; i < size_; ++i) {
  202. x_plus_delta[i] = std::max(x_plus_delta[i], lower_bounds_[i]);
  203. }
  204. }
  205. if (upper_bounds_.get() != nullptr) {
  206. for (int i = 0; i < size_; ++i) {
  207. x_plus_delta[i] = std::min(x_plus_delta[i], upper_bounds_[i]);
  208. }
  209. }
  210. return true;
  211. }
  212. std::string ToString() const {
  213. return StringPrintf(
  214. "{ this=%p, user_state=%p, state=%p, size=%d, "
  215. "constant=%d, index=%d, state_offset=%d, "
  216. "delta_offset=%d }",
  217. this,
  218. user_state_,
  219. state_,
  220. size_,
  221. is_set_constant_,
  222. index_,
  223. state_offset_,
  224. delta_offset_);
  225. }
  226. void EnableResidualBlockDependencies() {
  227. CHECK(residual_blocks_.get() == nullptr)
  228. << "Ceres bug: There is already a residual block collection "
  229. << "for parameter block: " << ToString();
  230. residual_blocks_ = std::make_unique<ResidualBlockSet>();
  231. }
  232. void AddResidualBlock(ResidualBlock* residual_block) {
  233. CHECK(residual_blocks_.get() != nullptr)
  234. << "Ceres bug: The residual block collection is null for parameter "
  235. << "block: " << ToString();
  236. residual_blocks_->insert(residual_block);
  237. }
  238. void RemoveResidualBlock(ResidualBlock* residual_block) {
  239. CHECK(residual_blocks_.get() != nullptr)
  240. << "Ceres bug: The residual block collection is null for parameter "
  241. << "block: " << ToString();
  242. CHECK(residual_blocks_->find(residual_block) != residual_blocks_->end())
  243. << "Ceres bug: Missing residual for parameter block: " << ToString();
  244. residual_blocks_->erase(residual_block);
  245. }
  246. // This is only intended for iterating; perhaps this should only expose
  247. // .begin() and .end().
  248. ResidualBlockSet* mutable_residual_blocks() { return residual_blocks_.get(); }
  249. double LowerBoundForParameter(int index) const {
  250. if (lower_bounds_.get() == nullptr) {
  251. return -std::numeric_limits<double>::max();
  252. } else {
  253. return lower_bounds_[index];
  254. }
  255. }
  256. double UpperBoundForParameter(int index) const {
  257. if (upper_bounds_.get() == nullptr) {
  258. return std::numeric_limits<double>::max();
  259. } else {
  260. return upper_bounds_[index];
  261. }
  262. }
  263. private:
  264. bool UpdatePlusJacobian() {
  265. if (manifold_ == nullptr) {
  266. return true;
  267. }
  268. // Update the Plus Jacobian. In some cases this is
  269. // wasted effort; if this is a bottleneck, we will find a solution
  270. // at that time.
  271. const int jacobian_size = Size() * TangentSize();
  272. InvalidateArray(jacobian_size, plus_jacobian_.get());
  273. if (!manifold_->PlusJacobian(state_, plus_jacobian_.get())) {
  274. LOG(WARNING) << "Manifold::PlusJacobian computation failed"
  275. "for x: "
  276. << ConstVectorRef(state_, Size()).transpose();
  277. return false;
  278. }
  279. if (!IsArrayValid(jacobian_size, plus_jacobian_.get())) {
  280. LOG(WARNING) << "Manifold::PlusJacobian computation returned "
  281. << "an invalid matrix for x: "
  282. << ConstVectorRef(state_, Size()).transpose()
  283. << "\n Jacobian matrix : "
  284. << ConstMatrixRef(
  285. plus_jacobian_.get(), Size(), TangentSize());
  286. return false;
  287. }
  288. return true;
  289. }
  290. double* user_state_ = nullptr;
  291. int size_ = -1;
  292. bool is_set_constant_ = false;
  293. Manifold* manifold_ = nullptr;
  294. // The "state" of the parameter. These fields are only needed while the
  295. // solver is running. While at first glance using mutable is a bad idea, this
  296. // ends up simplifying the internals of Ceres enough to justify the potential
  297. // pitfalls of using "mutable."
  298. mutable const double* state_ = nullptr;
  299. mutable std::unique_ptr<double[]> plus_jacobian_;
  300. // The index of the parameter. This is used by various other parts of Ceres to
  301. // permit switching from a ParameterBlock* to an index in another array.
  302. int index_ = -1;
  303. // The offset of this parameter block inside a larger state vector.
  304. int state_offset_ = -1;
  305. // The offset of this parameter block inside a larger delta vector.
  306. int delta_offset_ = -1;
  307. // If non-null, contains the residual blocks this parameter block is in.
  308. std::unique_ptr<ResidualBlockSet> residual_blocks_;
  309. // Upper and lower bounds for the parameter block. SetUpperBound
  310. // and SetLowerBound lazily initialize the upper_bounds_ and
  311. // lower_bounds_ arrays. If they are never called, then memory for
  312. // these arrays is never allocated. Thus for problems where there
  313. // are no bounds, or only one sided bounds we do not pay the cost of
  314. // allocating memory for the inactive bounds constraints.
  315. //
  316. // Upon initialization these arrays are initialized to
  317. // std::numeric_limits<double>::max() and
  318. // -std::numeric_limits<double>::max() respectively which correspond
  319. // to the parameter block being unconstrained.
  320. std::unique_ptr<double[]> upper_bounds_;
  321. std::unique_ptr<double[]> lower_bounds_;
  322. // Necessary so ProblemImpl can clean up the manifolds.
  323. friend class ProblemImpl;
  324. };
  325. } // namespace ceres::internal
  326. #include "ceres/internal/reenable_warnings.h"
  327. #endif // CERES_INTERNAL_PARAMETER_BLOCK_H_