Skip to main content

mathlink or wstp - Calling Mathematica from Visual Studio - example


I am trying to set up a connection between Visual Studio 10 and Mathematica 8.0 on a Windows system. I found several sources of information on the Wolfram website, but I get entirely lost. Most of the info deals with calling C++ functions from within Mathematica, but I would like it the other way round. I have quite a lot of functionality in my C++ code and I just call Mathematica to make some integral calculations.


Is there a good step-by-step approach to set up this connection? What files should I include in my project? What settings need to be changed? Does anyone have a complete example of pieces of code that gives insights in the general approach? That would be a great help.


I followed the standard instructions to add the lib, bin and include to the Visual Studio directories. More specifically, I was also wondering whether in this case I also need to make a .tm file? Should it be empty then? The example here looks very clear, but this is approached from a Linux/Mac perspective.



I hope someone could bring me some guidance.



Taking into account your suggestions, I have developed the following program, but there seems to be a problem with retrieving the packages. I have put the entire code here.


The program first sends some variables to Mathematica kernel and then it should perform the integral calculation. I am only interested in the final payoff that is calculated by Mathematica. While debugging, it seems to be that the program gets stuck in the removal of all the packages (in the loop just after sending the integral function). There is no problem with the function, I guess, because it performs the calculation perfectly when I just use Mathematica.


#include 
#include 

#include "mathlink.h"

static void init_and_openlink( int argc, char* argv[]);

static void error( MLINK lp);


MLENV ep = (MLENV)0;
MLINK lp = (MLINK)0;


int main(int argc, char* argv[]){
    int pkt;
    float payoff;


    init_and_openlink( argc, argv);

Now, I send the different parameters and the final function of interest:


    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);
        MLPutString(lp, "mue");
        MLPutReal(lp,0.25);
    MLEndPacket(lp);


    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);
        MLPutString(lp, "mui");
        MLPutReal(lp,0.25);
    MLEndPacket(lp);

    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);
        MLPutString(lp, "lambdae");
        MLPutReal(lp,0.25);

    MLEndPacket(lp);

    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);
        MLPutString(lp, "lambdai");
        MLPutReal(lp,0.25);
    MLEndPacket(lp);

    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);

        MLPutString(lp, "sigmap");
        MLPutReal(lp,0.05);
    MLEndPacket(lp);

    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);
        MLPutString(lp, "betae");
        MLPutReal(lp,0.05);
    MLEndPacket(lp);


    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);
        MLPutString(lp, "betai");
        MLPutReal(lp,0.05);
    MLEndPacket(lp);

    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);
        MLPutString(lp, "gammae");
        MLPutReal(lp,0.05);

    MLEndPacket(lp);

    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);
        MLPutString(lp, "gammai");
        MLPutReal(lp,0.10);
    MLEndPacket(lp);

    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);

        MLPutString(lp, "d");
        MLPutReal(lp,0.2);
    MLEndPacket(lp);

    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);
        MLPutString(lp, "experience1");
        MLPutInteger(lp,2);
    MLEndPacket(lp);


    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"SetDelayed",2);
        MLPutString(lp, "experience2");
        MLPutInteger(lp,3);
    MLEndPacket(lp);

    MLPutFunction( lp,"EnterTextPacket",1);
        MLPutString(lp,"expectedpayoff1bis[i1_, b1_,i2_,b2_] := NIntegrate[PDF[NormalDistribution[(1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100), (1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100)* Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience1])^2 + (gammai*Exp[-lambdai*i1])^2]], x]*(1 - CDF[NormalDistribution[(1 + betae*Exp[-mue*experience2] + betai*Exp[-mui*i2])*(1 + b2/100), (1 +  betae*Exp[-mue*experience2] + betai*Exp[-mui*i2])*(1 + b2/100)*Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience2])^2 + (gammai*Exp[-lambdai*i2])^2]], x]), {x, 0, 2}]*(NIntegrate[(x*PDF[NormalDistribution[(1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100), (1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100)* Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience1])^2 + (gammai*Exp[-lambdai*i1])^2]], x]*(1 - CDF[NormalDistribution[(1 + betae*Exp[-mue*experience2] + betai*Exp[-mui*i2])*(1 + b2/100), (1 + betae*Exp[-mue*experience2] + betai*Exp[-mui*i2])*(1 + b2/100)*Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience2])^2 + (gammai*Exp[-lambdai*i2])^2]], x])/NIntegrate[PDF[NormalDistribution[(1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100), (1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100)*Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience1])^2 + (gammai*Exp[-lambdai*i1])^2]], x]*(1 - CDF[NormalDistribution[(1 + betae*Exp[-mue*experience2] +betai*Exp[-mui*i2])*(1 + b2/100), (1 + betae*Exp[-mue*experience2] + betai*Exp[-mui*i2])*(1 + b2/100)*Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience2])^2 + (gammai*Exp[-lambdai*i2])^2]], x]), {x, 0, 2}]), {x, 0, 2}] - (1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1])) - d*i1/100) - (1 - d)*i1/100");
    MLEndPacket(lp);


