Attempt to compile boost unit test.
This commit is contained in:
Dongryeol Lee
@@ -4,6 +4,7 @@ cmake_minimum_required(VERSION 2.8)
|
||||
set(DIRS
|
||||
# compression
|
||||
# fast_multipole_method ## awaiting resolution of #10 (libint dependency)
|
||||
gp_regression
|
||||
# kde ## does not compile
|
||||
linear_algebra
|
||||
linear_regression
|
||||
|
||||
@@ -0,0 +1,11 @@
|
||||
cmake_minimum_required(VERSION 2.8)
|
||||
|
||||
# test executable
|
||||
add_executable(bilinear_form_test
|
||||
EXCLUDE_FROM_ALL
|
||||
bilinear_form_test.cc
|
||||
)
|
||||
# link dependencies of test executable
|
||||
target_link_libraries(bilinear_form_test
|
||||
fastlib
|
||||
)
|
||||
@@ -0,0 +1,128 @@
|
||||
/** @author Dongryeol Lee
|
||||
*
|
||||
* @file bilinear_form_estimator.h
|
||||
*/
|
||||
|
||||
#ifndef FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_BILINEAR_FORM_ESTIMATOR_H
|
||||
#define FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_BILINEAR_FORM_ESTIMATOR_H
|
||||
|
||||
#include <vector>
|
||||
#include "fastlib/la/matrix.h"
|
||||
#include "fastlib/la/uselapack.h"
|
||||
#include "linear_operator.h"
|
||||
|
||||
namespace fl {
|
||||
namespace ml {
|
||||
|
||||
class SquareRootTransformation {
|
||||
public:
|
||||
static double Transform(double val_in) {
|
||||
return sqrt(val_in);
|
||||
}
|
||||
};
|
||||
|
||||
class IdentityTransformation {
|
||||
public:
|
||||
static double Transform(double val_in) {
|
||||
return val_in;
|
||||
}
|
||||
};
|
||||
|
||||
class InverseTransformation {
|
||||
public:
|
||||
static double Transform(double val_in) {
|
||||
return 1.0 / val_in;
|
||||
}
|
||||
};
|
||||
|
||||
class LogTransformation {
|
||||
public:
|
||||
static double Transform(double val_in) {
|
||||
return log(val_in);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename TransformationType>
|
||||
class BilinearFormEstimator {
|
||||
|
||||
private:
|
||||
|
||||
class TridiagonalLinearOperator {
|
||||
private:
|
||||
const std::vector<double> *diagonal_entries_;
|
||||
const std::vector<double> *offdiagonal_entries_;
|
||||
|
||||
public:
|
||||
|
||||
int n_rows() const;
|
||||
|
||||
int n_cols() const;
|
||||
|
||||
double get(int row, int col) const;
|
||||
|
||||
const std::vector<double> *diagonal_entries() const;
|
||||
|
||||
TridiagonalLinearOperator(
|
||||
const std::vector<double> &diagonal_entries_in,
|
||||
const std::vector<double> &offdiagonal_entries_in);
|
||||
|
||||
int Apply(const Epetra_MultiVector &vecs,
|
||||
Epetra_MultiVector &prods) const;
|
||||
|
||||
void PrintDebug(const char *name = "", FILE *stream = stderr) const;
|
||||
};
|
||||
|
||||
private:
|
||||
|
||||
#ifdef EPETRA_MPI
|
||||
Epetra_MpiComm comm_;
|
||||
#else
|
||||
Epetra_SerialComm comm_;
|
||||
#endif
|
||||
|
||||
Anasazi::LinearOperator *linear_operator_;
|
||||
|
||||
const Epetra_Map *map_;
|
||||
|
||||
private:
|
||||
void AddExpert_(double scalar,
|
||||
const Epetra_MultiVector &source,
|
||||
Epetra_MultiVector *destination) const;
|
||||
|
||||
void Scale_(double scalar,
|
||||
const Epetra_MultiVector &source,
|
||||
Epetra_MultiVector *destination) const;
|
||||
|
||||
double Dot_(const Epetra_MultiVector &first_vec,
|
||||
const Epetra_MultiVector &second_vec) const;
|
||||
|
||||
template<typename LinearOperatorType>
|
||||
double ComputeQuadraticForm_(
|
||||
int num_iterations,
|
||||
const GenVector<double> &starting_vector,
|
||||
const LinearOperatorType &linear_operator_in,
|
||||
const Epetra_Map &map_in,
|
||||
int level,
|
||||
bool *broke_down);
|
||||
|
||||
public:
|
||||
|
||||
BilinearFormEstimator();
|
||||
|
||||
Anasazi::LinearOperator *linear_operator();
|
||||
|
||||
void Init(Anasazi::LinearOperator *linear_operator_in);
|
||||
|
||||
double Compute(
|
||||
const GenVector<double> &left_argument,
|
||||
const GenVector<double> &right_argument,
|
||||
bool naive_compute);
|
||||
|
||||
double Compute(const GenVector<double> &argument);
|
||||
|
||||
double NaiveCompute(const GenVector<double> &argument);
|
||||
};
|
||||
};
|
||||
};
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,356 @@
|
||||
/** @author Dongryeol Lee
|
||||
*
|
||||
* @file bilinear_form_estimator_dev.h
|
||||
*/
|
||||
|
||||
#ifndef FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_BILINEAR_FORM_ESTIMATOR_DEV_H
|
||||
#define FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_BILINEAR_FORM_ESTIMATOR_DEV_H
|
||||
|
||||
#include <vector>
|
||||
#include "fastlib/la/matrix.h"
|
||||
#include "fastlib/la/la.h"
|
||||
#include "fastlib/la/uselapack.h"
|
||||
#include "bilinear_form_estimator.h"
|
||||
|
||||
namespace fl {
|
||||
namespace ml {
|
||||
|
||||
template<typename TransformationType>
|
||||
void BilinearFormEstimator<TransformationType>::
|
||||
TridiagonalLinearOperator::PrintDebug(
|
||||
const char *name, FILE *stream) const {
|
||||
fprintf(stream, "----- MATRIX ------: %s\n", name);
|
||||
for (int r = 0; r < this->n_rows(); r++) {
|
||||
for (int c = 0; c < this->n_cols(); c++) {
|
||||
fprintf(stream, "%+3.3f ", this->get(r, c));
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
int BilinearFormEstimator<TransformationType>::
|
||||
TridiagonalLinearOperator::n_rows() const {
|
||||
return diagonal_entries_->size();
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
int BilinearFormEstimator<TransformationType>::
|
||||
TridiagonalLinearOperator::n_cols() const {
|
||||
return diagonal_entries_->size();
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
double BilinearFormEstimator<TransformationType>::
|
||||
TridiagonalLinearOperator::get(int row, int col) const {
|
||||
if (row == col) {
|
||||
return (*diagonal_entries_)[row];
|
||||
}
|
||||
else if (row == col + 1 || col == row + 1) {
|
||||
return (*offdiagonal_entries_)[ std::min(row, col)];
|
||||
}
|
||||
else {
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
const std::vector<double> *BilinearFormEstimator<TransformationType>::
|
||||
TridiagonalLinearOperator::diagonal_entries() const {
|
||||
|
||||
return diagonal_entries_;
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
BilinearFormEstimator<TransformationType>::
|
||||
TridiagonalLinearOperator::TridiagonalLinearOperator(
|
||||
const std::vector<double> &diagonal_entries_in,
|
||||
const std::vector<double> &offdiagonal_entries_in) {
|
||||
|
||||
diagonal_entries_ = &diagonal_entries_in;
|
||||
offdiagonal_entries_ = &offdiagonal_entries_in;
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
int BilinearFormEstimator<TransformationType>::
|
||||
TridiagonalLinearOperator::Apply(
|
||||
const Epetra_MultiVector &vecs,
|
||||
Epetra_MultiVector &prods) const {
|
||||
|
||||
prods.PutScalar(0);
|
||||
|
||||
for (int j = 0; j < diagonal_entries_->size(); j++) {
|
||||
|
||||
for (int k = 0; k < vecs.NumVectors(); k++) {
|
||||
|
||||
// Apply the diagonal entry.
|
||||
prods.Pointers()[k][j] += ((*diagonal_entries_)[j]) *
|
||||
vecs.Pointers()[k][j];
|
||||
|
||||
// Apply the lower diagonal entry.
|
||||
if (j > 0) {
|
||||
prods.Pointers()[k][j - 1] += ((*offdiagonal_entries_)[j - 1]) *
|
||||
vecs.Pointers()[k][j - 1];
|
||||
}
|
||||
|
||||
// Apply the upper diagonal entry.
|
||||
if (j < diagonal_entries_->size() - 1) {
|
||||
prods.Pointers()[k][j + 1] += ((*offdiagonal_entries_)[j]) *
|
||||
vecs.Pointers()[k][j + 1];
|
||||
}
|
||||
}
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
#ifdef EPETRA_MPI
|
||||
template<typename TransformationType>
|
||||
BilinearFormEstimator<TransformationType>::BilinearFormEstimator(): comm_(MPI_COMM_WORLD) {
|
||||
|
||||
}
|
||||
#else
|
||||
template<typename TransformationType>
|
||||
BilinearFormEstimator<TransformationType>::BilinearFormEstimator() {
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
template<typename TransformationType>
|
||||
Anasazi::LinearOperator *BilinearFormEstimator<TransformationType>::
|
||||
linear_operator() {
|
||||
return linear_operator_;
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
void BilinearFormEstimator<TransformationType>::Init(
|
||||
Anasazi::LinearOperator *linear_operator_in) {
|
||||
|
||||
linear_operator_ = linear_operator_in;
|
||||
map_ = &(linear_operator_->OperatorDomainMap());
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
double BilinearFormEstimator<TransformationType>::Dot_(
|
||||
const Epetra_MultiVector &first_vec,
|
||||
const Epetra_MultiVector &second_vec) const {
|
||||
|
||||
double dot_product = 0;
|
||||
for (int i = 0; i < first_vec.GlobalLength(); i++) {
|
||||
dot_product += first_vec.Pointers()[0][i] * second_vec.Pointers()[0][i];
|
||||
}
|
||||
return dot_product;
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
void BilinearFormEstimator<TransformationType>::AddExpert_(double scalar,
|
||||
const Epetra_MultiVector &source,
|
||||
Epetra_MultiVector *destination) const {
|
||||
for (int i = 0; i < source.GlobalLength(); i++) {
|
||||
destination->Pointers()[0][i] =
|
||||
destination->Pointers()[0][i] + scalar * source.Pointers()[0][i];
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
void BilinearFormEstimator<TransformationType>::Scale_(double scalar,
|
||||
const Epetra_MultiVector &source,
|
||||
Epetra_MultiVector *destination) const {
|
||||
for (int i = 0; i < source.GlobalLength(); i++) {
|
||||
destination->Pointers()[0][i] = scalar * source.Pointers()[0][i];
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
template<typename LinearOperatorType>
|
||||
double BilinearFormEstimator<TransformationType>::ComputeQuadraticForm_(
|
||||
int num_iterations,
|
||||
const GenVector<double> &starting_vector,
|
||||
const LinearOperatorType &linear_operator_in,
|
||||
const Epetra_Map &map_in,
|
||||
int level,
|
||||
bool *break_down) {
|
||||
|
||||
// The threshold for determining the convergence.
|
||||
const double convergence_threshold = 1e-7;
|
||||
|
||||
// If it is a 1 by 1 matrix, then apply the transformation.
|
||||
if (linear_operator_in.n_rows() == 1 && linear_operator_in.n_cols() == 1) {
|
||||
return TransformationType::Transform(linear_operator_in.get(0, 0));
|
||||
}
|
||||
|
||||
// The diagonal entries and the offdiagonal entries with the wrapper
|
||||
// class around it.
|
||||
std::vector<double> diagonal_entries;
|
||||
std::vector<double> offdiagonal_entries;
|
||||
TridiagonalLinearOperator tridiagonal_linear_operator(
|
||||
diagonal_entries, offdiagonal_entries);
|
||||
|
||||
// The basis vector in the previous iteration.
|
||||
Epetra_MultiVector previous_vector(map_in, 1);
|
||||
previous_vector.PutScalar(0);
|
||||
|
||||
// The basis vector in the current iteration.
|
||||
Epetra_MultiVector current_vector(map_in, 1);
|
||||
for (int i = 0; i < starting_vector.length(); i++) {
|
||||
current_vector.Pointers()[0][i] = starting_vector[i];
|
||||
}
|
||||
|
||||
// A temporary vector used for matrix-vector multiplication.
|
||||
Epetra_MultiVector residual_vector(map_in, 1);
|
||||
residual_vector.PutScalar(0);
|
||||
|
||||
// A temporary vector used for denoting the unit vector of varying
|
||||
// dimension.
|
||||
GenVector<double> unit_vector;
|
||||
unit_vector.Init(starting_vector.length());
|
||||
unit_vector.SetZero();
|
||||
unit_vector[0] = 1;
|
||||
|
||||
// The old bilinear estimate.
|
||||
double old_bilinear_estimate = std::numeric_limits<double>::max();
|
||||
|
||||
for (int j = 0; j < num_iterations; j++) {
|
||||
|
||||
linear_operator_in.Apply(current_vector, residual_vector);
|
||||
if (j > 0) {
|
||||
AddExpert_(- (offdiagonal_entries[j - 1]),
|
||||
previous_vector, &residual_vector);
|
||||
}
|
||||
|
||||
// The dot product with the residual and the current basis
|
||||
// vector. Compute the off-diagonal and the diagonal entries
|
||||
// in this iteration.
|
||||
double alpha_j = Dot_(residual_vector, current_vector);
|
||||
diagonal_entries.push_back(alpha_j);
|
||||
AddExpert_(- alpha_j, current_vector, &residual_vector);
|
||||
double beta_j_plus_one =
|
||||
sqrt(Dot_(residual_vector, residual_vector));
|
||||
|
||||
// Add in the offdiagonal entry.
|
||||
if (fabs(beta_j_plus_one) <= convergence_threshold) {
|
||||
*break_down = true;
|
||||
break;
|
||||
}
|
||||
offdiagonal_entries.push_back(beta_j_plus_one);
|
||||
|
||||
// Take the current tridiagonal decomposition and estimate the
|
||||
// bilinear form. It is essential that Epetra_Map is constructed
|
||||
// here right before the recursive call.
|
||||
GenVector<double> unit_vector_alias;
|
||||
unit_vector_alias.Alias(unit_vector.ptr(),
|
||||
diagonal_entries.size());
|
||||
Epetra_Map tridiagonal_linear_operator_map(
|
||||
unit_vector_alias.length(), 0, comm_);
|
||||
|
||||
bool subcase_break_down = false;
|
||||
double bilinear_estimate =
|
||||
ComputeQuadraticForm_(unit_vector_alias.length(), unit_vector_alias,
|
||||
tridiagonal_linear_operator,
|
||||
tridiagonal_linear_operator_map, level + 1,
|
||||
&subcase_break_down);
|
||||
|
||||
// Check whether it converged.
|
||||
if (subcase_break_down) {
|
||||
*break_down = true;
|
||||
break;
|
||||
}
|
||||
if (fabs(old_bilinear_estimate - bilinear_estimate) <=
|
||||
convergence_threshold && false) {
|
||||
*break_down = false;
|
||||
break;
|
||||
}
|
||||
else {
|
||||
old_bilinear_estimate = bilinear_estimate;
|
||||
}
|
||||
|
||||
// Update the previous vector and the next vector.
|
||||
for (int i = 0; i < current_vector.GlobalLength(); i++) {
|
||||
previous_vector.Pointers()[0][i] = current_vector.Pointers()[0][i];
|
||||
}
|
||||
Scale_(1.0 / beta_j_plus_one, residual_vector, ¤t_vector);
|
||||
}
|
||||
|
||||
return old_bilinear_estimate;
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
double BilinearFormEstimator<TransformationType>::Compute(
|
||||
const GenVector<double> &left_argument,
|
||||
const GenVector<double> &right_argument,
|
||||
bool naive_compute) {
|
||||
|
||||
// Use the formula: $u^T f(A) v = 0.25 * (y^T f(A) y - z^T f(A) z )
|
||||
// $ where $y = u + v$ and $z = u - v$.
|
||||
GenVector<double> sum, difference;
|
||||
la::AddInit(
|
||||
left_argument, right_argument, &sum);
|
||||
la::SubInit(
|
||||
right_argument, left_argument, &difference);
|
||||
double estimate =
|
||||
(naive_compute) ?
|
||||
0.25 * (NaiveCompute(sum) - NaiveCompute(difference)) :
|
||||
0.25 * (Compute(sum) - Compute(difference));
|
||||
return estimate;
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
double BilinearFormEstimator<TransformationType>::NaiveCompute(
|
||||
const GenVector<double> &argument) {
|
||||
|
||||
// The naive estimate to return.
|
||||
double naive_bilinear_estimate = 0;
|
||||
|
||||
// Compute the kernel matrix, and its eigendecomposition.
|
||||
GenMatrix<double> kernel_matrix;
|
||||
GenVector<double> eigenvalues;
|
||||
GenMatrix<double> eigenvectors;
|
||||
GenMatrix<double> eigenvectors_transposed;
|
||||
kernel_matrix.Init(linear_operator_->n_rows(), linear_operator_->n_cols());
|
||||
for (int j = 0; j < linear_operator_->n_cols(); j++) {
|
||||
for (int i = 0; i < linear_operator_->n_rows(); i++) {
|
||||
kernel_matrix.set(i, j, linear_operator_->get(i, j));
|
||||
}
|
||||
}
|
||||
la::SVDInit(
|
||||
kernel_matrix, &eigenvalues, &eigenvectors, &eigenvectors_transposed);
|
||||
eigenvalues.PrintDebug();
|
||||
|
||||
// Project the argument to the eigenspace.
|
||||
GenVector<double> projected_vector;
|
||||
la::MulInit(
|
||||
eigenvectors_transposed, argument, &projected_vector);
|
||||
for (int i = 0; i < eigenvalues.length(); i++) {
|
||||
naive_bilinear_estimate += math::Sqr(projected_vector[i]) *
|
||||
TransformationType::Transform(eigenvalues[i]);
|
||||
}
|
||||
return naive_bilinear_estimate;
|
||||
}
|
||||
|
||||
template<typename TransformationType>
|
||||
double BilinearFormEstimator<TransformationType>::Compute(
|
||||
const GenVector<double> &argument) {
|
||||
|
||||
// Pass in the normalized unit vector to the quadratic form
|
||||
// computation and correct it afterwards.
|
||||
GenVector<double> normalized_argument;
|
||||
double length = la::LengthEuclidean(argument);
|
||||
if (length > 0) {
|
||||
la::ScaleInit(
|
||||
1.0 / length, argument, &normalized_argument);
|
||||
}
|
||||
else {
|
||||
normalized_argument.Copy(argument);
|
||||
}
|
||||
bool break_down = false;
|
||||
return math::Sqr(length) * ComputeQuadraticForm_(
|
||||
normalized_argument.length(),
|
||||
normalized_argument,
|
||||
*linear_operator_,
|
||||
*map_,
|
||||
0,
|
||||
&break_down);
|
||||
}
|
||||
};
|
||||
};
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,84 @@
|
||||
/** @author Dongryeol Lee
|
||||
*
|
||||
* @file bilinear_form.test.cc
|
||||
*/
|
||||
|
||||
#undef BOOST_ALL_DYN_LINK
|
||||
#include "fastlib/fastlib.h"
|
||||
#include "boost/program_options.hpp"
|
||||
#include "boost/test/included/unit_test.hpp"
|
||||
#include "boost/mpl/map.hpp"
|
||||
#include "boost/mpl/if.hpp"
|
||||
#include "bilinear_form_estimator_dev.h"
|
||||
#include "log_determinant_dev.h"
|
||||
|
||||
#ifdef EPETRA_MPI
|
||||
#include "trilinos/Epetra_MpiComm.h"
|
||||
#else
|
||||
#include "trilinos/Epetra_SerialComm.h"
|
||||
#endif
|
||||
|
||||
namespace fl {
|
||||
namespace ml {
|
||||
namespace bilinear_form_test {
|
||||
class BilinearFormTestSuite : public boost::unit_test_framework::test_suite {
|
||||
|
||||
public:
|
||||
|
||||
class BilinearFormTest {
|
||||
public:
|
||||
|
||||
BilinearFormTest() {
|
||||
}
|
||||
|
||||
void RunTests() {
|
||||
|
||||
fprintf(stderr, "Running the tests:\n");
|
||||
|
||||
// Call MPI Finalize.
|
||||
MPI_Finalize();
|
||||
}
|
||||
};
|
||||
|
||||
public:
|
||||
|
||||
BilinearFormTestSuite()
|
||||
: boost::unit_test_framework::test_suite("Bilinear form test suite") {
|
||||
|
||||
// Create an instance of test.
|
||||
boost::shared_ptr<BilinearFormTest> instance(new BilinearFormTest());
|
||||
|
||||
// Create the test cases.
|
||||
boost::unit_test_framework::test_case* bilinear_form_test_case
|
||||
= BOOST_CLASS_TEST_CASE(
|
||||
&BilinearFormTest::RunTests, instance);
|
||||
// add the test cases to the test suite
|
||||
add(bilinear_form_test_case);
|
||||
}
|
||||
};
|
||||
};
|
||||
};
|
||||
};
|
||||
|
||||
boost::unit_test_framework::test_suite*
|
||||
init_unit_test_suite(int argc, char** argv) {
|
||||
|
||||
// Initialize MPI.
|
||||
#ifdef EPETRA_MPI
|
||||
MPI_Init(&argc, &argv);
|
||||
#endif
|
||||
|
||||
// create the top test suite
|
||||
boost::unit_test_framework::test_suite* top_test_suite
|
||||
= BOOST_TEST_SUITE("Bilinear form tests");
|
||||
|
||||
if (argc != 2) {
|
||||
NOTIFY("Wrong number of arguments for tree test. Expected test input files directory. Returning NULL.");
|
||||
return NULL;
|
||||
}
|
||||
|
||||
// add test suites to the top test suite
|
||||
std::string input_files_directory = argv[1];
|
||||
top_test_suite->add(new fl::ml::bilinear_form_test::BilinearFormTestSuite());
|
||||
return top_test_suite;
|
||||
}
|
||||
@@ -0,0 +1,166 @@
|
||||
#ifndef FASTLIB_CONTRIB_DONGRYEL_TRILINOS_WRAPPERS_KERNEL_LINEAR_OPERATOR_H
|
||||
#define FASTLIB_CONTRIB_DONGRYEL_TRILINOS_WRAPPERS_KERNEL_LINEAR_OPERATOR_H
|
||||
|
||||
#include "fastlib/math/fl_math.h"
|
||||
#include "fastlib/trilinos_wrappers/linear_operator.h"
|
||||
|
||||
namespace Anasazi {
|
||||
|
||||
template<bool dotproduct_selfcase_special>
|
||||
class DotProductTrait {
|
||||
public:
|
||||
template<typename KernelType, typename PointType>
|
||||
DotProductTrait(const KernelType *kernel,
|
||||
const PointType &first_point,
|
||||
const PointType &second_point,
|
||||
bool point_indices_are_same,
|
||||
double *dotproduct);
|
||||
};
|
||||
|
||||
template<>
|
||||
class DotProductTrait<true> {
|
||||
public:
|
||||
template<typename KernelType, typename PointType>
|
||||
DotProductTrait(const KernelType *kernel,
|
||||
const PointType &first_point,
|
||||
const PointType &second_point,
|
||||
bool point_indices_are_same,
|
||||
double *dotproduct) {
|
||||
*dotproduct = kernel->Dot(first_point, second_point,
|
||||
point_indices_are_same);
|
||||
}
|
||||
};
|
||||
|
||||
template<>
|
||||
class DotProductTrait<false> {
|
||||
public:
|
||||
template<typename KernelType, typename PointType>
|
||||
DotProductTrait(const KernelType *kernel,
|
||||
const PointType &first_point,
|
||||
const PointType &second_point,
|
||||
bool point_indices_are_same,
|
||||
double *dotproduct) {
|
||||
*dotproduct = kernel->Dot(first_point, second_point);
|
||||
}
|
||||
};
|
||||
|
||||
template < typename TableType, typename KernelType, bool do_centering,
|
||||
bool dotproduct_selfcase_special = false >
|
||||
class KernelLinearOperator: public LinearOperator {
|
||||
|
||||
private:
|
||||
|
||||
TableType *table_;
|
||||
|
||||
const KernelType *kernel_;
|
||||
|
||||
fl::data::MonolithicPoint<double> average_row_;
|
||||
|
||||
double average_;
|
||||
|
||||
public:
|
||||
|
||||
KernelLinearOperator(TableType &table_in,
|
||||
const KernelType &kernel_in,
|
||||
#ifdef EPETRA_MPI
|
||||
const Epetra_MpiComm &comm_in,
|
||||
#else
|
||||
const Epetra_SerialComm &comm_in,
|
||||
#endif
|
||||
const Epetra_Map &map_in) {
|
||||
|
||||
table_ = &table_in;
|
||||
kernel_ = &kernel_in;
|
||||
comm_ = &comm_in;
|
||||
map_ = &map_in;
|
||||
|
||||
if (do_centering) {
|
||||
average_row_.Init(table_in.n_entries());
|
||||
}
|
||||
average_ = 0;
|
||||
|
||||
if (do_centering) {
|
||||
|
||||
// Precompute the average. This is a naive way of computing it.
|
||||
for (int i = 0; i < table_in.n_entries(); i++) {
|
||||
double average_for_i_th_point = 0;
|
||||
for (int j = 0; j < table_in.n_entries(); j++) {
|
||||
average_for_i_th_point += kernel_value(i, j);
|
||||
}
|
||||
average_for_i_th_point /= ((double) table_in.n_entries());
|
||||
average_row_[i] = average_for_i_th_point;
|
||||
}
|
||||
for (int i = 0; i < table_in.n_entries(); i++) {
|
||||
average_ += average_row_[i];
|
||||
}
|
||||
average_ /= ((double) table_in.n_entries());
|
||||
}
|
||||
}
|
||||
|
||||
double centered_kernel_value(int row, int col) const {
|
||||
typename TableType::Dataset_t::Point_t row_point;
|
||||
typename TableType::Dataset_t::Point_t col_point;
|
||||
table_->get(row, &row_point);
|
||||
table_->get(col, &col_point);
|
||||
double dotproduct = 0;
|
||||
DotProductTrait<dotproduct_selfcase_special>(
|
||||
kernel_, row_point, col_point, row == col, &dotproduct);
|
||||
return dotproduct - average_row_[row] - average_row_[col] + average_;
|
||||
}
|
||||
|
||||
double kernel_value(int row, int col) const {
|
||||
typename TableType::Dataset_t::Point_t row_point;
|
||||
typename TableType::Dataset_t::Point_t col_point;
|
||||
table_->get(row, &row_point);
|
||||
table_->get(col, &col_point);
|
||||
double dotproduct = 0;
|
||||
DotProductTrait<dotproduct_selfcase_special>(
|
||||
kernel_, row_point, col_point, row == col, &dotproduct);
|
||||
return dotproduct;
|
||||
}
|
||||
|
||||
int Apply(
|
||||
const Epetra_MultiVector &vecs,
|
||||
Epetra_MultiVector &prods) const {
|
||||
|
||||
prods.PutScalar(0);
|
||||
|
||||
for (int j = 0; j < table_->n_entries(); j++) {
|
||||
for (int i = 0; i < table_->n_entries(); i++) {
|
||||
double pair_kernel_value = (do_centering) ?
|
||||
centered_kernel_value(i, j) : kernel_value(i, j);
|
||||
for (int k = 0; k < vecs.NumVectors(); k++) {
|
||||
prods.Pointers()[k][i] += pair_kernel_value * vecs.Pointers()[k][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
int n_rows() const {
|
||||
return table_->n_entries();
|
||||
}
|
||||
|
||||
int n_cols() const {
|
||||
return table_->n_entries();
|
||||
}
|
||||
|
||||
double get(int row, int col) const {
|
||||
return (do_centering) ? centered_kernel_value(row, col) : kernel_value(row, col);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename TableType, typename KernelType, bool do_centering>
|
||||
class OperatorTraits < double, Epetra_MultiVector,
|
||||
KernelLinearOperator<TableType, KernelType, do_centering> > {
|
||||
public:
|
||||
|
||||
static void Apply(const Epetra_Operator& Op,
|
||||
const Epetra_MultiVector& x,
|
||||
Epetra_MultiVector& y) {
|
||||
Op.Apply(x, y);
|
||||
}
|
||||
};
|
||||
};
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,129 @@
|
||||
|
||||
#ifndef FASTLIB_CONTRIB_DONGRYEL_FASTLIB_TRILINOS_WRAPPERS_LINEAR_OPERATOR_H
|
||||
#define FASTLIB_CONTRIB_DONGRYEL_FASTLIB_TRILINOS_WRAPPERS_LINEAR_OPERATOR_H
|
||||
|
||||
#undef F77_FUNC
|
||||
#undef LI
|
||||
#include "trilinos/AnasaziEpetraAdapter.hpp"
|
||||
#include "trilinos/AnasaziBasicEigenproblem.hpp"
|
||||
#include "trilinos/AnasaziBlockKrylovSchurSolMgr.hpp"
|
||||
#include "trilinos/AnasaziBasicSort.hpp"
|
||||
#include "trilinos/AztecOO.h"
|
||||
#include "trilinos/Epetra_BlockMap.h"
|
||||
#include "trilinos/Epetra_CrsMatrix.h"
|
||||
#include "trilinos/Epetra_DataAccess.h"
|
||||
#include "trilinos/Epetra_LinearProblem.h"
|
||||
#include "trilinos/Epetra_Map.h"
|
||||
#include "trilinos/Epetra_MultiVector.h"
|
||||
#include "trilinos/Epetra_Operator.h"
|
||||
|
||||
#ifdef EPETRA_MPI
|
||||
#include "trilinos/Epetra_MpiComm.h"
|
||||
#else
|
||||
#include "trilinos/Epetra_SerialComm.h"
|
||||
#endif
|
||||
|
||||
#include "trilinos/Epetra_Vector.h"
|
||||
#undef F77_FUNC
|
||||
#include "fastlib/la/matrix.h"
|
||||
#include <vector>
|
||||
|
||||
namespace Anasazi {
|
||||
|
||||
class LinearOperator: public virtual Epetra_Operator {
|
||||
|
||||
protected:
|
||||
|
||||
#ifdef EPETRA_MPI
|
||||
const Epetra_MpiComm *comm_;
|
||||
#else
|
||||
const Epetra_SerialComm *comm_;
|
||||
#endif
|
||||
|
||||
const Epetra_Map *map_;
|
||||
|
||||
public:
|
||||
|
||||
virtual ~LinearOperator() {
|
||||
}
|
||||
|
||||
LinearOperator() {
|
||||
comm_ = NULL;
|
||||
map_ = NULL;
|
||||
}
|
||||
|
||||
#ifdef EPETRA_MPI
|
||||
LinearOperator(const Epetra_MpiComm &comm_in,
|
||||
const Epetra_Map &map_in) {
|
||||
comm_ = &comm_in;
|
||||
map_ = &map_in;
|
||||
}
|
||||
#else
|
||||
LinearOperator(const Epetra_SerialComm &comm_in,
|
||||
const Epetra_Map &map_in) {
|
||||
comm_ = &comm_in;
|
||||
map_ = &map_in;
|
||||
}
|
||||
#endif
|
||||
|
||||
virtual int Apply(const Epetra_MultiVector &vec,
|
||||
Epetra_MultiVector &prod) const = 0;
|
||||
|
||||
int SetUseTranspose(bool use_transpose) {
|
||||
return -1;
|
||||
}
|
||||
|
||||
int ApplyInverse(const Epetra_MultiVector &X,
|
||||
Epetra_MultiVector &Y) const {
|
||||
return -1;
|
||||
}
|
||||
|
||||
double NormInf() const {
|
||||
return -1;
|
||||
}
|
||||
|
||||
const char *Label() const {
|
||||
return "Generic linear operator";
|
||||
}
|
||||
|
||||
bool UseTranspose() const {
|
||||
return false;
|
||||
}
|
||||
|
||||
bool HasNormInf() const {
|
||||
return false;
|
||||
}
|
||||
|
||||
const Epetra_Comm &Comm() const {
|
||||
return *comm_;
|
||||
}
|
||||
|
||||
const Epetra_Map &OperatorDomainMap() const {
|
||||
const Epetra_Map &map_reference = *map_;
|
||||
return map_reference;
|
||||
}
|
||||
|
||||
const Epetra_Map &OperatorRangeMap() const {
|
||||
const Epetra_Map &map_reference = *map_;
|
||||
return map_reference;
|
||||
}
|
||||
|
||||
void PrintDebug(const char *name = "", FILE *stream = stderr) const {
|
||||
fprintf(stream, "----- MATRIX ------: %s\n", name);
|
||||
for (int r = 0; r < this->n_rows(); r++) {
|
||||
for (int c = 0; c < this->n_cols(); c++) {
|
||||
fprintf(stream, "%+3.3f ", this->get(r, c));
|
||||
}
|
||||
fprintf(stream, "\n");
|
||||
}
|
||||
}
|
||||
|
||||
virtual int n_rows() const = 0;
|
||||
|
||||
virtual int n_cols() const = 0;
|
||||
|
||||
virtual double get(int row, int col) const = 0;
|
||||
};
|
||||
};
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,39 @@
|
||||
/** @author Dongryeol Lee
|
||||
*
|
||||
* @file log_determinant.h
|
||||
*/
|
||||
|
||||
#ifndef FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_LOG_DETERMINANT_H
|
||||
#define FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_LOG_DETERMINANT_H
|
||||
|
||||
#include "bilinear_form_estimator.h"
|
||||
#include "fastlib/la/matrix.h"
|
||||
|
||||
namespace fl {
|
||||
namespace ml {
|
||||
class LogDeterminant {
|
||||
|
||||
private:
|
||||
|
||||
fl::ml::BilinearFormEstimator<fl::ml::LogTransformation> bilinear_log_form_;
|
||||
|
||||
private:
|
||||
|
||||
void RandomVector_(GenVector<double> &v);
|
||||
|
||||
public:
|
||||
|
||||
LogDeterminant();
|
||||
|
||||
void Init(Anasazi::LinearOperator *linear_operator_in);
|
||||
|
||||
double MonteCarloCompute();
|
||||
|
||||
double Compute();
|
||||
|
||||
double NaiveCompute();
|
||||
};
|
||||
};
|
||||
};
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,85 @@
|
||||
/** @author Dongryeol Lee
|
||||
*
|
||||
* @file log_determinant_dev.h
|
||||
*/
|
||||
|
||||
#ifndef FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_LOG_DETERMINANT_DEV_H
|
||||
#define FASTLIB_CONTRIB_DONGRYEL_GP_REGRESSION_LOG_DETERMINANT_DEV_H
|
||||
|
||||
#include "fastlib/la/matrix.h"
|
||||
#include "bilinear_form_estimator_dev.h"
|
||||
#include "log_determinant.h"
|
||||
|
||||
namespace fl {
|
||||
namespace ml {
|
||||
LogDeterminant::LogDeterminant() {
|
||||
}
|
||||
|
||||
void LogDeterminant::Init(Anasazi::LinearOperator *linear_operator_in) {
|
||||
bilinear_log_form_.Init(linear_operator_in);
|
||||
}
|
||||
|
||||
void LogDeterminant::RandomVector_(GenVector<double> &v) {
|
||||
for (int i = 0; i < v.length(); i++) {
|
||||
v[i] = (math::Random() >= 0.5) ? 1 : -1;
|
||||
}
|
||||
}
|
||||
|
||||
double LogDeterminant::MonteCarloCompute() {
|
||||
|
||||
// A random vector for the samples.
|
||||
GenVector<double> random_vector;
|
||||
random_vector.Init(bilinear_log_form_.linear_operator()->n_rows());
|
||||
|
||||
double log_determinant = 0;
|
||||
const int num_samples = 100;
|
||||
|
||||
for (int i = 0; i < num_samples; i++) {
|
||||
RandomVector_(random_vector);
|
||||
log_determinant += bilinear_log_form_.Compute(random_vector);
|
||||
}
|
||||
log_determinant /= ((double) num_samples);
|
||||
return log_determinant;
|
||||
}
|
||||
|
||||
double LogDeterminant::NaiveCompute() {
|
||||
|
||||
double log_determinant = 0;
|
||||
GenMatrix<double> kernel_matrix;
|
||||
kernel_matrix.Init(bilinear_log_form_.linear_operator()->n_rows(),
|
||||
bilinear_log_form_.linear_operator()->n_cols());
|
||||
for (int j = 0; j < bilinear_log_form_.linear_operator()->n_cols(); j++) {
|
||||
for (int i = 0; i < bilinear_log_form_.linear_operator()->n_rows(); i++) {
|
||||
kernel_matrix.set(i, j, bilinear_log_form_.linear_operator()->get(i, j));
|
||||
}
|
||||
}
|
||||
GenVector<double> eigenvalues;
|
||||
la::SVDInit(kernel_matrix, &eigenvalues);
|
||||
|
||||
for (int i = 0; i < eigenvalues.length(); i++) {
|
||||
log_determinant += log(eigenvalues[i]);
|
||||
}
|
||||
return log_determinant;
|
||||
}
|
||||
|
||||
double LogDeterminant::Compute() {
|
||||
|
||||
// Do a naive for-loop over each row and apply $e_i^T log(A) e_i$.
|
||||
GenVector<double> i_th_unit_vector;
|
||||
i_th_unit_vector.Init(bilinear_log_form_.linear_operator()->n_rows());
|
||||
i_th_unit_vector.SetZero();
|
||||
double log_determinant = 0;
|
||||
|
||||
for (int i = 0; i < bilinear_log_form_.linear_operator()->n_rows(); i++) {
|
||||
i_th_unit_vector[i] = 1;
|
||||
if (i > 0) {
|
||||
i_th_unit_vector[i - 1] = 0;
|
||||
}
|
||||
log_determinant += bilinear_log_form_.Compute(i_th_unit_vector);
|
||||
}
|
||||
return log_determinant;
|
||||
}
|
||||
};
|
||||
};
|
||||
|
||||
#endif
|
||||
Reference in New Issue
Block a user