// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// example_tensor.cpp
// Tensor<T> as a multidimensional array: construction, shape/reshape/transpose,
// element-wise (Hadamard) product, contraction (= matmul), outer product, and
// Matrix interop.

#include <math/core/Tensor.hpp>
#include <math/core/matrix.hpp>
#include <iomanip>
#include <iostream>

using namespace sangi;

template <typename T>
static void print_tensor(const char* label, const Tensor<T>& t) {
    std::cout << "  " << label << " shape=[";
    for (std::size_t i = 0; i < t.rank(); ++i) std::cout << (i ? "," : "") << t.shape(i);
    std::cout << "] data=[";
    for (std::size_t i = 0; i < t.size(); ++i) std::cout << (i ? "," : "") << t.flat(i);
    std::cout << "]\n";
}

int main() {
    std::cout << std::setprecision(15);

    std::cout << "=== construction ===\n";
    Tensor<double> t1({2, 3});                                  // rank-2, zeros
    Tensor<double> t2({2, 3}, 1.0);                             // all ones
    Tensor<double> t3({2, 3}, {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}); // row-major fill
    print_tensor("zeros", t1);
    print_tensor("ones ", t2);
    print_tensor("fill ", t3);

    std::cout << "\n=== shape / reshape / transpose ===\n";
    auto r = t3.reshape({3, 2});      // same data, new shape
    auto tT = t3.transpose();         // reverse axes (2x3 -> 3x2 view)
    print_tensor("reshape", r);
    std::cout << "  transpose rank=" << tT.rank()
              << " shape={" << tT.shape(0) << "," << tT.shape(1) << "}\n";

    std::cout << "\n=== Hadamard (element-wise) product ===\n";
    Tensor<double> a({2, 2}, {1.0, 2.0, 3.0, 4.0});
    Tensor<double> b({2, 2}, {5.0, 6.0, 7.0, 8.0});
    print_tensor("a*b", hadamard(a, b));

    std::cout << "\n=== contraction = matrix product ===\n";
    // (3x4) . (4x5) contracting axis 1 of M with axis 0 of N -> (3x5)
    Tensor<double> M({3, 4}, 1.0);
    Tensor<double> N({4, 5}, 2.0);
    Tensor<double> C = contract(M, N, {{1, 0}});
    std::cout << "  C shape={" << C.shape(0) << "," << C.shape(1)
              << "} C.flat(0)=" << C.flat(0) << "  (= 4*1*2 = 8)\n";

    std::cout << "\n=== outer product ===\n";
    Tensor<double> u({3}, {1.0, 2.0, 3.0});
    Tensor<double> v({2}, {10.0, 20.0});
    print_tensor("u(x)v", outerProduct(u, v));

    std::cout << "\n=== Matrix interop ===\n";
    Matrix<double> mat(2, 2, 1.0);
    Tensor<double> tm(mat);                 // Matrix -> Tensor
    Matrix<double> back = tm.toMatrix();    // Tensor -> Matrix
    std::cout << "  rank=" << tm.rank() << " size=" << tm.size()
              << " back(1,1)=" << back(1, 1) << "\n";

    return 0;
}
