Value semantics

APM 50179 EP

The need for uvector 


class matrix
{
  public:

  // ....

  private:

    size_t m_nb_rows;
    size_t m_nb_cols;
    std::vector<double> m_data;
};

The need for uvector

std::vector

  • initializes all its values upon construction
  • even when these values will be modified just after
  • not efficient for intensive computation on large matrices

matrix matrix::sigmoid(const matrix& m)
{
  matrix res(m.nb_rows(), m.nb_cols());
  std::transform(m.m_data.begin(), m.m_data.end(), res.m_data.begin()
  [](double arg) { return 1. / (1. + std::exp(-arg)); });
}

Digression: heap allocation

Stack allocation


double m_data[N];
// N must be known during compilation
// N is not mutable, you cannot resize m_data;

Heap allocation


double* p_data = new double[size];
// size can be computed at runtime
delete[] p_data; // Never forget this
delete[] p_data; // Never do it twice on a non null pointer
p_data = nullptr;
delete[] p_data; // OK
>

Digression: heap allocation 


double* p1 = new double(size);
double* p2 = p1;
              

delete p1;
p2[2] // crash

Digression: heap allocation

The survival guide to dynamic allocation

  • Default initialize pointers to nullptr
  • "new" must have a "delete" counterpart
  • Always reset a deleted pointer to nullptr
  • Avoid dynamic allocation as much as possible

uvector - constructors 


class uvector
{
  public:

    uvector(std::size_t size = 0);
    uvector(std::size_t size, double value);

  private:

    double* p_data;
    std::size_t m_size;
};

uvector - constructors 


uvector::uvector(std::size_t size)
: p_data(nullptr), m_size(0)
{
  if(size != 0)
  {
    p_data = new double[size];
    m_size = size;
  }
}
              

uvector::uvector(std::size_t size, double value)
: uvector(size) // delegating constructor
{
  std::fill(p_data, p_data + size, value);
}

uvector - destructor 


class uvector
{
  public:

    uvector(std::size_t size = 0);
    uvector(std::size_t size, double value);

    ~uvector();
};
              

uvector::~uvector()
{
  delete[] p_data;
  p_data = nullptr;
  m_size = 0;
}

uvector - copy semantics 


class uvector
{
  public:

    uvector(std::size_t size = 0);
    uvector(std::size_t size, double value);

    ~uvector();

    uvector(const uvector&);
    uvector& operator=(const uvector&);
}

uvector - copy constructor 


uvector::uvector(const uvector& rhs)
: p_data(nullptr), m_size(0)
{
  if(rhs.m_size != 0)
  {
    p_data = new double[rhs.m_size];
    m_size = rhs.m_size;
    std::copy(rhs.p_data, rhs.p_data + m_size, p_data);
  }
}

uvector - assignment operator 


uvector& uvector::operator=(const uvector& rhs)
{
  delete[] p_data;
  p_data = new double[rhs.m_size];
  std::copy(rhs.p_data, rhs.p_data + rhs.m_size, p_data);
  m_size = rhs.m_size;
  return *this;
}
              

uvector v(2, 2.);
v = v;

uvector - assign operator

Copy and swap idiom


uvector& uvector::operator=(const uvector& rhs)
{
  double* tmp = new double[rhs.m_size];                // Always allocate a temporary
  std::copy(rhs.p_data, rhs.p_data + rhs.m_size, tmp); // copy
  std::swap(tmp, p_data);                              // then swap
  delete[] tmp;                                        // then delete the temporary
  m_size = rhs.m_size;
  return *this;
}
              

uvector& uvector::operator=(const uvector& rhs)
{
  double* tmp = new double[rhs.m_size];                // Always allocate a temporary
  std::copy(rhs.p_data, rhs.p_data + rhs.m_size, tmp); // copy
  tmp_size = rhs.m_size;                               // copy
  std::swap(tmp, p_data);                              // swap
  std::swap(tmp_size, m_size);                         // swap
  delete[] tmp;                                        // then delete the temporary
  return *this;
}

uvector - swap 


class uvector
{
  public:

    void swap(uvector& rhs);

  private:

    double* p_data;
    std::size_t m_size;
};

void swap(uvector& lhs, uvector& rhs);
              

void uvector::swap(uvector& rhs)
{
  using std::swap;
  swap(p_data, rhs.p_data);
  swap(m_size, rhs.m_size);
}

void swap(uvector& lhs, uvector& rhs)
{
  lhs.swap(rhs);
}

