fatigue_optimizer.cpp - main - starlv/Fatigue_plan_optimizer_cpp

fatigue_optimizer.cpp
1 167 строк43 KB
Agamirov L.V.
Add files via upload
22 май 2026, 18:40
22 май 2026, 18:401f16d79
// fatigue_optimizer.cpp
// TWO-STAGE FATIGUE TEST PLANNING OPTIMIZATION
// Stage 1: Theoretical optimization (based on Fisher information matrix)
// Stage 2: Monte Carlo refinement

#include <cmath>
#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <algorithm>
#include <random>
#include <chrono>
#include <sstream>
#include <map>
#include <set>
#include <functional>
#include <tuple>
#include "fatique_curve.cpp"

using namespace std;


// ====================== DATA STRUCTURES ======================

struct MaterialParams {
    double sigma_inf;       // fatigue limit, MPa
    double C;               // coefficient in fatigue curve equation
    double m;               // exponent
    double a;               // lgN dispersion parameter
    double alpha;           // exponent for dispersion
    vector<double> lgN_levels;   // lgN levels
    vector<double> N0_list;      // service lives for strength evaluation
    double C1;              // cost per specimen
    double C3;              // cost per hour of testing
    double f;               // loading frequency, Hz
    double eps;             // relative error
    double gamma;           // confidence half-width
    double t_quantile;      // Student's t-quantile
    int n_simulations;      // number of Monte Carlo simulations
    
    MaterialParams() : sigma_inf(250.0), C(800.0), m(0.12), a(0.08), alpha(0.8),
                       C1(1000.0), C3(1000.0), f(1000.0), eps(0.2), gamma(0.1),
                       t_quantile(2.0), n_simulations(700) {}
    
    double delta_target() const {
        return eps * gamma / t_quantile;
    }
    
    double sigma_from_N(double N) const {
        return sigma_inf + C * pow(N, -m);
    }
    
    double sigma_from_lgN(double lgN) const {
        return sigma_from_N(pow(10.0, lgN));
    }
};

struct Plan {
    vector<double> nu;          // distribution fractions
    int n;                      // total number of specimens
    double delta;               // achieved accuracy
    double cost;                // cost
    string description;         // description
    vector<int> n_distribution; // specimen distribution across levels
};

struct Variant {
    double beta;
    vector<double> nu;
    string description;
};

// ====================== FORWARD DECLARATIONS ======================
bool read_input(const string& filename, MaterialParams& params);
double calculate_cost(const MaterialParams& params, int n, const vector<double>& nu);
double theoretical_std_sigma_R(const MaterialParams& params, int n, const vector<double>& nu, double N0);
double theoretical_delta(const MaterialParams& params, int n, const vector<double>& nu, double N0);
pair<int, double> theoretical_min_n(const MaterialParams& params, const vector<double>& nu, double N0, int n_min, int n_max);
vector<Variant> generate_nu_variants(const MaterialParams& params);
Plan theoretical_optimization(const MaterialParams& params, double N0, ostream& out_file);
Plan monte_carlo_optimization(const MaterialParams& params, const Plan& theoretical_plan, double N0, ostream& out_file);

// ====================== READING INPUT DATA ======================

bool read_input(const string& filename, MaterialParams& params) {
    ifstream file(filename);
    if (!file.is_open()) {
        cerr << "Failed to open file: " << filename << endl;
        return false;
    }
    
    string line;
    while (getline(file, line)) {
        // Skip empty lines and comments
        if (line.empty() || line[0] == '#') continue;
        
        stringstream ss(line);
        string key;
        ss >> key;
        
        if (key == "sigma_inf") {
            ss >> params.sigma_inf;
        } else if (key == "C") {
            ss >> params.C;
        } else if (key == "m") {
            ss >> params.m;
        } else if (key == "a") {
            ss >> params.a;
        } else if (key == "alpha") {
            ss >> params.alpha;
        } else if (key == "lgN_levels") {
            params.lgN_levels.clear();
            double val;
            while (ss >> val) {
                params.lgN_levels.push_back(val);
            }
        } else if (key == "N0_list") {
            params.N0_list.clear();
            double val;
            while (ss >> val) {
                params.N0_list.push_back(val);
            }
        } else if (key == "C1") {
            ss >> params.C1;
        } else if (key == "C3") {
            ss >> params.C3;
        } else if (key == "f") {
            ss >> params.f;
        } else if (key == "eps") {
            ss >> params.eps;
        } else if (key == "gamma") {
            ss >> params.gamma;
        } else if (key == "t_quantile") {
            ss >> params.t_quantile;
        } else if (key == "n_simulations") {
            ss >> params.n_simulations;
        }
    }
    
    file.close();
    return true;
}

// ====================== COST CALCULATION ======================

double calculate_cost(const MaterialParams& params, int n, const vector<double>& nu) {
    vector<int> n_i(nu.size());
    int total = 0;
    for (size_t i = 0; i < nu.size(); ++i) {
        n_i[i] = (int)round(n * nu[i]);
        total += n_i[i];
    }
    
    // Adjust sum
    if (total != n) {
        int idx = max_element(n_i.begin(), n_i.end()) - n_i.begin();
        n_i[idx] += n - total;
    }
    
    // Ensure at least 1 specimen per level (if nu[i] > 0)
    for (size_t i = 0; i < nu.size(); ++i) {
        if (nu[i] > 1e-6 && n_i[i] == 0) {
            n_i[i] = 1;
            int idx = max_element(n_i.begin(), n_i.end()) - n_i.begin();
            n_i[idx]--;
        }
    }
    
    double total_cycles = 0.0;
    for (size_t i = 0; i < params.lgN_levels.size(); ++i) {
        total_cycles += n_i[i] * pow(10.0, params.lgN_levels[i]);
    }
    
    double time_hours = total_cycles / (params.f * 3600.0);
    return n + (params.C3 / params.C1) * time_hours;
}