I don't need any of the information retrieved from the above packages, so I try to get rid of everything. This is were the debugging is in a loop (When I pause the debugging, I get the screen with "no source available".


    while ((pkt=MLNextPacket(lp),pkt)&&pkt!=RETURNPKT) {
        MLNewPacket(lp);
        if (MLError(lp)) error(lp);
    }
    MLNewPacket(lp);

Afterwards I want to evaluate the function and the result needs to be printed on the screen. This is the only output I am interested in:


    MLPutFunction( lp,"EnterExpressionPacket",1);
        MLPutFunction(lp,"expectedpayoff1bis",4);

        MLPutInteger(lp,2);
        MLPutInteger(lp,10);
        MLPutInteger(lp,4);
        MLPutInteger(lp,12);
    MLEndPacket(lp);

    while( (pkt = MLNextPacket( lp), pkt) && pkt != RETURNPKT) {
        MLNewPacket( lp);
        if (MLError( lp)) error( lp);
    }

    MLGetReal( lp,&payoff);


    printf( "payoff: %f", payoff);


    MLPutFunction( lp, "Exit", 0);

    return 0;
}


The rest of the code is retrieved from the examples


static void error( MLINK lp){
    if( MLError( lp)){
        fprintf( stderr, "Error detected by MathLink: %s.\n",MLErrorMessage(lp));
    }
    else{
        fprintf( stderr, "Error detected by this program.\n");
    }
    exit(3);

}


static void deinit( void){
    if( ep) MLDeinitialize( ep);
}


static void closelink( void){
    if( lp) MLClose( lp);

}


static void init_and_openlink( int argc, char* argv[]){
    #if MLINTERFACE >= 3
    int err;
    #else
    long err;
    #endif /* MLINTERFACE >= 3 */


    ep =  MLInitialize( (MLParametersPointer)0);
    if( ep == (MLENV)0) exit(1);
    atexit( deinit);

    #if MLINTERFACE < 3
    lp = MLOpenArgv( ep, argv, argv + argc, &err);
    #else
    lp = MLOpenArgcArgv( ep, argc, argv, &err);
    #endif
    if(lp == (MLINK)0) exit(2);

    atexit( closelink); 
}

Answer



There are a few problems with your code. If you fix those up, as I did, your program will run fine.


First off, the reason you see "no source available" when you pause the program is probably that when you break, the program is down inside a MathLink function, and so it is complaining that it doesn't have access to the MathLink library source code. To debug your program, you generally set breakpoints in your code and step through it, and not try break in at an arbitrary time, which often will be inside some library that you did not write.


Back to the code. The most glaring error is that you send a bunch of evaluations to the kernel, but fail to read the packets that come back from each one. When you go to read what you think is your final result, you will actually be reading some earlier result that you never drained off the link. Every interaction with Mathematica will involve you sending a packet and then immediately draining off all the packets that come back. You will note that in my program I have written a discardResult() function that drains all packets up to, and including, the ReturnPacket containing the ignored result.


Don't use EnterExpressionPacket as the packet when sending things to the kernel. This engages the whole main loop, including assignment of In/Out, etc., and you get a different sequence of packets coming back. Better to use EvaluatePacket, from which you are guaranteed to get back a single ReturnPacket containing the result (you might also get other packets along the way, such as for messages and Print output).


Another error is that you are using MLPutString() to send your symbol names (mue, mui, lambdae, etc.) These are symbols, so use MLPutSymbol().


I also changed your use of SetDelayed to Set, just to be picky.


You used EnterTextPacket for the definition of expectedpayoff1bis so that you could send it as a string. Like EnterExpressionPacket, EnterTextPacket engages the main loop, and gives a different packet sequence in return, so avoid it. The best way to send Mathematica input in textual format is to use EvaluatePacket[ToExpression["string of code"]].



Here is my modified version of your program. It runs and gives a numerical answer for me.


#include 
#include 
#include "mathlink.h"

static void init_and_openlink( int argc, char* argv[]);
static void error( MLINK lp);
static void discardResult(MLINK lp);

MLENV ep = (MLENV)0;

MLINK lp = (MLINK)0;


