test_core_correctness.cpp

#include <iostream>
#include <cassert>
#include <cmath>
#include "allheader.h"
#include "network.h"

using namespace ml;
using namespace std;

// Tolerance for floating point comparisons
const double EPSILON = 1e-6;

bool approxEqual(double a, double b, double epsilon = EPSILON) {
    return std::abs(a - b) < epsilon;
}

void printMatrix(const Mat<double>& mat, const string& name) {
    cout << "\n" << name << " (" << mat.size().cy << "x" << mat.size().cx << "):" << endl;
    for (int i = 0; i < mat.size().cy; i++) {
        cout << "  [ ";
        for (int j = 0; j < mat.size().cx; j++) {
            cout << mat.getAt(i, j) << " ";
        }
        cout << "]" << endl;
    }
}

// TEST 1: Matrix multiplication with known values
bool test_matrix_multiplication() {
    cout << "\n=== TEST 1: Matrix Multiplication ===" << endl;

    // Test case: [1, 2] * [[3], [4]] should equal [11]
    //            (1x2)  *  (2x1)   =   (1x1)
    Mat<double> A(1, 2, 0);
    A.setAt(0, 0, 1.0);
    A.setAt(0, 1, 2.0);

    Mat<double> B(2, 1, 0);
    B.setAt(0, 0, 3.0);
    B.setAt(1, 0, 4.0);

    printMatrix(A, "A (input)");
    printMatrix(B, "B (weights)");

    // Hypothesis: A * B should give us [1*3 + 2*4] = [11]
    // Use false to let Mult transpose B from (2,1) to (1,2) for us
    Mat<double> result = Mult<double>(A, B, false);

    printMatrix(result, "A * B (result)");

    double expected = 1.0 * 3.0 + 2.0 * 4.0;  // = 11.0
    double actual = result.getAt(0, 0);

    cout << "Expected: " << expected << ", Actual: " << actual << endl;

    if (!approxEqual(actual, expected)) {
        cout << "FAIL: Matrix multiplication incorrect!" << endl;
        return false;
    }

    cout << "PASS: Matrix multiplication correct" << endl;
    return true;
}

// TEST 2: Matrix multiplication with larger matrices
bool test_matrix_multiplication_2x3() {
    cout << "\n=== TEST 2: Matrix Multiplication 2x3 ===" << endl;

    // Test: [1, 2, 3] * [[1], [2], [3]] = [1*1 + 2*2 + 3*3] = [14]
    //       (1x3)     *  (3x1)         =   (1x1)
    Mat<double> A(1, 3, 0);
    A.setAt(0, 0, 1.0);
    A.setAt(0, 1, 2.0);
    A.setAt(0, 2, 3.0);

    Mat<double> B(3, 1, 0);
    B.setAt(0, 0, 1.0);
    B.setAt(1, 0, 2.0);
    B.setAt(2, 0, 3.0);

    printMatrix(A, "A (1x3)");
    printMatrix(B, "B (3x1)");

    // Use false to let Mult transpose B for us
    Mat<double> result = Mult<double>(A, B, false);
    printMatrix(result, "Result");

    double expected = 1.0*1.0 + 2.0*2.0 + 3.0*3.0;  // = 14.0
    double actual = result.getAt(0, 0);

    cout << "Expected: " << expected << ", Actual: " << actual << endl;

    if (!approxEqual(actual, expected)) {
        cout << "FAIL: Matrix multiplication 2x3 incorrect!" << endl;
        return false;
    }

    cout << "PASS: Matrix multiplication 2x3 correct" << endl;
    return true;
}

