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

mat B = A + A;
mat C = A * B;
mat D = A % B;

cx_mat X(A,B);

B.zeros();
B.set_size(10,10);
B.zeros(5,6);

//
// fixed size matrices:

mat::fixed<5,6> F;
F.ones();

mat44 G;
G.randn();

cout << mat22().randu() << endl;

//
// constructing matrices from
// auxiliary (external) memory:

double aux_mem[24];
mat H(aux_mem, 4, 6, false);

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

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

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

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

vec a(5);  a = 123.0;

vec a(5);  a.fill(123.0);

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

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

cube::fixed<4,5,6> f;
f.ones();

// 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 vec field with 3 rows and 2 columns
field<vec> F(3,2);

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

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

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

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

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

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

sp_mat A(5,6);
sp_mat B(6,5);

A(0,0) = 1;
A(1,0) = 2;

B(0,0) = 3;
B(0,1) = 4;

sp_mat C = 2*B;

sp_mat D = A*C;


// batch insertion of two values at (5, 6) and (9, 9)
umat locations;
locations << 5 << 9 << endr
          << 6 << 9 << endr;

vec values;
values << 1.5 << 3.2 << endr;

sp_mat X(locations, values);

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

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

double* mem = A.colptr(2);

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

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

X.diag()  = randu<vec>(5);
X.diag() += 6;

mat X = ones<mat>(6,5);
vec v = linspace<vec>(10,15,6);

// add v to each column in X
X.each_col() += v;

// subtract v from columns 0 through to 3 in X
X.cols(0,3).each_col() -= v;

uvec indices(2);
indices(0) = 2;
indices(1) = 4;

// copy v to columns 2 and 4 in X
X.each_col(indices) = v;

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

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

mat A;

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

mat B = { 1, 2, 3, 4, 5, 6 };  // C++11 only
B.reshape(2,3);

cx_mat A( randu<mat>(4,4), randu<mat>(4,4) );

real(A).eval().save("A_real.dat", raw_ascii);
imag(A).eval().save("A_imag.dat", raw_ascii);

mat A(5,5);
A.eye();

mat B;
B.eye(5,5);

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

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

mat X = A.i();
mat Y = (A+A).i();

mat B = randu<mat>(4,1);
mat Z = A.i() * B; // automatically converted to Z=solve(A,B)

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

cout << A.in_range(0,0) << endl;  // true
cout << A.in_range(3,4) << endl;  // true
cout << A.in_range(4,5) << endl;  // false

mat A = randu<mat>(5,5);
cout << A.is_empty() << endl;

A.reset();
cout << A.is_empty() << endl;

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

B(1,1) = datum::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;

// C++11 only example
// need to include <random>

std::mt19937 engine;  // Mersenne twister random number engine

std::uniform_real_distribution<double> distr(0.0, 1.0);
  
mat A(4,5);
  
A.imbue( [&]() { return distr(engine); } );

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

// at column 2, insert a copy of B;
// A will now have 12 columns
A.insert_cols(2, B);

// at column 1, insert 5 zeroed columns;
// B will now have 7 columns
B.insert_cols(1, 5);

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();

vec v = randu<vec>(10);

cout << "min value is " << v.min() << endl;


uword  index;
double min_val = v.min(index);

cout << "index of min value is " << index << endl;


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

uword  row;
uword  col;
double min_val2 = A.max(row,col);

cout << "max value is at " << row << ',' << col << endl;

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();

// print a transposed version of A
A.t().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(4,5);
A.randu();

mat B;
B.randu(6,7);

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

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

A.reshape(5,4);

mat A = randu<mat>(4,5);
A.resize(7,6);

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

cx_mat C = zeros<cx_mat>(4,5);

C.set_real(A);
C.set_imag(B);

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

cx_mat C = cx_mat(A,B);

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

vec q;
q.set_size(100);

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

A.shed_row(2);
A.shed_cols(2,4);

mat A = randu<mat>(5,5);
cout << A.size() << endl;

A.clear();
cout << A.empty() << endl;

mat A = zeros<mat>(5,10);

A.submat(0,1,2,3) = randu<mat>(3,3);

// the following three statements
// access the same part of A
mat B = A.submat(0,1,2,3);
mat C = A.submat( span(0,2), span(1,3) );
mat D = A( span(0,2), span(1,3) );

// the following two statements
// access the same part of A
A.col(1)        = randu<mat>(5,1);
A(span::all, 1) = randu<mat>(5,1);

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

// get all elements of X that are greater than 0.5
vec q = X.elem( find(X > 0.5) );

// set four specific elements of X to 1
uvec indices;
indices << 2 << 3 << 6 << 8;

X.elem(indices) = ones<vec>(4);

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

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

A.subcube(0,0,1,  1,1,2) = randu<cube>(2,2,2);
A( span(0,1), span(0,1), span(1,2) ) = randu<cube>(2,2,2);

mat A = zeros<mat>(4,5);
mat B =  ones<mat>(6,7);

A.swap(B);

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

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

// C++11 only example

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

