Session 4#

Agenda#

  1. C++ Classes

  2. C++ Eigen

  3. C++ Armadillo

C++ Classes#

Simple Vector Class#

  1/*
  2 Simple demonstration of class to implement mathematical vectors
  3
  4 Useful References:
  5 * Discovering Modern C++ by Peter Gottschling (amongst many other C++ books)
  6 * https://en.cppreference.com/
  7*/
  8
  9// This is a header file, so we want to enure if only gets
 10// included once in the code so we don't end up with multiple
 11// redefinitions of the same functions in a project.
 12// On way to do this is to define a uniquely named variable
 13// using the prepocessor. Here, we name a variable "_SIMPLE_VECTOR_H_"
 14// If that variable is already defined, then nothing happens.
 15// Otherwise, we define it and the code defining the class is
 16// used.
 17#ifndef _SIMPLE_VECTOR_H_
 18#define _SIMPLE_VECTOR_H_
 19
 20// for C++ vectors ... these are arrays not math vectors
 21#include <vector>
 22
 23// Next we'll define a namespace. We don't technically have to
 24// do this, but we might want to use other libraries that have
 25// the same function or object names, so this will let us
 26// use both
 27namespace smu_demo_namespace
 28{
 29
 30    // We'll make a class named Vector
 31    class Vector
 32    {
 33
 34        // I like to define all of my member variables
 35        // at the top of a class. This is a style
 36        // preference
 37        //
 38        // in a class functions, variables, and objects
 39        // can be private, public, or protected
 40        //
 41        // private means only other members of this class
 42        // can modify or call
 43        //
 44        // protected means that other classes that this is
 45        // a sub class (also commonly called child) can modify
 46        // or call
 47        //
 48        // public means anything can modify or call.
 49        //
 50        // Generally speaking, all variables should either
 51        // be private or protected. The idea is that if you
 52        // need to access data in this class it should be
 53        // through a function and not directly.
 54
 55    private:
 56        // an array to contain the elements of our vector
 57        // its under the private heading because nothing outside
 58        // our class should be able to access it directly
 59        std::vector<double> m_data;
 60
 61        // Next we can declare the functions. For this, almost
 62        // everything will be public because we want to be
 63        // able to call them, but it is common to have private
 64        // or protected helper functions that you do not want
 65        // to be calld directly
 66
 67        // We'll have a private function to check sizes
 68
 69    public:
 70        // Constructors sets default arguments. The compiler automatically
 71        // makes one that doesn't do anything, but we'll want to initialize
 72        // our data with a constructor. We'll default to setting the vector
 73        // to 0 with the specified number of elements.
 74        // The name of the constructor always matches the name of the class
 75        Vector(const std::size_t &num_elements = 3);
 76
 77        // There is always a destructor as well to free memory. We don't need
 78        // it here, but it's always good to define a virtual destructor
 79        // in case someone uses this class for something else
 80        // virtual means it can be redfined.
 81        // The name is always ~ClassName
 82        virtual ~Vector(){};
 83
 84        // now we can define some vector functions, we might want to
 85        //
 86        // 1. add 2 vectors
 87        // 2. compute a dot product
 88        // 3. compute the norm
 89        // 4. access elements so we can modify values
 90        // 6. reset the vector a single value (e.g 0 or 1)
 91        // 7. normalize the vector
 92        // etc.
 93        // We won't implement every possible thing we might want,
 94        // but just a few to demonstrate the idea
 95
 96        // First lets define some functions to add vectors
 97        void add(const Vector &b);      
 98        void operator+(const Vector &b); 
 99
100        // for case 2 to work, we also need to define
101        // an = operator
102        void operator=(const Vector &b);
103
104        // compute a dot product
105        // we shouldn't modify the data so we add const
106        double dot(const Vector &b) const;
107
108        // compute the norm
109        double norm() const;
110
111        // functions to access elements / change values
112        // the first will let us change the value, the
113        // second will not
114        // The next 2 are the same ... but using the []
115        // operator
116        double &get(const std::size_t &index);
117        double get(const std::size_t &index) const;
118        double &operator[](const std::size_t &index);
119        double operator[](const std::size_t &index) const;
120
121        // function to set all the values of the vector to
122        // something
123        void setAllValues(const double &val);
124
125        // function to normalize the vector
126        void normalize();
127
128        // function to print the vector
129        void print() const;
130
131        // function to get the number of elements
132        std::size_t numElements() const;
133
134    private:
135        // function to check sizes
136        void check_sizes(const Vector &b) const;
137
138    }; // notice this ends with ;
139
140} // smu_demo_namespace
141
142// End the preprocessor header guard if / else
143// I like to add a comment after it refering to the variable it set
144#endif //_SIMPLE_VECTOR_H_

  1/*
  2 Simple demonstration of class to implement mathematical vectors
  3
  4 This file defines (implements) the classes declared in vector_class.hpp
  5
  6
  7 Useful References:
  8 * Discovering Modern C++ by Peter Gottschling (amongst many other C++ books)
  9 * https://en.cppreference.com/
 10*/
 11
 12#include <cmath>     // for std::sqrt
 13#include <iostream>  // for IO
 14#include <exception> // for throwing exceptions
 15#include "vector_class.hpp"
 16
 17// enter our namespace
 18namespace smu_demo_namespace
 19{
 20
 21    // Start by defining our constructor
 22    // Since all these functions are inside the class
 23    // we need to address them with ClassName::foo
 24    // Notice we set a default value in the declaration in the
 25    // header and omit it here
 26    Vector::Vector(const std::size_t &num_elements)
 27    {
 28        // initialize our array to be all zeros
 29        m_data.assign(num_elements, 0.0);
 30    }
 31
 32    // define our adding operators
 33    void Vector::add(const Vector &b)
 34    {
 35
 36        // use our private helper function to
 37        // make sure the sizes match
 38        check_sizes(b);
 39
 40        for (std::size_t i = 0; i < numElements(); ++i)
 41        {
 42            // the .at() checks bounds unless you turn of
 43            // debug flags
 44            m_data.at(i) += b[i];
 45        }
 46    }
 47
 48    void Vector::operator+(const Vector &b)
 49    {
 50        add(b);
 51    }
 52
 53    void Vector::operator=(const Vector &b)
 54    {
 55        check_sizes(b);
 56        for (std::size_t i = 0; i < numElements(); ++i)
 57        {
 58            m_data.at(i) = b[i];
 59        }
 60    }
 61
 62    double Vector::dot(const Vector &b) const
 63    {
 64        check_sizes(b);
 65        double dot_prod = 0;
 66        for (std::size_t i = 0; i < numElements(); ++i)
 67        {
 68            // the .at() checks bounds unless you turn of
 69            // debug flags
 70            dot_prod += m_data.at(i) * b[i];
 71        }
 72
 73        return dot_prod;
 74    }
 75
 76    double Vector::norm() const
 77    {
 78        double norm = 0;
 79        for (std::size_t i = 0; i < numElements(); ++i)
 80        {
 81            // the .at() checks bounds unless you turn of
 82            // debug flags
 83            norm += m_data.at(i) * m_data.at(i);
 84        }
 85
 86        return std::sqrt(norm);
 87    }
 88
 89    double &Vector::get(const std::size_t &index)
 90    {
 91        // the .at will check bounds to we don't need
 92        // to do our own check
 93        return m_data.at(index);
 94    }
 95
 96    double Vector::get(const std::size_t &index) const
 97    {
 98        // the .at will check bounds to we don't need
 99        // to do our own check
100        return m_data.at(index);
101    }
102
103    double &Vector::operator[](const std::size_t &index)
104    {
105        return m_data.at(index);
106    }
107    double Vector::operator[](const std::size_t &index) const
108    {
109        return m_data.at(index);
110    }
111
112    void Vector::setAllValues(const double &val)
113    {
114        m_data.assign(numElements(), val);
115    }
116
117    void Vector::normalize()
118    {
119        double norm = this->norm();
120        if (norm > 0)
121        {
122            for (std::size_t i = 0; i < numElements(); ++i)
123            {
124                // the .at() checks bounds unless you turn of
125                // debug flags
126                m_data.at(i) /= norm;
127            }
128        }
129        else
130        {
131            // TODO: maybe this should throw an error
132            throw std::logic_error("Norm is not positive");
133        }
134    }
135
136    void Vector::print() const
137    {
138        for (std::size_t i = 0; i < numElements() - 1; ++i)
139        {
140            std::cout << m_data.at(i) << ", ";
141        }
142        std::cout << m_data.back() << std::endl;
143    }
144
145    std::size_t Vector::numElements() const
146    {
147        return m_data.size();
148    }
149
150    void Vector::check_sizes(const Vector &b) const
151    {
152        if (numElements() != b.numElements())
153        {
154            throw std::range_error("The vectors are not the same length!");
155        }
156    }
157
158} // smu_demo_namespace

 1/*
 2 Simple demonstration of our vector class
 3
 4 This file runs the classes declared in
 5 vector_class.hpp and implemented in vector_class.cpp
 6
 7 Useful References:
 8 * Discovering Modern C++ by Peter Gottschling (amongst many other C++ books)
 9 * https://en.cppreference.com/
10*/
11
12#include <iostream>  // for IO
13#include <exception> // for throwing exceptions
14#include "vector_class.hpp"
15
16int main()
17{
18
19    // use our namespace
20    using namespace smu_demo_namespace;
21
22    // create 3 vectors
23    Vector a, b, c;
24
25    // set a to (1, 2, 3)
26    a[0] = 1;
27    a.get(1) = 2;
28    a[2] = 3;
29
30    std::cout << "a: " << std::endl;
31    a.print();
32
33    // set b to (4, 5, 6)
34    b.get(0) = 4;
35    b[1] = 5;
36    b.get(2) = 6;
37    std::cout << "b: " << std::endl;
38    b.print();
39
40    // print c
41    std::cout << "c: " << std::endl;
42    c.print();
43
44    // add b to a
45    a + b;
46    std::cout << "a: " << std::endl;
47    a.print();
48
49    c = a;
50    std::cout << "c: " << std::endl;
51    c.print();
52
53    for (int i = 0; i < c.numElements(); ++i)
54    {
55        c[i] = 1;
56    }
57    std::cout << "c: " << std::endl;
58    c.print();
59
60    std::cout << "c norm: " << c.norm() << std::endl;
61    std::cout << "normalized c:" << std::endl;
62    c.normalize();
63    c.print();
64
65    std::cout << "c dot c: " << c.dot(c) << std::endl;
66
67    Vector d(4);
68
69    // this should cause an error, but we can catch it
70    // by wrapping it in a try / catch block
71    try
72    {
73        c.dot(d);
74    }
75    catch (const std::exception &e)
76    {
77        // print out the error
78        // usually you stop the program after an error, but we'll keep going
79        std::cerr << e.what() << std::endl;
80    }
81
82    return 0;
83}

Class Inheritance#

  1/*
  2 Simple demonstration of classes and inheritance
  3
  4 Useful References:
  5 * Discovering Modern C++ by Peter Gottschling (amongst many other C++ books)
  6 * https://en.cppreference.com/
  7*/
  8
  9// This is a header file, so we want to enure if only gets
 10// included once in the code so we don't end up with multiple
 11// redefinitions of the same functions in a project.
 12// On way to do this is to define a uniquely named variable
 13// using the prepocessor. Here, we name a variable "_SIMPLE_INHERITANCE_H_"
 14// If that variable is already defined, then nothing happens.
 15// Otherwise, we define it and the code defining the class is
 16// used.
 17#ifndef _SIMPLE_INHERITANCE_H_
 18#define _SIMPLE_INHERITANCE_H_
 19
 20// Next we'll define a namespace. We don't technically have to
 21// do this, but we might want to use other libraries that have
 22// the same function or object names, so this will let us
 23// use both
 24namespace smu_demo_namespace
 25{
 26
 27    // Now we'll declare our first class in our smu_demo_namespace
 28    // A class is an object that can container variables, objects,
 29    // and functions.
 30    //
 31    // This is a simple class to define some properties of 2D shapes
 32    class ShapeProperties
 33    {
 34        // in a class functions, variables, and objects
 35        // can be private, public, or protected
 36        //
 37        // private means only other members of this class
 38        // can modify or call
 39        //
 40        // protected means that other classes that this is
 41        // a sub class (also commonly called child) can modify
 42        // or call
 43        //
 44        // public means anything can modify or call.
 45
 46        // we'll say that out shape has a few properties,
 47        // namely perimeter and area. These will be protected
 48        // because we're going to want anohter class to be
 49        // able to change them.
 50    protected:
 51        // everything below this is protected until we get
 52        // to another designation
 53        double m_area;
 54        double m_perimeter;
 55
 56        // we might want to be able to print this values out,
 57        // so we'll make some functions. In order to call these
 58        // functions outside of the class, we need to make them
 59        // public
 60    public:
 61        // the const after the function tells the compiler
 62        // that these functions will not modify any of the
 63        // variables or objects in this class. They are
 64        // only returning a copy of the values.
 65        double getArea() const;
 66        double getPerimeter() const;
 67
 68    }; // notice the the class ends with ;
 69
 70    // Now we have a shape class, but it's not really usable.
 71    // What shape is it? A triangle, square, circle ...
 72    // We can use the shape class above as a base to define
 73    // specific shapes
 74    class Triangle : public ShapeProperties
 75    {
 76        // a triangle has 3 sides so we'll want
 77        // variables to hold those lengths
 78    private:
 79        double m_side1;
 80        double m_side2;
 81        double m_side3;
 82
 83    public:
 84        // classes create objects using constructors
 85        // The compiler will make a default constructor with no
 86        // arguments (that doesn't set anything). So we'll make
 87        // our own.
 88        // First a default constructor with no arguments ... when
 89        // we define this, we'll set all the sides to 1 by default
 90        // Second, we'll make constructor that takes 3 sides lenghts
 91        //
 92        // The constructor is always named the same as the class
 93        Triangle();
 94        Triangle(const double &s1, const double &s2, const double &s3);
 95
 96        // we might want some helper functions to compute area and
 97        // perimeter. These only need to be called by other members
 98        // of the class, so we'll make them private
 99    private:
100        void computeArea();
101        void computePerimeter();
102    };
103
104    class Circle : public ShapeProperties
105    {
106        // we just need a radius
107        // we'll make this const so it can't be changed
108    private:
109        const double m_radius;
110
111    public:
112        // classes create objects using constructors
113        // The compiler will make a default constructor with no
114        // arguments (that doesn't set anything). So we'll make
115        // our own.
116        // First a default constructor with no arguments ... when
117        // we define this, we'll set the radius to 1 by default
118        // Second, we'll make constructor that takes 1 radius arg
119        //
120        // The constructor is always named the same as the class
121        Circle();
122        Circle(const double &r);
123
124        // we might want some helper functions to compute area and
125        // perimeter. These only need to be called by other members
126        // of the class, so we'll make them private
127    private:
128        void computeArea();
129        void computePerimeter();
130    };
131
132    class Rectangle : public ShapeProperties
133    {
134        // We need the 2 side lengths
135    private:
136        double m_side1;
137        double m_side2;
138
139    public:
140        // classes create objects using constructors
141        // The compiler will make a default constructor with no
142        // arguments (that doesn't set anything). So we'll make
143        // our own.
144        // First a default constructor with no arguments ... when
145        // we define this, we'll set both sides to 1 by default
146        // Second, we'll make constructor that takes the 2 side
147        // lengths as arguments
148
149        // The constructor is always named the same as the class
150        Rectangle();
151        Rectangle(const double &s1, const double &s2);
152
153        // we might want some helper functions to compute area and
154        // perimeter. These only need to be called by other members
155        // of the class, so we'll make them private
156    private:
157        void computeArea();
158        void computePerimeter();
159    };
160
161    // A square is a special case of a rectangle
162    class Square : public Rectangle
163    {
164
165    public:
166        // classes create objects using constructors
167        // The compiler will make a default constructor with no
168        // arguments (that doesn't set anything). So we'll make
169        // our own.
170        // First a default constructor with no arguments ... when
171        // we define this, we'll set the sides to 1 by default
172        // Second, we'll make constructor that takes the a side
173        // length as arguments
174
175        // The constructor is always named the same as the class
176        Square();
177        Square(const double &s);
178    };
179
180} // namespace smu_demo_namespace
181
182// End the preprocessor header guard if / else
183// I like to add a comment after it refering to the variable it set
184#endif //_SIMPLE_INHERITANCE_H_

  1/*
  2 Simple demonstration of classes and inheritance
  3
  4 Useful References:
  5 * Discovering Modern C++ by Peter Gottschling (amongst many other C++ books)
  6 * https://en.cppreference.com/
  7*/
  8
  9// This is a header file, so we want to enure if only gets
 10// included once in the code so we don't end up with multiple
 11// redefinitions of the same functions in a project.
 12// On way to do this is to define a uniquely named variable
 13// using the prepocessor. Here, we name a variable "_SIMPLE_INHERITANCE_H_"
 14// If that variable is already defined, then nothing happens.
 15// Otherwise, we define it and the code defining the class is
 16// used.
 17#ifndef _SIMPLE_INHERITANCE_H_
 18#define _SIMPLE_INHERITANCE_H_
 19
 20// Next we'll define a namespace. We don't technically have to
 21// do this, but we might want to use other libraries that have
 22// the same function or object names, so this will let us
 23// use both
 24namespace smu_demo_namespace
 25{
 26
 27    // Now we'll declare our first class in our smu_demo_namespace
 28    // A class is an object that can container variables, objects,
 29    // and functions.
 30    //
 31    // This is a simple class to define some properties of 2D shapes
 32    class ShapeProperties
 33    {
 34        // in a class functions, variables, and objects
 35        // can be private, public, or protected
 36        //
 37        // private means only other members of this class
 38        // can modify or call
 39        //
 40        // protected means that other classes that this is
 41        // a sub class (also commonly called child) can modify
 42        // or call
 43        //
 44        // public means anything can modify or call.
 45
 46        // we'll say that out shape has a few properties,
 47        // namely perimeter and area. These will be protected
 48        // because we're going to want anohter class to be
 49        // able to change them.
 50    protected:
 51        // everything below this is protected until we get
 52        // to another designation
 53        double m_area;
 54        double m_perimeter;
 55
 56        // we might want to be able to print this values out,
 57        // so we'll make some functions. In order to call these
 58        // functions outside of the class, we need to make them
 59        // public
 60    public:
 61        // the const after the function tells the compiler
 62        // that these functions will not modify any of the
 63        // variables or objects in this class. They are
 64        // only returning a copy of the values.
 65        double getArea() const;
 66        double getPerimeter() const;
 67
 68    }; // notice the the class ends with ;
 69
 70    // Now we have a shape class, but it's not really usable.
 71    // What shape is it? A triangle, square, circle ...
 72    // We can use the shape class above as a base to define
 73    // specific shapes
 74    class Triangle : public ShapeProperties
 75    {
 76        // a triangle has 3 sides so we'll want
 77        // variables to hold those lengths
 78    private:
 79        double m_side1;
 80        double m_side2;
 81        double m_side3;
 82
 83    public:
 84        // classes create objects using constructors
 85        // The compiler will make a default constructor with no
 86        // arguments (that doesn't set anything). So we'll make
 87        // our own.
 88        // First a default constructor with no arguments ... when
 89        // we define this, we'll set all the sides to 1 by default
 90        // Second, we'll make constructor that takes 3 sides lenghts
 91        //
 92        // The constructor is always named the same as the class
 93        Triangle();
 94        Triangle(const double &s1, const double &s2, const double &s3);
 95
 96        // we might want some helper functions to compute area and
 97        // perimeter. These only need to be called by other members
 98        // of the class, so we'll make them private
 99    private:
100        void computeArea();
101        void computePerimeter();
102    };
103
104    class Circle : public ShapeProperties
105    {
106        // we just need a radius
107        // we'll make this const so it can't be changed
108    private:
109        const double m_radius;
110
111    public:
112        // classes create objects using constructors
113        // The compiler will make a default constructor with no
114        // arguments (that doesn't set anything). So we'll make
115        // our own.
116        // First a default constructor with no arguments ... when
117        // we define this, we'll set the radius to 1 by default
118        // Second, we'll make constructor that takes 1 radius arg
119        //
120        // The constructor is always named the same as the class
121        Circle();
122        Circle(const double &r);
123
124        // we might want some helper functions to compute area and
125        // perimeter. These only need to be called by other members
126        // of the class, so we'll make them private
127    private:
128        void computeArea();
129        void computePerimeter();
130    };
131
132    class Rectangle : public ShapeProperties
133    {
134        // We need the 2 side lengths
135    private:
136        double m_side1;
137        double m_side2;
138
139    public:
140        // classes create objects using constructors
141        // The compiler will make a default constructor with no
142        // arguments (that doesn't set anything). So we'll make
143        // our own.
144        // First a default constructor with no arguments ... when
145        // we define this, we'll set both sides to 1 by default
146        // Second, we'll make constructor that takes the 2 side
147        // lengths as arguments
148
149        // The constructor is always named the same as the class
150        Rectangle();
151        Rectangle(const double &s1, const double &s2);
152
153        // we might want some helper functions to compute area and
154        // perimeter. These only need to be called by other members
155        // of the class, so we'll make them private
156    private:
157        void computeArea();
158        void computePerimeter();
159    };
160
161    // A square is a special case of a rectangle
162    class Square : public Rectangle
163    {
164
165    public:
166        // classes create objects using constructors
167        // The compiler will make a default constructor with no
168        // arguments (that doesn't set anything). So we'll make
169        // our own.
170        // First a default constructor with no arguments ... when
171        // we define this, we'll set the sides to 1 by default
172        // Second, we'll make constructor that takes the a side
173        // length as arguments
174
175        // The constructor is always named the same as the class
176        Square();
177        Square(const double &s);
178    };
179
180} // namespace smu_demo_namespace
181
182// End the preprocessor header guard if / else
183// I like to add a comment after it refering to the variable it set
184#endif //_SIMPLE_INHERITANCE_H_