// ====================== THEORETICAL ESTIMATION ======================

double theoretical_std_sigma_R(const MaterialParams& params, int n, 
                                const vector<double>& nu, double N0) {
    int n_levels = params.lgN_levels.size();
    double ln10 = log(10.0);
    
    double F11 = 0, F12 = 0, F13 = 0, F22 = 0, F23 = 0, F33 = 0;
    int active_levels = 0;
    
    for (int i = 0; i < n_levels; ++i) {
        double lgN = params.lgN_levels[i];
        double nu_i = nu[i];
        
        if (nu_i < 1e-6) continue;
        
        // Use fractional number of specimens (like Python)
        double n_i_double = n * nu_i;
        active_levels++;
        
        // Variance of lgN at this level
        double sigma_lgN = params.a * pow(lgN, params.alpha);
        double var_lgN = sigma_lgN * sigma_lgN;
        
        // Weight: n_i / var_lgN
        double w = n_i_double / var_lgN;
        
        // Stress at this level
        double sigma = params.sigma_from_lgN(lgN);
        double sigma_adj = sigma - params.sigma_inf;
        sigma_adj = max(sigma_adj, 1e-6);
        
        // Derivatives
        double d_dsigma_inf = (1.0 / params.m) * (1.0 / (ln10 * sigma_adj));
        double d_dC = (1.0 / params.m) * (1.0 / (ln10 * params.C));
        double d_dm = -(1.0 / (params.m * params.m)) * 
                      (log10(params.C) - log10(sigma_adj));
        
        // Add to Fisher information matrix
        F11 += w * d_dsigma_inf * d_dsigma_inf;
        F12 += w * d_dsigma_inf * d_dC;
        F13 += w * d_dsigma_inf * d_dm;
        F22 += w * d_dC * d_dC;
        F23 += w * d_dC * d_dm;
        F33 += w * d_dm * d_dm;
    }
    
    // Need at least 3 active levels for reliable estimation
    if (active_levels < 3) {
        return 0.0;  // Return 0.0 which will become 0.25 delta
    }
    
    // Check if matrix is non-zero
    if (F11 < 1e-12 && F22 < 1e-12 && F33 < 1e-12) {
        return 0.0;
    }
    
    // Construct 3x3 matrix
    double F[3][3] = {{F11, F12, F13}, {F12, F22, F23}, {F13, F23, F33}};
    
    // Compute condition number using eigenvalues approximation
    // For a 3x3 symmetric matrix, we can compute condition number as |lambda_max| / |lambda_min|
    
    // Compute characteristic polynomial coefficients: λ³ - trace*λ² + (sum of principal minors)*λ - det = 0
    double trace = F11 + F22 + F33;
    
    // Sum of principal minors (2x2 determinants)
    double sum_minors = (F11*F22 - F12*F12) + (F11*F33 - F13*F13) + (F22*F33 - F23*F23);
    
    double det = F11 * (F22 * F33 - F23 * F23) -
                 F12 * (F12 * F33 - F23 * F13) +
                 F13 * (F12 * F23 - F22 * F13);
    
    if (fabs(det) < 1e-12) {
        return 0.0;
    }
    
    // Estimate eigenvalues using quadratic approximation for the largest eigenvalue
    // Power iteration for max eigenvalue
    double v[3] = {1.0, 1.0, 1.0};
    for (int iter = 0; iter < 20; ++iter) {
        double w[3] = {0, 0, 0};
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                w[i] += F[i][j] * v[j];
            }
        }
        double norm = sqrt(w[0]*w[0] + w[1]*w[1] + w[2]*w[2]);
        for (int i = 0; i < 3; ++i) {
            v[i] = w[i] / norm;
        }
    }
    
    double lambda_max = 0;
    for (int i = 0; i < 3; ++i) {
        double val = 0;
        for (int j = 0; j < 3; ++j) {
            val += F[i][j] * v[j];
        }
        lambda_max += v[i] * val;
    }
    
    // Estimate min eigenvalue using inverse iteration
    // [F] * x = y, solve for x
    double det_inv = 1.0 / det;
    double invF[3][3] = {
        {(F22*F33 - F23*F23) * det_inv, (F13*F23 - F12*F33) * det_inv, (F12*F23 - F13*F22) * det_inv},
        {(F13*F23 - F12*F33) * det_inv, (F11*F33 - F13*F13) * det_inv, (F12*F13 - F11*F23) * det_inv},
        {(F12*F23 - F13*F22) * det_inv, (F12*F13 - F11*F23) * det_inv, (F11*F22 - F12*F12) * det_inv}
    };
    
    // Power iteration on inverse matrix to get smallest eigenvalue of F
    double v_inv[3] = {1.0, 1.0, 1.0};
    for (int iter = 0; iter < 20; ++iter) {
        double w[3] = {0, 0, 0};
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                w[i] += invF[i][j] * v_inv[j];
            }
        }
        double norm = sqrt(w[0]*w[0] + w[1]*w[1] + w[2]*w[2]);
        for (int i = 0; i < 3; ++i) {
            v_inv[i] = w[i] / norm;
        }
    }
    
    double lambda_min_inv = 0;
    for (int i = 0; i < 3; ++i) {
        double val = 0;
        for (int j = 0; j < 3; ++j) {
            val += invF[i][j] * v_inv[j];
        }
        lambda_min_inv += v_inv[i] * val;
    }
    
    double lambda_min = 1.0 / lambda_min_inv;
    double cond = lambda_max / lambda_min;
    
    // Python uses cond > 1e12 as threshold
    if (cond > 1e12) {
        return 0.0;  // Return 0.0 which will become 0.25 delta
    }
    
    // Invert 3x3 matrix for covariance
    double invF11 = (F22 * F33 - F23 * F23) / det;
    double invF12 = (F13 * F23 - F12 * F33) / det;
    double invF13 = (F12 * F23 - F13 * F22) / det;
    double invF22 = (F11 * F33 - F13 * F13) / det;
    double invF23 = (F12 * F13 - F11 * F23) / det;
    double invF33 = (F11 * F22 - F12 * F12) / det;
    
    // Check if inverse has negative diagonal elements (matrix not positive definite)
    if (invF11 < 0 || invF22 < 0 || invF33 < 0) {
        return 0.0;
    }
    
    // Derivatives of fatigue limit
    double term = pow(N0, -params.m);
    double dR_dsigma_inf = 1.0;
    double dR_dC = term;
    double dR_dm = -params.C * term * log(N0);
    
    // Variance of fatigue limit
    double var = dR_dsigma_inf * (dR_dsigma_inf * invF11 + dR_dC * invF12 + dR_dm * invF13) +
                 dR_dC * (dR_dsigma_inf * invF12 + dR_dC * invF22 + dR_dm * invF23) +
                 dR_dm * (dR_dsigma_inf * invF13 + dR_dC * invF23 + dR_dm * invF33);
    
    if (var < 0) var = 0;
    double result = sqrt(var);
    
    // Additional sanity check
    double sigma_R = params.sigma_from_N(N0);
    if (result < 0.1 || result > 100) {
        return 0.0;
    }
    
    return result;
}

