multifit_function.hpp

Go to the documentation of this file.

/*
 * $Id: multifit_function.hpp 233 2012-08-08 19:44:28Z jdl3 $
 * Copyright (C) 2010, 2011, 2012, 2020, 2021 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_MULTIFIT_FUNCTION_HPP
#define CCGSL_MULTIFIT_FUNCTION_HPP
 
#include<gsl/gsl_multifit_nlin.h>
#include"vector.hpp"
// For std::swap
#if __cplusplus <= 199711L
#include<algorithm>
#endif
 
namespace gsl {
  namespace multifit {
    class function : public gsl_multifit_function {
    public:
      struct Concept {
    virtual double f( gsl::vector const& x, gsl::vector& f ) = 0;
    virtual size_t size() const = 0;
    virtual size_t num_parameters() const = 0;
      };
#ifndef DOXYGEN_SKIP
    private:
      struct base_F {
    typedef int (*function_t)(gsl_vector const*,void*,gsl_vector*);
      };
      template<typename T>
      struct F : public base_F {
      public:
    F( T& t ) : t( t ){
      xv.wrap_gsl_vector_without_ownership( 0 );
      fv.wrap_gsl_vector_without_ownership( 0 );
    }
      public:
    static int fn( gsl_vector const* x, void* params, gsl_vector* fx ){
      // share gsl_vector* and gsl::vector contents.
      F<T>* ft = reinterpret_cast<F<T>*>( params );
      ft->xv.ccgsl_pointer = const_cast<gsl_vector*>( x );
      ft->fv.ccgsl_pointer = fx;
      return ft->t.f( ft->xv, ft->fv ); }
      private:
    T& t;
    gsl::vector xv;
    gsl::vector fv;
      };
      class Fn : public base_F {
      public:
    Fn( int (*const f)(gsl::vector const&,gsl::vector&), size_t n ) : f( f ), n( n ){
      xv.wrap_gsl_vector_without_ownership( 0 );
      fv.wrap_gsl_vector_without_ownership( 0 );
    }
    function_t function(){ return &fn; }
      private:
    int (*const f)(gsl::vector const&,gsl::vector&);
      public:
    size_t const n;
    static int fn( gsl_vector const* x, void* params, gsl_vector* fx ){
      // share gsl_vector* and gsl::vector contents.
      Fn* ft = reinterpret_cast<Fn*>( params );
      ft->xv.ccgsl_pointer = const_cast<gsl_vector*>( x );
      ft->fv.ccgsl_pointer = fx;
      return ft->f( ft->xv, ft->fv ); }
      private:
    gsl::vector xv;
    gsl::vector fv;
      };
      base_F* func;
      size_t* count;
#endif
    public:
      function(){
    func = 0;
    f = 0;
    n = 0;
    p = 0;
    params = 0;
    count = 0; // initially nullptr will do
      }
      template<typename T>
      void function_constructor( function&, T& );
      template<typename T>
      friend void function_constructor( function&, T& );
      explicit function( gsl_multifit_function& v ){
    // explicitly set function = 0 to indicate nothing to be deleted.
    func = 0;
    // copy
    f = v.f;
    n = v.n;
    p = v.p;
    params = v.params;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to function
      }
      template<typename T>
      function( T& t ){ function_constructor<T>( *this, t ); }
      function( int (* const f)(gsl::vector const&, gsl::vector&), size_t n ){
    func = 0;
    count = 0;
    func = new function::Fn( f, n );
    this->f = Fn::fn;
    this->n = n;
    this->p = p;
    params = func;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to function
      }
       
      // Copy and move constructors
      // copy constructor
      function( function const& v ) : func( v.func ), count( v.count ){
    f = v.f;
    n = v.n;
    p = v.p;
    params = v.params;
    if( count != 0 ) ++*count; // function is now shared.
      }
      // assignment operator
      function& operator=( function const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      delete func;
      delete count;
    }
    // Then copy
    func = v.func;
    f = v.f;
    n = v.n;
    p = v.p;
    params = v.params;
    // count = v.count;
    if( count != 0 ) ++*count; // function is now shared.
    return *this;
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      function( function&& v ) : func( v.func ), count( nullptr ){
    std::swap( f, v.f );
    std::swap( n, v.n );
    std::swap( p, v.p );
    std::swap( params, v.params );
    std::swap( count, v.count );
    v.func = nullptr;
      }
      function& operator=( function&& v ){
    std::swap( func, v.func );
    std::swap( f, v.f );
    std::swap( n, v.n );
    std::swap( p, v.p );
    std::swap( params, v.params );
    std::swap( count, v.count );
    return *this;
      }
#endif
      ~function(){ 
    // first, possibly delete anything pointed to by f
    if( count == 0 or --*count == 0 ){
      delete func;
      delete count;
    }
      }
    };
#ifndef DOXYGEN_SKIP
    template<typename T>
    inline void function_constructor( function& f, T& t ){
      f.func = 0;
      f.count = 0;
      // try constructing from T
      f.func = new function::F<T>( t );
      f.f = &function::F<T>::fn;
      f.n = t.size();
      f.p = t.num_parameters();
      f.params = f.func;
      // just plausibly we could fail to allocate count: no further action needed.
      f.count = new size_t;
      *f.count = 1; // initially there is just one reference to function
    }
    template<>
    inline void function_constructor( function& f, function& v ){
      f.func = v.func; f.count = v.count; f.f = v.f; f.n = v.n; f.p = v.p;
      f.params = v.params; if( f.count != 0 ) ++*f.count;
      /* function is now shared. */ }
    template<>
    inline void function_constructor( function& f, gsl_multifit_function& v ){
      f.f = v.f; f.n = v.n; f.p = v.p; f.params = v.params; }
#endif
    template<typename T>
    inline function make_function( T& t ){
      function fn( t );
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      return std::move( fn );
#else
      return fn;
#endif
    }
  }
}
 
#endif