Eigen Demo#

  1/*
  2 Simple demonstration of the Eigen Library
  3
  4 Useful References:
  5 * https://eigen.tuxfamily.org/dox/
  6 * Discovering Modern C++ by Peter Gottschling (amongst many other C++ books)
  7 * https://en.cppreference.com/
  8*/
  9
 10#include <iostream>       // IO
 11#include <Eigen/Dense>    // The Eigen dense linear algebra library
 12#include <Eigen/Geometry> // for rotation matrices
 13#include <cmath>          // for trig functions
 14
 15int main()
 16{
 17
 18    // Eigen is a templated linear algebra library. Templating allows
 19    // the library to change the size and type of various objects at
 20    // compile time and it also allows the compiler to do some clever
 21    // optimizations. For instance, it can automatically expand
 22    // expressiosn like
 23    //
 24    // d = a + b + c
 25    //
 26    // where the variables are all vectors (or matrices) to
 27    //
 28    // for (int i=0; i<n; ++i) {
 29    //   d[i] = a[i] + b[i] + c[i]
 30    // }
 31    //
 32    // This makes many linear algebra operations much more efficient than
 33    // traditional blas / lapack routines that may require multiple
 34    // evaluations
 35
 36    // For linear algebra, Eigen is built around a Matrix class
 37    // see https://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html
 38    // This class has the form
 39    //
 40    // template<typename Scalar_, int Rows_, int Cols_, int Options_, int MaxRows_, int MaxCols_>
 41    // class Eigen::Matrix<Scalar_, Rows_, Cols_, Options_, MaxRows_, MaxCols_>
 42    //
 43    // Scalar is the datatype, e.g. a doulbe
 44    // Rows and Cols are the number of rows and cols, respectively
 45    // options lets you set some things like row/col major, alignment, etc.
 46    //
 47    // Eigen provides some convenient type defs so you don't need to remember all that, e.g.
 48    // Matrix2d is a 2x2 square matrix of doubles (Matrix<double, 2, 2>)
 49    // Vector4f is a vector of 4 floats(Matrix<float, 4, 1>)
 50    // RowVector3i is a row vector of 3 ints(Matrix<int, 1, 3>)
 51    // MatrixXf is a dynamic size matrix of floats(Matrix<float, Dynamic, Dynamic>)
 52    // VectorXf is a dynamic size vector of floats(Matrix<float, Dynamic, 1>)
 53    // Matrix2Xf is a partially fixed size(dynamic - size) matrix of floats(Matrix<float, 2, Dynamic>)
 54    // MatrixX3d is a partially dynamic size(fixed - size) matrix of double(Matrix<double, Dynamic, 3>)
 55    // you can also make your own typedefs, e.g
 56    // this is the same as Vector3d
 57    // typedef Eigen::Matrix<double, 3, 1> Vector;
 58
 59    // Define Vectors a,b,c of length 3
 60    Eigen::Vector3d a, b, c;
 61
 62    // set a to (1,2,3)
 63    a << 1, 2, 3;
 64
 65    // set b to (4, 5, 6)
 66    b(0) = 4;
 67    b(1) = 5;
 68    b(2) = 6;
 69
 70    // set c = a + b
 71    c = a + b;
 72
 73    // print out a (it's a column vector, so transpose it so it's on 1 line)
 74    std::cout << "a: " << a.transpose() << std::endl;
 75    std::cout << "b: " << b.transpose() << std::endl;
 76    std::cout << "c: " << c.transpose() << std::endl;
 77
 78    // set c to ones
 79    c.setOnes();
 80    std::cout << "c: " << c.transpose() << std::endl;
 81    c.normalize();
 82    std::cout << "c: " << c.transpose() << std::endl;
 83    c = Eigen::Vector3d::Random(3);
 84    std::cout << "c: " << c.transpose() << std::endl;
 85    c.normalize();
 86    std::cout << "c: " << c.transpose() << std::endl;
 87
 88    // create a 3x3 matrix
 89    Eigen::Matrix3d A, B, C;
 90    double theta = 1.23423; // random angle in radians
 91    double omega = 3.854;
 92    double gamma = -2.863237;
 93
 94    // make A & B be a rotation matrix
 95    // rotates in xy-plane
 96    A << std::cos(theta), -std::sin(theta), 0,
 97        std::sin(theta), std::cos(theta), 0,
 98        0, 0, 1;
 99
100    // rotates in xy-plane
101    B.setIdentity();
102    B(0, 0) = std::cos(omega);
103    B(0, 2) = -std::sin(omega);
104    B(2, 0) = std::sin(omega);
105    B(2, 2) = std::cos(omega);
106
107    // rotates around the vector c
108    C = Eigen::AngleAxisd(gamma, c);
109
110    // create a new rotation matrix that is the combination
111    // of A, B, and C
112    auto R = A * B * C;
113
114    std::cout << "R:" << std::endl
115              << R << std::endl;
116
117    // rotate some vectors
118    std::cout << "a: " << a.transpose() << " (norm = " << a.norm() << ")" << std::endl;
119    std::cout << "b: " << b.transpose() << " (norm = " << b.norm() << ")" << std::endl;
120    std::cout << "c: " << c.transpose() << " (norm = " << c.norm() << ")" << std::endl;
121
122    a = R * a;
123    b = R * b;
124    c = R * c;
125
126    std::cout << "after rotating:" << std::endl;
127    std::cout << "\ta: " << a.transpose() << " (norm = " << a.norm() << ")" << std::endl;
128    std::cout << "\tb: " << b.transpose() << " (norm = " << b.norm() << ")" << std::endl;
129    std::cout << "\tc: " << c.transpose() << " (norm = " << c.norm() << ")" << std::endl;
130
131    // the inverse of rotation matrix is it's transpose ...
132    a = R.transpose() * a;
133    b = R.transpose() * b;
134    c = R.transpose() * c;
135
136    std::cout << "after rotating back:" << std::endl;
137    std::cout << "\ta: " << a.transpose() << " (norm = " << a.norm() << ")" << std::endl;
138    std::cout << "\tb: " << b.transpose() << " (norm = " << b.norm() << ")" << std::endl;
139    std::cout << "\tc: " << c.transpose() << " (norm = " << c.norm() << ")" << std::endl;
140
141    // We might want dynamic sizes, so we can do (X is dynamic, where we had 3 above)
142    const int n = 1987;
143    Eigen::MatrixXd M = Eigen::MatrixXd::Random(n, n);
144    Eigen::VectorXd x = Eigen::VectorXd::Random(n);
145    Eigen::VectorXd f(n);
146
147    // set f = Mx
148    f = M * x;
149
150    // Eigen has solvers built in
151    // see https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html
152    // We'll use the CompleteOrthogonalDecomposition method here. It's
153    // very accurate and reasonably fast for small to medium matrices
154    Eigen::VectorXd y = M.completeOrthogonalDecomposition().solve(f);
155
156    // if everything worked, we expect y = x (up to floating point rounding
157    // errors anyway). So lets check that:
158    double abs_error = (y - x).norm();
159    double rel_error = abs_error / x.norm();
160
161    std::cout << "Solve gave absolute error: " << abs_error
162              << " and relative error: " << rel_error << std::endl;
163
164    return 0;
165}

