IPS-PROD-TP

class: titlepage

.title[Introduction to Scientific Programming - TP :notebook:]

.subtitle[`IPS-PROD-TP` - N. Dubray - ENSIIE - 2017]

.row[
`$$\int_{-\infty}^\infty f(x)e^{-x^2} dx \simeq \sum_{i=0}^{n-1}w_if(x_i)$$`
]

.footnote[
[:book:](../toc/index.html)
]

---

layout: true
class: animated fadeIn numbers

.footnote[
`IPS-PROD-TP` - N. Dubray - ENSIIE - 2017 - [:book:](../toc/index.html)
]

---

# Project rules

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
]

.column.middle.grow[
## Rules:
* a "`AUTHORS`" file **must** be present and contain your name(s)
* the source **must** be written in `C++11`
* the source **must** be documented with `Doxygen`
* the compilation chain **must** use `GNU Make`
* the presentation **must** be at the location "`pres/index.html`"
* using branches is OK but **only master will be checked**
* the presentation **must** use `remark.js`
* mandatory unit tests **must** pass

## :moneybag: Bonuses:
* no warning at {compilation, bindings, tests, doc}
* project is implementing other methods
* source code is well written (comments, easy to read, etc...)
* commits are small and explicit
]

---

# IPS-PROD-TP - Week #1 - Initialization

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style A fill:#090 
]

.column.middle.grow[

## 1. Fork your project from IPS-DEV

 * create a new empty project in your gitlab and add files from the previous project
 * or fork your previous project in your gitlab

## 2. Adapt your code to pass the mandatory unit tests

 * change method names
 * adapt object hierarchy
 * change methods arguments
 * ...
]

---

# IPS-PROD-TP - Week #1 - Polynomials

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style A fill:#090 
]


.column.middle.grow[
## Hermite polynomials (physicists' version)

`$$
\begin{eqnarray*}
H_0(\zeta)&=&1\\
H_1(\zeta)&=&2\zeta\\
H_{n}(\zeta)&=&2\zeta H_{n-1}(\zeta)-2(n-1)H_{n-2}(\zeta)
\end{eqnarray*}
$$`

## Generalized Laguerre polynomials

`$$
\begin{eqnarray*}
L^{m}_0(\eta)&=&1\\
L^{m}_1(\eta)&=&1+m-\eta\\
L^{m}_{n}(\eta)&=&\left(2+\frac{m-1-\eta}{n}\right)L^{m}_{n-1}(\eta)-\left(1+\frac{m-1}{n}\right)L^{m}_{n-2}(\eta)
\end{eqnarray*}
$$`
]

---

# IPS-PROD-TP - Week #1 - Mandatory unit test

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style A fill:#090 
]

.column.middle.grow[
```C++
// Mandatory test #00 - Hermite and Laguerre polynomials

Poly poly;

arma::vec zVals, calcVals, targetVals;

zVals = {-3.1, -2.3, -1.0, -0.3, 0.1, 4.3, 9.2, 13.7};
poly.calcHermite(6, zVals); // compute Hermite polynomials for n in {0 ... 5}

calcVals   = poly.hermite(4); // n = 4
targetVals = {  1.02835360e+03,  2.05825600e+02, -2.00000000e+01,  7.80960000e+00,
                1.15216000e+01,  4.59456160e+03,  1.10572154e+05,  5.54643458e+05};
TS_ASSERT_DELTA(arma::norm(calcVals / targetVals - 1.0), 0.0, 1e-08);

calcVals   = poly.hermite(5); // n = 5
targetVals = { -4.76676832e+03, -3.88909760e+02,  8.00000000e+00, -3.17577600e+01,
                1.18403200e+01,  3.48375818e+04,  1.98557479e+06,  1.50339793e+07};
TS_ASSERT_DELTA(arma::norm(calcVals / targetVals - 1.0), 0.0, 1e-08);

zVals = {0.1, 0.3, 1.2, 1.8, 2.0, 2.5, 7.1, 11.1};
poly.calcLaguerre(6, 4, zVals); // compute generalized Laguerre polynomials for m in {0 ... 5} and n in {0 ... 3}

calcVals   = poly.laguerre(4, 2); // m = 4, n = 2
targetVals = {  14.405,  13.245,  8.52 ,  5.82 ,  5.,  3.125,  -2.395,  10.005};
TS_ASSERT_DELTA(arma::norm(calcVals / targetVals - 1.0), 0.0, 1e-08);

calcVals   = poly.laguerre(5, 3); // m = 5, n = 3
targetVals = { 53.23983333,  47.95550000,  27.87200000,  17.5880,
               14.66666667,   8.39583333,  -0.81183333,  10.1015};
TS_ASSERT_DELTA(arma::norm(calcVals / targetVals - 1.0), 0.0, 1e-08);
```
]

