Шилин Никита. Технология SEQ. Вычисление многомерных интегралов методом Монте-Карло. Вариант 12

tasks/shilin_n_monte_carlo_integration/common/include/common.hpp

@@ -0,0 +1,141 @@
+#pragma once
+
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <string>
+#include <tuple>
+#include <vector>
+
+#include "task/include/task.hpp"
+
+namespace shilin_n_monte_carlo_integration {
+
+enum class FuncType : std::uint8_t {
+  kConstant = 0,
+  kLinear = 1,
+  kProduct = 2,
+  kSumSquares = 3,
+  kSinProduct = 4,
+};
+
+class IntegrandFunction {
+ public:
+  static double Evaluate(FuncType func_type, const std::vector<double> &point) {
+    switch (func_type) {
+      case FuncType::kConstant:
+        return 1.0;
+      case FuncType::kLinear: {
+        double sum = 0.0;
+        for (double x : point) {
+          sum += x;
+        }
+        return sum;
+      }
+      case FuncType::kProduct: {
+        double prod = 1.0;
+        for (double x : point) {
+          prod *= x;
+        }
+        return prod;
+      }
+      case FuncType::kSumSquares: {
+        double sum = 0.0;
+        for (double x : point) {
+          sum += x * x;
+        }
+        return sum;
+      }
+      case FuncType::kSinProduct: {
+        double prod = 1.0;
+        for (double x : point) {
+          prod *= std::sin(x);
+        }
+        return prod;
+      }
+      default:
+        return 0.0;
+    }
+  }
+
+  static double AnalyticalIntegral(FuncType func_type, const std::vector<double> &lower,
+                                   const std::vector<double> &upper) {
+    size_t dim = lower.size();
+    switch (func_type) {
+      case FuncType::kConstant:
+        return ComputeVolume(lower, upper);
+      case FuncType::kLinear:
+        return ComputeLinearIntegral(lower, upper, dim);
+      case FuncType::kProduct:
+        return ComputeProductIntegral(lower, upper, dim);
+      case FuncType::kSumSquares:
+        return ComputeSumSquaresIntegral(lower, upper, dim);
+      case FuncType::kSinProduct:
+        return ComputeSinProductIntegral(lower, upper, dim);
+      default:
+        return 0.0;
+    }
+  }
+
+ private:
+  static double ComputeVolume(const std::vector<double> &lower, const std::vector<double> &upper) {
+    double vol = 1.0;
+    for (size_t i = 0; i < lower.size(); ++i) {
+      vol *= (upper[i] - lower[i]);
+    }
+    return vol;
+  }
+
+  static double VolumeExcludingDim(const std::vector<double> &lower, const std::vector<double> &upper, size_t exclude) {
+    double product = 1.0;
+    for (size_t j = 0; j < lower.size(); ++j) {
+      if (j != exclude) {
+        product *= (upper[j] - lower[j]);
+      }
+    }
+    return product;
+  }
+
+  static double ComputeLinearIntegral(const std::vector<double> &lower, const std::vector<double> &upper, size_t dim) {
+    double result = 0.0;
+    for (size_t i = 0; i < dim; ++i) {
+      double term = (upper[i] * upper[i] - lower[i] * lower[i]) / 2.0;
+      result += term * VolumeExcludingDim(lower, upper, i);
+    }
+    return result;
+  }
+
+  static double ComputeProductIntegral(const std::vector<double> &lower, const std::vector<double> &upper, size_t dim) {
+    double prod = 1.0;
+    for (size_t i = 0; i < dim; ++i) {
+      prod *= (upper[i] * upper[i] - lower[i] * lower[i]) / 2.0;
+    }
+    return prod;
+  }
+
+  static double ComputeSumSquaresIntegral(const std::vector<double> &lower, const std::vector<double> &upper,
+                                          size_t dim) {
+    double result = 0.0;
+    for (size_t i = 0; i < dim; ++i) {
+      double term = (upper[i] * upper[i] * upper[i] - lower[i] * lower[i] * lower[i]) / 3.0;
+      result += term * VolumeExcludingDim(lower, upper, i);
+    }
+    return result;
+  }
+
+  static double ComputeSinProductIntegral(const std::vector<double> &lower, const std::vector<double> &upper,
+                                          size_t dim) {
+    double prod = 1.0;
+    for (size_t i = 0; i < dim; ++i) {
+      prod *= (-std::cos(upper[i]) + std::cos(lower[i]));
+    }
+    return prod;
+  }
+};
+
+using InType = std::tuple<std::vector<double>, std::vector<double>, int, FuncType>;
+using OutType = double;
+using TestType = std::tuple<InType, std::string>;
+using BaseTask = ppc::task::Task<InType, OutType>;
+
+}  // namespace shilin_n_monte_carlo_integration