Armadillo Demo#

  1#include <armadillo>
  2#include <fmt/core.h>
  3#include <iostream>
  4
  5int main(int argc, char** argv) {
  6  std::cout << "Armadillo version: " << arma::arma_version::as_string() << std::endl;
  7  
  8  // construct a matrix according to given size and form of element initialisation
  9  arma::mat A(2,3, arma::fill::zeros);
 10  
 11  // .n_rows and .n_cols are read only
 12  fmt::print("A.n_rows: {}\nA.n_cols: {}\n", A.n_rows, A.n_cols);
 13  
 14  A(1,2) = 456.0;  // access an element (indexing starts at 0)
 15  A.print("A:");
 16  
 17  A = 5.0;         // scalars are treated as a 1x1 matrix
 18  A.print("A:");
 19  
 20  A.set_size(4,5); // change the size (data is not preserved)
 21  
 22  A.fill(5.0);     // set all elements to a specific value
 23  A.print("A:");
 24  
 25  A = { { 0.165300, 0.454037, 0.995795, 0.124098, 0.047084 },
 26        { 0.688782, 0.036549, 0.552848, 0.937664, 0.866401 },
 27        { 0.348740, 0.479388, 0.506228, 0.145673, 0.491547 },
 28        { 0.148678, 0.682258, 0.571154, 0.874724, 0.444632 },
 29        { 0.245726, 0.595218, 0.409327, 0.367827, 0.385736 } };
 30        
 31  A.print("A:");
 32  
 33  // determinant
 34  fmt::print("det(A): {}\n", det(A));
 35  
 36  // inverse
 37  inv(A).print("inv(A):");
 38  
 39  // save matrix as a text file
 40  A.save("A.txt", arma::raw_ascii);
 41  
 42  // load from file
 43  arma::mat B;
 44  B.load("A.txt");
 45  
 46  // submatrices
 47  B(arma::span(0,2), arma::span(3,4)).print("B( span(0,2), span(3,4) ):");
 48  B(0, 3, arma::size(3,2)).print("B( 0,3, size(3,2) ):");
 49
 50  B.row(0).print("B.row(0):");
 51  B.col(1).print("B.col(1):");
 52 
 53  // transpose
 54  std::cout << "B.t(): " << std::endl << B.t() << std::endl;
 55  
 56  // maximum from each column (traverse along rows)
 57  std::cout << "max(B): " << std::endl << max(B) << std::endl;
 58  
 59  // maximum from each row (traverse along columns)
 60  std::cout << "max(B,1): " << std::endl << max(B,1) << std::endl;
 61  
 62  // maximum value in B
 63  std::cout << "max(max(B)) = " << max(max(B)) << std::endl;
 64  
 65  // sum of each column (traverse along rows)
 66  std::cout << "sum(B): " << std::endl << sum(B) << std::endl;
 67  
 68  // sum of each row (traverse along columns)
 69  std::cout << "sum(B,1) =" << std::endl << sum(B,1) << std::endl;
 70  
 71  // sum of all elements
 72  std::cout << "accu(B): " << arma::accu(B) << std::endl;
 73  
 74  // trace = sum along diagonal
 75  std::cout << "trace(B): " << arma::trace(B) << std::endl;
 76  
 77  // generate the identity matrix
 78  arma::mat C = arma::eye<arma::mat>(4,4);
 79  
 80  // random matrix with values uniformly distributed in the [0,1] interval
 81  arma::mat D = arma::randu<arma::mat>(4,4);
 82  D.print("D:");
 83  
 84  // row vectors are treated like a matrix with one row
 85  arma::rowvec r = { 0.59119, 0.77321, 0.60275, 0.35887, 0.51683 };
 86  r.print("r:");
 87  
 88  // column vectors are treated like a matrix with one column
 89  arma::vec q = { 0.14333, 0.59478, 0.14481, 0.58558, 0.60809 };
 90  q.print("q:");
 91  
 92  // convert matrix to vector; data in matrices is stored column-by-column
 93  arma::vec v = arma::vectorise(A);
 94  v.print("v:");
 95  
 96  // dot or inner product
 97  std::cout << "as_scalar(r*q): " << arma::as_scalar(r*q) << std::endl;
 98  
 99  // outer product
100  std::cout << "q*r: " << std::endl << q*r << std::endl;
101  
102  // multiply-and-accumulate operation (no temporary matrices are created)
103  std::cout << "accu(A % B) = " << arma::accu(A % B) << std::endl;
104  
105  // example of a compound operation
106  B += 2.0 * A.t();
107  B.print("B:");
108  
109  // imat specifies an integer matrix
110  arma::imat AA = { { 1, 2, 3 },
111                    { 4, 5, 6 },
112                    { 7, 8, 9 } };
113  
114  arma::imat BB = { { 3, 2, 1 }, 
115                    { 6, 5, 4 },
116                    { 9, 8, 7 } };
117  
118  // comparison of matrices (element-wise); output of a relational operator is a umat
119  arma::umat ZZ = (AA >= BB);
120  ZZ.print("ZZ:");
121  
122  // cubes ("3D matrices")
123  arma::cube Q( B.n_rows, B.n_cols, 2 );
124  
125  Q.slice(0) = B;
126  Q.slice(1) = 2.0 * B;
127  
128  Q.print("Q:");
129  
130  // 2D field of matrices; 3D fields are also supported
131  arma::field<arma::mat> F(4,3); 
132  
133  for(arma::uword col=0; col < F.n_cols; ++col) {
134    for(arma::uword row=0; row < F.n_rows; ++row) {
135      F(row,col) = arma::randu<arma::mat>(2,3);  // each element in field<mat> is a matrix
136    }
137  }
138  
139  F.print("F:");
140  
141  return 0;
142}