uvector - assign operator

Copy and swap idiom


uvector& uvector::operator=(const uvector& rhs)
{
  uvector tmp(rhs);
  swap(*this, tmp);
  return *this;
}

rvalue reference 


uvector compute(const uvector& param)
{
  uvector res;
  // do some computation ...
  return res;
}

// Inefficient copy
uvector res = compute(my_huge_param);
              

uvector& compute(const uvector& param)
{
  uvector res;
  // do some computation ...
  return res; // Binding a temporary to a non const reference!
}

// dangling reference
uvector& res = compute(my_huge_param);

rvalue reference 


uvector compute(const uvector& param)
{
  uvector res;
  // do some computation ...
  return res;
}

// rvalue reference
uvector&& res = compute(my_huge_param);

lvalue vs rvalue

lvalue

  • can be on both side of an assignment
  • can have a name (or not)
  • lvalue references bind to lvalues

rvalue

  • can be on right hand side of an assignment only
  • does not have a name
  • rvalue references bind to rvalues

lvalue vs rvalue 


uvector& func(uvector& t)
{
  // do some stuff
  return t;
}

uvector t, something;
func(t) = something;
  • short: If it has a name, it is an lvalue
  • long: http://en.cppreference.com/w/cpp/language/value_category

lvalue vs rvalue 


uvector compute(const huge_param& param);

void function(const uvector& param);
void function(uvector&& param);
              

uvector m;
function(m);
function(compute(m));
uvector&& ref = compute(m);
function(ref);
              

uvector m;
function(m);                        // calls function(const uvector& param);
function(compute(m));               // calls function(uvector&& param);
uvector&& ref = compute(m);
function(ref);                      // calls function(const uvector& param);
              

function(static_cast<uvector&&>(ref));
function(std::move(uvector));

uvector - move semantic 


class uvector
{
  public:

    uvector(std::size_t = 0);
    uvector(std::size_t size, double value);
    ~uvector();

    uvector(const uvector&)
    uvector& operator=(const uvector& rhs);

    uvector(uvector&&);
    uvector& operator=(uvector&& rhs);
};

uvector - move constructor 


uvector::uvector(uvector&& rhs)
: p_data(std::move(p_data)), m_size(std::move(rhs.m_size))
{
}
              

uvector::uvector(uvector&& rhs)
: p_data(std::move(p_data)), m_size(std::move(rhs.m_size))
{
  rhs.p_data = nullptr;
  rhs.m_size = 0;
}

Move assign operator 


uvector& uvector::operator=(uvector&& rhs)
{
  using std::swap;
  swap(p_data, rhs.p_data);
  swap(m_size, rhs.m_size);
  return *this;
}

Auto-generation rules

  • The default constructor is auto-generated if there is no:
    • user-declared constructor
  • The default destructor is auto-generated if there is no:
    • user-declared destructor

Auto-generation rules

  • The copy constructor is auto-generated if there is no:
    • user-declared move constructor
    • user-declared move assignment operator
  • The copy assignment operator is auto-generated if there is no:
    • user-declared move constructor
    • user-declared move assignment operator

Auto-generation rules

  • The move constructor is auto-generated if there is no:
    • user-declared copy constructor
    • user-declared copy assignment operator
    • user-declared destructor
  • The move assignment operator is auto-generated if there is no:
    • user-declared copy constructor
    • user-declared copy assignment operator
    • user-declared destructor

Auto-generation rules

If you need to declare and implement any of

  • destructor
  • copy constructor
  • copy assignment operators
  • move constructor
  • move assignment operator

... you probably need to implement all off them

Genetic algorithm

A population of candidate solutions (called individuals or phenotypes) to an optimization problem is evolved toward better solutions. Typical steps are:

  • Initialization
  • Selection
  • Genetic operators (crossover and mutations)
  • Optional heuristics
  • Termination

Travelling salesman problem

Given a list of cities and the distance between each pair of cities, find the shortest route that visits each city exactly once and returns to the original city.

Travelling salesman problem

  • A Chromosome is a city. It can compute its distance to another city.
  • A phenotype is a circuit holding the cities. The fitness of the phenotype is the inverse of the total distance of the circuit.
  • A population is a set of circuits.

Travelling salesman problem

