Commit 02d18154 authored by Jack Poulson's avatar Jack Poulson

Adding support for DPP log-likelihood computation; Adding more timing data for...

Adding support for DPP log-likelihood computation; Adding more timing data for examples and writing TIFF for negative 2D Laplacian DPP.
parent fc3a3605
......@@ -133,6 +133,28 @@ void RunLDLTransposeFactorization(Int tile_size, Int block_size,
}
}
template <typename Field>
void TestFactorizations(Int matrix_size, Int num_rounds, Int tile_size,
Int block_size) {
BlasMatrix<Field> matrix;
Buffer<Field> extra_buffer;
for (Int round = 0; round < num_rounds; ++round) {
InitializeMatrix(matrix_size, &matrix);
RunCholeskyFactorization(tile_size, block_size, &matrix.view);
InitializeMatrix(matrix_size, &matrix);
RunLDLAdjointFactorization(tile_size, block_size, &matrix.view,
&extra_buffer);
InitializeMatrix(matrix_size, &matrix);
RunLDLTransposeFactorization(tile_size, block_size, &matrix.view,
&extra_buffer);
std::cout << std::endl;
}
}
} // anonymous namespace
int main(int argc, char** argv) {
......@@ -149,23 +171,22 @@ int main(int argc, char** argv) {
return 0;
}
BlasMatrix<Complex<double>> matrix;
Buffer<Complex<double>> extra_buffer;
std::cout << "Testing for float:" << std::endl;
TestFactorizations<float>(matrix_size, num_rounds, tile_size, block_size);
std::cout << std::endl;
for (Int round = 0; round < num_rounds; ++round) {
InitializeMatrix(matrix_size, &matrix);
RunCholeskyFactorization(tile_size, block_size, &matrix.view);
std::cout << "Testing for double:" << std::endl;
TestFactorizations<double>(matrix_size, num_rounds, tile_size, block_size);
std::cout << std::endl;
InitializeMatrix(matrix_size, &matrix);
RunLDLAdjointFactorization(tile_size, block_size, &matrix.view,
&extra_buffer);
std::cout << "Testing for Complex<float>:" << std::endl;
TestFactorizations<Complex<float>>(matrix_size, num_rounds, tile_size,
block_size);
std::cout << std::endl;
InitializeMatrix(matrix_size, &matrix);
RunLDLTransposeFactorization(tile_size, block_size, &matrix.view,
&extra_buffer);
std::cout << std::endl;
}
std::cout << "Testing for Complex<double>:" << std::endl;
TestFactorizations<Complex<double>>(matrix_size, num_rounds, tile_size,
block_size);
return 0;
}
......@@ -15,6 +15,51 @@
using catamari::Int;
#ifdef CATAMARI_HAVE_LIBTIFF
// Returns whether the given value exists within a sorted list.
template <typename Iterator, typename T>
bool IndexExists(Iterator beg, Iterator end, T value) {
auto iter = std::lower_bound(beg, end, value);
return iter != end && *iter == value;
}
// Writes out the sample to a TIFF file.
inline void WriteGridToTIFF(const std::string& filename, Int x_size, Int y_size,
const std::vector<Int>& sample, const Int box_size,
const catamari::Pixel& background_pixel,
const catamari::Pixel& active_pixel) {
const Int height = box_size * y_size;
const Int width = box_size * x_size;
const Int samples_per_pixel = 3;
std::vector<char> image(width * height * samples_per_pixel);
// Fill the background pixels.
for (Int i = 0; i < height; ++i) {
for (Int j = 0; j < width; ++j) {
catamari::PackPixel(i, j, width, background_pixel, &image);
}
}
for (Int y = 0; y < y_size; ++y) {
for (Int x = 0; x < x_size; ++x) {
const Int index = y + x * y_size;
const bool found = IndexExists(sample.begin(), sample.end(), index);
if (found) {
const Int i_offset = y * box_size;
const Int j_offset = x * box_size;
for (Int i = i_offset; i < i_offset + box_size; ++i) {
for (Int j = j_offset; j < j_offset + box_size; ++j) {
catamari::PackPixel(i, j, width, active_pixel, &image);
}
}
}
}
}
catamari::WriteTIFF(filename, height, width, samples_per_pixel, image);
}
#endif // ifdef CATAMARI_HAVE_LIBTIFF
// A list of properties to measure from DPP sampling.
struct Experiment {
// The number of seconds that elapsed during the object construction.
......@@ -120,10 +165,18 @@ inline void AsciiDisplaySample(Int x_size, Int y_size,
Experiment RunShifted2DNegativeLaplacianTest(
Int x_size, Int y_size, double diagonal_shift, double scale,
bool maximum_likelihood, bool analytical_ordering, bool print_sample,
bool ascii_display, char missing_char, char sampled_char, Int num_samples,
bool ascii_display, char missing_char, char sampled_char,
bool CATAMARI_UNUSED write_tiff, Int num_samples,
const catamari::DPPControl& dpp_control, bool print_progress) {
Experiment experiment;
#ifdef CATAMARI_HAVE_LIBTIFF
// TIFF configuration.
const catamari::Pixel background_pixel{char(255), char(255), char(255)};
const catamari::Pixel active_pixel{char(255), char(0), char(0)};
const Int box_size = 5;
#endif // ifdef CATAMARI_HAVE_LIBTIFF
// Read the matrix from file.
std::unique_ptr<catamari::CoordinateMatrix<double>> matrix =
Shifted2DNegativeLaplacian(x_size, y_size, diagonal_shift, scale,
......@@ -156,9 +209,17 @@ Experiment RunShifted2DNegativeLaplacianTest(
sample = dpp->Sample(maximum_likelihood);
const double sample_seconds = sample_timer.Stop();
std::cout << " sample took " << sample_seconds << " seconds." << std::endl;
const double log_likelihood = dpp->LogLikelihood();
std::cout << " log-likelihood was: " << log_likelihood << std::endl;
if (print_sample) {
quotient::PrintVector(sample, "sample", std::cout);
}
if (write_tiff) {
const std::string filename =
std::string("shifted_laplacian-") + std::to_string(run) + ".tif";
WriteGridToTIFF(filename, x_size, y_size, sample, box_size,
background_pixel, active_pixel);
}
if (ascii_display) {
AsciiDisplaySample(x_size, y_size, sample, missing_char, sampled_char);
}
......@@ -186,7 +247,9 @@ int main(int argc, char** argv) {
const bool print_sample = parser.OptionalInput<bool>(
"print_sample", "Print each DPP sample?", false);
const bool ascii_display = parser.OptionalInput<bool>(
"ascii_display", "Display sample in ASCII?", true);
"ascii_display", "Display sample in ASCII?", false);
const bool write_tiff =
parser.OptionalInput<bool>("write_tiff", "Write samples to TIFF?", true);
const char missing_char = parser.OptionalInput<char>(
"missing_char", "ASCII display for missing index.", ' ');
const char sampled_char = parser.OptionalInput<char>(
......@@ -222,7 +285,7 @@ int main(int argc, char** argv) {
const Experiment experiment = RunShifted2DNegativeLaplacianTest(
x_size, y_size, diagonal_shift, scale, maximum_likelihood,
analytical_ordering, print_sample, ascii_display, missing_char,
sampled_char, num_samples, dpp_control, print_progress);
sampled_char, write_tiff, num_samples, dpp_control, print_progress);
PrintExperiment(experiment);
return 0;
......
......@@ -106,6 +106,15 @@ std::vector<Int> DPP<Field>::Sample(bool maximum_likelihood) const {
}
}
template <class Field>
ComplexBase<Field> DPP<Field>::LogLikelihood() const {
if (is_supernodal_) {
return supernodal_dpp_->LogLikelihood();
} else {
return scalar_dpp_->LogLikelihood();
}
}
} // namespace catamari
#endif // ifndef CATAMARI_DPP_IMPL_H_
......@@ -44,6 +44,9 @@ class DPP {
// pivot is kept based upon which choice is most likely.
std::vector<Int> Sample(bool maximum_likelihood) const;
// Returns the log-likelihood of the last DPP sample.
ComplexBase<Field> LogLikelihood() const;
private:
// Whether or not a supernodal sampler was used. If it is true, only
// 'supernodal_dpp_' should be non-null, and vice versa.
......
......@@ -155,6 +155,19 @@ std::vector<Int> ScalarDPP<Field>::Sample(bool maximum_likelihood) const {
return UpLookingSample(maximum_likelihood);
}
template <typename Field>
ComplexBase<Field> ScalarDPP<Field>::LogLikelihood() const {
typedef ComplexBase<Field> Real;
const Int num_values = diagonal_factor_->values.Size();
Real log_likelihood = 0;
for (Int j = 0; j < num_values; ++j) {
log_likelihood += std::log(std::abs(diagonal_factor_->values[j]));
}
return log_likelihood;
}
} // namespace catamari
#endif // ifndef CATAMARI_SCALAR_DPP_IMPL_H_
......@@ -31,6 +31,9 @@ class ScalarDPP {
// pivot is kept based upon which choice is most likely.
std::vector<Int> Sample(bool maximum_likelihood) const;
// Returns the log-likelihood of the last sample.
ComplexBase<Field> LogLikelihood() const;
private:
typedef ComplexBase<Field> Real;
......
......@@ -79,6 +79,9 @@ class SupernodalDPP {
// pivot is kept based upon which choice is most likely.
std::vector<Int> Sample(bool maximum_likelihood) const;
// Returns the log-likelihood of the last sample.
ComplexBase<Field> LogLikelihood() const;
private:
typedef ComplexBase<Field> Real;
......
......@@ -162,6 +162,24 @@ void SupernodalDPP<Field>::AppendSupernodeSample(
}
}
template <typename Field>
ComplexBase<Field> SupernodalDPP<Field>::LogLikelihood() const {
typedef ComplexBase<Field> Real;
const Int num_supernodes = diagonal_factor_->blocks.Size();
Real log_likelihood = 0;
for (Int supernode = 0; supernode < num_supernodes; ++supernode) {
const ConstBlasMatrixView<Field> diag_block =
diagonal_factor_->blocks[supernode];
const Int supernode_size = diag_block.height;
for (Int j = 0; j < supernode_size; ++j) {
log_likelihood += std::log(std::abs(diag_block(j, j)));
}
}
return log_likelihood;
}
} // namespace catamari
#endif // ifndef CATAMARI_SUPERNODAL_DPP_COMMON_IMPL_H_
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment