whittaker.h

Go to the documentation of this file.

 
#ifndef __SABLIB_WHITTAKER_H__
#define __SABLIB_WHITTAKER_H__
 
#include <stdexcept>
#include <vector>
 
#include "../misc/diff.h"
 
namespace sablib {

const std::vector<double> Whittaker(
    const std::vector<double> & y, const std::vector<double> & w,
    const double lambda, const unsigned int s = 2
);

const std::vector<double> Whittaker(
    const std::vector<double> & y, const double lambda, const unsigned int s = 2
);

template <typename Derived>
const typename Derived::PlainObject
Whittaker(
    const Eigen::MatrixBase<Derived> & y,
    const Eigen::MatrixBase<Derived> & w,
    const Eigen::SparseMatrix<typename Derived::PlainObject::Scalar> & lambdaDTD
)
{
    // Although parameters are received as MatrixBase<Derived>, only vector classes are allowed.
    // Others will be rejected at compile time.
    static_assert(Derived::IsVectorAtCompileTime, "Error: y and w are not vector.");
 
    using Scalar = typename Derived::PlainObject::Scalar;
 
    Eigen::VectorXd z;
    Eigen::SparseMatrix<typename Derived::PlainObject::Scalar> W;
    Eigen::SimplicialCholesky< Eigen::SparseMatrix<Scalar> > solver;
 
    W = w.asDiagonal();
    solver.compute(W + lambdaDTD);
 
    if(solver.info() != Eigen::Success) {
        throw std::runtime_error("Whittaker(): solver calculation fails.");
    }
 
    z = solver.solve(W * y);
 
    return z;
}
 
}; // namespace sablib
 
#endif // __SABLIB_WHITTAKER_H__