Real_Function.h

#ifndef QNANO_NEW_REAL_FUNCTION_DEFINED_H
#define QNANO_NEW_REAL_FUNCTION_DEFINED_H

#include <iostream>
#include <fstream>
#include <cstdlib>

template <typename T>
class Function_Template{
public:
  virtual T fct(double x)=0;

  inline T operator()(double x){ 
    return fct(x);
  }

  virtual void plot(const std::string &filename, double xbegin, double xend, size_t Nintervals){
    std::ofstream ofs(filename.c_str());
    double h=(xend-xbegin)/(Nintervals);
    for(size_t i=0; i<Nintervals+1; i++){
      double x=xbegin+h*i;
      ofs<<x<<" "<<fct(x)<<std::endl;
    }
  }
};


typedef Function_Template<double> Real_Function;

class: public Real_Function{
public: double fct(double x){return 0.;}
} Real_Function_Zero;

class: public Real_Function{
public: double fct(double x){return x;}
} Real_Function_Identity;



class Multivar{
public:
  size_t sz;
  double *d;
  
  inline size_t size()const{
    return sz;
  }
  void resize(size_t newsz){
    sz=newsz;
    delete[] d;
    d=new double[sz];
  }
  void set_zero(){
    for(size_t i=0; i<sz; i++)d[i]=0;
  }
  double access_check(size_t i){
    if(i<0){std::cerr<<"Error accessing Multivar: i<0"<<std::endl; exit(1);}
    if(i>=sz){std::cerr<<"Error accessing Multivar: i>="<<sz<<std::endl; exit(1);}
    return d[i];
  }
  double operator()(size_t i)const{
    return d[i];
  }
  double &operator()(size_t i){
    return d[i];
  }
  double &operator[](size_t i){
    return d[i];
  }
  double operator[](size_t i)const{
    return d[i];
  }


  inline void add(const Multivar &other){
    for(size_t i=0; i<sz; i++)d[i]+=other.d[i];
  }
  inline void subtract(const Multivar &other){
    for(size_t i=0; i<sz; i++)d[i]-=other.d[i];
  }
  inline void add_scaled(double scale, const Multivar &other){
    for(size_t i=0; i<sz; i++)d[i]+=scale*other.d[i];
  }
  inline void scale(double scale){
    for(size_t i=0; i<sz; i++)d[i]*=scale;
  }

  void copy(const Multivar &other){
    resize(other.sz);
    for(size_t i=0; i<sz; i++)d[i]=other.d[i];
  }
  Multivar(const Multivar &other){ d=new double[0]; copy(other);}
  Multivar & operator=(const Multivar &other){ copy(other);return *this;}

  Multivar(size_t sz_): sz(sz_){
    d=new double[sz];
    set_zero();
  }
  Multivar(){ 
    d=new double[0];
  }
  ~Multivar(){
    delete[] d;
  }

};


template <typename T> 
class Multivar_Function_Template{
public:

  virtual T fct(const Multivar &v)=0;

  inline T operator()(const Multivar &v){
    return fct(v);
  }

};
typedef Multivar_Function_Template<double> Multivar_Real_Function;



template <typename T> 
class Function_Template_from_Multivar: public Function_Template<T>{
public:
  Multivar_Function_Template<T> *mvfct;
  Multivar mv_start;
  Multivar mv_dir;
  

  virtual T fct(double x){
    Multivar mv(mv_start);
    mv.add_scaled(x, mv_dir);
    return mvfct->fct(mv);
  }
  Function_Template_from_Multivar(Multivar_Function_Template<T>*  mvfct_, const Multivar &mv_start_, const  Multivar &mv_dir_) : mvfct(mvfct_), mv_start(mv_start_), mv_dir(mv_dir_){
  }
  
};
typedef Function_Template_from_Multivar<double> Real_Function_from_Multivar;



template <typename T> 
class Function_Template_along_Axis: public Function_Template<T>{
public:
  Multivar_Function_Template<T> *mvfct;
  size_t axis;
  Multivar mv_ref;
  
  virtual T fct(double x){
    mv_ref(axis)=x;
    return mvfct->fct(mv_ref);
  }
  Function_Template_along_Axis(Multivar_Function_Template<T>&  mvfct_, size_t axis_, const Multivar &mv_ref_): mvfct(&mvfct_), axis(axis_), mv_ref(mv_ref_){
  }
};
typedef Function_Template_along_Axis<double> Real_Function_along_Axis;

#endif