vectorfunction.h

Go to the documentation of this file.

#ifndef MATH_VECTOR_FUNCTION_H
#define MATH_VECTOR_FUNCTION_H

#include "function.h"
#include "vector.h"
#include "matrix.h"
#include <memory>
#include <vector>
#include <assert.h>

namespace Math {


class ScalarFieldDirectionalFunction : public RealFunction
{
public:
  ScalarFieldDirectionalFunction(ScalarFieldFunction& _f,const Vector& x,const Vector& n,bool ref=false);
  virtual std::string Label() const;
  virtual void PreEval(Real t);
  virtual Real Eval(Real t);
  virtual Real Deriv(Real t);
  virtual Real Deriv2(Real t);
  virtual const Vector& LastEvaluation() const { return tmp; }

private:
  ScalarFieldFunction* f;
  Vector x;
  Vector n;
  Vector tmp;
};

class ScalarFieldProjectionFunction : public RealFunction
{
public:
  ScalarFieldProjectionFunction(ScalarFieldFunction& f,const Vector& x,int i);
  virtual std::string Label() const;
  virtual void PreEval(Real t);
  virtual Real Eval(Real t);
  virtual Real Deriv(Real t);
  virtual Real Deriv2(Real t);

private:
  ScalarFieldFunction* f;
  Real xi0;
  int i;
  Vector tmp;
};

class LinearScalarFieldFunction : public ScalarFieldFunction
{
public:
  LinearScalarFieldFunction(const Vector& _a,Real _b,bool ref=false) {
    if(ref) {
      a.setRef(_a);
      b = _b;
    }
    else {
      a=_a;
      b=_b;
    }
  }

  virtual std::string Label() const { return "<a^t*x+b>"; }
  virtual Real Eval(const Vector& x) { return a.dot(x)+b; }
  virtual void Gradient(const Vector& x,Vector& grad) { grad = a; }
  virtual Real Gradient_i(const Vector& x,int i) { return a(i); }
  virtual void Hessian(const Vector& x,Matrix& H) { H.setZero(); }
  virtual Real Hessian_ij(const Vector& x,int i,int j) { return 0; }

  Vector a;
  Real b;
};

class QuadraticScalarFieldFunction : public ScalarFieldFunction 
{
 public:
  QuadraticScalarFieldFunction(const Matrix& _A,const Vector& _b,Real _c,bool ref=false) {
    if(ref) {
      A.setRef(_A);
      b.setRef(_b);
      c = _c;
    }
    else {
      A=_A;
      b=_b;
      c=_c;
    }
  }
  virtual std::string Label() const { return "[0.5*x'Ax+bx+c]"; }
  virtual void PreEval(const Vector& x) { A.mul(x,Ax); }
  virtual Real Eval(const Vector& x) { return Half*x.dot(Ax) + x.dot(b) + c; }
  virtual void Gradient(const Vector& x,Vector& grad) 
  { grad.add(Ax,b); }
  virtual Real Gradient_i(const Vector& x,int i)
  { return Ax(i) + b(i); }
  virtual Real DirectionalDeriv(const Vector& x,const Vector& h)
  { return Ax.dot(h) + b.dot(h); }
  virtual void Hessian(const Vector& x,Matrix& H) { H = A; }
  virtual Real Hessian_ij(const Vector& x,int i,int j) { return A(i,j); }
  virtual Real DirectionalDeriv2(const Vector& x,const Vector& h)
  { Real sum = 0;
    for(int i=0;i<A.m;i++) sum += h(i)*A.dotRow(i,h);
    return sum;
  }

  Matrix A;
  Vector b;
  Real c;
  //temp
  Vector Ax;
};

struct NormSquaredScalarFieldFunction : public ScalarFieldFunction
{
  virtual std::string Label() const { return "<|x|^2>"; }
  virtual Real Eval(const Vector& x) { return x.normSquared(); }
  virtual void Gradient(const Vector& x,Vector& grad) { grad.add(x,x); }
  virtual Real Gradient_i(const Vector& x,int i) { return Two*x(i); }
  virtual void Hessian(const Vector& x,Matrix& H) { H.setIdentity(); H.inplaceMul(Two); }
  virtual Real Hessian_ij(const Vector& x,int i,int j) { return Two*Delta(i,j); }
};

struct NormScalarFieldFunction : public ScalarFieldFunction
{
  NormScalarFieldFunction(Real _degree=Two) : degree(_degree) {}
  virtual std::string Label() const;
  virtual void PreEval(const Vector& x);
  virtual Real Eval(const Vector& x);
  virtual void Gradient(const Vector& x,Vector& grad);
  virtual Real Gradient_i(const Vector& x,int i);
  virtual void Hessian(const Vector& x,Matrix& H);
  virtual Real Hessian_ij(const Vector& x,int i,int j);

