Difference between revisions of "ApCoCoA-1:MatlabInterface"

From ApCoCoAWiki
m (added link to Compile instructions)
Line 34: Line 34:
 
#include <vector>
 
#include <vector>
 
using std::vector;
 
using std::vector;
extern "C" __declspec(dllexport) int PreProcess(double
 
*tolerance, int rows, int cols, double *OrigArray, double
 
*PreProcessedArray, int aType)
 
{
 
// Variables
 
vector<double> toler(cols);
 
vector<ApproxPts::ApproxPt> OrigPts(rows,
 
  
ApproxPts::ApproxPt(cols));
+
extern "C" __declspec(dllexport) int PreProcess(double *tolerance, int rows, int cols, double *OrigArray, double *PreProcessedArray, int aType) {
vector<ApproxPts::ApproxPt> PreprocessedPts;
+
// Variables
// Define toler parameter
+
vector<double> toler(cols);
for (int i=0; i < cols; i++)
+
vector<ApproxPts::ApproxPt> OrigPts(rows, ApproxPts::ApproxPt(cols));
{toler[i]=tolerance[i];}
+
vector<ApproxPts::ApproxPt> PreprocessedPts;
// Rearrange original points to fit needed parameter structure
+
// Build toler parameter
for (int i=0; i < rows; ++i)
+
for (int i=0; i < cols; i++)
for (int j=0; j < cols; ++j)
+
  {toler[i]=tolerance[i];}
{OrigPts[i][j]=OrigArray[i*cols+j];}
+
// Copy original data to vector
// Call Preprocess function in ApCoCoALib
+
for (int i=0; i < rows; ++i) {
switch ( aType ) {
+
  for (int j=0; j < cols; ++j){
case 1 :  
+
    OrigPts[i][j]=OrigArray[j*rows+i];
ApproxPts::PreprocessGridAlgm(PreprocessedPts, OrigPts, toler);
+
  }
break;
+
}
case 2 :  
+
// Call Preprocess function in ApCoCoALib
ApproxPts::PreprocessAggrAlgm(PreprocessedPts, OrigPts, toler);
+
switch ( aType ) {
break;
+
  case 1 :  
case 3:
+
    ApproxPts::PreprocessGridAlgm(PreprocessedPts, OrigPts, toler);
                ApproxPts::PreprocessSubdivAlgm(PreprocessedPts, OrigPts, toler);
+
    break;
break;
+
  case 2 :  
}
+
    ApproxPts::PreprocessAggrAlgm(PreprocessedPts, OrigPts, toler);
// Rearrange preprocessed points
+
    break;
for (int i=0; i < PreprocessedPts.size(); i++)
+
  case 3:
for (int j=0; j < cols; j++)
+
    ApproxPts::PreprocessSubdivAlgm(PreprocessedPts, OrigPts, toler);
{PreProcessedArray[i*cols+j] = PreprocessedPts[i][j];}
+
    break;
// Return number of calculated rows/points
+
}
return PreprocessedPts.size();
+
// Rearrange preprocessed points
}</cpp>
+
int rowsRes = PreprocessedPts.size();
 +
for (int i=0; i < rowsRes; i++) {
 +
  for (int j=0; j < cols; j++) {
 +
  PreProcessedArray[j*rowsRes+i] = PreprocessedPts[i][j];
 +
  }
 +
}
 +
// Return number of calculated rows/points
 +
return rowsRes;
 +
}</cpp>
  
 
and PreProcessing.c
 