//=====================================================================

double theoretical_delta(const MaterialParams& params, int n,
                          const vector<double>& nu, double N0) {
    double true_val = params.sigma_from_N(N0);
    double std_theo = theoretical_std_sigma_R(params, n, nu, N0);
    
    if (true_val <= 0) return 0.25;
    if (std_theo <= 1e-6) return 0.25;  // Если вернулся 0.0 (сингулярная матрица)
    
    double delta = std_theo / true_val;
    
    // Bound reasonable values
    if (delta < 0.001) return 0.25;
    if (delta > 0.5) return 0.25;
    
    return delta;
}

//======================================================================

pair<int, double> theoretical_min_n(const MaterialParams& params,
                                     const vector<double>& nu, double N0,
                                     int n_min=1, int n_max=300) {
    double delta_target = params.delta_target();
    
    // Count levels with significant nu
    int levels_with_nu = 0;
    for (double nu_i : nu) {
        if (nu_i > 1e-6) levels_with_nu++;
    }
    
    // Minimum n should be at least number of active levels
    // and at least 3 for fitting
    int min_required = max(levels_with_nu, 3);
    n_min = max(n_min, min_required);
    
    // Start searching from n_min
    for (int n = n_min; n <= min(n_min + 30, n_max); ++n) {
        double delta = theoretical_delta(params, n, nu, N0);
        if (delta <= delta_target && delta > 0 && delta < 0.1) {
            // Found a solution, try to decrease
            for (int n2 = max(n_min, n - 10); n2 <= n; ++n2) {
                double delta2 = theoretical_delta(params, n2, nu, N0);
                if (delta2 <= delta_target && delta2 > 0 && delta2 < 0.1) {
                    return {n2, delta2};
                }
            }
            return {n, delta};
        }
    }
    
    // Binary search for larger n
    int left = n_min, right = n_max;
    int best_n = n_max;
    double best_delta = theoretical_delta(params, n_max, nu, N0);
    
    while (left <= right) {
        int mid = (left + right) / 2;
        double delta = theoretical_delta(params, mid, nu, N0);
        if (delta <= delta_target && delta > 0 && delta < 0.1) {
            best_n = mid;
            best_delta = delta;
            right = mid - 1;
        } else {
            left = mid + 1;
        }
    }
    
    return {best_n, best_delta};
}



// ====================== GENERATION OF NU VARIANTS ======================

vector<Variant> generate_nu_variants(const MaterialParams& params) {
    vector<Variant> variants;
    int n_levels = params.lgN_levels.size();
    double lgN_min = params.lgN_levels[0];
    double lgN_max = params.lgN_levels.back();
    double eps = 0.01;
    
    // 1. Uniform distribution
    vector<double> nu_uniform(n_levels, 1.0 / n_levels);
    variants.push_back({0.0, nu_uniform, "Uniform"});
    
    // 2. Right shift (power law)
    vector<double> p_vals = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 7.0, 10.0};
    for (double p : p_vals) {
        vector<double> weights(n_levels);
        double total = 0.0;
        for (int i = 0; i < n_levels; ++i) {
            double lgN_norm = (params.lgN_levels[i] - lgN_min) / (lgN_max - lgN_min);
            weights[i] = pow(lgN_norm + eps, p);
            total += weights[i];
        }
        vector<double> nu(n_levels);
        for (int i = 0; i < n_levels; ++i) {
            nu[i] = weights[i] / total;
        }
        string desc = "Right shift p=" + to_string(p);
        variants.push_back({p, nu, desc});
    }
    
    // 3. Left shift
    for (double p : p_vals) {
        vector<double> weights(n_levels);
        double total = 0.0;
        for (int i = 0; i < n_levels; ++i) {
            double lgN_norm = (lgN_max - params.lgN_levels[i]) / (lgN_max - lgN_min);
            weights[i] = pow(lgN_norm + eps, p);
            total += weights[i];
        }
        vector<double> nu(n_levels);
        for (int i = 0; i < n_levels; ++i) {
            nu[i] = weights[i] / total;
        }
        string desc = "Left shift p=" + to_string(p);
        variants.push_back({-p, nu, desc});
    }
    
    // 4. Symmetric configurations for 5 levels
    if (n_levels == 5) {
        vector<vector<double>> sym_configs = {
            {0.3, 0.1, 0.1, 0.1, 0.3}, {0.25, 0.15, 0.1, 0.15, 0.25},
            {0.2, 0.2, 0.2, 0.2, 0.2}, {0.15, 0.25, 0.2, 0.25, 0.15},
            {0.1, 0.3, 0.2, 0.3, 0.1}, {0.05, 0.35, 0.2, 0.35, 0.05}
        };
        vector<string> sym_desc = {
            "Symmetric (0.3,0.1,0.1,0.1,0.3)", "Symmetric (0.25,0.15,0.1,0.15,0.25)",
            "Symmetric (0.2,0.2,0.2,0.2,0.2)", "Symmetric (0.15,0.25,0.2,0.25,0.15)",
            "Symmetric (0.1,0.3,0.2,0.3,0.1)", "Symmetric (0.05,0.35,0.2,0.35,0.05)"
        };
        for (size_t j = 0; j < sym_configs.size(); ++j) {
            double total = 0;
            for (double v : sym_configs[j]) total += v;
            vector<double> nu(n_levels);
            for (int i = 0; i < n_levels; ++i) nu[i] = sym_configs[j][i] / total;
            variants.push_back({0.0, nu, sym_desc[j]});
        }
    }
    
    // 5. Extreme configurations for 5 levels
    if (n_levels == 5) {
        vector<vector<double>> ext_configs = {
            {0.05, 0.05, 0.1, 0.3, 0.5}, {0.5, 0.3, 0.1, 0.05, 0.05}
        };
        vector<string> ext_desc = {
            "Extreme right (0.05,0.05,0.1,0.3,0.5)",
            "Extreme left (0.5,0.3,0.1,0.05,0.05)"
        };
        for (size_t j = 0; j < ext_configs.size(); ++j) {
            double total = 0;
            for (double v : ext_configs[j]) total += v;
            vector<double> nu(n_levels);
            for (int i = 0; i < n_levels; ++i) nu[i] = ext_configs[j][i] / total;
            variants.push_back({0.0, nu, ext_desc[j]});
        }
    }
    
    return variants;
}

// ====================== STAGE 1: THEORETICAL OPTIMIZATION ======================

Plan theoretical_optimization(const MaterialParams& params, double N0, ostream& out_file) {
    cout << "\n  Theoretical optimization for N0=" << scientific << setprecision(0) << N0 << "..." << endl;
    out_file << "\n  Theoretical optimization for N0=" << scientific << setprecision(0) << N0 << "..." << endl;
    
    vector<Variant> variants = generate_nu_variants(params);
    cout << "    Generating distribution variants: " << variants.size() << " pcs." << endl;
    
    double delta_target = params.delta_target();
    
    struct Result {
        Variant var;
        int n;
        double delta;
        double cost;
        vector<int> n_dist;
    };
    vector<Result> results;
    
    int total_variants = variants.size();
    int current_var = 0;
    
    for (const auto& var : variants) {
        current_var++;
        if (current_var % 5 == 0 || current_var == total_variants) {
            cout << "\r    Processing variants: " << current_var << "/" << total_variants << flush;
        }
        
        auto [n_opt, delta_opt] = theoretical_min_n(params, var.nu, N0, 1, 300);
        double cost = calculate_cost(params, n_opt, var.nu);
        
        // Specimen distribution (round to integers)
        vector<int> n_dist(var.nu.size());
        int total = 0;
        for (size_t i = 0; i < var.nu.size(); ++i) {
            n_dist[i] = (int)round(n_opt * var.nu[i]);
            if (n_dist[i] < 0) n_dist[i] = 0;
            total += n_dist[i];
        }
        
        // Ensure at least 1 specimen per level with nu[i] > 0
        for (size_t i = 0; i < var.nu.size(); ++i) {
            if (var.nu[i] > 1e-6 && n_dist[i] == 0) {
                n_dist[i] = 1;
                total++;
            }
        }
        
        // Adjust total to match n_opt
        if (total != n_opt) {
            int idx = max_element(n_dist.begin(), n_dist.end()) - n_dist.begin();
            n_dist[idx] += n_opt - total;
            if (n_dist[idx] < 0) n_dist[idx] = 0;
        }
        
        results.push_back({var, n_opt, delta_opt, cost, n_dist});
    }
    cout << endl;
    
    out_file << "\n    " << left << setw(40) << "Variant" 
        << right << setw(4) << "n" 
        << setw(10) << "delta" 
        << setw(10) << "cost"
        << " nu1...nun" << endl;
    out_file << "    " << string(40, '-') << " " << string(4, '-') 
        << " " << string(10, '-') << " " << string(10, '-') << " " << string(50, '-') << endl;
    
    for (const auto& res : results) {
        string status = (res.delta <= delta_target && res.delta > 0) ? "OK" : "NO";
        out_file << "    " << left << setw(40) << res.var.description.substr(0, 38)
            << right << setw(4) << res.n
            << fixed << setw(10) << setprecision(5) << res.delta
            << setw(10) << setprecision(1) << res.cost << "   [";
        for (size_t i = 0; i < res.var.nu.size(); ++i) {
            out_file << fixed << setprecision(3) << res.var.nu[i];
            if (i < res.var.nu.size() - 1) out_file << ", ";
        }
        out_file << "] " << status << endl;
    }
    
    // Select best plan
    int best_idx = 0;
    vector<int> valid_indices;
    for (size_t i = 0; i < results.size(); ++i) {
        if (results[i].delta <= delta_target && results[i].delta > 0 && results[i].delta < 0.5) {
            valid_indices.push_back(i);
        }
    }
    
    if (!valid_indices.empty()) {
        best_idx = valid_indices[0];
        for (size_t idx : valid_indices) {
            if (results[idx].cost < results[best_idx].cost) {
                best_idx = idx;
            }
        }
        cout << "\n    ★ Best theoretical plan (accuracy achieved):" << endl;
        out_file << "\n    ★ Best theoretical plan (accuracy achieved):" << endl;
    } else {
        for (size_t i = 0; i < results.size(); ++i) {
            if (results[i].delta < results[best_idx].delta && results[i].delta > 0) {
                best_idx = i;
            }
        }
        cout << "\n    ★ Best theoretical plan (accuracy not achieved, minimal delta):" << endl;
        out_file << "\n    ★ Best theoretical plan (accuracy not achieved, minimal delta):" << endl;
    }
    
    Result& best = results[best_idx];
    
    cout << "      " << best.var.description << endl;
    cout << "      nu = [";
    for (size_t i = 0; i < best.var.nu.size(); ++i) {
        cout << fixed << setprecision(4) << best.var.nu[i];
        if (i < best.var.nu.size() - 1) cout << ", ";
    }
    cout << "]" << endl;
    cout << "      n = " << best.n << " -> distribution: [";
    for (size_t i = 0; i < best.n_dist.size(); ++i) {
        cout << best.n_dist[i];
        if (i < best.n_dist.size() - 1) cout << ", ";
    }
    cout << "]" << endl;
    cout << "      delta = " << fixed << setprecision(5) << best.delta 
        << " (target <= " << delta_target << ")" << endl;
    cout << "      cost = " << fixed << setprecision(1) << best.cost << endl;
    
    out_file << "      " << best.var.description << endl;
    out_file << "      nu = [";
    for (size_t i = 0; i < best.var.nu.size(); ++i) {
        out_file << fixed << setprecision(4) << best.var.nu[i];
        if (i < best.var.nu.size() - 1) out_file << ", ";
    }
    out_file << "]" << endl;
    out_file << "      n = " << best.n << " -> distribution: [";
    for (size_t i = 0; i < best.n_dist.size(); ++i) {
        out_file << best.n_dist[i];
        if (i < best.n_dist.size() - 1) out_file << ", ";
    }
    out_file << "]" << endl;
    out_file << "      delta = " << fixed << setprecision(5) << best.delta 
        << " (target <= " << delta_target << ")" << endl;
    out_file << "      cost = " << fixed << setprecision(1) << best.cost << endl;
    
    Plan plan;
    plan.nu = best.var.nu;
    plan.n = best.n;
    plan.delta = best.delta;
    plan.cost = best.cost;
    plan.description = best.var.description;
    plan.n_distribution = best.n_dist;
    return plan;
}

// ====================== STAGE 2: MONTE CARLO REFINEMENT ======================

class MonteCarloEstimator {
private:
    const MaterialParams& params;
    double N0;
    int n_simulations;
    mt19937 rng;
    
   

vector<vector<double>> generate_experiment(int n, const vector<double>& nu, int seed) {
    mt19937 local_rng(seed + 10000);
    normal_distribution<> normal(0.0, 1.0);
    
    // ТОЧНО КАК В test.cpp
    vector<int> n_i(nu.size());
    for (size_t i = 0; i < nu.size(); ++i) {
        n_i[i] = max(1, (int)round(n * nu[i]));
    }
    
    int total = 0;
    for (int v : n_i) total += v;
    if (total != n) {
        int idx = max_element(n_i.begin(), n_i.end()) - n_i.begin();
        n_i[idx] += n - total;
    }
    
    vector<vector<double>> data;
    for (size_t i = 0; i < params.lgN_levels.size(); ++i) {
        if (n_i[i] <= 0) continue;
        
        double lgN_mean = params.lgN_levels[i];
        double sigma_lgN = params.a * pow(lgN_mean, params.alpha);
        double sigma_level = params.sigma_from_lgN(lgN_mean);
        
        for (int j = 0; j < n_i[i]; ++j) {
            double lgN_sample = lgN_mean + sigma_lgN * normal(local_rng);
            data.push_back({sigma_level, lgN_sample});
        }
    }
    return data;
}
 





double estimate_strength_from_data(const vector<vector<double>>& data, double N0) {
    if (data.size() < 4) return -1.0;
    
    map<double, vector<double>> level_data;
    for (const auto& point : data) {
        level_data[point[0]].push_back(point[1]);
    }
    
    vector<double> sigma_vals, lgN_vals;
    vector<int> ni_vals;
    vector<double> slgN_vals;
    
    for (auto& kv : level_data) {
        double sigma = kv.first;
        vector<double>& samples = kv.second;
        int ni = samples.size();
        
        double mean = 0.0;
        for (double v : samples) mean += v;
        mean /= ni;
        
        double var = 0.0;
        for (double v : samples) var += (v - mean) * (v - mean);
        double std = sqrt(var / max(1, ni - 1));
        
        if (std < 0.01) std = 0.01;
        
        sigma_vals.push_back(sigma);
        lgN_vals.push_back(mean);
        ni_vals.push_back(ni);
        slgN_vals.push_back(std);
    }
    
    double sigma_inf_est = 250.0;
    double C_est = 800.0;
    double m_est = 0.12;
    double Q;
    
    int iterations = fatique_curve_fit(1, sigma_vals, ni_vals, lgN_vals, slgN_vals,
                                       sigma_inf_est, C_est, m_est, Q);
    
    if (iterations > 0 && sigma_inf_est > 0 && C_est > 0 && m_est > 0) {
        //sigma_inf_est = max(200.0, min(350.0, sigma_inf_est));
        //C_est = max(400.0, min(1600.0, C_est));
        //m_est = max(0.08, min(0.20, m_est));
        return sigma_inf_est + C_est * pow(N0, -m_est);
    }
    
    return -1.0;
}






    
public:
    MonteCarloEstimator(const MaterialParams& p, double n0, int sims)
        : params(p), N0(n0), n_simulations(sims) {
        random_device rd;
        rng.seed(rd());
    }
    
    double compute_delta(int n, const vector<double>& nu) {
        cout << "      Computing delta for n=" << n << "..." << endl;
        
        vector<double> estimates;
        estimates.reserve(n_simulations);
        
        for (int seed = 0; seed < n_simulations; ++seed) {
            auto data = generate_experiment(n, nu, seed);
            double est = estimate_strength_from_data(data, N0);
            
            if (est > 0) {
                estimates.push_back(est);
            }
            
            if ((seed + 1) % 100 == 0 || seed + 1 == n_simulations) {
                cout << "\r        Progress: " << (seed + 1) << "/" << n_simulations 
                     << ", valid estimates: " << estimates.size() << flush;
            }
        }
        cout << endl;
        
        if (estimates.size() < 20) {
            cout << "        WARNING: Only " << estimates.size() << " valid estimates (need 20)" << endl;
            return 0.25;
        }
        
        double mean = 0.0;
        for (double v : estimates) mean += v;
        mean /= estimates.size();
        
        double var = 0.0;
        for (double v : estimates) var += (v - mean) * (v - mean);
        var /= estimates.size();
        
        double delta = sqrt(var) / mean;
        cout << "        mean = " << mean << ", std = " << sqrt(var) << ", delta = " << delta << endl;
        
        return delta;
    }
};




