Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
Tempo2PhaseModel.cpp
1 /*
2  * The MIT License (MIT)
3  *
4  * Copyright (c) 2016 The SKA organisation (except where explicity noted below)
5  *
6  * Permission is hereby granted, free of charge, to any person obtaining a copy
7  * of this software and associated documentation files (the "Software"), to deal
8  * in the Software without restriction, including without limitation the rights
9  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10  * copies of the Software, and to permit persons to whom the Software is
11  * furnished to do so, subject to the following conditions:
12  *
13  * The above copyright notice and this permission notice shall be included in all
14  * copies or substantial portions of the Software.
15  *
16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22  * SOFTWARE.
23  */
24 #include "Tempo2PhaseModel.h"
25 #include "panda/Error.h"
26 #include "panda/CaseInsensitiveString.h"
27 #include <boost/algorithm/string/trim.hpp>
28 #include <boost/filesystem.hpp>
29 #include <boost/filesystem/fstream.hpp>
30 #include <cmath>
31 #include <limits>
32 
33 namespace {
34 
35 /*
36  * All functions in this anonymous namespace are complete rewrites of the
37  * tempo2 predictor code.
38  * copyright details of the original code can be found in:
39  *
40  * https://bitbucket.org/psrsoft/tempo2/overview
41  */
42 
43 /* Basic Chebyshev polynomial */
44 class Cheby2D
45 {
46  public:
47  Cheby2D()
48  : _coeff(nullptr)
49  , _nx(0U)
50  , _ny(0U)
51  {
52  }
53 
54  ~Cheby2D()
55  {
56  free(_coeff);
57  }
58 
59  void resize(unsigned x, unsigned y);
60 
61  // return the x, y coefficient
62  inline
63  long double& coeff(unsigned x, unsigned y) {
64  return _coeff[y*_nx+x];
65  }
66 
67  // evaluate the chebyshev polynomial
68  long double evaluate(long double x, long double y) const;
69 
70  // calculate the derived Cheby2D
71  void construct_x_derivative(Cheby2D& result) const;
72 
73  protected:
74  long double evaluate_1d(long double x, int index) const;
75 
76  private:
77  long double* _coeff;
78  unsigned _nx, _ny;
79 
80 };
81 
82 void Cheby2D::resize(unsigned x, unsigned y)
83 {
84  if(_coeff) {
85  free(_coeff);
86  _coeff = nullptr;
87  }
88  _nx = x;
89  _ny = y;
90  if(_nx && _ny)
91  _coeff = static_cast<long double*>(malloc(_nx*_ny*sizeof(long double)));
92 }
93 
94 long double Cheby2D::evaluate_1d(long double x, int index) const
95 {
96  // Basic 1D Clenshaw's formula applied to a row of coeffs
97  long double di = 0.0L;
98  long double di1 = 0.0L;
99  long double di2;
100  long double* coeffcurr = _coeff + (index + 1) * _nx - 1;
101  for(std::size_t i = _nx - 1; i > 0; --i) {
102  di2 = di1;
103  di1 = di;
104  di = 2.0L * x * di1 - di2 + *(coeffcurr--); // i.e. + c_ij
105  }
106  // compute di0, special case
107  return x * di - di1 + 0.5L * (*coeffcurr--);
108 }
109 
110 long double Cheby2D::evaluate(long double x, long double y) const
111 {
112  // Clenshaw's formula applied twice. The cheby polynomials
113  // in x are evaluated to provide the coefficients for the
114  // cheby polynomial in y
115 
116  long double y_coeff, dj, dj1, dj2;
117 
118  dj = dj1 = 0.0L;
119  for (int j = _ny - 1; j > 0; --j) {
120  y_coeff = evaluate_1d(x, j);
121  dj2 = dj1;
122  dj1 = dj;
123  dj = 2.0L * y * dj1 - dj2 + y_coeff;
124  }
125  return y * dj - dj1 + 0.5L * evaluate_1d(x, 0);
126 }
127 
128 void Cheby2D::construct_x_derivative(Cheby2D& derived_cheby) const
129 {
130  unsigned i=0;
131  assert(derived_cheby._nx == _nx && derived_cheby._ny == _ny);
132  assert(_nx>2);
133 
134  for (unsigned iy=0; iy < _ny; ++iy)
135  {
136  unsigned ioff = i+_nx;
137  derived_cheby._coeff[ioff-1] = 0.0;
138  derived_cheby._coeff[ioff-2] = 2 * (_nx-1) * _coeff[ioff-1];
139  unsigned ix=_nx-1;
140  do
141  {
142  derived_cheby._coeff[i+ix]= derived_cheby._coeff[i+ix+2] + 2 * (ix+1) * _coeff[i+ix+1];
143  } while(ix-- != 0);
144  i+=_nx;
145  }
146 }
147 
148 /* Chebyshev predictive phase model */
149 struct ChebyModel
150 {
151  public:
152  void load(boost::filesystem::ifstream& is, unsigned& line_count)
153  {
154  std::string line;
155  std::istringstream iss;
156  unsigned ix=0;
157  // lambda to skip whitespace/comments
158  auto next_line = [&]() -> bool {
159  while( std::getline(is, line) ) {
160  ++line_count;
161  boost::trim_left(line);
162  if(line.empty() || line.front() == '#') {
163  continue; // ignore comments
164  }
165  iss.str(line);
166  iss.clear();
167  return true;
168  }
169  return false;
170  };
171  auto parse_error = [&](std::string const& msg)
172  {
173  ska::panda::Error e(msg);
174  e << " on line " << line_count;
175  throw e;
176  };
177  unsigned nx=0;
178  unsigned ny=0;
179  ska::panda::CaseInsensitiveString keyword;
180  if(!next_line()) return;
181  iss >> keyword;
182  if(keyword == "ChebyModel") {
183  iss >> keyword;
184  if(keyword != "BEGIN") {
185  parse_error("expecting ChebyModel BEGIN");
186  }
187  }
188  while( next_line() ) {
189  iss >> keyword;
190  if(keyword == "ChebyModel") {
191  iss >> keyword;
192  if(keyword == "BEGIN") {
193  parse_error("unexpected ChebyModel BEGIN");
194  continue;
195  }
196  if(keyword == "END") {
197  cheby.construct_x_derivative(frequency_cheby);
198  return;
199  }
200  }
201  if(keyword == "PSRNAME") {
202  iss >> psrname;
203  continue;
204  }
205  if(keyword == "SITENAME") {
206  iss >> sitename;
207  continue;
208  }
209  if(keyword == "TIME_RANGE")
210  {
211  iss >> mjd_start;
212  iss >> mjd_end;
213  continue;
214  }
215  if(keyword == "FREQ_RANGE")
216  {
217  iss >> freq_start;
218  iss >> freq_end;
219  continue;
220  }
221  if(keyword == "DISPERSION_CONSTANT")
222  {
223  iss >> dispersion_constant;
224  continue;
225  }
226  if(keyword == "NCOEFF_TIME")
227  {
228  iss >> nx;
229  continue;
230  }
231  if (keyword == "NCOEFF_FREQ")
232  {
233  iss >> ny;
234  continue;
235  }
236  if (keyword == "COEFFS")
237  {
238  if(ix==0) {
239  if(!(nx&&ny))
240  parse_error("unexpected COEFFS keyword before NCOEFF_TIME or NCOEFF_FREQ");
241  cheby.resize(nx,ny);
242  frequency_cheby.resize(nx,ny);
243  }
244  unsigned iy=0;
245  while(iy<ny) {
246  if(iss.eof()) {
247  next_line();
248  }
249  iss >> cheby.coeff(ix, iy);
250  ++iy;
251  }
252  ++ix;
253  continue;
254  }
255  parse_error(std::string("unrecognised keyword :") + keyword.c_str());
256  } // while
257  }
258 
259  public:
260  std::string psrname;
261  std::string sitename;
262  long double mjd_start, mjd_end;
263  long double freq_start, freq_end;
264  long double dispersion_constant; // phase = polyco + d_c/f^2
265  Cheby2D cheby;
266  Cheby2D frequency_cheby;
267 };
268 
269 
270 class ChebyModelSet {
271  public:
272  ~ChebyModelSet() {
273  for(auto ptr : _models) {
274  delete ptr;
275  }
276  }
277  std::vector<ChebyModel*>& models() { return _models; };
278 
279  // the number of models is the set
280  std::size_t size() const { return _models.size(); }
281 
282  // return the evaluation of phase from the closest model mathcing the params
283  // N.B. ignoring DM in the model
284  long double phase(long double mjd, long double freq) const
285  {
286  ChebyModel const* model = get_nearest_model(mjd);
287  if(model) {
288  // linear fit across the model range
289  return model->cheby.evaluate(
290  -1.0L + 2.0L * (mjd - model->mjd_start)/(model->mjd_end - model->mjd_start)
291  ,
292  -1.0L + 2.0L * (freq - model->freq_start)/(model->freq_end - model->freq_start));
293  //+ model->dispersion_constant/(freq * freq);
294  }
295  return -1.0;
296  }
297 
298  void load(boost::filesystem::ifstream& is)
299  {
300  std::string line;
301  unsigned line_number=0;
302  ska::panda::CaseInsensitiveString keyword;
303  unsigned model_number;
304  while( std::getline(is, line) ) {
305  ++line_number;
306  boost::trim_left(line);
307  if(line.empty() || line.front() == '#') {
308  continue; // ignore comments
309  }
310  std::istringstream iss(line);
311  iss >> keyword >> model_number;
312  if(keyword != "ChebyModelSet") {
313  throw ska::panda::Error("expecting a ChebyModelSet format");
314  }
315  _models.reserve(model_number);
316  for(unsigned num =0; num < model_number; ++num)
317  {
318  ChebyModel* model = new ChebyModel();
319  _models.push_back(model);
320  model->load(is, line_number);
321  }
322  }
323  }
324 
325  protected:
327  ChebyModel const* get_nearest_model(long double mjd) const
328  {
329  ChebyModel const* rv_model = nullptr;
330  long double min_offset=std::numeric_limits<long double>::max(), offset;
331  for(ChebyModel const* model : _models) {
332  // check model for matching times
333  if(mjd < model->mjd_start || mjd > model->mjd_end) {
334  continue;
335  }
336  offset = std::abs((model->mjd_start+model->mjd_end)/2 - mjd);
337  if( offset < min_offset ) {
338  min_offset = offset;
339  rv_model = model;
340  }
341  }
342  return rv_model;
343  }
344 
345  private:
346  std::vector<ChebyModel*> _models;
347 };
348 
349 } // namespace
350 
351 namespace ska {
352 namespace cheetah {
353 namespace generators {
354 
355 struct T2Predictor {
356  public:
357  T2Predictor()
358  : _modelset(nullptr)
359  {
360  }
361  ~T2Predictor() {
362  delete _modelset;
363  }
364 
365  long double phase(long double mjd, long double freq) // freq in MHz
366  {
367  return _modelset->phase(mjd, freq);
368  }
369 
370  // load in model from file
371  void load(boost::filesystem::path const& filename) {
372  boost::filesystem::ifstream is(filename);
373  if(!is.is_open()) {
374  panda::Error e("unable to open file: ");
375  e << filename;
376  throw e;
377  }
378  if(_modelset == nullptr) _modelset = new ChebyModelSet;
379  _modelset->load(is);
380  }
381 
382  private:
383  ChebyModelSet* _modelset;
384 
385 };
386 
387 Tempo2PhaseModel::Tempo2PhaseModel(Tempo2PhaseModelConfig const& config)
388  : _pred(new T2Predictor)
389 {
390  load(config.filename());
391 }
392 
393 void Tempo2PhaseModel::load(boost::filesystem::path const& tempo2_predictor_file)
394 {
395  // check the file exists
396  if(!boost::filesystem::is_regular_file(tempo2_predictor_file))
397  {
398  panda::Error e("file does not exist:");
399  e << tempo2_predictor_file;
400  throw e;
401  }
402  PANDA_LOG << "tempo2 predictor model loading from file " << tempo2_predictor_file;
403 
404  _pred->load(tempo2_predictor_file);
405 
406 }
407 
408 Tempo2PhaseModel::~Tempo2PhaseModel()
409 {
410 }
411 
412 double Tempo2PhaseModel::operator()(utils::ModifiedJulianClock::time_point const& mjd, boost::units::quantity<data::MegaHertz, double> const mhz) const
413 {
414  long double const mjd_dbl = mjd.time_since_epoch().count();
415  return _pred->phase(mjd_dbl, mhz.value());
416 }
417 
418 } // namespace generators
419 } // namespace cheetah
420 } // namespace ska
configuration parameters for a tempo2 predictor model
Some limits and constants for FLDO.
Definition: Brdz.h:35
boost::filesystem::path const & filename() const
return the tempo2 file containing the model parameters