Introduction to C++

Hello world

%%file hello.cpp

#include <iostream>

int main() {
    std::cout << "Hello, world!" << std::endl;
}

Overwriting hello.cpp

Compilation

%%bash

g++ hello.cpp -o hello.exe

Execution

%%bash

./hello.exe

Hello, world!

Namespaces

Just like Python, C++ has namespaces that allow us to build large libraries without worrying about name collisions. In the Hello world program, we used the explicit name std::cout indicating that cout is a member of the standard workspace. We can also use the using keyword to import selected functions or classes from a namespace.

using std::cout;

int main()
{
    cout << "Hello, world!\n";
}

For small programs, we sometimes import the entire namespace for convenience, but this may cause namespace collisions in larger programs.

using namespace std;

int main()
{
    cout << "Hello, world!\n";
}

You can easily create your own namespace.

namespace sta_663 {
    const double pi=2.14159;

    void greet(string name) {
        cout << "\nTraditional first program\n";
        cout << "Hello, " << name << "\n";
    }
}

int main()
{
    cout << "\nUsing namespaces\n";
    string name = "Tom";
    cout << sta_663::pi << "\n";
    sta_663::greet(name);
}

Types

%%file dtypes.cpp

#include <iostream>
#include <complex>

using std::cout;

int main() {
    // Boolean
    bool a = true, b = false;

    cout << "and            " << (a and b) << "\n";
    cout << "&&             " << (a && b) << "\n";
    cout << "or             " << (a or b) << "\n";
    cout << "||             " << (a || b) << "\n";
    cout << "not            " << not (a or b) << "\n";
    cout << "!              " << !(a or b) << "\n";

    // Integral numbers
    cout << "char           " << sizeof(char) << "\n";
    cout << "short int      " << sizeof(short int) << "\n";
    cout << "int            " << sizeof(int) << "\n";
    cout << "long           " << sizeof(long) << "\n";

    // Floating point numbers
    cout << "float          " << sizeof(float) << "\n";
    cout << "double         " << sizeof(double) << "\n";
    cout << "long double    " << sizeof(long double) << "\n";
    cout << "complex double " << sizeof(std::complex<double>) << "\n";

    // Characters and strings
    char c = 'a'; // Note single quotes
    char word[] = "hello"; // C char arrays
    std::string s = "hello"; // C++ string

    cout << c << "\n";
    cout << word << "\n";
    cout << s << "\n";
}

Overwriting dtypes.cpp

%%bash

g++ dtypes.cpp -o dtypes.exe
./dtypes.exe

and            0
&&             0
or             1
||             1
not            0
!              0
char           1
short int      2
int            4
long           8
float          4
double         8
long double    16
complex double 16
a
hello
hello

Type conversions

%%file type.cpp
#include <iostream>
using std::cout;
using std::string;
using std::stoi;

int main() {
    char c = '3'; // A char is an integer type
    string s = "3"; // A string is not an integer type
    int i = 3;
    float f = 3;
    double d = 3;

    cout << c << "\n";
    cout << i << "\n";
    cout << f << "\n";
    cout << d << "\n";

    cout << "c + i is " << c + i << "\n";
    cout << "c + i is " << c - '0' + i << "\n";
    cout << "s + i is " << stoi(s) + i << "\n"; // Use std::stod to convert to double
}

Overwriting type.cpp

%%bash

g++ -o type.exe type.cpp -std=c++14

%%bash

./type.exe

3
3
3
3
c + i is 54
c + i is 6
s + i is 6

Command line inputs

%%file command_args.cpp
#include <iostream>
#include <exception>
using std::cout;
using std::stoi;

int main(int argc, char *argv[]) {
    for (int i=0; i<argc; i++) {
        cout << i << ": " << argv[i];
        try {
            stoi(argv[i]);
            cout << " is an integer\n";
        } catch (std::exception& e) {
            cout << " is not an integer\n";
        }
    }
}

Overwriting command_args.cpp

%%bash

g++ -o command_args.exe command_args.cpp -std=c++14

%%bash

./command_args.exe  1 2 hello goodbye

0: ./command_args.exe is not an integer
1: 1 is an integer
2: 2 is an integer
3: hello is not an integer
4: goodbye is not an integer

Header, implementation and driver files

Header file(s)

%%file func02.hpp

double add(double x, double y);
double mult(double x, double y);

Overwriting func02.hpp

Implementation file(s)

%%file func02.cpp

double add(double x, double y) {
    return x + y;
}

double mult(double x, double y) {
    return x * y;
}

Overwriting func02.cpp

Driver program

%%file test_func02.cpp

#include <iostream>
#include "func02.hpp"

int main() {
    double a = 3;
    double b = 4;

    std::cout << add(a, b) << std::endl;
    std::cout << mult(a, b) << std::endl;

}