and PreProcessing.c
<c> #include "mex.h"
+
<c>#include "mex.h"
#include "matrix.h"
+
#include "matrix.h"
__declspec(dllexport) int PreProcess(double *tolerance, int rows, int cols,  
+
__declspec(dllimport) int PreProcess(double *tolerance, int rows, int cols, double *OrigArray, double *PreProcessedArray, int aType);
                    double *OrigArray, double *PreProcessedArray, int aType);
 
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const
 
  
mxArray *prhs[])
+
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
+
{
// Variables
+
// Variables
mxArray *Toler, *Rows, *Cols, *OrigPts, *AlgType;
+
mxArray *Toler, *OrigPts, *AlgType;
double *preProcessed, *origData, *result, *tol;
+
double *preProcessed, *origData, *result, *tol;
int rowLen, colLen, strLen, numPP, nElPP, nEl, i, algo;
+
int rowLen, colLen, strLen, numPP, nElPP, nEl, i, algo;
/*  
+
/*  
Check for proper number of input arguments.  
+
Check for proper number of input arguments.  
Required are:  
+
Required are:  
1. Tooler (tolerance vector, for each dimension)
+
  1. Tooler (tolerance vector, for each dimension)
2. Rows (of original data matrix - points of data)
+
  4. OrigPoints (Pointer for data, 1xn matrix)
3. Cols (of original data matrix - dimension of data)
+
  5. AlgType (1 = Grid, 2 = Aggr, 3 = Subdiv)
4. OrigPoints (Pointer for data, 1xn matrix)
+
*/     
5. AlgType (1 = Grid, 2 = Aggr, 3 = Subdiv)
+
if (nrhs != 3)  
*/     
+
  {mexErrMsgTxt("3 input arguments required.");}
if (nrhs != 5)  
+
// Get input data
{mexErrMsgTxt("5 input arguments required.");}
+
Toler = prhs[0];
// Get input data
+
OrigPts = prhs[1];
Toler = prhs[0];
+
AlgType = prhs[2];
Rows = prhs[1];
+
tol=mxGetPr(Toler);
Cols = prhs[2];
+
origData = mxGetPr(OrigPts);
OrigPts = prhs[3];
+
rowLen = (int) mxGetM(OrigPts);
AlgType = prhs[4];
+
colLen = (int) mxGetN(OrigPts);
tol=mxGetPr(Toler);
+
algo = mxGetScalar(AlgType);
rowLen = (int) mxGetScalar(Rows);
+
// Prepare array for result
colLen = (int) mxGetScalar(Cols);
+
nEl = rowLen*colLen;
origData = mxGetPr(OrigPts);
+
preProcessed = malloc(nEl * sizeof(double));
algo = mxGetScalar(AlgType);
+
// Run calculation
// Prepare array for result
+
numPP = PreProcess(tol, rowLen, colLen, origData, preProcessed, algo);
nEl = rowLen*colLen;
+
nElPP = numPP*colLen;
preProcessed = malloc(nEl * sizeof(double));
+
// Return preprocessed points
// Run calculation
+
plhs[0] = mxCreateDoubleMatrix(numPP, colLen, mxREAL);
numPP = PreProcess(tol, rowLen, colLen, origData, preProcessed, algo);
+
result = mxGetPr(plhs[0]);
nElPP = numPP*colLen;
+
for (i=0; i < nElPP; ++i)
// Return preprocessed points
+
{
plhs[0] = mxCreateDoubleMatrix(1, nElPP, mxREAL);
+
  result[i] = preProcessed[i];
for (i=0; i < nElPP; i++)
+
}
{
+
// clean up
result = mxGetPr(plhs[0]);
+
free(preProcessed);
result[i] = preProcessed[i];
+
}</c>
}
 
// clean up
 
free(preProcessed);
 
}</c>
 
  
 
We use '''mex.bat''' to compile the second file into a mex plugin (as '''C''' code) while we use '''cl.exe''' to compile the first file into a DLL.
 
We use '''mex.bat''' to compile the second file into a mex plugin (as '''C''' code) while we use '''cl.exe''' to compile the first file into a DLL.

Revision as of 11:46, 29 November 2007

MatlabToolbox Interface

The ApCoCoALib is written in C++ and therefore the potential number of users is limited to people who are comfortable writing C++ code. But many users of Computer Algebra systems or other software might benefit from the functionality provided by the CoCoALib. One such piece of software is Matlab which among other things allows the user to integrate C, Fortran and C++ libraries. The integration is not too difficult for libraries written in C or Fortran but there are severe restrictions for libraries written in C++. Due to that fact we are currently working on the necessary bits to provide the functionality of the CoCoALib in MatLab.

We want to provide the MatLabToolbox as binary blob as well as in source. Since the CoCoALib runs on a wide number of 32 and 64 bit platforms we expect to provide binaries for

  • x86 & x86-64 Linux
  • x86 & x86-64 Windows
  • 32bit PPC MacOSX
  • Sparc Solaris

This obviously depends on access to appropriate hardware and MatLab binaries.

Downloads

Caution: The current source codes are prepared to compile with MSVC 2005 only. The use of any other compiler but MSVC 2005 will lead to spectacular link failures since one cannot link C++ binaries produced by MSVC 2005 with for example MSVC 2003. You also need to have the runtime of MSVC 2005 installed on your system to make this work. In case you experience problems you can always ask questions in the ApCoCoA forum.

The process of building the MatlabInterface is described at Compile on Windows. Links for downloads are also offered there.

Design

C++ with exceptions and RTTI vs. MatLab

The design is quite complex due to the fact that integrating C++ code into a mex file has severe restrictions, most significantly:

  • Limitation to g++ 3.2.2 on Linux. Using g++ 3.2.1 or 3.2.3 might lead to crashes or unexpected behaviour in general. While this might have been ok two years ago most current systems now ship with g++ 4.0 or 4.1.
  • The Microsoft compiler does disable exception handling or RTTI per default. Mex also explicitly disables exceptions in C++ code in the mexopts.bat file. While one obviously is free to change that we have been unable to create a mex plugin compiled with exceptions and RTTI. Since those two features are essential to the CoCoALib a different way had to be found.

