// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2023 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
//   this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
//   this list of conditions and the following disclaimer in the documentation
//   and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
//   used to endorse or promote products derived from this software without
//   specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
//         keir@google.com (Keir Mierle)
//
// The Problem object is used to build and hold least squares problems.

#ifndef CERES_PUBLIC_PROBLEM_H_
#define CERES_PUBLIC_PROBLEM_H_

#include <array>
#include <cstddef>
#include <map>
#include <memory>
#include <set>
#include <vector>

#include "ceres/context.h"
#include "ceres/internal/disable_warnings.h"
#include "ceres/internal/export.h"
#include "ceres/internal/port.h"
#include "ceres/types.h"

namespace ceres {

class CostFunction;
class EvaluationCallback;
class LossFunction;
class Manifold;
class Solver;
struct CRSMatrix;

namespace internal {
class Preprocessor;
class ProblemImpl;
class ParameterBlock;
class ResidualBlock;
}  // namespace internal

// A ResidualBlockId is an opaque handle clients can use to remove residual
// blocks from a Problem after adding them.
using ResidualBlockId = internal::ResidualBlock*;

// A class to represent non-linear least squares problems. Such
// problems have a cost function that is a sum of error terms (known
// as "residuals"), where each residual is a function of some subset
// of the parameters. The cost function takes the form
//
//    N    1
//   SUM  --- loss( || r_i1, r_i2,..., r_ik ||^2  ),
//   i=1   2
//
// where
//
//   r_ij is residual number i, component j; the residual is a function of some
//        subset of the parameters x1...xk. For example, in a structure from
//        motion problem a residual might be the difference between a measured
//        point in an image and the reprojected position for the matching
//        camera, point pair. The residual would have two components, error in x
//        and error in y.
//
//   loss(y) is the loss function; for example, squared error or Huber L1
//           loss. If loss(y) = y, then the cost function is non-robustified
//           least squares.
//
// This class is specifically designed to address the important subset of
// "sparse" least squares problems, where each component of the residual depends
// only on a small number of parameters, even though the total number of
// residuals and parameters may be very large. This property affords tremendous
// gains in scale, allowing efficient solving of large problems that are
// otherwise inaccessible.
//
// The canonical example of a sparse least squares problem is
// "structure-from-motion" (SFM), where the parameters are points and cameras,
// and residuals are reprojection errors. Typically a single residual will
// depend only on 9 parameters (3 for the point, 6 for the camera).
//
// To create a least squares problem, use the AddResidualBlock() and
// AddParameterBlock() methods, documented below. Here is an example least
// squares problem containing 3 parameter blocks of sizes 3, 4 and 5
// respectively and two residual terms of size 2 and 6:
//
//   double x1[] = { 1.0, 2.0, 3.0 };
//   double x2[] = { 1.0, 2.0, 3.0, 5.0 };
//   double x3[] = { 1.0, 2.0, 3.0, 6.0, 7.0 };
//
//   Problem problem;
//
//   problem.AddResidualBlock(new MyUnaryCostFunction(...), nullptr, x1);
//   problem.AddResidualBlock(new MyBinaryCostFunction(...), nullptr, x2, x3);
//
// Please see cost_function.h for details of the CostFunction object.
class CERES_EXPORT Problem {
 public:
  struct CERES_EXPORT Options {
    // These flags control whether the Problem object owns the CostFunctions,
    // LossFunctions, and Manifolds passed into the Problem.
    //
    // If set to TAKE_OWNERSHIP, then the problem object will delete the
    // corresponding object on destruction. The destructor is careful to delete
    // the pointers only once, since sharing objects is allowed.
    Ownership cost_function_ownership = TAKE_OWNERSHIP;
    Ownership loss_function_ownership = TAKE_OWNERSHIP;
    Ownership manifold_ownership = TAKE_OWNERSHIP;

    // If true, trades memory for faster RemoveResidualBlock() and
    // RemoveParameterBlock() operations.
    //
    // By default, RemoveParameterBlock() and RemoveResidualBlock() take time
    // proportional to the size of the entire problem. If you only ever remove
    // parameters or residuals from the problem occasionally, this might be
    // acceptable. However, if you have memory to spare, enable this option to
    // make RemoveParameterBlock() take time proportional to the number of
    // residual blocks that depend on it, and RemoveResidualBlock() take (on
    // average) constant time.
    //
    // The increase in memory usage is two-fold: an additional hash set per
    // parameter block containing all the residuals that depend on the parameter
    // block; and a hash set in the problem containing all residuals.
    bool enable_fast_removal = false;

