igl::mosek Namespace Reference

Classes
struct	MosekData

Functions
template<typename DerivedV , typename DerivedEle , typename Derivedb , typename Derivedbc , typename DerivedW >
IGL_INLINE bool	bbw (const Eigen::PlainObjectBase< DerivedV > &V, const Eigen::PlainObjectBase< DerivedEle > &Ele, const Eigen::PlainObjectBase< Derivedb > &b, const Eigen::PlainObjectBase< Derivedbc > &bc, igl::BBWData &data, igl::mosek::MosekData &mosek_data, Eigen::PlainObjectBase< DerivedW > &W)
IGL_INLINE MSKrescodee	mosek_guarded (const MSKrescodee r)
IGL_INLINE bool	mosek_linprog (const Eigen::VectorXd &c, const Eigen::SparseMatrix< double > &A, const Eigen::VectorXd &lc, const Eigen::VectorXd &uc, const Eigen::VectorXd &lx, const Eigen::VectorXd &ux, Eigen::VectorXd &x)
IGL_INLINE bool	mosek_linprog (const Eigen::VectorXd &c, const Eigen::SparseMatrix< double > &A, const Eigen::VectorXd &lc, const Eigen::VectorXd &uc, const Eigen::VectorXd &lx, const Eigen::VectorXd &ux, const MSKenv_t &env, Eigen::VectorXd &x)
template<typename Index , typename Scalar >
IGL_INLINE bool	mosek_quadprog (const Index n, std::vector< Index > &Qi, std::vector< Index > &Qj, std::vector< Scalar > &Qv, const std::vector< Scalar > &c, const Scalar cf, const Index m, std::vector< Scalar > &Av, std::vector< Index > &Ari, const std::vector< Index > &Acp, const std::vector< Scalar > &lc, const std::vector< Scalar > &uc, const std::vector< Scalar > &lx, const std::vector< Scalar > &ux, MosekData &mosek_data, std::vector< Scalar > &x)
IGL_INLINE bool	mosek_quadprog (const Eigen::SparseMatrix< double > &Q, const Eigen::VectorXd &c, const double cf, const Eigen::SparseMatrix< double > &A, const Eigen::VectorXd &lc, const Eigen::VectorXd &uc, const Eigen::VectorXd &lx, const Eigen::VectorXd &ux, MosekData &mosek_data, Eigen::VectorXd &x)

Function Documentation

bbw()

template<typename DerivedV , typename DerivedEle , typename Derivedb , typename Derivedbc , typename DerivedW >

IGL_INLINE bool igl::mosek::bbw	(	const Eigen::PlainObjectBase< DerivedV > &	V,
		const Eigen::PlainObjectBase< DerivedEle > &	Ele,
		const Eigen::PlainObjectBase< Derivedb > &	b,
		const Eigen::PlainObjectBase< Derivedbc > &	bc,
		igl::BBWData &	data,
		igl::mosek::MosekData &	mosek_data,
		Eigen::PlainObjectBase< DerivedW > &	W
	)

{
  using namespace std;
  using namespace Eigen;
  assert(!data.partition_unity && "partition_unity not implemented yet");
  // number of domain vertices
  int n = V.rows();
  // number of handles
  int m = bc.cols();
  // Build biharmonic operator
  Eigen::SparseMatrix<typename DerivedV::Scalar> Q;
  harmonic(V,Ele,2,Q);
  W.derived().resize(n,m);
  // No linear terms
  VectorXd c = VectorXd::Zero(n);
  // No linear constraints
  SparseMatrix<typename DerivedW::Scalar> A(0,n);
  VectorXd uc(0,1),lc(0,1);
  // Upper and lower box constraints (Constant bounds)
  VectorXd ux = VectorXd::Ones(n);
  VectorXd lx = VectorXd::Zero(n);
  // Loop over handles
  for(int i = 0;i<m;i++)
  {
    if(data.verbosity >= 1)
    {
      cout<<"BBW: Computing weight for handle "<<i+1<<" out of "<<m<<
        "."<<endl;
    }
    VectorXd bci = bc.col(i);
    VectorXd Wi;
    // impose boundary conditions via bounds
    slice_into(bci,b,ux);
    slice_into(bci,b,lx);
    bool r = mosek_quadprog(Q,c,0,A,lc,uc,lx,ux,mosek_data,Wi);
    if(!r)
    {
      return false;
    }
    W.col(i) = Wi;
  }
#ifndef NDEBUG
    const double min_rowsum = W.rowwise().sum().array().abs().minCoeff();
    if(min_rowsum < 0.1)
    {
      cerr<<"bbw.cpp: Warning, minimum row sum is very low. Consider more "
        "active set iterations or enforcing partition of unity."<<endl;
    }
#endif
 
  return true;
}

