dA.cpp

#include "dA.h"
#include <algorithm>
#include <cmath>
#include <iostream>
#include <numeric>
#include <random>

// dA_params实现
dA_params::dA_params(size_t n_visible, size_t n_hidden, double lr,
                     double corruption_level, size_t gracePeriod,
                     double hiddenRatio)
    : n_visible(n_visible), n_hidden(n_hidden), lr(lr),
      corruption_level(corruption_level), gracePeriod(gracePeriod),
      hiddenRatio(hiddenRatio) {}

// dA实现
dA::dA(const dA_params &params) : params(params), n(0), rng(1234) {
    // 如果指定了hiddenRatio，根据可见层计算隐藏层大小
    if (params.hiddenRatio > 0.0) {
        this->params.n_hidden = static_cast<size_t>(
            std::ceil(params.n_visible * params.hiddenRatio));
    }
    // 确保n_hidden至少为1，否则后续操作无法进行
    if (this->params.n_hidden == 0) {
        this->params.n_hidden = 1;
    }

    // 初始化归一化边界
    norm_max.resize(params.n_visible, -std::numeric_limits<double>::infinity());
    norm_min.resize(params.n_visible, std::numeric_limits<double>::infinity());

    // 初始化权重矩阵 W (均匀分布)
    double a = 1.0 / params.n_visible;
    W.resize(params.n_visible, std::vector<double>(this->params.n_hidden));

    std::uniform_real_distribution<double> dist(-a, a);
    for (size_t i = 0; i < params.n_visible; ++i) {
        for (size_t j = 0; j < this->params.n_hidden; ++j) {
            W[i][j] = dist(rng);
        }
    }

    // 初始化偏置向量 (全0)
    hbias.resize(this->params.n_hidden, 0.0);
    vbias.resize(params.n_visible, 0.0);

    // 计算W的转置
    W_prime.resize(this->params.n_hidden,
                   std::vector<double>(params.n_visible));
    for (size_t i = 0; i < params.n_visible; ++i) {
        for (size_t j = 0; j < this->params.n_hidden; ++j) {
            W_prime[j][i] = W[i][j];
        }
    }
}

std::vector<double> dA::get_corrupted_input(const std::vector<double> &input,
                                            double corruption_level) {
    assert(corruption_level < 1.0);

    std::vector<double> corrupted = input;
    std::bernoulli_distribution dist(1.0 - corruption_level);

    for (auto &val : corrupted) {
        val *= dist(rng) ? 1.0 : 0.0;
    }

    return corrupted;
}

std::vector<double> dA::get_hidden_values(const std::vector<double> &input) {
    // 检查n_hidden是否为0，如果为0，返回空向量
    if (params.n_hidden == 0) {
        return std::vector<double>();
    }

    std::vector<double> hidden(params.n_hidden, 0.0);

    // 计算 input · W + hbias
    for (size_t j = 0; j < params.n_hidden; ++j) {
        for (size_t i = 0; i < params.n_visible; ++i) {
            hidden[j] += input[i] * W[i][j];
        }
        hidden[j] += hbias[j];
    }

    // 应用sigmoid激活函数
    for (auto &val : hidden) {
        val = utils::sigmoid(val);
    }

    return hidden;
}

std::vector<double>
dA::get_reconstructed_input(const std::vector<double> &hidden) {
    // 检查n_hidden是否为0，如果为0，返回与输入大小相同的零向量
    if (params.n_hidden == 0) {
        return std::vector<double>(params.n_visible, 0.0);
    }

    std::vector<double> reconstructed(params.n_visible, 0.0);

    // 计算 hidden · W' + vbias
    for (size_t i = 0; i < params.n_visible; ++i) {
        for (size_t j = 0; j < params.n_hidden; ++j) {
            reconstructed[i] += hidden[j] * W_prime[j][i];
        }
        reconstructed[i] += vbias[i];
    }

    // 应用sigmoid激活函数
    for (auto &val : reconstructed) {
        val = utils::sigmoid(val);
    }

    return reconstructed;
}

