blob: b18a1c6a1e35a5294b3545bdb835ae3ad4fbc3b4 [file] [log] [blame]
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2023 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: vitus@google.com (Michael Vitus)
#include "ceres/parallel_for.h"
#include <cmath>
#include <condition_variable>
#include <mutex>
#include <numeric>
#include <random>
#include <thread>
#include <tuple>
#include <vector>
#include "ceres/context_impl.h"
#include "ceres/internal/config.h"
#include "ceres/parallel_vector_ops.h"
#include "glog/logging.h"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
namespace ceres::internal {
using testing::ElementsAreArray;
using testing::UnorderedElementsAreArray;
// Tests the parallel for loop computes the correct result for various number of
// threads.
TEST(ParallelFor, NumThreads) {
ContextImpl context;
context.EnsureMinimumThreads(/*num_threads=*/2);
const int size = 16;
std::vector<int> expected_results(size, 0);
for (int i = 0; i < size; ++i) {
expected_results[i] = std::sqrt(i);
}
for (int num_threads = 1; num_threads <= 8; ++num_threads) {
std::vector<int> values(size, 0);
ParallelFor(&context, 0, size, num_threads, [&values](int i) {
values[i] = std::sqrt(i);
});
EXPECT_THAT(values, ElementsAreArray(expected_results));
}
}
// Tests parallel for loop with ranges
TEST(ParallelForWithRange, NumThreads) {
ContextImpl context;
context.EnsureMinimumThreads(/*num_threads=*/2);
const int size = 16;
std::vector<int> expected_results(size, 0);
for (int i = 0; i < size; ++i) {
expected_results[i] = std::sqrt(i);
}
for (int num_threads = 1; num_threads <= 8; ++num_threads) {
std::vector<int> values(size, 0);
ParallelFor(
&context, 0, size, num_threads, [&values](std::tuple<int, int> range) {
auto [start, end] = range;
for (int i = start; i < end; ++i) values[i] = std::sqrt(i);
});
EXPECT_THAT(values, ElementsAreArray(expected_results));
}
}
// Tests the parallel for loop with the thread ID interface computes the correct
// result for various number of threads.
TEST(ParallelForWithThreadId, NumThreads) {
ContextImpl context;
context.EnsureMinimumThreads(/*num_threads=*/2);
const int size = 16;
std::vector<int> expected_results(size, 0);
for (int i = 0; i < size; ++i) {
expected_results[i] = std::sqrt(i);
}
for (int num_threads = 1; num_threads <= 8; ++num_threads) {
std::vector<int> values(size, 0);
ParallelFor(
&context, 0, size, num_threads, [&values](int thread_id, int i) {
values[i] = std::sqrt(i);
});
EXPECT_THAT(values, ElementsAreArray(expected_results));
}
}
// Tests nested for loops do not result in a deadlock.
TEST(ParallelFor, NestedParallelForDeadlock) {
ContextImpl context;
context.EnsureMinimumThreads(/*num_threads=*/2);
// Increment each element in the 2D matrix.
std::vector<std::vector<int>> x(3, {1, 2, 3});
ParallelFor(&context, 0, 3, 2, [&x, &context](int i) {
std::vector<int>& y = x.at(i);
ParallelFor(&context, 0, 3, 2, [&y](int j) { ++y.at(j); });
});
const std::vector<int> results = {2, 3, 4};
for (const std::vector<int>& value : x) {
EXPECT_THAT(value, ElementsAreArray(results));
}
}
// Tests nested for loops do not result in a deadlock for the parallel for with
// thread ID interface.
TEST(ParallelForWithThreadId, NestedParallelForDeadlock) {
ContextImpl context;
context.EnsureMinimumThreads(/*num_threads=*/2);
// Increment each element in the 2D matrix.
std::vector<std::vector<int>> x(3, {1, 2, 3});
ParallelFor(&context, 0, 3, 2, [&x, &context](int thread_id, int i) {
std::vector<int>& y = x.at(i);
ParallelFor(&context, 0, 3, 2, [&y](int thread_id, int j) { ++y.at(j); });
});
const std::vector<int> results = {2, 3, 4};
for (const std::vector<int>& value : x) {
EXPECT_THAT(value, ElementsAreArray(results));
}
}
TEST(ParallelForWithThreadId, UniqueThreadIds) {
// Ensure the hardware supports more than 1 thread to ensure the test will
// pass.
const int num_hardware_threads = std::thread::hardware_concurrency();
if (num_hardware_threads <= 1) {
LOG(ERROR)
<< "Test not supported, the hardware does not support threading.";
return;
}
ContextImpl context;
context.EnsureMinimumThreads(/*num_threads=*/2);
// Increment each element in the 2D matrix.
std::vector<int> x(2, -1);
std::mutex mutex;
std::condition_variable condition;
int count = 0;
ParallelFor(&context,
0,
2,
2,
[&x, &mutex, &condition, &count](int thread_id, int i) {
std::unique_lock<std::mutex> lock(mutex);
x[i] = thread_id;
++count;
condition.notify_all();
condition.wait(lock, [&]() { return count == 2; });
});
EXPECT_THAT(x, UnorderedElementsAreArray({0, 1}));
}
// Helper function for partition tests
bool BruteForcePartition(
int* costs, int start, int end, int max_partitions, int max_cost);
// Basic test if MaxPartitionCostIsFeasible and BruteForcePartition agree on
// simple test-cases
TEST(GuidedParallelFor, MaxPartitionCostIsFeasible) {
std::vector<int> costs, cumulative_costs, partition;
costs = {1, 2, 3, 5, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0};
cumulative_costs.resize(costs.size());
std::partial_sum(costs.begin(), costs.end(), cumulative_costs.begin());
const auto dummy_getter = [](const int v) { return v; };
// [1, 2, 3] [5], [0 ... 0, 7, 0, ... 0]
EXPECT_TRUE(MaxPartitionCostIsFeasible(0,
costs.size(),
3,
7,
0,
cumulative_costs.data(),
dummy_getter,
&partition));
EXPECT_TRUE(BruteForcePartition(costs.data(), 0, costs.size(), 3, 7));
// [1, 2, 3, 5, 0 ... 0, 7, 0, ... 0]
EXPECT_TRUE(MaxPartitionCostIsFeasible(0,
costs.size(),
3,
18,
0,
cumulative_costs.data(),
dummy_getter,
&partition));
EXPECT_TRUE(BruteForcePartition(costs.data(), 0, costs.size(), 3, 18));
// Impossible since there is item of cost 7
EXPECT_FALSE(MaxPartitionCostIsFeasible(0,
costs.size(),
3,
6,
0,
cumulative_costs.data(),
dummy_getter,
&partition));
EXPECT_FALSE(BruteForcePartition(costs.data(), 0, costs.size(), 3, 6));
// Impossible
EXPECT_FALSE(MaxPartitionCostIsFeasible(0,
costs.size(),
2,
10,
0,
cumulative_costs.data(),
dummy_getter,
&partition));
EXPECT_FALSE(BruteForcePartition(costs.data(), 0, costs.size(), 2, 10));
}
// Randomized tests for MaxPartitionCostIsFeasible
TEST(GuidedParallelFor, MaxPartitionCostIsFeasibleRandomized) {
std::vector<int> costs, cumulative_costs, partition;
const auto dummy_getter = [](const int v) { return v; };
// Random tests
const int kNumTests = 1000;
const int kMaxElements = 32;
const int kMaxPartitions = 16;
const int kMaxElCost = 8;
std::mt19937 rng;
std::uniform_int_distribution<int> rng_N(1, kMaxElements);
std::uniform_int_distribution<int> rng_M(1, kMaxPartitions);
std::uniform_int_distribution<int> rng_e(0, kMaxElCost);
for (int t = 0; t < kNumTests; ++t) {
const int N = rng_N(rng);
const int M = rng_M(rng);
int total = 0;
costs.clear();
for (int i = 0; i < N; ++i) {
costs.push_back(rng_e(rng));
total += costs.back();
}
cumulative_costs.resize(N);
std::partial_sum(costs.begin(), costs.end(), cumulative_costs.begin());
std::uniform_int_distribution<int> rng_seg(0, N - 1);
int start = rng_seg(rng);
int end = rng_seg(rng);
if (start > end) std::swap(start, end);
++end;
int first_admissible = 0;
for (int threshold = 1; threshold <= total; ++threshold) {
const bool bruteforce =
BruteForcePartition(costs.data(), start, end, M, threshold);
if (bruteforce && !first_admissible) {
first_admissible = threshold;
}
const bool binary_search =
MaxPartitionCostIsFeasible(start,
end,
M,
threshold,
start ? cumulative_costs[start - 1] : 0,
cumulative_costs.data(),
dummy_getter,
&partition);
EXPECT_EQ(bruteforce, binary_search);
EXPECT_LE(partition.size(), M + 1);
// check partition itself
if (binary_search) {
ASSERT_GT(partition.size(), 1);
EXPECT_EQ(partition.front(), start);
EXPECT_EQ(partition.back(), end);
const int num_partitions = partition.size() - 1;
EXPECT_LE(num_partitions, M);
for (int j = 0; j < num_partitions; ++j) {
int total = 0;
for (int k = partition[j]; k < partition[j + 1]; ++k) {
EXPECT_LT(k, end);
EXPECT_GE(k, start);
total += costs[k];
}
EXPECT_LE(total, threshold);
}
}
}
}
}
TEST(GuidedParallelFor, PartitionRangeForParallelFor) {
std::vector<int> costs, cumulative_costs, partition;
const auto dummy_getter = [](const int v) { return v; };
// Random tests
const int kNumTests = 1000;
const int kMaxElements = 32;
const int kMaxPartitions = 16;
const int kMaxElCost = 8;
std::mt19937 rng;
std::uniform_int_distribution<int> rng_N(1, kMaxElements);
std::uniform_int_distribution<int> rng_M(1, kMaxPartitions);
std::uniform_int_distribution<int> rng_e(0, kMaxElCost);
for (int t = 0; t < kNumTests; ++t) {
const int N = rng_N(rng);
const int M = rng_M(rng);
int total = 0;
costs.clear();
for (int i = 0; i < N; ++i) {
costs.push_back(rng_e(rng));
total += costs.back();
}
cumulative_costs.resize(N);
std::partial_sum(costs.begin(), costs.end(), cumulative_costs.begin());
std::uniform_int_distribution<int> rng_seg(0, N - 1);
int start = rng_seg(rng);
int end = rng_seg(rng);
if (start > end) std::swap(start, end);
++end;
int first_admissible = 0;
for (int threshold = 1; threshold <= total; ++threshold) {
const bool bruteforce =
BruteForcePartition(costs.data(), start, end, M, threshold);
if (bruteforce) {
first_admissible = threshold;
break;
}
}
EXPECT_TRUE(first_admissible != 0 || total == 0);
partition = PartitionRangeForParallelFor(
start, end, M, cumulative_costs.data(), dummy_getter);
ASSERT_GT(partition.size(), 1);
EXPECT_EQ(partition.front(), start);
EXPECT_EQ(partition.back(), end);
const int num_partitions = partition.size() - 1;
EXPECT_LE(num_partitions, M);
for (int j = 0; j < num_partitions; ++j) {
int total = 0;
for (int k = partition[j]; k < partition[j + 1]; ++k) {
EXPECT_LT(k, end);
EXPECT_GE(k, start);
total += costs[k];
}
EXPECT_LE(total, first_admissible);
}
}
}
// Recursively try to partition range into segements of total cost
// less than max_cost
bool BruteForcePartition(
int* costs, int start, int end, int max_partitions, int max_cost) {
if (start == end) return true;
if (start < end && max_partitions == 0) return false;
int total_cost = 0;
for (int last_curr = start + 1; last_curr <= end; ++last_curr) {
total_cost += costs[last_curr - 1];
if (total_cost > max_cost) break;
if (BruteForcePartition(
costs, last_curr, end, max_partitions - 1, max_cost))
return true;
}
return false;
}
// Tests if guided parallel for loop computes the correct result for various
// number of threads.
TEST(GuidedParallelFor, NumThreads) {
ContextImpl context;
context.EnsureMinimumThreads(/*num_threads=*/2);
const int size = 16;
std::vector<int> expected_results(size, 0);
for (int i = 0; i < size; ++i) {
expected_results[i] = std::sqrt(i);
}
std::vector<int> costs, cumulative_costs;
for (int i = 1; i <= size; ++i) {
int cost = i * i;
costs.push_back(cost);
if (i == 1) {
cumulative_costs.push_back(cost);
} else {
cumulative_costs.push_back(cost + cumulative_costs.back());
}
}
for (int num_threads = 1; num_threads <= 8; ++num_threads) {
std::vector<int> values(size, 0);
ParallelFor(
&context,
0,
size,
num_threads,
[&values](int i) { values[i] = std::sqrt(i); },
cumulative_costs.data(),
[](const int v) { return v; });
EXPECT_THAT(values, ElementsAreArray(expected_results));
}
}
TEST(ParallelAssign, D2MulX) {
const int kVectorSize = 1024 * 1024;
const int kMaxNumThreads = 8;
const double kEpsilon = 1e-16;
const Vector D_full = Vector::Random(kVectorSize * 2);
const ConstVectorRef D(D_full.data() + kVectorSize, kVectorSize);
const Vector x = Vector::Random(kVectorSize);
const Vector y_expected = D.array().square() * x.array();
ContextImpl context;
context.EnsureMinimumThreads(kMaxNumThreads);
for (int num_threads = 1; num_threads <= kMaxNumThreads; ++num_threads) {
Vector y_observed(kVectorSize);
ParallelAssign(
&context, num_threads, y_observed, D.array().square() * x.array());
// We might get non-bit-exact result due to different precision in scalar
// and vector code. For example, in x86 mode mingw might emit x87
// instructions for scalar code, thus making bit-exact check fail
EXPECT_NEAR((y_expected - y_observed).squaredNorm(),
0.,
kEpsilon * y_expected.squaredNorm());
}
}
TEST(ParallelAssign, SetZero) {
const int kVectorSize = 1024 * 1024;
const int kMaxNumThreads = 8;
ContextImpl context;
context.EnsureMinimumThreads(kMaxNumThreads);
for (int num_threads = 1; num_threads <= kMaxNumThreads; ++num_threads) {
Vector x = Vector::Random(kVectorSize);
ParallelSetZero(&context, num_threads, x);
CHECK_EQ(x.squaredNorm(), 0.);
}
}
} // namespace ceres::internal