Skip to content
/ matrix Public

A modern, C++20-native, single-file header-only dense 2D matrix library.

82 stars 8 forks Branches Tags Activity
Star
Notifications You must be signed in to change notification settings

fengwang/matrix

BranchesTags

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

374 Commits
.github/workflows
.github/workflows
 
 
examples
examples
 
 
images
images
 
 
tests
tests
 
 
.gitignore
.gitignore
 
 
Makefile
Makefile
 
 
ReadMe.md
ReadMe.md
 
 
TODO
TODO
 
 
matrix.hpp
matrix.hpp
 
 

Repository files navigation

A modern, C++20-native, single-file header-only dense 2D matrix library.

CI

Contents

Example usage

including the header file

//your_source_code.cpp
#include "matrix.hpp"

typical compile and link command

g++ -o your_exe_file your_source_code.cpp -std=c++20 -O2 -pthread -lstdc++fs

Please note std::thread is not enabled by default. If you prefer multi-thread mode, pass -DPARALLEL option to compiler, and add necessary link options. std::filesystem is used, make sure corresponding library option is passed during link time (-lstdc++fs for g++).

Variadic macro __VA_OPT__ is used. It is officially supported since c++20(link1, link2), so the compiler must be compatible with c++20.

basic

creating matrices

  • generating matrix

    • creating a matrix of size 12 X 34:
    feng::matrix<double> m{ 12, 34 };
    • creating a random matrix of size 12 X 34, in range [0.0, 1.0]:
    auto rand = feng::rand<double>(12, 34);
    • creating a matrix of size 12 X 34, with all elements to be 0:
    auto zero = feng::zeros<double>(12, 34);
    • creating a matrix of size 12 X 34, with all elements to be 1:
    auto one = feng::ones<double>(12, 34);
    • and ones_like:
    auto another_one = feng::ones_like( one );
    • creating a matrix of size 12 X 34, with all elements to be uninitialized:
    auto one = feng::empty<double>(12, 34);
    • converting value types
    auto one_fp32 = feng::empty<float>(12, 34);
auto one_fp64 = one_fp32.astype<double>(); //
    • creating a matrix of size 1 X (12*34), with all elements from 0 to 12x34, then reshape to 12 X 34:
    auto one = feng::arange<double>(12*34);
one.reshape( 12, 34 );
    • loading matrix from a local txt file ./mat.txt, the delimiter can be either of , ,, \t or ;, the line end is \n:
    feng::matrix<double> mat;
mat.load_txt( './mat.txt' );
    • loading matrix from an opencv matrix instance with interface from_opencv
    cv::Mat M( 2, 2, CV_8UC3, cv::Scalar(0, 0, 255) );
std::cout << "OPENCV matrix:\n" << M << std::endl;
feng::matrix<std::uint32_t> mat;
mat.from_opencv( M );
std::cout << "Converted to feng::matrix:\n" << mat << std::endl;

    This will produce output like:

    OPENCV matrix:
[  0,   0, 255,   0,   0, 255;
   0,   0, 255,   0,   0, 255]
Converted to feng::matrix:
0       0       255     0       0       255
0       0       255     0       0       255

    However, to compile and link the code above, make sure to

    1. define the opencv guard by passing -DOPENCV option to the compiler (g++),
    2. tell the compile where to find the opencv header files, for example, passing pkg-config --cflags opencv4 to the compiler (g++), and
    3. tell the linker which libraries to link against, for example, passing pkg-config --libs opencv4 to the linker (g++). Or
    4. compile and link in a single command, such as
    g++ -o ./test_test  -std=c++20 -Wall -Wextra -fmax-errors=1 -Ofast -flto=auto  -funroll-all-loops -pipe -march=native -DPARALLEL -DOPENCV `pkg-config --cflags opencv4` -Wno-deprecated-enum-enum-conversion `pkg-config --libs opencv4`  -pthread -lstdc++fs -Wl,--gc-sections -flto tests/test.cc

    And to convert a matrix instance to opencv matrix:

    cv::Mat m = mat.to_opencv( 3 );

    in whcih the parameter 3 is for the image channels. By default this parameter is 1, and up to 4 channels are allowed.

  • others

create, row, col, size, shape, clear

feng::matrix<double> m{ 64, 256 };
m.save_as_bmp( "./images/0002_create.bmp" );

assert( m.row() == 64 );
assert( m.col() == 256 );
assert( m.size() == m.row() * m.col() );

auto const [r,c] = m.shape();
assert( r == m.row() );
assert( c == m.col() );

m.clear();
assert( 0 == m.row() );
assert( 0 == m.col() );

create

element access using operator [] or operator ()

feng::matrix<double> m{ 64, 256 };
for ( auto r = 12; r != 34; ++r )
    for ( auto c = 34; c != 45; ++c )
        m[r][c] = 1.0;

for ( auto r = 34; r != 45; ++r )
    for ( auto c = 123; c != 234; ++c )
        m(r, c) = -1.0;

m.save_as_bmp( "./images/0019_create.bmp" );

create

range-based for access

feng::matrix<double> m{ 64, 256 };
int starter = 0;
double const keys[] = { 1, 2, 3, 4, 5, 6, 7, 8 };
for ( auto& x : m )
{
    int val = starter++ & 0x7;
    x = keys[val];
}
m.save_as_bmp( "./images/0000_access.bmp" );

access

copying, resizing and reshaping

feng::matrix<double> m{ 64, 256 };
for ( auto r = 12; r != 34; ++r )
    for ( auto c = 12; c != 34; ++c )
        m[r][c] = 1.0;
m.save_as_bmp( "./images/0020_create.bmp" );

created matrix m:

create

feng::matrix<double> n = m; //copying
n.save_as_bmp( "./images/0021_create.bmp" );

copied matrix n:

create

n.resize( 63, 244 );
n.save_as_bmp( "./images/0022_create.bmp" );

resized matrix n:

create

m.reshape( m.col(), m.row() );
m.save_as_bmp( "./images/0023_create.bmp" );

reshaped matrix m:

create

matrix slicing

feng::matrix<double> m{ 64, 256 };
std::fill( m.upper_diag_begin(1), m.upper_diag_end(1), 1.0 );
std::fill( m.diag_begin(), m.diag_end(), 1.0 );
std::fill( m.lower_diag_begin(1), m.lower_diag_end(1), 1.0 );
m.save_as_bmp( "./images/0000_slicing.bmp" );

matrix slicing

feng::matrix<double> n{ m, 0, 32, 0, 64 };
n.save_as_bmp( "./images/0001_slicing.bmp" );

matrix slicing

feng::matrix<double> p{ m, {16, 48}, {0, 64} };
p.save_as_bmp( "./images/0002_slicing.bmp" );

matrix slicing

meshgrid

meshgrid returns 2-D grid coordinates based on the coordinates contained in interger x and y.

auto const& [X, Y] = feng::meshgrid( 3, 5 );
std::cout << X << std::endl;
std::cout << Y << std::endl;

This will produce

0       1       2
0       1       2
0       1       2
0       1       2
0       1       2

0       0       0
1       1       1
2       2       2
3       3       3
4       4       4

while the code below

auto const& [X, Y] = feng::meshgrid( 384, 512 );
X.save_as_bmp( "./images/0000_meshgrid_x.bmp", "grey" );
Y.save_as_bmp( "./images/0000_meshgrid_y.bmp", "grey" );

generates two images

meshgrid x

meshgrid y

