mat A = rand<mat>(5,5);
double x = A(1,2);
mat B = A + A;
mat C = A * B;
mat D = A % B;
B.zeros();
B.set_size(10,10);
B.zeros(10,10);

cx_mat X(A,B);

mat A(5,5);
A = 123.0;

mat A(5,5);
A.fill(123.0);

mat A = 123.0 * ones<mat>(5,5);

colvec x(10);
colvec y = zeros<colvec>(10,1);

mat    A = rand<mat>(10,10);
colvec z = A.col(5); // extract a column vector

colvec q(5);
q = 123.0;

colvec q(5);
q.fill(123.0);

colvec q = 123.0 * ones<colvec>(5,1);

rowvec x(10);
rowvec y = zeros<mat>(1,10);

mat A = rand<mat>(10,10);
rowvec z = A.row(5); // extract a row vector

rowvec r(5);
r = 123.0;

rowvec r(5);
r.fill(123.0);

rowvec r = 123.0 * ones<rowvec>(1,5);

cube x(1,2,3);
cube y = rand<cube>(4,5,6);

mat A = y.slice(1);  // extract a slice from the cube
                     // (each slice is a matrix)

mat B = rand<mat>(4,5);
y.slice(2) = B;     // set a slice in the cube

cube q = y + y;     // cube addition
cube r = y % y;     // element-wise cube multiplication

mat A = rand<mat>(5,6);
mat B;
B.copy_size(A);

cout << B.n_rows << endl;
cout << B.n_cols << endl;

mat    X = rand<mat>(5,5);
colvec a = X.diag();
colvec b = X.diag(1);
colvec c = X.diag(-2);
X.diag() = rand<colvec>(5);

mat A = rand<mat>(10,10);
A(9,9) = 123.0;
double x = A.at(9,9);
double y = A[99];

colvec p = rand<colvec>(10,1);
p(9) = 123.0;
double z = p[9];

mat A(5,5);
A.fill(123.0);

mat A = rand<mat>(5,5);
mat B = rand<mat>(5,5);

B(1,1) = 0.0 / 0.0;  // create a NaN

cout << A.is_finite() << endl;
cout << B.is_finite() << endl;

colvec q(5);
cout << "q has " << q.n_cols << " columns" << endl;

mat A = rand<mat>(5,10);
A.ones();      // sets all elements to one
A.ones(10,20); // sets the size to 10 columns and 20 rows
               // followed by setting all elements to one

mat A = rand<mat>(5,10);
mat B = rand<mat>(5,10);
mat C = rand<mat>(10,5);

mat P = A + B;
mat Q = A - B;
mat R = -B;
mat S = A / 123.0;
mat T = A % B;
mat U = A * C;

// V is constructed without temporaries
mat V = A + B + A + B;

imat AA = "1 2 3; 4 5 6; 7 8 9;";
imat BB = "3 2 1; 6 5 4; 9 8 7;";

// compare elements
umat ZZ = (AA >= BB);

mat A = rand<mat>(5,5);
mat B = rand<mat>(6,6);

A.print();

// "B =" is the optional header line
B.print("B =");

cout << A << endl;
cout << "B =" << endl << B << endl;

mat A = rand(5,5);

cout.precision(11);
cout.setf(ios::fixed);

A.raw_print(cout, "A =");

mat A = rand<mat>(5, 5);
A.reset();

mat A = rand<mat>(5,5);

// default save format is arma_binary
A.save("A1.mat");
A.save("A2.mat", arma_ascii);

mat B;
// automatically detect format type
B.load("A1.mat");

mat C;
// force loading in the arma_ascii format
C.load("A2.mat", arma_ascii);

mat A;
A.set_size(5,10);

colvec q;
q.set_size(100);

cube A = rand<cube>(2,3,4);
mat B = A.slice(1);
A.slice(0) = rand<mat>(2,3);
A.subcube(0,0,1,  1,1,2) = rand<cube>(2,2,2);
A.slice(0)(1,2) = 99.0;

mat A = zeros<mat>(5,10);
mat B = A.submat(0,1,2,3);
A.submat(0,1,2,3) = rand<mat>(3,3);
A.row(1) = rand<mat>(1,10);

mat X = rand<mat>(5,5);
X.swap_rows(0,4);

