cubature.hpp

Go to the documentation of this file.

#ifndef BORDERBASIX_CUBATURE_HPP
#define BORDERBASIX_CUBATURE_HPP
#include <map>
#include <cstdlib>
#include <borderbasix/mdebug.hpp>
#include <borderbasix/mpoly.hpp>
#include <borderbasix/solversdp_csdp.h>
#include <borderbasix/solversdp_sdpa.h>
#include <borderbasix/mpoly/matrix_of.hpp>
#include <tensor-decomposition/Decomposer.hpp>
#define MMX_DEBUG_COUT

#define TMPL template<class T>

template<class T>
struct cubature {

    typedef Monom<double, dynamicexp<'x'> > Monomial;
    typedef Dlex<Monomial> Ordering;
    typedef MPoly<std::vector<Monomial>, Ordering > Polynomial;

    cubature();

    void set_input(const Polynomial &p);

    void obj_def(int i)        { m_init_p=i;}
  
    T&   obj_mat(int i, int j) { return m_P[m_s*i+j]; }

    void run(int k, list<Polynomial> list2);

    void set_threshold(double eps) {m_eps_reconstr=eps;}

    void set_sdpsolver(solversdp* slv) {m_slvsdp=slv;}

    void set_degree  (int deg) {m_degree=deg;}

    void set_numeq   (int g) {m_numeq=g;}

    void set_numneq   (int h) {m_numneq=h;}

    void set_verbose  (int v) {m_verbose=v;}


private:

    void init_rnd_obj(unsigned s);
    void init_idt_obj(unsigned s);


    Polynomial m_input;
    solversdp* m_slvsdp;
    T*         m_P;
    unsigned   m_s;
    int        m_init_p;
    double     m_eps_reconstr;
    int        m_verbose;


public:
    std::vector<T> weight;
    std::vector<T> point;
    unsigned       dim;
    int        m_degree;
    int        m_numeq;
    int        m_numneq;
    T&    P_ini(int i, int j);

};

TMPL
cubature<T>::cubature(): m_slvsdp(new solversdp_csdp),
    m_eps_reconstr(1e-5), m_verbose(0) {
}

TMPL
void cubature<T>::init_rnd_obj(unsigned s) {
    T* tmp=new T[s*s];
    m_P=new T[s*s];
    m_s=s;

    for(unsigned j=0;j<s*s;j++) {
        m_P[j]=(T)0.;
        tmp[j]=double(rand())/RAND_MAX;
        //std::cout<<"tmp"<<tmp[j]<<std::endl;
    }


    for(unsigned k=0;k<s;k++)
        for(unsigned l=0;l<s;l++)
            for(unsigned n=0;n<s;n++)
                m_P[k*s+l] += tmp[k*s+n]*tmp[l*s+n];

}

TMPL
void cubature<T>::init_idt_obj(unsigned s) {

    m_P=new T[s*s];
    m_s=s;

    for(unsigned k=0;k<s;k++)
        for(unsigned l=0;l<s;l++) {
            if(k==l)
                m_P[s*k+l]=(T)1.;
            else
                m_P[s*k+l]=(T)0.;
        }
}


TMPL
void cubature<T>::set_input(const Polynomial& p) {
    m_input = p;
}

