fit.hpp

Go to the documentation of this file.

/*
 * $Id: fit.hpp 281 2012-08-26 14:35:20Z jdl3 $
 * Copyright (C) 2011, 2012, 2020 John D Lamb
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
 
#ifndef CCGSL_FIT_HPP
#define CCGSL_FIT_HPP
 
#include<cmath>
#include<exception>
#include<gsl/gsl_fit.h>
#include"vector.hpp"
 
namespace gsl {
  namespace fit {
    inline int linear( double const* x, size_t const xstride, double const* y,
               size_t const ystride,
               size_t const n, double* c0, double* c1, double* cov00, double* cov01,
               double* cov11, double* sumsq ){
      return gsl_fit_linear( x, xstride, y, ystride, n, c0, c1, cov00, cov01, cov11, sumsq ); }
    
    inline int wlinear( double const* x, size_t const xstride, double const* w,
            size_t const wstride,
            double const* y, size_t const ystride, size_t const n, double* c0,
            double* c1,
            double* cov00, double* cov01, double* cov11, double* chisq ){
      return gsl_fit_wlinear( x, xstride, w, wstride, y, ystride, n, c0, c1, cov00, cov01,
                  cov11, chisq ); }
 
    inline int linear_est( double const x, double const c0, double const c1, double const cov00,
               double const cov01, double const cov11, double* y, double* y_err ){
      return gsl_fit_linear_est( x, c0, c1, cov00, cov01, cov11, y, y_err ); }
    
    inline int mul( double const* x, size_t const xstride, double const* y, size_t const ystride,
            size_t const n, double* c1, double* cov11, double* sumsq ){
      return gsl_fit_mul( x, xstride, y, ystride, n, c1, cov11, sumsq ); }
    
    inline int wmul( double const* x, size_t const xstride, double const* w,
             size_t const wstride,
             double const* y, size_t const ystride, size_t const n,
             double* c1, double* cov11, double* chisq ){
      return gsl_fit_wmul( x, xstride, w, wstride, y, ystride, n, c1, cov11, chisq ); }
 
    inline int mul_est( double const x, double const c1, double const cov11, double* y, double* y_err ){
      return gsl_fit_mul_est( x, c1, cov11, y, y_err ); } 
 
    // Versions without stride arguments
 
    inline int linear( double const* x, double const* y,
               size_t const n, double* c0, double* c1, double* cov00,
               double* cov01, double* cov11,
               double* sumsq ){
      return gsl_fit_linear( x, 1, y, 1, n, c0, c1, cov00, cov01, cov11, sumsq ); }
    
    inline int wlinear( double const* x, double const* w,
            double const* y, size_t const n, double* c0, double* c1,
            double* cov00, double* cov01, double* cov11, double* chisq ){
      return gsl_fit_wlinear( x, 1, w, 1, y, 1, n, c0, c1, cov00, cov01, cov11, chisq ); }
 
    inline int mul( double const* x, double const* y,
            size_t const n, double* c1, double* cov11, double* sumsq ){
      return gsl_fit_mul( x, 1, y, 1, n, c1, cov11, sumsq ); }
    
    inline int wmul( double const* x, double const* w,
             double const* y, size_t const n, double* c1, double* cov11,
             double* chisq ){
      return gsl_fit_wmul( x, 1, w, 1, y, 1, n, c1, cov11, chisq ); }
 
    // CCGSL additions
    
#ifndef DOXYGEN_SKIP
    inline int linear( gsl::vector const& x, gsl::vector const& y, double* c0, double* c1,
               double* cov00, double* cov01, double* cov11, double* sumsq,
               size_t const n = 0 ){
      return linear( x.data(), x.get()->stride, y.data(), y.get()->stride,
             n == 0 ? y.size() / y.get()->stride : n, c0, c1, cov00,
             cov01, cov11, sumsq ); }
 
    inline int wlinear( gsl::vector const& x, gsl::vector const& w,
            gsl::vector const& y, double* c0, double* c1, double* cov00,
            double* cov01,
            double* cov11, double* chisq, size_t const n = 0 ){
      return wlinear( x.data(), x.get()->stride, w.data(), w.get()->stride,
              y.data(), y.get()->stride, n == 0 ? y.size() / y.get()->stride : n, c0, c1,
              cov00, cov01, cov11, chisq ); }
    
    inline int mul( gsl::vector const& x, gsl::vector const& y, double* c1, double* cov11,
            double* sumsq, size_t const n = 0 ){
      return mul( x.data(), x.get()->stride, y.data(), y.get()->stride,
          n == 0 ? y.size() / y.get()->stride : n, c1, cov11, sumsq ); }
 
    inline int wmul( gsl::vector const& x, gsl::vector const& w,
             gsl::vector const& y, double* c1, double* cov11, double* chisq,
             size_t const n = 0 ){
      return wmul( x.data(), x.get()->stride, w.data(), w.get()->stride,
           y.data(), y.get()->stride, n == 0 ? y.size() / y.get()->stride : n,
           c1, cov11, chisq ); }
#endif /* DOXYGEN_SKIP */
 
    // Template versions
