func.h

Go to the documentation of this file.

// Copyright (C) 2006 Ivo Nowak and Stefan Vigerske
// All Rights Reserved.
// This code is published under the Common Public License.
//
// Author: Stefan Vigerske

#ifndef FUNC_H
#define FUNC_H

#include "standard.h"
#include "usermatrix.h"

class Func;
class Decomposition;

class SparsityInfo {
        friend class Decomposition;
        private:
                void add_nonlinear(int index) {
                        nonlinear->insert(pair<int, NonlinearVariable>(index, NonlinearVariable()));
                }

                void add_quadratic(int index, double coeff_lin, double coeff_quad);

                void add_linear(int index, double coeff) {
                        linear->insert(pair<int, LinearVariable>(index, LinearVariable(coeff)));
                }

                void add_sparsity_pattern(int index1, int index2, double coeff);

        public:
                class LinearVariable {
                        public:
                                double coeff;
                                LinearVariable(double coeff_)
                                : coeff(coeff_)
                                { }

                                friend ostream& operator<<(ostream& out, const LinearVariable& var) {
                                        out << var.coeff;
                                        return out;
                                }
                };

                class NonlinearVariable {
                        public:
                };

                class QuadraticVariable {
                        public:
                                double coeff_lin, coeff_quad;

                                QuadraticVariable(double coeff_lin_, double coeff_quad_)
                                : coeff_lin(coeff_lin_), coeff_quad(coeff_quad_)
                                { }

                                friend ostream& operator<<(ostream& out, const QuadraticVariable& var) {
                                        out << var.coeff_lin << ',' << var.coeff_quad;
                                        return out;
                                }
                };

                class NonlinearConnection {
                        public:
                                double coeff;

                                NonlinearConnection(double coeff_)
                                : coeff(coeff_)
                                { }

                                friend ostream& operator<<(ostream& out, const NonlinearConnection& nc) {
                                        out << nc.coeff;
                                        return out;
                                }
                };

                Pointer<map<int, SparsityInfo::LinearVariable> > linear;
                Pointer<map<int, SparsityInfo::NonlinearVariable> > nonlinear;

                Pointer<map<int, SparsityInfo::QuadraticVariable> > quadratic;

                Pointer<map<pair<int, int>, SparsityInfo::NonlinearConnection> > sparsity_pattern;

                SparsityInfo(int init_level=0);

                SparsityInfo(const SparsityInfo& si);

                SparsityInfo(const SparsityInfo& si1, const SparsityInfo& si2);

                friend ostream& operator<<(ostream& out, const SparsityInfo& si);

                bool compute_sparsity_pattern(const Func& f, const vector<dvector>& sample_set);

//              void compute_sparsity_pattern(const UserMatrix& A);

                void add(const UserVector<double>& b);
                void add(const UserVector<double>& b, const ivector& block);
                void add(const UserMatrix& A);
                void add(const UserMatrix& A, const ivector& block);
                void add(const SparsityInfo& si);
                void add(const SparsityInfo& si, const ivector& block);

                int nr_var(bool linear_=true, bool nonlinear_=true) {
                        return (linear_ ? linear->size() : 0)+(nonlinear_ ? nonlinear->size() : 0);
                }
};

class VariableIterator_Type {
        protected:
                bool linear;
                bool nonlinear;
                bool quadratic;

        public:
                typedef enum { LINEAR=0, NONLINEAR=1, QUADRATIC=3, END=4 } VarType;

                VariableIterator_Type(bool linear_=true, bool nonlinear_=true, bool quadratic_=false)
                : linear(linear_), nonlinear(nonlinear_), quadratic(quadratic_)
                { }

                virtual ~VariableIterator_Type() { }

                virtual int operator()() const=0;

                virtual double coeff_lin() const=0;

                virtual double coeff_quad() const=0;

                virtual VarType type() const=0;

                virtual void operator++()=0;

                virtual operator bool() const=0;
};

class VariableIterator : public VariableIterator_Type {
        private:
                const SparsityInfo& sparsity;
                map<int, SparsityInfo::LinearVariable>::iterator it_linear;
                map<int, SparsityInfo::NonlinearVariable>::iterator it_nonlinear;
                map<int, SparsityInfo::QuadraticVariable>::iterator it_quadratic;
                VarType whereami;

                void init();
        public:
                VariableIterator(const Func& f_, bool linear_=true, bool nonlinear_=true, bool quadratic_=false);
                VariableIterator(const SparsityInfo& sparsity_, bool linear_=true, bool nonlinear_=true, bool quadratic_=false);

