-
Notifications
You must be signed in to change notification settings - Fork 4
Design
This documents holds designs descriptions and decisions.
- compilation: CMake
- unit-testing: googletest
- documentation: doxygen for now. Later on, we might want to complement it with breathe
- Exploration crap in exploration sub-directory
- likelihood source in likelihood sub-directory
- tests in test sub-directory of relevant source, e.g. likelihood/tests/.
General types, macros, and defines all go in DCProgsConfig.h.in in the root of the directory. If we need better separation at some point, lets think about it then.
Coding styles:
- tabwidth is 2
- column width is 100
- with the exception of
$Q$ , variables are all lower case - classes start with a capital and are camel-cased
- method and variables are lowercase, with _ as word separation
- private/protected members end with an "_"
- input variables in a method start with an "_"
In order to make exceptions thrown explicitely by DCProgs easily catch-able, it is best that all thrown exceptions derive from a single root. There are a lot of libraries out there for that purpose, notably boost::Exceptions. However, in order to keep dependencies minimal, it might be best to not go there just yet.
The exception hierarcy is currently in errors.h in the likelihood directory. There is no more structure than deriving from root:
namespace DCProgs {
namespace errors {
//! All explicit DCProgs exception derive from this, for easy catching,
class Root : public std::exception {};
//! Math (convergence, domain, etc) error
class Math : public Root { };
//! Math error which carries a message.
class Mass : public Math { ... };
//! Input error to a math problem
class Domain : public Math, virtual public std::domain_error { ... }
}
}
Please add exceptions to that file. Improve if you want a more explicit hierarchy.
The whole deal should go through a function of the form:
//! Returns the likelihood of a burst
//! \param[in] _gmat: Input extended or ideal G matrix defining the markov process.
//! \param[in] _series: series with open/shut times for a single burst
template<class T_GMATRIX, class T> t_real
likelihood(T_GMATRIX const &_gmat, Eigen::DenseBase<T> const &_series) {
// The folloing interface should exist
t_rmatrix X = _gmat.laplace_aa(...)
t_rmatrix X = _gmat.laplace_fa(...)
t_rmatrix X = _gmat.laplace_af(...)
t_rmatrix X = _gmat.laplace_ff(...)
}This function yields the likelihood. It is currently unimplemented. It acts on one of two objects: (i) the ideal transition matrix of the Markov process, (ii) the missed-even equivalent. These objects are defined with the following interface.
namespace DCProgs {
//! \brief Ideal transition matrix of open and shut intervals
//! \details Given a transition matrix $Q$ it is possible to figure out the evolution of any given
//! system.
class IdealG : protected StateMatrix {
public:
IdealG() : StateMatrix() {}
//! \brief Constructor with parameters.
//! \details Calls set method with input parameters.
//! \param[in] _matrix: Any matrix or matrix expression from Eigen. Will become the transition
//! matrix. Diagonal elements are transformed as explain in set(). Open
//! states should be in the top rows.
//! \param[in] _nopen: Number of open states.
//! \throws errors::Domain if input has incorrect values or size.
template<class T>
IdealG(Eigen::DenseBase<T> const &_matrix, t_int _nopen);
//! \brief Sets Q matrix and the number of open states.
//! \details Enforces \f[Q{ii} = -\sum_{j\neqi} Q{ij}]\f.
//! It is expected that open states are the top rows [0, _nopen].
void set(t_rmatrix const &_Q, t_int const &_nopen);
//! Gets Q matrix.
t_rmatrix const & get_Q() const { return this->matrix; }
//! Gets the number of open states
t_int const & get_nopen() const { return this->nopen; }
//! Open to open transitions.
auto aa(t_real t) const -> decltype( t_rmatrix::Zero(nopen, nopen) );
//! Shut to shut transitions.
t_rmatrix ff(t_real t) const;
//! Shut to open transitions.
auto fa(t_real t) const -> decltype( (t*StateMatrix::ff()).exp()*StateMatrix::fa() );
//! Open to shut transitions.
auto af(t_real t) const -> decltype( (t*StateMatrix::aa()).exp()*StateMatrix::af() );
//! Laplace transform of open to open transitions.
auto laplace_aa(t_real s) const -> decltype(aa(0)) { return this->aa(0); }
//! Laplace transform of shut to shut transitions.
auto laplace_ff(t_real s) const -> decltype(ff(0)) { return this->ff(0); }
//! Laplace transform of shut to open transitions.
t_rmatrix laplace_fa(t_real s) const;
//! Open to shut transitions.
t_rmatrix laplace_af(t_real t) const;
};
}The main interface are the laplace_... functions that will be called to actually compute the
likelihoods. The class derives from StateMatrix, a tuple (more or less) which holds the matrix of
transition rates for the mechanism, and a the number of open states. The open states are always at
the top of the matrix. In any case, StateMatrix is defined as:
namespace DCProgs {
//! \brief State matrix that can be partitioned into open/shut states.
//! \details In practice, this is a two tuple with some helper functions to get corners.
struct StateMatrix {
//! Number of open states.
t_int nopen;
//! The matrix itself.
t_rmatrix matrix;
//! Constructor
StateMatrix() : matrix(0,0), nopen(0) {}
//! Constructor
template<class T>
StateMatrix (Eigen::DenseBase<T> const &_c, t_int _nopen = 0)
: matrix(_c), nopen(_nopen) {}
//! Open to open transitions.
auto aa() const -> decltype( matrix.topLeftCorner(nopen, nopen) );
//! Open to shut transitions.
auto af() const -> decltype( matrix.topRightCorner(nopen, nopen) );
//! Shut to open transitions.
auto fa() const -> decltype( matrix.bottomLeftCorner(nopen, nopen) );
//! Shut to shut transitions.
auto ff() const -> decltype( matrix.bottomRightCorner(nopen, nopen) );
};
}NOTE: Currently only in branch features/recursion.
This feature is implemented in likelihood/recursion.h. The objective is to offer a way to compute
the
template<class T, class T_ZERO>
typename T::t_element recursion_formula( T & _C, t_int _i, t_int _m, t_int _l,
T_ZERO const &_zero );This approach means that we can test the implementation of the recursion using less complex object
(e.g. where each
The first template argument, T & _C, is an object with a specific interface. In practice, it
stands for the whole set of _C can somehow compute prior _C must provide is as follows:
class T {
//! Type of the element
typedef t_element;
//! Returns (prior) element in recursion
t_element operator()(t_int _i, t_int _j, t_int _m);
//! \brief Returns D objects, e.g. \f$A_{iAF}e^{Q_{FF}\tau}Q_{FA}\f$.
auto getD(t_int _i) const;
//! Returns specific eigenvalue of \f$Q\f$.
auto get_eigval(t_int _i) const;
//! Returns number of eigenvalues.
t_int nbeigval(t_int _i) const;
};The operator() should return
The second template argument to recursion_formula should be a
functor returning an appropriate zero: e.g. a matrix (or Eigen expression) if T::t_element is a
matrix, or simply the real number zero if T:t_element is a real number.
- calculation of the equilibrium vectors in equilibrium.h. There are currently two interfaces, one for c++, and one that can be accesed from python via ctypes, as done in exploration/Time series.ipynb (branch features/time_filter).