mat A = randu<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 = randu<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 = randu<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 = randu<cube>(4,5,6);

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

mat B = randu<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

// 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");

// 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) = randu<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");

mat X(4,5);
cout << "X has " << X.n_cols << " columns" << endl;

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

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

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

colvec a = X.diag();
colvec b = X.diag(1);
colvec c = X.diag(-2);

X.diag() = randu<colvec>(5);

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

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

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

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

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

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

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

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

mat A = randu<mat>(1,5);
mat B = randu<mat>(5,1);
mat C = randu<mat>(5,5);

cout << A.is_vec() << endl;
cout << B.is_vec() << endl;
cout << C.is_vec() << endl;

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


mat::iterator a = X.begin();
mat::iterator b = X.end();

for(mat::iterator i=a; i!=b; ++i)
  {
  cout << *i << endl;
  }


mat::col_iterator c = X.begin_col(1);  // start of column 1
mat::col_iterator d = X.end_col(3);    // end of column 3

for(mat::col_iterator i=c; i!=d; ++i)
  {
  cout << *i << endl;
  (*i) = 123.0;
  }

cube X = randu<cube>(2,3,4);


cube::iterator a = X.begin();
cube::iterator b = X.end();

for(cube::iterator i=a; i!=b; ++i)
  {
  cout << *i << endl;
  }


cube::slice_iterator c = X.begin_slice(1);  // start of slice 1
cube::slice_iterator d = X.end_slice(2);    // end of slice 2

for(cube::slice_iterator i=c; i!=d; ++i)
  {
  cout << *i << endl;
  (*i) = 123.0;
  }

      mat A = randu<mat>(5,5);
const mat B = randu<mat>(5,5);

      double* A_mem = A.memptr();
const double* B_mem = B.memptr();

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

mat A = randu<mat>(5,10);
mat B = randu<mat>(5,10);
mat C = randu<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 = randu<mat>(5,5);
mat B = randu<mat>(6,6);

A.print();

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

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

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

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

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

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

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

A.save("A1.mat");  // default save format is arma_binary
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);


// example of saving/loading using a stream
std::stringstream s;
A.save(s);

mat D;
D.load(s);


// example of testing for success
mat E;
bool status = E.load("A2.mat");

if(status == true)
  {
  cout << "loaded okay" << endl;
  }
else
  {
  cout << "problem with loading" << endl;
  }

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

colvec q;
q.set_size(100);

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

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

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

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

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 = randu<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 = randu<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 = randu<mat>(4,4);
mat B = randu<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);

colvec q = linspace<colvec>(10, 20, 5);
rowvec r = linspace<rowvec>(10, 20, 5);

mat X = linspace<mat>(10, 20, 5);

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 = randu<colvec>(5);
rowvec r = randu<rowvec>(5);
mat    A = randu<mat>(5,6);
cube   C = randu<cube>(5,6,7);

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

vec A = randu<vec>(5);
mat X = toeplitz(A);

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 = randu<mat>(5,5);
mat B = abs(A); 

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

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

mat X = eps(A);

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

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

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

mat B = randu<mat>(5,5);
double y = accu(A % B);

// operator % performs element-wise multiplication,
// hence accu(A % B) is a "multiply-and-accumulate"
// operation

rowvec r = randu<rowvec>(5);
colvec q = randu<colvec>(5);
mat    X = randu<mat>(5,5);

// examples of some expressions
// for which optimised implementations exist

double a = as_scalar(r*q);
double b = as_scalar(r*X*q);
double c = as_scalar(r*diagmat(X)*q);
double d = as_scalar(r*inv(diagmat(X))*q);

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

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

double x = dot(a,b);

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

double val;
double sign;

log_det(val, sign, A);

colvec q = randu<colvec>(5);
double x = norm(q, 2);
double y = norm(q, "inf");

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

u32 r = rank(A);

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

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

colvec d = diagvec(A);

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

mat A = randu<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 = randu<colvec>(10,1);
double x = prod(q);

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

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

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

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

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

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

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

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

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

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

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

vec A = randu<vec>(128) - 0.5;
vec B = randu<vec>(128) - 0.5;

vec C = conv(A,B);

mat  A = randu<mat>(5,5);
fmat B = conv_to<fmat>::from(A);

typedef std::vector<double> stdvec;

stdvec x(3);
x[0] = 0.0; x[1] = 1.0;  x[2] = 2.0;

colvec y = conv_to< colvec >::from(x);
stdvec z = conv_to< stdvec >::from(y); 

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

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

mat R = cor(X,Y);

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

mat C = cov(X,Y);

vec a = randu<vec>(3);
vec b = randu<vec>(3);

vec c = cross(a,b);

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

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

mat  A  = randu<mat>(5,5);
mat  B  = randu<mat>(5,5);
uvec q1 = find(A > B);
uvec q2 = find(A > 0.5);
uvec q3 = find(A > 0.5, 3, "last");

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

mat B = fliplr(A);
mat C = flipud(A);

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

cx_mat A = randu<cx_mat>(5,5);

mat B = imag(A);
mat C = real(A);

mat A = randu<mat>(4,5);
mat B = randu<mat>(4,6);
mat C = randu<mat>(6,5);

mat X = join_rows(A,B);
mat Y = join_cols(A,C);

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

mat K = kron(A,B);

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

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

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

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

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

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

mat R = chol(Y);

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

vec eigval;
mat eigvec;

eig_sym(eigval, eigvec, B);

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

cx_vec eigval;
cx_mat eigvec;

eig_gen(eigval, eigvec, A);

mat A = randu<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 in C are set to zero.
mat C = inv(diagmat(A));

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

mat L;
mat U;
mat P;

lu(L, U, P, A);

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

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

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

mat coeff;
mat score;
vec latent;
vec tsquared;

princomp(coeff, score, latent, tsquared, A);

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

mat coeff;
vec latent;
vec explained; 
princomp_cov(coeff, latent, explained, A);

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

qr(Q,R,X);

mat A = randu<mat>(5,5);
vec b = randu<vec>(5);
mat B = randu<mat>(5,5);

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

vec x2;
bool status = solve(x2, A, b);

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

mat U;
vec s;
mat V;
bool status = svd(U,s,V,X);

cout << "2.0 * pi = " << 2.0 * math::pi() << endl;

cout << "log_max for floats = ";
cout << fmath::log_max() << endl;

cout << "log_max for doubles = ";
cout << math::log_max() << endl;

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