References igl::harmonic(), mosek_quadprog(), Eigen::PlainObjectBase< Derived >::resize(), Eigen::PlainObjectBase< Derived >::rows(), and igl::slice_into().

Here is the call graph for this function:

mosek_guarded()

IGL_INLINE MSKrescodee igl::mosek::mosek_guarded

(

const MSKrescodee

r

)

{
  using namespace std;
  if(r != MSK_RES_OK)
  {
    /* In case of an error print error code and description. */      
    char symname[MSK_MAX_STR_LEN];
    char desc[MSK_MAX_STR_LEN];
    MSK_getcodedesc(r,symname,desc);
    cerr<<"MOSEK ERROR ("<<r<<"): "<<symname<<" - '"<<desc<<"'"<<endl;
  }
  return r;
}

Referenced by mosek_linprog(), mosek_linprog(), and mosek_quadprog().

Here is the caller graph for this function:

mosek_linprog() [1/2]

IGL_INLINE bool igl::mosek::mosek_linprog	(	const Eigen::VectorXd &	c,
		const Eigen::SparseMatrix< double > &	A,
		const Eigen::VectorXd &	lc,
		const Eigen::VectorXd &	uc,
		const Eigen::VectorXd &	lx,
		const Eigen::VectorXd &	ux,
		const MSKenv_t &	env,
		Eigen::VectorXd &	x
	)

{
  // following http://docs.mosek.com/7.1/capi/Linear_optimization.html
  using namespace std;
  // number of constraints
  const int m = A.rows();
  // number of variables
  const int n = A.cols();
 
 
  vector<double> vAv;
  vector<int> vAri,vAcp;
  int nr;
  harwell_boeing(A,nr,vAv,vAri,vAcp);
 
  MSKtask_t task;
  // Create the optimization task
  mosek_guarded(MSK_maketask(env,m,n,&task));
  // no threads
  mosek_guarded(MSK_putintparam(task,MSK_IPAR_NUM_THREADS,1));
  if(m>0)
  {
    // Append 'm' empty constraints, the constrainst will initially have no
    // bounds
    mosek_guarded(MSK_appendcons(task,m));
  }
  mosek_guarded(MSK_appendvars(task,n));
 
  
  const auto & key = [](const double lxj, const double uxj) ->
    MSKboundkeye
  {
    MSKboundkeye k = MSK_BK_FR;
    if(isfinite(lxj) && isfinite(uxj))
    {
      if(lxj == uxj)
      {
        k = MSK_BK_FX;
      }else{
        k = MSK_BK_RA;
      }
    }else if(isfinite(lxj))
    {
      k = MSK_BK_LO;
    }else if(isfinite(uxj))
    {
      k = MSK_BK_UP;
    }
    return k;
  };
 
  // loop over variables
  for(int j = 0;j<n;j++)
  {
    if(c.size() > 0)
    {
      // Set linear term c_j in the objective
      mosek_guarded(MSK_putcj(task,j,c(j)));
    }
 
    // Set constant bounds on variable j
    const double lxj = lx.size()>0?lx[j]:-numeric_limits<double>::infinity();
    const double uxj = ux.size()>0?ux[j]: numeric_limits<double>::infinity();
    mosek_guarded(MSK_putvarbound(task,j,key(lxj,uxj),lxj,uxj));
 
    if(m>0)
    {
      // Input column j of A
      mosek_guarded(
        MSK_putacol(
          task,
          j,
          vAcp[j+1]-vAcp[j],
          &vAri[vAcp[j]],
          &vAv[vAcp[j]])
        );
    }
  }
  // loop over constraints
  for(int i = 0;i<m;i++)
  {
    // Set constraint bounds for row i
    const double lci = lc.size()>0?lc[i]:-numeric_limits<double>::infinity();
    const double uci = uc.size()>0?uc[i]: numeric_limits<double>::infinity();
    mosek_guarded(MSK_putconbound(task,i,key(lci,uci),lci,uci));
  }
 
  // Now the optimizer has been prepared
  MSKrescodee trmcode;
  // run the optimizer
  mosek_guarded(MSK_optimizetrm(task,&trmcode));
  // Get status
  MSKsolstae solsta;
  MSK_getsolsta (task,MSK_SOL_ITR,&solsta);
  bool success = false;
  switch(solsta)
  {
    case MSK_SOL_STA_OPTIMAL:   
    case MSK_SOL_STA_NEAR_OPTIMAL:
      x.resize(n);
      /* Request the basic solution. */ 
      MSK_getxx(task,MSK_SOL_BAS,x.data()); 
      success = true;
      break;
    case MSK_SOL_STA_DUAL_INFEAS_CER:
    case MSK_SOL_STA_PRIM_INFEAS_CER:
    case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
    case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:  
      //printf("Primal or dual infeasibility certificate found.\n");
      break;
    case MSK_SOL_STA_UNKNOWN:
      //printf("The status of the solution could not be determined.\n");
      break;
    default:
      //printf("Other solution status.");
      break;
  }
  MSK_deletetask(&task);
  return success;
}