                virtual int operator()() const;

                virtual double coeff_lin() const;
                virtual double coeff_quad() const;

                virtual VarType type() const { return whereami; }

                virtual void operator++();

                virtual operator bool() const { return (whereami!=END); }
};

class Func {
                friend ostream& operator << (ostream& out, const Func& a)
                {  a.print(out); return out; };

                friend class VariableIterator;
                friend class MinusFunc;
                friend class SumFunc;

        protected:
                int dim_;

                Pointer<SparsityInfo> sparsity;

                virtual SparsityInfo& get_sparsity() { assert(sparsity); return *sparsity; }

        public:
                Pointer<ostream> out_func_p;

#define out_func if (out_func_p) (*out_func_p)

      Pointer<ostream> out_func_log_p;
#define out_func_log if (out_func_log_p) (*out_func_log_p)

      // methods:
      Func(int n=0, Pointer<ostream> out_func_p_=out_out_p, Pointer<ostream> out_func_log_p_=out_log_p)
      : dim_(n), out_func_p(out_func_p_), out_func_log_p(out_func_log_p_)
      { };

      virtual ~Func() { }

      int dim() const { return dim_; };

                        virtual bool sparsity_available() const { return sparsity; }

                        virtual const SparsityInfo& get_sparsity() const { assert(sparsity); return *sparsity; }

                        virtual bool compute_sparsity_pattern(const vector<dvector>& sample_set) { return get_sparsity().compute_sparsity_pattern(*this, sample_set); }

      virtual double eval(const UserVector<double>& x) const=0;

      virtual void grad(UserVector<double>& g, const UserVector<double>& x) const=0;

      virtual dvector grad(const dvector& x) const {
        dvector g(x.dim());
        grad(g, x);
        return g;
      }

      virtual void HessMult(UserVector<double>& y, const UserVector<double>& x, const UserVector<double>& z) const=0;

      virtual dvector HessMult(const dvector& x, const dvector& z) const {
        dvector y(x.size());
        HessMult(y, x, z);
        return y;
      }

      virtual int valgrad(double& val, UserVector<double>& g, const UserVector<double>& x) const {
        val=eval(x);
        grad(g,x);
        return 0;
      };

      virtual void add_grad(UserVector<double>& g, const UserVector<double>& x) const {
        Pointer<UserVector<double> > gr(g.getemptycopy());
        grad(*gr, x);
        g+=*gr;
      };

      virtual void add_valgrad(double& val, UserVector<double>& g, const UserVector<double>& x) const {
        val+=eval(x);
        add_grad(g,x);
      };

      virtual double min_eig_hess(const UserVector<double>& x, Param* param=NULL) const;

//      virtual dvector get_convex_shift(const UserVector<double>& lower, const UserVector<double>& upper, int lev, double offset) const;

#ifdef FILIB_AVAILABLE

                        virtual bool is_interval_compliant() const { return false; }

                        virtual interval<double> eval(const IntervalVector& x) const {
                                IntervalVector g(x.dim());
                                interval<double> val;
                                valgrad(val, g, x);
                                return val;
                        };

                        virtual void grad(IntervalVector& g, const IntervalVector& x) const {
                                interval<double> val;
                                valgrad(val, g, x);
                        };

                        virtual int valgrad(interval<double>& val, IntervalVector& g, const IntervalVector& x) const {
                                out_err << "Interval evaluation not implemented for this class." << endl;
                                exit(-1);
                                return 1;
                        }
#endif
                        typedef enum { CONVEX=1, CONCAVE=2, LINEAR=CONVEX | CONCAVE, UNKNOWN=4, INDEFINITE=12 } CurvatureType;
                        
                        friend ostream& operator<<(ostream& out, const CurvatureType& ct);

                        virtual void set_curvature(CurvatureType ct)=0;
                        virtual CurvatureType get_curvature() const=0;

                        static CurvatureType add_curvatures(double a1, CurvatureType ct1, double a2, CurvatureType ct2);
                        static CurvatureType mult_curvature(double a, CurvatureType ct);

      virtual void print(ostream &out) const
      {  out << "Func: dim=" << dim() << endl;  };
};

// -----------------------------------------------------------------

class HessMatrix : public UserMatrix {
  protected:
    Pointer<const UserVector<double> > x0;

    Pointer<const Func> f;

  public:
    HessMatrix(Pointer<const Func> f_, Pointer<const UserVector<double> > x0_)
    : UserMatrix(f_->dim()), f(f_), x0(x0_)
    { }