// TEST 3: Forward propagation with known weights
bool test_forward_propagation_simple() {
    cout << "\n=== TEST 3: Forward Propagation Simple ===" << endl;

    // Create a simple 2-input, 2-hidden, 1-output network
    // We'll set specific weights to verify correctness
    Layer<double>* input = new Layer<double>(2, "input");
    Layer<double>* hidden = new Layer<double>(2, "hidden");
    Layer<double>* output = new Layer<double>(1, "output");

    Network<double>* net = new Network<double>();
    net->setInputLayer(input);
    net->setOutputLayer(output);

    net->connect(input, hidden);
    net->connect(hidden, output);
    net->init();

    // Set specific weights for input->hidden
    // Input layer has 2 inputs + 1 bias = 3 outputs
    // Hidden layer has 2 nodes
    // So weights should be (2, 3) - 2 rows (one per hidden node), 3 cols (one per input+bias)
    Mat<double> weights_ih = input->getWeights(hidden);
    cout << "\nInitial input->hidden weights size: " << weights_ih.size().cy << "x" << weights_ih.size().cx << endl;

    // Set simple weights: all 1.0 for first hidden node, all 0.5 for second
    Mat<double> new_weights_ih(2, 3, 0);
    new_weights_ih.setAt(0, 0, 1.0);  // hidden node 0, input 0
    new_weights_ih.setAt(0, 1, 1.0);  // hidden node 0, input 1
    new_weights_ih.setAt(0, 2, 1.0);  // hidden node 0, bias
    new_weights_ih.setAt(1, 0, 0.5);  // hidden node 1, input 0
    new_weights_ih.setAt(1, 1, 0.5);  // hidden node 1, input 1
    new_weights_ih.setAt(1, 2, 0.5);  // hidden node 1, bias
    input->setWeights(hidden, new_weights_ih);

    printMatrix(new_weights_ih, "Input->Hidden weights");

    // Set weights for hidden->output
    // Hidden layer has 2 nodes + 1 bias = 3 outputs
    // Output layer has 1 node
    // So weights should be (1, 3)
    Mat<double> new_weights_ho(1, 3, 0);
    new_weights_ho.setAt(0, 0, 1.0);  // output node 0, hidden 0
    new_weights_ho.setAt(0, 1, 1.0);  // output node 0, hidden 1
    new_weights_ho.setAt(0, 2, 1.0);  // output node 0, bias
    hidden->setWeights(output, new_weights_ho);

    printMatrix(new_weights_ho, "Hidden->Output weights");

    // Feed input [1.0, 1.0]
    Mat<double> inputData(1, 2, 0);
    inputData.setAt(0, 0, 1.0);
    inputData.setAt(0, 1, 1.0);

    printMatrix(inputData, "Input data");

    Mat<double> result = net->feed(inputData);

    printMatrix(result, "Network output");

    // Manual calculation:
    // Input layer receives [1.0, 1.0]
    // Input layer activates with sigmoid, then adds bias
    // Activated input: [sigmoid(1.0), sigmoid(1.0), 1.0] = [0.731, 0.731, 1.0]
    //
    // Hidden layer input = activated_input * weights_ih^T
    // For hidden node 0: 0.731*1.0 + 0.731*1.0 + 1.0*1.0 = 2.462
    // For hidden node 1: 0.731*0.5 + 0.731*0.5 + 1.0*0.5 = 1.231
    // Hidden activation: [sigmoid(2.462), sigmoid(1.231)]

    double sig1 = 1.0 / (1.0 + exp(-1.0));
    cout << "\nManual calculation:" << endl;
    cout << "sigmoid(1.0) = " << sig1 << endl;

    delete net;
    delete output;
    delete hidden;
    delete input;

    cout << "PASS: Forward propagation completed (manual verification needed)" << endl;
    return true;
}