    // By default, Ceres performs a variety of safety checks when constructing
    // the problem. There is a small but measurable performance penalty to these
    // checks, typically around 5% of construction time. If you are sure your
    // problem construction is correct, and 5% of the problem construction time
    // is truly overhead you want to avoid, then you can set
    // disable_all_safety_checks to true.
    //
    // WARNING: Do not set this to true, unless you are absolutely sure of what
    // you are doing.
    bool disable_all_safety_checks = false;

    // A Ceres global context to use for solving this problem. This may help to
    // reduce computation time as Ceres can reuse expensive objects to create.
    // The context object can be nullptr, in which case Ceres may create one.
    //
    // Ceres does NOT take ownership of the pointer.
    Context* context = nullptr;

    // Using this callback interface, Ceres can notify you when it is about to
    // evaluate the residuals or Jacobians. With the callback, you can share
    // computation between residual blocks by doing the shared computation in
    // EvaluationCallback::PrepareForEvaluation() before Ceres calls
    // CostFunction::Evaluate(). It also enables caching results between a pure
    // residual evaluation and a residual & Jacobian evaluation.
    //
    // Problem DOES NOT take ownership of the callback.
    //
    // NOTE: Evaluation callbacks are incompatible with inner iterations. So
    // calling Solve with Solver::Options::use_inner_iterations = true on a
    // Problem with a non-null evaluation callback is an error.
    EvaluationCallback* evaluation_callback = nullptr;
  };

  // The default constructor is equivalent to the invocation
  // Problem(Problem::Options()).
  Problem();
  explicit Problem(const Options& options);

  // Move constructor and assignment.
  Problem(Problem&&);
  Problem& operator=(Problem&&);

  Problem(const Problem&) = delete;
  Problem& operator=(const Problem&) = delete;

  ~Problem();

  // Add a residual block to the overall cost function. The cost function
  // carries with its information about the sizes of the parameter blocks it
  // expects. The function checks that these match the sizes of the parameter
  // blocks listed in parameter_blocks. The program aborts if a mismatch is
  // detected. loss_function can be nullptr, in which case the cost of the term
  // is just the squared norm of the residuals.
  //
  // The user has the option of explicitly adding the parameter blocks using
  // AddParameterBlock. This causes additional correctness checking; however,
  // AddResidualBlock implicitly adds the parameter blocks if they are not
  // present, so calling AddParameterBlock explicitly is not required.
  //
  // The Problem object by default takes ownership of the cost_function and
  // loss_function pointers (See Problem::Options to override this behaviour).
  // These objects remain live for the life of the Problem object. If the user
  // wishes to keep control over the destruction of these objects, then they can
  // do this by setting the corresponding enums in the Options struct.
  //
  // Note: Even though the Problem takes ownership of cost_function and
  // loss_function, it does not preclude the user from re-using them in another
  // residual block. The destructor takes care to call delete on each
  // cost_function or loss_function pointer only once, regardless of how many
  // residual blocks refer to them.
  //
  // Example usage:
  //
  //   double x1[] = {1.0, 2.0, 3.0};
  //   double x2[] = {1.0, 2.0, 5.0, 6.0};
  //   double x3[] = {3.0, 6.0, 2.0, 5.0, 1.0};
  //
  //   Problem problem;
  //
  //   problem.AddResidualBlock(new MyUnaryCostFunction(...), nullptr, x1);
  //   problem.AddResidualBlock(new MyBinaryCostFunction(...), nullptr, x2, x1);
  //
  // Add a residual block by listing the parameter block pointers directly
  // instead of wrapping them in a container.
  template <typename... Ts>
  ResidualBlockId AddResidualBlock(CostFunction* cost_function,
                                   LossFunction* loss_function,
                                   double* x0,
                                   Ts*... xs) {
    const std::array<double*, sizeof...(Ts) + 1> parameter_blocks{{x0, xs...}};
    return AddResidualBlock(cost_function,
                            loss_function,
                            parameter_blocks.data(),
                            static_cast<int>(parameter_blocks.size()));
  }