arange

arange<Type>([start, ]stop, [step, ])

Return evenly spaced row matrix within a given interval.

auto m = feng::arange<double>( 256*256 );
m.reshape( 256, 256 );
m.save_as_bmp( "./images/0000_arange.bmp" );

arange 256X256

clip

For an normal matrix m

    feng::matrix<double> m{ 64, 256 };
    std::generate( m.begin(), m.end(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
    m.save_as_bmp( "./images/0000_clip.bmp" );

clip0

it can be transformed to range [0, 1] by applying sin on it

    m = feng::sin(m);
    m.save_as_bmp( "./images/0001_clip.bmp" );

clip1

then this matrix can be clipped to range [0.1, 0.9]

    auto const& cm0 = feng::clip( 0.1, 0.9 )( m );
    cm0.save_as_bmp( "./images/0002_clip.bmp" );

clip2

or even to range [0.4, 0.6]

    auto const& cm1 = feng::clip( 0.4, 0.6 )( m );
    cm1.save_as_bmp( "./images/0003_clip.bmp" );

clip3

iterations

elementwise apply

feng::matrix<double> m{ 64, 256 };
std::generate( m.begin(), m.end(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0000_apply.bmp" );

before apply:

apply

m.apply( [](auto& x) { x = std::sin(x); } );
m.save_as_bmp( "./images/0001_apply.bmp" );

after apply:

apply

iteration from head to tail

feng::matrix<double> m{ 64, 256 };
std::generate( m.begin(), m.end(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0000_create.bmp" );

create

iteration from tail to head

feng::matrix<double> m{ 64, 256 };
std::generate( m.rbegin(), m.rend(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0006_create.bmp" );

create

iteration through a selected row

feng::matrix<double> m{ 64, 256 };
std::generate( m.row_begin(17), m.row_end(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0001_create.bmp" );

create

reverse iteration through a selected row

feng::matrix<double> m{ 64, 256 };
std::generate( m.row_rbegin(17), m.row_rend(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0003_create.bmp" );

create

iteration through a selected column

feng::matrix<double> m{ 64, 256 };
std::generate( m.col_begin(17), m.col_end(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0004_create.bmp" );

create

reverse iteration through a selected column

feng::matrix<double> m{ 64, 256 };
std::generate( m.col_rbegin(17), m.col_rend(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0003_create.bmp" );

create

iteration through diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.diag_begin(), m.diag_end(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0011_create.bmp" );

create

reverse iteration through diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.diag_rbegin(), m.diag_rend(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0012_create.bmp" );

create

iteration through upper diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.upper_diag_begin(17), m.upper_diag_end(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0007_create.bmp" );

create

reverse iteration through upper diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.upper_diag_rbegin(17), m.upper_diag_rend(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0008_create.bmp" );

create

iteration through lower diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.lower_diag_begin(17), m.lower_diag_end(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0009_create.bmp" );

create

reverse iteration through lower diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.lower_diag_rbegin(17), m.lower_diag_rend(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0010_create.bmp" );

create

iteration through anti diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.anti_diag_begin(), m.anti_diag_end(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0017_create.bmp" );

create

reverse iteration through anti diagonal

feng::matrix<double> m{ 64, 256 };
iag_rbegin(), m.anti_diag_rend(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0018_create.bmp" );

create

iterator through upper anti diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.upper_anti_diag_begin(17), m.upper_anti_diag_end(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0013_create.bmp" );

create

reverse iteration through upper anti diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.upper_anti_diag_rbegin(17), m.upper_anti_diag_rend(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0014_create.bmp" );

create

iteration through lower anti diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.lower_anti_diag_begin(17), m.lower_anti_diag_end(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0015_create.bmp" );

create

reverse iteration through lower anti diagonal

feng::matrix<double> m{ 64, 256 };
std::generate( m.lower_anti_diag_rbegin(17), m.lower_anti_diag_rend(17),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0016_create.bmp" );

create

functions

clone -- matrix slicing

feng::matrix<double> m{ 64, 256 };
std::fill( m.diag_begin(), m.diag_end(), 1.1 );
m.save_as_bmp( "./images/0000_clone.bmp" );

matrix m:

clone

auto n = m.clone( 0, 32, 0, 64 );
n.save_as_bmp( "./images/0001_clone.bmp" );

m slicing of [0:32, 0:64]:

clone

n.clone( m, 32, 64, 0, 64 );
n.save_as_bmp( "./images/0002_clone.bmp" );

m slicing of [32:64, 0:64]:

clone

astype

For a normal matrix such like

feng::matrix<double> m{ 64, 256 };
std::generate( m.begin(), m.end(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return std::sin(init); }; }() );
m.save_as_bmp( "./images/0000_astype.bmp" );

There are many colors in its visualization.

astype_0

astype can convert it to 3 colors:

m = m * 2.0;
auto const& mm = m.astype<int>();
mm.save_as_bmp( "./images/0001_astype.bmp" );

astype_1

data -- raw memory access

feng::matrix<double> m{ 64, 256 };
m.save_as_bmp( "./images/0000_data.bmp" );

original matrix

data

auto ptr = m.data();
for (  auto idx = 0UL; idx != m.size(); ++idx )
    ptr[idx] = std::sin( idx*idx*0.1 );
m.save_as_bmp( "./images/0001_data.bmp" );

after modification

data

det -- matrix determinant

feng::matrix<double> m{ 128, 128 };
std::generate( m.diag_begin(), m.diag_end(), [](){ double x = 0.9; return [x]() mutable { x+= 0.156; return x; }(); } );
double det1 = m.det();
double det2 = std::accumulate( m.diag_begin(), m.diag_end(), 1.0, []( double x, double y ){ return x*y; } );
std::cout << det1 << "\t:\t" << det2 << std::endl;

generated output is

1069.00941294551	:	1069.0094129455

operator divide-equal

auto m = feng::rand<double>( 197, 197 );
auto n = m;
n /= 2.0;
m /= n;
m.save_as_bmp( "images/0000_divide_equal.bmp" );

divide equal

matrix inverse

auto const& m = feng::rand<double>( 128, 128 );
auto const& n = m.inverse();
auto const& identity = m * n;
identity.save_as_bmp( "./images/0000_inverse.bmp" );

matrix inverse

save matrix to images with colormap

Here we demonstrate how to save matrix to images with specified colormap. There are 18 builtin colormaps:

  • autumn
  • bluehot
  • bone
  • cool
  • copper
  • default
  • gray
  • hotblue
  • hot
  • hsv
  • jet
  • lines
  • obscure
  • parula
  • pink
  • spring
  • summer
  • winter
  • gray

First we load the matrix from a '.txt' file

feng::matrix<double> m;
m.load_txt( "./images/Lenna.txt" );

Then we can save this matrix to a '.bmp' file with default colormap:

m.save_as_bmp( "./images/0000_save_with_colormap_default.bmp" );

The default image looks like:

colormap-default

m.save_as_bmp( "./images/0000_save_with_colormap_parula.bmp", "parula" );

The parula image looks like:

colormap-parula

m.save_as_bmp( "./images/0000_save_with_colormap_bluehot.bmp", "bluehot" );

The bluehot image looks like:

colormap-bluehot

m.save_as_bmp( "./images/0000_save_with_colormap_hotblue.bmp", "hotblue" );

The hotblue image looks like:

colormap-hotblue

m.save_as_bmp( "./images/0000_save_with_colormap_jet.bmp", "jet" );

The jet image looks like:

colormap-jet

m.save_as_bmp( "./images/0000_save_with_colormap_obscure.bmp", "obscure" );

The obscure image looks like:

colormap-obscure

m.save_as_bmp( "./images/0000_save_with_colormap_gray.bmp", "gray" );

The gray image looks like:

colormap-gray

m.save_as_bmp( "./images/0000_save_with_colormap_hsv.bmp", "hsv" );

The hsv image looks like:

colormap-hsv

m.save_as_bmp( "./images/0000_save_with_colormap_hot.bmp", "hot" );

The hot image looks like:

colormap-hot

m.save_as_bmp( "./images/0000_save_with_colormap_cool.bmp", "cool" );

The cool image looks like:

colormap-cool

m.save_as_bmp( "./images/0000_save_with_colormap_spring.bmp", "spring" );

The spring image looks like:

colormap-spring

m.save_as_bmp( "./images/0000_save_with_colormap_summer.bmp", "summer" );

The summer image looks like:

colormap-summer

m.save_as_bmp( "./images/0000_save_with_colormap_autumn.bmp", "autumn" );

The autumn image looks like:

colormap-autumn

m.save_as_bmp( "./images/0000_save_with_colormap_winter.bmp", "winter" );

The winter image looks like:

colormap-winter

m.save_as_bmp( "./images/0000_save_with_colormap_bone.bmp", "bone" );

The bone image looks like:

colormap-bone

m.save_as_bmp( "./images/0000_save_with_colormap_copper.bmp", "copper" );

The copper image looks like:

colormap-copper

m.save_as_bmp( "./images/0000_save_with_colormap_pink.bmp", "pink" );

The pink image looks like:

colormap-pink

m.save_as_bmp( "./images/0000_save_with_colormap_lines.bmp", "lines" );

The lines image looks like:

colormap-lines

For sparse data such as particles, it is highly recommended to use colormap bluehot:

feng::matrix<double> m;
m.load_txt( "./images/star.txt" );
m.save_as_bmp( "./images/0001_star_bluehot.bmp", "bluehot" );

The starry image generate is demonstrated below:

colormap-starry-bluehot

save load

To load an image from a txt file, we can use .load_txt method:

feng::matrix<double> m;
m.load_txt( "./images/Lenna.txt" );
m.save_as_txt( "./images/0000_save_load.txt" );
m.save_as_binary( "./images/0000_save_load.bin" );
m.save_as_bmp( "./images/0000_save_load.bmp" );

The image loaded is

image saved

feng::matrix<double> n;
n.load_txt( "./images/0000_save_load.txt" );
n.save_as_bmp( "./images/0001_save_load.bmp" );

image saved

n.load_binary( "./images/0000_save_load.bin" );
n.save_as_pgm( "./images/0002_save_load.pgm" );

image saved

load npy

Loading a matrix created by numpy is straightforward:

feng::matrix<double> mat;
mat.load_npy( "./images/64.npy");

save load bmp

To load an image from a bmp file, we can use feng::load_bmp function, which will return an oject of type std::optional<std::array<feng::matrix<std::uint8_t>,3>>:

std::optional<std::array<matrix<std::uint8_t>,3>> load_bmp( std::string const& file_path ) {...}

We first generate an image of Lenna:

feng::matrix<double> m;
m.load_txt( "./images/Lenna.txt" );
m.save_as_bmp( "./images/Lenna.bmp", "gray" );

which looks like:

Lenna

Then we can try to load it directly:

auto const& mat_3 = feng::load_bmp( "./images/Lenna.bmp" );

if successfull, ``mat_3with hold a channel-first image. To accessmat_3` we need to verify it is accessible first by:

if ( mat_3 )
{

Then we can visualize its red channel:

(*mat_3)[0].save_as_bmp( "./images/0001_save_load_julia_red.bmp", "gray" );

red channel

green channel:

(*mat_3)[1].save_as_bmp( "./images/0001_save_load_julia_green.bmp", "gray" );

green channel

and blue channel:

(*mat_3)[2].save_as_bmp( "./images/0001_save_load_julia_blue.bmp", "gray" );

blue channel

save png

It is possible to save matrix as a png file the same way as bmp

m.save_as_png( "./images/0000_save_with_colormap_default.png" );
m.save_as_png( "./images/0000_save_with_colormap_parula.png", "parula" );

png

operator minus equal

    feng::matrix<double> image;
    image.load_txt( "images/Lenna.txt" );
    image.save_as_bmp("images/0000_minus_equal.bmp", "gray");

image minus equal

    double const min = *std::min_element( image.begin(), image.end() );
    image -= min;
    image.save_as_bmp("images/0001_minus_equal.bmp", "jet");

image minus equal

    image -= image;
    image.save_as_bmp("images/0002_minus_equal.bmp");

image minus equal

plot

plot is an alias name of save_as_bmp:

feng::matrix<double> m;
m.load_txt( "./images/Lenna.txt" );
m.plot( "./images/0000_plot_default.bmp" );

default_plot

m.plot( "./images/0000_plot_jet.bmp", "jet" );

default_jet

Juliet set

Having plot/save_as_bmp method implemented, it is convenient to plot a juliet set to an image. The polynomial function is: f_c(z) = z^n + c Suppose we have the ranges of z, the x and y coordinates in a 2D image, the power n and the constant c, we can generate an image like this

auto make_julia_set( std::complex<double> const& lower_left, std::complex<double> const& upper_right, std::complex<double>const& cval, unsigned long const dim = 1024, unsigned long const iterations = 1024, unsigned long const powers = 2 )
{
    std::complex<double> spacing_ratio{ (std::real(upper_right)-std::real(lower_left)) / static_cast<double>(dim),
                                        (std::imag(upper_right)-std::imag(lower_left)) / static_cast<double>(dim) };

    feng::matrix<std::complex<double>> cmat{ dim, dim };
    for ( auto r = 0UL; r != dim; ++r )
        for ( auto c = 0UL; c != dim; ++c )
            cmat[r][c] = std::complex<double>{ static_cast<double>(r), static_cast<double>(c) };

    auto const& converter = [&spacing_ratio, &lower_left, &upper_right, iterations, powers, &cval]( std::complex<double>& c )
    {
        std::complex<double> z = lower_left + std::complex<double>{ std::real(spacing_ratio)*std::real(c), std::imag(spacing_ratio)*std::imag(c) };
        c = std::complex<double>{ static_cast<double>(iterations), 0 };
        for ( auto idx = 0UL; idx != iterations; ++idx )
        {
            z = std::pow( z, powers ) + cval;
            if ( std::abs(z) > 2.0 )
            {
                c = std::complex<double>{ static_cast<double>(idx), 0 };
                break;
            }
        }
    };

    cmat.apply( converter );
    return feng::real(cmat);
}

where the lower_left and upper_right defines the size of the canvas, the cval is for the constant c, the dim gives the sampling density, the iterations gives the maximum value in the sampled positions, and the powers for the maximum order of the polinomial.

Then we are able to generate a series of the fractional images:

unsigned long const n = 32;
for ( unsigned r = 0; r != n; ++r )
    for ( unsigned c = 0; c != n; ++c )
    {
        std::complex<double> zc{ double(r)/n*1.8-0.9, double(c)/n*1.8-0.9 };
        auto&&  mat = make_julia_set( std::complex<double>{-1.5, -1.0}, std::complex<double>{1.5, 1.0}, zc, 1024, 1024, 4 );
        std::string file_name = std::string{"./images/julia_set_4/0001_julia_set_"} + std::to_string(r) + std::string{"-"} + std::to_string(c) + std::string{".bmp"};
        mat.save_as_bmp( file_name, "tealhot" );
    }

A typical result image looks like this:

juliet_set

operator multiply equal

auto m = feng::rand( 127, 127 );
m *= m.inverse();
m.save_as_bmp("images/0001_multiply_equal.bmp");

image multiply equal

operator plus equal

    feng::matrix<double> image;
    image.load_txt( "images/Lenna.txt" );
    image.save_as_bmp("images/0000_plus_equal.bmp", "gray");

image plus equal

    double const mn = *std::min_element( image.begin(), image.end() );
    double const mx = *std::max_element( image.begin(), image.end() );
    image = (image - mn)/(mx - mn);

    auto const& noise = feng::rand<double>( image.row(), image.col() );
    image += 0.1*noise;
    image.save_as_bmp("images/0001_plus_equal.bmp", "gray");

image plus equal

operator prefix

auto const& m = feng::random<double>( 127, 127 );
auto const& pp = +m;
auto const& pm = -m;
auto const& shoule_be_zero = pp + pm;
shoule_be_zero.save_as_bmp("images/0000_prefix.bmp");

image prefix

elementwise mathematical functions

Most common mathematical functions are supported as elementwise matrix operator. Here only sin and sinh are demonstrated

elementwise sin

feng::matrix<double> m{ 64, 256 };
std::generate( m.begin(), m.end(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init; }; }() );
m.save_as_bmp( "./images/0000_sin.bmp" );

sin image orig

m = feng::sin(m);
m.save_as_bmp( "./images/0001_sin.bmp" );

sin image after

elementwise sinh

feng::matrix<double> m{ 64, 256 };
std::generate( m.begin(), m.end(),  [](){ double init = 0.0; return [init]() mutable { init += 0.1; return init/500.0; }; }() );
m.save_as_bmp( "./images/0000_sinh.bmp" );

sinh image orig

m = feng::sinh(m);
m.save_as_bmp( "./images/0001_sinh.bmp" );

sinh image after

common functions

eye function

auto const& m = feng::eye<double>( 128, 128 );
m.save_as_bmp( "./images/0000_eye.bmp" );

eye image

linspace

auto const& m = feng::linspace<double>( 1, 10, 10 );
std::cout << "linspace<double>(1, 10, 10):\n" << m << std::endl;

gives out an array of size 1 x 10:

linspace(1, 10, 10): 1 2 3 4 5 6 7 8 9 10

auto const& m = feng::linspace<double>( 1, 10, 10, false );
std::cout << "linspace<double>(1, 10, 10, false):\n" << m << std::endl;

gives out an array of size 1 x 9:

linspace(1, 10, 10, false): 1 1.89999999999999991 2.79999999999999982 3.69999999999999973 4.59999999999999964 5.5 6.40000000000000036 7.30000000000000071 8.20000000000000107 9.10000000000000142

And the prototype of linspace is:

matrix<T> linspace( T start, T stop, const std::uint_least64_t num = 50ULL, bool end_point=true )

magic function

Calling magic method is quite straightforward:

std::cout << "Magic 3\n" << feng::magic( 3 ) << std::endl;
std::cout << "Magic 4\n" << feng::magic( 4 ) << std::endl;
std::cout << "Magic 5\n" << feng::magic( 5 ) << std::endl;
std::cout << "Magic 6\n" << feng::magic( 6 ) << std::endl;

This will produce a series of magic matrices:


Magic 3
 8      1       6
3       5       7
4       9       2

Magic 4
 16     3       2       13
5       10      11      8
9       6       7       12
4       15      14      1

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

Magic 6
 32     29      4       1       24      21
30      31      2       3       22      23
12      9       17      20      28      25
10      11      18      19      26      27
13      16      33      36      8       5
14      15      34      35      6       7

Also we can expand it a bit to do a better visualization:

unsigned long n = 38;
unsigned long pixs = 16;

auto const& mat = feng::magic( n );

feng::matrix<double> v_mat( n*pixs, n*pixs );

for ( auto r = 0UL; r != n; ++r )
    for ( auto c = 0UL; c != n; ++c )
        for ( auto rr = 0UL; rr != pixs; ++rr )
            for ( auto cc = 0UL; cc != pixs; ++cc )
                v_mat[r*pixs+rr][c*pixs+cc] = mat[r][c];

v_mat.save_as_bmp("./images/0001_magic.bmp");

This produces an image looks like:

magic image

matrix convolution

feng::matrix<double> m;
m.load_txt( "./images/Lenna.txt" );
m.save_as_bmp( "./images/0000_conv.bmp", "gray" );

convolution 1

Full 2D matrix convolution with zero-paddings is given by conv or conv2 :

feng::matrix<double> filter{3, 3, {0.0, 1.0, 0.0,
                                   1.0,-4.0, 1.0,
                                   0.0, 1.0, 0.0}};
auto const& edge = feng::conv( m, filter );
edge.save_as_bmp( "./images/0001_conv.bmp", "gray" );

convolution 2

The convolution has three modes: valid, same, and full

The valid mode gives out convolution result without zero-paddings:

auto const& edge_valid = feng::conv( m, filter, "valid" );
edge_valid.save_as_bmp( "./images/0001_conv_valid.bmp", "gray" );

convolution valid

The same mode returns the central part of the convolution result with zero-paddings, of the same size as the larger matrix passed to function conv

auto const& edge_same = feng::conv( m, filter, "same" );
edge_same.save_as_bmp( "./images/0001_conv_same.bmp", "gray" );

convolution same

full mode is the default mode:

auto const& edge_full = feng::conv( m, filter, "full" );
edge_full.save_as_bmp( "./images/0001_conv_full.bmp", "gray" );

convolution full

make view function

feng::matrix<double> m;
m.load_txt( "./images/Lenna.txt" );
m.save_as_bmp( "./images/0000_make_view.bmp" );

make view 1

auto const[r,c] = m.shape();
auto const& v = feng::make_view( m, {r>>2, (r>>2)*3}, {c>>2, (c>>2)*3} );

v.save_as_bmp( "./images/0001_make_view.bmp" );

make view 2

And creating new matrix from a view

auto new_matrix{v};
new_matrix.save_as_bmp( "./images/0002_make_view.bmp" );

make view 3

A matrix view has several methods, row(), col(), shape() and operator[]():

feng::matrix<double> n{ v.row(), v.col() }; // row() and col() of a matrix view
for ( auto r = 0UL; r != n.row(); ++r )
    for ( auto c = 0UL; c != n.col(); ++c )
        n[r][c] = v[r][c]; // accessing matrix elements using operator [], read-only
n.save_as_bmp( "./images/0003_make_view.bmp", "gray" );

make view 4

singular value decomposition

We first Load an image from hardisk and normalize it to range [0,1]

// load
feng::matrix<double> m;
m.load_txt( "./images/Teacher.txt" );
// normalize
auto const mx = *std::max_element( m.begin(), m.end() );
auto const mn = *std::min_element( m.begin(), m.end() );
m = ( m - mn ) / ( mx - mn + 1.0e-10 );
// take a snapshot
m.save_as_bmp( "./images/0000_singular_value_decomposition.bmp", "gray" );

This image looks like:

svd_1

Then we add some white noise to remove possible singularity in it:

// adding noise
auto const[r, c] = m.shape();
m += feng::rand<double>( r, c );
// record noisy matrix
m.save_as_bmp( "./images/0001_singular_value_decomposition.bmp", "gray" );

The noisy image now looks like

svd_2

We execute Singular Value Decomposition by calling function std::optional<std::tuple<matrix, matrix, matrix>> singular_value_decomposition( matrix const& ), or svd

// execute svd
auto const& svd = feng::singular_value_decomposition( m );

If the svd is successfully, we can verify the accuricy by reconstructing the noisy image by matrix multiplications

// check svd result
if (svd) // case successful
{
	// extracted svd result matrices, u, v w
	auto const& [u, v, w] = (*svd);
	// try to reconstruct matrix using  u * v * w'
	auto const& m_ = u * v * (w.transpose());
	// record reconstructed matrix
	m_.save_as_bmp( "./images/0002_singular_value_decomposition.bmp", "gray" );

The reconstructed image looks like:

svd_3

One interesting application of SVD is data compression. In the code above, we use full rank to restore the original noisy image. However, we can select only 1/2 or 1/4 or even less ranks to approximate the original image. The code below demonstrates how:

	auto dm = std::min( r, c );
	auto factor = 2UL;
	while ( dm >= factor )
	{
		auto new_dm = dm / factor;

		feng::matrix<double> const new_u{ u, std::make_pair(0UL, r), std::make_pair(0UL, new_dm) };
		feng::matrix<double> const new_v{ v, std::make_pair(0UL, new_dm), std::make_pair(0UL, new_dm) };
		feng::matrix<double> const new_w{ w, std::make_pair(0UL, c), std::make_pair(0UL, new_dm) };

		auto const& new_m = new_u * new_v * new_w.transpose();

		new_m.save_as_bmp( "./images/0003_singular_value_decomposition_"+std::to_string(new_dm)+".bmp", "gray" );

		factor *= 2UL;
	}
}
else
{
	std::cout << "Failed to execute Singular Value Decomposition for this matrix!\n";
}

When using 256 ranks, the reconstructed image lookes like:

svd_4

When using 128 ranks, the reconstructed image lookes like:

svd_4

When using 64 ranks, the reconstructed image lookes like:

svd_4

When using 32 ranks, the reconstructed image lookes like:

svd_4

When using 16 ranks, the reconstructed image lookes like:

svd_4

When using 8 ranks, the reconstructed image lookes like:

svd_4

When using 4 ranks, the reconstructed image lookes like:

svd_4

When using 2 ranks, the reconstructed image lookes like:

svd_4

When using 1 ranks, the reconstructed image lookes like:

svd_4

pooling

We are able to pooling an image with function pooling( matrix, dim_row, dim_col, option ), where option can be either of mean, max, or min, if no option provided, then mean is applied.

For an normal image

feng::matrix<double> m;
m.load_txt( "./images/Lenna.txt" );
m.save_as_bmp( "./images/0000_pooling.bmp", "gray" );

Lenna

A 2X2 mean pooling looks like:

auto const& pooling_2 = feng::pooling( m, 2 );
pooling_2.save_as_bmp( "./images/0000_pooling_2.bmp", "gray" );

Lenna pooling 2

And a 4X4 pooling is

auto const& pooling_4 = feng::pooling( m, 4 );
pooling_4.save_as_bmp( "./images/0000_pooling_4.bmp", "gray" );

Lenna pooling 4

For hyterdyne pooling of 2X4

auto const& pooling_2_4 = feng::pooling( m, 2, 4 );
pooling_2_4.save_as_bmp( "./images/0000_pooling_2_4.bmp", "gray" );

Lenna pooling 2 4

And 4X2

auto const& pooling_4_2 = feng::pooling( m, 4, 2 );
pooling_4_2.save_as_bmp( "./images/0000_pooling_4_2.bmp", "gray" );

Lenna pooling 4_2

Also min pooling is possible

auto const& pooling_min = feng::pooling( m, 2, "min" );
pooling_min.save_as_bmp( "./images/0000_pooling_2_min.bmp", "gray" );

Lenna pooling min

And max pooling

auto const& pooling_max = feng::pooling( m, 2, "max" );
pooling_max.save_as_bmp( "./images/0000_pooling_2_max.bmp", "gray" );

Lenna pooling max

gauss jordan elimination

for a random matrix

auto const& m = feng::rand<double>( 64, 128);
m.save_as_bmp( "./images/0000_gauss_jordan_elimination.bmp", "gray" );

gauss_jordan_elimination_0

auto const& n = feng::gauss_jordan_elimination( m ); //<- also `feng::rref(m);`, alias name from Matlab

if (n)
    (*n).save_as_bmp( "./images/0001_gauss_jordan_elimination.bmp", "gray" );
else
    std::cout << "Failed to execute Gauss-Jordan Elimination for matrix m.\n";

after applying Gauss Jordan elimination, the matrix is reduced to a form of

gauss_jordan_elimination_1

lu decomposition

We load a Lena from harddisk file ./images/Lenna.txt using function load_txt(std::string const&):

    // initial matrix
    feng::matrix<double> m;
    m.load_txt( "./images/Lenna.txt" );
    m.save_as_bmp( "./images/0000_lu_decomposition.bmp", "gray" );

The loaded image lookes like below:

lu_0

Then we scale this image to range [0,1] and add some uniform random noise to it (to remove singularity of the original image)

    // adding noise
    double mn = *std::min_element( m.begin(), m.end() );
    double mx = *std::max_element( m.begin(), m.end() );
    m = (m-mn) / (mx - mn + 1.0e-10);
    auto const& [row, col] = m.shape();
    m += feng::rand<double>( row, col );
    m.save_as_bmp( "./images/0001_lu_decomposition.bmp", "gray" );

The noised image lookes like:

lu_1

we can do LU decomposition simply with std::optional<std::tuple<matrix, matrix>>lu_decomposition( matrix const& ) function

    // lu decomposition
    auto const& lu = feng::lu_decomposition( m );
    if (lu)
    {
        auto const& [l, u] = lu.value();
        l.save_as_bmp( "./images/0002_lu_decomposition.bmp", "jet" );

the result of the LU decompositon is a Maybe monad of a tuple of two matrices, i.e., std::optional<std::tuple<feng::matrix<Type, Allocator>, feng::matrix<Type, Allocator>>>, therefor, we need to check its value before using it.

Then we can draw the lower matrix L

lu_2

        u.save_as_bmp( "./images/0003_lu_decomposition.bmp", "jet" );

And the upper matrix U

lu_3

        auto const& reconstructed = l * u;
        reconstructed.save_as_bmp( "./images/0004_lu_decomposition.bmp", "gray" );

We can multiply L and U back to see if the decomposition is correct or not.

lu_4

    }
    else
    {
        std::cout << "Error: Failed to execute lu decomposition for matrix m!\n";
    }

A typical use of LU Decomposition is to solve an equation in the form of Ax=b, this is done by calling auto const& x = feng::lu_solver(A,b), and, again, the returned value is a Maybe monad, std::optional<matrix> lu_solver( matrix const&, matrix const& ), therefore we need to check its value before using it.

    auto const X = feng::rand<double>( row, 1 );
    auto const b = m * X;
    auto const& ox = feng::lu_solver( m, b );

    if (ox)
    {
        auto const& diff = ox.value() - X;
        auto const mae = std::sqrt( std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0) / diff.size() );
        std::cout << "mean absolute error for lu solver is " << mae << "\n";
    }
    else
    {
        std::cout << "Error: Failed to solve equation with lu solver!\n";
    }

And we can also evaluate the solver's accuracy with the mean absolute value error; the output as small as

mean absolute error for lu solver is 2.34412875135465e-10

License

Copyright <2018> <Feng Wang>

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Dependency

This library only depends on a C++-20 standard compiler.

Installation

This is a single-file header-only library. Put matrix.hpp directly into the project source tree or somewhere reachable from your project.

Synopsis

Building tests and examples

Simple execute make or make test or make example at the root folder.

Notes and references

Design

Requirements

In the table below, M denotes a matrix of type T using allocator of type A, m is a value of type A, a and b denote values of type M, u denotes an identifier, r denotes a non-constant value of type M, and rv denotes a non-const rvalue of type M.

Expression ReturnType Operational Sematics Assertion, pre-/post-condition Complexity
M::value_type T T is Erasable from M compile time
M::reference T& compile time
M::const_reference T cosnt& compile time
M::difference_type signed integer type identical to the difference type of M::iterator and M::const_iterator compile time
M::size_type unsigned integer type any type that can represent non-negative value of difference_type compile time
M::allocator_type A M::allocator_type::value_type is identical to M::value_type compile time
M::iterator iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements, convertible to M::const_iterator compile time
M::reverse_iterator iterator type whose value type is T reverse_iterator<iterator> compile time
M::const_iterator constant iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements compile time
M::const_reverse_iterator constant iterator type whose value type is T reverse_iterator<const_iterator> compile time
M::row_iterator iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements, convertible to M::const_row_iterator compile time
M::reverse_row_iterator iterator type whose value type is T reverse_iterator<row_iterator> compile time
M::const_row_iterator constant iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements compile time
M::const_reverse_row_iterator constant iterator type whose value type is T reverse_iterator<const_row_iterator> compile time
M::col_iterator iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements, convertible to M::const_col_iterator compile time
M::reverse_col_iterator iterator type whose value type is T reverse_iterator<col_iterator> compile time
M::const_col_iterator constant iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements compile time
M::const_reverse_col_iterator constant iterator type whose value type is T reverse_iterator<const_col_iterator> compile time
M::diag_iterator iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements, convertible to M::const_diag_iterator compile time
M::reverse_diag_iterator iterator type whose value type is T reverse_iterator<diag_iterator> compile time
M::const_diag_iterator constant iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements compile time
M::const_reverse_diag_iterator constant iterator type whose value type is T reverse_iterator<const_diag_iterator> compile time
M::diag_iterator iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements, convertible to M::const_diag_iterator compile time
M::reverse_diag_iterator iterator type whose value type is T reverse_iterator<diag_iterator> compile time
M::const_diag_iterator constant iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements compile time
M::const_reverse_diag_iterator constant iterator type whose value type is T reverse_iterator<const_diag_iterator> compile time
M::anti_diag_iterator iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements, convertible to M::const_anti_diag_iterator compile time
M::reverse_anti_diag_iterator iterator type whose value type is T reverse_iterator<anti_diag_iterator> compile time
M::const_anti_diag_iterator constant iterator type whose value type is T any iterator category that meets the RandomAccessIterator requirements compile time
M::const_reverse_anti_diag_iterator constant iterator type whose value type is T reverse_iterator<const_anti_diag_iterator> compile time

Header <matrix> synopsis

namespace xxx
{

    template < typename T, class Allocator = allocator<T> >
    struct matrix
    {
    	//types:
        typedef T 								                        value_type;
        typedef Allocator                                               allocator_type;
        typedef value_type& 					                        reference;
        typedef value_type const&                                       const_reference;
        typedef implementation-defined                                  iterator;
        typedef implementation-defined                                  const_iterator;
        typedef implementation-defined                                  row_iterator;
        typedef implementation-defined                                  const_row_iterator;
        typedef implementation-defined                                  col_iterator;
        typedef implementation-defined                                  const_col_iterator;
        typedef implementation-defined                                  diag_iterator;
        typedef implementation-defined                                  const_diag_iterator;
        typedef implementation-defined                                  anti_diag_iterator;
        typedef implementation-defined                                  const_anti_diag_iterator;
        typedef implementation-defined                                  size_type;
        typedef implementation-defined                                  difference_type;
        typedef typename allocator_trait<allocator_type>::pointer       pointer;
        typedef typename allocator_trait<allocator_type>::const_pointer const_pointer;
        typedef std::reverse_iterator<iterator>                         reverse_iterator;
        typedef std::reverse_iterator<const_iterator>                   const_reverse_iterator;
        typedef std::reverse_iterator<row_iterator>                     row_reverse_iterator;
        typedef std::reverse_iterator<const_row_iterator>               const_row_reverse_iterator;
        typedef std::reverse_iterator<col_iterator>                     col_reverse_iterator;
        typedef std::reverse_iterator<const_col_iterator>               const_col_reverse_iterator;
        typedef std::reverse_iterator<diag_iterator>                    diag_reverse_iterator;
        typedef std::reverse_iterator<const_diag_iterator>              const_diag_reverse_terator;
        typedef std::reverse_iterator<anti_diag_iterator>               anti_diag_reverse_iterator;
        typedef std::reverse_iterator<const_anti_diag_iterator>         const_anti_diag_reverse_terator;


        // construct, copy and destroy
        matrix() noexcept;
        explicit matrix ( allocator_type const& ) noexcept;
        explicit matrix ( size_type row_, size_type col_, allocator_type const& = Allocator() );
        matrix( matrix const& );
        matrix( matrix&& ) noexcept;
        matrix( matrix const&, allocator_type const& );
        matrix( matrix&&, allocator_type const& );
        ~matrix();

        matrix& operator = ( matrix const& );
        matrix& operator = ( matrix && ) noexcept( allocator_traits<allocator_type>::propagate_on_container_move_assignment::value || allocator_traits<allocator_type>::is_always_equal::value );

        allocator_type get_allocator() const noexcept;


        //iterators
        iterator                                begin() noexcept;
        const_iterator                          begin() const noexcept;
        iterator                                end() noexcept;
        const_iterator                          end() const noexcept;
        row_iterator                            row_begin(size_type) noexcept;
        const_row_iterator                      row_begin(size_type) const noexcept;
        row_iterator                            row_end(size_type) noexcept;
        const_row_iterator                      row_end(size_type) const noexcept;
        col_iterator                            col_begin(size_type) noexcept;
        const_col_iterator                      col_begin(size_type) const noexcept;
        col_iterator                            col_end(size_type) noexcept;
        const_col_iterator                      col_end(size_type) const noexcept;
        diag_iterator                           diag_begin(difference_type) noexcept;
        const_diag_iterator                     diag_begin(difference_type) const noexcept;
        diag_iterator                           diag_end(difference_type) noexcept;
        const_diag_iterator                     diag_end(difference_type) const noexcept;

        //reverse iterators
        reverse_iterator                        rbegin() noexcept;
        const_reverse_iterator                  rbegin() const noexcept;
        reverse_iterator                        rend() noexcept;
        const_reverse_iterator                  rend() const noexcept;
        row_reverse_iterator                    row_rbegin(size_type) noexcept;
        const_row_reverse_iterator              row_rbegin(size_type) const noexcept;
        row_reverse_iterator                    row_rend(size_type) noexcept;
        const_row_reverse_iterator              row_rend(size_type) const noexcept;
        col_reverse_iterator                    col_rbegin(size_type) noexcept;
        const_col_reverse_iterator              col_rbegin(size_type) const noexcept;
        col_reverse_iterator                    col_rend(size_type) noexcept;
        const_col_reverse_iterator              col_rend(size_type) const noexcept;
        diag_reverse_iterator                   diag_rbegin(difference_type) noexcept;
        const_diag_reverse_iterator             diag_rbegin(difference_type) const noexcept;
        diag_reverse_iterator                   diag_rend(difference_type) noexcept;
        const_diag_reverse_iterator             diag_rend(difference_type) const noexcept;

        //const iterators
        const_iterator                          cbegin() const noexcept;
        const_iterator                          cend() const noexcept;
        const_row_iterator                      row_cbegin(size_type) const noexcept;
        const_row_iterator                      row_cend(size_type) const noexcept;
        const_col_iterator                      col_cbegin(size_type) const noexcept;
        const_col_iterator                      col_cend(size_type) const noexcept;
        const_diag_iterator                     diag_cbegin(difference_type) const noexcept;
        const_diag_iterator                     diag_cend(difference_type) const noexcept;
        const_reverse_iterator                  crbegin() const noexcept;
        const_reverse_iterator                  crend() const noexcept;
        const_row_reverse_iterator              row_crbegin(size_type) const noexcept;
        const_row_reverse_iterator              row_crend(size_type) const noexcept;
        const_col_reverse_iterator              col_crbegin(size_type) const noexcept;
        const_col_reverse_iterator              col_crend(size_type) const noexcept;
        const_diag_reverse_iterator             diag_crbegin(difference_type) const noexcept;
        const_diag_reverse_iterator             diag_crend(difference_type) const noexcept;

        //capacity and shape
        size_type                               size() const noexcept;
        size_type                               row() const noexcept;
        size_type                               col() const noexcept;
        void                                    resize( size_type row_, size_type col_ );
        void                                    resize( size_type row_, size_type col_, value_type const& value_ );
        void                                    reshape( size_type row_, size_type col_ );
        void                                    transpose();


        //element access
        iterator                                operator[](size_type row_ );
        const_iterator                          operator[](size_type row_ ) const; //TODO
        reference                               operator()( size_type row_, size_type col_ );
        const_reference                         operator()( size_type row_, size_type col_ ) const;
        reference                               at( size_type row_, size_type col_ );
        const_reference                         at( size_type row_, size_type col_ ) const;


        //data access
        pointer                                 data() noexcept;
        const_pointer                           data() const noexcept;

        //modifiers
        void                                    clean() noexcept;
        void                                    swap( matrix& ) noexcept( allocator_traits<allocator_type>::propagate_on_container_move_assignment::value || allocator_traits<allocator_type>::is_always_equal::value );

        //loader, saver and ploter
        void                                    load( string const& file_name_ );
        void                                    save_as( string const& file_name_ );
        void                                    plot( string const& file_name_ ) const;
        void                                    plot( string const& file_name_, string const& builtin_color_scheme_name_ ) const;

        //unary operators
        matrix const                            operator+() const;
        matrix const                            operator-() const;
        matrix const                            operator~() const;
        matrix const                            operator~() const; //bool only

        //computed assignment
        //TODO: return optional<matrix>?
        matrix&                                 operator*=( matrix const& );
        matrix&                                 operator/=( matrix const& );
        matrix&                                 operator+=( matrix const& );
        matrix&                                 operator-=( matrix const& );
        matrix&                                 operator*=( value_type const );
        matrix&                                 operator/=( value_type const );
        matrix&                                 operator%=( value_type const );
        matrix&                                 operator+=( value_type const );
        matrix&                                 operator-=( value_type const );
        matrix&                                 operator^=( value_type const );
        matrix&                                 operator&=( value_type const );
        matrix&                                 operator|=( value_type const );
        matrix&                                 operator<<=( value_type const );
        matrix&                                 operator>>=( value_type const );


        //basic numeric operations
        value_type                              det() const;
        value_type                              tr() const;


    };

    //building functions
    template< typename T, typename A > matrix<T,A> const make_eye( size_type row_, size_type col_, A const& alloc_ );
    template< typename T > matrix<T,std::allocator<T>> const make_eye( size_type row_, size_type col_ );
    template< typename T, typename A > matrix<T,A> const make_zeros( size_type row_, size_type col_, A const& alloc_ );
    template< typename T > matrix<T,std::allocator<T>> const make_zeros( size_type row_, size_type col_ );
    template< typename T, typename A > matrix<T,A> const make_ones( size_type row_, size_type col_, A const& alloc_ );
    template< typename T > matrix<T,std::allocator<T>> const make_ones( size_type row_, size_type col_ );
    template< typename T, typename A > matrix<T,A> const make_diag( matrix<T,A> const );
    template< typename T, typename A > matrix<T,A> const make_triu( matrix<T,A> const );
    template< typename T, typename A > matrix<T,A> const make_tril( matrix<T,A> const );
    template< typename T, typename A > matrix<T,A> const make_rand( size_type row_, size_type col_, A const& alloc_ );
    template< typename T > matrix<T,std::allocator<T>> const make_rand( size_type row_, size_type col_ );
    template< typename T, typename A > matrix<T,A> const make_hilb( size_type n_, A const& alloc_ );
    template< typename T > matrix<T,std::allocator<T>> const make_hilb( size_type n_ );
    template< typename T, typename A > matrix<T,A> const make_magic( size_type n_, A const& alloc_ );
    template< typename T > matrix<T,std::allocator<T>> const make_magic( size_type n_ );
    template< typename T, typename A, typename Input_Itor_1, typename Input_Iterator_2 > matrix<T, A> const make_toeplitz( Input_Iterator_1 begin_, Input_Iterator_1 end_, Input_Iterator_2 begin_2_, A const alloc_ );
    template< typename T, typename A, typename Input_Itor > matrix<T, A> const make_toeplitz( Input_Iterator begin_, Input_Iterator end_, A const alloc_ );
    template< typename T, typename Input_Itor_1, typename Input_Iterator_2 > const matrix<T, std::allocator<T> > make_toeplitz( Input_Iterator_1 begin_, Input_Iterator_1 end_, Input_Iterator_2 begin_2_ ):
    template< typename T, typename Input_Itor > matrix<T, std::allocator<T>> const make_toeplitz( Input_Iterator begin_, Input_Iterator end_ );
    template< typename T, typename A > matrix<T,A> const make_horizontal_cons( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<T,A> const make_vertical_cons( matrix<T,A> const&, matrix<T,A> const& );



    //binary operation
    //TODO: return optional<matrix>?
    template< typename T, typename A > matrix<T,A> const operator * ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<T,A> const operator * ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<T,A> const operator * ( T const&, matrix<T,A> const& );

    template< typename T, typename A > matrix<T,A> const operator / ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<T,A> const operator / ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<T,A> const operator / ( T const&, matrix<T,A> const& );

    template< typename T, typename A > matrix<T,A> const operator + ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<T,A> const operator + ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<T,A> const operator + ( T const&, matrix<T,A> const& );

    template< typename T, typename A > matrix<T,A> const operator - ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<T,A> const operator - ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<T,A> const operator - ( T const&, matrix<T,A> const& );

    template< typename T, typename A > matrix<T,A> const operator % ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<T,A> const operator ^ ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<T,A> const operator & ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<T,A> const operator | ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<T,A> const operator << ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<T,A> const operator >> ( matrix<T,A> const&, T const& );

    template< typename T, typename A > std::ostream& operator << ( std::ostream&, matrix<T,A> const& );
    template< typename T, typename A > std::istream& operator << ( std::istream&, matrix<T,A> const& );

    //logical
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator == ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator == ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator == ( T const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator != ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator != ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator != ( T const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator >= ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator >= ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator >= ( T const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator <= ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator <= ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator <= ( T const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator > ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator > ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator > ( T const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator < ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator < ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator < ( T const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator && ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator && ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator && ( T const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator || ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator || ( matrix<T,A> const&, T const& );
    template< typename T, typename A > matrix<bool,std::allocator_traits<A>::rebind_alloc<bool> > const operator || ( T const&, matrix<T,A> const& );


    //numeric functions
    template< typename T, typename A, typename F > matrix<T,A> const element_wise_apply ( matrix<T,A> const&, F const& f_ );

    //Linear Equations

    //mldivide 	Solve systems of linear equations Ax = B for x
    template< typename T, typename A, typename F > matrix<T,A> const mldivide ( matrix<T,A> const&, matrix<T,A> const& );
    //mrdivide 	Solve systems of linear equations xA = B for x
    template< typename T, typename A, typename F > matrix<T,A> const mrdivide ( matrix<T,A> const&, matrix<T,A> const& );
    //linsolve 	Solve linear system of equations
    template< typename T, typename A, typename F > matrix<T,A> const mrdivide ( matrix<T,A> const&, matrix<T,A> const& );
    //inv 	Matrix inverse
    template< typename T, typename A, typename F > matrix<T,A> const inv ( matrix<T,A> const& );
    //pinv 	Moore-Penrose pseudoinverse of matrix
    template< typename T, typename A, typename F > matrix<T,A> const pinv ( matrix<T,A> const& );
    //lscov 	Least-squares solution in presence of known covariance
    template< typename T, typename A, typename F > matrix<T,A> const lscov ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A, typename F > matrix<T,A> const lscov ( matrix<T,A> const&, matrix<T,A> const&, matrix<T,A> const& );
    //lsqnonneg 	Solve nonnegative linear least-squares problem
    template< typename T, typename A, typename F > matrix<T,A> const lsqnonneg ( matrix<T,A> const&, matrix<T,A> const& );
    template< typename T, typename A, typename F > matrix<T,A> const lsqnonneg ( matrix<T,A> const&, matrix<T,A> const&, matrix<T,A> const& );
    //sylvester 	Solve Sylvester equation AX + XB = C for X
    template< typename T, typename A, typename F > matrix<T,A> const sylvester ( matrix<T,A> const&, matrix<T,A> const&, matrix<T,A> const& );

    //Eigenvalues ans Singular Values

    //eig 	Eigenvalues and eigenvectors
    template< typename T, typename A, typename F > std::tuple<matrix<T,A>,matrix<T,A>> const eig( matrix<T,A> const& );
    //eigs 	Subset of eigenvalues and eigenvectors -- TODO
    //balance 	Diagonal scaling to improve eigenvalue accuracy
    //svd 	Singular value decomposition
    template< typename T, typename A, typename F > std::tuple<matrix<T,A>,matrix<T,A>,matrix<T,A>> const svd( matrix<T,A> const& );
    //svds 	Subset of singular values and vectors -- TODO
    //gsvd 	Generalized singular value decomposition
    template< typename T, typename A, typename F > std::tuple<matrix<T,A>,matrix<T,A>,matrix<T,A>,matrix<T,A>> const gsvd( matrix<T,A> const&, matrix<T,A> const& );
    //ordeig 	Eigenvalues of quasitriangular matrices
    template< typename T, typename A, typename F > matrix<T,A> const ordeig( matrix<T,A> const& );
    template< typename T, typename A, typename F > matrix<T,A> const ordeig( matrix<T,A> const&, matrix<T,A> const& );



}

//linear algebra
/*
ordqz 	Reorder eigenvalues in QZ factorization
ordschur 	Reorder eigenvalues in Schur factorization
polyeig 	Polynomial eigenvalue problem
qz 	QZ factorization for generalized eigenvalues
hess 	Hessenberg form of matrix
schur 	Schur decomposition
rsf2csf 	Convert real Schur form to complex Schur form
cdf2rdf 	Convert complex diagonal form to real block diagonal form
lu 	LU matrix factorization
ldl 	Block LDL' factorization for Hermitian indefinite matrices
chol 	Cholesky factorization
cholupdate 	Rank 1 update to Cholesky factorization
qr 	Orthogonal-triangular decomposition
qrdelete 	Remove column or row from QR factorization
qrinsert 	Insert column or row into QR factorization
qrupdate 	Rank 1 update to QR factorization
planerot 	Givens plane rotation
transpose 	Transpose vector or matrix
ctranspose 	Complex conjugate transpose
mtimes 	Matrix Multiplication
mpower 	Matrix power
sqrtm 	Matrix square root
expm 	Matrix exponential
logm 	Matrix logarithm
funm 	Evaluate general matrix function
kron 	Kronecker tensor product
cross 	Cross product
dot 	Dot product
bandwidth 	Lower and upper matrix bandwidth
tril 	Lower triangular part of matrix
triu 	Upper triangular part of matrix
isbanded 	Determine if matrix is within specific bandwidth
isdiag 	Determine if matrix is diagonal
ishermitian 	Determine if matrix is Hermitian or skew-Hermitian
issymmetric 	Determine if matrix is symmetric or skew-symmetric
istril 	Determine if matrix is lower triangular
istriu 	Determine if matrix is upper triangular
norm 	Vector and matrix norms
normest 	2-norm estimate
cond 	Condition number with respect to inversion
condest 	1-norm condition number estimate
rcond 	Reciprocal condition number
condeig 	Condition number with respect to eigenvalues
det 	Matrix determinant
null 	Null space
orth 	Orthonormal basis for range of matrix
rank 	Rank of matrix
rref 	Reduced row echelon form (Gauss-Jordan elimination)
trace 	Sum of diagonal elements
subspace 	Angle between two subspaces
cosm
sinm
tanm
ctanm
acosm
asinm
atanm
atan2m
coshm
sinhm
tanhm
acosh
asinh
atanh
expm -- expm(A) and expm(A, B)
logm
log10m
exp2m
log2m
logm
powm
sqrtm
cbrtm
hypotm
erfm
erfcm
tgammam
lgammam


//trigonometric
sin
cos
tan
cot
sec
csc

arcsin
arccos
arctan
arccot
arcsec
arccsc


//hyperbolic
sinh
cosh
tanh
coth
sech
csch
arsinh
arcosh
artanh
arsech
arcsch
arcoth
*/


//Fourier Analysis and Filtering
/*
fft
fftshift
ifft
ifftshift
conv
filter
ss2tf
*/

About

A modern, C++20-native, single-file header-only dense 2D matrix library.

Topics

c-plus-plus matrix numeric header-only single-file cpp20

Resources

Readme
Activity

Stars

82 stars

Watchers

3 watching

Forks

8 forks
Report repository

Releases 3

v20230402 Latest
Apr 2, 2023
+ 2 releases

Packages

No packages published

Contributors 3

  •  
  •  
  •  

Languages