// example_vector.cpp -- Vector API examples
// All code snippets are subroutines called from main()
//
// Build:
//   cl /std:c++latest /EHsc /O2 /MD /utf-8 /I<calx>/include examples\example_vector.cpp /MACHINE:X64

#include <math/core/vector.hpp>
#include <iostream>
#include <iomanip>
#include <cmath>

using namespace calx;

// ============================================================
// Constructors
// ============================================================

void demo_constructors() {
    std::cout << "--- constructors ---" << std::endl;

    // Default (empty)
    Vector<double> empty;
    std::cout << "Vector<double>()       size = " << empty.size() << std::endl;

    // Size (zero-initialized)
    Vector<double> zeros(5);
    std::cout << "Vector<double>(5)      = [";
    for (size_t i = 0; i < zeros.size(); ++i) std::cout << (i ? ", " : "") << zeros[i];
    std::cout << "]" << std::endl;

    // Size + value
    Vector<double> ones(4, 1.0);
    std::cout << "Vector<double>(4, 1.0) = [";
    for (size_t i = 0; i < ones.size(); ++i) std::cout << (i ? ", " : "") << ones[i];
    std::cout << "]" << std::endl;

    // Initializer list
    Vector<double> v{1.0, 2.0, 3.0, 4.0, 5.0};
    std::cout << "Vector{1,2,3,4,5}      = [";
    for (size_t i = 0; i < v.size(); ++i) std::cout << (i ? ", " : "") << v[i];
    std::cout << "]" << std::endl;
}

void demo_static_vector() {
    std::cout << "--- StaticVector ---" << std::endl;

    StaticVector<double, 3> a{1.0, 2.0, 3.0};
    StaticVector<double, 3> b{4.0, 5.0, 6.0};
    std::cout << "a = [" << a[0] << ", " << a[1] << ", " << a[2] << "]" << std::endl;
    std::cout << "b = [" << b[0] << ", " << b[1] << ", " << b[2] << "]" << std::endl;
    std::cout << "dot(a, b)   = " << dot(a, b) << std::endl;
    auto c = cross(a, b);
    std::cout << "cross(a, b) = [" << c[0] << ", " << c[1] << ", " << c[2] << "]" << std::endl;
}

// ============================================================
// Element access
// ============================================================

void demo_element_access() {
    std::cout << "--- element access ---" << std::endl;

    Vector<double> v{10.0, 20.0, 30.0, 40.0, 50.0};
    std::cout << "v = [";
    for (size_t i = 0; i < v.size(); ++i) std::cout << (i ? ", " : "") << v[i];
    std::cout << "]" << std::endl;
    std::cout << "v[0]     = " << v[0] << std::endl;
    std::cout << "v[4]     = " << v[4] << std::endl;
    std::cout << "v.at(2)  = " << v.at(2) << std::endl;
    std::cout << "v.front()= " << v.front() << std::endl;
    std::cout << "v.back() = " << v.back() << std::endl;
    std::cout << "v.size() = " << v.size() << std::endl;
}

// ============================================================
// Arithmetic (expression templates)
// ============================================================

static void print_vec(const char* label, const Vector<double>& v) {
    std::cout << label << " = [";
    for (size_t i = 0; i < v.size(); ++i) std::cout << (i ? ", " : "") << v[i];
    std::cout << "]" << std::endl;
}

void demo_arithmetic() {
    std::cout << "--- arithmetic (expression templates) ---" << std::endl;

    Vector<double> a{1.0, 2.0, 3.0};
    Vector<double> b{4.0, 5.0, 6.0};
    print_vec("a", a);
    print_vec("b", b);

    // Expression templates: no temporaries until assignment
    Vector<double> c = a + b;
    print_vec("a + b", c);

    Vector<double> d = a - b;
    print_vec("a - b", d);

    Vector<double> e = 2.0 * a;
    print_vec("2 * a", e);

    Vector<double> f = a + 2.0 * b;
    print_vec("a + 2*b", f);
}

void demo_compound_assignment() {
    std::cout << "--- compound assignment ---" << std::endl;

    Vector<double> v{1.0, 2.0, 3.0};
    print_vec("v", v);
    v += Vector<double>{10.0, 20.0, 30.0};
    print_vec("v += {10,20,30}", v);
    v *= 0.5;
    print_vec("v *= 0.5", v);
}

// ============================================================
// Dot product, norms
// ============================================================

void demo_dot_norm() {
    std::cout << "--- dot product & norms ---" << std::endl;

    Vector<double> a{1.0, 2.0, 3.0};
    Vector<double> b{4.0, 5.0, 6.0};
    print_vec("a", a);
    print_vec("b", b);

    std::cout << "dot(a, b)    = " << dot(a, b) << std::endl;
    std::cout << "a.dot(b)     = " << a.dot(b) << std::endl;
    std::cout << "norm(a)      = " << norm(a) << "  (L2)" << std::endl;
    std::cout << "norm_l1(a)   = " << norm_l1(a) << std::endl;
    std::cout << "norm_linf(a) = " << norm_linf(a) << std::endl;
    std::cout << "norm_lp(a,3) = " << norm_lp(a, 3.0) << std::endl;
}

// ============================================================
// Normalization
// ============================================================

