10ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Ceres Solver - A fast non-linear least squares minimizer
20ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
30ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// http://code.google.com/p/ceres-solver/
40ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//
50ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Redistribution and use in source and binary forms, with or without
60ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// modification, are permitted provided that the following conditions are met:
70ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//
80ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// * Redistributions of source code must retain the above copyright notice,
90ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//   this list of conditions and the following disclaimer.
100ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// * Redistributions in binary form must reproduce the above copyright notice,
110ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//   this list of conditions and the following disclaimer in the documentation
120ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//   and/or other materials provided with the distribution.
130ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// * Neither the name of Google Inc. nor the names of its contributors may be
140ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//   used to endorse or promote products derived from this software without
150ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//   specific prior written permission.
160ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//
170ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
180ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
190ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
200ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
210ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
220ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
230ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
240ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
250ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
260ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
270ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// POSSIBILITY OF SUCH DAMAGE.
280ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//
290ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Author: sameeragarwal@google.com (Sameer Agarwal)
300ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
310ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include <limits>
320ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include <vector>
330ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "ceres/block_random_access_sparse_matrix.h"
340ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "ceres/internal/eigen.h"
350ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "glog/logging.h"
360ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "gtest/gtest.h"
370ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
380ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongnamespace ceres {
390ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongnamespace internal {
400ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
410ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus KongTEST(BlockRandomAccessSparseMatrix, GetCell) {
420ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  vector<int> blocks;
430ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  blocks.push_back(3);
440ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  blocks.push_back(4);
450ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  blocks.push_back(5);
460ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  const int num_rows = 3 + 4 + 5;
470ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
480ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  set< pair<int, int> > block_pairs;
490ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  int num_nonzeros = 0;
500ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  block_pairs.insert(make_pair(0, 0));
510ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  num_nonzeros += blocks[0] * blocks[0];
520ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
530ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  block_pairs.insert(make_pair(1, 1));
540ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  num_nonzeros += blocks[1] * blocks[1];
550ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
560ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  block_pairs.insert(make_pair(1, 2));
570ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  num_nonzeros += blocks[1] * blocks[2];
580ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
590ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  block_pairs.insert(make_pair(2, 0));
600ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  num_nonzeros += blocks[2] * blocks[0];
610ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
620ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  BlockRandomAccessSparseMatrix m(blocks, block_pairs);
630ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  EXPECT_EQ(m.num_rows(), num_rows);
640ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  EXPECT_EQ(m.num_cols(), num_rows);
650ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
660ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  for (set<pair<int, int> >::const_iterator it = block_pairs.begin();
670ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong       it != block_pairs.end();
680ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong       ++it) {
690ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    const int row_block_id = it->first;
700ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    const int col_block_id = it->second;
710ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    int row;
720ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    int col;
730ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    int row_stride;
740ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    int col_stride;
750ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    CellInfo* cell =  m.GetCell(row_block_id, col_block_id,
760ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong                                &row, &col,
770ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong                                &row_stride, &col_stride);
780ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    EXPECT_TRUE(cell != NULL);
790ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    EXPECT_EQ(row, 0);
800ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    EXPECT_EQ(col, 0);
810ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    EXPECT_EQ(row_stride, blocks[row_block_id]);
820ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    EXPECT_EQ(col_stride, blocks[col_block_id]);
830ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
840ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // Write into the block
850ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    MatrixRef(cell->values, row_stride, col_stride).block(
860ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong        row, col, blocks[row_block_id], blocks[col_block_id]) =
870ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong        (row_block_id + 1) * (col_block_id +1) *
880ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong        Matrix::Ones(blocks[row_block_id], blocks[col_block_id]);
890ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  }
900ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
910ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  const TripletSparseMatrix* tsm = m.matrix();
920ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  EXPECT_EQ(tsm->num_nonzeros(), num_nonzeros);
930ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  EXPECT_EQ(tsm->max_num_nonzeros(), num_nonzeros);
940ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
950ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  Matrix dense;
960ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  tsm->ToDenseMatrix(&dense);
970ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
980ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  double kTolerance = 1e-14;
990ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1000ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // (0,0)
1010ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  EXPECT_NEAR((dense.block(0, 0, 3, 3) - Matrix::Ones(3, 3)).norm(),
1020ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong              0.0,
1030ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong              kTolerance);
1040ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // (1,1)
1050ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  EXPECT_NEAR((dense.block(3, 3, 4, 4) - 2 * 2 * Matrix::Ones(4, 4)).norm(),
1060ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong              0.0,
1070ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong              kTolerance);
1080ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // (1,2)
1090ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  EXPECT_NEAR((dense.block(3, 3 + 4, 4, 5) - 2 * 3 * Matrix::Ones(4, 5)).norm(),
1100ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong              0.0,
1110ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong              kTolerance);
1120ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // (2,0)
1130ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  EXPECT_NEAR((dense.block(3 + 4, 0, 5, 3) - 3 * 1 * Matrix::Ones(5, 3)).norm(),
1140ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong              0.0,
1150ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong              kTolerance);
1160ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1170ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // There is nothing else in the matrix besides these four blocks.
1180ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  EXPECT_NEAR(dense.norm(), sqrt(9. + 16. * 16. + 36. * 20. + 9. * 15.),
1190ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong              kTolerance);
1200ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong}
1210ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1220ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// IntPairToLong is private, thus this fixture is needed to access and
1230ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// test it.
1240ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongclass BlockRandomAccessSparseMatrixTest : public ::testing::Test {
1250ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong public:
1260ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  virtual void SetUp() {
1270ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    vector<int> blocks;
1280ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    blocks.push_back(1);
1290ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    set< pair<int, int> > block_pairs;
1300ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    block_pairs.insert(make_pair(0, 0));
1310ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    m_.reset(new BlockRandomAccessSparseMatrix(blocks, block_pairs));
1320ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  }
1330ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1340ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  void CheckIntPair(int a, int b) {
1350ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    int64 value = m_->IntPairToLong(a, b);
1360ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    EXPECT_GT(value, 0) << "Overflow a = " << a << " b = " << b;
1370ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    EXPECT_GT(value, a) << "Overflow a = " << a << " b = " << b;
1380ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    EXPECT_GT(value, b) << "Overflow a = " << a << " b = " << b;
1390ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  }
1400ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1410ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong private:
1420ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  scoped_ptr<BlockRandomAccessSparseMatrix> m_;
1430ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong};
1440ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1450ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus KongTEST_F(BlockRandomAccessSparseMatrixTest, IntPairToLongOverflow) {
1460ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  CheckIntPair(numeric_limits<int>::max(), numeric_limits<int>::max());
1470ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong}
1480ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1490ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong}  // namespace internal
1500ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong}  // namespace ceres
151