    HessMatrix(const Func& f_, const UserVector<double>& x0_)
    : UserMatrix(f_.dim()), f(&f_, false), x0(&x0_, false)
    { }

    void MultV(UserVector<double> &y, const UserVector<double> &x) const {
      f->HessMult(y, *x0, x);
    }

                using UserMatrix::MultV;

    void print(ostream &out) const {
      out << "HessMatrix for UserVector<double>: " << *x0 << "function: " << *f;
    }

};

// ----------------------------------------------------------------------

class MinusFunc : public Func {
  protected:
    Pointer<Func> f;

                SparsityInfo& get_sparsity() { return sparsity ? *sparsity : f->get_sparsity(); }

  public:
    MinusFunc(Pointer<Func> f_, Pointer<ostream> out_func_p_=out_out_p, Pointer<ostream> out_func_log_p_=out_log_p)
    : Func(f_ ? f_->dim() : 0, out_func_p_, out_func_log_p_), f(f_)
    { assert(f_ != NULL);
    }

                virtual bool sparsity_available() const { return sparsity || f->sparsity_available(); }

                virtual const SparsityInfo& get_sparsity() const { return sparsity ? *sparsity : ((const Func*)(Func*)f)->get_sparsity(); }

                virtual bool compute_sparsity_pattern(const vector<dvector>& sample_set);

    double eval(const UserVector<double>& x) const {
      return - f->eval(x);
    }

    void grad(UserVector<double>& g, const UserVector<double>& x) const {
      f->grad(g,x);
      g*=-1;
    }
                
                dvector grad(const dvector& x) const {
                        return -f->grad(x);
                }

    void HessMult(UserVector<double>& y, const UserVector<double>& x, const UserVector<double>& z) const {
      f->HessMult(y, x, z);
      y*=-1;
    }
                
                dvector HessMult(const dvector& x, const dvector& z) const {
                        return -f->HessMult(x,z);
                }

#ifdef FILIB_AVAILABLE
                bool is_interval_compliant() const { return f->is_interval_compliant(); }

                interval<double> eval(const IntervalVector& x) const {
                        return -f->eval(x);
                };

                void grad(IntervalVector& g, const IntervalVector& x) const {
                        f->grad(g,x);
                        g*=-1;
                };

                int valgrad(interval<double>& val, IntervalVector& g, const IntervalVector& x) const {
                        int ret=f->valgrad(val, g, x);
                        val*=-1;
                        g*=-1;
                        return ret;
                }

                using Func::valgrad;

#endif

                virtual void set_curvature(CurvatureType ct) {
                        switch(ct) {
                                case Func::CONVEX: f->set_curvature(Func::CONCAVE); break;
                                case Func::CONCAVE: f->set_curvature(Func::CONVEX); break;
                                default: f->set_curvature(ct); // UNKNOWN, UNDEFINED or LINEAR
                        }
                }

                virtual CurvatureType get_curvature() const {
                        Func::CurvatureType ct(f->get_curvature());
                        switch(ct) {
                                case Func::CONVEX: return Func::CONCAVE;
                                case Func::CONCAVE: return Func::CONVEX;
                        }
                        return ct;
                }

    void print(ostream& out) const {
      out << "MinusFunc: function: " << *f;
    }

};

// ----------------------------------------------------------------------

class SumFunc: public Func {
  private:
        Pointer<Func> f;
        Pointer<Func> g;
        double a,b;

                Func::CurvatureType curv_type;

                SparsityInfo& get_sparsity() {
                        assert(sparsity_available());
                        if (sparsity) return *sparsity;
                        return f->get_sparsity(); // this gives me only a const-member
                }

  public:
        SumFunc(Pointer<Func> f_, Pointer<Func> g_, double a_=1., double b_=1., Pointer<ostream> out_func_p_=out_out_p, Pointer<ostream> out_func_log_p_=out_log_p)
    : Func(f_ ? f_->dim() : (g_ ? g_->dim() : 0), out_func_p_, out_func_log_p_), f(f_), g(g_), a(a_), b(b_),
                  curv_type(add_curvatures(f_ ? a_ : 0., f_ ? f_->get_curvature() : Func::LINEAR, g_ ? b_ : 0., g_ ? g_->get_curvature() : Func::LINEAR))
    { assert(f || g);
                        if (!f) { f=g; a=b; g=NULL; } // move to front
                        if (f && g && f->sparsity_available() && g->sparsity_available())
                                sparsity=new SparsityInfo(((const Func*)(Func*)f)->get_sparsity(), ((const Func*)(Func*)g)->get_sparsity());
                };

