added multi strat functions

Author: Alexander Robbins

include/lfmc/estimator/control_variate_estimator.hpp

@@ -2,7 +2,8 @@
 
 #include "lfmc/core/types.hpp"
 #include "lfmc/estimator/estimator.hpp"
-
+#include <cmath>
+#include <limits>
 #include <expected>
 #include <string>
 #include <vector>
@@ -18,6 +19,7 @@ class ControlVariateEstimator : public Estimator {
     double sum_y = 0.0;  // Sum of control variate values
     double sum_xy = 0.0; // Sum of products of payoffs and control variate
     double sum_yy = 0.0; // Sum of squares of control variate values
+    double sum_xx = 0.0;
 
     double control_expectation_ = 0.0; // Expected value of control variate (known analytically)
 
@@ -29,6 +31,11 @@ class ControlVariateEstimator : public Estimator {
     bool converged() const override;
     std::expected<double, std::string> result() const override;
     std::expected<void, std::string> merge(Estimator const& other) override;
+
+    double mean() const noexcept;
+    double variance() const noexcept;
+    double std_error() const noexcept;
+    std::size_t sample_count() const noexcept;
 };
 
 } // namespace lfmc

include/lfmc/estimator/monte_carlo_estimator.hpp

@@ -12,13 +12,19 @@ namespace lfmc {
 class MonteCarloEstimator : public Estimator {
   private:
     double sum = 0.0;
+    double sum_sq = 0.0;
     std::size_t count = 0;
 
   public:
     std::expected<void, std::string> add_payoffs(const std::vector<Payoffs>& payoffs) override;
     bool converged() const override;
     std::expected<double, std::string> result() const override;
     std::expected<void, std::string> merge(Estimator const& other) override;
+
+    double mean() const noexcept;
+    double variance() const noexcept;
+    double std_error() const noexcept;
+    std::size_t sample_count() const noexcept;
 };
 
 } // namespace lfmc

include/lfmc/path_generator/path_generator.hpp

@@ -27,8 +27,8 @@ class PathGenerator {
 
         std::vector<Path> paths;
         for (const auto& norm : normals) {
-            Path path(steps + 1);
-
+            Path path;
+            path.reserve(steps + 1);  
             double t = 0.0;
             double x = process_.initial();
             path.push_back(x);

include/lfmc/pipeline/pipeline.hpp

@@ -49,6 +49,15 @@ template <StochasticProcess SP, NumericalScheme<SP> NS> class Pipeline {
 
         return estimator_->result();
     }
+    
+    const Estimator* get_estimator() const noexcept {
+        return estimator_.get();
+    }
+
+    Estimator* get_estimator() noexcept {
+        return estimator_.get();
+    }
+
 };
 
 } // namespace lfmc

include/lfmc/strategy/multi_strategy_runner.hpp

@@ -0,0 +1,77 @@
+#pragma once
+
+#include "lfmc/strategy/strategy_metrics.hpp"
+#include "lfmc/strategy/strategy_runner.hpp"
+
+#include <algorithm>
+#include <memory>
+#include <thread>
+#include <vector>
+
+namespace lfmc {
+
+template <StochasticProcess SP, NumericalScheme<SP> NS>
+class MultiStrategyRunner {
+private:
+    std::vector<std::unique_ptr<StrategyRunner<SP, NS>>> strategies;
+    std::vector<std::thread> threads;
+
+public:
+    void add_strategy(std::unique_ptr<StrategyRunner<SP, NS>> strategy) {
+        strategies.push_back(std::move(strategy));
+    }
+
+    // Run all strategies for a warmup period
+    std::expected<void, std::string> run_warmup(size_t steps, double T, 
+                                                  size_t warmup_iterations) {
+        threads.clear();
+        threads.reserve(strategies.size());
+
+        for (auto& strategy : strategies) {
+            threads.emplace_back([&strategy, steps, T, warmup_iterations]() {
+                for (size_t i = 0; i < warmup_iterations; ++i) {
+                    auto result = strategy->run_batch(steps, T);
+                    if (!result) {
+                        // Error occurred - need this to get log
+                        return;
+                    }
+                }
+            });
+        }
+
+        for (auto& thread : threads) {
+            if (thread.joinable()) {
+                thread.join();
+            }
+        }
+
+        return {};
+    }
+
+    // Get current metricss
+    std::vector<StrategyMetrics> get_all_metrics() const {
+        std::vector<StrategyMetrics> metrics;
+        metrics.reserve(strategies.size());
+        
+        for (const auto& strategy : strategies) {
+            metrics.push_back(strategy->get_metrics());
+        }
+        
+        return metrics;
+    }
+
+    // Find best strat
+    size_t best_strategy_index() const {
+        auto metrics = get_all_metrics();
+        if (metrics.empty()) {
+            return 0;
+        }
+        return std::min_element(metrics.begin(), metrics.end()) - metrics.begin();
+    }
+
+    size_t strategy_count() const noexcept { 
+        return strategies.size(); 
+    }
+};
+
+} // namespace lfmc