// add 123 to every element
A.transform( [](double val) { return (val + 123.0); } );

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(uword 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(uword 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(uword 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() / (sd.t() * sd);

wall_clock timer;

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

timer.tic();
for(uword 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);

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

vec  v = ones<vec>(10);
uvec u = ones<uvec>(11);
mat  A = ones<mat>(5,6);
cube Q = ones<cube>(5,6,7);

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

vec  v = randu<vec>(5);
mat  A = randu<mat>(5,6);
cube Q = randu<cube>(5,6,7);

mat A = randu<mat>(2, 3);

mat B = repmat(A, 4, 5);

sp_mat A = speye<sp_mat>(5,5);

sp_mat A = sprandu<sp_mat>(100, 200, 0.1);

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

vec  v = zeros<vec>(10);
uvec u = zeros<uvec>(11);
mat  A = zeros<mat>(5,6);
cube Q = 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 B = eps(A);

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

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

mat44  B = randu<mat>(4,4);

double y = det(B);       // use fast algorithm by default
double z = det(B, true); // use slow algorithm

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

double x = dot(a,b);

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

double val;
double sign;

log_det(val, sign, A);

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

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

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

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

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

vec V = randu<vec>(10);
mat X = randu<mat>(5,5);


// status1 will be set to true if vector V has all non-zero elements
bool status1 = all(V);

// status2 will be set to true if vector V has all elements greater than 0.5
bool status2 = all(V > 0.5);

// status3 will be set to true if matrix X has all elements greater than 0.6;
// note the use of vectorise()
bool status3 = all(vectorise(X) > 0.6);

// generate a row vector indicating which columns of X have all elements greater than 0.7
umat A = all(X > 0.7);

vec V = randu<vec>(10);
mat X = randu<mat>(5,5);


// status1 will be set to true if vector V has any non-zero elements
bool status1 = any(V);

// status2 will be set to true if vector V has any elements greater than 0.5
bool status2 = any(V > 0.5);

// status3 will be set to true if matrix X has any elements greater than 0.6;
// note the use of vectorise()
bool status3 = any(vectorise(X) > 0.6);

// generate a row vector indicating which columns of X have elements greater than 0.7
umat A = any(X > 0.7);

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 = cumsum(A);

vec x = randu<vec>(10);
vec y = cumsum(x);

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

vec  v  = randn<vec>(1000); // Gaussian distribution

uvec h1 = hist(v, 11);
uvec h2 = hist(v, linspace<vec>(-2,2,11));

vec  v = randn<vec>(1000);  // Gaussian distribution

uvec h = histc(v, linspace<vec>(-2,2,11));

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

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

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 = resize(A, 7, 6);

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

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

vec  q       = randu<vec>(10);
uvec indices = sort_index(q);

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

mat B = symmatu(A);
mat C = symmatl(A);

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

mat B = trans(A);
mat C = A.t();    // equivalent to trans(A), but more compact

mat A = randu<mat>(5,5);
mat U = trimatu(A);
mat L = trimatl(A);

// tell the inv() function to look only
// at the upper triangular part
mat X = inv( trimatu(U) );

mat X;
X << 1 << 2 << endr
  << 2 << 3 << endr;

mat Y = unique(X);

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

vec v = vectorise(A);

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

mat R = chol(Y);

mat A = randu<mat>(50,50);
mat B = A.t()*A;  // generate a symmetric matrix

vec eigval;
mat eigvec;

eig_sym(eigval, eigvec, B);        // use standard algorithm by default

eig_sym(eigval, eigvec, B, "dc");  // use "divide & conquer" algorithm

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

cx_vec eigval;
cx_mat eigvec;

eig_gen(eigval, eigvec, A);

   vec X = randu<vec>(100);
cx_vec Y = fft(X, 128);

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


// tell inv() to look only at the upper triangular part of A
mat D = inv( trimatu(A) );


// tell inv() that AA is a symmetric positive definite matrix
mat AA = A*A.t();
mat E  = inv( sympd(AA) );


mat44 F = randu<mat>(4,4);

mat G = inv(F);       // use fast algorithm by default
mat H = inv(F, true); // use slow algorithm

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

mat L, U, P;

lu(L, U, P, A);

mat B = P.t()*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 X = randu<mat>(5,5);
mat Q, R;

qr(Q,R,X);

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

qr_econ(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;
solve(x2, A, b);

// tell solve() to look only at the upper triangular part of A
mat Y = solve( trimatu(A), B );


mat44 C = randu<mat>(4,4);
mat44 D = randu<mat>(4,4);

mat E = solve(C, D);       // use fast algorithm by default
mat F = solve(C, D, true); // use slow algorithm

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

mat U;
vec s;
mat V;

svd(U,s,V,X);        // use standard algorithm by default

svd(U,s,V,X, "dc");  // use "divide & conquer" algorithm

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

mat U;
vec s;
mat V;
svd_econ(U, s, V, X, 'l');

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

mat X1 = syl(A, B, C);

mat X2;
syl(X2, A, B, C);

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

B(1,1) = datum::nan;

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

cout << is_finite( 0.123456789 ) << endl;
cout << is_finite( datum::nan  ) << endl;
cout << is_finite( datum::inf  ) << endl;

// print "hello" to the current err1 stream
get_stream_err1() << "hello" << endl;

// change the err2 stream to be a file
ofstream f("my_log.txt");
set_stream_err2(f);

// trying to invert a singular matrix
// will print a message to the err2 stream
// and throw an exception
mat X = zeros<mat>(5,5);
mat Y = inv(X);

// disable messages being printed to the err2 stream
std::ostream nullstream(0);
set_stream_err2(nullstream);

cout << "2.0 * pi = " << 2.0 * datum::pi << endl;

cout << "speed of light = " << datum::c_0 << endl;

cout << "log_max for floats = ";
cout << fdatum::log_max << endl;

cout << "log_max for doubles = ";
cout << datum::log_max << endl;