tasks/shilin_n_monte_carlo_integration/info.json

@@ -0,0 +1,9 @@
+{
+  "student": {
+    "first_name": "Никита",
+    "group_number": "3823Б1ПР1",
+    "last_name": "Шилин",
+    "middle_name": "Дмитриевич",
+    "task_number": "1"
+  }
+}

tasks/shilin_n_monte_carlo_integration/seq/include/ops_seq.hpp

@@ -0,0 +1,29 @@
+#pragma once
+
+#include <vector>
+
+#include "shilin_n_monte_carlo_integration/common/include/common.hpp"
+#include "task/include/task.hpp"
+
+namespace shilin_n_monte_carlo_integration {
+
+class ShilinNMonteCarloIntegrationSEQ : public BaseTask {
+ public:
+  static constexpr ppc::task::TypeOfTask GetStaticTypeOfTask() {
+    return ppc::task::TypeOfTask::kSEQ;
+  }
+  explicit ShilinNMonteCarloIntegrationSEQ(const InType &in);
+
+ private:
+  bool ValidationImpl() override;
+  bool PreProcessingImpl() override;
+  bool RunImpl() override;
+  bool PostProcessingImpl() override;
+
+  std::vector<double> lower_bounds_;
+  std::vector<double> upper_bounds_;
+  int num_points_{};
+  FuncType func_type_{};
+};
+
+}  // namespace shilin_n_monte_carlo_integration

tasks/shilin_n_monte_carlo_integration/seq/src/ops_seq.cpp

@@ -0,0 +1,90 @@
+#include "shilin_n_monte_carlo_integration/seq/include/ops_seq.hpp"
+
+#include <cstddef>
+#include <vector>
+
+#include "shilin_n_monte_carlo_integration/common/include/common.hpp"
+
+namespace shilin_n_monte_carlo_integration {
+
+ShilinNMonteCarloIntegrationSEQ::ShilinNMonteCarloIntegrationSEQ(const InType &in) {
+  SetTypeOfTask(GetStaticTypeOfTask());
+  GetInput() = in;
+  GetOutput() = 0.0;
+}
+
+bool ShilinNMonteCarloIntegrationSEQ::ValidationImpl() {
+  const auto &[lower, upper, n, func_type] = GetInput();
+  if (lower.size() != upper.size() || lower.empty()) {
+    return false;
+  }
+  if (n <= 0) {
+    return false;
+  }
+  for (size_t i = 0; i < lower.size(); ++i) {
+    if (lower[i] >= upper[i]) {
+      return false;
+    }
+  }
+  if (func_type < FuncType::kConstant || func_type > FuncType::kSinProduct) {
+    return false;
+  }
+  constexpr size_t kMaxDimensions = 10;
+  return lower.size() <= kMaxDimensions;
+}
+
+bool ShilinNMonteCarloIntegrationSEQ::PreProcessingImpl() {
+  const auto &[lower, upper, n, func_type] = GetInput();
+  lower_bounds_ = lower;
+  upper_bounds_ = upper;
+  num_points_ = n;
+  func_type_ = func_type;
+  return true;
+}
+
+bool ShilinNMonteCarloIntegrationSEQ::RunImpl() {
+  auto dimensions = static_cast<int>(lower_bounds_.size());
+
+  // Quasi-random Kronecker sequence: alpha_d = fractional part of sqrt(prime_d)
+  const std::vector<double> alpha = {
+      0.41421356237309504,  // frac(sqrt(2))
+      0.73205080756887729,  // frac(sqrt(3))
+      0.23606797749978969,  // frac(sqrt(5))
+      0.64575131106459059,  // frac(sqrt(7))
+      0.31662479035539984,  // frac(sqrt(11))
+      0.60555127546398929,  // frac(sqrt(13))
+      0.12310562561766059,  // frac(sqrt(17))
+      0.35889894354067355,  // frac(sqrt(19))
+      0.79583152331271838,  // frac(sqrt(23))
+      0.38516480713450403   // frac(sqrt(29))
+  };
+
+  std::vector<double> current(dimensions, 0.5);
+  std::vector<double> point(dimensions);
+  double sum = 0.0;
+
+  for (int i = 0; i < num_points_; ++i) {
+    for (int di = 0; di < dimensions; ++di) {
+      current[di] += alpha[di];
+      if (current[di] >= 1.0) {
+        current[di] -= 1.0;
+      }
+      point[di] = lower_bounds_[di] + ((upper_bounds_[di] - lower_bounds_[di]) * current[di]);
+    }
+    sum += IntegrandFunction::Evaluate(func_type_, point);
+  }
+
+  double volume = 1.0;
+  for (int di = 0; di < dimensions; ++di) {
+    volume *= (upper_bounds_[di] - lower_bounds_[di]);
+  }
+
+  GetOutput() = volume * sum / static_cast<double>(num_points_);
+  return true;
+}
+
+bool ShilinNMonteCarloIntegrationSEQ::PostProcessingImpl() {
+  return true;
+}
+
+}  // namespace shilin_n_monte_carlo_integration

