overhaul of the function class

examples/exCompiledDoubleFunction.cpp

@@ -17,14 +17,14 @@ int main(int argc, char* argv[])
     }
     source = argv[1];
     
-    CompiledDoubleFunction f(1, source);
+    CompiledDoubleFunction f(source, 1);
     
     cout << "-- Program:" << endl << f << endl;
     cout << "-- Values:" << endl;
     for (double x = 0.0; x < 10.0; x += 0.5)
     {
         cout << "f(" << right << setw(6) << strFrom<double>(x) << ")= "
-             << f(x) << endl;
+             << f(&x) << endl;
     }
     
     return EXIT_SUCCESS;

examples/exDerivative.cpp

@@ -23,7 +23,7 @@ int main(int argc, char *argv[])
     maxOrder = strTo<Index>(argv[2]);
     x        = strTo<double>(argv[3]);
     
-    CompiledDoubleFunction f(1, source);
+    DoubleFunction    f = compile(source, 1);
     CentralDerivative df(f);
     
     for (Index i = 1; i <= 4; ++i)
@@ -32,7 +32,7 @@ int main(int argc, char *argv[])
         for (Index j = 0; j <= maxOrder; ++j)
         {
             df.setOrder(j, i);
-            cout << "d^" << j << "f(" << x << ")= " << df(x) << endl;
+            cout << "d^" << j << "f(" << x << ")= " << df(&x) << endl;
         }
     }
     

examples/exFit.cpp

@@ -15,10 +15,10 @@ const double exactPar[2] = {0.5,5.0}, dx = 10.0/static_cast<double>(nPoint);
 int main(void)
 {
     // generate fake data
-    XYStatData          data(nPoint, 1, 1);
+    XYStatData  data(nPoint, 1, 1);
-    RandGen             rg;
+    RandGen     rg;
-    double              x_k, y_k;
+    double      x_k, y_k;
-    CompiledDoubleModel f(1, 2, "return p_1*exp(-x_0*p_0);");
+    DoubleModel f = compile("return p_1*exp(-x_0*p_0);", 1, 2);
 
     for (Index k = 0; k < nPoint; ++k)
     {

examples/exIntegrator.cpp

@@ -20,9 +20,9 @@ int main(int argc, char* argv[])
     xMin   = strTo<double>(argv[2]);
     xMax   = strTo<double>(argv[3]);
     
-    CompiledDoubleFunction f(1, source);
+    DoubleFunction    f = compile(source, 1);
-    GslQagsIntegrator      integrator;
+    GslQagsIntegrator integrator;
-    double                 result;
+    double            result;
     
     result = integrator(f, xMin, xMax);
     cout << "function integral on [" << xMin << ", " << xMax << "] = ";

examples/exMin.cpp

@@ -17,12 +17,15 @@ int main(int argc, char* argv[])
     }
     source = argv[1];
     
-    CompiledDoubleFunction f(1, source);
+    DoubleFunction  f = compile(source, 1);
-    MinuitMinimizer minimizer(1);
+    MinuitMinimizer minimize(1);
+    DVec init(1);
     double min;
-    
+
-    minimizer.setVerbosity(Minimizer::Verbosity::Debug);
+    init(0) = 0.1;
-    min = minimizer(f)(0);
+    minimize.setInit(init);
+    minimize.setVerbosity(Minimizer::Verbosity::Debug);
+    min = minimize(f)(0);
     cout << "function minimum = " << min << endl;
     
     return EXIT_SUCCESS;

examples/exPlot.cpp

@@ -20,11 +20,10 @@ int main(void)
     }
     p << PlotRange(Axis::x, -5.0, 5.0) << PlotRange(Axis::y, -5.0, 20.0);
     p << Color("rgb 'blue'") << PlotFunction(StdMath::tgamma, -5, 5);
-    p << PlotFunction(CompiledDoubleFunction(1, "return cos(x_0)^2;"), -5, 5);
+    p << PlotFunction(compile("return cos(x_0)^2;", 1), -5, 5);
     p << Color("rgb 'brown'") << PlotCommand("x**3");
     p << PlotCommand("x**2");