  // Add a residual block by providing a vector of parameter blocks.
  // Ownership of cost_function and loss_function is handled as described above.
  ResidualBlockId AddResidualBlock(
      CostFunction* cost_function,
      LossFunction* loss_function,
      const std::vector<double*>& parameter_blocks);

  // Add a residual block by providing a pointer to the parameter block array
  // and the number of parameter blocks.
  // Ownership of cost_function and loss_function is handled as described above.
  ResidualBlockId AddResidualBlock(CostFunction* cost_function,
                                   LossFunction* loss_function,
                                   double* const* const parameter_blocks,
                                   int num_parameter_blocks);

  // Add a parameter block with appropriate size to the problem. Repeated calls
  // with the same arguments are ignored. Repeated calls with the same double
  // pointer but a different size will result in a crash.
  void AddParameterBlock(double* values, int size);

  // Add a parameter block with appropriate size and Manifold to the
  // problem. It is okay for manifold to be nullptr.
  //
  // Repeated calls with the same arguments are ignored. Repeated calls
  // with the same double pointer but a different size result in a crash
  // (unless Solver::Options::disable_all_safety_checks is set to true).
  //
  // Repeated calls with the same double pointer and size but different Manifold
  // is equivalent to calling SetManifold(manifold), i.e., any previously
  // associated Manifold object will be replaced with the manifold.
  void AddParameterBlock(double* values, int size, Manifold* manifold);

  // Remove a parameter block from the problem. The Manifold of the parameter
  // block, if it exists, will persist until the deletion of the problem
  // (similar to cost/loss functions in residual block removal). Any residual
  // blocks that depend on the parameter are also removed, as described above
  // in RemoveResidualBlock().
  //
  // If Problem::Options::enable_fast_removal is true, then the removal is fast
  // (almost constant time). Otherwise, removing a parameter block will incur a
  // scan of the entire Problem object.
  //
  // WARNING: Removing a residual or parameter block will destroy the implicit
  // ordering, rendering the Jacobian or residuals returned from the solver
  // uninterpretable. If you depend on the evaluated Jacobian, do not use
  // remove! This may change in a future release.
  void RemoveParameterBlock(const double* values);

  // Remove a residual block from the problem. Any parameters that the residual
  // block depends on are not removed. The cost and loss functions for the
  // residual block will not get deleted immediately; won't happen until the
  // problem itself is deleted.
  //
  // WARNING: Removing a residual or parameter block will destroy the implicit
  // ordering, rendering the Jacobian or residuals returned from the solver
  // uninterpretable. If you depend on the evaluated Jacobian, do not use
  // remove! This may change in a future release.
  void RemoveResidualBlock(ResidualBlockId residual_block);

  // Hold the indicated parameter block constant during optimization.
  void SetParameterBlockConstant(const double* values);

  // Allow the indicated parameter block to vary during optimization.
  void SetParameterBlockVariable(double* values);

  // Returns true if a parameter block is set constant, and false otherwise. A
  // parameter block may be set constant in two ways: either by calling
  // SetParameterBlockConstant or by associating a Manifold with a zero
  // dimensional tangent space with it.
  bool IsParameterBlockConstant(const double* values) const;

  // Set the Manifold for the parameter block. Calling SetManifold with nullptr
  // will clear any previously set Manifold for the parameter block.
  //
  // Repeated calls will result in any previously associated Manifold object to
  // be replaced with the manifold.
  //
  // The manifold is owned by the Problem by default (See Problem::Options to
  // override this behaviour).
  //
  // It is acceptable to set the same Manifold for multiple parameter blocks.
  void SetManifold(double* values, Manifold* manifold);

  // Get the Manifold object associated with this parameter block.
  //
  // If there is no Manifold object associated then nullptr is returned.
  const Manifold* GetManifold(const double* values) const;

  // Returns true if a Manifold is associated with this parameter block, false
  // otherwise.
  bool HasManifold(const double* values) const;

  // Set the lower/upper bound for the parameter at position "index".
  void SetParameterLowerBound(double* values, int index, double lower_bound);
  void SetParameterUpperBound(double* values, int index, double upper_bound);

  // Get the lower/upper bound for the parameter at position "index". If the
  // parameter is not bounded by the user, then its lower bound is
  // -std::numeric_limits<double>::max() and upper bound is
  // std::numeric_limits<double>::max().
  double GetParameterLowerBound(const double* values, int index) const;
  double GetParameterUpperBound(const double* values, int index) const;