References Eigen::SparseMatrix< _Scalar, _Options, _StorageIndex >::cols(), igl::harwell_boeing(), mosek_guarded(), and Eigen::SparseMatrix< _Scalar, _Options, _StorageIndex >::rows().

Here is the call graph for this function:

mosek_linprog() [2/2]

IGL_INLINE bool igl::mosek::mosek_linprog	(	const Eigen::VectorXd &	c,
		const Eigen::SparseMatrix< double > &	A,
		const Eigen::VectorXd &	lc,
		const Eigen::VectorXd &	uc,
		const Eigen::VectorXd &	lx,
		const Eigen::VectorXd &	ux,
		Eigen::VectorXd &	x
	)

{
  // variables for mosek task, env and result code
  MSKenv_t env;
  // Create the MOSEK environment
  mosek_guarded(MSK_makeenv(&env,NULL));
  // initialize mosek environment
#if MSK_VERSION_MAJOR <= 7
  mosek_guarded(MSK_initenv(env));
#endif
  const bool ret = mosek_linprog(c,A,lc,uc,lx,ux,env,x);
  MSK_deleteenv(&env);
  return ret;
}

References mosek_guarded(), and mosek_linprog().

Referenced by mosek_linprog().

Here is the call graph for this function:

Here is the caller graph for this function:

mosek_quadprog() [1/2]

IGL_INLINE bool igl::mosek::mosek_quadprog	(	const Eigen::SparseMatrix< double > &	Q,
		const Eigen::VectorXd &	c,
		const double	cf,
		const Eigen::SparseMatrix< double > &	A,
		const Eigen::VectorXd &	lc,
		const Eigen::VectorXd &	uc,
		const Eigen::VectorXd &	lx,
		const Eigen::VectorXd &	ux,
		MosekData &	mosek_data,
		Eigen::VectorXd &	x
	)

{
  using namespace Eigen;
  using namespace std;
 
  typedef int Index;
  typedef double Scalar;
  // Q should be square
  assert(Q.rows() == Q.cols());
  // Q should be symmetric
#ifdef EIGEN_HAS_A_BUG_AND_FAILS_TO_LET_ME_COMPUTE_Q_MINUS_Q_TRANSPOSE
  assert( (Q-Q.transpose()).sum() < FLOAT_EPS);
#endif
  // Only keep lower triangular part of Q
  SparseMatrix<Scalar> QL;
  //QL = Q.template triangularView<Lower>();
  QL = Q.triangularView<Lower>();
  VectorXi Qi,Qj;
  VectorXd Qv;
  find(QL,Qi,Qj,Qv);
  vector<Index> vQi = matrix_to_list(Qi);
  vector<Index> vQj = matrix_to_list(Qj);
  vector<Scalar> vQv = matrix_to_list(Qv);
 
  // Convert linear term
  vector<Scalar> vc = matrix_to_list(c);
 
  assert(lc.size() == A.rows());
  assert(uc.size() == A.rows());
  // Convert A to harwell boeing format
  vector<Scalar> vAv;
  vector<Index> vAr,vAc;
  Index nr;
  harwell_boeing(A,nr,vAv,vAr,vAc);
 
  assert(lx.size() == Q.rows());
  assert(ux.size() == Q.rows());
  vector<Scalar> vlc = matrix_to_list(lc);
  vector<Scalar> vuc = matrix_to_list(uc);
  vector<Scalar> vlx = matrix_to_list(lx);
  vector<Scalar> vux = matrix_to_list(ux);
 
  vector<Scalar> vx;
  bool ret = mosek_quadprog<Index,Scalar>(
    Q.rows(),vQi,vQj,vQv,
    vc,
    cf,
    nr,
    vAv, vAr, vAc,
    vlc,vuc,
    vlx,vux,
    mosek_data,
    vx);
  list_to_matrix(vx,x);
  return ret;
}