int main(int argc, char* argv[])
{
    int pkt;
    double payoff;

    init_and_openlink( argc, argv);


    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "mue");
        MLPutReal(lp,0.25);
    MLEndPacket(lp);
    discardResult(lp);

    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "mui");

        MLPutReal(lp,0.25);
    MLEndPacket(lp);
    discardResult(lp);

    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "lambdae");
        MLPutReal(lp,0.25);
    MLEndPacket(lp);
    discardResult(lp);


    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "lambdai");
        MLPutReal(lp,0.25);
    MLEndPacket(lp);
    discardResult(lp);

    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);

        MLPutSymbol(lp, "sigmap");
        MLPutReal(lp,0.05);
    MLEndPacket(lp);
    discardResult(lp);

    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "betae");
        MLPutReal(lp,0.05);
    MLEndPacket(lp);

    discardResult(lp);

    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "betai");
        MLPutReal(lp,0.05);
    MLEndPacket(lp);
    discardResult(lp);

    MLPutFunction( lp,"EvaluatePacket",1);

        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "gammae");
        MLPutReal(lp,0.05);
    MLEndPacket(lp);
    discardResult(lp);

    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "gammai");
        MLPutReal(lp,0.10);

    MLEndPacket(lp);
    discardResult(lp);

    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "d");
        MLPutReal(lp,0.2);
    MLEndPacket(lp);
    discardResult(lp);


    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "experience1");
        MLPutInteger(lp,2);
    MLEndPacket(lp);
    discardResult(lp);

    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"Set",2);
        MLPutSymbol(lp, "experience2");

        MLPutInteger(lp,3);
    MLEndPacket(lp);
    discardResult(lp);

    MLPutFunction( lp,"EvaluatePacket",1);
      MLPutFunction( lp,"ToExpression",1);
        MLPutString(lp,"expectedpayoff1bis[i1_, b1_,i2_,b2_] := NIntegrate[PDF[NormalDistribution[(1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100), (1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100)* Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience1])^2 + (gammai*Exp[-lambdai*i1])^2]], x]*(1 - CDF[NormalDistribution[(1 + betae*Exp[-mue*experience2] + betai*Exp[-mui*i2])*(1 + b2/100), (1 +  betae*Exp[-mue*experience2] + betai*Exp[-mui*i2])*(1 + b2/100)*Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience2])^2 + (gammai*Exp[-lambdai*i2])^2]], x]), {x, 0, 2}]*(NIntegrate[(x*PDF[NormalDistribution[(1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100), (1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100)* Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience1])^2 + (gammai*Exp[-lambdai*i1])^2]], x]*(1 - CDF[NormalDistribution[(1 + betae*Exp[-mue*experience2] + betai*Exp[-mui*i2])*(1 + b2/100), (1 + betae*Exp[-mue*experience2] + betai*Exp[-mui*i2])*(1 + b2/100)*Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience2])^2 + (gammai*Exp[-lambdai*i2])^2]], x])/NIntegrate[PDF[NormalDistribution[(1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100), (1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1]))*(1 + b1/100)*Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience1])^2 + (gammai*Exp[-lambdai*i1])^2]], x]*(1 - CDF[NormalDistribution[(1 + betae*Exp[-mue*experience2] +betai*Exp[-mui*i2])*(1 + b2/100), (1 + betae*Exp[-mue*experience2] + betai*Exp[-mui*i2])*(1 + b2/100)*Sqrt[sigmap^2 + (gammae*Exp[-lambdae*experience2])^2 + (gammai*Exp[-lambdai*i2])^2]], x]), {x, 0, 2}]), {x, 0, 2}] - (1 + (betae*Exp[-mue*experience1] + betai*Exp[-mui*i1])) - d*i1/100) - (1 - d)*i1/100");
    MLEndPacket(lp);
    discardResult(lp);


    MLPutFunction( lp,"EvaluatePacket",1);
        MLPutFunction(lp,"expectedpayoff1bis",4);
        MLPutInteger(lp,2);
        MLPutInteger(lp,10);
        MLPutInteger(lp,4);
        MLPutInteger(lp,12);
    MLEndPacket(lp);

    while( (pkt = MLNextPacket( lp), pkt) && pkt != RETURNPKT) {
        MLNewPacket( lp);

        if (MLError( lp)) error( lp);
    }
    MLGetReal( lp,&payoff);

    printf( "payoff: %f", payoff);

    return 0;
}



static void discardResult(MLINK lp) {
    int pkt;
    while( (pkt = MLNextPacket( lp), pkt) && pkt != RETURNPKT) {
        MLNewPacket( lp);
        if (MLError( lp)) error( lp);
    }
    MLNewPacket(lp);
}



static void error( MLINK lp){
    if( MLError( lp)){
        fprintf( stderr, "Error detected by MathLink: %s.\n",MLErrorMessage(lp));
    }
    else{
        fprintf( stderr, "Error detected by this program.\n");
    }
    exit(3);
}



static void deinit( void){
    if( ep) MLDeinitialize( ep);
}


static void closelink( void){
    if( lp) MLClose( lp);
}



static void init_and_openlink( int argc, char* argv[]){
    #if MLINTERFACE >= 3
    int err;
    #else
    long err;
    #endif /* MLINTERFACE >= 3 */

    ep =  MLInitialize( (MLParametersPointer)0);
    if( ep == (MLENV)0) exit(1);

    atexit( deinit);

    #if MLINTERFACE < 3
    lp = MLOpenArgv( ep, argv, argv + argc, &err);
    #else
    lp = MLOpenArgcArgv( ep, argc, argv, &err);
    #endif
    if(lp == (MLINK)0) exit(2);
    atexit( closelink); 
}

Comments

Post a Comment

Popular posts from this blog

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]
Post a Comment
Read more

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1.
Post a Comment
Read more