  // Number of parameter blocks in the problem. Always equals
  // parameter_blocks().size() and parameter_block_sizes().size().
  int NumParameterBlocks() const;

  // The size of the parameter vector obtained by summing over the sizes of all
  // the parameter blocks.
  int NumParameters() const;

  // Number of residual blocks in the problem. Always equals
  // residual_blocks().size().
  int NumResidualBlocks() const;

  // The size of the residual vector obtained by summing over the sizes of all
  // of the residual blocks.
  int NumResiduals() const;

  // The size of the parameter block.
  int ParameterBlockSize(const double* values) const;

  // The dimension of the tangent space of the Manifold for the parameter block.
  // If there is no Manifold associated with this parameter block, then
  // ParameterBlockTangentSize = ParameterBlockSize.
  int ParameterBlockTangentSize(const double* values) const;

  // Is the given parameter block present in this problem or not?
  bool HasParameterBlock(const double* values) const;

  // Fills the passed parameter_blocks vector with pointers to the parameter
  // blocks currently in the problem. After this call, parameter_block.size() ==
  // NumParameterBlocks.
  void GetParameterBlocks(std::vector<double*>* parameter_blocks) const;

  // Fills the passed residual_blocks vector with pointers to the residual
  // blocks currently in the problem. After this call, residual_blocks.size() ==
  // NumResidualBlocks.
  void GetResidualBlocks(std::vector<ResidualBlockId>* residual_blocks) const;

  // Get all the parameter blocks that depend on the given residual block.
  void GetParameterBlocksForResidualBlock(
      const ResidualBlockId residual_block,
      std::vector<double*>* parameter_blocks) const;

  // Get the CostFunction for the given residual block.
  const CostFunction* GetCostFunctionForResidualBlock(
      const ResidualBlockId residual_block) const;

  // Get the LossFunction for the given residual block. Returns nullptr
  // if no loss function is associated with this residual block.
  const LossFunction* GetLossFunctionForResidualBlock(
      const ResidualBlockId residual_block) const;

  // Get all the residual blocks that depend on the given parameter block.
  //
  // If Problem::Options::enable_fast_removal is true, then getting the residual
  // blocks is fast and depends only on the number of residual
  // blocks. Otherwise, getting the residual blocks for a parameter block will
  // incur a scan of the entire Problem object.
  void GetResidualBlocksForParameterBlock(
      const double* values,
      std::vector<ResidualBlockId>* residual_blocks) const;

  // Options struct to control Problem::Evaluate.
  struct EvaluateOptions {
    // The set of parameter blocks for which evaluation should be
    // performed. This vector determines the order that parameter blocks occur
    // in the gradient vector and in the columns of the Jacobian matrix. If
    // parameter_blocks is empty, then it is assumed to be equal to the vector
    // containing ALL the parameter blocks. Generally speaking the parameter
    // blocks will occur in the order in which they were added to the
    // problem. But, this may change if the user removes any parameter blocks
    // from the problem.
    //
    // NOTE: This vector should contain the same pointers as the ones used to
    // add parameter blocks to the Problem. These parameter blocks should NOT
    // point to new memory locations. Bad things will happen otherwise.
    std::vector<double*> parameter_blocks;

    // The set of residual blocks to evaluate. This vector determines the order
    // in which the residuals occur, and how the rows of the Jacobian are
    // ordered. If residual_blocks is empty, then it is assumed to be equal to
    // the vector containing ALL the residual blocks. Generally speaking the
    // residual blocks will occur in the order in which they were added to the
    // problem. But, this may change if the user removes any residual blocks
    // from the problem.
    std::vector<ResidualBlockId> residual_blocks;

    // Even though the residual blocks in the problem may contain loss
    // functions, setting apply_loss_function to false will turn off the
    // application of the loss function to the output of the cost function. This
    // is of use for example if the user wishes to analyse the solution quality
    // by studying the distribution of residuals before and after the solve.
    bool apply_loss_function = true;

    int num_threads = 1;
  };