// TEST 4: Verify layer input/output dimensions
bool test_layer_dimensions() {
    cout << "\n=== TEST 4: Layer Dimensions ===" << endl;

    Layer<double>* layer1 = new Layer<double>(3, "layer1");  // 3 input nodes
    Layer<double>* layer2 = new Layer<double>(2, "layer2");  // 2 input nodes

    cout << "Layer1 - Input nodes: " << layer1->getNumInputNodes()
         << ", Output nodes: " << layer1->getNumOutputNodes() << endl;
    cout << "Layer2 - Input nodes: " << layer2->getNumInputNodes()
         << ", Output nodes: " << layer2->getNumOutputNodes() << endl;

    layer1->connect(layer2);
    layer1->initWeights(layer2);

    Mat<double> weights = layer1->getWeights(layer2);
    cout << "Weights matrix size: " << weights.size().cy << " rows x "
         << weights.size().cx << " cols" << endl;

    // Hypothesis: Weights should be (layer2 input nodes, layer1 output nodes)
    // layer1 has 3 inputs + 1 bias = 4 outputs
    // layer2 has 2 inputs
    // So weights should be (2, 4)

    bool dimensionsCorrect = (weights.size().cy == layer2->getNumInputNodes() &&
                             weights.size().cx == layer1->getNumOutputNodes());

    if (!dimensionsCorrect) {
        cout << "FAIL: Weight dimensions incorrect!" << endl;
        cout << "Expected: (" << layer2->getNumInputNodes() << ", "
             << layer1->getNumOutputNodes() << ")" << endl;
        cout << "Actual: (" << weights.size().cy << ", " << weights.size().cx << ")" << endl;
        delete layer1;
        delete layer2;
        return false;
    }

    delete layer1;
    delete layer2;

    cout << "PASS: Layer dimensions correct" << endl;
    return true;
}

// TEST 5: Verify bias is added correctly
bool test_bias_addition() {
    cout << "\n=== TEST 5: Bias Addition ===" << endl;

    // Create a simple test matrix
    Mat<double> mat(1, 2, 0);
    mat.setAt(0, 0, 0.5);
    mat.setAt(0, 1, 0.7);

    printMatrix(mat, "Original matrix");

    // Add bias
    pushBiasCol<double>(mat);

    printMatrix(mat, "After adding bias");

    // Verify dimensions
    if (mat.size().cx != 3) {
        cout << "FAIL: Bias not added correctly. Expected 3 columns, got "
             << mat.size().cx << endl;
        return false;
    }

    // Verify bias value
    if (!approxEqual(mat.getAt(0, 2), 1.0)) {
        cout << "FAIL: Bias value incorrect. Expected 1.0, got "
             << mat.getAt(0, 2) << endl;
        return false;
    }

    cout << "PASS: Bias addition correct" << endl;
    return true;
}

// TEST 6: Test actual network computation with hand-verified values
bool test_network_computation_manual() {
    cout << "\n=== TEST 6: Network Computation (Hand-Verified) ===" << endl;

    // Create a minimal network: 1 input -> 1 output
    Layer<double>* input = new Layer<double>(1, "input");
    Layer<double>* output = new Layer<double>(1, "output");

    Network<double>* net = new Network<double>();
    net->setInputLayer(input);
    net->setOutputLayer(output);
    net->connect(input, output);
    net->init();

    // Set specific weights
    // Input has 1 node + 1 bias = 2 outputs
    // Output has 1 node
    // Weights: (1, 2)
    Mat<double> weights(1, 2, 0);
    weights.setAt(0, 0, 2.0);  // weight for input
    weights.setAt(0, 1, 0.5);  // weight for bias
    input->setWeights(output, weights);

    printMatrix(weights, "Weights");

    // Feed input [0.0]
    Mat<double> inputData(1, 1, 0);
    inputData.setAt(0, 0, 0.0);

    printMatrix(inputData, "Input");

    Mat<double> result = net->feed(inputData);

    printMatrix(result, "Output");

    // Manual calculation:
    // 1. Input layer receives [0.0]
    // 2. Input is activated: sigmoid(0.0) = 0.5
    // 3. Bias is added: [0.5, 1.0]
    // 4. Weighted sum for output: 0.5*2.0 + 1.0*0.5 = 1.0 + 0.5 = 1.5
    // 5. Output activation: sigmoid(1.5) ≈ 0.8176

    double sigmoid0 = 1.0 / (1.0 + exp(-0.0));  // = 0.5
    double weighted_sum = sigmoid0 * 2.0 + 1.0 * 0.5;  // = 1.5
    double expected_output = 1.0 / (1.0 + exp(-weighted_sum));  // sigmoid(1.5)

    cout << "\nManual calculation:" << endl;
    cout << "sigmoid(0.0) = " << sigmoid0 << endl;
    cout << "Weighted sum = " << weighted_sum << endl;
    cout << "Expected output = sigmoid(" << weighted_sum << ") = " << expected_output << endl;
    cout << "Actual output = " << result.getAt(0, 0) << endl;

    bool passed = approxEqual(result.getAt(0, 0), expected_output);

    delete net;
    delete output;
    delete input;

    if (!passed) {
        cout << "FAIL: Network computation does not match hand calculation!" << endl;
        return false;
    }

    cout << "PASS: Network computation matches hand calculation" << endl;
    return true;
}