References Eigen::SparseMatrix< _Scalar, _Options, _StorageIndex >::cols(), igl::find(), igl::FLOAT_EPS, igl::harwell_boeing(), igl::list_to_matrix(), igl::matrix_to_list(), Eigen::SparseMatrix< _Scalar, _Options, _StorageIndex >::rows(), Eigen::SparseMatrixBase< Derived >::transpose(), and Eigen::SparseMatrixBase< Derived >::triangularView().

Here is the call graph for this function:

mosek_quadprog() [2/2]

template<typename Index , typename Scalar >

IGL_INLINE bool igl::mosek::mosek_quadprog	(	const Index	n,
		std::vector< Index > &	Qi,
		std::vector< Index > &	Qj,
		std::vector< Scalar > &	Qv,
		const std::vector< Scalar > &	c,
		const Scalar	cf,
		const Index	m,
		std::vector< Scalar > &	Av,
		std::vector< Index > &	Ari,
		const std::vector< Index > &	Acp,
		const std::vector< Scalar > &	lc,
		const std::vector< Scalar > &	uc,
		const std::vector< Scalar > &	lx,
		const std::vector< Scalar > &	ux,
		MosekData &	mosek_data,
		std::vector< Scalar > &	x
	)

• Directs the log stream to the 'printstr' function. *‍/ / Little function mosek needs in order to know how to print to std out