It has also been reported that migrating C++ code from R13 to R14 has usually not been as easy as one could hope. This is obviously only anecdotal evidence, as each of the google groups on the other hands provides plenty of people seeking help.

Using C wrapper around the ApCoCoALib

The solution is to use extern "C" functions inside a second DLL. That second DLL can be compiled with exceptions and RTTI. For a demonstration look at the following two code snippets demonstrating the integration of the PreProcessing function:

ApCoCoALib.cpp: <cpp>#include "CoCoA/library.H" using namespace CoCoA;

  1. include <vector>

using std::vector;

extern "C" __declspec(dllexport) int PreProcess(double *tolerance, int rows, int cols, double *OrigArray, double *PreProcessedArray, int aType) { // Variables vector<double> toler(cols); vector<ApproxPts::ApproxPt> OrigPts(rows, ApproxPts::ApproxPt(cols)); vector<ApproxPts::ApproxPt> PreprocessedPts; // Build toler parameter for (int i=0; i < cols; i++)

 {toler[i]=tolerance[i];}

// Copy original data to vector for (int i=0; i < rows; ++i) {

 for (int j=0; j < cols; ++j){
   OrigPts[i][j]=OrigArray[j*rows+i];
 }

} // Call Preprocess function in ApCoCoALib switch ( aType ) {

 case 1 : 
   ApproxPts::PreprocessGridAlgm(PreprocessedPts, OrigPts, toler);
   break;
 case 2 : 
   ApproxPts::PreprocessAggrAlgm(PreprocessedPts, OrigPts, toler);
   break;
 case 3:
   ApproxPts::PreprocessSubdivAlgm(PreprocessedPts, OrigPts, toler);
   break;

} // Rearrange preprocessed points int rowsRes = PreprocessedPts.size(); for (int i=0; i < rowsRes; i++) {

 for (int j=0; j < cols; j++) {

PreProcessedArray[j*rowsRes+i] = PreprocessedPts[i][j];

 }

} // Return number of calculated rows/points return rowsRes; }</cpp>

and PreProcessing.c <c>#include "mex.h"

  1. include "matrix.h"

__declspec(dllimport) int PreProcess(double *tolerance, int rows, int cols, double *OrigArray, double *PreProcessedArray, int aType);

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])

{ // Variables mxArray *Toler, *OrigPts, *AlgType; double *preProcessed, *origData, *result, *tol; int rowLen, colLen, strLen, numPP, nElPP, nEl, i, algo; /* Check for proper number of input arguments. Required are:

 1. Tooler (tolerance vector, for each dimension)
 4. OrigPoints (Pointer for data, 1xn matrix)
 5. AlgType (1 = Grid, 2 = Aggr, 3 = Subdiv)
  • /

if (nrhs != 3)

 {mexErrMsgTxt("3 input arguments required.");}

// Get input data Toler = prhs[0]; OrigPts = prhs[1]; AlgType = prhs[2]; tol=mxGetPr(Toler); origData = mxGetPr(OrigPts); rowLen = (int) mxGetM(OrigPts); colLen = (int) mxGetN(OrigPts); algo = mxGetScalar(AlgType); // Prepare array for result nEl = rowLen*colLen; preProcessed = malloc(nEl * sizeof(double)); // Run calculation numPP = PreProcess(tol, rowLen, colLen, origData, preProcessed, algo); nElPP = numPP*colLen; // Return preprocessed points plhs[0] = mxCreateDoubleMatrix(numPP, colLen, mxREAL); result = mxGetPr(plhs[0]); for (i=0; i < nElPP; ++i) {

 result[i] = preProcessed[i];

} // clean up free(preProcessed); }</c>

We use mex.bat to compile the second file into a mex plugin (as C code) while we use cl.exe to compile the first file into a DLL.

Functions provided

Obviously the goal is to provide every available function in the ApCoCoALib to users of Matlab. Due to our own needs at the moment the first functions to work will be:

  • PreProcessing
  • Approximate Buchberger-Möller
  • Syzygy computation for ideals

PreProcessing

Due to the current implementation of PreProcessing.c the data needs to be a vector but most of the data we work with in Matlab is in matrix form. Thus the following code might be usefull to use the PreProcessing function in your Matlab codes:

PrepPreProcessing.m <matlab> function result = PrepPreProcessing(tolerance, data, algtype)

   [nRows, nCols] = size(data);
   transferedData = reshape(data',1,nRows*nCols);
   resultPP = PreProcessing(tolerance,nRows,nCols,transferedData,algtype);
   nRows2 = length(resultPP)/nCols;
   result = reshape(resultPP,nCols,nRows2)';</matlab>