void demo_normalize() {
    std::cout << "--- normalization ---" << std::endl;

    Vector<double> v{3.0, 4.0};
    print_vec("v", v);
    std::cout << "norm(v) = " << norm(v) << std::endl;

    Vector<double> u = normalized(v);
    print_vec("normalized(v)", u);
    std::cout << "norm(normalized(v)) = " << norm(u) << std::endl;

    v.normalize();  // in-place
    print_vec("v after normalize()", v);
}

// ============================================================
// Block operations (head, tail, segment)
// ============================================================

void demo_block_ops() {
    std::cout << "--- block operations (head/tail/segment) ---" << std::endl;

    Vector<double> v{10.0, 20.0, 30.0, 40.0, 50.0};
    print_vec("v", v);

    // Read through view
    auto h = v.head(3);
    std::cout << "v.head(3) = [";
    for (size_t i = 0; i < h.size(); ++i) std::cout << (i ? ", " : "") << h[i];
    std::cout << "]" << std::endl;

    auto t = v.tail(2);
    std::cout << "v.tail(2) = [";
    for (size_t i = 0; i < t.size(); ++i) std::cout << (i ? ", " : "") << t[i];
    std::cout << "]" << std::endl;

    auto s = v.segment(1, 3);
    std::cout << "v.segment(1,3) = [";
    for (size_t i = 0; i < s.size(); ++i) std::cout << (i ? ", " : "") << s[i];
    std::cout << "]" << std::endl;

    // Write through view
    v.head(2) = Vector<double>{99.0, 88.0};
    print_vec("after v.head(2) = {99,88}", v);
}

// ============================================================
// BLAS Level-1: axpy, axpby
// ============================================================

void demo_blas() {
    std::cout << "--- BLAS Level-1: axpy / axpby ---" << std::endl;

    Vector<double> x{1.0, 2.0, 3.0};
    Vector<double> y{10.0, 20.0, 30.0};
    print_vec("x", x);
    print_vec("y", y);

    axpy(2.0, x, y);  // y += 2*x
    print_vec("after axpy(2, x, y): y", y);

    Vector<double> result(3);
    axpby(0.5, x, 0.1, y, result);  // result = 0.5*x + 0.1*y
    print_vec("axpby(0.5, x, 0.1, y)", result);
}

// ============================================================
// VectorMap (zero-copy view of external memory)
// ============================================================

void demo_vector_map() {
    std::cout << "--- VectorMap (zero-copy view) ---" << std::endl;

    double data[] = {1.0, 2.0, 3.0, 4.0};
    VectorMap<double> m(data, 4);
    std::cout << "data[] = {1, 2, 3, 4}" << std::endl;
    std::cout << "map[2] = " << m[2] << std::endl;

    m[2] = 99.0;
    std::cout << "after map[2] = 99: data[2] = " << data[2] << std::endl;

    // Convert to Vector (copy)
    Vector<double> v = m;
    print_vec("Vector from map", v);
}

// ============================================================
// Range-based for loop & iterators
// ============================================================

void demo_iterators() {
    std::cout << "--- iterators & range-for ---" << std::endl;

    Vector<double> v{1.0, 4.0, 9.0, 16.0};
    std::cout << "v = [";
    for (size_t i = 0; i < v.size(); ++i) std::cout << (i ? ", " : "") << v[i];
    std::cout << "]" << std::endl;
    std::cout << "sqrt of each element: [";
    bool first = true;
    for (double x : v) {
        if (!first) std::cout << ", ";
        std::cout << std::sqrt(x);
        first = false;
    }
    std::cout << "]" << std::endl;
}

// ============================================================
// Cross product (StaticVector<T,3>)
// ============================================================

void demo_cross_product() {
    std::cout << "--- cross product (3D) ---" << std::endl;

    StaticVector<double, 3> i{1.0, 0.0, 0.0};
    StaticVector<double, 3> j{0.0, 1.0, 0.0};
    std::cout << "i = [" << i[0] << ", " << i[1] << ", " << i[2] << "]" << std::endl;
    std::cout << "j = [" << j[0] << ", " << j[1] << ", " << j[2] << "]" << std::endl;
    auto k = cross(i, j);
    std::cout << "i x j = [" << k[0] << ", " << k[1] << ", " << k[2] << "]  (= k)" << std::endl;
    auto ji = cross(j, i);
    std::cout << "j x i = [" << ji[0] << ", " << ji[1] << ", " << ji[2] << "]  (= -k)" << std::endl;
}

// ============================================================
// main
// ============================================================

int main() {
    std::cout << "=== calx Vector API Examples ===" << std::endl << std::endl;

    demo_constructors();
    demo_static_vector();
    std::cout << std::endl;

    demo_element_access();
    std::cout << std::endl;

    demo_arithmetic();
    demo_compound_assignment();
    std::cout << std::endl;

    demo_dot_norm();
    std::cout << std::endl;

    demo_normalize();
    std::cout << std::endl;

    demo_block_ops();
    std::cout << std::endl;

    demo_blas();
    std::cout << std::endl;

    demo_vector_map();
    std::cout << std::endl;

    demo_iterators();
    demo_cross_product();

    std::cout << std::endl;
    std::cout << "=== All examples completed ===" << std::endl;
    return 0;
}