double dA::train(const std::vector<double> &x) {
    n++;

    // 更新归一化边界
    for (size_t i = 0; i < params.n_visible; ++i) {
        if (x[i] > norm_max[i]) {
            norm_max[i] = x[i];
        }
        if (x[i] < norm_min[i]) {
            norm_min[i] = x[i];
        }
    }

    // 0-1归一化
    std::vector<double> x_normalized(params.n_visible);
    for (size_t i = 0; i < params.n_visible; ++i) {
        x_normalized[i] =
            (x[i] - norm_min[i]) / (norm_max[i] - norm_min[i] + 1e-16);
    }

    // 添加噪声 (如果需要)
    std::vector<double> tilde_x;
    if (params.corruption_level > 0.0) {
        tilde_x = get_corrupted_input(x_normalized, params.corruption_level);
    } else {
        tilde_x = x_normalized;
    }

    // 确保n_hidden大于0
    if (params.n_hidden == 0) {
        // 如果n_hidden为0，返回一个默认误差值
        return 0.0;
    }

    // 前向传播
    std::vector<double> y = get_hidden_values(tilde_x);
    std::vector<double> z = get_reconstructed_input(y);

    // 计算误差
    std::vector<double> L_h2(params.n_visible);
    for (size_t i = 0; i < params.n_visible; ++i) {
        L_h2[i] = x_normalized[i] - z[i];
    }

    // 计算隐藏层梯度
    std::vector<double> L_h1(params.n_hidden, 0.0);
    for (size_t j = 0; j < params.n_hidden; ++j) {
        for (size_t i = 0; i < params.n_visible; ++i) {
            L_h1[j] += L_h2[i] * W[i][j];
        }
        L_h1[j] *= y[j] * (1.0 - y[j]);
    }

    // 更新权重和偏置
    for (size_t i = 0; i < params.n_visible; ++i) {
        for (size_t j = 0; j < params.n_hidden; ++j) {
            W[i][j] += params.lr * (tilde_x[i] * L_h1[j] + L_h2[i] * y[j]);
            // 同时更新W_prime (保持W'是W的转置)
            W_prime[j][i] = W[i][j];
        }
        vbias[i] += params.lr * L_h2[i];
    }

    for (size_t j = 0; j < params.n_hidden; ++j) {
        hbias[j] += params.lr * L_h1[j];
    }

    // 计算并返回RMSE
    double sum_squared_error = 0.0;
    for (size_t i = 0; i < params.n_visible; ++i) {
        sum_squared_error += L_h2[i] * L_h2[i];
    }

    return std::sqrt(sum_squared_error / params.n_visible);
}

std::vector<double> dA::reconstruct(const std::vector<double> &x) {
    // 确保n_hidden大于0
    if (params.n_hidden == 0) {
        // 如果n_hidden为0，返回与输入大小相同的零向量
        return std::vector<double>(params.n_visible, 0.0);
    }

    return get_reconstructed_input(get_hidden_values(x));
}

double dA::execute(const std::vector<double> &x) {
    if (inGrace()) {
        return 0.0;
    }

    // 确保n_hidden大于0
    if (params.n_hidden == 0) {
        // 如果n_hidden为0，返回一个默认误差值
        return 0.0;
    }

    // 0-1归一化
    std::vector<double> x_normalized(params.n_visible);
    for (size_t i = 0; i < params.n_visible; ++i) {
        x_normalized[i] =
            (x[i] - norm_min[i]) / (norm_max[i] - norm_min[i] + 1e-16);
    }

    // 重构输入
    std::vector<double> z = reconstruct(x_normalized);

    // 计算MSE
    double sum_squared_error = 0.0;
    for (size_t i = 0; i < params.n_visible; ++i) {
        double diff = x_normalized[i] - z[i];
        sum_squared_error += diff * diff;
    }

    return std::sqrt(sum_squared_error / params.n_visible);
}

bool dA::inGrace() const { return n < params.gracePeriod; }