Solutions¶

#include <iostream>
#include <Eigen/Dense>
#include <print>

#include "utils_print.h"

int main()
{
    using Eigen::RowVector;
    using Eigen::VectorXi;
    
    RowVector<double, 5> v1{1.0, 2.0, 3.0, 4.0, 5.0};
    VectorXi v2{{1, 2, 3, 4, 5}};

    utils::print("v1", v1);
    utils::print("v2", v2);

    std::printf("v1[0] = %.1f\n", v1[0]);
    std::printf("v2[0] = %d\n", v2[0]);
}

#include <iostream>
#include <print>
#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::Matrix3d;
    using Eigen::Matrix;
    using std::cout;

    Matrix3d A;
    A.setZero();

    // A = Matrix3d::Zero();

    utils::print("A", A);

    Matrix<double, 2, 4> B;

    B << 1.0, 2.0, 3.0, 4.0,
         5.0, 6.0, 7.0, 8.0;

    utils::print("B", B);
}

#include <iostream>
#include <cmath>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::RowVector;

    // -5 -4 -3 -2 -1 0 1 2 3 4 5
    RowVector<double, 11> v1;

    v1.setLinSpaced(11, -5.0, 5.0);
    // v1 = Vector<double, 11>::LinSpaced(11, -5.0, 5.0);

    utils::print("v1", v1);

    auto v2 = v1.cwiseAbs().cwiseSqrt();
    auto v2a = v1.array().abs().sqrt().matrix();

    utils::print("v2", v2);
    utils::print("v2a", v2a);

    auto v3 = v1;

    for (auto i = 0; i < v3.size(); ++i)
        v3[i] = exp(v3[i]);

    utils::print("v3", v3);

    std::println("v1.sum() = {}", v1.sum());
    std::println("v1.prod() = {}", v1.prod());
    std::println("v1.mean() = {}", v1.mean());
    std::println("v1.maxCoeff() = {}", v1.maxCoeff());
    std::println("v1.minCoeff() = {}", v1.minCoeff());
}

#include <iostream>
#include <cmath>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::RowVector;
    using std::cout;

    RowVector<double, 10> v1;
    v1.setRandom();

    utils::print("v1", v1);

    RowVector<double, 10> v2;
    v2.setRandom();

    utils::print("v2", v2);

    utils::print("v1 + v2", v1 + v2);
    utils::print("v1 - v2", v1 - v2);
    utils::print("v1.cwiseProduct(v2)", v1.cwiseProduct(v2));
    std::println("v1.dot(v2) = {}", v1.dot(v2));
}

#include <iostream>
#include <cmath>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::Matrix3d;
    using std::cout;

    Matrix3d A;
    A.setRandom();

    utils::print("A", A);

    Matrix3d B;
    B.setRandom();

    utils::print("B", B);

    utils::print("A + B", A + B);
    utils::print("A - B", A - B);
    utils::print("A.array() * B.array()", A.array() * B.array());
    utils::print("A * B", A * B);
}

#include <iostream>
#include <cmath>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::Matrix3d;
    using Eigen::Vector3d;
    using std::cout;

    Matrix3d A;
    A.setRandom();

    utils::print("A", A);

    Vector3d b; // ColumnVector [3 x 1]
    b.setRandom();

    utils::print("b", b);
    utils::print("A * b", A * b);
}

#include <iostream>
#include <cmath>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::Matrix4d;
    using Eigen::Matrix3d;
    using Eigen::MatrixXd;
    using Eigen::FullPivLU;
    using std::cout;

    Matrix4d A;
    A.setRandom();

    utils::print("A", A);

    std::println("det(A) = {}", A.determinant());

    utils::print("A^(-1)", A.inverse());
    utils::print("A^T", A.transpose());
    utils::print("A * A^(-1)", A * A.inverse());

    if ((A * A.inverse()).isIdentity())
        std::println("A * A^(-1) is identity");
    else
        std::println("A * A^(-1) is not identity");

    Matrix3d B;

    B << 1, 2, 3, 4, 5, 6, 7, 8, 9;

    auto det = B.determinant();
    bool isSingular = std::abs(det) < 1e-10;

    if (isSingular)
        std::println("B is singular");
    else
        std::println("B is not singular");

    FullPivLU<MatrixXd> lu(B);

    bool isSingular_lu = !lu.isInvertible();

    if (isSingular_lu)
        std::println("B is singular (LU)");
    else
        std::println("B is not singular (LU)");
}

#include <iostream>
#include <cmath>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::Matrix4d;
    using Eigen::MatrixXd;
    using Eigen::Matrix2d;

    Matrix4d A;

    A << 1, 2, 3, 4, 
         5, 6, 7, 8, 
         9, 10, 11, 12, 
         13, 14, 15, 16;

    utils::print("A = ", A);

    MatrixXd B = A.block(0, 0, 2, 2);

    utils::print("B = ", B);

    MatrixXd C = A.block(2, 2, 2, 2);

    utils::print("C = ", C);

    Matrix2d D;

    D << 42, 42, 
         42, 42;

    utils::print("D = ", D);

    A.block(2, 2, 2, 2) = D;

    utils::print("A = ", A);
}

