PHP Classes

File: matrix.cpp

Recommend this page to a friend!
  Classes of Cuthbert Martin Lwinga   PHP Matrix Extension Class   matrix.cpp   Download  
File: matrix.cpp
Role: Auxiliary data
Content type: text/plain
Description: Auxiliary data
Class: PHP Matrix Extension Class
Perform matrix calculations using a PHP extension
Author: By
Last change:
Date: 1 month ago
Size: 12,322 bytes
 

Contents

Class file image Download
#include "matrix.h" #include <stdexcept> #include <iostream> #include <thread> #include <vector> #include <cmath> #include <functional> Matrix::Matrix(int rows, int cols, double value) : data(rows, cols) { data.setConstant(value); } Matrix::Matrix(const Eigen::MatrixXd& other) : data(other) {} Matrix::Matrix(const std::vector<std::vector<double>>& inputData) { int rows = inputData.size(); int cols = rows > 0 ? inputData[0].size() : 0; data.resize(rows, cols); // Check if all rows have the same number of columns for (int i = 0; i < rows; ++i) { if (inputData[i].size() != cols) { throw std::invalid_argument("All rows must have the same number of columns"); } } // Log the number of rows and columns std::cout << "Initializing matrix with " << rows << " rows and " << cols << " columns using multi-threading." << std::endl; // Determine the number of threads int numThreads = std::min(20, static_cast<int>(std::thread::hardware_concurrency())); int rowsPerThread = rows / numThreads; // Define a lambda function to initialize a range of elements auto initRange = [&](int startRow, int endRow) { for (int i = startRow; i < endRow; ++i) { for (int j = 0; j < cols; ++j) { data(i, j) = inputData[i][j]; } } // Debug: Log each thread's range std::cout << "Thread completed rows " << startRow << " to " << endRow << std::endl; }; // Create and run threads std::vector<std::thread> threads; for (int i = 0; i < numThreads; ++i) { int startRow = i * rowsPerThread; int endRow = (i == numThreads - 1) ? rows : (i + 1) * rowsPerThread; threads.emplace_back(initRange, startRow, endRow); } // Wait for all threads to finish for (auto& thread : threads) { if (thread.joinable()) { thread.join(); } } std::cout << "Matrix initialization complete." << std::endl; } void addRange(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B, Eigen::MatrixXd& C, int startRow, int endRow) { C.block(startRow, 0, endRow - startRow, A.cols()) = A.block(startRow, 0, endRow - startRow, A.cols()) + B.block(startRow, 0, endRow - startRow, B.cols()); } void subtractRange(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B, Eigen::MatrixXd& C, int startRow, int endRow) { C.block(startRow, 0, endRow - startRow, A.cols()) = A.block(startRow, 0, endRow - startRow, A.cols()) - B.block(startRow, 0, endRow - startRow, B.cols()); } Matrix Matrix::add(const Matrix& other) const { if (data.rows() != other.data.rows() || data.cols() != other.data.cols()) { throw std::invalid_argument("Matrix dimensions must match for addition"); } // Enable Eigen's multi-threading Eigen::setNbThreads(std::thread::hardware_concurrency()); Matrix result(data.rows(), data.cols()); result.data = data + other.data; return result; } Matrix Matrix::subtract(const Matrix& other) const { if (data.rows() != other.data.rows() || data.cols() != other.data.cols()) { throw std::invalid_argument("Matrix dimensions must match for subtraction"); } // Enable Eigen's multi-threading Eigen::setNbThreads(std::thread::hardware_concurrency()); Matrix result(data.rows(), data.cols()); result.data = data - other.data; return result; } Matrix Matrix::log() const { Matrix result(data.rows(), data.cols()); result.data = data.array().log(); return result; } Matrix Matrix::sqrt() const { Matrix result(data.rows(), data.cols()); result.data = data.array().sqrt(); return result; } Matrix Matrix::exp(double base) const { Matrix result(data.rows(), data.cols()); result.data = data.array().pow(base); return result; } Matrix Matrix::sum(int axis) const { if (axis == -1) { // Sum all elements and return a 1x1 matrix double total_sum = data.sum(); return Matrix(1, 1, total_sum); } else if (axis == 0) { // Sum along columns (resulting in a row vector) Eigen::RowVectorXd col_sum = data.colwise().sum(); Matrix result(1, col_sum.size()); result.data = col_sum; return result; } else if (axis == 1) { // Sum along rows (resulting in a column vector) Eigen::VectorXd row_sum = data.rowwise().sum(); Matrix result(row_sum.size(), 1); result.data = row_sum; return result; } else { throw std::invalid_argument("Invalid axis for sum. Axis must be -1, 0, or 1."); } } Matrix Matrix::inverse() const { if (data.rows() != data.cols()) { throw std::invalid_argument("Matrix must be square to calculate its inverse."); } return Matrix(data.inverse().eval()); } double Matrix::determinant() const { if (data.rows() != data.cols()) { throw std::invalid_argument("Matrix must be square to calculate its determinant."); } return data.determinant(); } std::pair<Matrix, Matrix> Matrix::eigen() const { if (data.rows() != data.cols()) { throw std::invalid_argument("Matrix must be square to calculate its eigenvalues and eigenvectors."); } Eigen::EigenSolver<Eigen::MatrixXd> solver(data); Matrix eigenvalues(1, data.rows()); Matrix eigenvectors(data.rows(), data.cols()); eigenvalues.data = solver.eigenvalues().real(); eigenvectors.data = solver.eigenvectors().real(); return std::make_pair(eigenvalues, eigenvectors); } Matrix Matrix::operator*(double scalar) const { Matrix result(data.rows(), data.cols()); result.data = data * scalar; return result; } Matrix Matrix::transpose() const { Matrix result(data.cols(), data.rows()); result.data = data.transpose(); return result; } void dotProductRange(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B, Eigen::MatrixXd& C, int startRow, int endRow) { C.block(startRow, 0, endRow - startRow, B.cols()) = A.block(startRow, 0, endRow - startRow, A.cols()) * B; } Matrix Matrix::dot(const Matrix& other) const { if (data.cols() != other.data.rows()) { throw std::invalid_argument("Matrix dimensions must be compatible for dot product"); } Matrix result(data.rows(), other.data.cols()); int numThreads = std::min(20, static_cast<int>(std::thread::hardware_concurrency())); // Use at most 20 threads or the number of available hardware threads int rowsPerThread = data.rows() / numThreads; auto dataRef = std::cref(data); auto otherDataRef = std::cref(other.data); auto resultDataRef = std::ref(result.data); std::vector<std::thread> threads; for (int i = 0; i < numThreads; ++i) { int startRow = i * rowsPerThread; int endRow = (i == numThreads - 1) ? data.rows() : (i + 1) * rowsPerThread; threads.push_back(std::thread(dotProductRange, dataRef, otherDataRef, resultDataRef, startRow, endRow)); } for (auto& t : threads) { t.join(); } return result; } Matrix Matrix::addScalar(const double other) const { Matrix result(data.rows(), data.cols()); result.data = data.array() + other; return result; } Matrix Matrix::subtractScalar(const double other) const { Matrix result(data.rows(), data.cols()); result.data = data.array() - other; return result; } // other functions std::vector<int> Matrix::argmax(int axis) const { std::vector<int> indices; if (axis == 0) { // Max along columns indices.resize(data.cols()); auto find_max_in_col = [&](int start, int end) { for (int j = start; j < end; ++j) { int maxIndex = 0; for (int i = 1; i < data.rows(); ++i) { if (data(i, j) > data(maxIndex, j)) { maxIndex = i; } } indices[j] = maxIndex; } }; int num_threads = std::thread::hardware_concurrency(); std::vector<std::thread> threads(num_threads); int cols_per_thread = data.cols() / num_threads; int remaining_cols = data.cols() % num_threads; int start = 0; for (int i = 0; i < num_threads; ++i) { int end = start + cols_per_thread + (remaining_cols > 0 ? 1 : 0); if (remaining_cols > 0) remaining_cols--; threads[i] = std::thread(find_max_in_col, start, end); start = end; } for (auto& th : threads) { th.join(); } } else if (axis == 1) { // Max along rows indices.resize(data.rows()); auto find_max_in_row = [&](int start, int end) { for (int i = start; i < end; ++i) { int maxIndex = 0; for (int j = 1; j < data.cols(); ++j) { if (data(i, j) > data(i, maxIndex)) { maxIndex = j; } } indices[i] = maxIndex; } }; int num_threads = std::thread::hardware_concurrency(); std::vector<std::thread> threads(num_threads); int rows_per_thread = data.rows() / num_threads; int remaining_rows = data.rows() % num_threads; int start = 0; for (int i = 0; i < num_threads; ++i) { int end = start + rows_per_thread + (remaining_rows > 0 ? 1 : 0); if (remaining_rows > 0) remaining_rows--; threads[i] = std::thread(find_max_in_row, start, end); start = end; } for (auto& th : threads) { th.join(); } } else { throw std::invalid_argument("Axis must be 0 or 1"); } return indices; } Matrix Matrix::clip(double min_val, double max_val) const { Matrix result(data.rows(), data.cols()); result.data = data.array().min(max_val).max(min_val); return result; } // Operator overloads for +, -, * Matrix Matrix::operator/(const Matrix& other) const { if (data.rows() != other.getRows() || data.cols() != other.getCols()) { throw std::invalid_argument("Matrices must have the same dimensions for division."); } Matrix result(data.rows(), data.cols()); auto divide_row = [&](int start, int end) { for (int i = start; i < end; ++i) { for (int j = 0; j < data.cols(); ++j) { if (other(i, j) == 0) { throw std::invalid_argument("Division by zero encountered in matrix division."); } result(i, j) = data(i, j) / other(i, j); } } }; int num_threads = std::thread::hardware_concurrency(); std::vector<std::thread> threads(num_threads); int rows_per_thread = data.rows() / num_threads; int remaining_rows = data.rows() % num_threads; int start = 0; for (int i = 0; i < num_threads; ++i) { int end = start + rows_per_thread + (remaining_rows > 0 ? 1 : 0); if (remaining_rows > 0) remaining_rows--; threads[i] = std::thread(divide_row, start, end); start = end; } for (auto& th : threads) { th.join(); } return result; } Matrix Matrix::operator/(double scalar) const { if (scalar == 0) { throw std::invalid_argument("Division by zero encountered in scalar division."); } Matrix result(data.rows(), data.cols()); auto divide_row = [&](int start, int end) { for (int i = start; i < end; ++i) { for (int j = 0; j < data.cols(); ++j) { result(i, j) = data(i, j) / scalar; } } }; int num_threads = std::thread::hardware_concurrency(); std::vector<std::thread> threads(num_threads); int rows_per_thread = data.rows() / num_threads; int remaining_rows = data.rows() % num_threads; int start = 0; for (int i = 0; i < num_threads; ++i) { int end = start + rows_per_thread + (remaining_rows > 0 ? 1 : 0); if (remaining_rows > 0) remaining_rows--; threads[i] = std::thread(divide_row, start, end); start = end; } for (auto& th : threads) { th.join(); } return result; } void Matrix::display() const { std::cout << data << std::endl; }