                virtual bool sparsity_available() const {
                        if (sparsity) return true;
                        if (g) return false; // f and g defined, but no sparsity set
                        return f->sparsity_available();
                }

                virtual const SparsityInfo& get_sparsity() const {
                        assert(sparsity_available());
                        if (sparsity) return *sparsity;
                        return f->get_sparsity();
                }

                virtual bool compute_sparsity_pattern(const vector<dvector>& sample_set);

                double eval(const UserVector<double>& x) const {
                        return a * (f ? f->eval(x) : 0) + b * (g ? g->eval(x) : 0);
                }

                void grad(UserVector<double>& y, const UserVector<double>& x) const;
                
                using Func::grad;

                int valgrad(double& val, UserVector<double>& y, const UserVector<double>& x) const;
                
                void HessMult(UserVector<double>& y, const UserVector<double>& x, const UserVector<double>& z) const;
                
                using Func::HessMult;

#ifdef FILIB_AVAILABLE
                bool is_interval_compliant() const { return (f ? f->is_interval_compliant() : true) && (g ? g->is_interval_compliant() : true); }

                interval<double> eval(const IntervalVector& x) const {
                        return a*(f ? f->eval(x) : 0)+ b*(g ? g->eval(x) : 0);
                };

                void grad(IntervalVector& y, const IntervalVector& x) const;

                int valgrad(interval<double>& val, IntervalVector& y, const IntervalVector& x) const;
#endif

                Func::CurvatureType get_curvature() const { return curv_type; }
                void set_curvature(Func::CurvatureType ct) { curv_type=ct; }

                void print(ostream& out) const {
                        out << "SumFunc a*f(x)+b*g(x): a=" << a << " b=" << b;
                        if (f) out << endl << *f; else out << " f=NULL";
                        if (g) out << endl << *g; else out << " g=NULL" << endl;
                }
                
};


// ----------------------------------------------------------------------

class SepFunc : public Func {
   public:
      vector<ivector> block;

      SepFunc(int n, Pointer<ostream> out_func_p_=out_out_p, Pointer<ostream> out_func_log_p_=out_log_p)
      : Func(n, out_func_p_, out_func_log_p_), block(1)
      { block[0].resize(n);
        for (int i=0; i<n; i++) block[0][i]=i;
      };

      SepFunc(const vector<ivector>& block_, Pointer<ostream> out_func_p_=out_out_p, Pointer<ostream> out_func_log_p_=out_log_p)
      : Func(0, out_func_p_, out_func_log_p_), block(block_)
      { set_dim();
      };

                        SepFunc(const SepFunc& f)
                        : Func(f.dim(), f.out_func_p, f.out_func_log_p), block(f.block)
                        { }

      void set_dim() {
        dim_=0;
        for (int i=0; i<block.size(); i++) dim_+=block[i].size();
      }

      virtual void print(ostream &out) const {
        out << "SepFunc: dim=" << dim() << endl;
        for (int i=0; i<block.size(); i++)
          out << "block " << i << ": " << block[i];
      };
};

// ----------------------------------------------------------------------

class SepQcFunc: public SepFunc {
        friend class Decomposition;
        private:
                vector<Pointer<SparsityInfo> > sparsity_block;
        public:
        class VariableIterator : public VariableIterator_Type {
                private:
                        const SepQcFunc& f;
                        map<int, SparsityInfo::LinearVariable>::iterator it_linear;
                        map<int, SparsityInfo::NonlinearVariable>::iterator it_nonlinear;
                        map<int, SparsityInfo::QuadraticVariable>::iterator it_quadratic;
                        VarType whereami;
                        int blocknr;
                        void find_start();

                public:
                        VariableIterator(const SepQcFunc& f_, bool linear_=true, bool nonlinear_=true, bool quadratic_=false);

                        virtual int bnr() const { return blocknr; }
                        virtual int bindex() const;
                        virtual int operator()() const { return f.block[blocknr][bindex()]; }

                        virtual double coeff_lin() const;
                        virtual double coeff_quad() const;

                        virtual VarType type() const { return whereami; }

                        virtual void operator++();

                        virtual operator bool() const { return (whereami!=END); }
        };
        typedef enum { CONSTANT, LINEAR, QUADRATIC, NONQUAD } ftype;
        friend ostream& operator<<(ostream& out, const ftype& ft);

      vector<Pointer<UserMatrix> > A;

      vector<Pointer<UserVector<double> > > b;

      vector<Pointer<Func> > s;

                        vector<Func::CurvatureType> curv_type;

      double c;