Overwriting test_func02.cpp

Compilation

%%bash

g++ test_func02.cpp func02.cpp -o test_func02.exe

Execution

%%bash

./test_func02.exe

7
12

Using make

%%file Makefile

test_func02.exe: test_func02.o func02.o
     g++ -o test_func02.exe test_func02.o func02.o

test_func02.o: test_func02.cpp func02.hpp
     g++ -c test_func02.cpp

func02.o: func02.cpp
     g++ -c func02.cpp

Overwriting Makefile

Compilation

%%bash

make

g++ -c test_func02.cpp
g++ -c func02.cpp
g++ -o test_func02.exe test_func02.o func02.o

Execution

%%bash

./test_func02.exe

7
12

A more flexible Makefile

%%file Makefile2

CC=g++
CFLAGS=-Wall -std=c++14

test_func02.exe: test_func02.o func02.o
     $(CC) $(CFLAGS) -o test_func02.exe test_func02.o func02.o

test_func02.o: test_func02.cpp func02.hpp
     $(CC) $(CFLAGS) -c test_func02.cpp

func02.o: func02.cpp
     $(CC) $(CFLAGS) -c func02.cpp

Overwriting Makefile2

Compilation

Note that no re-compilation occurs!

%%bash

make -f Makefile2

make: `test_func02.exe' is up to date.

Execution

%%bash

./test_func02.exe

7
12

Input and output

%%file data.txt
9 6

Overwriting data.txt

%%file io.cpp

#include <fstream>
#include "func02.hpp"

int main() {
    std::ifstream fin("data.txt");
    std::ofstream fout("result.txt");

    double a, b;

    fin >> a >> b;
    fin.close();

    fout << add(a, b) << std::endl;
    fout << mult(a, b) << std::endl;
    fout.close();
}

Overwriting io.cpp

%%bash

g++ io.cpp -o io.exe func02.cpp

%%bash

./io.exe

! cat result.txt

15
54

Arrays

%%file array.cpp

#include <iostream>
using std::cout;
using std::endl;

int main() {

    int N = 3;
    double counts[N];

    counts[0] = 1;
    counts[1] = 3;
    counts[2] = 3;

    double avg = (counts[0] + counts[1] + counts[2])/3;

    cout << avg << endl;
}

Overwriting array.cpp

%%bash

g++ -o array.exe array.cpp

%%bash

./array.exe

2.33333

Loops

%%file loop.cpp

#include <iostream>
using std::cout;
using std::endl;
using std::begin;
using std::end;

int main()
{
    int x[] = {1, 2, 3, 4, 5};

    cout << "\nTraditional for loop\n";
    for (int i=0; i < sizeof(x)/sizeof(x[0]); i++) {
        cout << i << endl;
    }

    cout << "\nUsing iterators\n";
    for (auto it=begin(x); it != end(x); it++) {
        cout << *it << endl;
    }

    cout << "\nRanged for loop\n\n";
    for (auto &i : x) {
        cout << i << endl;
    }
}

Overwriting loop.cpp

%%bash

g++ -o loop.exe loop.cpp -std=c++14

%%bash

./loop.exe

Traditional for loop
0
1
2
3
4

Using iterators
1
2
3
4
5

Ranged for loop

1
2
3
4
5

Function arguments

%%file func_arg.cpp

#include <iostream>
using std::cout;
using std::endl;

// Value parameter
void f1(int x) {
    x *= 2;
    cout << "In f1    : x=" << x << endl;
}

// Reference parameter
void f2(int &x) {
    x *= 2;
    cout << "In f2    : x=" << x << endl;
}

/* Note
If you want to avoid side effects
but still use references to avoid a copy operation
use a const refernece like this to indicate that x cannot be changed

void f2(const int &x)
*/

/* Note
Raw pointers are prone to error and
generally avoided in modern C++
See unique_ptr and shared_ptr
*/

// Raw pointer parameter
void f3(int *x) {
    *x *= 2;
    cout << "In f3    : x=" << *x << endl;
}

int main() {
    int x = 1;

    cout << "Before f1: x=" << x << "\n";
    f1(x);
    cout << "After f1 : x=" << x << "\n";

    cout << "Before f2: x=" << x << "\n";
    f2(x);
    cout << "After f2 : x=" << x << "\n";

    cout << "Before f3: x=" << x << "\n";
    f3(&x);
    cout << "After f3 : x=" << x << "\n";
}

Overwriting func_arg.cpp

%%bash

c++ -o func_arg.exe func_arg.cpp --std=c++14

%%bash

./func_arg.exe

Before f1: x=1
In f1    : x=2
After f1 : x=1
Before f2: x=1
In f2    : x=2
After f2 : x=2
Before f3: x=2
In f3    : x=4
After f3 : x=4

Arrays, pointers and dynamic memory

We generally avoid using raw pointers in C++, but this is standard in C and you should at least understand what is going on.

In C++, we typically use smart pointers, STL containers or convenient array constructs provided by libraries such as Eigen and Armadillo.

Pointers and addresses

%%file p01.cpp

#include <iostream>

using std::cout;

int main() {
    int x = 23;
    int *xp;
    xp = &x;

    cout << "x                    " << x << "\n";
    cout << "Address of x         " << &x << "\n";
    cout << "Pointer to x         " << xp << "\n";
    cout << "Value at pointer to x " << *xp << "\n";
}

Overwriting p01.cpp

%%bash

g++ -o p01.exe p01.cpp -std=c++14
./p01.exe

x                    23
Address of x         0x7ffee96d66dc
Pointer to x         0x7ffee96d66dc
Value at pointer to x 23

Arrays

%%file p02.cpp

#include <iostream>

using std::cout;
using std::begin;
using std::end;

int main() {
    int xs[] = {1,2,3,4,5};

    int ys[3];
    for (int i=0; i<5; i++) {
        ys[i] = i*i;
    }

    for (auto x=begin(xs); x!=end(xs); x++) {
        cout << *x << " ";
    }
    cout << "\n";

    for (auto x=begin(ys); x!=end(ys); x++) {
        cout << *x << " ";
    }
    cout << "\n";
}

Overwriting p02.cpp

%%bash

g++ -o p02.exe p02.cpp -std=c++14
./p02.exe

9 16 3 4 5
0 1 4

Dynamic memory

%%file p03.cpp

#include <iostream>

using std::cout;
using std::begin;
using std::end;

int main() {

    // declare memory
    int *z = new int; // single integer
    *z = 23;

    // Allocate on heap
    int *zs = new int[3]; // array of 3 integers
    for (int i=0; i<3; i++) {
        zs[i] = 10*i;
    }

    cout << *z << "\n";

    for (int i=0; i < 3; i++) {
        cout << zs[i] << " ";
    }
    cout << "\n";

    // need for manual management of dynamically assigned memory
    delete z;
    delete[] zs;
}

Overwriting p03.cpp

%%bash

g++ -o p03.exe p03.cpp -std=c++14
./p03.exe

23
0 10 20

Pointer arithmetic

%%file p04.cpp

#include <iostream>

using std::cout;
using std::begin;
using std::end;

int main() {
    int xs[] = {100,200,300,400,500,600,700,800,900,1000};

    cout << xs << ": " << *xs  << "\n";
    cout << &xs << ": " << *xs  << "\n";
    cout << &xs[3] << ": " << xs[3] << "\n";
    cout << xs+3 << ": " << *(xs+3)  << "\n";
}

Overwriting p04.cpp

%%bash

g++ -o p04.exe p04.cpp
./p04.exe

0x7ffee28046b0: 100
0x7ffee28046b0: 100
0x7ffee28046bc: 400
0x7ffee28046bc: 400

C style dynamic memory for jagged array (“matrix”)

%%file p05.cpp

#include <iostream>

using std::cout;
using std::begin;
using std::end;

int main() {
    int m = 3;
    int n = 4;
    int **xss = new int*[m]; // assign memory for m pointers to int
    for (int i=0; i<m; i++) {
        xss[i] = new int[n]; // assign memory for array of n ints
        for (int j=0; j<n; j++) {
            xss[i][j] = i*10 + j;
        }
    }

    for (int i=0; i<m; i++) {
        for (int j=0; j<n; j++) {
            cout << xss[i][j] << "\t";
        }
        cout << "\n";
    }

    // Free memory
    for (int i=0; i<m; i++) {
        delete[] xss[i];
    }
    delete[] xss;
}

Overwriting p05.cpp

%%bash

g++ -o p05.exe p05.cpp
./p05.exe

0       1       2       3
10      11      12      13
20      21      22      23

Functions

%%file func01.cpp

#include <iostream>

double add(double x, double y) {
    return x + y;
}

double mult(double x, double y) {
    return x * y;
}

int main() {
    double a = 3;
    double b = 4;

    std::cout << add(a, b) << std::endl;
    std::cout << mult(a, b) << std::endl;

}

Overwriting func01.cpp

%%bash

g++ -o func01.exe func01.cpp  -std=c++14
./func01.exe

7
12

Function parameters

%%file func02.cpp

#include <iostream>

double* add(double *x, double *y, int n) {
    double *res = new double[n];

    for (int i=0; i<n; i++) {
        res[i] = x[i] + y[i];
    }
    return res;
}

int main() {
    double a[] = {1,2,3};
    double b[] = {4,5,6};

    int n = 3;
    double *c = add(a, b, n);

    for (int i=0; i<n; i++) {
        std::cout << c[i] << " ";
    }
    std::cout << "\n";

    delete[] c; // Note difficulty of book-keeping when using raw pointers!
}

Overwriting func02.cpp

%%bash

g++ -o func02.exe func02.cpp  -std=c++14
./func02.exe

5 7 9

%%file func03.cpp

#include <iostream>
using std::cout;

// Using value
void foo1(int x) {
    x = x + 1;
}


// Using pointer
void foo2(int *x) {
    *x = *x + 1;
}

// Using ref
void foo3(int &x) {
    x = x + 1;
}

int main() {
    int x = 0;

    cout << x << "\n";
    foo1(x);
    cout << x << "\n";
    foo2(&x);
    cout << x << "\n";
    foo3(x);
    cout << x << "\n";
}

Overwriting func03.cpp

%%bash

g++ -o func03.exe func03.cpp  -std=c++14
./func03.exe

0
0
1
2

Generic programming with templates

%%file template.cpp

#include <iostream>

template<typename T>
T add(T a, T b) {
    return a + b;
}

int main() {
    int m =2, n =3;
    double u = 2.5, v = 4.5;

    std::cout << add(m, n) << std::endl;
    std::cout << add(u, v) << std::endl;
}

Overwriting template.cpp

%%bash

g++ -o template.exe template.cpp

%%bash

./template.exe

5
7

Anonymous functions

%%file lambda.cpp

#include <iostream>
using std::cout;
using std::endl;

int main() {

    int a = 3, b = 4;
    int c = 0;

    // Lambda function with no capture
    auto add1 = [] (int a, int b) { return a + b; };
    // Lambda function with value capture
    auto add2 = [c] (int a, int b) { return c * (a + b); };
    // Lambda funciton with reference capture
    auto add3 = [&c] (int a, int b) { return c * (a + b); };

    // Change value of c after function definition
    c += 5;

    cout << "Lambda function\n";
    cout << add1(a, b) <<  endl;
    cout << "Lambda function with value capture\n";
    cout << add2(a, b) <<  endl;
    cout << "Lambda function with reference capture\n";
    cout << add3(a, b) <<  endl;

}

Overwriting lambda.cpp

%%bash

c++ -o lambda.exe lambda.cpp --std=c++14

%%bash

./lambda.exe

Lambda function
7
Lambda function with value capture
0
Lambda function with reference capture
35

Function pointers

%%file func_pointer.cpp

#include <iostream>
#include <vector>
#include <functional>

using std::cout;
using std::endl;
using std::function;
using std::vector;

int main()
{
    cout << "\nUsing generalized function pointers\n";
    using func = function<double(double, double)>;

    auto f1 = [](double x, double y) { return x + y; };
    auto f2 = [](double x, double y) { return x * y; };
    auto f3 = [](double x, double y) { return x + y*y; };

    double x = 3, y = 4;

    vector<func> funcs = {f1, f2, f3,};

    for (auto& f : funcs) {
        cout << f(x, y) << "\n";
    }
}

Overwriting func_pointer.cpp

%%bash

g++ -o func_pointer.exe func_pointer.cpp -std=c++14

%%bash

./func_pointer.exe

Using generalized function pointers
7
12
19

Standard template library (STL)

%%file stl.cpp

#include <iostream>
#include <vector>
#include <map>
#include <unordered_map>

using std::vector;
using std::map;
using std::unordered_map;
using std::string;
using std::cout;
using std::endl;

struct Point{
    int x;
    int y;

    Point(int x_, int y_) :
      x(x_), y(y_) {};
};

int main() {
    vector<int> v1 = {1,2,3};
    v1.push_back(4);
    v1.push_back(5);

    cout << "Vecotr<int>" << endl;
    for (auto n: v1) {
        cout << n << endl;
    }
    cout << endl;

    vector<Point> v2;
    v2.push_back(Point(1, 2));
    v2.emplace_back(3,4);

    cout <<  "Vector<Point>" << endl;
    for (auto p: v2) {
        cout << "(" << p.x << ", " << p.y << ")" << endl;
    }
    cout << endl;

    map<string, int> v3 = {{"foo", 1}, {"bar", 2}};
    v3["hello"] = 3;
    v3.insert({"goodbye", 4});

    // Note the a C++ map is ordered
    // Note using (traditional) iterators instead of ranged for loop
    cout << "Map<string, int>" << endl;
    for (auto iter=v3.begin(); iter != v3.end(); iter++) {
        cout << iter->first << ": " << iter->second << endl;
    }
    cout << endl;

    unordered_map<string, int> v4 = {{"foo", 1}, {"bar", 2}};
    v4["hello"] = 3;
    v4.insert({"goodbye", 4});

    // Note the unordered_map is similar to Python' dict.'
    // Note using ranged for loop with const ref to avoid copying or mutation
    cout << "Unordered_map<string, int>" << endl;
    for (const auto& i: v4) {
        cout << i.first << ": " << i.second << endl;
    }
    cout << endl;
}

Overwriting stl.cpp

%%bash

g++ -o stl.exe stl.cpp -std=c++14

%%bash

./stl.exe

Vecotr<int>
1
2
3
4
5

Vector<Point>
(1, 2)
(3, 4)

Map<string, int>
bar: 2
foo: 1
goodbye: 4
hello: 3

Unordered_map<string, int>
goodbye: 4
hello: 3
bar: 2
foo: 1

STL algorithms

%%file stl_algorithm.cpp

#include <vector>
#include <iostream>
#include <numeric>

using std::cout;
using std::endl;
using std::vector;
using std::begin;
using std::end;

int main() {
    vector<int> v(10);

    // iota is somewhat like range
    std::iota(v.begin(), v.end(), 1);

    for (auto i: v) {
        cout << i << " ";
    }
    cout << endl;

    // C++ version of reduce
    cout << std::accumulate(begin(v), end(v), 0) << endl;

    // Accumulate with lambda
    cout << std::accumulate(begin(v), end(v), 1, [](int a, int b){return a * b; }) << endl;
}

Overwriting stl_algorithm.cpp

%%bash

g++ -o stl_algorithm.exe stl_algorithm.cpp -std=c++14

%%bash

./stl_algorithm.exe

1 2 3 4 5 6 7 8 9 10
55
3628800

Random numbers

%%file random.cpp

#include <iostream>
#include <random>
#include <functional>

using std::cout;
using std::random_device;
using std::mt19937;
using std::default_random_engine;
using std::uniform_int_distribution;
using std::poisson_distribution;
using std::student_t_distribution;
using std::bind;

// start random number engine with fixed seed
// Note default_random_engine may give differnet values on different platforms
// default_random_engine re(1234);

// or
// Using a named engine will work the same on differnt platforms
// mt19937 re(1234);

// start random number generator with random seed
random_device rd;
mt19937 re(rd());

uniform_int_distribution<int> uniform(1,6); // lower and upper bounds
poisson_distribution<int> poisson(30); // rate
student_t_distribution<double> t(10); // degrees of freedom

int main()
{
    cout << "\nGenerating random numbers\n";

    auto runif = bind (uniform, re);
    auto rpois = bind(poisson, re);
    auto rt = bind(t, re);

    for (int i=0; i<10; i++) {
        cout << runif() << ", " << rpois() <<  ", " << rt() << "\n";

    }
}

Overwriting random.cpp

%%bash

g++ -o random.exe random.cpp -std=c++14

%%bash

./random.exe

Generating random numbers
2, 28, -0.189118
6, 40, 0.271387
3, 27, 1.87313
6, 38, -0.110778
1, 27, 0.918611
6, 30, 1.07302
4, 34, 0.886518
2, 34, -0.550586
1, 23, 1.86504
4, 30, 0.178585

Statistics

A nicer library for working with probability distributions. Show integration with Armadillo. Integration with Eigen is also possible.

import os

if not os.path.exists('./stats'):
    ! git clone https://github.com/kthohr/stats.git

%%file stats.cpp
#define STATS_USE_ARMA
#include <iostream>
#include "stats.hpp"

using std::cout;
using std::endl;

// set seed for randome engine to 1776
std::mt19937_64 engine(1776);

int main() {
    // evaluate the normal PDF at x = 1, mu = 0, sigma = 1
    double dval_1 = stats::dnorm(1.0,0.0,1.0);

    // evaluate the normal PDF at x = 1, mu = 0, sigma = 1, and return the log value
    double dval_2 = stats::dnorm(1.0,0.0,1.0,true);

    // evaluate the normal CDF at x = 1, mu = 0, sigma = 1
    double pval = stats::pnorm(1.0,0.0,1.0);

    // evaluate the Laplacian quantile at p = 0.1, mu = 0, sigma = 1
    double qval = stats::qlaplace(0.1,0.0,1.0);

    // draw from a t-distribution dof = 30
    double rval = stats::rt(30);

    // matrix output
    arma::mat beta_rvs = stats::rbeta<arma::mat>(5,5,3.0,2.0);
    // matrix input
    arma::mat beta_cdf_vals = stats::pbeta(beta_rvs,3.0,2.0);

    cout << "evaluate the normal PDF at x = 1, mu = 0, sigma = 1" << endl;
    cout << dval_1 << endl;
    cout << "evaluate the normal PDF at x = 1, mu = 0, sigma = 1, and return the log value" << endl;
    cout << dval_2 << endl;
    cout << "evaluate the normal CDF at x = 1, mu = 0, sigma = 1" << endl;
    cout << pval << endl;
    cout << "evaluate the Laplacian quantile at p = 0.1, mu = 0, sigma = 1" << endl;
    cout << qval << endl;
    cout << "draw from a t-distribution dof = 30" << endl;
    cout << rval << endl;
    cout << "draws from a beta distribuiotn to populate Armadillo matrix" << endl;
    cout << beta_rvs << endl;
    cout << "evaluaate CDF for beta draws" << endl;
    cout << beta_cdf_vals << endl;
}

Overwriting stats.cpp

%%bash

g++ -std=c++11 -I./stats/include stats.cpp -o stats.exe

%%bash

./stats.exe

evaluate the normal PDF at x = 1, mu = 0, sigma = 1
0.241971
evaluate the normal PDF at x = 1, mu = 0, sigma = 1, and return the log value
-1.41894
evaluate the normal CDF at x = 1, mu = 0, sigma = 1
0.841345
evaluate the Laplacian quantile at p = 0.1, mu = 0, sigma = 1
-1.60944
draw from a t-distribution dof = 30
-1.80271
draws from a beta distribuiotn to populate Armadillo matrix
   0.7117   0.6939   0.3855   0.2690   0.3303
   0.5459   0.5328   0.6257   0.3173   0.5513
   0.8201   0.6763   0.8341   0.3584   0.7111
   0.8470   0.5821   0.6299   0.7157   0.5238
   0.6388   0.8524   0.8932   0.5235   0.6822

evaluaate CDF for beta draws
   0.6723   0.6410   0.1629   0.0622   0.1085
   0.3843   0.3632   0.5200   0.0974   0.3931
   0.8493   0.6098   0.8691   0.1346   0.6711
   0.8865   0.4445   0.5274   0.6793   0.3490
   0.5431   0.8935   0.9409   0.3485   0.6203

Numerics

Using Armadillo

%%file test_arma.cpp

#include <iostream>
#include <armadillo>

using std::cout;
using std::endl;

int main()
{
    using namespace arma;

    vec u = linspace<vec>(0,1,5);
    vec v = ones<vec>(5);
    mat A = randu<mat>(4,5); // uniform random deviates
    mat B = randn<mat>(4,5); // normal random deviates

    cout << "\nVecotrs in Armadillo\n";
    cout << u << endl;
    cout << v << endl;
    cout << u.t() * v << endl;

    cout << "\nRandom matrices in Armadillo\n";
    cout << A << endl;
    cout << B << endl;
    cout << A * B.t() << endl;
    cout << A * v << endl;

    cout << "\nQR in Armadillo\n";
    mat Q, R;
    qr(Q, R, A.t() * A);
    cout << Q << endl;
    cout << R << endl;
}

Overwriting test_arma.cpp

%%bash

g++ -o test_arma.exe test_arma.cpp -std=c++14 -larmadillo

%%bash

./test_arma.exe

Vecotrs in Armadillo
        0
   0.2500
   0.5000
   0.7500
   1.0000

   1.0000
   1.0000
   1.0000
   1.0000
   1.0000

   2.5000


Random matrices in Armadillo
   7.8264e-06   5.3277e-01   6.7930e-01   8.3097e-01   6.7115e-01
   1.3154e-01   2.1896e-01   9.3469e-01   3.4572e-02   7.6982e-03
   7.5561e-01   4.7045e-02   3.8350e-01   5.3462e-02   3.8342e-01
   4.5865e-01   6.7886e-01   5.1942e-01   5.2970e-01   6.6842e-02

  -0.7649   1.2041  -0.7020   1.1862   0.8284
   1.7313   0.0937   1.6814  -0.8631   0.9426
   0.1454  -0.6920   0.2742   0.6810   0.2091
   0.7032   0.2610   0.1752   1.3165  -0.5897

   1.7063   1.1075   0.5238   0.9563
  -0.4457   1.7973   0.1490   0.3544
  -0.4095   2.2726   0.2990   0.4551
   0.7857   1.3368   0.1140   1.2487

   2.7142
   1.3275
   1.6230
   2.2535


QR in Armadillo
  -0.6776   0.6377   0.3253   0.0944  -0.1395
  -0.3188  -0.4409   0.3521   0.4177   0.6368
  -0.5524  -0.2366  -0.7774   0.1504  -0.1092
  -0.2443  -0.5816   0.4071  -0.1709  -0.6380
  -0.2727  -0.0688  -0.0099  -0.8745   0.3950

  -1.1785e+00  -1.3394e+00  -2.1015e+00  -1.3527e+00  -1.0229e+00
            0  -8.3421e-01  -9.7607e-01  -9.9515e-01  -5.3245e-01
            0            0  -4.6336e-01   7.6793e-02  -4.0376e-03
            0            0            0  -2.0273e-01  -3.2745e-01
            0            0            0            0   1.1102e-16

Using Eigen

%%file test_eigen.cpp
#include <iostream>
#include <fstream>
#include <random>
#include <Eigen/Dense>
#include <functional>

using std::cout;
using std::endl;
using std::ofstream;

using std::default_random_engine;
using std::normal_distribution;
using std::bind;

// start random number engine with fixed seed
default_random_engine re{12345};

normal_distribution<double> norm(5,2); // mean and standard deviation
auto rnorm = bind(norm, re);

int main()
{
    using namespace Eigen;

    VectorXd x1(6);
    x1 << 1, 2, 3, 4, 5, 6;
    VectorXd x2 = VectorXd::LinSpaced(6, 1, 2);
    VectorXd x3 = VectorXd::Zero(6);
    VectorXd x4 = VectorXd::Ones(6);
    VectorXd x5 = VectorXd::Constant(6, 3);
    VectorXd x6 = VectorXd::Random(6);

    double data[] = {6,5,4,3,2,1};
    Map<VectorXd> x7(data, 6);

    VectorXd x8 = x6 + x7;

    MatrixXd A1(3,3);
    A1 << 1 ,2, 3,
          4, 5, 6,
          7, 8, 9;
    MatrixXd A2 = MatrixXd::Constant(3, 4, 1);
    MatrixXd A3 = MatrixXd::Identity(3, 3);

    Map<MatrixXd> A4(data, 3, 2);

    MatrixXd A5 = A4.transpose() * A4;
    MatrixXd A6 = x7 * x7.transpose();
    MatrixXd A7 = A4.array() * A4.array();
    MatrixXd A8 = A7.array().log();
    MatrixXd A9 = A8.unaryExpr([](double x) { return exp(x); });
    MatrixXd A10 = MatrixXd::Zero(3,4).unaryExpr([](double x) { return rnorm(); });

    VectorXd x9 = A1.colwise().norm();
    VectorXd x10 = A1.rowwise().sum();

    MatrixXd A11(x1.size(), 3);
    A11 << x1, x2, x3;

    MatrixXd A12(3, x1.size());
    A12 << x1.transpose(),
          x2.transpose(),
          x3.transpose();

    JacobiSVD<MatrixXd> svd(A10, ComputeThinU | ComputeThinV);


    cout << "x1: comman initializer\n" << x1.transpose() << "\n\n";
    cout << "x2: linspace\n" << x2.transpose() << "\n\n";
    cout << "x3: zeors\n" << x3.transpose() << "\n\n";
    cout << "x4: ones\n" << x4.transpose() << "\n\n";
    cout << "x5: constant\n" << x5.transpose() << "\n\n";
    cout << "x6: rand\n" << x6.transpose() << "\n\n";
    cout << "x7: mapping\n" << x7.transpose() << "\n\n";
    cout << "x8: element-wise addition\n" << x8.transpose() << "\n\n";

    cout << "max of A1\n";
    cout << A1.maxCoeff() << "\n\n";
    cout << "x9: norm of columns of A1\n" << x9.transpose() << "\n\n";
    cout << "x10: sum of rows of A1\n" << x10.transpose() << "\n\n";

    cout << "head\n";
    cout << x1.head(3).transpose() << "\n\n";
    cout << "tail\n";
    cout << x1.tail(3).transpose() << "\n\n";
    cout << "slice\n";
    cout << x1.segment(2, 3).transpose() << "\n\n";

    cout << "Reverse\n";
    cout << x1.reverse().transpose() << "\n\n";

    cout << "Indexing vector\n";
    cout << x1(0);
    cout << "\n\n";

    cout << "A1: comma initilizer\n";
    cout << A1 << "\n\n";
    cout << "A2: constant\n";
    cout << A2 << "\n\n";
    cout << "A3: eye\n";
    cout << A3 << "\n\n";
    cout << "A4: mapping\n";
    cout << A4 << "\n\n";
    cout << "A5: matrix multiplication\n";
    cout << A5 << "\n\n";
    cout << "A6: outer product\n";
    cout << A6 << "\n\n";
    cout << "A7: element-wise multiplication\n";
    cout << A7 << "\n\n";
    cout << "A8: ufunc log\n";
    cout << A8 << "\n\n";
    cout << "A9: custom ufucn\n";
    cout << A9 << "\n\n";
    cout << "A10: custom ufunc for normal deviates\n";
    cout << A10 << "\n\n";
    cout << "A11: np.c_\n";
    cout << A11 << "\n\n";
    cout << "A12: np.r_\n";
    cout << A12 << "\n\n";

    cout << "2x2 block startign at (0,1)\n";
    cout << A1.block(0,1,2,2) << "\n\n";
    cout << "top 2 rows of A1\n";
    cout << A1.topRows(2) << "\n\n";
    cout << "bottom 2 rows of A1";
    cout << A1.bottomRows(2) << "\n\n";
    cout << "leftmost 2 cols of A1";
    cout << A1.leftCols(2) << "\n\n";
    cout << "rightmost 2 cols of A1";
    cout << A1.rightCols(2) << "\n\n";

    cout << "Diagonal elements of A1\n";
    cout << A1.diagonal() << "\n\n";
    A1.diagonal() = A1.diagonal().array().square();
    cout << "Transforming diagonal eelemtns of A1\n";
    cout << A1 << "\n\n";

    cout << "Indexing matrix\n";
    cout << A1(0,0) << "\n\n";

    cout << "singular values\n";
    cout << svd.singularValues() << "\n\n";

    cout << "U\n";
    cout << svd.matrixU() << "\n\n";

    cout << "V\n";
    cout << svd.matrixV() << "\n\n";
}

Overwriting test_eigen.cpp

import os
if not os.path.exists('./eigen'):
    ! hg clone https://bitbucket.org/eigen/eigen/

%%bash

g++ -o test_eigen.exe test_eigen.cpp -std=c++14 -I./eigen

%%bash

./test_eigen.exe

x1: comman initializer
1 2 3 4 5 6

x2: linspace
  1 1.2 1.4 1.6 1.8   2

x3: zeors
0 0 0 0 0 0

x4: ones
1 1 1 1 1 1

x5: constant
3 3 3 3 3 3

x6: rand
 -0.999984  -0.736924   0.511211 -0.0826997  0.0655345  -0.562082

x7: mapping
6 5 4 3 2 1

x8: element-wise addition
 5.00002  4.26308  4.51121   2.9173  2.06553 0.437918

max of A1
9

x9: norm of columns of A1
8.12404 9.64365  11.225

x10: sum of rows of A1
 6 15 24

head
1 2 3

tail
4 5 6

slice
3 4 5

Reverse
6 5 4 3 2 1

Indexing vector
1

A1: comma initilizer
1 2 3
4 5 6
7 8 9

A2: constant
1 1 1 1
1 1 1 1
1 1 1 1

A3: eye
1 0 0
0 1 0
0 0 1

A4: mapping
6 3
5 2
4 1

A5: matrix multiplication
77 32
32 14

A6: outer product
36 30 24 18 12  6
30 25 20 15 10  5
24 20 16 12  8  4
18 15 12  9  6  3
12 10  8  6  4  2
 6  5  4  3  2  1

A7: element-wise multiplication
36  9
25  4
16  1

A8: ufunc log
3.58352 2.19722
3.21888 1.38629
2.77259       0

A9: custom ufucn
36  9
25  4
16  1

A10: custom ufunc for normal deviates
5.17237 3.73572 6.29422 6.55268
5.33713 3.88883 1.93637 4.39812
8.22086 6.94502 6.36617  6.5961

A11: np.c_
  1   1   0
  2 1.2   0
  3 1.4   0
  4 1.6   0
  5 1.8   0
  6   2   0

A12: np.r_
  1   2   3   4   5   6
  1 1.2 1.4 1.6 1.8   2
  0   0   0   0   0   0

2x2 block startign at (0,1)
2 3
5 6

top 2 rows of A1
1 2 3
4 5 6

bottom 2 rows of A14 5 6
7 8 9

leftmost 2 cols of A11 2
4 5
7 8

rightmost 2 cols of A12 3
5 6
8 9

Diagonal elements of A1
1
5
9

Transforming diagonal eelemtns of A1
 1  2  3
 4 25  6
 7  8 81

Indexing matrix
1

singular values
19.5001
2.80674
1.29869

U
  -0.5585   0.75124  0.351733
-0.406818 -0.617595  0.673109
-0.722895  -0.23284 -0.650544

V
-0.564245 -0.471949 0.0490777
-0.445586 -0.431952 -0.451578
-0.456672  0.730484 -0.480641
-0.523957  0.238904  0.750102

Check SVD

import numpy as np

A10 = np.array([
    [5.17237, 3.73572, 6.29422, 6.55268],
    [5.33713, 3.88883, 1.93637, 4.39812],
    [8.22086, 6.94502, 6.36617,  6.5961]
])

U, s, Vt = np.linalg.svd(A10, full_matrices=False)

s

array([19.50007376,  2.80674189,  1.29869186])

U

array([[-0.55849978, -0.75124103, -0.3517313 ],
       [-0.40681745,  0.61759344, -0.67311062],
       [-0.72289526,  0.2328417 ,  0.65054376]])

Vt.T

array([[-0.56424535,  0.47194895, -0.04907563],
       [-0.44558625,  0.43195279,  0.45157518],
       [-0.45667231, -0.73048295,  0.48064272],
       [-0.52395657, -0.23890509, -0.75010267]])