{
  // I J V vectors of Q should all be same length
  assert(Qv.size() == Qi.size());
  assert(Qv.size() == Qj.size());
  // number of columns in linear constraint matrix must be ≤ number of
  // variables
  assert( (int)Acp.size() == (n+1));
  // linear bound vectors must be size of number of constraints or empty
  assert( ((int)lc.size() == m) || ((int)lc.size() == 0));
  assert( ((int)uc.size() == m) || ((int)uc.size() == 0));
  // constant bound vectors must be size of number of variables or empty
  assert( ((int)lx.size() == n) || ((int)lx.size() == 0));
  assert( ((int)ux.size() == n) || ((int)ux.size() == 0));
 
  // allocate space for solution in x
  x.resize(n);
 
  // variables for mosek task, env and result code
  MSKenv_t env;
  MSKtask_t task;
 
  // Create the MOSEK environment
#if MSK_VERSION_MAJOR >= 7
  mosek_guarded(MSK_makeenv(&env,NULL));
#elif MSK_VERSION_MAJOR == 6
  mosek_guarded(MSK_makeenv(&env,NULL,NULL,NULL,NULL));
#endif
  //const auto & printstr = [](void *handle, char str[])
  //{
  //  printf("%s",str);
  //}
  //mosek_guarded(MSK_linkfunctoenvstream(env,MSK_STREAM_LOG,NULL,printstr));
  // initialize mosek environment
#if MSK_VERSION_MAJOR <= 7
  mosek_guarded(MSK_initenv(env));
#endif
  // Create the optimization task
  mosek_guarded(MSK_maketask(env,m,n,&task));
  verbose("Creating task with %ld linear constraints and %ld variables...\n",m,n);
  //mosek_guarded(MSK_linkfunctotaskstream(task,MSK_STREAM_LOG,NULL,printstr));
  // Give estimate of number of variables
  mosek_guarded(MSK_putmaxnumvar(task,n));
  if(m>0)
  {
    // Give estimate of number of constraints
    mosek_guarded(MSK_putmaxnumcon(task,m));
    // Give estimate of number of non zeros in A
    mosek_guarded(MSK_putmaxnumanz(task,Av.size()));
  }
  // Give estimate of number of non zeros in Q
  mosek_guarded(MSK_putmaxnumqnz(task,Qv.size()));
  if(m>0)
  {
    // Append 'm' empty constraints, the constrainst will initially have no
    // bounds
#if MSK_VERSION_MAJOR >= 7
    mosek_guarded(MSK_appendcons(task,m));
#elif MSK_VERSION_MAJOR == 6
    mosek_guarded(MSK_append(task,MSK_ACC_CON,m));
#endif
  }
  // Append 'n' variables
#if MSK_VERSION_MAJOR >= 7
  mosek_guarded(MSK_appendvars(task,n));
#elif MSK_VERSION_MAJOR == 6
  mosek_guarded(MSK_append(task,MSK_ACC_VAR,n));
#endif
  // add a contant term to the objective
  mosek_guarded(MSK_putcfix(task,cf));
 
  // loop over variables
  for(int j = 0;j<n;j++)
  {
    if(c.size() > 0)
    {
      // Set linear term c_j in the objective
      mosek_guarded(MSK_putcj(task,j,c[j]));
    }
 
    // Set constant bounds on variable j
    if(lx[j] == ux[j])
    {
      mosek_guarded(MSK_putbound(task,MSK_ACC_VAR,j,MSK_BK_FX,lx[j],ux[j]));
    }else
    {
      mosek_guarded(MSK_putbound(task,MSK_ACC_VAR,j,MSK_BK_RA,lx[j],ux[j]));
    }
 
    if(m>0)
    {
      // Input column j of A
#if MSK_VERSION_MAJOR >= 7
      mosek_guarded(
        MSK_putacol(
          task,
          j,
          Acp[j+1]-Acp[j],
          &Ari[Acp[j]],
          &Av[Acp[j]])
        );
#elif MSK_VERSION_MAJOR == 6
      mosek_guarded(
        MSK_putavec(
          task,
          MSK_ACC_VAR,
          j,
          Acp[j+1]-Acp[j],
          &Ari[Acp[j]],
          &Av[Acp[j]])
        );
#endif
    }
  }
 
  // loop over constraints
  for(int i = 0;i<m;i++)
  {
    // put bounds  on constraints
    mosek_guarded(MSK_putbound(task,MSK_ACC_CON,i,MSK_BK_RA,lc[i],uc[i]));
  }
 
  // Input Q for the objective (REMEMBER Q SHOULD ONLY BE LOWER TRIANGLE)
  mosek_guarded(MSK_putqobj(task,Qv.size(),&Qi[0],&Qj[0],&Qv[0]));
 
  // Set up task parameters
  for(
    std::map<MSKiparame,int>::iterator pit = mosek_data.intparam.begin();
    pit != mosek_data.intparam.end(); 
    pit++)
  {
    mosek_guarded(MSK_putintparam(task,pit->first,pit->second));
  }
  for(
    std::map<MSKdparame,double>::iterator pit = mosek_data.douparam.begin();
    pit != mosek_data.douparam.end(); 
    pit++)
  {
    mosek_guarded(MSK_putdouparam(task,pit->first,pit->second));
  }
 
  // Now the optimizer has been prepared
  MSKrescodee trmcode;
  // run the optimizer
  mosek_guarded(MSK_optimizetrm(task,&trmcode));
 
  //MSK_solutionsummary(task,MSK_STREAM_LOG);
 
  // Get status of solution
  MSKsolstae solsta;
#if MSK_VERSION_MAJOR >= 7
  MSK_getsolsta (task,MSK_SOL_ITR,&solsta);
#elif MSK_VERSION_MAJOR == 6
  MSK_getsolutionstatus(task,MSK_SOL_ITR,NULL,&solsta);
#endif
 
  bool success = false;
  switch(solsta)
  {
    case MSK_SOL_STA_OPTIMAL:   
    case MSK_SOL_STA_NEAR_OPTIMAL:
      MSK_getsolutionslice(task,MSK_SOL_ITR,MSK_SOL_ITEM_XX,0,n,&x[0]);
      //printf("Optimal primal solution\n");
      //for(size_t j=0; j<n; ++j)
      //{
      //  printf("x[%ld]: %g\n",j,x[j]);
      //}
      success = true;
      break;
    case MSK_SOL_STA_DUAL_INFEAS_CER:
    case MSK_SOL_STA_PRIM_INFEAS_CER:
    case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
    case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:  
      //printf("Primal or dual infeasibility certificate found.\n");
      break;
    case MSK_SOL_STA_UNKNOWN:
      //printf("The status of the solution could not be determined.\n");
      break;
    default:
      //printf("Other solution status.");
      break;
  }
 
  MSK_deletetask(&task);
  MSK_deleteenv(&env);
 
  return success;
}

References igl::mosek::MosekData::douparam, igl::mosek::MosekData::intparam, mosek_guarded(), and verbose.

Referenced by bbw().

Here is the call graph for this function:

Here is the caller graph for this function: