// Ceres Solver - A fast non-linear least squares minimizer

// Copyright 2017 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)

#ifndef CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_

#define CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_

#include <vector>

#include "ceres/block_sparse_matrix.h"

#include "ceres/compressed_row_sparse_matrix.h"

#include "ceres/internal/scoped_ptr.h"

namespace ceres {

namespace internal {

// This class is used to repeatedly compute the inner product

//

//   result = m' * m

//

// where the sparsity structure of m remains constant across calls.

//

// Upon creation, the class computes and caches information needed to

// compute v, and then uses it to efficiently compute the product

// every time InnerProductComputer::Compute is called.

//

// See sparse_normal_cholesky_solver.cc for example usage.

//

// Note that the result matrix is a block upper or lower-triangular

// matrix, i.e., it will contain entries in the upper or lower

// triangular part of the matrix corresponding to the block that occur

// along its diagonal.

//

// This is not a problem as sparse linear algebra libraries can ignore

// these entries with ease and the space used is minimal/linear in the

// size of the matrices.

class InnerProductComputer {

 public:

  // Factory

  //

  // m is the input matrix

  //

  // Since m' * m is a symmetric matrix, we only compute half of the

  // matrix and the value of storage_type which must be

  // UPPER_TRIANGULAR or LOWER_TRIANGULAR determines which half is

  // computed.

  //

  // The user must ensure that the matrix m is valid for the life time

  // of this object.

  static InnerProductComputer* Create(

      const BlockSparseMatrix& m,

      CompressedRowSparseMatrix::StorageType storage_type);

  // This factory method allows the user control over range of row

  // blocks of m that should be used to compute the inner product.

  //

  // a = m(start_row_block : end_row_block, :);

  // result = a' * a;

  static InnerProductComputer* Create(

      const BlockSparseMatrix& m,

      int start_row_block,

      int end_row_block,

      CompressedRowSparseMatrix::StorageType storage_type);

  // Update result_ to be numerically equal to m' * m.

  void Compute();

  // Accessors for the result containing the inner product.

  //

  // Compute must be called before accessing this result for

  // the first time.

  const CompressedRowSparseMatrix& result() const { return *result_; }

  CompressedRowSparseMatrix* mutable_result() const { return result_.get(); }

 private:

  // A ProductTerm is a term in the block inner product of a matrix

  // with itself.

  struct ProductTerm {

    ProductTerm(const int row, const int col, const int index)

        : row(row), col(col), index(index) {}

    bool operator<(const ProductTerm& right) const {

      if (row == right.row) {

        if (col == right.col) {

          return index < right.index;

        }

        return col < right.col;

      }

      return row < right.row;

    }

    int row;

    int col;

    int index;

  };

  InnerProductComputer(const BlockSparseMatrix& m,

                       int start_row_block,

                       int end_row_block);

  void Init(CompressedRowSparseMatrix::StorageType storage_type);

  CompressedRowSparseMatrix* CreateResultMatrix(

      const CompressedRowSparseMatrix::StorageType storage_type,

      int num_nonzeros);

  int ComputeNonzeros(const std::vector<ProductTerm>& product_terms,

                      std::vector<int>* row_block_nnz);

  void ComputeOffsetsAndCreateResultMatrix(

      const CompressedRowSparseMatrix::StorageType storage_type,

      const std::vector<ProductTerm>& product_terms);

  const BlockSparseMatrix& m_;

  const int start_row_block_;

  const int end_row_block_;

  scoped_ptr<CompressedRowSparseMatrix> result_;

  // For each term in the inner product, result_offsets_ contains the

  // location in the values array of the result_ matrix where it

  // should be stored.

  //

  // This is the principal look up table that allows this class to

  // compute the inner product fast.

  std::vector<int> result_offsets_;

};

}  // namespace internal

}  // namespace ceres

#endif  // CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_