gpu_postprocessor.cu - master - nik2703/boundary_approx

gpu_postprocessor.cu
309 строк11 KB
nik2703
upload files
08 окт 2025, 19:39
08 окт 2025, 19:39ccf6b5d
// gpu_postprocessor.cu
// Финальная реализация Шага 2.3: Постобработка и валидация на GPU (БЕЗ THRUST)
#include "gpu_postprocessor.h"
#include "gpu_segmentation.h"
#include <cuda_runtime.h>
#include <cmath>
#include <iostream>
#include <vector>
#include <algorithm>
#include <functional>

// ==========================
// 1. ВСПОМОГАТЕЛЬНЫЕ ФУНКЦИИ (CPU)
// ==========================
void generateMonomialPowers(std::vector<std::vector<int>>& powers, int num_vars, int max_degree) {
    powers.clear();
    std::vector<int> current(num_vars, 0);
    std::function<void(int, int)> generate = [&](int var_idx, int remaining_degree) {
        if (var_idx == num_vars - 1) {
            current[var_idx] = remaining_degree;
            powers.push_back(current);
            return;
        }
        for (int d = 0; d <= remaining_degree; ++d) {
            current[var_idx] = d;
            generate(var_idx + 1, remaining_degree - d);
        }
    };
    for (int total_deg = 0; total_deg <= max_degree; ++total_deg) {
        generate(0, total_deg);
    }
}

// ==========================
// 2. CUDA ЯДРА
// ==========================
__device__ float evaluatePolynomial(
    const float* x_prime,
    const int* monomial_powers,
    const float* coefficients,
    int proj_dim,
    int n_basis,
    int point_idx
) {
    float result = 0.0f;
    const float* x = x_prime + point_idx * proj_dim;
    for (int b = 0; b < n_basis; ++b) {
        float monomial_val = 1.0f;
        const int* powers = monomial_powers + b * proj_dim;
        for (int d = 0; d < proj_dim; ++d) {
            int exp = powers[d];
            if (exp == 0) continue;
            float base = x[d];
            float power_result = 1.0f;
            while (exp > 0) {
                if (exp & 1) power_result *= base;
                base *= base;
                exp >>= 1;
            }
            monomial_val *= power_result;
        }
        result += coefficients[b] * monomial_val;
    }
    return result;
}

__global__ void computeResidualsKernel(
    const float* x_prime,
    const float* y_values,
    const int* monomial_powers,
    const float* coefficients,
    int proj_dim,
    int n_basis,
    int n,
    float* residuals
) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= n) return;

    float y_pred = evaluatePolynomial(x_prime, monomial_powers, coefficients, proj_dim, n_basis, idx);
    float y_true = y_values[idx];
    float residual = y_true - y_pred;
    residuals[idx] = residual * residual;
}

__global__ void computeGradientNormsKernel(
    const float* x_prime,
    const int* monomial_powers,
    const float* coefficients,
    int proj_dim,
    int n_basis,
    int n,
    float* grad_norms_sq
) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= n) return;

    const float* x = x_prime + idx * proj_dim;
    float norm_sq = 0.0f;

    for (int var_idx = 0; var_idx < proj_dim; ++var_idx) {
        float partial = 0.0f;
        for (int b = 0; b < n_basis; ++b) {
            const int* powers = monomial_powers + b * proj_dim;
            int exp = powers[var_idx];
            if (exp == 0) continue;
            float monomial_val = exp * coefficients[b];
            for (int d = 0; d < proj_dim; ++d) {
                int p = (d == var_idx) ? (exp - 1) : powers[d];
                if (p == 0) continue;
                float base = x[d];
                float power_result = 1.0f;
                while (p > 0) {
                    if (p & 1) power_result *= base;
                    base *= base;
                    p >>= 1;
                }
                monomial_val *= power_result;
            }
            partial += monomial_val;
        }
        norm_sq += partial * partial;
    }
    grad_norms_sq[idx] = norm_sq;
}

__global__ void reduceSumKernel(float* data, int n) {
    extern __shared__ float sdata[];
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;

    sdata[tid] = (i < n) ? data[i] : 0.0f;
    __syncthreads();

    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }

    if (tid == 0) data[blockIdx.x] = sdata[0];
}

__global__ void reduceMaxKernel(float* data, int n) {
    extern __shared__ float sdata[];
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;

    sdata[tid] = (i < n) ? data[i] : 0.0f;
    __syncthreads();

    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] = fmaxf(sdata[tid], sdata[tid + s]);
        }
        __syncthreads();
    }

    if (tid == 0) data[blockIdx.x] = sdata[0];
}

// ==========================
// 3. РЕАЛИЗАЦИЯ МЕТОДА КЛАССА
// ==========================
GPUPostProcessor::ValidationResult GPUPostProcessor::validateApproximation(
    const GPUSegment& segment,
    const ApproximationResult& result,
    int poly_degree
) {
    ValidationResult validation;
    validation.is_valid = true;
    validation.rss = 0.0f;
    validation.max_gradient_norm = 0.0f;
    validation.warning_message = "";

    if (segment.points.empty() || result.coefficients.empty()) {
        validation.is_valid = false;
        validation.warning_message = "Пустой сегмент или отсутствуют коэффициенты.";
        return validation;
    }

    size_t n = segment.points.size();
    size_t full_dim = segment.points[0].coords.size();
    size_t proj_dim = full_dim - 1;
    int n_basis = result.coefficients.size();

    // 1. Подготовка данных на GPU
    float *d_x_prime, *d_y_values, *d_monomial_powers, *d_coefficients;
    float *d_residuals, *d_grad_norms_sq;

    // Host data preparation
    std::vector<float> h_x_prime(n * proj_dim);
    std::vector<float> h_y_values(n);
    for (size_t i = 0; i < n; ++i) {
        const auto& p = segment.points[i];
        int k = 0;
        for (size_t d = 0; d < full_dim; ++d) {
            if (static_cast<int>(d) != segment.approx_axis) {
                h_x_prime[i * proj_dim + k] = p.coords[d];
                k++;
            }
        }
        h_y_values[i] = p.coords[segment.approx_axis];
    }

    std::vector<std::vector<int>> monomial_powers_vec;
    generateMonomialPowers(monomial_powers_vec, proj_dim, poly_degree);
    std::vector<float> h_monomial_powers(n_basis * proj_dim);
    for (int b = 0; b < n_basis; ++b) {
        for (int d = 0; d < proj_dim; ++d) {
            h_monomial_powers[b * proj_dim + d] = static_cast<float>(monomial_powers_vec[b][d]);
        }
    }

    // Device memory allocation
    cudaMalloc(&d_x_prime, n * proj_dim * sizeof(float));
    cudaMalloc(&d_y_values, n * sizeof(float));
    cudaMalloc(&d_monomial_powers, n_basis * proj_dim * sizeof(float));
    cudaMalloc(&d_coefficients, n_basis * sizeof(float));
    cudaMalloc(&d_residuals, n * sizeof(float));
    cudaMalloc(&d_grad_norms_sq, n * sizeof(float));

    // Copy data to device
    cudaMemcpy(d_x_prime, h_x_prime.data(), n * proj_dim * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_y_values, h_y_values.data(), n * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_monomial_powers, h_monomial_powers.data(), n_basis * proj_dim * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_coefficients, result.coefficients.data(), n_basis * sizeof(float), cudaMemcpyHostToDevice);

    // 2. Launch kernels
    int blockSize = 256;
    int gridSize = (n + blockSize - 1) / blockSize;

    computeResidualsKernel<<<gridSize, blockSize>>>(
        d_x_prime, d_y_values, (int*)d_monomial_powers, d_coefficients,
        proj_dim, n_basis, n, d_residuals
    );

    computeGradientNormsKernel<<<gridSize, blockSize>>>(
        d_x_prime, (int*)d_monomial_powers, d_coefficients,
        proj_dim, n_basis, n, d_grad_norms_sq
    );

    cudaDeviceSynchronize();

    // 3. Reduce RSS
    float total_rss = 0.0f;
    if (n > 1) {
        float *d_temp = d_residuals;
        int current_n = n;
        while (current_n > 1) {
            int grid = (current_n + blockSize - 1) / blockSize;
            reduceSumKernel<<<grid, blockSize, blockSize * sizeof(float)>>>(d_temp, current_n);
            current_n = grid;
            d_temp = d_residuals;
        }
        cudaMemcpy(&total_rss, d_residuals, sizeof(float), cudaMemcpyDeviceToHost);
    } else {
        cudaMemcpy(&total_rss, d_residuals, sizeof(float), cudaMemcpyDeviceToHost);
    }

    // 4. Reduce Max Gradient Norm
    float max_grad_norm_sq = 0.0f;
    if (n > 1) {
        float *d_temp = d_grad_norms_sq;
        int current_n = n;
        while (current_n > 1) {
            int grid = (current_n + blockSize - 1) / blockSize;
            reduceMaxKernel<<<grid, blockSize, blockSize * sizeof(float)>>>(d_temp, current_n);
            current_n = grid;
            d_temp = d_grad_norms_sq;
        }
        cudaMemcpy(&max_grad_norm_sq, d_grad_norms_sq, sizeof(float), cudaMemcpyDeviceToHost);
    } else {
        cudaMemcpy(&max_grad_norm_sq, d_grad_norms_sq, sizeof(float), cudaMemcpyDeviceToHost);
    }

    validation.rss = total_rss;
    validation.max_gradient_norm = sqrtf(max_grad_norm_sq);

    // 5. Validation
    float avg_rss_per_point = total_rss / n;
    if (avg_rss_per_point > RSS_HIGH_THRESHOLD) {
        validation.warning_message += "Высокая остаточная ошибка. ";
    }
    if (avg_rss_per_point < RSS_LOW_THRESHOLD) {
        validation.warning_message += "Возможно переобучение. ";
    }
    if (validation.max_gradient_norm > GRADIENT_THRESHOLD) {
        validation.warning_message += "Обнаружен высокий градиент (возможны ложные экстремумы). ";
    }

    std::cout << "Валидация сегмента (" << segment.cluster_i << ", " << segment.cluster_j << "):\n"
              << "  RSS: " << validation.rss << "\n"
              << "  Max Gradient Norm: " << validation.max_gradient_norm << "\n";
    if (!validation.warning_message.empty()) {
        std::cout << "  Предупреждения: " << validation.warning_message << "\n";
    }

    // Cleanup
    cudaFree(d_x_prime);
    cudaFree(d_y_values);
    cudaFree(d_monomial_powers);
    cudaFree(d_coefficients);
    cudaFree(d_residuals);
    cudaFree(d_grad_norms_sq);

    return validation;
}