parallelism_usage_demo.cpp

#include <diffeq.hpp>
#include <vector>
#include <iostream>
#include <chrono>
#include <random>
#include <thread>
#include <future>
#include <algorithm>
#include <execution>
#include <ranges>
#include <numeric>
 
void modern_parallel_integration_example() {
    std::cout << "\n=== Modern Parallel Integration Example ===" << std::endl;
    
    // Example: Parallel ODE integration for multiple initial conditions
    std::vector<std::vector<double>> initial_conditions;
    for (int i = 0; i < 100; ++i) {
        initial_conditions.push_back({static_cast<double>(i), 0.0});
    }
    
    // Simple exponential decay: dy/dt = -0.1 * y
    auto system = [](double t, const std::vector<double>& y, std::vector<double>& dydt) {
        dydt[0] = -0.1 * y[0];
        dydt[1] = -0.2 * y[1];
    };
    
    std::cout << "Integrating " << initial_conditions.size() << " initial conditions in parallel..." << std::endl;
    
    auto start_time = std::chrono::high_resolution_clock::now();
    
    // Use standard C++20 parallel execution
    std::for_each(std::execution::par, initial_conditions.begin(), initial_conditions.end(),
        [&](std::vector<double>& state) {
            auto integrator = std::make_unique<diffeq::RK45Integrator<std::vector<double>>>(system);
            integrator->step(state, 0.01); // Single integration step
        });
    
    auto end_time = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end_time - start_time);
    
    std::cout << "[PASS] Parallel integration completed in " << duration.count() << " us!" << std::endl;
    std::cout << "Result for initial condition 10: [" << initial_conditions[10][0] 
              << ", " << initial_conditions[10][1] << "]" << std::endl;
}
 
namespace robotics_control {
 
// Robot arm dynamics (simplified)
struct RobotArmSystem {
    double mass = 1.0;
    double length = 0.5;
    double damping = 0.1;
    double gravity = 9.81;
    
    void operator()(double t, const std::vector<double>& state, std::vector<double>& derivative) {
        // state = [angle, angular_velocity]
        // Simple pendulum with control input
        double angle = state[0];
        double angular_vel = state[1];
        
        // Control input (placeholder - would come from feedback controller)
        double control_torque = -2.0 * angle - 0.5 * angular_vel;  // PD controller
        
        derivative[0] = angular_vel;
        derivative[1] = -(gravity / length) * std::sin(angle) - 
                        (damping / (mass * length * length)) * angular_vel + 
                        control_torque / (mass * length * length);
    }
};
 
void demonstrate_realtime_control() {
    std::cout << "\n=== Robotics Control System with Modern Async ===" << std::endl;
    
    // Setup multiple control systems (e.g., different robot joints)
    std::vector<std::vector<double>> joint_states;
    for (int i = 0; i < 6; ++i) {  // 6-DOF robot arm
        joint_states.push_back({0.1 * i, 0.0}); // [angle, angular_velocity]
    }
    
    auto start_time = std::chrono::high_resolution_clock::now();
    
    // Parallel control loop using standard C++ parallelism
    std::for_each(std::execution::par, joint_states.begin(), joint_states.end(),
        [](std::vector<double>& state) {
            RobotArmSystem system;
            auto integrator = std::make_unique<diffeq::RK45Integrator<std::vector<double>>>(system);
            integrator->step(state, 0.001); // 1ms control timestep
        });
    
    auto end_time = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end_time - start_time);
    
    std::cout << "Parallel control completed in " << duration.count() << " us" << std::endl;
    std::cout << "Average per joint: " << duration.count() / joint_states.size() << " us" << std::endl;
    
    // Advanced async approach using diffeq async capabilities
    std::cout << "\n--- Advanced Async Approach ---" << std::endl;
    
    // Create integrator for robot dynamics
    auto robot_system = RobotArmSystem{};
    auto integrator = std::make_unique<diffeq::RK45Integrator<std::vector<double>>>(robot_system);
    
    // Initial state: [angle=0.1 rad, angular_velocity=0]
    std::vector<double> state = {0.1, 0.0};
    double dt = 0.001;  // 1ms timestep for 1kHz control
    double simulation_time = 1.0;
    
    std::cout << "Running real-time robot control simulation..." << std::endl;
    std::cout << "Target frequency: 1kHz (1ms timestep)" << std::endl;
    
    auto advanced_start_time = std::chrono::high_resolution_clock::now();
    
    // Simulate real-time control loop with async execution
    std::vector<std::future<void>> control_futures;
    
    for (double t = 0.0; t < simulation_time; t += dt) {
        // Execute control step asynchronously
        auto control_future = std::async(std::launch::async, [&, t]() {
            integrator->step(state, dt);
        });
        
        // Simulate sensor reading (parallel)
        auto sensor_future = std::async(std::launch::async, [&]() {
            // Placeholder for sensor data processing
            std::this_thread::sleep_for(std::chrono::microseconds(50));
            return state[0];  // Return angle measurement
        });
        
        // Wait for both to complete
        control_future.wait();
        double measured_angle = sensor_future.get();
        
        // Output every 100ms
        if (static_cast<int>(t * 1000) % 100 == 0) {
            std::cout << "t=" << t << "s, angle=" << state[0] 
                      << " rad, measured=" << measured_angle << " rad" << std::endl;
        }
    }
    
    auto advanced_end_time = std::chrono::high_resolution_clock::now();
    auto advanced_duration = std::chrono::duration_cast<std::chrono::milliseconds>(advanced_end_time - advanced_start_time);
    