TMPL
void cubature<T>::run(int k, list<Polynomial> list2) {

    std::map<Monomial,T> coeff;
    std::map<Monomial,int> mat;

    Monomial m;


    mdebug()<<"Initial moments:"<<m_input;
    for(Polynomial::const_iterator it= m_input.begin();it!= m_input.end();it++) {
        m =*it;
        m.SetCoeff(1);
        coeff[m]=it->GetCoeff();
    }


    int i = 0, j = 0, l = 0, i2,j2,contneq=0,l2=0,dim2i=0,dim2j=0;
    int c = 0;
    int dimbloq,u;
    //int d = Degree(this->m_input);
    this->dim = leading_var(this->m_input)+1;
    //mdebug()<<"d"<<d<<"dim"<<this->dim;
    mdebug()<<"nv:"<<this->dim;

    std::vector<Monomial> new_mon;

    Polynomial moments,mnq;

    Polynomial lm2,lm = matrix_of<Polynomial>::choose(this->dim,k);



    mdebug()<<"lm:"<<lm<<"lm.size()"<<lm.size();

    // Initialize matrix P with random coefficients and
    // the objective function <X, P^t P, X> = trace(P X P^t).
    if(this->m_init_p==0)
        this->init_rnd_obj(lm.size());
    else
        this->init_idt_obj(lm.size());


    std::map<Monomial,int>::const_iterator idx;

    mdebug()<<"after initialize matrix P";

    // Hankel matrix indexed by the monomials of lm
    for(Polynomial::const_iterator it= lm.begin();it!= lm.end();it++) {
        j=i;
        for(Polynomial::const_iterator jt= it; jt!= lm.end();jt++) {
            m=(*it)*(*jt);
            //std::cout<<"mon"<<m<<std::endl;
            // If degree <=d then coefficients in M0.
            if(Degree(m)<=m_degree) {
                m_slvsdp->coeff(0,i,j)=coeff[m];
            } else {

                idx=mat.find(m);
                if(idx!=mat.end()) {
                    l = idx->second+1;
                } else {
                    mat[m]=c;
                    new_mon.push_back(m);
                    c++;
                    l=c;
                }

                m_slvsdp->coeff(l,i,j)=(T)1;

                m_slvsdp->obj_coeff(l)+=obj_mat(i,j);

                if(i!=j) {
                    //                  //  m_slvsdp->coeff(l,j,i)=(T)1;
                    m_slvsdp->obj_coeff(l)+=obj_mat(j,i);
                }
            }

            mdebug()<<"l"<<l<<",i:"<<i<<",j:"<<j<<",value"<<m_slvsdp->coeff(l,i,j);
            j++;
        }
        i++;
    }

    // Dimension of the Hankel matrix
    m_slvsdp->dim()=lm.size();


    //Localizing Matrices
    dim2i=lm.size();
    dim2j=lm.size();
    if (m_numneq !=0)
      for(list<Polynomial>::const_iterator itneq=list2.begin();itneq!=list2.end();itneq++)
      {
          contneq++;
          std::cout<<k<<Degree(*itneq)<<std::endl;
          dimbloq=(2*k-Degree(*itneq))/2;

          //dimbloq=u-floor(Degree(*itneq)/2);
          std::cout<<"dimbloq"<<dimbloq<<std::endl;
          lm2=matrix_of<Polynomial>::choose(this->dim,dimbloq);
          //std::cout<<"lm2"<<lm2<<std::endl;
          m_slvsdp->dim(contneq)=lm2.size();
          i2=0;
          j2=0;
          for(Polynomial::const_iterator itneq2= lm2.begin();itneq2!= lm2.end();itneq2++) {
              j2=i2;
              for(Polynomial::const_iterator jtneq2= itneq2; jtneq2!= lm2.end();jtneq2++) {
                  mnq=(*itneq)*(*itneq2)*(*jtneq2);
                  for(Polynomial::const_iterator itmnq= mnq.begin(); itmnq!= mnq.end();itmnq++) {

                      if(Degree(*itmnq)<=m_degree) {
                          m_slvsdp->coeff(0,contneq,i2,j2)+=(itmnq->GetCoeff())*coeff[*itmnq];
                      } else {


                          l2 = mat[*itmnq]+1;
                          m_slvsdp->coeff(l2,contneq,i2,j2)=itmnq->GetCoeff();


                      }


                  }
                  mdebug()<<"l"<<l2<<"bloque"<<contneq<<",i:"<<i2<<",j:"<<j2<<",value"<<m_slvsdp->coeff(l2,contneq,i2,j2);
                  j2++;
              }
              i2++;

          }
          std::cout<<"dim2i"<<dim2i<<std::endl;
          std::cout<<"lm2 size"<<lm2.size()<<std::endl;
          dim2i+=lm2.size();
          dim2j+=lm2.size();
      }



    mdebug()<<"After hankel matrices construction";

    ;

    // Solve min <X,P^t P> s.t. X = M0 + y1 M1 + ...
    m_slvsdp->run();


    // Construct the extended moment sequence from the solution
    // of the SDP problem.
    moments=m_input;


    for(unsigned i=0;i<m_slvsdp->obj_size();i++) {
        new_mon[i].SetCoeff(m_slvsdp->output[i]);
        moments+=new_mon[i];

    }


    mdebug()<<"Extended moments:"<<moments;

    mmx::Decomposer<Polynomial> decomp(moments,this->m_eps_reconstr,this->m_verbose,-1,-1,this->m_s+1,0);
    decomp.main_algo();


    // Get the weights
    weight.resize(decomp.lambdas.size());
    for(unsigned i=0;i<decomp.lambdas.size();i++)
        weight[i] = decomp.lambdasopt[i].real();


    // Get the points

    point.resize(decomp.lambdasopt.size()*this->dim);
    for(unsigned i=0;i<decomp.lambdasopt.size();i++) {
        for(unsigned j=0;j<decomp.linear_forms.nbrow();j++)
            point[i*this->dim+j]=decomp.linear_forms(j,i).real();


    }

    m_slvsdp->clear();

}

#undef TMPL
#endif// BORDERBASIX_CUBATURE_HPP