PhaseGradients.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2004 - 2016 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// psrchive/More/MEAL/MEAL/PhaseGradients.h
 
#ifndef __MEAL_PhaseGradients_H
#define __MEAL_PhaseGradients_H
 
#include "MEAL/Univariate.h"
#include "MEAL/Complex2.h"
#include "MEAL/Parameters.h"
 
namespace MEAL {
 
 
  template<typename T>
  class PhaseGradients : public Univariate<T>
  {
 
  public:
 
    PhaseGradients (unsigned ngradient = 0);
 
    PhaseGradients (const PhaseGradients&);
 
    PhaseGradients& operator = (const PhaseGradients&);
 
    PhaseGradients* clone () const;
 
    unsigned get_ngradient () const;
 
    void set_igradient (unsigned igradient);
 
    unsigned get_igradient () const;
 
    void set_offset (unsigned igradient, double offset);
 
    double get_offset (unsigned igradient) const;
 
    void add_gradient ();
 
    void remove_gradient ();
 
    void resize (unsigned ngradient);
 
    // ///////////////////////////////////////////////////////////////////
    //
    // Function implementation
    //
    // ///////////////////////////////////////////////////////////////////
 
    std::string get_name () const;
 
  protected:
 
    void calculate (typename T::Result& result, std::vector<typename T::Result>* gradient);
   
    Parameters parameters;
 
    std::vector<double> offsets;
 
    unsigned igradient;
 
  };
 
}
 
 
template<typename T>
MEAL::PhaseGradients<T>::PhaseGradients (unsigned ncoef)
  : parameters (this)
{
  resize (ncoef);
  igradient = 0;
}
 
template<typename T>
MEAL::PhaseGradients<T>::PhaseGradients (const PhaseGradients& copy)
  : parameters (this)
{
  operator = (copy);
}
 
template<typename T>
MEAL::PhaseGradients<T>& 
MEAL::PhaseGradients<T>::operator = (const PhaseGradients& copy)
{
  if (&copy == this)
    return *this;
 
  parameters = copy.parameters;
  offsets = copy.offsets;
  igradient = copy.igradient;
 
  this->set_evaluation_changed();
 
  return *this;
}
 
template<typename T>
MEAL::PhaseGradients<T>* MEAL::PhaseGradients<T>::clone () const
{
  return new PhaseGradients( *this );
}
 
template<typename T>
void MEAL::PhaseGradients<T>::set_igradient (unsigned i)
{
  if (i == igradient)
    return;
 
  if (i >= get_ngradient())
    throw Error (InvalidParam, "MEAL::PhaseGradients<T>::set_igradient",
                 "invalid index=%d >= ngradient=%d", i, get_ngradient());
 
  igradient = i;
  this->set_evaluation_changed ();
}
 
template<typename T>
unsigned MEAL::PhaseGradients<T>::get_igradient () const
{
  return igradient;
}
 
template<typename T>
void MEAL::PhaseGradients<T>::set_offset (unsigned igradient, double offset)
{
#if 0
  cerr << "MEAL::PhaseGradients<T>::set_offset"
    " i=" << igradient << " offset=" << offset << endl;
#endif
  offsets.at(igradient) = offset;
}
 
template<typename T>
double MEAL::PhaseGradients<T>::get_offset (unsigned igradient) const
{
  return offsets.at(igradient);
}
 
template<typename T>
void MEAL::PhaseGradients<T>::resize (unsigned ngradient)
{
  unsigned current = get_ngradient();
 
  if (current == ngradient)
    return;
 
  parameters.resize( ngradient );
  offsets.resize (ngradient);
 
  for (unsigned i=current; i<ngradient; i++)
  {
    parameters.set_name (i, "phi_" + tostring(i));
    offsets[i] = 0.0;
  }
 
  if (ngradient)
    set_igradient (ngradient-1);
}
 
template<typename T>
void MEAL::PhaseGradients<T>::add_gradient ()
{
  resize ( get_ngradient() + 1 );
}
 
template<typename T>
void MEAL::PhaseGradients<T>::remove_gradient ()
{
  resize ( get_ngradient() - 1 );
}
 
template<typename T>
unsigned MEAL::PhaseGradients<T>::get_ngradient () const
{
  return parameters.get_nparam();
}
 
template<typename T>
std::string MEAL::PhaseGradients<T>::get_name () const
{
  return "PhaseGradients";
}
 
template<typename T>
void MEAL::PhaseGradients<T>::calculate (typename T::Result& result,
                                         std::vector<typename T::Result>* grad)
{
  double x = this->get_abscissa();
 
  double phase = (this->get_param(igradient) + this->get_offset(igradient)) * x;
 
  if (this->verbose)
    std::cerr << "MEAL::PhaseGradients<T>::calculate phase=" << phase
              << std::endl;
 
  double cos_phase = ::cos(phase);
  double sin_phase = ::sin(phase);
 
  result = std::complex<double>(cos_phase, sin_phase);
 
  if (grad)
  {
    for (unsigned i=0; i<grad->size(); i++)
      (*grad)[i] = 0.0;
 
    (*grad)[igradient] = x * std::complex<double>(-sin_phase, cos_phase);
    
    if (this->verbose)
      std::cerr << "MEAL::PhaseGradients<T>::calculate gradient" << std::endl
                << "   " << (*grad)[igradient] << std::endl;
  }
}
 
#endif // __MEAL_PhaseGradients_H