  // Evaluate Problem. Any of the output pointers can be nullptr. Which residual
  // blocks and parameter blocks are used is controlled by the EvaluateOptions
  // struct above.
  //
  // Note 1: The evaluation will use the values stored in the memory locations
  // pointed to by the parameter block pointers used at the time of the
  // construction of the problem. i.e.,
  //
  //   Problem problem;
  //   double x = 1;
  //   problem.AddResidualBlock(new MyCostFunction, nullptr, &x);
  //
  //   double cost = 0.0;
  //   problem.Evaluate(Problem::EvaluateOptions(), &cost,
  //                    nullptr, nullptr, nullptr);
  //
  // The cost is evaluated at x = 1. If you wish to evaluate the problem at x =
  // 2, then
  //
  //   x = 2;
  //   problem.Evaluate(Problem::EvaluateOptions(), &cost,
  //                    nullptr, nullptr, nullptr);
  //
  // is the way to do so.
  //
  // Note 2: If no Manifolds are used, then the size of the gradient vector (and
  // the number of columns in the Jacobian) is the sum of the sizes of all the
  // parameter blocks. If a parameter block has a Manifold, then it contributes
  // "TangentSize" entries to the gradient vector (and the number of columns in
  // the Jacobian).
  //
  // Note 3: This function cannot be called while the problem is being solved,
  // for example it cannot be called from an IterationCallback at the end of an
  // iteration during a solve.
  //
  // Note 4: If an EvaluationCallback is associated with the problem, then its
  // PrepareForEvaluation method will be called every time this method is called
  // with new_point = true.
  bool Evaluate(const EvaluateOptions& options,
                double* cost,
                std::vector<double>* residuals,
                std::vector<double>* gradient,
                CRSMatrix* jacobian);

  // Evaluates the residual block, storing the scalar cost in *cost, the
  // residual components in *residuals, and the Jacobians between the parameters
  // and residuals in jacobians[i], in row-major order.
  //
  // If residuals is nullptr, the residuals are not computed.
  //
  // If jacobians is nullptr, no Jacobians are computed. If jacobians[i] is
  // nullptr, then the Jacobian for that parameter block is not computed.
  //
  // It is not okay to request the Jacobian w.r.t a parameter block that is
  // constant.
  //
  // The return value indicates the success or failure. Even if the function
  // returns false, the caller should expect the output memory locations to have
  // been modified.
  //
  // The returned cost and Jacobians have had robustification and Manifold
  // applied already; for example, the Jacobian for a 4-dimensional quaternion
  // parameter using the "QuaternionParameterization" is num_residuals by 3
  // instead of num_residuals by 4.
  //
  // apply_loss_function as the name implies allows the user to switch the
  // application of the loss function on and off.
  //
  // If an EvaluationCallback is associated with the problem, then its
  // PrepareForEvaluation method will be called every time this method is called
  // with new_point = true. This conservatively assumes that the user may have
  // changed the parameter values since the previous call to evaluate / solve.
  // For improved efficiency, and only if you know that the parameter values
  // have not changed between calls, see
  // EvaluateResidualBlockAssumingParametersUnchanged().
  bool EvaluateResidualBlock(ResidualBlockId residual_block_id,
                             bool apply_loss_function,
                             double* cost,
                             double* residuals,
                             double** jacobians) const;

  // Same as EvaluateResidualBlock except that if an EvaluationCallback is
  // associated with the problem, then its PrepareForEvaluation method will be
  // called every time this method is called with new_point = false.
  //
  // This means, if an EvaluationCallback is associated with the problem then it
  // is the user's responsibility to call PrepareForEvaluation before calling
  // this method if necessary, i.e. iff the parameter values have been changed
  // since the last call to evaluate / solve.
  //
  // This is because, as the name implies, we assume that the parameter blocks
  // did not change since the last time PrepareForEvaluation was called (via
  // Solve, Evaluate or EvaluateResidualBlock).
  bool EvaluateResidualBlockAssumingParametersUnchanged(
      ResidualBlockId residual_block_id,
      bool apply_loss_function,
      double* cost,
      double* residuals,
      double** jacobians) const;

  // Returns reference to the options with which the Problem was constructed.
  const Options& options() const;

  // Returns pointer to Problem implementation
  internal::ProblemImpl* mutable_impl();

 private:
  std::unique_ptr<internal::ProblemImpl> impl_;
};

}  // namespace ceres

#include "ceres/internal/reenable_warnings.h"

#endif  // CERES_PUBLIC_PROBLEM_H_