  Real degree;
  Real norm;
};

class MinimumScalarFieldFunction : public ScalarFieldFunction
{
public:
  virtual std::string Label() const { return "<min_i xi>"; }
  virtual void PreEval(const Vector& x);
  virtual Real Eval(const Vector& x);
  virtual void Gradient(const Vector& x,Vector& grad);
  virtual Real Gradient_i(const Vector& x,int i);
  virtual void Hessian(const Vector& x,Matrix& H) { H.setZero(); }
  virtual Real Hessian_ij(const Vector& x,Matrix& H) { return Zero; }

  Real xmin;
  int index;
};

class MaximumScalarFieldFunction : public ScalarFieldFunction
{
public:
  virtual std::string Label() const { return "<max_i xi>"; }
  virtual void PreEval(const Vector& x);
  virtual Real Eval(const Vector& x);
  virtual void Gradient(const Vector& x,Vector& grad);
  virtual Real Gradient_i(const Vector& x,int i);
  virtual void Hessian(const Vector& x,Matrix& H) { H.setZero(); }
  virtual Real Hessian_ij(const Vector& x,Matrix& H) { return Zero; }

  Real xmax;
  int index;
};

class ComposeScalarFieldFunction : public ScalarFieldFunction
{
public:
  ComposeScalarFieldFunction(std::shared_ptr<RealFunction> _f,std::shared_ptr<ScalarFieldFunction> _g) : f(_f), g(_g) {}
  virtual std::string Label() const;
  virtual void PreEval(const Vector& x);
  virtual Real Eval(const Vector& x) { return f->Eval(gx); }
  virtual void Gradient(const Vector& x,Vector& grad);
  virtual Real Gradient_i(const Vector& x,int i);
  virtual void Hessian(const Vector& x,Matrix& H);
  virtual Real Hessian_ij(const Vector& x,int i,int j); 
 
  std::shared_ptr<RealFunction> f;
  std::shared_ptr<ScalarFieldFunction> g;
  //temp
  Real gx;
};




class Vector3FieldFunction : public VectorFieldFunction
{
public:
  virtual ~Vector3FieldFunction() {}
  virtual void Curl(const Vector& x,Vector& curl)
  { 
    assert(x.n==3);
    curl.resize(3);
    curl(0)=Jacobian_ij(x,2,1)-Jacobian_ij(x,1,2);
    curl(1)=Jacobian_ij(x,0,2)-Jacobian_ij(x,2,0);
    curl(2)=Jacobian_ij(x,1,0)-Jacobian_ij(x,0,1);
  }
  virtual int NumDimensions() const { return 3; }
};

class LinearVectorFieldFunction : public VectorFieldFunction
{
 public:
  LinearVectorFieldFunction(const Matrix& _A,const Vector& _b,bool ref=false) {
    if(ref) {
      A.setRef(_A);
      b.setRef(_b);
    }
    else {
      A=_A;
      b=_b;
    }
  }
  virtual std::string Label() const { return "[A*x+b]"; }
  virtual int NumDimensions() const { return A.m; }
  virtual void Eval(const Vector& x,Vector& v) { A.mul(x,v); v+=b;  }
  virtual Real Eval_i(const Vector& x,int i) { return A.dotRow(i,x)+b(i); }
  virtual Real Jacobian_ij(const Vector& x,int i,int j) { return A(i,j); }
  virtual void Jacobian_i(const Vector& x,int i,Vector& Ji) { A.getRowCopy(i,Ji); }
  virtual void Jacobian_j(const Vector& x,int j,Vector& Jj) { A.getColCopy(j,Jj); }
  virtual void Jacobian(const Vector& x,Matrix& J) { J=A; }
  virtual void DirectionalDeriv(const Vector& x,const Vector& h,Vector& v) { A.mul(h,v); }
  virtual void Hessian_i(const Vector& x,int i,Matrix& Hi) { Hi.setZero(); }
  virtual Real Hessian_ijk(const Vector& x,int i,int j,int k) { return 0; }

  Matrix A;
  Vector b;
};

class Compose_SF_VF_Function : public ScalarFieldFunction
{
public:
  Compose_SF_VF_Function(std::shared_ptr<ScalarFieldFunction> _f,std::shared_ptr<VectorFieldFunction> _g) : f(_f), g(_g) {}
  virtual std::string Label() const;
  virtual void PreEval(const Vector& x);
  virtual Real Eval(const Vector& x) { return f->Eval(gx); }
  virtual void Gradient(const Vector& x,Vector& grad);
  virtual Real Gradient_i(const Vector& x,int i);
  virtual void Hessian(const Vector& x,Matrix& H);
  virtual Real Hessian_ij(const Vector& x,int i,int j);
  