#include <iostream>
#include <cmath>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::Matrix;
    using Eigen::VectorXd;
    using std::cout;

    Matrix<double, 3, 4> A;

    A.setRandom();

    utils::print("A", A);

    utils::print("A.row(1)", A.row(1));
    utils::print("A.col(2)", A.col(2));

    A.row(1) << 1, 2, 3, 4;

    utils::print("A", A);

    VectorXd temp = A.col(2);
    A.col(2) = A.col(3);
    A.col(3) = temp;

    utils::print("A", A);

    A.col(2).swap(A.col(3));

    utils::print("A", A);
}

#include <iostream>
#include <cmath>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::Matrix;
    using Eigen::Vector;

    Matrix<double, 2, 3> A;

    A << 1, 2, 
         3, 4, 
         5, 6;

    utils::print("A", A);

    Matrix<double, 3, 2> B = A.reshaped(3, 2);

    utils::print("B", B);

    Vector<double, 6> v;

    v << 1, 2, 3, 4, 5, 6;

    utils::print("v", v);

    Matrix<double, 2, 3> C = v.reshaped(2, 3);

    utils::print("C", C);
}

#include <iostream>
#include <cmath>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::Matrix2d;
    using Eigen::MatrixXd;

    Matrix2d A;

    A << 1, 2, 
         3, 4;

    utils::print("A", A);

    Matrix2d B;

    B << 5, 6, 
         7, 8;

    utils::print("B", B);

    MatrixXd C(A.rows(), A.cols() + B.cols());

    C << A, B;

    utils::print("C", C);

    MatrixXd D(A.rows() + B.rows(), A.cols());

    D << A, B;

    utils::print("D", D);
}

#include <iostream>
#include <cmath>
#include <vector>
#include <print>

#include <Eigen/Dense>
#include "utils_print.h"

int main()
{
    using Eigen::Matrix;
    using Eigen::Matrix2d;
    using Eigen::MatrixXd;
    using Eigen::VectorXd;
    using Eigen::Index;

    Matrix<double, 5, 5> A;

    A <<  1,  2,  3,  4,  5, 
          6,  7,  8,  9, 10, 
         11, 12, 13, 14, 15, 
         16, 17, 18, 19, 20, 
         21, 22, 23, 24, 25;

    utils::print("A", A);

    std::vector<Index> rows = {1, 3, 4};
    std::vector<Index> cols = {0, 2};

    MatrixXd B = A(rows, cols);

    utils::print("B", B);

    utils::print("diagonal of A", A.diagonal());

    auto n = A.rows();

    VectorXd antiDiagonal(n);

    for (int i = 0; i < n; i++)
        antiDiagonal(i) = A(i, n - i - 1);

    utils::print("anti-diagonal of A", antiDiagonal);

    auto tiledRows = 4;
    auto tiledCols = 3;

    Matrix2d C;

    C << 1, 2, 3, 4;

    MatrixXd D(tiledRows * C.rows(), tiledCols * C.cols());

    for (auto i = 0; i < tiledRows; i++)
        for (auto j = 0; j < tiledCols; j++)
            D.block(i * C.rows(), j * C.cols(), C.rows(), C.cols()) = C;

    utils::print("D", D);

    // Demonstrate different formatting options
    std::print("--- Different Format Examples ---\n\n");
    
    Matrix2d E;
    E << 1.23456789, 2.3456789,
         3.456789,   4.56789;
    
    utils::print("E (Default)", E);
    utils::print("E (Compact)", E, PrintFormat::Compact);
    utils::print("E (Python)", E, PrintFormat::Python);
    utils::print("E (Clean)", E, PrintFormat::Clean);
    utils::print("E (Octave)", E, PrintFormat::Octave);
    utils::print("E (CSV)", E, PrintFormat::CSV);

    // Custom format with specific precision
    Eigen::IOFormat CustomFmt(2, 0, ", ", "\n", "[", "]");
    utils::print("E (Custom 2 decimals)", E, CustomFmt);
}

#include <iostream>
#include <cmath>
#include <vector>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::MatrixXd;
    using Eigen::VectorXd;

    MatrixXd A(10, 10);

    A << MatrixXd::Random(10, 10);

    utils::print("A", A);

    VectorXd b(10);

    b << VectorXd::Random(10);

    utils::print("b", b);

    VectorXd x = A.fullPivLu().solve(b);

    utils::print("The solution is:", x);

    std::println("The relative error is: {}", (A * x - b).norm() / b.norm());

    std::println("A * x = {}", A * x);
}

#include <iostream>
#include <cmath>
#include <vector>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

