ʻOhana
Population structure, admixture history, and selection using learning methods.
jade.controller.hpp
1 /* -------------------------------------------------------------------------
2  Ohana
3  Copyright (c) 2015-2020 Jade Cheng (\___/)
4  Jade Cheng <info@jade-cheng.com> (='.'=)
5  ------------------------------------------------------------------------- */
6 
7 #ifndef JADE_CONTROLLER_HPP__
8 #define JADE_CONTROLLER_HPP__
9 
10 #include "jade.likelihood.hpp"
11 #include "jade.settings.hpp"
12 #include "jade.simplex.hpp"
13 #include "jade.stopwatch.hpp"
14 
15 namespace jade
16 {
17  ///
18  /// A template for a class that performs the Nelder-Mead optimization. This
19  /// class implements shared features of the two concrete controllers of this
20  /// module, one that uses a Newick tree as input and another that does not.
21  ///
22  template <typename TValue>
24  {
25  public:
26  /// The value type.
27  typedef TValue value_type;
28 
29  /// The likelihood type.
31 
32  /// The matrix type.
34 
35  /// The options type.
37 
38  /// The settings type.
40 
41  /// The simplex type.
43 
44  /// The container type for the simplex.
46 
47  /// The exit condition for the simplex.
49 
50  /// The loc arguments for the simplex.
52 
53  ///
54  /// Reclaims resources, if any, used by the instance.
55  ///
56  inline virtual ~basic_controller()
57  {
58  }
59 
60  ///
61  /// Computes the objective function by decoding the specified
62  /// Nelder- Mead parameters into a covariance matrix, performing a
63  /// Cholskey square root, calculating the determinant, computing the
64  /// inverse, and calling the likelihood function.
65  ///
66  /// \return The negation of the log likelihood.
67  ///
69  const container_type & params) ///< The parameters.
70  {
71  static const auto inf = std::numeric_limits<value_type>::max();
72 
73  //
74  // The LAPACK routines need only the lower triangle stored in
75  // the matrix. If these routines are successful, the lower
76  // triangle is copied to the upper triangle before calculating
77  // the likelihood.
78  //
79  if (!_decode_lower(_c, params))
80  return inf;
81 
82  //
83  // Variance and covariance values must be positive; if one is not,
84  // return infinity to reject these parameters.
85  //
86  for (size_t row = 0; row < _rk; row++)
87  for (size_t col = 0; col <= row; col++)
88  if (_c(row, col) <= value_type(0))
89  return inf;
90 
91  //
92  // Compute the inverse and the log of the determinant. If this
93  // fails, the matrix is not positive semidefinite; return Infinity
94  // to indicate these are unacceptable parameters.
95  //
96  value_type log_c_det;
97  if (!_c.invert(log_c_det))
98  return inf;
99 
100  //
101  // The Nelder-Mead algorithm minimizes the objective function, so
102  // return the negation of the log-likelihood function.
103  //
104  return -_likelihood(_c, log_c_det);
105  }
106 
107  ///
108  /// Writes results to standard output and files.
109  ///
110  virtual void emit_results(
111  const options_type & opts, ///< The options.
112  const simplex_type & simplex, ///< The simplex.
113  const exit_condition_type &) ///< The context.
114  {
115  //
116  // Decode the final set of parameters (the optimized matrix).
117  //
118  _decode_lower(_c, simplex.get_vertex());
119  _c.copy_lower_to_upper();
120 
121  //
122  // Display the log-likelihood, which is the negative of the last
123  // value from the objective function.
124  //
125  std::cout
126  << "\nlog likelihood = "
127  << -simplex.get_objval()
128  << std::endl;
129 
130  //
131  // Display or save the covariance matrix.
132  //
133  if (opts.is_cout_specified())
134  {
135  const auto & cout = opts.get_cout();
136  std::cout << "Writing C matrix to " << cout << std::endl;
137 
138  std::ofstream out (cout);
139  if (!out.good())
140  throw error() << "failed to create matrix '" << cout << "'";
141 
143  out << _c;
144  }
145  else
146  {
147  std::cout << "[C Matrix]\n";
149  std::cout << _c << std::endl;
150  }
151  }
152 
153  ///
154  /// \return A reference to the covariance matrix.
155  ///
156  inline const matrix_type & get_c() const
157  {
158  return _c;
159  }
160 
161  ///
162  /// \return The rooted K value for this instance.
163  ///
164  inline size_t get_rk() const
165  {
166  return _rk;
167  }
168 
169  ///
170  /// Creates and returns the initial set of parameters for the Nelder-
171  /// Mead algorithm.
172  ///
173  /// \return The initial parameters for the Nelder-Mead algorithm.
174  ///
176 
177  ///
178  /// Logs information about the specified context for one iteration of
179  /// the Nelder-Mead algorithm.
180  ///
182  const log_args_type & log_args) ///< The Nelder-Mead arguments.
183  {
184  const auto lle = -log_args.simplex->get_objval();
185 
186  const auto dlle = log_args.iteration == 1
187  ? value_type(0)
188  : lle - _lle;
189 
190  std::ostringstream line;
191  line << log_args.iteration
192  << std::fixed << std::setprecision(6)
193  << '\t' << _iteration_time;
195  line << '\t' << dlle << '\t' << lle;
196  std::cout << line.str() << std::endl;
197 
198  _lle = lle;
199  _iteration_time = stopwatch();
200  }
201 
202  protected:
203  ///
204  /// Initializes a new instance of the class based on the specified
205  /// settings.
206  ///
208  const settings_type & settings) ///< The algorithm settings.
209  : _rk (settings.get_rf().get_height())
210  , _c (settings.get_c())
211  , _lle (0)
212  , _likelihood (settings.get_rf(), settings.get_mu())
213  , _iteration_time ()
214  {
215  assert(_rk > 0);
216  assert(_c.is_size(_rk, _rk));
217  }
218 
219  ///
220  /// Decodes the specified Nelder-Mead container and stores the result
221  /// into the lower triangle, including the diagonal, of the covariance
222  /// matrix.
223  /// \return True if successful; otherwise, false.
224  ///
225  virtual bool _decode_lower(
226  matrix_type & dst, ///< The covariance matrix.
227  const container_type & src) ///< The Nelder-Mead container.
228  = 0;
229 
230  private:
231  size_t _rk; // The rooted K value; i.e. K - 1.
232  matrix_type _c; // The covariance matrix.
233  value_type _lle; // The log-likelihood value.
234  likelihood_type _likelihood; // The likelihood functor.
235  stopwatch _iteration_time; // The time elapsed per iteration.
236  };
237 }
238 
239 #endif // JADE_CONTROLLER_HPP__
jade::basic_controller::basic_controller
basic_controller(const settings_type &settings)
Initializes a new instance of the class based on the specified settings.
Definition: jade.controller.hpp:207
jade::basic_controller::init_parameters
virtual container_type init_parameters()=0
Creates and returns the initial set of parameters for the Nelder- Mead algorithm.
jade::basic_controller
A template for a class that performs the Nelder-Mead optimization. This class implements shared featu...
Definition: jade.controller.hpp:24
jade::basic_controller::compute_objfunc
value_type compute_objfunc(const container_type &params)
Computes the objective function by decoding the specified Nelder- Mead parameters into a covariance m...
Definition: jade.controller.hpp:68
jade::basic_options< value_type >
jade::basic_controller::settings_type
basic_settings< value_type > settings_type
The settings type.
Definition: jade.controller.hpp:39
jade::basic_simplex::log_args
Arguments passed to the logging function.
Definition: jade.simplex.hpp:115
jade::basic_options::is_cout_specified
bool is_cout_specified() const
Definition: nemeco/jade.options.hpp:182
jade::basic_controller::_decode_lower
virtual bool _decode_lower(matrix_type &dst, const container_type &src)=0
Decodes the specified Nelder-Mead container and stores the result into the lower triangle,...
jade::basic_controller::log_iteration
void log_iteration(const log_args_type &log_args)
Logs information about the specified context for one iteration of the Nelder-Mead algorithm.
Definition: jade.controller.hpp:181
jade::basic_controller::exit_condition_type
simplex_type::exit_condition_type exit_condition_type
The exit condition for the simplex.
Definition: jade.controller.hpp:48
jade::basic_stopwatch< std::chrono::high_resolution_clock >
jade::basic_matrix< value_type >::set_high_precision
static void set_high_precision(std::ostream &out)
Sets the specified stream to scientific notation with a high precision.
Definition: jade.matrix.hpp:1154
jade::basic_controller::options_type
basic_options< value_type > options_type
The options type.
Definition: jade.controller.hpp:36
jade::basic_controller::likelihood_type
basic_likelihood< value_type > likelihood_type
The likelihood type.
Definition: jade.controller.hpp:30
jade::basic_simplex
A template for a class that implements the Nelder-Mead Simplex Method, which attempts to minimize an ...
Definition: jade.simplex.hpp:21
jade::basic_controller::get_rk
size_t get_rk() const
Definition: jade.controller.hpp:164
jade::basic_simplex::container_type
std::vector< value_type > container_type
The container type.
Definition: jade.simplex.hpp:27
jade::basic_likelihood< value_type >
jade::basic_matrix::is_size
bool is_size(const basic_matrix &other) const
Definition: jade.matrix.hpp:896
jade::basic_controller::container_type
simplex_type::container_type container_type
The container type for the simplex.
Definition: jade.controller.hpp:45
jade::basic_controller::log_args_type
simplex_type::log_args log_args_type
The loc arguments for the simplex.
Definition: jade.controller.hpp:51
jade::basic_controller::emit_results
virtual void emit_results(const options_type &opts, const simplex_type &simplex, const exit_condition_type &)
Writes results to standard output and files.
Definition: jade.controller.hpp:110
jade::basic_simplex::get_objval
const value_type & get_objval() const
Definition: jade.simplex.hpp:714
jade::basic_matrix::invert
bool invert(value_type &log_det)
Computes and stores the inverse of this matrix using the Cholesky square root method....
Definition: jade.matrix.hpp:796
jade::basic_settings
A template for a class encapsulating the settings provided to the optimizer.
Definition: cpax/jade.settings.hpp:22
jade::basic_simplex::get_vertex
const container_type & get_vertex() const
Definition: jade.simplex.hpp:734
jade::basic_controller::value_type
TValue value_type
The value type.
Definition: jade.controller.hpp:27
jade::basic_simplex::exit_condition::type
type
The exit condition.
Definition: jade.simplex.hpp:38
jade::basic_controller::~basic_controller
virtual ~basic_controller()
Reclaims resources, if any, used by the instance.
Definition: jade.controller.hpp:56
jade::basic_controller::matrix_type
basic_matrix< value_type > matrix_type
The matrix type.
Definition: jade.controller.hpp:33
jade::basic_matrix::copy_lower_to_upper
void copy_lower_to_upper()
Copies the elements in the lower triangle of the square matrix to the corresponding upper triangle el...
Definition: jade.matrix.hpp:266
jade::basic_error
A template for a class representing an exception thrown from this namespace.
Definition: jade.error.hpp:20
jade::basic_matrix< value_type >
jade::basic_controller::get_c
const matrix_type & get_c() const
Definition: jade.controller.hpp:156
jade::basic_simplex::log_args::simplex
const basic_simplex * simplex
The simplex instance.
Definition: jade.simplex.hpp:118
jade::basic_controller::simplex_type
basic_simplex< value_type > simplex_type
The simplex type.
Definition: jade.controller.hpp:42
jade::basic_options::get_cout
const std::string & get_cout() const
Definition: nemeco/jade.options.hpp:104
jade::basic_simplex::log_args::iteration
size_t iteration
The completed iteration.
Definition: jade.simplex.hpp:116