Plan monte_carlo_optimization(const MaterialParams& params, const Plan& theoretical_plan,
                               double N0, ostream& out_file) {
    cout << "\n  Stage 2: Monte Carlo refinement (simulations=" << params.n_simulations << ")..." << endl;
    out_file << "\n  Stage 2: Monte Carlo refinement (simulations=" << params.n_simulations << ")..." << endl;
    
    MonteCarloEstimator mc(params, N0, params.n_simulations);
    double delta_target = params.delta_target();
    
    cout << "    Estimating accuracy of theoretical plan..." << endl;
    int n_current = theoretical_plan.n;
    double delta_current = mc.compute_delta(n_current, theoretical_plan.nu);
    cout << "    Theoretical plan: n=" << n_current << ", delta=" << fixed << setprecision(5) << delta_current << endl;
    out_file << "    Theoretical plan: n=" << n_current << ", delta=" << fixed << setprecision(5) << delta_current << endl;
    
    int n_min = 4;
    int n_max = 500;
    
    // БИНАРНЫЙ ПОИСК
    if (delta_current > delta_target) {
        cout << "    Accuracy not achieved (delta=" << delta_current << "), increasing n..." << endl;
        out_file << "    Accuracy not achieved, increasing n..." << endl;
        
        int left = n_current;
        int right = n_max;
        int best_n = n_max;
        double best_delta = mc.compute_delta(n_max, theoretical_plan.nu);
        
        // Бинарный поиск минимального n, при котором delta <= delta_target
        while (left <= right) {
            int mid = (left + right) / 2;
            cout << "    Checking n=" << mid << "..." << endl;
            double delta = mc.compute_delta(mid, theoretical_plan.nu);
            
            if (delta <= delta_target && delta > 0) {
                best_n = mid;
                best_delta = delta;
                right = mid - 1;  // пробуем меньше
                cout << "      OK (delta=" << delta << "), trying smaller" << endl;
            } else {
                left = mid + 1;   // нужно больше
                cout << "      NOT OK (delta=" << delta << "), need larger" << endl;
            }
        }
        
        n_current = best_n;
        delta_current = best_delta;
        cout << "    Found n = " << n_current << ", delta = " << delta_current << endl;
        out_file << "    Found n = " << n_current << ", delta = " << delta_current << endl;
        
        if (delta_current > delta_target) {
            cout << "    WARNING: Even at n=" << n_max << " accuracy not achieved!" << endl;
            out_file << "    WARNING: Even at n=" << n_max << " accuracy not achieved!" << endl;
        }
        
    } else if (delta_current < delta_target * 0.99) {
        cout << "    Accuracy excessive (delta=" << delta_current << "), trying to decrease n..." << endl;
        out_file << "    Accuracy excessive, trying to decrease n..." << endl;
        
        int left = n_min;
        int right = n_current;
        int best_n = n_min;
        double best_delta = mc.compute_delta(n_min, theoretical_plan.nu);
        
        // Бинарный поиск минимального n, при котором delta <= delta_target
        while (left <= right) {
            int mid = (left + right) / 2;
            cout << "    Checking n=" << mid << "..." << endl;
            double delta = mc.compute_delta(mid, theoretical_plan.nu);
            
            if (delta <= delta_target && delta > 0) {
                best_n = mid;
                best_delta = delta;
                right = mid - 1;  // пробуем ещё меньше
                cout << "      OK (delta=" << delta << "), trying smaller" << endl;
            } else {
                left = mid + 1;   // нужно больше
                cout << "      NOT OK (delta=" << delta << "), need larger" << endl;
            }
        }
        
        n_current = best_n;
        delta_current = best_delta;
        cout << "    Reduced to n = " << n_current << ", delta = " << delta_current << endl;
        out_file << "    Reduced to n = " << n_current << ", delta = " << delta_current << endl;
    }
    
    // Финальное распределение
    vector<int> n_dist(theoretical_plan.nu.size());
    int total = 0;
    for (size_t i = 0; i < theoretical_plan.nu.size(); ++i) {
        n_dist[i] = (int)round(n_current * theoretical_plan.nu[i]);
        if (n_dist[i] < 0) n_dist[i] = 0;
        total += n_dist[i];
    }
    
    for (size_t i = 0; i < theoretical_plan.nu.size(); ++i) {
        if (theoretical_plan.nu[i] > 1e-6 && n_dist[i] == 0) {
            n_dist[i] = 1;
            total++;
        }
    }
    
    if (total != n_current && total > 0) {
        int idx = max_element(n_dist.begin(), n_dist.end()) - n_dist.begin();
        n_dist[idx] += n_current - total;
        if (n_dist[idx] < 0) n_dist[idx] = 0;
    }
    
    double cost = calculate_cost(params, n_current, theoretical_plan.nu);
    
    cout << "\n    ★ Final plan (Monte Carlo):" << endl;
    cout << "      nu = [";
    for (size_t i = 0; i < theoretical_plan.nu.size(); ++i) {
        cout << fixed << setprecision(4) << theoretical_plan.nu[i];
        if (i < theoretical_plan.nu.size() - 1) cout << ", ";
    }
    cout << "]" << endl;
    cout << "      n = " << n_current << " -> distribution: [";
    for (size_t i = 0; i < n_dist.size(); ++i) {
        cout << n_dist[i];
        if (i < n_dist.size() - 1) cout << ", ";
    }
    cout << "]" << endl;
    cout << "      delta = " << fixed << setprecision(5) << delta_current 
        << " (target <= " << delta_target << ")" << endl;
    cout << "      cost = " << fixed << setprecision(1) << cost << endl;
    
    out_file << "\n    ★ Final plan (Monte Carlo):" << endl;
    out_file << "      nu = [";
    for (size_t i = 0; i < theoretical_plan.nu.size(); ++i) {
        out_file << fixed << setprecision(4) << theoretical_plan.nu[i];
        if (i < theoretical_plan.nu.size() - 1) out_file << ", ";
    }
    out_file << "]" << endl;
    out_file << "      n = " << n_current << " -> distribution: [";
    for (size_t i = 0; i < n_dist.size(); ++i) {
        out_file << n_dist[i];
        if (i < n_dist.size() - 1) out_file << ", ";
    }
    out_file << "]" << endl;
    out_file << "      delta = " << fixed << setprecision(5) << delta_current 
        << " (target <= " << delta_target << ")" << endl;
    out_file << "      cost = " << fixed << setprecision(1) << cost << endl;
    
    Plan final_plan;
    final_plan.nu = theoretical_plan.nu;
    final_plan.n = n_current;
    final_plan.delta = delta_current;
    final_plan.cost = cost;
    final_plan.description = theoretical_plan.description;
    final_plan.n_distribution = n_dist;
    return final_plan;
}












// ====================== MAIN FUNCTION ======================