---

# IPS-PROD-TP - Weeks #2-#5 - Project description

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style B fill:#090
    style C fill:#090
    style D fill:#090
    style E fill:#090
]

.column.middle.grow[
## Summary
During this 6-week long project, the local density of a nuclear system will be calculated and plotted.

.vspace[]

## Nuclear local density
Definition:
`$$\rho(\mathbf{r})\equiv \sum_a \sum_b \rho_{ab}\psi_a(\mathbf{r})\psi^*_b(\mathbf{r})$$`
]

---

# IPS-PROD-TP - Weeks #2-#5 - Project description

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style B fill:#090
    style C fill:#090
    style D fill:#090
    style E fill:#090
]


.column.middle.grow[
## Basis functions

They are defined as
`$$
\psi_{m,n,n_z}(r_\perp, \theta, z)
        \equiv
    Z(z, n_z)
    .
    R(r_\perp, m, n)
    .
         e^{im\theta}
$$`

with

`$$
Z(z, n_z)
\equiv
\phi_{n_z}(z)
=
         \frac{1}{\sqrt{b_z}}
         \frac{1}{\sqrt{2^{n_z} \sqrt{\pi}n_z!}}
         e^{-\frac{z^2}{2b_z^2}}H_{n_z}\left(\frac{z}{b_z}\right)
$$`

and
`$$
R(r_\perp, m, n)
         \equiv
         \frac{1}{b_{\perp}\sqrt{\pi}}
         \sqrt{\frac{n!}{(n+|m|)!}}
         e^{-\frac{r_{\perp}^2}{2b_{\perp}^2}}
         \left(\frac{r_{\perp}}{b_{\perp}}\right)^{|m|}
         L_n^{|m|}\left(\frac{r_{\perp}^2}{b_{\perp}^2}\right).
$$`

Parameters: `\(b_\perp, b_z\)`.
]

---

# IPS-PROD-TP - Weeks #2-#5 - Project description

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style B fill:#090
    style C fill:#090
    style D fill:#090
    style E fill:#090
]


.column.middle.grow[
## Basis truncation

We define
`$$
n_z^\textrm{max}(i) \equiv (N+2).Q^\frac{2}{3}+\frac{1}{2}-i.Q
$$`

and
`$$
m^\textrm{max} \equiv \textrm{sup}\left\{i:n_z^\textrm{max}(i)\ge 1\right\}.
$$`

The quantum number values are the **integers** verifying
`$$\begin{eqnarray}
0 &\le m \le& m^\textrm{max}-1\\
0 &\le n \le& \frac{1}{2}(m^\textrm{max}-m-1)\\
0 &\le n_z \le& n_z^\textrm{max}(m+2n+1)-1.
\end{eqnarray}$$`

Parameters: `\(N, Q\)`.
]


---

# IPS-PROD-TP - Weeks #2-#5 - Mandatory unit test

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style B fill:#090
    style C fill:#090
    style D fill:#090
    style E fill:#090
]

.column.middle.grow[
```C++
// Mandatory test #01 - Basis truncation

//     br = 1.935801664793151, bz = 2.829683956491218, N = 14, Q = 1.3
Basis basis(1.935801664793151,      2.829683956491218,     14,     1.3);

TS_ASSERT_EQUALS(basis.mMax, 14);

arma::ivec nMax = {7, 7, 6, 6, 5, 5, 4, 4, 3, 3, 2, 2, 1, 1};

TS_ASSERT(!arma::any(basis.nMax - nMax));

arma::imat n_zMax = {{18, 15, 13, 10, 7, 5, 2}, 
                     {16, 14, 11,  9, 6, 3, 1}, 
                     {15, 13, 10,  7, 5, 2, 0}, 
                     {14, 11,  9,  6, 3, 1, 0}, 
                     {13, 10,  7,  5, 2, 0, 0}, 
                     {11,  9,  6,  3, 1, 0, 0}, 
                     {10,  7,  5,  2, 0, 0, 0}, 
                     { 9,  6,  3,  1, 0, 0, 0}, 
                     { 7,  5,  2,  0, 0, 0, 0}, 
                     { 6,  3,  1,  0, 0, 0, 0}, 
                     { 5,  2,  0,  0, 0, 0, 0}, 
                     { 3,  1,  0,  0, 0, 0, 0}, 
                     { 2,  0,  0,  0, 0, 0, 0}, 
                     { 1,  0,  0,  0, 0, 0, 0}};

// check if matrices are equal
TS_ASSERT(!arma::any(arma::any(basis.n_zMax - n_zMax)));
```
]

---

# IPS-PROD-TP - Weeks #2-#5 - Mandatory unit test

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style B fill:#090
    style C fill:#090
    style D fill:#090
    style E fill:#090
]

.column.middle.grow[
```C++
// Mandatory test #02 - Basis r-functions

//     br = 1.935801664793151, bz = 2.829683956491218, N = 14, Q = 1.3
Basis basis(1.935801664793151,      2.829683956491218,     14,     1.3);

arma::vec r = {3.1, 2.3, 1.0, 0.0, 0.1, 4.3, 9.2, 13.7};

arma::vec res00 = { 8.08521235111303e-02,
                    1.43887615825118e-01,
                    2.55045100912706e-01,
                    2.91450097294984e-01,
                    2.91061479407116e-01,
                    2.47240792330589e-02,
                    3.63004153921473e-06,
                    3.87659726026123e-12};

TS_ASSERT_DELTA(arma::norm(basis.rPart(r, 0, 0) - res00), 0.0, 1e-15);

arma::vec res82 = { 5.87858442372438e-02,
                    1.35240488413384e-02,
                    4.06810074575519e-05,
                    0.00000000000000e+00,
                    4.92817669085478e-13,
                    8.52011998934850e-02,
                    5.20525909328609e-02,
                    1.44615166152252e-05};

TS_ASSERT_DELTA(arma::norm(basis.rPart(r, 8, 2) - res82), 0.0, 1e-15);
```
]

---

# IPS-PROD-TP - Weeks #2-#5 - Mandatory unit test

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style B fill:#090
    style C fill:#090
    style D fill:#090
    style E fill:#090
]

.column.middle.grow[
```C++
// Mandatory test #03 - Basis z-functions

//     br = 1.935801664793151, bz = 2.829683956491218, N = 14, Q = 1.3
Basis basis(1.935801664793151,      2.829683956491218,     14,     1.3);

arma::vec z = {-10.1, -8.4, -1.0, 0.0, 0.1, 4.3, 9.2, 13.7};

arma::vec res00 = { 7.64546544834383e-04,
                    5.44886272162148e-03,
                    4.19492564268520e-01,
                    4.46522724110539e-01,
                    4.46243982300708e-01,
                    1.40736821086932e-01,
                    2.26186220733178e-03,
                    3.62929640195959e-06};

TS_ASSERT_DELTA(arma::norm(basis.zPart(z, 0) - res00), 0.0, 1e-15);

arma::vec res15 = {-9.48674551049192e-02,
                   -1.40338701953237e-03,
                    1.85620628040096e-01,
                   -0.00000000000000e+00,
                   -3.93028470685214e-02,
                   -1.79526868763440e-01,
                    2.15604096600475e-01,
                    2.44977220882127e-01};

TS_ASSERT_DELTA(arma::norm(basis.zPart(z, 15) - res15), 0.0, 1e-15);
```
]

---

# IPS-PROD-TP - Weeks #2-#5 - Nuclear local density

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style B fill:#090
    style C fill:#090
    style D fill:#090
    style E fill:#090
]

.column.middle.grow[
## Direct algorithm
```C++
arma::mat result = arma::zeros(nbR, nbZ); // number of points on r- and z- axes
for (int m = 0; m < basis.mMax; m++)
{
  for (int n = 0; n < basis.nMax(m); n++)
  {
    for (int n_z = 0; n_z < basis.n_zMax(m, n); n_z++)
    {
      for (int mp = 0; mp < basis.mMax; mp++)
      {
        for (int np = 0; np < basis.nMax(mp); np++)
        {
          for (int n_zp = 0; n_zp < basis.n_zMax(mp, np); n_zp++)
          {
            arma::mat funcA = basis.basisFunc( m,  n,  n_z, zVals, rVals);
            arma::mat funcB = basis.basisFunc(mp, np, n_zp, zVals, rVals);
            result += funcA % funcB * rho(m, n, n_z, mp, np, n_zp); // mat += mat % mat * double
          }
        }
      }
    }
  }
}
```
]