tasks/shilin_n_monte_carlo_integration/settings.json

@@ -0,0 +1,6 @@
+{
+  "tasks": {
+    "seq": "enabled"
+  },
+  "tasks_type": "threads"
+}

tasks/shilin_n_monte_carlo_integration/tests/functional/main.cpp

@@ -0,0 +1,92 @@
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <string>
+#include <tuple>
+#include <vector>
+
+#include "shilin_n_monte_carlo_integration/common/include/common.hpp"
+#include "shilin_n_monte_carlo_integration/seq/include/ops_seq.hpp"
+#include "util/include/func_test_util.hpp"
+#include "util/include/util.hpp"
+
+namespace shilin_n_monte_carlo_integration {
+
+class ShilinNRunFuncTestsThreads : public ppc::util::BaseRunFuncTests<InType, OutType, TestType> {
+ public:
+  static std::string PrintTestParam(const TestType &test_param) {
+    return std::get<1>(test_param);
+  }
+
+ protected:
+  void SetUp() override {
+    TestType params = std::get<static_cast<std::size_t>(ppc::util::GTestParamIndex::kTestParams)>(GetParam());
+    input_data_ = std::get<0>(params);
+  }
+
+  bool CheckTestOutputData(OutType &output_data) final {
+    const auto &[lower, upper, n, func_type] = input_data_;
+    double expected = IntegrandFunction::AnalyticalIntegral(func_type, lower, upper);
+
+    double volume = 1.0;
+    for (size_t i = 0; i < lower.size(); ++i) {
+      volume *= (upper[i] - lower[i]);
+    }
+    double epsilon = std::max(volume * 10.0 / std::sqrt(static_cast<double>(n)), 1e-2);
+    return std::abs(output_data - expected) <= epsilon;
+  }
+
+  InType GetTestInputData() final {
+    return input_data_;
+  }
+
+ private:
+  InType input_data_;
+};
+
+namespace {
+
+TEST_P(ShilinNRunFuncTestsThreads, MonteCarloIntegration) {
+  ExecuteTest(GetParam());
+}
+
+const std::array<TestType, 8> kTestParam = {{
+    std::make_tuple(std::make_tuple(std::vector<double>{0.0}, std::vector<double>{1.0}, 10000, FuncType::kLinear),
+                    "kLinear_1D"),
+    std::make_tuple(std::make_tuple(std::vector<double>{0.0}, std::vector<double>{2.0}, 10000, FuncType::kSumSquares),
+                    "kSumSquares_1D"),
+    std::make_tuple(
+        std::make_tuple(std::vector<double>{0.0, 0.0}, std::vector<double>{1.0, 1.0}, 50000, FuncType::kConstant),
+        "kConstant_2D"),
+    std::make_tuple(
+        std::make_tuple(std::vector<double>{0.0, 0.0}, std::vector<double>{1.0, 1.0}, 50000, FuncType::kLinear),
+        "kLinear_2D"),
+    std::make_tuple(
+        std::make_tuple(std::vector<double>{0.0, 0.0}, std::vector<double>{1.0, 1.0}, 50000, FuncType::kProduct),
+        "kProduct_2D"),
+    std::make_tuple(
+        std::make_tuple(std::vector<double>{0.0, 0.0}, std::vector<double>{1.0, 1.0}, 50000, FuncType::kSinProduct),
+        "kSinProduct_2D"),
+    std::make_tuple(std::make_tuple(std::vector<double>{0.0, 0.0, 0.0}, std::vector<double>{1.0, 1.0, 1.0}, 100000,
+                                    FuncType::kLinear),
+                    "kLinear_3D"),
+    std::make_tuple(std::make_tuple(std::vector<double>{0.0, 0.0, 0.0}, std::vector<double>{1.0, 1.0, 1.0}, 100000,
+                                    FuncType::kProduct),
+                    "kProduct_3D"),
+}};
+
+const auto kTestTasksList = ppc::util::AddFuncTask<ShilinNMonteCarloIntegrationSEQ, InType>(
+    kTestParam, PPC_SETTINGS_shilin_n_monte_carlo_integration);
+
+const auto kGtestValues = ppc::util::ExpandToValues(kTestTasksList);
+
+const auto kPerfTestName = ShilinNRunFuncTestsThreads::PrintFuncTestName<ShilinNRunFuncTestsThreads>;
+
+INSTANTIATE_TEST_SUITE_P(MonteCarloTests, ShilinNRunFuncTestsThreads, kGtestValues, kPerfTestName);
+
+}  // namespace
+
+}  // namespace shilin_n_monte_carlo_integration