int main() {
    cout << "======================================================================" << endl;
    cout << "TWO-STAGE FATIGUE TEST PLANNING OPTIMIZATION" << endl;
    cout << "Stage 1: Theoretical optimization (based on Fisher information matrix)" << endl;
    cout << "Stage 2: Monte Carlo refinement" << endl;
    cout << "======================================================================" << endl;
    cout << endl;
    
    // Redirect output to file
    ofstream out_file("optimization_result.out");
    if (!out_file.is_open()) {
        cerr << "Failed to create output file" << endl;
        return 1;
    }
    
    // Read parameters
    cout << "Reading parameters from input.inp..." << endl;
    MaterialParams params;
    if (!read_input("input.inp", params)) {
        out_file << "ERROR: Failed to read input.inp file" << endl;
        cerr << "ERROR: Failed to read input.inp file" << endl;
        return 1;
    }
    cout << "Parameters loaded." << endl << endl;
    
    out_file << "======================================================================" << endl;
    out_file << "TWO-STAGE FATIGUE TEST PLANNING OPTIMIZATION" << endl;
    out_file << "Stage 1: Theoretical optimization (based on Fisher information matrix)" << endl;
    out_file << "Stage 2: Monte Carlo refinement" << endl;
    out_file << "======================================================================" << endl;
    
    out_file << "\nProblem parameters:" << endl;
    out_file << "  sigma_inf = " << params.sigma_inf << " MPa" << endl;
    out_file << "  C = " << params.C << ", m = " << params.m << endl;
    out_file << "  a = " << params.a << ", alpha = " << params.alpha << endl;
    out_file << "  lgN levels (" << params.lgN_levels.size() << " pcs): ";
    for (double v : params.lgN_levels) out_file << v << " ";
    out_file << endl;
    out_file << "  C3/C1 = " << params.C3 / params.C1 << endl;
    out_file << "  eps = " << params.eps << ", gamma = " << params.gamma 
             << ", t = " << params.t_quantile << endl;
    out_file << "  delta_target = " << fixed << setprecision(5) << params.delta_target() << endl;
    out_file << "  f = " << params.f << " Hz" << endl;
    out_file << "  n_simulations = " << params.n_simulations << endl;
    
    map<string, pair<Plan, Plan>> all_results;
    
    int n0_count = params.N0_list.size();
    int current_n0 = 0;
    
    for (double N0 : params.N0_list) {
        current_n0++;
        cout << "\n============================================================" << endl;
        cout << "OPTIMIZATION FOR N0 = " << scientific << setprecision(0) << N0;
        cout << " [" << current_n0 << "/" << n0_count << "]" << endl;
        cout << "============================================================" << endl;
        
        out_file << "\n======================================================================" << endl;
        out_file << "OPTIMIZATION FOR N0 = " << scientific << setprecision(0) << N0 << endl;
        out_file << "======================================================================" << endl;
        
        auto start = chrono::high_resolution_clock::now();
        
        cout << "\n[1/2] Theoretical optimization..." << endl;
        Plan theoretical = theoretical_optimization(params, N0, out_file);
        cout << "      [OK] Theoretical optimization completed" << endl;
        
        cout << "\n[2/2] Monte Carlo refinement (simulations=" << params.n_simulations << ")..." << endl;
        Plan final_plan = monte_carlo_optimization(params, theoretical, N0, out_file);
        cout << "      [OK] Monte Carlo refinement completed" << endl;
        
        auto end = chrono::high_resolution_clock::now();
        double elapsed = chrono::duration<double>(end - start).count();
        cout << "\n      Execution time: " << fixed << setprecision(1) << elapsed << " sec" << endl;
        
        out_file << "\n    Execution time: " << elapsed << " sec" << endl;
        
        string key = "N0_" + to_string((long long)N0);
        all_results[key] = {theoretical, final_plan};
    }
    
    // Print summary table
    cout << "\n\n======================================================================" << endl;
    cout << "SUMMARY TABLE OF RESULTS" << endl;
    cout << "======================================================================" << endl;
    cout << "\n" << left << setw(15) << "N0"
         << setw(15) << "Stage"
         << setw(6) << "n"
         << setw(10) << "delta"
         << setw(10) << "cost"
         << "Distribution" << endl;
    cout << string(80, '-') << endl;
    
    out_file << "\n======================================================================" << endl;
    out_file << "SUMMARY TABLE OF RESULTS" << endl;
    out_file << "======================================================================" << endl;
    out_file << "\n" << left << setw(15) << "N0"
             << setw(15) << "Stage"
             << setw(6) << "n"
             << setw(10) << "delta"
             << setw(10) << "cost"
             << "Distribution (n1...nn)" << endl;
    out_file << string(80, '-') << endl;
    
    for (const auto& [key, plans] : all_results) {
        const Plan& theo = plans.first;
        const Plan& mc = plans.second;
        
        string theo_dist = "";
        for (size_t i = 0; i < theo.n_distribution.size(); ++i) {
            theo_dist += to_string(theo.n_distribution[i]);
            if (i < theo.n_distribution.size() - 1) theo_dist += ", ";
        }
        
        string mc_dist = "";
        for (size_t i = 0; i < mc.n_distribution.size(); ++i) {
            mc_dist += to_string(mc.n_distribution[i]);
            if (i < mc.n_distribution.size() - 1) mc_dist += ", ";
        }
        
        // Console output
        cout << left << setw(15) << key
             << setw(15) << "Theoretical"
             << right << setw(6) << theo.n
             << fixed << setw(10) << setprecision(5) << theo.delta
             << setw(10) << setprecision(1) << theo.cost
             << " [" << theo_dist << "]" << endl;
        
        cout << left << setw(15) << ""
             << setw(15) << "Monte Carlo"
             << right << setw(6) << mc.n
             << fixed << setw(10) << setprecision(5) << mc.delta
             << setw(10) << setprecision(1) << mc.cost
             << " [" << mc_dist << "]" << endl;
        cout << endl;
        
        // File output
        out_file << left << setw(15) << key
                 << setw(15) << "Theoretical"
                 << right << setw(6) << theo.n
                 << fixed << setw(10) << setprecision(5) << theo.delta
                 << setw(10) << setprecision(1) << theo.cost
                 << " [" << theo_dist << "]" << endl;
        
        out_file << left << setw(15) << ""
                 << setw(15) << "Monte Carlo"
                 << right << setw(6) << mc.n
                 << fixed << setw(10) << setprecision(5) << mc.delta
                 << setw(10) << setprecision(1) << mc.cost
                 << " [" << mc_dist << "]" << endl;
        out_file << endl;
    }
    
    cout << "\n======================================================================" << endl;
    cout << "OPTIMIZATION COMPLETED" << endl;
    cout << "Results saved to file: optimization_result.out" << endl;
    cout << "======================================================================" << endl;
    
    out_file << "\n======================================================================" << endl;
    out_file << "RESULTS SAVED TO FILE: optimization_result.out" << endl;
    out_file << "======================================================================" << endl;
    
    out_file.close();
    
    return 0;
}