---

# IPS-PROD-TP - Weeks #2-#5 - Nuclear local density

.mermaid[
    graph LR
    A("Week #1<br/>Initialization") --> B
    B("Week #2<br/>Project") --> C
    C("Week #3<br/>Project") --> D
    D("Week #4<br/>Project") --> E
    E("Week #5<br/>Project") --> F("Week #6<br/>Presentation")
    style B fill:#090
    style C fill:#090
    style D fill:#090
    style E fill:#090
]

.column.middle.grow[
## Work to be done

* calculate Hermite and Generalized Laguerre polynomials
* calculate z- and r- functions
* define the basis quantum number values
* calculate the nuclear local density with the direct algorithm
* calculate the nuclear local density with an **optimized algorithm** (see below)
* implement the mandatory unit tests
* plot the nuclear local density in the (x, z) plane (with Matplotlib, Plotly.js...)
* plot the nuclear local density in the (x, y, z) space (with POV-Ray, Blender...)
* present your results using the [presentation template](../template_ensiie/index.html)

## Possible optimizations of the direct algorithm

* use `\(\rho_{mnn_z,m'n'n'_z} = \delta_{mm'}\rho_{mnn_z,mn'n'_z}\)`
* use `\(\rho_{mnn_z,m'n'n'_z} = \rho_{m'n'n'_z,mnn_z}\)`
* optimize loops by using partial sums and pre-calculations
]

---

# POV-Ray - `df3` writer

```C++
std::string cubeToDf3(const arma::cube &m)
{
  std::stringstream ss(std::stringstream::out | std::stringstream::binary);
  int nx = m.n_rows;
  int ny = m.n_cols;
  int nz = m.n_slices;
  ss.put(nx >> 8);
  ss.put(nx & 0xff);
  ss.put(ny >> 8);
  ss.put(ny & 0xff);
  ss.put(nz >> 8);
  ss.put(nz & 0xff);
  double theMin = 0.0;
  double theMax = m.max();
  for (uint k = 0; k < m.n_slices; k++)
  {
    for (uint j = 0; j < m.n_cols; j++)
    {
      for (uint i = 0; i < m.n_rows; i++)
      {
        uint v = 255 * (fabs(m(i, j, k)) - theMin) / (theMax - theMin);
        ss.put(v);
      }
    }
  }
  return ss.str();
}
```
where `m` is an `arma::cube` representing the density to be plotted.
**Mandatory dimensions for `m`**:

.hcenter[
| Axis |    x   |    y   |    z   |
| :--: | :----: | :----: | :----: |
| min  | -10 fm | -10 fm | -20 fm |
| max  |  10 fm |  10 fm |  20 fm |
| nbp  |  32    |  32    |  64    |
]

---

# Blender - `raw` writer

```C++
std::string cubeToRaw(const arma::cube &m)
{
  std::stringstream ss(std::stringstream::out | std::stringstream::binary);
  double theMin = 0.0;
  double theMax = m.max();

  for (uint k = 0; k < m.n_slices; k++)
  {
    for (uint j = 0; j < m.n_cols; j++)
    {
      for (uint i = 0; i < m.n_rows; i++)
      {
        uint v = 255 * (fabs(m(i, j, k)) - theMin) / (theMax - theMin);
        ss.put(v);
      }
    }
  }

  return ss.str();
}
```
where `m` is an `arma::cube` representing the density to be plotted.
**Mandatory dimensions for `m`**:

.hcenter[
| Axis |    x   |    y   |    z   |
| :--: | :----: | :----: | :----: |
| min  | -10 fm | -10 fm | -20 fm |
| max  |  10 fm |  10 fm |  20 fm |
| nbp  |  32    |  32    |  64    |
]

---

class: middle

# Density file to be loaded

The file containing the values of the `\(\rho\)` matrix can be found [here](files/rho.arma).
To read it, use the `.load()` method of the `arma::mat` class:

```C++
arma::mat rho;
rho.load("rho.arma", arma::arma_ascii);
```

The corresponding basis has the following parameters:

```C++
//     br = 1.935801664793151, bz = 2.829683956491218, N = 14, Q = 1.3
Basis basis(1.935801664793151,      2.829683956491218,     14,     1.3);
```