include/lfmc/strategy/strategy_factory.hpp

@@ -0,0 +1,79 @@
+#pragma once
+
+#include "lfmc/estimator/control_variate_estimator.hpp"
+#include "lfmc/estimator/monte_carlo_estimator.hpp"
+#include "lfmc/numerical_scheme/euler_maruyama.hpp"
+#include "lfmc/path_generator/path_generator.hpp"
+#include "lfmc/payoff/control_variate_payoffs.hpp"
+#include "lfmc/payoff/payoff.hpp"
+#include "lfmc/pipeline/pipeline.hpp"
+#include "lfmc/random_source/antithetic_random_source.hpp"
+#include "lfmc/random_source/pseudo_random_source.hpp"
+#include "lfmc/strategy/strategy_runner.hpp"
+#include "lfmc/stochastic_process/stochastic_process.hpp"
+
+#include <memory>
+#include <random>
+#include <string>
+
+namespace lfmc {
+
+template <StochasticProcess SP, NumericalScheme<SP> NS>
+class StrategyFactory {
+public:
+    static std::unique_ptr<StrategyRunner<SP, NS>>
+    create_pseudo_random_strategy(SP process, NS scheme, 
+                                   std::unique_ptr<Payoff> payoff,
+                                   unsigned seed = std::random_device{}()) {
+        auto rs = std::make_unique<PseudoRandomSource>(seed);
+        auto pg = std::make_unique<PathGenerator<SP, NS>>(process, scheme);
+        auto est = std::make_unique<MonteCarloEstimator>();
+
+        Pipeline<SP, NS> pipeline(std::move(rs), std::move(pg), 
+                                  std::move(payoff), std::move(est));
+
+        return std::make_unique<StrategyRunner<SP, NS>>(
+            std::move(pipeline), "PseudoRandom");
+    }
+
+    static std::unique_ptr<StrategyRunner<SP, NS>>
+    create_antithetic_strategy(SP process, NS scheme,
+                               std::unique_ptr<Payoff> payoff,
+                               unsigned seed = std::random_device{}()) {
+        auto rs = std::make_unique<AntitheticRandomSource>(seed);
+        auto pg = std::make_unique<PathGenerator<SP, NS>>(process, scheme);
+        auto est = std::make_unique<MonteCarloEstimator>();
+
+        Pipeline<SP, NS> pipeline(std::move(rs), std::move(pg),
+                                  std::move(payoff), std::move(est));
+
+        return std::make_unique<StrategyRunner<SP, NS>>(
+            std::move(pipeline), "Antithetic");
+    }
+
+    static std::unique_ptr<StrategyRunner<SP, NS>>
+    create_control_variate_strategy(SP process, NS scheme,
+                                    std::unique_ptr<Payoff> target_payoff,
+                                    std::unique_ptr<Payoff> control_payoff,
+                                    double control_expectation,
+                                    unsigned seed = std::random_device{}()) {
+        auto rs = std::make_unique<PseudoRandomSource>(seed);
+        auto pg = std::make_unique<PathGenerator<SP, NS>>(process, scheme);
+        auto po = std::make_unique<ControlVariatePayoff>(
+            std::move(target_payoff), std::move(control_payoff));
+        auto est = std::make_unique<ControlVariateEstimator>(control_expectation);
+
+        Pipeline<SP, NS> pipeline(std::move(rs), std::move(pg),
+                                  std::move(po), std::move(est));
+
+        return std::make_unique<StrategyRunner<SP, NS>>(
+            std::move(pipeline), "ControlVariate");
+    }
+
+    // Add more factory methods as we make/implement more VR strategies:
+    // - create_quasi_monte_carlo_strategy()
+    // - create_importance_sampling_strategy()
+    // etc.
+};
+
+} // namespace lfmc

include/lfmc/strategy/strategy_metrics.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <cstddef>
+#include <string>
+
+namespace lfmc {
+
+struct StrategyMetrics {
+    std::string name;
+    double mean;
+    double variance;
+    double std_error;
+    std::size_t samples;
+    long long elapsed_ms;
+
+    bool operator<(const StrategyMetrics& other) const noexcept {
+        return variance < other.variance;
+    }
+};
+
+} // namespace lfmc

include/lfmc/strategy/strategy_runner.hpp