// TEST 7: Test weight matrix multiplication orientation
bool test_weight_multiplication_orientation() {
    cout << "\n=== TEST 7: Weight Multiplication Orientation ===" << endl;

    // This test verifies that weights are multiplied in the correct order
    // activated_input (1, n+bias) * weights^T (n+bias, m) = output (1, m)

    Mat<double> activated_input(1, 3, 0);  // 2 nodes + 1 bias
    activated_input.setAt(0, 0, 1.0);
    activated_input.setAt(0, 1, 2.0);
    activated_input.setAt(0, 2, 1.0);  // bias

    // Weights for 2 output nodes: (2, 3)
    Mat<double> weights(2, 3, 0);
    weights.setAt(0, 0, 1.0);  // node0: w0
    weights.setAt(0, 1, 0.0);  // node0: w1
    weights.setAt(0, 2, 0.0);  // node0: bias
    weights.setAt(1, 0, 0.0);  // node1: w0
    weights.setAt(1, 1, 1.0);  // node1: w1
    weights.setAt(1, 2, 0.0);  // node1: bias

    printMatrix(activated_input, "Activated input (1x3)");
    printMatrix(weights, "Weights (2x3)");

    // Multiply: activated_input * weights^T
    // Use true because we want to treat each row of weights as a target node
    // (this matches how the network code uses it)
    Mat<double> result = Mult<double>(activated_input, weights, true);

    printMatrix(result, "Result (should be 1x2)");

    // Expected result:
    // node0: 1.0*1.0 + 2.0*0.0 + 1.0*0.0 = 1.0
    // node1: 1.0*0.0 + 2.0*1.0 + 1.0*0.0 = 2.0
    // Result: [1.0, 2.0]

    if (result.size().cy != 1 || result.size().cx != 2) {
        cout << "FAIL: Result dimensions incorrect!" << endl;
        return false;
    }

    if (!approxEqual(result.getAt(0, 0), 1.0) || !approxEqual(result.getAt(0, 1), 2.0)) {
        cout << "FAIL: Weight multiplication values incorrect!" << endl;
        cout << "Expected: [1.0, 2.0], Got: [" << result.getAt(0, 0) << ", "
             << result.getAt(0, 1) << "]" << endl;
        return false;
    }

    cout << "PASS: Weight multiplication orientation correct" << endl;
    return true;
}

int main() {
    cout << "========================================" << endl;
    cout << "CORE NETWORK CORRECTNESS TESTS" << endl;
    cout << "Testing with scientific rigor and evidence" << endl;
    cout << "========================================" << endl;

    int passed = 0;
    int total = 0;

    total++; if (test_matrix_multiplication()) passed++;
    total++; if (test_matrix_multiplication_2x3()) passed++;
    total++; if (test_bias_addition()) passed++;
    total++; if (test_layer_dimensions()) passed++;
    total++; if (test_weight_multiplication_orientation()) passed++;
    total++; if (test_network_computation_manual()) passed++;
    total++; if (test_forward_propagation_simple()) passed++;

    cout << "\n========================================" << endl;
    cout << "RESULTS: " << passed << "/" << total << " tests passed" << endl;
    cout << "========================================" << endl;

    return (passed == total) ? 0 : 1;
}