#ifndef DOXYGEN_SKIP
    template<typename ARRAY>
    inline int linear( ARRAY const& x, ARRAY const& y, double* c0, double* c1,
               double* cov00, double* cov01, double* cov11, double* sumsq,
               size_t const n = 0 ){
      return linear( x.data(), 1, y.data(), 1,
             n == 0 ? y.size() : n, c0, c1, cov00, cov01, cov11, sumsq ); }
 
    template<typename ARRAY>
    inline int wlinear( ARRAY const& x, ARRAY const& w,
            ARRAY const& y, double* c0, double* c1, double* cov00, double* cov01,
            double* cov11, double* chisq, size_t const n = 0 ){
      return wlinear( x.data(), 1, w.data(), 1, y.data(), 1, n == 0 ? y.size() : n, c0, c1,
              cov00, cov01, cov11, chisq ); }
    
    template<typename ARRAY>
    inline int mul( ARRAY const& x, ARRAY const& y, double* c1, double* cov11,
            double* sumsq, size_t const n = 0 ){
      return mul( x.data(), 1, y.data(), 1, n == 0 ? y.size() : n, c1, cov11, sumsq ); }
 
    template<typename ARRAY>
    inline int wmul( ARRAY const& x, ARRAY const& w,
             ARRAY const& y, double* c1, double* cov11, double* chisq,
             size_t const n = 0 ){
      return wmul( x.data(), 1, w.data(), 1, y.data(), 1, n == 0 ? y.size() : n,
           c1, cov11, chisq ); }
#endif
    // Template versions with no pointers to doubles
    template<typename ARRAY>
    inline int linear( ARRAY const& x, ARRAY const& y, double& c0, double& c1,
               double& cov00, double& cov01, double& cov11, double& sumsq,
               size_t const n = 0 ){
      if( 0 == n ){
    if( x.size() != y.size() )
      throw std::length_error("gsl::fit::linear(): x and y must have same size.");
      } else {
    if( x.size() != n)
      throw std::length_error("gsl::fit::linear(): if n ≠ 0, you must have x.size() = n.");
    if( y.size() != n)
      throw std::length_error("gsl::fit::linear(): if n ≠ 0, you must have y.size() = n.");
      }
      return linear( x.data(), 1, y.data(), 1,
             n == 0 ? y.size() : n, &c0, &c1, &cov00, &cov01, &cov11, &sumsq ); }
 
    template<typename ARRAY>
    inline int wlinear( ARRAY const& x, ARRAY const& w,
            ARRAY const& y, double& c0, double& c1, double& cov00, double& cov01,
            double& cov11, double& chisq, size_t const n = 0 ){
      if( 0 == n ){
    if( x.size() != y.size() )
      throw std::length_error("gsl::fit::linear(): x and y must have same size.");
    if( w.size() != y.size() )
      throw std::length_error("gsl::fit::linear(): w and y must have same size.");
      } else {
    if( x.size() != n)
      throw std::length_error("gsl::fit::linear(): if n ≠ 0, you must have x.size() = n.");
    if( y.size() != n)
      throw std::length_error("gsl::fit::linear(): if n ≠ 0, you must have y.size() = n.");
    if( w.size() != n)
      throw std::length_error("gsl::fit::linear(): if n ≠ 0, you must have w.size() = n.");
      }
      return wlinear( x.data(), 1, w.data(), 1, y.data(), 1, n == 0 ? y.size() : n, &c0, &c1,
              &cov00, &cov01, &cov11, &chisq ); }
    
    template<typename ARRAY>
    inline int mul( ARRAY const& x, ARRAY const& y, double& c1, double& cov11,
            double& sumsq, size_t const n = 0 ){
      if( 0 == n ){
    if( x.size() != y.size() )
      throw std::length_error("gsl::fit::linear(): x and y must have same size.");
      } else {
    if( x.size() != n)
      throw std::length_error("gsl::fit::linear(): if n ≠ 0, you must have x.size() = n.");
    if( y.size() != n)
      throw std::length_error("gsl::fit::linear(): if n ≠ 0, you must have y.size() = n.");
      }
      return mul( x.data(), 1, y.data(), 1, n == 0 ? y.size() : n, &c1, &cov11, &sumsq ); }
 
    template<typename ARRAY>
    inline int wmul( ARRAY const& x, ARRAY const& w,
             ARRAY const& y, double& c1, double& cov11, double& chisq,
             size_t const n = 0 ){
      if( 0 == n ){
    if( x.size() != y.size() )
      throw std::length_error("gsl::fit::linear(): x and y must have same size.");
    if( w.size() != y.size() )
      throw std::length_error("gsl::fit::linear(): w and y must have same size.");
      } else {
    if( x.size() != n)
      throw std::length_error("gsl::fit::linear(): if n ≠ 0, you must have x.size() = n.");
    if( y.size() != n)
      throw std::length_error("gsl::fit::linear(): if n ≠ 0, you must have y.size() = n.");
    if( w.size() != n)
      throw std::length_error("gsl::fit::linear(): if n ≠ 0, you must have w.size() = n.");
      }
      return wmul( x.data(), 1, w.data(), 1, y.data(), 1, n == 0 ? y.size() : n,
           &c1, &cov11, &chisq ); }
    
  }
}
#endif