EnergyMinimization.cpp

Go to the documentation of this file.

#include "EnergyMinimization.h"
#include "CirclePattern.h"
#ifdef _WIN32
#ifndef MSKDLL 
#define MSKDLL 
#endif
#endif
#include "mosek.h"   /* Include the MOSEK definition file. */
#include "util.h"
#include <fstream>
#include <iomanip>

using namespace std;


static int MSKAPI nlspar(void *nlhandle, int  *numgrdobjnz, int  *grdobjsub, int  i,
                            int  *convali, int  *grdconinz, int  *grdconisub, int  yo,
                            int  numycnz, int  ycsub[], int  maxnumhesnz,
                            int  *numhesnz, int  *hessubi, int  *hessubj) 
{

   Mesh * mesh = (Mesh *)nlhandle;
   Mesh::HalfEdgeIterator hIter = mesh->halfEdgeIterator();
   int NumInteriorEdges = mesh->numEdges() - mesh->numBoundary();

   if ( numgrdobjnz ) {
      numgrdobjnz[0] = mesh->numFaces();

      if ( grdobjsub ) {
          for (int i = 0; i <  mesh->numFaces(); i++)
             grdobjsub[i] = i;
      }
   }

   if ( convali )
      convali[0] = 0;    /* Zero because no nonlinear expression in the constraint. */

   if ( grdconinz )
      grdconinz[0] = 0;  /* Zero because no nonlinear expression in the constraint. */

   if ( numhesnz ) {
      if (yo)
          numhesnz[0] = NumInteriorEdges + mesh->numFaces();   
      else
          numhesnz[0] = 0;
   }

   if ( maxnumhesnz ) {
      check_error (maxnumhesnz < (NumInteriorEdges + mesh->numFaces()), "Not enough space was allocated for hessian.");
      if (yo) {
          if  (hessubi && hessubj) 
             CirclePattern::buildSparsityForHessian(mesh, hessubi, hessubj);
      }
   }
   return ( MSK_RES_OK );
} /* nlspar */


static int MSKAPI nleval(void   *nlhandle, double *xx, double yo, double *yc,
                            double *objval, int    *numgrdobjnz, int    *grdobjsub,
                            double *grdobjval, int    numi, int    *subi, double *conval,
                            int    *grdconptrb, int    *grdconptre, int    *grdconsub,
                            double *grdconval, double *grdlag, int    maxnumhesnz,
                            int    *numhesnz, int    *hessubi, int    *hessubj, 
                            double *hesval)
{
   Mesh * mesh = (Mesh *)nlhandle;
   int NumInteriorEdges = mesh->numEdges() - mesh->numBoundary();

   if ( objval) {
      objval[0] = CirclePattern::evaluateEnergy(mesh, xx);   
   }
   if ( numgrdobjnz ) {
      numgrdobjnz[0] = mesh->numFaces();                     
      for (int k = 0; k < mesh->numFaces(); k++) 
          grdobjsub[k] = k;
      CirclePattern::computeGradient(mesh, xx, grdobjval);
   }
   if ( conval )
      for(int k=0; k<numi; ++k)
          conval[k] = 0.0;

   if ( grdlag ) 
      CirclePattern::computeGradient(mesh, xx, grdlag);

   if ( maxnumhesnz ) {
      if ( yo==0.0 ) 
          numhesnz[0] = 0;
      else {                                                
          check_error (maxnumhesnz < (NumInteriorEdges + mesh->numFaces()), 
             "Not enough space was allocated for hessian.");
          numhesnz[0] = NumInteriorEdges + mesh->numFaces();

          if ( hessubi && hessubj && hesval ) {
             CirclePattern::buildSparsityForHessian(mesh, hessubi, hessubj);
             CirclePattern::computeHessian(mesh, xx, hesval);
          }
      }
   }
   return ( MSK_RES_OK );
} /* nleval */

EnergyMinimization::EnergyMinimization(Mesh * _mesh, const MSKenv_t & env) {
   mesh = _mesh;

   NUMVAR = mesh->numFaces();
   NUMCON = 1;
   NUMANZ = NUMVAR;

   int r = MSK_RES_OK;
   r = MSK_maketask(env, NUMCON, NUMVAR, &task);
   check_error (r != MSK_RES_OK, "Cannot make MOSEK task for radii energy minimization.");
   r = MSK_linkfiletotaskstream(task, MSK_STREAM_LOG, "_RadiiMinMosekOut.txt", 0);
   check_error (r != MSK_RES_OK, "Cannot link MOSEK \"_RadiiMinMosekOut.txt\" output file.");
}

void EnergyMinimization::doSolve() {
   allocateMosekArrays();
   setDataStructres();
   minimize();   
}

void EnergyMinimization::setDataStructres() {
   bkc[0] = MSK_BK_FX;
   blc[0] = 0.0;
   buc[0] = 0.0;

   for (int i = 0; i < NUMVAR; i++){
      c[i]     = 1.0;
      ptrb[i]  = i;  ptre[i] = i+1;
      sub[i]   = 0;  val[i]  = 1.0;

      bkx[i]   = MSK_BK_FR;
      blx[i]   = -MSK_INFINITY;
      bux[i]   = MSK_INFINITY;
   }
}

void EnergyMinimization::minimize() {
   int r = MSK_RES_OK;

   r = MSK_inputdata(task, NUMCON, NUMVAR, NUMCON, NUMVAR, c, 0.0, ptrb, ptre, sub, val, bkc, blc, buc, bkx, blx, bux);
   check_error ( r != MSK_RES_OK, 
      "Could not input constraints for radii minimization.\n Check  \"RadiiMinMosekOut.txt\" for output." );

   r = MSK_putnlfunc(task, mesh, nlspar, nleval);
   check_error ( r != MSK_RES_OK, 
      "Could not input nonlinear function for radii minimization.\n Check  \"RadiiMinMosekOut.txt\" for output." );

   r = MSK_optimize(task);
   check_error ( r != MSK_RES_OK, 
      "Could not solve for radii minimization.\n Check  \"RadiiMinMosekOut.txt\" for output." );

   r = MSK_getsolutionslice(task, MSK_SOL_ITR, MSK_SOL_ITEM_XX, 0, NUMVAR, xx);
   check_error ( r != MSK_RES_OK, "Could not get radii solution.\n Check  \"RadiiMinMosekOut.txt\" for output." );

   MSK_solutionsummary(task, MSK_STREAM_ERR);

   Mesh::FaceIterator fIter = mesh->faceIterator();
   for (int i = 0; !fIter.end(); fIter++) { 
      fIter.face()->r = exp(xx[i++]);
   }
}

void EnergyMinimization::allocateMosekArrays() {
   bkc = new int[NUMCON];
   bkx = new int[NUMVAR];
   ptrb = new int[NUMVAR];
   ptre = new int[NUMVAR];
   sub = new int[NUMANZ];
   blc = new double[NUMCON];
   buc = new double[NUMCON];
   c = new double[NUMVAR];
   blx = new double[NUMVAR];
   bux = new double[NUMVAR];
   val = new double[NUMANZ];
   xx = new double[NUMVAR];
}

void EnergyMinimization::clearMosekArrays() {
   if (task) {

      delete [] bkc;
      delete [] bkx;
      delete [] ptrb;
      delete [] ptre;
      delete [] sub;
      delete [] blc;
      delete [] buc;
      delete [] c;
      delete [] blx;
      delete [] bux;
      delete [] val;
      MSK_deletetask(&task);
      task = NULL;
   }
}

---