    std::cout << "Simulation completed in " << advanced_duration.count() << "ms" << std::endl;
    std::cout << "Average loop time: " << advanced_duration.count() / (simulation_time / dt) << "ms" << std::endl;
    std::cout << "Final robot state: angle=" << state[0] << " rad, velocity=" << state[1] << " rad/s" << std::endl;
}
 
} // namespace robotics_control
 
namespace stochastic_research {
 
// Geometric Brownian Motion for financial modeling
struct GeometricBrownianMotion {
    double mu = 0.05;    // drift
    double sigma = 0.2;  // volatility
    
    void operator()(double t, const std::vector<double>& state, std::vector<double>& derivative) {
        derivative[0] = mu * state[0];  // deterministic part
    }
    
    void diffusion(double t, const std::vector<double>& state, std::vector<double>& noise_coeff) {
        noise_coeff[0] = sigma * state[0];  // stochastic part
    }
};
 
void demonstrate_monte_carlo_simulation() {
    std::cout << "\n=== Stochastic Process Research with Modern Parallelism ===" << std::endl;
    
    const int num_simulations = 10000;
    const double initial_price = 100.0;
    const double dt = 0.01;
    const double t_final = 1.0;
    
    std::cout << "Running " << num_simulations << " Monte Carlo simulations..." << std::endl;
    
    auto start_time = std::chrono::high_resolution_clock::now();
    
    // Use standard C++20 parallel execution for Monte Carlo
    std::vector<double> final_prices(num_simulations);
    
    std::for_each(std::execution::par, 
        std::views::iota(0, num_simulations).begin(),
        std::views::iota(0, num_simulations).end(),
        [&](int i) {
            std::mt19937 rng(i);  // Seed with simulation index
            std::normal_distribution<double> normal(0.0, 1.0);
            
            // Define SDE parameters
            const double mu = 0.05;    // drift coefficient
            const double sigma = 0.2;  // volatility coefficient
            
            // Simple Euler-Maruyama implementation
            std::vector<double> state = {initial_price};
            for (double t = 0.0; t < t_final; t += dt) {
                double dW = normal(rng) * std::sqrt(dt);
                state[0] += mu * state[0] * dt + sigma * state[0] * dW;
            }
            
            final_prices[i] = state[0];
        });
    
    auto end_time = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time);
    
    // Calculate statistics
    double mean_price = std::reduce(final_prices.begin(), final_prices.end()) / num_simulations;
    double variance = 0.0;
    for (double price : final_prices) {
        variance += (price - mean_price) * (price - mean_price);
    }
    variance /= num_simulations;
    double std_dev = std::sqrt(variance);
    
    std::cout << "Monte Carlo simulation completed in " << duration.count() << "ms" << std::endl;
    std::cout << "Mean final price: " << mean_price << std::endl;
    std::cout << "Standard deviation: " << std_dev << std::endl;
    std::cout << "Theoretical mean: " << initial_price * std::exp(0.05) << std::endl;
}
 
} // namespace stochastic_research
 
namespace hardware_benchmark {
 
struct ExponentialDecay {
    double k = 1.0;
    
    void operator()(double t, const std::vector<double>& state, std::vector<double>& derivative) {
        derivative[0] = -k * state[0];
    }
};
 
void benchmark_hardware_targets() {
    std::cout << "\n=== Hardware Performance Benchmarking ===" << std::endl;
    
    const int num_integrations = 1000;  // Reduced from 10000 to 1000 for faster execution
    const double dt = 0.01;
    const double t_final = 1.0;
    
    auto system = ExponentialDecay{};
    auto integrator = std::make_unique<diffeq::RK45Integrator<std::vector<double>>>(system);
    std::vector<double> initial_state = {1.0};
    
    // Sequential execution
    auto seq_start = std::chrono::high_resolution_clock::now();
    std::vector<std::vector<double>> seq_states(num_integrations, initial_state);
            for (int i = 0; i < num_integrations; ++i) {
            integrator->integrate(seq_states[i], dt, t_final);
        }
    auto seq_end = std::chrono::high_resolution_clock::now();
    auto seq_duration = std::chrono::duration_cast<std::chrono::milliseconds>(seq_end - seq_start);
    
    // Parallel execution
    auto par_start = std::chrono::high_resolution_clock::now();
    std::vector<std::vector<double>> par_states(num_integrations, initial_state);
    std::for_each(std::execution::par, par_states.begin(), par_states.end(),
        [&](std::vector<double>& state) {
            auto local_integrator = std::make_unique<diffeq::RK45Integrator<std::vector<double>>>(system);
            local_integrator->integrate(state, dt, t_final);
        });
    auto par_end = std::chrono::high_resolution_clock::now();
    auto par_duration = std::chrono::duration_cast<std::chrono::milliseconds>(par_end - par_start);
    
    std::cout << "Sequential execution: " << seq_duration.count() << "ms" << std::endl;
    std::cout << "Parallel execution: " << par_duration.count() << "ms" << std::endl;
    std::cout << "Speedup: " << static_cast<double>(seq_duration.count()) / par_duration.count() << "x" << std::endl;
    std::cout << "Hardware concurrency: " << std::thread::hardware_concurrency() << " threads" << std::endl;
}
 
} // namespace hardware_benchmark
 
void demonstrate_all_parallelism_features() {
    std::cout << "\n=== Complete Parallelism Feature Demonstration ===" << std::endl;
    
    modern_parallel_integration_example();
    robotics_control::demonstrate_realtime_control();
    stochastic_research::demonstrate_monte_carlo_simulation();
    hardware_benchmark::benchmark_hardware_targets();
    
    std::cout << "\n=== All Examples Completed Successfully! ===" << std::endl;
}
 
int main() {
    std::cout << "Modern DiffeQ Parallelism Examples" << std::endl;
    std::cout << "===================================" << std::endl;
    
    demonstrate_all_parallelism_features();
    
    return 0;
}