Selection

  • Roulette wheel (probability to be selected proportional to fitness)
  • Rank (probability to be selected proportional to rank)
  • Tournament (best of small random groups
  • Elitism (ensures best phenotypes get into the enxt generation)

Travelling salesman problem

Write a select method that implements a tournament selection:

  • The method should accept a population and a tournament size arguments.
  • The tournament size should be a data member of the algorithm.

Travelling salesman problem

Write a mutate method:

  • A mutation is a permutation of two cities in the circuit.
  • The mutation method should accept a circuit and a mutation rate as arguments.
  • The mutation rate should be a data member of the algorithm.

Travelling salesman problem

Write a crossover method:

  • It copies the phenotype of a parent between two indices (two-points crossover)
  • It copies the other cities from the other parent (two-points crossover
  • It ensure it contains all the cities, and no duplicates (consistency)

Travelling salesman problem

Write the evolve method:

  • It accepts a population as an argument
  • It applies elitism for the fittest circuit
  • It creates a new population thanks to cross over
  • It applies mutation to the newly created population

Laplacian 1d

We want to solve the following equation $$ \left\{ \begin{array}{l} -u''(x) = 1 \; \text{for} \; x \in [0, 1] \\ u(0) = 0 \; \text{and} \; u(1) = 0. \end{array} \right. $$ The exact solution is $u(x) = \frac{x - x^2}{2}$.

Laplacian 1d

We use the approximation $$ u''(x) \approx \frac{u(x-h) - 2u(x) + u(x+h)}{h^2} $$

where $h = \frac{1}{N}$ and $N$ is a positive integer.

We replace $u$ by a discretization $u_i = u(i/N)$.

Laplacian 1d

Then the discretized problem becomes $$ \frac{2u_i-u_{i-1}-u_{i+1}}{h^2} = 1 \; \text{for} \; i = 1, \cdots , N-1, $$ and $$ u_0 = u_N = 0. $$

Laplacian 1d

The problem can be rewritten as $$ Ax = 1 $$ where $$ A = \frac{1}{h^2} \left( \begin{array}{ccccc} 2 & -1 & 0 & \cdots & 0 \\ -1 & 2 & -1 & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & -1 & 2 & -1 \\ 0 & \cdots & 0 & -1 & 2 \\ \end{array} \right) $$

Laplacian 1d

Using the class matrix introduced in the previous lecture, create the matrix $A$

Laplacian 1d

Create the vector $x$ where all the entries will be equal to $0$

Create the vector $b$ where all the entries will be equal to $1$

Laplacian 1d

There are many ways to solve this linear system: LU, GMRES, CG, ...

In the following, we will implement the CG (conjugate gradient) algorithm.

conjugate gradient algorithm

AXPY

Implement the following algorithm $$ y = \alpha x + y. $$

AYPX

Implement the following algorithm $$ y = x + \alpha y. $$
Implement the scalar product and a $L^2$-norm.
Implement the matrix-vector product.
Implement the conjugate gradient using the previous functions.

CSR (Compressed Sparse Row)

The matrix $A$ of the laplacian contains a lot of zero.

Then the matrix-vector product computes a lot of zero.

The CSR format stores only the non-zero elements.

CSR (Compressed Sparse Row)

This format is composed of three arrays
  • an array val containing the non-zero elements of the matrix
  • an array row_ptr of the size $n_{rows} + 1$ containing the cumlulative number of non-zero elements of the row $i$. It is defined by a recursive relation
    • row_ptr[0] = 0
    • row_ptr[i] = row_ptr[i - 1] + number of non-zero elements in the row $i-1$
  • an array col containing the column indices of the non-zero elements

CSR: an example

$$ A = \left( \begin{array}{cccc} 2 & 0 & 0 & 1 \\ 0 & -2 & 0 & 0 \\ 0 & 0 & 3 & 1 \\ -4 & 0 & 0 & 1 \\ \end{array} \right) $$
  • val = [2, 1, -2, 3, 1, -4, 1]
  • row_ptr = [0, 2, 3, 5, 7]
  • col = [0, 3, 1, 2, 3, 0, 3]

CSR

Create a class csr with the three arrays and the methods nb_rows and nb_cols as described in the previous class matrix.

CSR

Create a method find_element(i, j) which returns a pointer to the entry in val if the element exists, nullptr otherwise.

CSR

Create a method insert_element(i, j) which insert an element in (i, j) and return an iterator to this new entry in val.

CSR

Create the operator()(i, j) and initialize the 1d laplacian as done with the class matrix.

CSR

Implement the matrix-vector product and test the conjugate gradient.