data_analysis.zig

//! Zigen Example: PCA Data Analysis
//! Principal Component Analysis on iris-like flower measurements.

const std = @import("std");
const Zigen = @import("zigen");

pub fn main() !void {
    std.debug.print(
        \\
        \\  ╔═══════════════════════════════════════════════╗
        \\  ║    Zigen — PCA Data Analysis                  ║
        \\  ╚═══════════════════════════════════════════════╝
        \\
        \\
    , .{});

    // ── Dataset ───────────────────────────────────────────────────────
    std.debug.print("━━ Dataset: Flower Measurements (10 samples) ━━━━━━━━━\n\n", .{});

    const data = [10][4]f32{
        .{ 5.1, 3.5, 1.4, 0.2 }, // Setosa
        .{ 4.9, 3.0, 1.4, 0.2 },
        .{ 5.4, 3.9, 1.7, 0.4 },
        .{ 4.6, 3.4, 1.4, 0.3 },
        .{ 5.0, 3.4, 1.5, 0.2 },
        .{ 7.0, 3.2, 4.7, 1.4 }, // Versicolor
        .{ 6.4, 3.2, 4.5, 1.5 },
        .{ 6.9, 3.1, 4.9, 1.5 },
        .{ 5.5, 2.3, 4.0, 1.3 },
        .{ 6.5, 2.8, 4.6, 1.5 },
    };

    for (data, 0..) |row, i| {
        std.debug.print("  {d:>2}: [{d:.1}, {d:.1}, {d:.1}, {d:.1}]\n", .{
            i + 1, row[0], row[1], row[2], row[3],
        });
    }
    std.debug.print("\n", .{});

    // ── Step 1: Compute mean ──────────────────────────────────────────
    std.debug.print("━━ Step 1: Feature Means ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    var mean = [4]f32{ 0, 0, 0, 0 };
    for (data) |row| {
        for (0..4) |j| mean[j] += row[j];
    }
    for (0..4) |j| mean[j] /= 10.0;

    std.debug.print("  Mean: [{d:.3}, {d:.3}, {d:.3}, {d:.3}]\n\n", .{
        mean[0], mean[1], mean[2], mean[3],
    });

    // ── Step 2: Covariance matrix ─────────────────────────────────────
    std.debug.print("━━ Step 2: Covariance Matrix ━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    var cov_data: [4][4]f32 = undefined;
    for (0..4) |i| {
        for (0..4) |j| {
            var s: f32 = 0;
            for (data) |row| s += (row[i] - mean[i]) * (row[j] - mean[j]);
            cov_data[i][j] = s / 9.0;
        }
    }

    const cov = Zigen.Matrix4f.fromArray(cov_data);
    for (0..4) |i| {
        std.debug.print("  | {d:>7.4} {d:>7.4} {d:>7.4} {d:>7.4} |\n", .{
            cov.at(i, 0), cov.at(i, 1), cov.at(i, 2), cov.at(i, 3),
        });
    }
    std.debug.print("\n", .{});

    // ── Step 3: Eigenvalues (PCA) ─────────────────────────────────────
    std.debug.print("━━ Step 3: Principal Components ━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    const eig = Zigen.SelfAdjointEigenSolver(f32, 4).compute(cov);
    const ev = eig.eigenvalues();

    const total_var = ev.at(0) + ev.at(1) + ev.at(2) + ev.at(3);

    std.debug.print("  Component | Eigenvalue | Variance\n", .{});
    var i: usize = 4;
    while (i > 0) {
        i -= 1;
        const pct = ev.at(i) / total_var * 100.0;
        std.debug.print("  PC{d}       |   {d:>7.4}  |  {d:>5.1}%\n", .{
            4 - i, ev.at(i), pct,
        });
    }

    const top2_var = ev.at(3) + ev.at(2);
    std.debug.print("\n  Top 2 PCs explain {d:.1}% of variance\n\n", .{
        top2_var / total_var * 100.0,
    });

    // ── Step 4: Principal direction ───────────────────────────────────
    std.debug.print("━━ Step 4: First Principal Direction ━━━━━━━━━━━━━━━━━━\n\n", .{});

    const evecs = eig.eigenvectors();
    const pc1 = evecs.col(3);
    std.debug.print("  PC1 = [{d:.4}, {d:.4}, {d:.4}, {d:.4}]\n\n", .{
        pc1.at(0), pc1.at(1), pc1.at(2), pc1.at(3),
    });

    std.debug.print("PCA analysis complete!\n", .{});
}