decomp.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 DECOMP_H
#define DECOMP_H

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

class DecompGraph {
  friend ostream& operator<<(ostream& out, const DecompGraph& g);

  public:
                class Edge;
                class Node {
                        friend class DecompGraph;
                        private:
                                int set_component(int comp);
                        public:
                                map<int, set<Edge>::iterator> adj;
                                int weight;
                                int component;

                                Node(int weight_=0)
                                : weight(weight_), component(-1)
                                { }
                };

                class Edge {
                        public:
                                map<int, Node>::iterator node1, node2;
                                mutable int weight;

                                Edge(int weight_=0)
                                : weight(weight_)
                                { }

                                Edge(map<int, Node>::iterator node1_, map<int, Node>::iterator node2_, int weight_=0)
                                : node1(node1_), node2(node2_), weight(weight_)
                                { }

//                              map<int, Node>::iterator neighbour(int index) const {
//                                      return (index==node1->first) ? node2 : node1;
//                              }

                                Node& neighbour(const Node& node) const {
                                        return (&node==&node1->second) ? node2->second : node1->second;
                                }

                                bool operator<(const Edge& e) const {
                                        return pair<int,int>(node1->first, node2->first)<pair<int,int>(e.node1->first, e.node2->first);
                                }
                };

                map<int, Node> nodes;
                set<Edge> edges;

                int nrcomp;

    DecompGraph()
                : nrcomp(0)
                { }

                DecompGraph(const SparsityInfo& si);

                virtual ~DecompGraph() { }

                map<int, Node>::iterator add_node(int index, int weight=0);

                set<Edge>::iterator add_edge(map<int, Node>::iterator node1, map<int, Node>::iterator node2, int weight=0);
                set<Edge>::iterator add_edge(int node1, int node2, int weight=0);

        int n() const { return nodes.size(); }

        int m() const { return edges.size(); }

                void compute_connected_components();

                void get_component_members(vector<list<int> >& members);

                void compute_partition(int nparts);
};

class SplitFunc: public Func {
        private:
                const Pointer<Func> orig;

                ivector indices;

    mutable dvector point;

                Func::CurvatureType curv_type;

        public:
    vector<bool> ignore;

                SplitFunc(const Pointer<Func>& orig_, const ivector& indices_, const dvector& point_, const vector<int>& ignore_=vector<int>(0))
                : Func(indices_.dim()), orig(orig_), indices(indices_), point(point_), curv_type(orig_->get_curvature()), ignore(indices_.dim(), false)
                { assert(orig!=NULL);
                        for (vector<int>::const_iterator it(ignore_.begin()); it!=ignore_.end(); it++) ignore[*it]=true;
                }

                SplitFunc(const Pointer<Func>& orig_, const ivector& indices_, const dvector& point_, const set<int>& ignore_, Pointer<SparsityInfo> si=NULL)
                : Func(indices_.dim()), orig(orig_), indices(indices_), point(point_), curv_type(orig_->get_curvature()), ignore(indices_.dim(), false)
                { assert(orig!=NULL);
                        for (set<int>::const_iterator it(ignore_.begin()); it!=ignore_.end(); ++it) ignore[*it]=true;
                        sparsity=si;
                }

                void set_point(const UserVector<double>& x) const {
                        for (int i=0; i<indices.size(); i++)
                                if (!ignore[i]) point[indices[i]]=x(i);
                }

                double eval(const UserVector<double>& x) const {
                        set_point(x);
                        return orig->eval(point);
                }

        
                void grad(UserVector<double>& g, const UserVector<double>& x) const {
                        set_point(x);
                        dvector biggrad(orig->dim());
                        orig->grad(biggrad, point);
                        g=biggrad(indices);
                        for (int i=0; i<indices.size(); i++)
                                if (ignore[i]) g[i]=0.;
                }

                using Func::grad;

                void HessMult(UserVector<double>& y, const UserVector<double>& x, const UserVector<double>& z) const {
                        set_point(z);
                        
                        dvector big_x(point);
                        for (int i=0; i<indices.size(); i++)
                                if (!ignore[i]) big_x.SetElement(indices[i], x(i));
                        
                        dvector bighess(orig->dim());
                        orig->HessMult(bighess, big_x, point);
                        y=bighess(indices);
                }

                using Func::HessMult;

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

                        void set_point(IntervalVector& bigx, const IntervalVector& x) const {
                                for (int i=0; i<indices.size(); i++)
                                        if (!ignore[i]) bigx[indices[i]]=x(i);
                        }

                        interval<double> eval(const IntervalVector& x) const {
                                IntervalVector bigx(point);
                                set_point(bigx, x);
                                return orig->eval(bigx);
                        }

                        void grad(IntervalVector& g, const IntervalVector& x) const {
                                IntervalVector bigx(point);
                                set_point(bigx, x);
                                IntervalVector biggrad(bigx.dim());
                                orig->grad(biggrad, bigx);
                                for (int i=0; i<indices.size(); i++)
                                        g[i] = ignore[i] ? interval<double>(0.) : biggrad(indices[i]);
                        }

                        int valgrad(interval<double>& val, IntervalVector& y, const IntervalVector& x) const {
                                IntervalVector bigx(point);
                                set_point(bigx, x);
                                IntervalVector biggrad(bigx.dim());
                                int ret=orig->valgrad(val, biggrad, bigx);
                                for (int i=0; i<indices.size(); i++)
                                        y[i] = ignore[i] ? interval<double>(0.) : biggrad(indices[i]);
                                return ret;
                        }

                using Func::valgrad;
#endif

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

                void print(ostream& out) const;
};


class SplittingScheme2 {
public:
        vector<list<pair<int, int> > > new_pos;

        vector<ivector> newblock;

        SplittingScheme2(const DecompGraph& graph);
};

class Decomposition {
        private:
                bool replace_if_quadratic;

                int avg_blocksize;

                Pointer<Param> param;

        public:
                Pointer<set<int> > I_var;

                Decomposition(Pointer<Param> param_=NULL)
                : param(param_), replace_if_quadratic(param_ ? param_->get_i("Decomposition replace quadratic", 1) : 1),
                  avg_blocksize(param_ ? param_->get_i("Decomposition average blocksize", 5) : 5)
                { }

        
        void set_sparsity_pattern(MinlpProblem& prob, const vector<vector<dvector> >& sample_set);

        
        Pointer<MinlpProblem> decompose(MinlpProblem& prob, vector<vector<dvector> >& sample_set);

        Pointer<MinlpProblem> decompose(MinlpProblem& prob, vector<Pointer<SplittingScheme2> >& ss);
        void decompose(SepQcFunc& f, int block_offset, const SepQcFunc& old_f, int k, const SplittingScheme2& ss, Pointer<dvector> primal);

        Pointer<SplittingScheme2> get_splittingscheme(MinlpProblem& prob, int blocknr);

        Pointer<SepQcFunc> decompose(Pointer<Func> old_f, const dvector& primal);
};

#endif

---