@@ -0,0 +1,72 @@
+#pragma once
+
+#include "lfmc/estimator/estimator.hpp"
+#include "lfmc/estimator/monte_carlo_estimator.hpp"
+#include "lfmc/estimator/control_variate_estimator.hpp"
+#include "lfmc/pipeline/pipeline.hpp"
+#include "lfmc/strategy/strategy_metrics.hpp"
+#include "lfmc/timing/timing.hpp"
+
+#include <atomic>
+#include <memory>
+#include <string>
+
+namespace lfmc {
+
+template <StochasticProcess SP, NumericalScheme<SP> NS>
+class StrategyRunner {
+private:
+    Pipeline<SP, NS> pipeline;
+    std::string name;
+    Timer timer;
+    
+    std::atomic<double> current_mean{0.0};
+    std::atomic<double> current_variance{std::numeric_limits<double>::max()};
+    std::atomic<size_t> samples_processed{0};
+
+public:
+    StrategyRunner(Pipeline<SP, NS> pipeline, std::string name)
+        : pipeline(std::move(pipeline)), name(std::move(name)) {}
+
+    // Run a single batch of simulations and update metrics I made in the other metrics file
+    std::expected<void, std::string> run_batch(size_t steps, double T) {
+        auto result = pipeline.run(steps, T);
+        if (!result) {
+            return std::unexpected(result.error());
+        }
+        
+        
+        auto* estimator = pipeline.get_estimator();
+        
+        // had to implement this, it wasn't in the base class, but need it to update the metrics
+        if (auto* mc_est = dynamic_cast<MonteCarloEstimator*>(estimator)) {
+            current_mean.store(mc_est->mean(), std::memory_order_relaxed);
+            current_variance.store(mc_est->variance(), std::memory_order_relaxed);
+            samples_processed.store(mc_est->sample_count(), std::memory_order_relaxed);
+        }
+        // Try ControlVariateEstimator, same as above
+        else if (auto* cv_est = dynamic_cast<ControlVariateEstimator*>(estimator)) {
+            current_mean.store(cv_est->mean(), std::memory_order_relaxed);
+            current_variance.store(cv_est->variance(), std::memory_order_relaxed);
+            samples_processed.store(cv_est->sample_count(), std::memory_order_relaxed);
+        }
+        
+        return {};
+    }
+
+    StrategyMetrics get_metrics() const noexcept {
+        return StrategyMetrics{
+            .name = name,
+            .mean = current_mean.load(std::memory_order_relaxed),
+            .variance = current_variance.load(std::memory_order_relaxed),
+            .std_error = std::sqrt(current_variance.load(std::memory_order_relaxed) / 
+                                   samples_processed.load(std::memory_order_relaxed)),
+            .samples = samples_processed.load(std::memory_order_relaxed),
+            .elapsed_ms = timer.elapsedMilliseconds()
+        };
+    }
+
+    const std::string& get_name() const noexcept { return name; }
+};
+
+} // namespace lfmc

src/estimator/control_variate_estimator.cpp

@@ -16,6 +16,7 @@ ControlVariateEstimator::add_payoffs(const std::vector<Payoffs>& payoffs) {
         sum_y += y;
         sum_xy += x * y;
         sum_yy += y * y;
+        sum_xx += x * x;
         ++count;
     }
 
@@ -61,8 +62,53 @@ std::expected<void, std::string> ControlVariateEstimator::merge(Estimator const&
     sum_y += other_estimator->sum_y;
     sum_xy += other_estimator->sum_xy;
     sum_yy += other_estimator->sum_yy;
+    sum_xx += other_estimator->sum_xx;
 
     return {};
 }
 
+double ControlVariateEstimator::mean() const noexcept {
+    if (count == 0) return 0.0;
+    
+    double n = static_cast<double>(count);
+    double mean_x = sum_x / n;
+    double mean_y = sum_y / n;
+    
+    double cov_xy = (sum_xy / n) - mean_x * mean_y;
+    double var_y = (sum_yy / n) - mean_y * mean_y;
+    
+    if (var_y == 0.0) return mean_x;
+    
+    double beta = cov_xy / var_y;
+    return mean_x - beta * (mean_y - control_expectation_);
+}
+
+double ControlVariateEstimator::variance() const noexcept {
+    if (count < 2) return std::numeric_limits<double>::max();
+    
+    double n = static_cast<double>(count);
+    double mean_x = sum_x / n;
+    double mean_y = sum_y / n;
+    
+    double var_x = (sum_xx / n) - mean_x * mean_x;
+    double var_y = (sum_yy / n) - mean_y * mean_y;
+    double cov_xy = (sum_xy / n) - mean_x * mean_y;
+    
+    if (var_y == 0.0) return var_x;
+    
+    double beta = cov_xy / var_y;
+    
+    // Var(X - β(Y - E[Y])) = Var(X) + β²Var(Y) - 2βCov(X,Y)
+    return var_x + beta * beta * var_y - 2.0 * beta * cov_xy;
+}
+
+double ControlVariateEstimator::std_error() const noexcept {
+    if (count < 2) return std::numeric_limits<double>::max();
+    return std::sqrt(variance() / static_cast<double>(count));
+}
+
+std::size_t ControlVariateEstimator::sample_count() const noexcept {
+    return count;
+}
+
 } // namespace lfmc