:arrow_right: check that this basis contains 374 vectors !

The quantum numbers of the basis vector corresponding to the `\(i\)`-th row (or column) of the `\(\rho\)` matrix can be calculated with:

```C++
uint i = 0;
for (int m = 0; m < basis.mMax; m++)
for (int n = 0; n < basis.nMax(m); n++)
for (int n_z = 0; n_z < basis.n_zMax(m, n); n_z++)
{
  std::cout << "Basis vector " << i << ": m=" << m << " n=" << n << " n_z=" << n_z << std::endl;
  i++;
}
```

:arrow_right: You can check the symmetries of the `\(\rho\)` matrix:
* `\(\rho\)` should be symmetric : `\(\rho = \rho^t\)`
* `\(\rho\)` should be `\(m-\)`diagonal : `\(\rho(a, b) = \delta_{m_a, m_b} \rho(a, b)\)`

---

class: middle

# Presentation

## Directives

* the deadline for the projects is **01/09/2018 at 8:00pm**
* at the deadline, the branch `master` will be pulled
* check that your repository is **at the right address** in the following slides
* the slides **must be** in your repository (branch `master`, location `pres/index.html`)
* do not bring your slides for the presentation, they will already be on the screen
* if you want to plug a laptop to show something else than the slides, **check if it works before**
* be on time (5 min before the scheduled time)
* presentation format: **15 min + 5 min for questions**

.row[
.column.w48[
## :warning: Mandatory content :warning:

* explain the problem you have solved
* explain the associated physics content
* describe your project
* **show your results**
* **show the mandatory unit tests**
]
.column.w48[
## Possible content

* show some parts of your code
* show some benchmarks
* show some documentation
* show some of your `git` history
* show your `Makefile` organization
]
]

---

class: middle

# Presentations - 1st part

.hcenter[
| Presentation time | Students                                         | Gitlab account                                                                                                                         |
| :---------------: | :----------------------------------------------: | :------------------------------------------------------------------------------------------------------------------------------------: |
|  9:00             | Younès **Benallal**</br> Thomas **Dilasser**     | [younes.benallal](http://pedago-etu.ensiie.fr/thomas.dilasser/IPS-project)                                                             |
|  9:25             | Eléna **Bouchaud** </br> Hanna **Ben Jedidia**   | [elena.bouchaud-bougouin](http://pedago-etu.ensiie.fr/hanna.benjedidia/Bouchaud_BenJedidia_IPS)                                        |
|  9:50             | Yufei **Wang** </br> Fatma **Bali**              | [yufei.wang](http://pedago-etu.ensiie.fr/yufei.wang/IPS-PRO)                                                                           |
| 10:15             | Alexis **Puddu** </br> Domitille **Ailleret**    | [domitille.ailleret](http://pedago-etu.ensiie.fr/alexis.puddu/IPS/)                                                                    |
| 10:40             | Dinghao **Li** </br> Tingting **Wang**           | [tingting.wang](http://pedago-etu.ensiie.fr/tingting.wang/projet-IPS-prod.git)                                                         |
]

---

class: middle

# Presentations - 2nd part

.hcenter[
| Presentation time | Students                                         | Gitlab account                                                                                                                         |
| :---------------: | :----------------------------------------------: | :------------------------------------------------------------------------------------------------------------------------------------: |
| 11:05             | Aymeric **Marboeuf** </br> Pierre **Prevelle**   | [pierre.prevelle](http://pedago-etu.ensiie.fr/pierre.prevelle/prod)                                                                    |
| 11:30             | Louis **Melchior** </br> Tàzio **Gennuso**       | [tazio.gennuso](http://pedago-etu.ensiie.fr/louis.melchior/projet2IPS)                                                                 |
| 11:55             | Yael **Halimi** </br> Van man **Nguyen**         | [vanman.nguyen](http://pedago-etu.ensiie.fr/yael.halimi/Projet_IPS)                                                                    |
| 12:20             | Abdourahmane **Afifi** </br> Ruide **Song**      | [ruide.song](http://pedago-etu.ensiie.fr/ruide.song/ips-prod)                                                                          |
| 12:45             | Bylel **Trabelsi** </br> Théotime **Donnenfeld** | [theotime.donnenfeld](http://pedago-etu.ensiie.fr/theotime.donnenfeld/ProjetIPS2)                                                      |
]