<a name="RSF"></a>
# RBComb simulation framework
This project aims at developing a framework that can be used to simulate any project
that is to be undertaken on the RBComb platform. It is designed in modular fashion,
such that it is flexible enough to adapt easily to any given situation. The list of
contents is the following:

1. [ RBComb simulation framework. ](#RSF)
2. [ Implementing a new project. ](#INP)
2. [ Classes. ](#Classes)
   - [ Vec2. ](#vec2)
   - [ Diagonalizer. ](#diagonalizer)
   - [ DrumParameters. ](#drumparams)
   - [ DrumVariables. ](#drumvars)
   - [ Drum. ](#drum)
   - [ Force. ](#force)
   - [ Driver. ](#driver)
   - [ Coupler. ](#coupler)
   - [ LatticeGenerator. ](#latticegenerator)
     - [ Neighbour Ordering Convention. ](#noc)
     - [ Inexistent Neighbour Convention. ](#inc)
   - [ MatrixElementCalculator. ](#mec)
   - [ RK4Buffer. ](#rk4b)
   - [ Rk4Stepper. ](#rk4s)
   - [ SystemParameters. ](#sysparams)
   - [ Grabber. ](#grabber)
   - [ System. ](#system)



<a name="INP"></a>
## Implementing a new project
In order to implement a new project, follow these steps:

0. Create a new folder in `projects`.
   - Create subfolders `bin`, `include`, `lib` and `results`
   - Copy the `lib` contents into `projects/newproj/lib`
1. Identify force on a single drum
   - Depends on static parameters (w, m, a, etc.) and dynamic parameters and variables (x, v, etc.).
   - Coupling and driving parameters are dynamic.
2. Write custom versions of classes `DrumParameters` (everything that is static) and `DrumVariables`
(everything that is dynamic) to accomodate these. Inspiration in `lib/drum_variables.hpp` and
`lib/drum_parameters.hpp`.
4. Write custom `Grabber` that extracts the relevant data (inspiration in `include/grabber.hpp`)
3. Write custom child of `Force` to calculate the force.
4. Write custom children of `Coupler` and `Driver` to correctly update the drive and couplings.
5. If desired, write a custom child of `LatticeGenerator` to generate a custom lattice,
or use the generator `projects/braidingTightBinding/include/rbcomb_generator_braid.hpp` to
create the RBComb.
6. If desired, write a child of `MatrixElementCalculator` to calculate matrix elements
for the dynamic matrix.
7. Check the unit tests (`unit_tests`) or other projects to see how a `main.cpp` could be
constructed.

Hints for organizing and running simulations
1. Use a Makefile.
   - The framework should be compiled with c++ 2a.
2. Organize executables in `bin`
3. Store data and plots in subfolders of `results`



<a name="Classes"></a>
## Classes
The code is structured in an object oriented approach. The classes that likely will
not need to be adapted for a new situation are found in the `lib` folder. They are
described in the following. Note that qualifiers, references and the like are discarded
where it improves legibility. Consult the source files for more information.

<a name="vec2"></a>
### `Vec2` (vec2.hpp), 2-vector utility class
1. Template arguments
   - `value_t`: type of vector entries
2. Members
   - Access
     - `value_t x()`
       - returns x entry
     - `value_t y()`
       - returns y entry
     - `Vec2 normalized()`
       - returns normalized version of vector
     - `value_t r()`
       - returns length
     - `value_t phi()`
       - returns angle (`std::atan2` version of it)
   - Member functions
     - `value_t r_wrt(Vec2)`
       - returns length with origin at argument
     - `value_t phi_wrt(Vec2)`
       - returns angle with origin at argument
     - `value_t norm()`
       - returns norm
     - `value_t norm_sq()`
       - returns square of norm
   - Modifiers
     - `Vec2 normalize()`
       - normalizes the vector and returns it
     - `Vec2 rotate(Vec2, value_t)`
       - rotates the vector and returns it
   - Supported Operators, All of these work as one would expect
     - `*` with `Vec2` (inner product) and `value_t`
     - `/` with `value_t`
     - `+, -` with `Vec2`
     - All versions of `op=` of the above
     - `[]` with `std::size_t`
     - `<<` with `std::ostream`

<a name="diagonalizer"></a>
### `Diagonalizer` (diagonalizer.hpp), class to diagonalize symmetric Matrices
1. Member functions
   - `std::vector<double> ev(std::vector<double> mat, size_t N)`
     - returns eigenvalues of the symmetric matrix mat of linear size N
     - throws upon diagonalization failure
   - `std::pair<std::vector<double>, std::vector<double> > evv(std::vector<double> mat, size_t N)`
     - returns pair of (eigenvalues, eigenvectors) of the symmetric matrix mat of size N
     - throws upon diagonalization failure
2. Further developments
   - Only finding eigenvectors and -values in a certain range may be added later on

<a name="drumparams"></a>
### `DrumParameters`(drum_parameters.hpp), contains data members characterizing the static state of a drum

<a name="drumvars"></a>
### `DrumVariables` (drum_variables.hpp), contains data members characterizing the dynamic state of a drum

<a name="drum"></a>
### `Drum` (drum.hpp), represents a single drum top resonator
1. Template arguments
   - `value_t`: Scalar type
   - `params_t`: Drum parameters container type
   - `vars_t`: Drum variables container type
   - `sbuffer_t`: Stepper buffer container type
2. Access
   - `params_t get_parameters()`
     - returns the parameters, const and reference versions implemented
   - `vars_t get_variables()`
     - returns the variables, const and reference versions implemented
   - `sbuffer_t get_sbuffer()`
     - returns the stepper buffer, const and reference versions implemented
3. Modifiers
   - `void set_coupling_0(value_t)`
     - Sets coupling 0
   - `void set_coupling_1(value_t)`
     - Sets coupling 1
   - `void set_coupling_2(value_t)`
     - Sets coupling 2
   - `void set_drive(value_t)`
     - Sets central electrode coupling
4. Description
   - A drum is described by a set of (static) parameters (stiffness, mass, x-y position, etc),
which are to be stored in a container of type `params_t`. The variables (displacement, velocity,
electrode charges, etc.) are stored in a container of type `vars_t`. Example classes for
these two types are `lib/drum_parameters.hpp` and `lib/drum_variables.hpp`.
However, these containers likely need to be adapted to the situation at hand.
When time evolving, the stepper will use the container of
type `sbuffer_t` to store its intermediate results. Note that the default constructor of
this class is `delete`'d. It should be constructed from an object of type `params_t`.
5. Further developments
   - Abstract interfaces for `params_t` and `vars_t` could be added, but they would
be trivial.

<a name="force"></a>
### `Force` (force.hpp), force functional
1. Template arguments
   - `value_t`: Scalar type
   - `params_t`: Drum parameters type
   - `vars_t`: Drum variables type
   - `buffer_t`: Stepper buffer type
2. Virtual functions
   - `value_t operator()(drum_t drum, drum_t n1, drum_t n2, drum_t n3, value_t time)`
     - Returns force on `drum` at `time`, given its three neighbours `n1`, `n2`, `n3`
3. Description
   - This interface is a guide to complete implementation of a force functional. Any force
functional should derive from this class, but the child type should then be used in
order to avoid the vtable penalty.
   - The type `drum_t` is a `Drum` with the given
template arguments. Typically, this functional would make heavy use of the `Drum`
access members `get_parameters()` and `get_variables()`. The `time` argument of the
functional exists to fit special cases as well. The file `include/force_simple.hpp`
showcases how a real force functional could be written.

<a name="driver"></a>
### `Driver` (driver.hpp), calculate drive of drums
1. Template arguments
   - `value_t`: Scalar type
   - `drum_t`: Drum type
2. Virtual functions
   - `void precompute(value_t t_end, value_t dt, std::vector<drum_t> drum_vec)`
     - Called once at begin of `System` lifetime
   - `void step(value_t dt)`
     - Move in time by `dt`
   - `value_t operator()(size_t drum_index)`
     - Returns drive of drum `drum_index` (wrt `drum_vec`) at current time
3. Description
   - This interface is a guide to complete implementation of a drive calculation class. Any driver
class should derive from this class, but the child type should then be used in
order to avoid the vtable penalty.
   - The purpose of this class is to set the drive of each drum at specific times.
   - In the `precompute` function, this class is passed all information it could need about
the system. Hence it can in principle precompute all values for all drums and all times
of the simulation.
   - The member `step` is called to inform the `Driver` that time is advanced by the passed argument.
Note that an rk4 scheme advances time in steps of `dt/2`.
   - The functional should return the current drive on the drum with index passed as argument.
   - An example implementation of a `Driver` is shown in `include/driver_simple.hpp`.

<a name="coupler"></a>
### `Coupler` (coupler.hpp), calculate couplings between drums
1. Template arguments
   - `value_t`: Scalar type
   - `drum_t`: Drum type
2. Virtual functions
   - `void precompute(value_t t_end, value_t dt, std::vector<drum_t> drum_vec)`
     - Called once at begin of `System` lifetime
   - `void step(value_t dt)`
     - Move in time by `dt`
   - `value_t operator()(size_t drum_index, size_t neighbour_index)`
     - Returns coupling between drums `drum_index` and `neighbour_index` (wrt `drum_vec`) at current time
3. Description
   - This interface is a guide to complete implementation of a coupling calculation class. Any coupler
class should derive from this class, but the child type should then be used in
order to avoid the vtable penalty.

   - The purpose of this class is to set the coupling of each neighbouring pair of drums at specific times.

   - In the `precompute` function, this class is passed all information it could need about
the system. Hence it can in principle precompute all values for all drums and all times
of the simulation.

   - The member `step` is called to inform the `Coupler` that time is advanced by the passed argument.
Note that an rk4 scheme advances time in steps of `dt/2`.

   - The functional should return the current coupling between the two drums with indices passed as arguments.

   - An example implementation of a `Coupler` is shown in `include/coupler_simple.hpp`.

<a name="latticegenerator"></a>
### `LatticeGenerator` (lattice_generator.hpp), generates drum lattices
1. Template arguments
   - `value_t`: Scalar type
   - `params_t`: Drum parameters type
   - `vars_t`: Drum variables type
   - `sbuffer_t`: Stepper buffer type
2. Virtual functions
   - `std::pair<std::vector<drum_t>, std::vector<int> > operator()(params_t)`
     - Takes a `params_t`
     - Returns a pair that characterizes the generated lattice
       - a vector of drums `ds`
       - an adjacency vector of vectors `adj`, such that `ds[i]` and `ds[adj[i][0]]` are neighbours
     - All drums have the same `params_t`, except that the `position` members differ.
3. Description
   - An example child of the `LatticeGenerator` is shown in the file `include/rbcomb_generator.hpp`.
4. Further developments
   - In the future, there may be another overload for the functional. For example, it could either take an
`std::vector<params_t>` or an additional random number generator to construct the drums
differently.
5. Dependents
   - None
6. Typical dependencies
   - `params_t`
     - Existence of `position` member
<a name="noc"></a>
#### Neighbour ordering convention
An important note is the __convention of neighbour ordering__. Each drum has neighbours 0 thru 3.
For drums in different sublattices, these neighbours are:
- Sublattice 'A':
  - 0, `adj[0]`: straight down
  - 1, `adj[1]`: top left
  - 2, `adj[2]`: top right
- Sublattice 'B':
  - 0, `adj[0]`: straight up
  - 1, `adj[1]`: bottom right
  - 2, `adj[2]`: bottom left
Here _adj[]_ signifies the adjacency list of the given drum. Similarly, the couplings
_t0_ thru _t2_ in objects of type `params_t` should also respect this ordering. More
generally, whenever neighbours of a specific drum are ordered in some fashion, they
are assumed to respect the above convention. Note that with this convention,
neighbours see each other as the same neighbour index (the i-th neighbour of j
  sees j as its i-th neighbour). __Never violate this convention__.
<a name="inc"></a>
#### Inexistent neighbour convention
Another __convention__ concerns __inexistent neighbours__. For that purpose, this class should append
an auxiliary drum to the end of the drum vector. If neighbour i of a drum does not exist
(due to boundary, for example), the corresponding neighbour index will point to the
auxiliary drum, i.e. it will show an index `drum_vec.size()`. All couplings of this auxiliary
drum are to be kept at 0. This condition can then be applied in force calculation to avoid branching,
if one uses the couplings of the neighbours instead of the considered drum.

<a name="mec"></a>
### `MatrixElementCalculator` (matrix_element_calculator.hpp), calculates matrix elements
1. Template arguments
   - `value_t`: Scalar type
   - `params_t`: Drum parameters type
   - `vars_t`: Drum variables type
   - `drum_t`: Drum type
2. Virtual functions
   - `value_t operator()(size_t index, std::vector<drum_t>)`
     - Returns the diagonal element (index, index).
   - `value_t operator()(size_t index1, size_t index2, std::vector<drum_t>)`
     - Returns the coupling element (index1, index2), where these are each others neighbour 0
   - `value_t operator()(size_t index1, size_t index2, std::vector<drum_t>, int)`
     - Returns the coupling element (index1, index2), where these are each others neighbour 1
   - `value_t operator()(size_t index1, size_t index2, std::vector<drum_t>, int, int)`
     - Returns the coupling element (index1, index2), where these are each others neighbour 2
3. Description
   - The overload of the functional is done to avoid branching. When building the matrix, one calls
the individual functions correctly to accomodate the correct neighbours.
4. Dependents
   - `System`
     - Existence of the virtual functions
5. Typical dependencies
   - `params_t`
     - Class semantics
   - `vars_t`
     - Class semantics
   - `drum_t`
     - `get_parameters()`
     - `get_variables()`

<a name="rk4b"></a>
### `RK4Buffer` (rk4_buffer.hpp), holds Rk4Stepper intermediate results

<a name="rk4s"></a>
### `Rk4Stepper` (rk4_stepper.hpp), performs timesteps using rk4 scheme
1. Template arguments
   - `value_t`: Scalar type
   - `params_t`: Drum parameters container type
   - `vars_t`: Drum variables container type
   - `buffer_t`: Stepper buffer container type
   - `force_t`: Force functional type
2. Member functions
   - `void step_1(force_t, std::vector<drum_t>, std::vector<std::vector<int> >, value_t dt, value_t time)`
   - `void step_2(force_t, std::vector<drum_t>, std::vector<std::vector<int> >, value_t dt, value_t time)`
   - `void step_3(force_t, std::vector<drum_t>, std::vector<std::vector<int> >, value_t dt, value_t time)`
     - All of the above perform one step of a timestep, between successive steps certain
     other updates need to be taken care of.
     - Arguments: Force functional, Drum vector, Adjacency vector, time step, start time of current step

<a name="sysparams"></a>
### `SystemParameters` (system_parameters.hpp), holds system parameters
1. Template arguments
   - `coupler_t`: Coupler type
   - `driver_t`: Driver type
2. Public data members
   - `coupler`: The coupler_t object of the system
   - `driver`: The driver_t object of the system
   - `adjacency_vector`: A `std::vector<std::vector<int> >` representing the adjacency vector

<a name="grabber"></a>
### `Grabber` (grabber.hpp), grabs and saves data
1. Template arguments
   - `value_t`: Scalar type
   - `drum_t`: Drum type
2. Member functions
   - `void init(value_t t_end, value_t dt, std::vector<drum_t>, std::vector<std::vector<int> >)`
     - Called upon `System` lifetime start.
   - `bool grab(std::vector<drum_t>, value_t time)`
     - Grabs data from the drums.
     - May use current time to decide if to grab data.
     - Returns `true` if data was grabbed.
   - `bool save()`
     - Saves data to the specified files.
     - Is called by `System`.
     - Returns `true` on success.
3. Description
   - This class needs to be re-implemented whenever either `DrumParameters` or
   `DrumVariables` changes.
   - One should also implement a helper struct, as seen in `include/grabber.hpp`

<a name="system"></a>
### `System` (system.hpp), holds all parts together
1. Template arguments
   - `value_t`: Scalar type
   - `drum_t`: Drum type
   - `grabber_t`: Data exfiltrator type
   - `sysparams_t`: System parameters type
   - `force_t`: Force functional type
   - `coupler_t`: Coupler type
   - `driver_t`: Driver type
   - `stepper_t`: Stepper type
   - `matelecalc_t`: Matrix element calculator type
2. Member functions
   - `void simulate()`
     - Runs a simulation with the set parameters.
   - `void step()`
     - Performs one timestep with the set parameters.
   - `void reset_time()`
     - Resets the simulation time to zero.
   - `void set_step(value_t)`
     - Sets the timestep.
   - `bool save()`
     - Calls `grabber_t::save()`.
   - `std::vector<value_t> get_matrix(matelecalc_t)`
     - Returns the dynamic matrix as provided by the argument.