nmrLinearRegression.h

Go to the documentation of this file.

/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-    */
/* ex: set filetype=cpp softtabstop=4 shiftwidth=4 tabstop=4 cindent expandtab: */
 
/*
  Author(s):  Peter Kazanzides
  Created on: 2016-04-29
 
  (C) Copyright 2016 Johns Hopkins University (JHU), All Rights Reserved.
 
--- begin cisst license - do not edit ---
 
This software is provided "as is" under an open source license, with
no warranty.  The complete license can be found in license.txt and
http://www.cisst.org/cisst/license.txt.
 
--- end cisst license ---
 
*/
 

 
#pragma once
 
#ifndef _nmrLinearRegression_h
#define _nmrLinearRegression_h
 
#include <vector>
#include <cisstCommon/cmnTypeTraits.h>
#include <cisstVector/vctFixedSizeVectorTypes.h>
#include <cisstVector/vctDynamicVectorTypes.h>
#include <cisstVector/vctForwardDeclarations.h>
 
// Always the last file to include!
// (Not really needed since everything is templated)
#include <cisstNumerical/nmrExport.h>

 
template <class _elementType>
class nmrLinearRegressionSolver
{
public:
    typedef _elementType ElementType;
    typedef typename cmnTypeTraits<_elementType>::VaArgPromotion SummationType;

    nmrLinearRegressionSolver(_elementType tol = cmnTypeTraits<_elementType>::DefaultTolerance)
        : numpts(0), Sx(0), Sy(0), Sxx(0), Sxy(0), Syy(0), tolerance(tol) {}
    virtual ~nmrLinearRegressionSolver() {}

    size_t NumPoints() const { return numpts; }

    _elementType GetTolerance() const { return tolerance; }

    void SetTolerance(_elementType tol) { tolerance = tol; }

    virtual void Clear()
    { numpts = 0; Sx = Sy = Sxx = Sxy = Syy = 0; }

    virtual bool Sample(const _elementType &x, const _elementType &y);

    virtual bool Sample(const vctFixedSizeVector<_elementType, 2> &in)
    { return Sample(in.X(), in.Y()); }

    virtual bool Sample(const vctFixedSizeConstVectorRef<_elementType, 2, 1> &in)
    { return Sample(in.X(), in.Y()); }

    virtual bool Sample(const vctDynamicVector<_elementType> &x,
                        const vctDynamicVector<_elementType> &y)
    {   return SampleVector<vctDynamicVectorOwner<_elementType> >(x,y); }
 
    virtual bool Sample(const vctDynamicConstVectorRef<_elementType> &x,
                        const vctDynamicConstVectorRef<_elementType> &y)
    {   return SampleVector<vctDynamicVectorRefOwner<_elementType> >(x,y); }

    virtual bool Sample(const std::vector<_elementType> &x, const std::vector<_elementType> &y)
    {
        // C++11 introduces std::vector::data, which would allow x.data() instead of &x.front()
        return SampleVector(vctDynamicConstVectorRef<_elementType>(x.size(), &x.front()),
                            vctDynamicConstVectorRef<_elementType>(y.size(), &y.front()));
    }

    virtual bool EstimateAsFractions(SummationType &slope_num, SummationType &yint_num,
                                     SummationType &denom, SummationType *tse_num = 0);

    virtual bool Estimate(_elementType &slope, _elementType &yint, _elementType *mse = 0);
 
protected:
    size_t numpts;                // Number of data points
    SummationType Sx, Sy;         // Sums of X and Y values
    SummationType Sxx, Sxy, Syy;  // Sums of X*X, X*Y, and Y*Y
    SummationType tolerance;      // Tolerance for avoiding divide by 0

    virtual bool ComputeSums(const _elementType &x, const _elementType &y);

    template <class _vectorOwnerType>
    bool SampleVector(const vctDynamicConstVectorBase<_vectorOwnerType, _elementType> &x,
                      const vctDynamicConstVectorBase<_vectorOwnerType, _elementType> &y)
    {
        if (x.size() != y.size()) return false;
        numpts += x.size();
        Sx += x.SumOfElements();
        Sy += y.SumOfElements();
        Sxx += x.DotProduct(x);
        Syy += y.DotProduct(y);
        Sxy += x.DotProduct(y);
        return true;
    }
 
};
 
template <class _elementType>
bool nmrLinearRegressionSolver<_elementType>::ComputeSums(const _elementType &x, const _elementType &y)
{
    Sx += x;
    Sy += y;
    Sxx += x*x;
    Sxy += x*y;
    Syy += y*y;
    return true;
}
 
template <class _elementType>
bool nmrLinearRegressionSolver<_elementType>::Sample(const _elementType &x, const _elementType &y)
{
    // ComputeSums always returns true, but in future could check for overflow
    bool ret = ComputeSums(x,y);
    if (ret) numpts++;
    return ret;
}
 
template <class _elementType>
bool nmrLinearRegressionSolver<_elementType>::EstimateAsFractions(SummationType &slope_num,
                                                                  SummationType &yint_num, SummationType &denom,
                                                                  SummationType *tse_num)
{
    denom = static_cast<SummationType>(numpts)*Sxx - Sx*Sx;
    slope_num = static_cast<SummationType>(numpts)*Sxy - Sx*Sy;
    yint_num  = Sxx*Sy - Sx*Sxy;
    // Terms in equation for tse_num are grouped to try to reduce likelihood of fixed point overflow
    if (tse_num)
        *tse_num = (Syy*denom - static_cast<SummationType>(numpts)*Sxy*Sxy) - Sy*(yint_num - Sx*Sxy);
    return true;
}
 
