// 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 "ceres/sparse_cholesky.h"

#include <numeric>

#include <vector>

#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 "ceres/internal/scoped_ptr.h"

#include "ceres/random.h"

#include "glog/logging.h"

#include "gtest/gtest.h"

namespace ceres {

namespace internal {

BlockSparseMatrix* CreateRandomFullRankMatrix(const int num_col_blocks,

                                              const int min_col_block_size,

                                              const int max_col_block_size,

                                              const double block_density) {

  // Create a random matrix

  BlockSparseMatrix::RandomMatrixOptions options;

  options.num_col_blocks = num_col_blocks;

  options.min_col_block_size = min_col_block_size;

  options.max_col_block_size = max_col_block_size;

  options.num_row_blocks = 2 * num_col_blocks;

  options.min_row_block_size = 1;

  options.max_row_block_size = max_col_block_size;

  options.block_density = block_density;

  scoped_ptr<BlockSparseMatrix> random_matrix(

      BlockSparseMatrix::CreateRandomMatrix(options));

  // Add a diagonal block sparse matrix to make it full rank.

  Vector diagonal = Vector::Ones(random_matrix->num_cols());

  scoped_ptr<BlockSparseMatrix> block_diagonal(

      BlockSparseMatrix::CreateDiagonalMatrix(

          diagonal.data(), random_matrix->block_structure()->cols));

  random_matrix->AppendRows(*block_diagonal);

  return random_matrix.release();

}

bool ComputeExpectedSolution(const CompressedRowSparseMatrix& lhs,

                             const Vector& rhs,

                             Vector* solution) {

  Matrix eigen_lhs;

  lhs.ToDenseMatrix(&eigen_lhs);

  if (lhs.storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {

    Matrix full_lhs = eigen_lhs.selfadjointView<Eigen::Upper>();

    Eigen::LLT<Matrix, Eigen::Upper> llt =

        eigen_lhs.selfadjointView<Eigen::Upper>().llt();

    if (llt.info() != Eigen::Success) {

      return false;

    }

    *solution = llt.solve(rhs);

    return (llt.info() == Eigen::Success);

  }

  Matrix full_lhs = eigen_lhs.selfadjointView<Eigen::Lower>();

  Eigen::LLT<Matrix, Eigen::Lower> llt =

      eigen_lhs.selfadjointView<Eigen::Lower>().llt();

  if (llt.info() != Eigen::Success) {

    return false;

  }

  *solution = llt.solve(rhs);

  return (llt.info() == Eigen::Success);

}

void SparseCholeskySolverUnitTest(

    const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,

    const OrderingType ordering_type,

    const bool use_block_structure,

    const int num_blocks,

    const int min_block_size,

    const int max_block_size,

    const double block_density) {

  scoped_ptr<SparseCholesky> sparse_cholesky(SparseCholesky::Create(

      sparse_linear_algebra_library_type, ordering_type));

  const CompressedRowSparseMatrix::StorageType storage_type =

      sparse_cholesky->StorageType();

  scoped_ptr<BlockSparseMatrix> m(CreateRandomFullRankMatrix(

      num_blocks, min_block_size, max_block_size, block_density));

  scoped_ptr<InnerProductComputer> inner_product_computer(

      InnerProductComputer::Create(*m, storage_type));

  inner_product_computer->Compute();

  CompressedRowSparseMatrix* lhs = inner_product_computer->mutable_result();

  if (!use_block_structure) {

    lhs->mutable_row_blocks()->clear();

    lhs->mutable_col_blocks()->clear();

  }

  Vector rhs = Vector::Random(lhs->num_rows());

  Vector expected(lhs->num_rows());

  Vector actual(lhs->num_rows());

  EXPECT_TRUE(ComputeExpectedSolution(*lhs, rhs, &expected));

  std::string message;

  EXPECT_EQ(sparse_cholesky->FactorAndSolve(

                lhs, rhs.data(), actual.data(), &message),

            LINEAR_SOLVER_SUCCESS);

  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;

}

typedef ::testing::tuple<SparseLinearAlgebraLibraryType, OrderingType, 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) == AMD ? "AMD" : "NATURAL") << "_"

     << (::testing::get<2>(param) ? "UseBlockStructure" : "NoBlockStructure");

  return ss.str();

}

class SparseCholeskyTest : public ::testing::TestWithParam<Param> {};

TEST_P(SparseCholeskyTest, FactorAndSolve) {

  SetRandomState(2982);

  const int kMinNumBlocks = 1;

  const int kMaxNumBlocks = 10;

  const int kNumTrials = 10;

  const int kMinBlockSize = 1;

  const int kMaxBlockSize = 5;

  for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks;

       ++num_blocks) {

    for (int trial = 0; trial < kNumTrials; ++trial) {

      const double block_density = std::max(0.1, RandDouble());

      Param param = GetParam();

      SparseCholeskySolverUnitTest(::testing::get<0>(param),

                                   ::testing::get<1>(param),

                                   ::testing::get<2>(param),

                                   num_blocks,

                                   kMinBlockSize,

                                   kMaxBlockSize,

                                   block_density);

    }

  }

}

#ifndef CERES_NO_SUITESPARSE

INSTANTIATE_TEST_CASE_P(SuiteSparseCholesky,

                        SparseCholeskyTest,

                        ::testing::Combine(::testing::Values(SUITE_SPARSE),

                                           ::testing::Values(AMD, NATURAL),

                                           ::testing::Values(true, false)),

                        ParamInfoToString);

#endif

#ifndef CERES_NO_CXSPARSE

INSTANTIATE_TEST_CASE_P(CXSparseCholesky,

                        SparseCholeskyTest,

                        ::testing::Combine(::testing::Values(CX_SPARSE),

                                           ::testing::Values(AMD, NATURAL),

                                           ::testing::Values(true, false)),

                        ParamInfoToString);

#endif

#ifdef CERES_USE_EIGEN_SPARSE

INSTANTIATE_TEST_CASE_P(EigenSparseCholesky,

                        SparseCholeskyTest,

                        ::testing::Combine(::testing::Values(EIGEN_SPARSE),

                                           ::testing::Values(AMD, NATURAL),

                                           ::testing::Values(true, false)),

                        ParamInfoToString);

#endif

}  // namespace internal

}  // namespace ceres