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

#include <memory>
#include "ceres/subset_preconditioner.h"
#include "Eigen/Dense"
#include "Eigen/SparseCore"
#include "ceres/block_sparse_matrix.h"
#include "ceres/compressed_row_sparse_matrix.h"
#include "ceres/inner_product_computer.h"
#include "ceres/internal/eigen.h"
#include "glog/logging.h"
#include "gtest/gtest.h"

namespace ceres {
namespace internal {

// TODO(sameeragarwal): Refactor the following two functions out of
// here and sparse_cholesky_test.cc into a more suitable place.
template <int UpLoType>
bool SolveLinearSystemUsingEigen(const Matrix& lhs,
                                 const Vector rhs,
                                 Vector* solution) {
  Eigen::LLT<Matrix, UpLoType> llt = lhs.selfadjointView<UpLoType>().llt();
  if (llt.info() != Eigen::Success) {
    return false;
  }
  *solution = llt.solve(rhs);
  return (llt.info() == Eigen::Success);
}

// Use Eigen's Dense Cholesky solver to compute the solution to a
// sparse linear system.
bool ComputeExpectedSolution(const CompressedRowSparseMatrix& lhs,
                             const Vector& rhs,
                             Vector* solution) {
  Matrix dense_triangular_lhs;
  lhs.ToDenseMatrix(&dense_triangular_lhs);
  if (lhs.storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
    Matrix full_lhs = dense_triangular_lhs.selfadjointView<Eigen::Upper>();
    return SolveLinearSystemUsingEigen<Eigen::Upper>(full_lhs, rhs, solution);
  }
  return SolveLinearSystemUsingEigen<Eigen::Lower>(
      dense_triangular_lhs, rhs, solution);
}

typedef ::testing::tuple<SparseLinearAlgebraLibraryType, bool> Param;

std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
  Param param = info.param;
  std::stringstream ss;
  ss << SparseLinearAlgebraLibraryTypeToString(::testing::get<0>(param)) << "_"
     << (::testing::get<1>(param) ? "Diagonal" : "NoDiagonal");
  return ss.str();
}

class SubsetPreconditionerTest : public ::testing::TestWithParam<Param> {
 protected:
  virtual void SetUp() {
    BlockSparseMatrix::RandomMatrixOptions options;
    options.num_col_blocks = 4;
    options.min_col_block_size = 1;
    options.max_col_block_size = 4;
    options.num_row_blocks = 8;
    options.min_row_block_size = 1;
    options.max_row_block_size = 4;
    options.block_density = 0.9;

    m_.reset(BlockSparseMatrix::CreateRandomMatrix(options));
    start_row_block_ = m_->block_structure()->rows.size();

    // Ensure that the bottom part of the matrix has the same column
    // block structure.
    options.col_blocks = m_->block_structure()->cols;
    b_.reset(BlockSparseMatrix::CreateRandomMatrix(options));
    m_->AppendRows(*b_);

    // Create a Identity block diagonal matrix with the same column
    // block structure.
    diagonal_ = Vector::Ones(m_->num_cols());
    block_diagonal_.reset(BlockSparseMatrix::CreateDiagonalMatrix(
        diagonal_.data(), b_->block_structure()->cols));

    // Unconditionally add the block diagonal to the matrix b_,
    // because either it is either part of b_ to make it full rank, or
    // we pass the same diagonal matrix later as the parameter D. In
    // either case the preconditioner matrix is b_' b + D'D.
    b_->AppendRows(*block_diagonal_);
    inner_product_computer_.reset(InnerProductComputer::Create(
        *b_, CompressedRowSparseMatrix::UPPER_TRIANGULAR));
    inner_product_computer_->Compute();
  }

  std::unique_ptr<BlockSparseMatrix> m_;
  std::unique_ptr<BlockSparseMatrix> b_;
  std::unique_ptr<BlockSparseMatrix> block_diagonal_;
  std::unique_ptr<InnerProductComputer> inner_product_computer_;
  std::unique_ptr<Preconditioner> preconditioner_;
  Vector diagonal_;
  int start_row_block_;
};

TEST_P(SubsetPreconditionerTest, foo) {
  Param param = GetParam();
  Preconditioner::Options options;
  options.subset_preconditioner_start_row_block = start_row_block_;
  options.sparse_linear_algebra_library_type = ::testing::get<0>(param);
  preconditioner_.reset(new SubsetPreconditioner(options, *m_));

  const bool with_diagonal = ::testing::get<1>(param);
  if (!with_diagonal) {
    m_->AppendRows(*block_diagonal_);
  }

  EXPECT_TRUE(
      preconditioner_->Update(*m_, with_diagonal ? diagonal_.data() : NULL));

  // Repeatedly apply the preconditioner to random vectors and check
  // that the preconditioned value is the same as one obtained by
  // solving the linear system directly.
  for (int i = 0; i < 5; ++i) {
    CompressedRowSparseMatrix* lhs = inner_product_computer_->mutable_result();
    Vector rhs = Vector::Random(lhs->num_rows());
    Vector expected(lhs->num_rows());
    EXPECT_TRUE(ComputeExpectedSolution(*lhs, rhs, &expected));

    Vector actual(lhs->num_rows());
    preconditioner_->RightMultiply(rhs.data(), actual.data());

    Matrix eigen_lhs;
    lhs->ToDenseMatrix(&eigen_lhs);
    EXPECT_NEAR((actual - expected).norm() / actual.norm(),
                0.0,
                std::numeric_limits<double>::epsilon() * 10)
        << "\n"
        << eigen_lhs << "\n"
        << expected.transpose() << "\n"
        << actual.transpose();
  }
}

#ifndef CERES_NO_SUITESPARSE
INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithSuiteSparse,
                        SubsetPreconditionerTest,
                        ::testing::Combine(::testing::Values(SUITE_SPARSE),
                                           ::testing::Values(true, false)),
                        ParamInfoToString);
#endif

#ifndef CERES_NO_CXSPARSE
INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithCXSparse,
                        SubsetPreconditionerTest,
                        ::testing::Combine(::testing::Values(CX_SPARSE),
                                           ::testing::Values(true, false)),
                        ParamInfoToString);
#endif

#ifdef CERES_USE_EIGEN_SPARSE
INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithEigenSparse,
                        SubsetPreconditionerTest,
                        ::testing::Combine(::testing::Values(EIGEN_SPARSE),
                                           ::testing::Values(true, false)),
                        ParamInfoToString);
#endif

}  // namespace internal
}  // namespace ceres