mat A = rand<mat>(5,10);
A.zeros();      // sets all elements to zero
A.zeros(10,20); // sets the size to 10 columns and 20 rows
                // followed by setting all elements to zero

// create a colvec field with 3 rows and 2 columns
field<colvec> F(3,2);

// access components of the field
F(0,0) = colvec(5);
F(1,1) = rand<colvec>(6);
F(2,0).set_size(7);

// access element 1 of the colvec stored at 2,0
double x = F(2,0)(1);

// copy rows
F.row(0) = F.row(2);

// extract a row of colvecs from F
field<colvec> G = F.row(1);

// print the field to the standard output
G.print("G =");

// save the field to a binary file
G.save("colvec_field");


// create a field of strings
field<std::string> S(3,2);

S(0,0) = "hello";
S(1,0) = "there";

// string fields can be saved as plain text files
S.save("string_field");

running_stat<double> stats;

for(u32 i=0; i<10000; ++i)
  {
  double sample = double(rand())/RAND_MAX;
  stats(sample);
  }

cout << "mean = " << stats.mean() << endl;
cout << "var  = " << stats.var()  << endl;
cout << "min  = " << stats.min()  << endl;
cout << "max  = " << stats.max()  << endl;

running_stat_vec<rowvec::elem_type> stats;

rowvec sample;

for(u32 i=0; i<10000; ++i)
  {
  sample = rand<rowvec>(5);
  stats(sample);
  }

cout << "mean = " << stats.mean() << endl;
cout << "var  = " << stats.var()  << endl;
cout << "min  = " << stats.min()  << endl;
cout << "max  = " << stats.max()  << endl;

//
//

running_stat_vec<rowvec::elem_type> more_stats(true);

for(u32 i=0; i<20; ++i)
  {
  sample = rand<rowvec>(3);
  
  sample(1) -= sample(0);
  sample(2) += sample(1);
  
  more_stats(sample);
  }

cout << "covariance matrix = " << endl;
cout << more_stats.cov() << endl;

rowvec sd = more_stats.stddev();

cout << "correlations = " << endl;
cout << more_stats.cov() / (trans(sd) * sd);

wall_clock timer;

mat A = rand<mat>(4,4);
mat B = rand<mat>(4,4);
mat C;

timer.tic();
for(u32 i=0; i<100000; ++i)
  C = A + B + A + B;

double n_secs = timer.toc();
cout << "took " << n_secs << " seconds" << endl;

mat A = eye<mat>(5,5);
mat B = 123.0 * eye<mat>(5,5);

double start = 10.0;
double end   = 20.0;
u32    N     = 5;

mat    X     = linspace(start, end, N);
mat    Y     = linspace(start, end, N, 1);

colvec q     = linspace(start, end, N);
rowvec r     = linspace(start, end, N, 1);

mat A = ones<mat>(5,6);
mat B = 123.0 * ones<mat>(5,6);
cube C = 123.0 * ones<cube>(5,6,7);

colvec q = rand<colvec>(5);
rowvec r = rand<rowvec>(5);
mat A = rand<mat>(5,6);
cube C = rand<cube>(5,6,7);

mat A = rand<mat>(2, 3);
mat B = repmat(A, 4, 5);

colvec q = zeros<colvec>(5);
rowvec r = zeros<rowvec>(5);
mat A = zeros<mat>(5,6);
cube C = zeros<cube>(5,6,7);

mat A = rand<mat>(5,5);
mat B = abs(A); 

cx_mat X = rand<cx_mat>(5,5);
mat    Y = abs(X);

mat A = rand<mat>(4,5);

mat X = eps(A);

mat X = rand<mat>(5,5);
mat Y = exp(X);

mat X = rand<mat>(5,5);
mat Y = cos(X);

mat A = rand<mat>(5,5);
double x = accu(A);

// % performs element-wise multiplication, hence
// accu(A % B) is a "multiply-and-accumulate" operation
mat B = rand<mat>(5,5);
double y = accu(A % B);

rowvec r = rand<rowvec>(5);
colvec q = rand<colvec>(5);
mat    A = rand<mat>(5,5);

double x = as_scalar(r*q);
double y = as_scalar(r*A*q);

mat A = rand<mat>(5,5);
double x = det(A);

colvec a = rand<colvec>(10);
colvec b = rand<colvec>(10);

double x = dot(a,b);

colvec q = rand<colvec>(5);
double x = norm(q,2);