-    p << PlotFunction(TabFunction(x, y), 0.,
+    p << PlotFunction(interpolate(x, y), 0., nPoint/2.0 - 1.1);
-                      static_cast<double>(nPoint)/2.0 - 1.1);
     cout << p << endl;
     p.display();
     

examples/exRootFinder.cpp

@@ -8,8 +8,8 @@ using namespace Latan;
 int main(void)
 {
     constexpr double a = 1., b = 10.;
-    DoubleFunction f1(2, [a](const double *x){return a*(1.-x[0]);});
+    DoubleFunction f1([a](const double *x){return a*(1.-x[0]);}, 2);
-    DoubleFunction f2(2, [b](const double *x){return b*(x[1]-x[0]*x[0]);});
+    DoubleFunction f2([b](const double *x){return b*(x[1]-x[0]*x[0]);}, 2);
     vector<DoubleFunction *> system = {&f1, &f2};
     GslHybridRootFinder solve(2);
     DVec init(2), x;

lib/Chi2Function.cpp

@@ -29,8 +29,7 @@ using namespace Latan;
  ******************************************************************************/
 // constructors ////////////////////////////////////////////////////////////////
 Chi2Function::Chi2Function(const XYStatData &data)
-: DoubleFunction(0)
+: data_(data)
-, data_(data)
 , buffer_(new Chi2FunctionBuffer)
 {
     resizeBuffer();
@@ -299,3 +298,22 @@ double Chi2Function::operator()(const double *arg) const
     
     return res;
 }
+
+// DoubleFunction factory //////////////////////////////////////////////////////
+DoubleFunction Chi2Function::makeFunction(const bool makeHardCopy) const
+{
+    DoubleFunction res;
+
+    if (makeHardCopy)
+    {
+        Chi2Function copy(*this);
+
+        res.setFunction([copy](const double *p){return copy(p);}, getNPar());
+    }
+    else
+    {
+        res.setFunction([this](const double *p){return (*this)(p);}, getNPar());
+    }
+
+    return res;
+}

lib/Chi2Function.hpp

@@ -32,7 +32,7 @@ BEGIN_LATAN_NAMESPACE
 // forward declaration of XYStatData
 class XYStatData;
 
-class Chi2Function: public DoubleFunction
+class Chi2Function: public DoubleFunctionFactory
 {
 private:
     struct Chi2FunctionBuffer
@@ -57,10 +57,9 @@ public:
             void  setModel(const std::vector<const DoubleModel *> &modelVector);
             void  requestInit(void) const;
     // function call
-    using DoubleFunction::operator();
+    double operator()(const double *arg) const;
-protected:
+    // factory
-    // function call
+    virtual DoubleFunction makeFunction(const bool makeHardCopy = true) const;
-    virtual double operator()(const double *arg) const;
 private:
     // access
     void resizeBuffer(void) const;

lib/CompiledFunction.cpp

@@ -28,12 +28,12 @@ using namespace Latan;
  *                   Compiled double function implementation                  *
  ******************************************************************************/
 // constructors ////////////////////////////////////////////////////////////////
-CompiledDoubleFunction::CompiledDoubleFunction(const unsigned nArg)
+CompiledDoubleFunction::CompiledDoubleFunction(const Index nArg)
-: DoubleFunction(nArg, nullptr)
+: nArg_(nArg)
 {}
 
-CompiledDoubleFunction::CompiledDoubleFunction(const unsigned nArg,
+CompiledDoubleFunction::CompiledDoubleFunction(const string &code,
-                                               const string &code)
+                                               const Index nArg)
 : CompiledDoubleFunction(nArg)
 {
     setCode(code);
@@ -60,7 +60,7 @@ void CompiledDoubleFunction::compile(void) const
     if (!*isCompiled_)
     {
         varAddress_->clear();
-        for (Index i = 0; i < getNArg(); ++i)
+        for (Index i = 0; i < nArg_; ++i)
         {
             varAddress_->push_back(context_->addVariable("x_" + strFrom(i)));
         }
@@ -75,7 +75,7 @@ double CompiledDoubleFunction::operator()(const double *arg) const
     double result;
     
     compile();
-    for (unsigned int i = 0; i < getNArg(); ++i)
+    for (unsigned int i = 0; i < nArg_; ++i)
     {
         context_->setVariable((*varAddress_)[i], arg[i]);
     }
@@ -102,3 +102,30 @@ ostream & Latan::operator<<(ostream &out, CompiledDoubleFunction &f)
     
     return out;
 }
+
+// DoubleFunction factory //////////////////////////////////////////////////////
+DoubleFunction CompiledDoubleFunction::makeFunction(const bool makeHardCopy)
+                                                    const
+{
+    DoubleFunction res;
+
+    if (makeHardCopy)
+    {
+        CompiledDoubleFunction copy(*this);
+
+        res.setFunction([copy](const double *p){return copy(p);}, nArg_);
+    }
+    else
+    {
+        res.setFunction([this](const double *p){return (*this)(p);}, nArg_);
+    }
+    
+    return res;
+}
+
+DoubleFunction Latan::compile(const string code, const Index nArg)
+{
+    CompiledDoubleFunction compiledFunc(code, nArg);
+
+    return compiledFunc.makeFunction();
+}

lib/CompiledFunction.hpp

@@ -29,27 +29,30 @@ BEGIN_LATAN_NAMESPACE
 /******************************************************************************
  *                      compiled double function class                        *
  ******************************************************************************/
-class CompiledDoubleFunction: public DoubleFunction
+class CompiledDoubleFunction: public DoubleFunctionFactory
 {
 public:
     // constructors
-    explicit CompiledDoubleFunction(const unsigned nArg);
+    explicit CompiledDoubleFunction(const Index nArg);
-    CompiledDoubleFunction(const unsigned nArg, const std::string &code);
+    CompiledDoubleFunction(const std::string &code, const Index nArg);
     // destructor
     virtual ~CompiledDoubleFunction(void) = default;
     // access
     std::string getCode(void);
     void        setCode(const std::string &code);
     // function call
-    using DoubleFunction::operator();
+    double operator()(const double *arg) const;
-    virtual double operator()(const double *arg) const;
     // IO
     friend std::ostream & operator<<(std::ostream &out,
                                      CompiledDoubleFunction &f);
+    // factory
+    virtual DoubleFunction makeFunction(const bool makeHardCopy = true) const;
 private:
+
     // compile
     void compile(void) const;
 private:
+    Index                                      nArg_;
     std::string                                code_;
     std::shared_ptr<MathInterpreter>           interpreter_;
     std::shared_ptr<RunContext>                context_;
@@ -59,6 +62,9 @@ private:
 
 std::ostream & operator<<(std::ostream &out, CompiledDoubleFunction &f);
 
+// DoubleFunction factory
+DoubleFunction compile(const std::string code, const Index nArg);
+
 END_LATAN_NAMESPACE
 
 #endif // Latan_CompiledFunction_hpp_

lib/CompiledModel.cpp

@@ -28,14 +28,13 @@ using namespace Latan;
  *                 CompiledDoubleModel implementation                         *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-CompiledDoubleModel::CompiledDoubleModel(const unsigned nArg,
+CompiledDoubleModel::CompiledDoubleModel(const Index nArg, const Index nPar)
-                                         const unsigned nPar)
+: nArg_(nArg)
-: DoubleModel(nArg, nPar, nullptr)
+, nPar_(nPar)
 {}
 
-CompiledDoubleModel::CompiledDoubleModel(const unsigned nArg,
+CompiledDoubleModel::CompiledDoubleModel(const string &code, const Index nArg,
-                                         const unsigned nPar,
+                                         const Index nPar)
-                                         const string &code)
 : CompiledDoubleModel(nArg, nPar)
 {
     setCode(code);
@@ -64,11 +63,11 @@ void CompiledDoubleModel::compile(void) const
     {
         varAddress_->clear();
         parAddress_->clear();
-        for (Index i = 0; i < getNArg(); ++i)
+        for (Index i = 0; i < nArg_; ++i)
         {
             varAddress_->push_back(context_->addVariable("x_" + strFrom(i)));
         }
-        for (Index j = 0; j < getNPar(); ++j)
+        for (Index j = 0; j < nPar_; ++j)
         {
             parAddress_->push_back(context_->addVariable("p_" + strFrom(j)));
         }
@@ -84,11 +83,11 @@ double CompiledDoubleModel::operator()(const double *arg,
     double result;
     
     compile();
-    for (unsigned int i = 0; i < getNArg(); ++i)
+    for (unsigned int i = 0; i < nArg_; ++i)
     {
         context_->setVariable((*varAddress_)[i], arg[i]);
     }
-    for (unsigned int j = 0; j < getNPar(); ++j)
+    for (unsigned int j = 0; j < nPar_; ++j)
     {
         context_->setVariable((*parAddress_)[j], par[j]);
     }
@@ -115,3 +114,32 @@ ostream & Latan::operator<<(std::ostream &out, CompiledDoubleModel &m)
     
     return out;
 }
+
+// DoubleModel factory /////////////////////////////////////////////////////////
+DoubleModel CompiledDoubleModel::makeModel(const bool makeHardCopy) const
+{
+    DoubleModel res;
+
+    if (makeHardCopy)
+    {
+        CompiledDoubleModel copy(*this);
+
+        res.setFunction([copy](const double *x, const double *p)
+                        {return copy(x, p);}, nArg_, nPar_);
+    }
+    else
+    {
+        res.setFunction([this](const double *x, const double *p)
+                        {return (*this)(x, p);}, nArg_, nPar_);
+    }
+
+    return res;
+}
+
+DoubleModel Latan::compile(const std::string &code, const Index nArg,
+                           const Index nPar)
+{
+    CompiledDoubleModel compiledModel(code, nArg, nPar);
+
+    return compiledModel.makeModel();
+}

lib/CompiledModel.hpp

@@ -29,28 +29,30 @@ BEGIN_LATAN_NAMESPACE
 /******************************************************************************
  *                     compiled double model class                            *
  ******************************************************************************/
-class CompiledDoubleModel: public DoubleModel
+class CompiledDoubleModel: public DoubleModelFactory
 {
 public:
     // constructor
-    explicit CompiledDoubleModel(const unsigned nArg, const unsigned nPar);
+    CompiledDoubleModel(const Index nArg, const Index nPar);
-    CompiledDoubleModel(const unsigned nArg, const unsigned nPar,
+    CompiledDoubleModel(const std::string &code, const Index nArg,
-                        const std::string &code);
+                        const Index nPar);
     // destructor
     virtual ~CompiledDoubleModel(void) = default;
     // access
     std::string getCode(void);
     void        setCode(const std::string &code);
     // function call
-    using DoubleModel::operator();
+    double operator()(const double *arg, const double *par) const;
-    virtual double operator()(const double *arg, const double *par) const;
     // IO
     friend std::ostream & operator<<(std::ostream &out,
                                      CompiledDoubleModel &f);
+    // factory
+    DoubleModel makeModel(const bool makeHardCopy = true) const;
 private:
     // compile
     void compile(void) const;
 private:
+    Index                                      nArg_, nPar_;
     std::string                                code_;
     std::shared_ptr<MathInterpreter>           interpreter_;
     std::shared_ptr<RunContext>                context_;
@@ -60,6 +62,10 @@ private:
 
 std::ostream & operator<<(std::ostream &out, CompiledDoubleModel &f);
 
+// DoubleModel factory
+DoubleModel compile(const std::string &code, const Index nArg,
+                    const Index nPar);
+
 END_LATAN_NAMESPACE
 
 #endif // Latan_CompiledModel_hpp_

lib/Derivative.cpp

@@ -32,15 +32,14 @@ using namespace Math;
 // constructor /////////////////////////////////////////////////////////////////
 Derivative::Derivative(const DoubleFunction &f, const Index dir,
                        const double step)
-: DoubleFunction(f.getNArg())
+: f_(f)
-, f_(&f)
 , dir_(dir)
 , step_(step)
 , buffer_(new DVec(f.getNArg()))
 {}
 
 Derivative::Derivative(const DoubleFunction &f, const Index dir,
-                       const Index order, const DVec point, const double step)
+                       const Index order, const DVec &point, const double step)
 : Derivative(f, dir, step)
 {
     setOrderAndPoint(order, point);
@@ -64,10 +63,10 @@ double Derivative::getStep(void) const
 
 void Derivative::setFunction(const DoubleFunction &f)
 {
-    f_ = &f;
+    f_ = f;
 }
 
-void Derivative::setOrderAndPoint(const Index order, const DVec point)
+void Derivative::setOrderAndPoint(const Index order, const DVec &point)
 {
     if (order >= point.size())
     {
@@ -132,20 +131,47 @@ void Derivative::makeCoefficients(void)
 // function call ///////////////////////////////////////////////////////////////
 double Derivative::operator()(const double *x) const
 {
-    ConstMap<DVec> xMap(x, (*f_).getNArg());
+    ConstMap<DVec> xMap(x, f_.getNArg());
     double         res = 0.;
     
     *buffer_ = xMap;
     FOR_VEC(point_, i)
     {
         (*buffer_)(dir_) = x[dir_] + point_(i)*step_;
-        res += coefficient_[i]*(*f_)(*buffer_);
+        res += coefficient_[i]*f_(*buffer_);
     }
     res /= pow(step_, order_);
     
     return res;
 }
 
+// function factory ////////////////////////////////////////////////////////////
+DoubleFunction Derivative::makeFunction(const bool makeHardCopy) const
+{
+    DoubleFunction res;
+
+    if (makeHardCopy)
+    {
+        Derivative copy(*this);
+
+        res.setFunction([copy](const double *x){return copy(x);}, f_.getNArg());
+    }
+    else
+    {
+        res.setFunction([this](const double *x){return (*this)(x);},
+                        f_.getNArg());
+    }
+
+    return res;
+}
+
+DoubleFunction Latan::derivative(const DoubleFunction &f, const Index dir,
+                                 const Index order, const DVec point,
+                                 const double step)
+{
+    return Derivative(f, dir, order, point, step).makeFunction();
+}
+
 /******************************************************************************
  *                     CentralDerivative implementation                       *
  ******************************************************************************/
@@ -188,3 +214,11 @@ void CentralDerivative::tuneStep(void)
     
     setStep(step);
 }
+
+// function factory ////////////////////////////////////////////////////////////
+DoubleFunction Latan::centralDerivative(const DoubleFunction &f,
+                                        const Index dir, const Index order,
+                                        const Index precOrder)
+{
+    return CentralDerivative(f, dir, order, precOrder).makeFunction();
+}

lib/Derivative.hpp

@@ -28,14 +28,14 @@ BEGIN_LATAN_NAMESPACE
 /******************************************************************************
  *                              Derivative                                    *
  ******************************************************************************/
-class Derivative: public DoubleFunction
+class Derivative: public DoubleFunctionFactory
 {
 public:
     static constexpr double defaultStep = 1.0e-2;
 public:
     // constructor
     Derivative(const DoubleFunction &f, const Index dir, const Index order,
-               const DVec point, const double step = defaultStep);
+               const DVec &point, const double step = defaultStep);
     // destructor
     virtual ~Derivative(void) = default;
     // access
@@ -43,11 +43,12 @@ public:
     Index  getOrder(void) const;
     double getStep(void) const;
     void   setFunction(const DoubleFunction &f);
-    void   setOrderAndPoint(const Index order, const DVec point);
+    void   setOrderAndPoint(const Index order, const DVec &point);
     void   setStep(const double step);
     // function call
-    using DoubleFunction::operator();
+    double operator()(const double *x) const;
-    virtual double operator()(const double *x) const;
+    // function factory
+    virtual DoubleFunction makeFunction(const bool makeHardCopy = true) const;
 protected:
     // constructor
     Derivative(const DoubleFunction &f, const Index dir,
@@ -55,14 +56,18 @@ protected:
 private:
     void makeCoefficients(void);
 private:
-    const DoubleFunction  *f_;
+    DoubleFunction        f_;
     Index                 dir_, order_;
     double                step_;
     DVec                  point_, coefficient_;
     std::shared_ptr<DVec> buffer_;
 };
 
-class CentralDerivative: private Derivative
+DoubleFunction derivative(const DoubleFunction &f, const Index dir,
+                          const Index order, const DVec point,
+                          const double step = Derivative::defaultStep);
+
+class CentralDerivative: public Derivative
 {
 public:
     static const Index defaultPrecOrder = 2;
@@ -89,6 +94,11 @@ private:
     Index precOrder_;
 };
 
+DoubleFunction centralDerivative(const DoubleFunction &f, const Index dir = 0,
+                                 const Index order = 1,
+                                 const Index precOrder =
+                                    CentralDerivative::defaultPrecOrder);
+
 END_LATAN_NAMESPACE
 
 #endif // Latan_Derivative_hpp_

lib/Function.cpp

@@ -27,7 +27,7 @@ using namespace Latan;
  *                        DoubleFunction implementation                       *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-DoubleFunction::DoubleFunction(const Index nArg, const vecFunc &f)
+DoubleFunction::DoubleFunction(const vecFunc &f, const Index nArg)
 : buffer_(new DVec)
 {
     setFunction(f, nArg);
@@ -100,6 +100,28 @@ double DoubleFunction::operator()(void) const
     return (*this)(nullptr);
 }
 
+// bind ////////////////////////////////////////////////////////////////////////
+DoubleFunction DoubleFunction::bind(const Index argIndex, const double val)
+{
+    Index            nArg = getNArg();
+    shared_ptr<DVec> buf(new DVec(nArg));
+    DoubleFunction   copy(*this), bindFunc;
+
+    auto func = [copy, buf, argIndex, val](const double *arg)
+    {
+         ConstMap<DVec> argMap(arg, buf->size());
+
+        *buf = argMap;
+        (*buf)(argIndex) = val;
+
+        return copy(*buf);
+    };
+
+    bindFunc.setFunction(func, nArg - 1);
+
+    return bindFunc;
+}
+
 // arithmetic operators ////////////////////////////////////////////////////////
 DoubleFunction DoubleFunction::operator-(void) const
 {

lib/Function.hpp

@@ -39,21 +39,22 @@ private:
     typedef std::function<double(const double *)> vecFunc;
 public:
     // constructor
-    explicit DoubleFunction(const Index nArg = 0,
+    explicit DoubleFunction(const vecFunc &f = nullptr, const Index nArg = 0);
-                            const vecFunc &f = nullptr);
     // destructor
     virtual ~DoubleFunction(void) = default;
     // access
     virtual Index getNArg(void) const;
     void    setFunction(const vecFunc &f, const Index nArg);
     // function call
-    virtual double operator()(const double *arg) const;
+    double operator()(const double *arg) const;
     double operator()(const DVec &arg) const;
     double operator()(const std::vector<double> &arg) const;
     double operator()(std::stack<double> &arg) const;
     double operator()(void) const;
     template <typename... Ts>
     double operator()(const double arg0, const Ts... args) const;
+    // bind
+    DoubleFunction bind(const Index argIndex, const double val);
     // arithmetic operators
     DoubleFunction   operator-(void) const;
     DoubleFunction & operator+=(const DoubleFunction &f);
@@ -180,6 +181,20 @@ DSample DoubleFunctionSample::operator()(const double arg0,
     return (*this)(arg);
 }
 
+/******************************************************************************
+ *                       DoubleFunctionFactory class                          *
+ ******************************************************************************/
+class DoubleFunctionFactory
+{
+public:
+    // constructor
+    DoubleFunctionFactory(void) = default;
+    // destructor
+    virtual ~DoubleFunctionFactory(void) = default;
+    // factory
+    virtual DoubleFunction makeFunction(const bool makeHardCopy) const = 0;
+};
+
 END_LATAN_NAMESPACE
 
 #endif // Latan_Function_hpp_

lib/Math.cpp

@@ -29,11 +29,11 @@ using namespace Latan;
 
 #define DEF_STD_FUNC_1ARG(name) \
 auto name##VecFunc = [](const double arg[1]){return (name)(arg[0]);};\
-DoubleFunction STDMATH_NAMESPACE::name(1, name##VecFunc);
+DoubleFunction STDMATH_NAMESPACE::name(name##VecFunc, 1);
 
 #define DEF_STD_FUNC_2ARG(name) \
 auto name##VecFunc = [](const double arg[2]){return (name)(arg[0], arg[1]);};\
-DoubleFunction STDMATH_NAMESPACE::name(2, name##VecFunc);
+DoubleFunction STDMATH_NAMESPACE::name(name##VecFunc, 2);
 
 // Trigonometric functions
 DEF_STD_FUNC_1ARG(cos)

lib/Model.cpp

@@ -29,10 +29,10 @@ using namespace Latan;
  *                           Model implementation                             *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-DoubleModel::DoubleModel(const Index nArg, const Index nPar, const vecFunc &f)
+DoubleModel::DoubleModel(const vecFunc &f, const Index nArg, const Index nPar)
 : size_(new ModelSize)
 {
-    setFunction(nArg, nPar, f);
+    setFunction(f, nArg, nPar);
 }
 
 // access //////////////////////////////////////////////////////////////////////
@@ -46,8 +46,8 @@ Index DoubleModel::getNPar(void) const
     return size_->nPar;
 }
 
-void DoubleModel::setFunction(const Index nArg, const Index nPar,
+void DoubleModel::setFunction(const vecFunc &f, const Index nArg,
-                              const vecFunc &f)
+                              const Index nPar)
 {
     size_->nArg = nArg;
     size_->nPar = nPar;
@@ -95,25 +95,35 @@ double DoubleModel::operator()(const double *data, const double *par) const
 // model bind //////////////////////////////////////////////////////////////////
 DoubleFunction DoubleModel::fixArg(const DVec &arg) const
 {
-    auto modelWithVec = [this](const DVec &x, const double *p)
+    DoubleModel copy(*this);
-                        {return (*this)(x.data(), p);};
+
+    auto modelWithVec = [copy](const DVec &x, const double *p)
+    {
+        return copy(x.data(), p);
+    };
     auto modelBind    = bind(modelWithVec, arg, _1);
 
-    return DoubleFunction(getNPar(), modelBind);
+    return DoubleFunction(modelBind, getNPar());
 }
 
 DoubleFunction DoubleModel::fixPar(const DVec &par) const
 {
-    auto modelWithVec = [this](const double *x, const DVec &p)
+    DoubleModel copy(*this);
-                              {return (*this)(x, p.data());};
+
+    auto modelWithVec = [copy](const double *x, const DVec &p)
+    {
+        return copy(x, p.data());
+    };
     auto modelBind    = bind(modelWithVec, _1, par);
 
-    return DoubleFunction(getNArg(), modelBind);
+    return DoubleFunction(modelBind, getNArg());
 }
 
 DoubleFunction DoubleModel::toFunction(void) const
 {
-    auto func = [this](const double *x){return (*this)(x, x + getNArg());};
+    DoubleModel copy(*this);
 
-    return DoubleFunction(getNArg() + getNPar(), func);
+    auto func = [copy](const double *x){return copy(x, x + copy.getNArg());};
+
+    return DoubleFunction(func, getNArg() + getNPar());
 }

lib/Model.hpp

@@ -38,20 +38,20 @@ private:
     struct ModelSize{Index nArg, nPar;};
 public:
     // constructor
-    DoubleModel(const Index nArg = 0, const Index nPar = 0,
+    DoubleModel(const vecFunc &f = nullptr, const Index nArg = 0,
-                const vecFunc &f = nullptr);
+                const Index nPar = 0);
     // destructor
     virtual ~DoubleModel(void) = default;
     // access
     virtual Index getNArg(void) const;
     virtual Index getNPar(void) const;
-    void  setFunction(const Index nArg = 0, const Index nPar = 0,
+    void  setFunction(const vecFunc &f, const Index nArg,
-                      const vecFunc &f = nullptr);
+                      const Index nPar);
     // function call
     double operator()(const DVec &data, const DVec &par) const;
     double operator()(const std::vector<double> &data,
                       const std::vector<double> &par) const;
-    virtual double operator()(const double *data, const double *par) const;
+    double operator()(const double *data, const double *par) const;
     // bind
     DoubleFunction fixArg(const DVec &arg) const;
     DoubleFunction fixPar(const DVec &par) const;
@@ -64,6 +64,20 @@ private:
     vecFunc                    f_;
 };
 
+/******************************************************************************
+ *                       base class for model factories                       *
+ ******************************************************************************/
+class DoubleModelFactory
+{
+public:
+    // constructor
+    DoubleModelFactory(void) = default;
+    // destructor
+    virtual ~DoubleModelFactory(void) = default;
+    // factory
+    virtual DoubleModel makeModel(const bool makeHardCopy) const = 0;
+};
+
 END_LATAN_NAMESPACE
 
 #endif // Latan_Model_hpp_

lib/TabFunction.cpp

@@ -27,10 +27,6 @@ using namespace Latan;
  *                       TabFunction implementation                           *
  ******************************************************************************/
 // constructors ////////////////////////////////////////////////////////////////
-TabFunction::TabFunction(void)
-: DoubleFunction(1)
-{}
-
 TabFunction::TabFunction(const DVec &x, const DVec &y)
 : TabFunction()
 {
@@ -85,3 +81,33 @@ double TabFunction::operator()(const double *arg) const
     
     return y_a + (x - x_a)*(y_b - y_a)/(x_b - x_a);
 }
+
+// DoubleFunction factory //////////////////////////////////////////////////////
+DoubleFunction TabFunction::makeFunction(const bool makeHardCopy) const
+{
+    DoubleFunction res;
+
+    if (makeHardCopy)
+    {
+        TabFunction copy(*this);
+
+        res.setFunction([copy](const double *x){return copy(x);}, 1);
+    }
+    else
+    {
+        res.setFunction([this](const double *x){return (*this)(x);}, 1);
+    }
+
+    return res;
+}
+
+DoubleFunction Latan::interpolate(const DVec &x, const DVec &y)
+{
+    return TabFunction(x,y).makeFunction();
+}
+
+DoubleFunction Latan::interpolate(const XYStatData &data, const Index i,
+                                  const Index j)
+{
+    return TabFunction(data, i, j).makeFunction();
+}

lib/TabFunction.hpp

@@ -31,11 +31,11 @@ BEGIN_LATAN_NAMESPACE
  *                      tabulated function: 1D only                           *
  ******************************************************************************/
 
-class TabFunction: public DoubleFunction
+class TabFunction: public DoubleFunctionFactory
 {
 public:
     // constructors
-    TabFunction(void);
+    TabFunction(void) = default;
     TabFunction(const DVec &x, const DVec &y);
     TabFunction(const XYStatData &data, const Index i = 0, const Index j = 0);
     // destructor
@@ -44,12 +44,17 @@ public:
     void setData(const DVec &x, const DVec &y);
     void setData(const XYStatData &data, const Index i = 0, const Index j = 0);
     // function call
-    using DoubleFunction::operator();
+    double operator()(const double *arg) const;
-    virtual double operator()(const double *arg) const;
+    // factory
+    virtual DoubleFunction makeFunction(const bool makeHardCopy = true) const;
 private:
     std::map<double, double> value_;
 };
 
+DoubleFunction interpolate(const DVec &x, const DVec &y);
+DoubleFunction interpolate(const XYStatData &data, const Index i = 0,
+                           const Index j = 0);
+
 END_LATAN_NAMESPACE
 
 #endif // Latan_TabFunction_hpp_

lib/XYStatData.cpp

@@ -249,10 +249,11 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
     minimizer.setInit(fullInit);
     
     // fit
-    FitResult result;
+    DoubleFunction chi2 = chi2_.makeFunction(false);
+    FitResult      result;
     
-    result        = minimizer(chi2_);
+    result        = minimizer(chi2);
-    result.chi2_  = chi2_(result);
+    result.chi2_  = chi2(result);
     result.nDof_  = chi2_.getNDof();
     result.model_.resize(modelVector.size());
     for (unsigned int j = 0; j < modelVector.size(); ++j)