      SepQcFunc(int n, Pointer<ostream> out_func_p_=out_out_p, Pointer<ostream> out_func_log_p_=out_log_p)
      : SepFunc(n, out_func_p_, out_func_log_p_), A(1), b(1), c(0), s(1), curv_type(1, Func::UNKNOWN), sparsity_block(1)
      { }

      SepQcFunc(const vector<ivector>& block_, Pointer<ostream> out_func_p_=out_out_p, Pointer<ostream> out_func_log_p_=out_log_p)
      : SepFunc(block_, out_func_p_, out_func_log_p_), A(block_.size()), b(block_.size()), c(0), s(block_.size()), curv_type(block_.size(), Func::UNKNOWN), sparsity_block(block_.size())
      { }

      SepQcFunc(Pointer<UserMatrix> A_, Pointer<UserVector<double> > b_=NULL, Pointer<Func> s_=NULL, double c_=0, Pointer<SparsityInfo> sparsity_=NULL, Pointer<ostream> out_func_p_=out_out_p, Pointer<ostream> out_func_log_p_=out_log_p)
      : SepFunc(A_ ? A_->dim() : b_ ? b_->dim() : s_->dim(), out_func_p_, out_func_log_p_), A(1, A_), b(1, b_), s(1, s_), c(c_),
                          curv_type(1, A_ ? Func::UNKNOWN : (s_ ? s_->get_curvature() : Func::LINEAR)), sparsity_block(1, sparsity_)
      { sparsity=sparsity_; }

      SepQcFunc(const SepQcFunc &f, bool minus_func);

                        SepQcFunc(const SepQcFunc& f, const UserVector<double>& point, const int degree=2);

                        virtual bool sparsity_available(int k) const {
                                if (sparsity_block[k]) return true;
                                return false;
                        }

                        virtual bool sparsity_available() const { return sparsity; }

                        virtual SparsityInfo& get_sparsity(int k) {
                                return *sparsity_block[k];
                        }

                        virtual const SparsityInfo& get_sparsity(int k) const {
                                return *sparsity_block[k];
                        }

                        virtual SparsityInfo& get_sparsity() { assert(sparsity); return *sparsity; }

                        using Func::get_sparsity;

                        virtual void set_sparsity(int k, Pointer<SparsityInfo> si) {
                                sparsity_block[k]=si;
                        }
                        virtual void set_sparsity();

                        bool compute_sparsity(int block_nr, const vector<dvector>& sample_set, bool replace_if_quadratic=false);

      virtual ftype type(int k) {
        if (s[k]) return NONQUAD;
        if (A[k]) return QUADRATIC;
                                if (b[k]) return LINEAR;
        return CONSTANT;
      }

                        void add_var(int index, int bnum);

      void addmult(const SepQcFunc& f, double a_=1., double b_=1.);

      double eval(const UserVector<double>& x) const;

      double eval(const UserVector<double>& x, int k) const {
        return (A[k] ? A[k]->xAx(x) : 0) + (b[k] ? 2* (*b[k]*x) : 0) + (s[k] ? s[k]->eval(x) : 0);
      }

      void grad(UserVector<double>& g, const UserVector<double>& x) const;

                        using Func::grad;

                        void grad(UserVector<double>& g, const UserVector<double>& x, int bnum) const;

                        dvector grad(const UserVector<double>& x, int bnum) const {
                                dvector g(x.dim()); grad(g,x,bnum); return g;
                        }

      int valgrad(double& val, UserVector<double>& g, const UserVector<double>& x) const;

      void HessMult(UserVector<double>& y, const UserVector<double>& x, const UserVector<double>& z) const;

      void HessMult(UserVector<double>& y, const UserVector<double>& x, const UserVector<double>& z, int k) const;

                        using Func::HessMult;

#ifdef FILIB_AVAILABLE
                        bool is_interval_compliant() const;

                        interval<double> eval(const IntervalVector& x) const;

                        interval<double> eval(const IntervalVector& x, int k) const;

                        void grad(IntervalVector& g, const IntervalVector& x) const;

                        void grad(IntervalVector& g, const IntervalVector& x, int k) const;

                        int valgrad(interval<double>& val, IntervalVector& y, const IntervalVector& x) const;
#endif

                        Func::CurvatureType get_curvature(int k) const { return curv_type[k]; }
                        Func::CurvatureType get_curvature() const;
                        void set_curvature(int k, Func::CurvatureType ct) { curv_type[k]=ct; }
                        void set_curvature(Func::CurvatureType ct);

      void print(ostream &out) const;

                        void print(ostream& out, vector<Pointer<char> > var_names) const;

};

#endif // FUNC_H

---