mat A = rand<mat>(4,5);

u32 r = rank(A);

mat A = rand<mat>(5,5);
double x = trace(A);

colvec q = rand<colvec>(10,1);
double x = max(q);

mat A = rand<mat>(10,10);
rowvec b = max(A);

// same result as max(A)
// the 0 explicitly indicates "traverse across rows"
rowvec c = max(A,0); 

// the 1 explicitly indicates "traverse across columns"
colvec d = max(A,1);

// find the overall maximum value
double y = max(max(A));

colvec q = rand<colvec>(10,1);
double x = prod(q);

mat A = rand<mat>(10,10);
rowvec b = prod(A);

// same result as prod(A)
// the 0 explicitly indicates "traverse across rows"
rowvec c = prod(A,0);

// 1 explicitly indicates "traverse across columns"
colvec d = prod(A,1);

// find the overall product
double y = prod(prod(A));

colvec q = rand<colvec>(10,1);
double x = sum(q);

mat A = rand<mat>(10,10);
rowvec b = sum(A);

// same result as sum(A)
// the 0 explicitly indicates "traverse across rows"
rowvec c = sum(A,0);

// 1 explicitly indicates "traverse across columns"
colvec d = sum(A,1);

// find the overall sum
double y = sum(sum(A));

mat A    = rand<mat>(5,5);
mat B    = mean(A);
mat C    = var(A);
double m = mean(mean(A));

colvec q = rand<colvec>(5);
double v = var(q);

mat  A = rand<mat>(5,5);

fmat B = conv_to<fmat>::from(A);

cx_mat X = rand<cx_mat>(5,5);
cx_mat Y = conj(X);

mat X = rand<mat>(4,5);
mat Y = rand<mat>(4,5);

mat R = cor(X,Y);

mat X = rand<mat>(4,5);
mat Y = rand<mat>(4,5);

mat C = cov(X,Y);

mat A = rand<mat>(5,5);
mat B = diagmat(A);
mat C = A*diagmat(A);

rowvec q = rand<rowvec>(5);
colvec r = rand<colvec>(5);
mat X = diagmat(q)*diagmat(r);

cx_mat A = rand<cx_mat>(5,10);
cx_mat B = htrans(A);

cx_mat A = rand<cx_mat>(5,5);
mat B = imag(A);
mat C = real(A);

mat A = rand<mat>(4,5);
mat B = rand<mat>(5,4);

mat K = kron(A,B);

mat A = rand<mat>(10, 5);
mat B = reshape(A, 5, 10);

mat A = rand<mat>(4,5);
mat B = shuffle(A);

mat A = rand<mat>(10,10);
mat B = sort(A);

colvec  q       = rand<colvec>(10);
ucolvec indices = sort_index(q);

mat A = rand<mat>(5,10);
mat B = trans(A);

mat X = rand<mat>(5,5);
mat Y = trans(X)*X;

mat R = chol(Y);

mat A = rand<mat>(10,10);
mat B = trans(A)*A;  // generate a symmetric matrix

colvec eigval;
mat    eigvec;

eig_sym(eigval, eigvec, B);

mat A = rand<mat>(10,10);

cx_colvec eigval;
cx_mat    eigvec;

eig_gen(eigval, eigvec, A);

mat A = rand<mat>(5,5);
mat B = inv(A);

// Diagonal elements in C are set to the reciprocal
// of the corresponding elements in A. Off-diagonal
// elements of C are set to zero.
mat C = inv(diagmat(A));

mat A = rand<mat>(10,10);

mat L;
mat U;
mat P;

lu(L, U, P, A);

mat B = trans(P)*L*U;

mat A = rand<mat>(4,5);

mat B = pinv(A);

mat X = rand<mat>(5,5);
mat Q, R;

qr(Q,R,X);

mat A    = rand<mat>(5,5);
colvec b = rand<colvec>(5);
mat B    = rand<mat>(5,5);

colvec x = solve(A,b);
mat X    = solve(A,B);

mat X = rand<mat>(5,5);

mat U;
colvec s;
mat V;
svd(U,s,V,X);

cout << "2.0 * pi = " << 2.0 * math::pi() << endl;
cout << "log_max, when using floats  = " << fmath::log_max() << endl;
cout << "log_max, when using doubles = " <<  math::log_max() << endl;

cout << "speed of light = " << phy::c_0() << endl;