template <class _elementType>
bool nmrLinearRegressionSolver<_elementType>::Estimate(_elementType &slope, _elementType &yint,
                                                       _elementType *mse)
{
    SummationType slope_num, yint_num, denom, tse_num;
    EstimateAsFractions(slope_num, yint_num, denom, &tse_num);
 
    // Following check is to avoid division by 0 or by very small number.
    // Note that for fixed point numbers, the default tolerance is 0 so we
    // also check whether we are equal to the tolerance.
    if ((denom >= -tolerance) && (denom <= tolerance))
        return false;
 
    slope = static_cast<_elementType>(slope_num/denom);
    yint  = static_cast<_elementType>(yint_num/denom);
    if (mse)
        *mse = static_cast<_elementType>(tse_num/(denom*static_cast<SummationType>(numpts)));
    return true;
}

 
template <class _vectorType>
bool nmrLinearRegression(const _vectorType &x, const _vectorType &y,
                         typename _vectorType::value_type &slope, typename _vectorType::value_type &yint,
                         typename _vectorType::value_type *mse = 0,
                         typename _vectorType::value_type tolerance = cmnTypeTraits<typename _vectorType::value_type>::DefaultTolerance)
{
    nmrLinearRegressionSolver<typename _vectorType::value_type> solver(tolerance);
    bool ret = solver.Sample(x,y);
    if (ret)
       ret = solver.Estimate(slope, yint, mse);
    return ret;
}

 
template <class _elementType>
class nmrLinearRegressionWindowSolver : public nmrLinearRegressionSolver<_elementType>
{
protected:
    typedef nmrLinearRegressionSolver<_elementType> BaseClass;
    typedef typename BaseClass::SummationType SummationType;
    typedef typename std::vector<vctFixedSizeVector<_elementType, 2> > WindowType;
    typedef typename std::vector<vctFixedSizeVector<_elementType, 2> >::iterator WindowTypeIterator;
    WindowType window;
    WindowTypeIterator iter;
public:
    nmrLinearRegressionWindowSolver(size_t length) : window(length)
    { iter = window.begin(); }
    ~nmrLinearRegressionWindowSolver() {}
 
    virtual size_t WindowLength() const { return window.size(); }
 
    bool Sample(const _elementType &x, const _elementType &y)
    {
        // Update the window
        *iter++ = vctFixedSizeVector<_elementType, 2>(x,y);
        if (iter == window.end())
           iter = window.begin();
        // Increment numpts up to size of window
        if (this->numpts < window.size())
            this->numpts++;
        return true;
    }
 
    // The methods that take vector inputs are not currently implemented, as they
    // are less likely to be useful with a moving window.
    bool Sample(const vctDynamicVector<_elementType> & CMN_UNUSED(x),
                const vctDynamicVector<_elementType> & CMN_UNUSED(y))
    { return false; }
    bool Sample(const vctDynamicConstVectorRef<_elementType> & CMN_UNUSED(x),
                const vctDynamicConstVectorRef<_elementType> & CMN_UNUSED(y))
    { return false; }
    bool Sample(const std::vector<_elementType> & CMN_UNUSED(x),
                const std::vector<_elementType> & CMN_UNUSED(y))
    { return false; }

    virtual bool Recalculate()
    {
        size_t num = this->numpts;
        BaseClass::Clear();
        for (size_t i = 0; i < num; i++)
            BaseClass::ComputeSums(window[i].X(), window[i].Y());
        this->numpts = num;
        return true;
    }
 
    bool EstimateAsFractions(SummationType &slope_num, SummationType &yint_num,
                             SummationType &denom, SummationType *tse_num = 0)
    {
        Recalculate();
        return BaseClass::EstimateAsFractions(slope_num, yint_num, denom, tse_num);
    }
 
};

 
template <class _elementType>
class nmrLinearRegressionWindowRecursiveSolver : public nmrLinearRegressionWindowSolver<_elementType>
{
    typedef nmrLinearRegressionWindowSolver<_elementType> BaseClass;
    typedef typename BaseClass::SummationType SummationType;
    typedef typename BaseClass::ElementType ElementType;
 
public:
    nmrLinearRegressionWindowRecursiveSolver(size_t length) : BaseClass(length) {}
    ~nmrLinearRegressionWindowRecursiveSolver() {}

    bool Sample(const _elementType &x, const _elementType &y)
    {
        if (this->numpts >= this->window.size()) {
            // First, remove the oldest sample
            ElementType old_x = this->iter->X();
            ElementType old_y = this->iter->Y();
            this->Sx -= old_x;
            this->Sy -= old_y;
            this->Sxx -= old_x*old_x;
            this->Sxy -= old_x*old_y;
            this->Syy -= old_y*old_y;
        }
        // Add the new value
        bool ret = BaseClass::ComputeSums(x,y);
        // Update the window (and increment numpts if needed)
        if (ret)
           ret = BaseClass::Sample(x,y);
        return ret;
    }

    bool EstimateAsFractions(SummationType &slope_num, SummationType &yint_num,
                             SummationType &denom, SummationType *tse_num = 0)
    {
        return nmrLinearRegressionSolver<_elementType>::EstimateAsFractions(slope_num, yint_num, denom, tse_num);
    }
 
};
 
#endif