  std::shared_ptr<ScalarFieldFunction> f;
  std::shared_ptr<VectorFieldFunction> g;
  //temp
  Vector gx;
  Vector gradf;
  Matrix Jg;
};

class Compose_VF_VF_Function : public ScalarFieldFunction
{
public:
  Compose_VF_VF_Function(std::shared_ptr<VectorFieldFunction> _f,std::shared_ptr<VectorFieldFunction> _g) : f(_f), g(_g) {}
  virtual std::string Label() const;
  virtual void PreEval(const Vector& x);
  virtual void Eval(const Vector& x,Vector& v) { f->Eval(gx,v); }
  virtual Real Eval_i(const Vector& x,int i) { return f->Eval_i(gx,i); }
  virtual void Jacobian(const Vector& x,Matrix& J);
  virtual void Jacobian_i(const Vector& x,int i,Vector& Ji);
  virtual void Jacobian_j(const Vector& x,int j,Vector& Jj);
 
  std::shared_ptr<VectorFieldFunction> f;
  std::shared_ptr<VectorFieldFunction> g;
  Vector gx;
  Matrix Jg,Jf;
};

class ComponentVectorFieldFunction : public VectorFieldFunction
{
public:
  virtual std::string Label() const;
  virtual std::string Label(int i) const;
  virtual int NumDimensions() const;
  virtual void PreEval(const Vector& x);
  virtual void Eval(const Vector& x, Vector& v);
  virtual Real Eval_i(const Vector& x,int i);
  virtual void Jacobian(const Vector& x,Matrix& J);
  virtual void Jacobian_i(const Vector& x,int i,Vector& Ji);
  virtual void DirectionalDeriv(const Vector& x,const Vector& h,Vector& v);
  virtual void Hessian_i(const Vector& x,int i,Matrix& Hi);

  std::vector<std::shared_ptr<ScalarFieldFunction> > functions;
};

class VectorFieldProjectionFunction : public ScalarFieldFunction
{
public:
  VectorFieldProjectionFunction(VectorFieldFunction& f,int i);
  virtual std::string Label() const;
  virtual void PreEval(const Vector& x) { f->PreEval(x); }
  virtual Real Eval(const Vector& x);
  virtual void Gradient(const Vector& x,Vector& grad);

private:
  VectorFieldFunction* f;
  int i;
};


class CompositeVectorFieldFunction : public VectorFieldFunction
{
public:
  CompositeVectorFieldFunction();
  CompositeVectorFieldFunction(std::shared_ptr<VectorFieldFunction> f1,std::shared_ptr<VectorFieldFunction>& f2);
  CompositeVectorFieldFunction(const std::vector<std::shared_ptr<VectorFieldFunction> >& fs);
  virtual std::string Label() const;
  virtual std::string Label(int i) const;
  virtual int NumDimensions() const;
  virtual void PreEval(const Vector& x);
  virtual void Eval(const Vector& x, Vector& v);
  virtual Real Eval_i(const Vector& x,int i);
  virtual void Jacobian(const Vector& x,Matrix& J);
  virtual void Jacobian_i(const Vector& x,int i,Vector& Ji);
  virtual void DirectionalDeriv(const Vector& x,const Vector& h,Vector& v);
  virtual void Hessian_i(const Vector& x,int i,Matrix& Hi);

  int GetFunction(int &i) const;

  std::vector<std::shared_ptr<VectorFieldFunction> > functions;
};

class IndexedVectorFieldFunction : public VectorFieldFunction
{
public:
  virtual std::string Label() const;
  virtual std::string Label(int i) const;
  virtual int NumDimensions() const;
  virtual void PreEval(const Vector& x);
  virtual void Eval(const Vector& x, Vector& v);
  virtual Real Eval_i(const Vector& x,int i);
  virtual void Jacobian(const Vector& x,Matrix& J);
  virtual void Jacobian_i(const Vector& x,int i,Vector& Ji);
  virtual void DirectionalDeriv(const Vector& x,const Vector& h,Vector& v);
  virtual void Hessian_i(const Vector& x,int i,Matrix& Hi);

  std::shared_ptr<VectorFieldFunction> function;
  std::vector<int> xindices,findices;

  //temp
  Vector vf,xsub;
  Matrix Jf;
};

class SliceVectorFieldFunction : public VectorFieldFunction
{
public:
  virtual std::string Label() const;
  virtual std::string Label(int i) const;
  virtual int NumDimensions() const;
  virtual void PreEval(const Vector& x);
  virtual void Eval(const Vector& x, Vector& v);
  virtual Real Eval_i(const Vector& x,int i);
  virtual void Jacobian(const Vector& x,Matrix& J);
  virtual void Jacobian_i(const Vector& x,int i,Vector& Ji);
  virtual void DirectionalDeriv(const Vector& x,const Vector& h,Vector& v);
  virtual void Hessian_i(const Vector& x,int i,Matrix& Hi);

  std::shared_ptr<VectorFieldFunction> function;
  Vector x0;
  std::vector<int> xindices;

  //temp
  Vector vf,xfull;
  Matrix Jf;
};


} //namespace Math

#endif