int main()
{
    using Eigen::Matrix3d;
    using Eigen::SelfAdjointEigenSolver;

    Matrix3d A;

    A << 1, 2, 7, 
         2, 5, 6, 
         7, 6, 10;

    utils::print("A", A);

    SelfAdjointEigenSolver<Eigen::Matrix3d> es(A);

    if (es.info() != Eigen::Success)
    {
        std::println("Eigen solver failed.");
        return 1;
    }

    utils::print("The eigenvalues of A are:", es.eigenvalues());
    utils::print("The matrix of eigenvectors, V, is:", es.eigenvectors());

    auto v1 = es.eigenvectors().col(0);
    auto v2 = es.eigenvectors().col(1);
    auto v3 = es.eigenvectors().col(2);

    auto lambda1 = es.eigenvalues()(0);
    auto lambda2 = es.eigenvalues()(1);
    auto lambda3 = es.eigenvalues()(2);

    utils::print("The first eigenvector is:", v1);
    utils::print("The second eigenvector is:", v2);
    utils::print("The third eigenvector is:", v3);

    std::println();

    std::println("The first eigenvalue is: {}", lambda1);
    std::println("The second eigenvalue is: {}", lambda2);
    std::println("The third eigenvalue is: {}\n", lambda3);

    Matrix3d I = Matrix3d::Identity();

    std::println("(A-lambda1 * I) * v1 = {}", (A - lambda1 * I) * v1);
    std::println("(A-lambda2 * I) * v2 = {}", (A - lambda2 * I) * v2);
    std::println("(A-lambda3 * I) * v3 = {}", (A - lambda3 * I) * v3);
}

#include <cmath>
#include <vector>
#include <numbers>
#include <print>

#include <Eigen/Dense>

#include "utils_print.h"

#include <img/Image.h>

Eigen::MatrixXd adjustBrightnessContrast(const Eigen::MatrixXd &image, double alpha, double beta)
{
    // Convert to array for element-wise operations
    Eigen::ArrayXXd result = image.array();

    // Apply contrast (multiply)
    result = alpha * result;

    // Apply brightness (add)
    result = result + beta;

    // Clamp values to valid pixel range [0, 1]
    result = result.max(0.0).min(1.0);

    return result.matrix();
}

Eigen::MatrixXd createGaussianKernel(int size, double sigma)
{
    if (size % 2 == 0)
    {
        size++; // Make sure kernel size is odd
    }

    Eigen::MatrixXd kernel(size, size);
    int center = size / 2;
    double sum = 0.0;

    // Fill the kernel with Gaussian values
    for (int y = 0; y < size; y++)
    {
        for (int x = 0; x < size; x++)
        {
            int dx = x - center;
            int dy = y - center;
            double exponent = -(dx * dx + dy * dy) / (2 * sigma * sigma);
            kernel(y, x) = std::exp(exponent) / (2 * std::numbers::pi * sigma * sigma);
            sum += kernel(y, x);
        }
    }

    // Normalize the kernel so its sum equals 1
    kernel /= sum;

    return kernel;
}

Eigen::MatrixXd applyGaussianFilter(const Eigen::MatrixXd &image, int kernelSize, double sigma)
{
    Eigen::MatrixXd result = image;
    Eigen::MatrixXd kernel = createGaussianKernel(kernelSize, sigma);

    int kernelRadius = kernelSize / 2;

    for (int y = kernelRadius; y < image.rows() - kernelRadius; y++)
    {
        for (int x = kernelRadius; x < image.cols() - kernelRadius; x++)
        {
            double sum = 0.0;
            for (int ky = 0; ky < kernelSize; ky++)
            {
                for (int kx = 0; kx < kernelSize; kx++)
                {
                    int imgX = x + kx - kernelRadius;
                    int imgY = y + ky - kernelRadius;
                    sum += image(imgY, imgX) * kernel(ky, kx);
                }
            }
            result(y, x) = sum;
        }
    }
    return result;
}

img::ImageGd copyToImage(const Eigen::MatrixXd &mat)
{
    img::ImageGd image;
    const int rows = static_cast<int>(mat.rows());
    const int cols = static_cast<int>(mat.cols());
    
    image.resize(rows, cols);
    for (int i = 0; i < rows; i++)
        for (int j = 0; j < cols; j++)
            image(i, j) = mat(i, j);

    return image;
}

int main()
{
    using std::cout;
    using Eigen::MatrixXd;
    using img::ImageGd;

    ImageGd image;

    std::println("Loading image...");

    if (!load("images/half-moon-986269.png", image))
    {
        std::println("Error loading image");
        return 1;
    }

    // Create a matrix with the image data

    auto mat = image.as_matrix();

    std::println("Adjusting brightness contrast...");

    MatrixXd brighterImageMat = adjustBrightnessContrast(mat, 1.0, 0.5);

    // Create a new image with the modified data

    std::println("Converting to image...");

    auto brighterImage = copyToImage(brighterImageMat);

    // Save the new image

    std::println("Saving image...");

    if (!save("images/brighter-half-moon-986269.png", brighterImage))
    {
        std::println("Error saving image");
        return 1;
    }

    std::println("Applying Gaussian filter...");

    MatrixXd blurredImageMat = applyGaussianFilter(mat, 10, 1.0);

    // Create a new image with the modified data

    std::println("Converting to image...");

    auto blurredImage = copyToImage(blurredImageMat);

    // Save the new image

    std::println("Saving image...");

    if (!save("images/blurred-half-moon-986269.png", blurredImage))
    {
        std::println("Error saving image");
        return 1;
    }
}