学习笔记
首发于学习笔记
基于WSTP的C++与Mathematica的混合编程及其应用

基于WSTP的C++与Mathematica的混合编程及其应用

基于WSTP的C++与Mathematica的混合编程及其应用

  • 本文严禁任何个人或网站转载

摘要

本文主要探讨在Wolfram Symbolic Transfer Protocol (简称为WSTP,其前身就是MathLink)下的C++与Mathematica的混合编程。文中会详细介绍编译器的配置以及通过几个实例附上源代码来帮助初学者快速将这个方法应用到自己的科学研究中。并且最后我会给出如何利用欧洲核子中心(CERN)的Minuit2拟合程序包来拟合多参数函数,当然其前提是将C++与Mathematica交互进行使用。

引言

Mathematica作为Wolfram Research的一款科学计算工具,与MatLab、Maple并称为三大数学软件。它很好地结合了数值与符号计算引擎、编程语言、文本、图形及声音等系统的处理,并且也提供有大量与外部程序进行交互的接口。在数学、物理、工程、经济、地理等领域有着极为重要的应用。尤其是它强大的符号计算能力更是让数学和物理学家从基础、复杂又重复的体力劳动般的公式推导中解脱了出来。

Mathematica本身提供有大量的数学函数,几乎可以满足各种数学要求。比如积分、微分、求极限、求和、解方程、对矩阵的各种操作等,可谓是科研工作者的人工智能助手。比如说,在粒子物理中,高阶的越来越复杂的费曼图的计算就意味着越来越复杂的符号计算,学过《量子场论》的大概都明白树图阶的康普顿散射真正动笔计算起来也还是很复杂很麻烦的。那么面对复杂的符号计算及积分,人手工来做显然无能为力。这个时候Mathematica可就能帮上大忙,比如只要载入高能包(HighEnergyPhysics)便可以进行复杂的洛伦兹指标的收缩或者伽马矩阵的求迹。

尽管Mathematica很强大很智能,但它并不是万能的。有时候在特定的科学领域,有一些特定的程序包是由C或者C++这种高级编译型语言所写的。当然,如果你不嫌浪费时间的话完全可以将C/C++的程序包翻译成Wolfram版本的程序包,但显然没人愿意这么做。

而C++作为面向对象的编程语言,在编程语言排行榜上几乎一直稳居前五,其编程时的强大的自由度使其自发明一来就一直具有强大的生命力。但是,对于数学及物理工作者而言,这并非一个很好的语言,因为它能提供的数学库函数中就只有最基本的基础函数而已,如果想从C++的数学库函数中调用一个求积分的或者求导数的那是不可能的,除非你自己用C++写一套可以求积分和导数的算法。

那么如果将两者结合起来的话岂不是会更加强大?答案是肯定的。

当然上面提到的问题在大多数情况下都不会发生,因为很多问题单独利用Mathematica或者C/C++就能解决,除非上面所提到的特定领域的特定程序是用C/C++写的,那么这个时候将两者结合起来将十分重要。比如在本文后面将着重介绍的CERN的Minuit2拟合程序包,就是用C++写的,后面我会详细给出此拟合程序的用法。

本文的结构按照如下方式进行安排:引言之后我会给出编译器的详细配置步骤,接着介绍几个简例并同时附上C++与Mathematica的源代码。最后重点介绍如何利用CERN的Minuit2拟合程序在将C++与Mathematica连接起来的情况下进行多参数拟合。

配置及简例

  • 配置

本文所使用的工具为:Mathematica11.0中文版与Visual Studio 2013中文版。

Mathematica安装包和Visual Studio 2013的安装包大家在Wolfram官网和Visual Studio官网上都可以下载到。先将下载好的Mathematica安装包和Visual Studio 2013的安装包进行安装,假设它们分别安装在D盘上的Mathematica11VisualStudio文件夹下,然后进行如下步骤(自Mathematica10.0以后,安装的时候Mathematica并不会自动将WSTP的动态链接库文件放到系统目录下,所以需要我们自己配置):

先是Mathematica的配置,其步骤如下:

  1. 打开D:\Mathematica11\SystemFiles\Links\WSTP\DeveloperKit你将会发现两个名为WindowsWindows-x86-64的文件夹,它们分别对应于32位系统与64位的操作系统。
  2. 接着打开\Windows\SystemAdditions,将其中的wstp32i4.dll文件拷贝到C:\Windows\System32目录下;再打开\Windows-x86-64\SystemAdditions,将其中的wstp64i4.dll拷贝到C:\Windows\SysWOW64目录下 (有时候这两个动态链接库文件在系统文件夹下的位置可能是反的,也就是如果到时候程序源代码提示“无法解析的外部符号”与“无法解析的外部命令”时,将两个的位置调换一下,我的电脑就是这种情况)。

接着是Visual Studio的配置,其步骤如下:
  1. D:\Mathematica11\SystemFiles\Links\WSTP\DeveloperKit\Windows\CompilerAdditions文件夹下面的wstp32i4m.lib文件拷贝到D:\VisualStudio\VC\lib文件夹下;将D:\Mathematica11\SystemFiles\Links\WSTP\DeveloperKit\Windows-x86-64\CompilerAdditions文件夹下面的wstp64i4m.lib文件拷贝到D:\VisualStudio\VC\lib\amd64文件夹下 (当然,也可以不用这么做,你可以将wstp32i4m.libwstp64i4m.lib放在任何你想放的文件夹下,只是在编译器中你需要将对应的路径加进去,至于怎么加,请自行百度Visual Studio中.lib库文件搜索路径的添加方法)。
  2. 除此之外,你需要将\Windows\CompilerAdditions下的wstp.h的头文件拷贝到D:\VisualStudio\VC\include文件夹下(当然,也可以不用这么做,你可以将wstp.h放在任何你想放的文件夹下,只是在编译器中你需要将对应的路径加进去,至于怎么加,请自行百度Visual Studio中.h头文件搜索路径的添加方法)。
  3. 当在Visual Studio中新建好一个工程之后,你需要在编译器的项目/属性/配置属性/连接器/输入/附加依赖项中输入wstp32i4m.lib,并与后面的“kernel32.lib;user32.lib;”用“;”隔开。

上面提到的Visual Studio的这种配置方法只是为了方便,将库文件和头文件分别放在Visual Studio的默认路径下,省去了每次新建新工程时都要配置一次的烦恼。

关于配置的更详细的资料请看:WSTP Development in C (Windows)

关于WSTP的详细介绍请看:WSTP and External Program Communication

配置好之后,下面就可以让C++与Mathematica进行通讯了。

  • 简例

  1. 从C++输入一个长方体的长(Length)、宽(Width)和高(Height)三个参数,将这三个参数传递给Mathematica,然后Mathematica计算好之后再传给C++,最后从屏幕上输出结果。

其对应的C++代码如下:

#include"wstp.h"
#include<iostream>
#include<math.h>
#include<cmath>
#include<vector>
using namespace std;

class Volume
{
public:
Volume(double len, double wid, double hei) :Length(len), Width(wid), Height(hei)
{
int argc = 4;
char *argv[5] = { "-linkname", "Volume", "-linkmode", "connect", NULL };
Link = WSOpen(argc, argv);
if ((WSLINK)0 == Link)
{
cout << " Unable to create the link..." << endl;
}
else
{
cout << "Link is Created successfully..." << endl;
}
cout << "Length=" << Length << "   " << "Width=" << Width << "   " << "Height=" << Height << endl;
}
~Volume(){};
double getLength() const{ return Length; }
double getWidth() const{ return Width; }
double getHeight() const{ return Height; }
vector<double>LWH()
{
lwh.push_back(getLength());
lwh.push_back(getWidth());
lwh.push_back(getHeight());
return lwh;
}
void ParametersTransfer(vector<double>&lwh)
{
WSPutReal64List(Link, &lwh[0], lwh.size());
}
virtual double operator()()
{
WSGetReal64(Link, &volume);
return volume;
}
virtual double CVolume()const
{
return getLength()*getWidth()*getHeight();
}
void LinkClose(){ WSClose(Link); cout << "The link has been shuted down..." << endl; }
private:
double Length;
double Width;
double Height;
WSLINK Link;
vector<double> lwh;
double volume;
};

void main()
{
WSENV env;
env = WSInitialize((WSEnvironmentParameter)0);
if ((WSENV)0 == env)
{
cout << "Unable to initialize environment..." << endl;
}
else
{
cout << "The environment is initialized successfully..." << endl;
}
Volume Cuboid(2.1,3.2,4.5);
cout << "The volume calculated with C language is:" << Cuboid.CVolume() << endl;
Cuboid.ParametersTransfer(Cuboid.LWH());
cout << "The volume calculated with Wolfram language is:" << Cuboid() << endl;
Cuboid.LinkClose();
WSDeinitialize(env);
system("pause");
}

至于对应的Mathematica代码,则简单地多,其对应的代码如下:

Link = LinkCreate["Volume", LinkProtocol -> "SharedMemory", LinkMode -> Listen];
{Len, Wid, Hei} = LinkRead[Link]; 
LinkWrite[Link, Len*Wid*Hei]; 
Print["Volume=", Len*Wid*Hei]


其对应的C++的输出为:

The environment is initialized successfully...
Link is Created successfully...
Length=2.1   Width=3.2   Height=4.5
The volume calculated with C language is:30.24
The volume calculated with Wolfram language is:30.24
The link has been shuted down...

对应的Mathematica的输出为:

Volume=30.24

2. 从C++输入f(x)=\sin^6 x,然后调用Mathematica求解此函数的积分,并返回结果。

其对应的C++代码如下:

#include"wstp.h"
#include<iostream>
#include<math.h>
#include<cmath>
#include<stdlib.h>
#include<string>
using namespace std;

class Integration
{
public:
Integration()
{
int argc = 4;
char *argv[5] = { "-linkname", "Integration", "-linkmode", "connect", NULL };
Link = WSOpen(argc, argv);
if ((WSLINK)0 == Link)
{
cout << " Unable to create the link..." << endl;
}
else
{
cout << "Link is Created successfully..." << endl;
}
}
~Integration(){}
void ExpToMathematica()
{
//向Mathematica发送表达式,使其计算符号积分
WSPutFunction(Link, "EvaluatePacket", 1);
WSPutFunction(Link, "ToExpression", 1);
WSPutString(Link, Expression);
WSEndPacket(Link);
}
void SymbExpToC()
{
//得到Mathematica计算后的符号积分结果,并输出
WSNewPacket(Link);
WSGetString(Link, &String);
cout << "The symbol form integration" << "'" << Expression << "'" << " calculated with Wolfram Language is:" << endl;
cout << String << endl;
WSReleaseString(Link, String);
		WSEndPacket(Link);
}
void SendxToMathematica(double x)
{
//向Mathematica发送要计算的x的值
WSNewPacket(Link);
WSPutReal64(Link, x);
WSEndPacket(Link);
}
void GetValueFromMathematica()
{
//得到Mathematica计算后的数值结果,并输出
WSNewPacket(Link);
WSGetReal64(Link, &Result);
cout << "The numerical form integration result calculated with Wolfram Language is:" << Result << endl;
WSEndPacket(Link);
}
void LinkClose()
{
WSClose(Link);
cout << "The link has been shuted down..." << endl;
}
private:
WSLINK Link;
const char *Expression = "f[x_]=Integrate[Sin[x]^6,x]";
const char *String;
double Result = 0.;
};

void main()
{
WSENV env;
env = WSInitialize((WSEnvironmentParameter)0);
if ((WSENV)0==env)
{ 
cout << "Unable to initialize environment..." << endl;
}
else
{
cout << "The environment is initialized successfully..." << endl;
}
Integration Intg;
Intg.ExpToMathematica();
Intg.SymbExpToC();
Intg.SendxToMathematica(3.2);
Intg.GetValueFromMathematica();
Intg.LinkClose();
WSDeinitialize(env);
system("pause");
}

其对应的Mathematica代码如下:

Link = LinkCreate["Integration", LinkProtocol -> "SharedMemory", 
   LinkMode -> Listen];
SymbolForm = LinkRead[Link];
LinkWrite[Link, ToString[f[x] // CForm]];
x = LinkRead[Link];
LinkWrite[Link, f[x]]

其对应的C++的输出为:

The environment is initialized successfully...
Link is Created successfully...
The symbol form integration'f[x_]=Integrate[Sin[x]^6,x]' calculated with Wolfram Language is:
(5*x)/16. - (15*Sin(2*x))/64. + (3*Sin(4*x))/64. - Sin(6*x)/192.
The numerical form integration result calculated with Wolfram Language is:0.981748
The link has been shuted down...

如果上面的代码你看的云里雾里,但工作中又要用到时,建议你看:

WSTP的API函数的详述:WSTP API-Wolfram Language Documentation

以及C++相关的教程

Minuit2的简介及应用

  • Minuit2的简介

Minuit2的前身即为Minuit,是一款数值最小化的计算程序,原始版本是CERN的物理学家Fred.James在20世纪70年代用Fortran语言写的,在21世纪初,Fortran版本被翻译成了更加流行的且面向对象的编程语言C++。最新的版本可以在Minuit官网的下载页面Minuit上看到。

更专业更详细的介绍请看:Minuit

下面介绍Minuit2的配置及应用:

首先从官网上下载下来最新版名为Minuit2-5.34.14.tar的压缩文件,将其解压后对于我们有用的文件为:incsrc以及test文件夹下的内容,其中最重要的是incsrc,因为test文件夹下面是一些使用的例子。

其次,我们需要将src文件下面的后缀为.cxx的原文件全部编译成库文件,具体做法是在Visual Studio上建立Win32项目,然后选择静态库,并取消预编译头选项。建好之后将src文件夹下的.cxx文件全部添加到资源文件(如果你想运行官网上mnusersguide.pdf这个帮助文档最后关于高斯函数的拟合程序的话,那你还需要将test/MnSim文件夹下面的GaussDataGen.cxxGaussFcn.cxx一起加到资源文件中)。


在库文件的生成过程中可能会出现以下问题

GenAlgoOptions.cxx 的原文件有错误,建议直接移除。

xutility.cpp 的原文件可能会警告代号为“error c4996”的错误,这是是编译器导致的,可以在预处理器定义中加入“_CRT_SECURE_NO_WARNINGS”即可。

如果编译顺利通过,那么假设我们生成了一个名为MinuitLib.lib的库文件,然后将它同样拷贝至D:\VisualStudio\VC\lib文件夹下。并在以后新建的工程中,在项目/属性/配置属性/连接器/输入/附加依赖项中输入MinuitLib.lib。
  • Minuit2的应用

假如说实验上测得一系列含有误差的数据,其从图上画出来如下所示:

然后我们的模型函数为f(x)=ax^2+bx+c,也就是说来求当\chi^2最小时候的二次函数的系数a,b,c.

其拟合程序的C++源代码如下:

#include<iostream>
#include<cmath>
#include<fstream>
#include<vector>
#include<string>
#include <cassert>
#include "wstp.h"
#include "Minuit2/FCNBase.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnPrint.h"
using namespace ROOT::Minuit2;
class FunctionFit : public FCNBase
{
public:
FunctionFit() :fErrorDef(1.)
{
env = WSInitialize((WSEnvironmentParameter)0);
if ((WSENV)0 == env)
{
std::cout << "Unable to initialize environment..." << std::endl;
}
else
{
std::cout << "The environment is initialized successfully..." << std::endl;
}
int argc = 4;
char *argv[5] = { "-linkname", "QuadraticFunctionFit", "-linkmode", "connect", NULL };
Link = WSOpen(argc, argv);
if ((WSLINK)0 == Link)
{
std::cout << " Unable to create the link..." << std::endl;
}
else
{
std::cout << "Link is Created successfully..." << std::endl;
}
WSGetInteger32(Link, &NumOfParas);
WSGetInteger32(Link, &dof);
std::cout << "The Number of Parameters is:" << NumOfParas << std::endl;
std::cout << "The Value of d.o.f is:" << dof << std::endl;
}
~FunctionFit(){}
virtual double Up()const { return fErrorDef; }
void SetErrorDef(double ErrorDef){ fErrorDef = ErrorDef; }
virtual double operator()(const std::vector<double>&Paras)const
{
double ReturnChiSqure;
assert(Paras.size() == NumOfParas);
std::for_each(Paras.begin(), Paras.end(), [](double item)
{
std::cout << item << "   " << std::endl;
});
WSPutReal64List(Link, &Paras[0], NumOfParas);
WSGetReal64(Link, &ReturnChiSqure);
std::cout << "chiSqure/d.o.f=" << ReturnChiSqure / dof << std::endl;
return ReturnChiSqure;
}
void LinkCloseAndDeinitialize(){ WSClose(Link); WSDeinitialize(env); std::cout << "The link has been shuted down..." << std::endl; }
private:
WSENV env;
WSLINK Link;
double fErrorDef;
int NumOfParas;
int dof;
};
void main()
{
FunctionFit QuadraticFunctionFit;
MnUserParameters UserParameters;
UserParameters.Add("a", 1.0, 1.0);
UserParameters.Add("b", 1.0, 1.0);
UserParameters.Add("c", 1.0, 1.0);
MnMigrad mnmigrad(QuadraticFunctionFit, UserParameters, 0);
FunctionMinimum Min = mnmigrad();
std::cout << Min.UserParameters() << std::endl;
std::ofstream Minout("Minout.txt");
Minout <<"The total number of function calls during the minimization is:"<< Min.NFcn() << std::endl;
Minout << Min.UserParameters() << std::endl;
Minout.close();
QuadraticFunctionFit.LinkCloseAndDeinitialize();
system("pause");
} 

对应的Mathematica源代码如下:

Needs["ErrorBarPlots`"]
f[x_] := a x^2 + b x + c;
ExpData = {{{-8, -220}, ErrorBar[10]}, {{-7, -148}, 
    ErrorBar[10]}, {{-6, -100}, ErrorBar[10]}, {{-5, -45}, 
    ErrorBar[10]}, {{-4, -6}, ErrorBar[10]}, {{-3, 20}, 
    ErrorBar[10]}, {{-2, 55}, ErrorBar[10]}, {{-1, 78}, 
    ErrorBar[10]}, {{0, 106}, ErrorBar[10]}, {{1, 115}, 
    ErrorBar[10]}, {{2, 125}, ErrorBar[10]}, {{3, 125}, 
    ErrorBar[10]}, {{4, 115}, ErrorBar[10]}, {{5, 108}, 
    ErrorBar[10]}, {{6, 78}, ErrorBar[10]}, {{7, 50}, 
    ErrorBar[10]}, {{8, 20}, ErrorBar[10]}, {{9, -5}, 
    ErrorBar[10]}, {{10, -45}, ErrorBar[10]}, {{11, -95}, 
    ErrorBar[10]}, {{12, -150}, ErrorBar[10]}, {{13, -220}, 
    ErrorBar[10]}};
DataPoints = 22; NumOfParas = 3; dof = DataPoints - NumOfParas;
Link = LinkCreate["QuadraticFunctionFit", 
   LinkProtocol -> "SharedMemory", LinkMode -> Listen];
LinkWrite[Link, NumOfParas]; LinkWrite[Link, dof];
Counter = 0;
While[Links["QuadraticFunctionFit"] != {},
 Counter = Counter + 1;
 {a, b, c} = LinkRead[Link];
 chiSqure = 
  Sum[(f[ExpData[[i, 1, 1]]] - ExpData[[i, 1, 2]])^2/10^2, {i, 1, 22}];
 Print["a=", a, "   ", "b=", b, "   ", "c=", c, "      ", 
  "\!\(\*SuperscriptBox[\(\[Chi]\), \(2\)]\)/d.o.f=", chiSqure/dof, 
  "   ", "计算进度为:", Counter];
 LinkWrite[Link, chiSqure];
 ]
LinkClose[Link];
ExpPlot = ErrorListPlot[ExpData, PlotTheme -> "Scientific"];
ModelPlot = Plot[-3.0451 x^2 + 15.2396 x + 101.7254, {x, -8, 13}];
Show[ExpPlot, ModelPlot]

下面给出拟合结果:

其中,橙色的带有误差棒的点为实验值,而蓝色曲线则为拟合出来的曲线。

对应的参数值为(C++的输出):

The total number of function calls during the minimization is:44
# ext. ||   Name    ||   type  ||     Value     ||  Error +/- 
   0   ||         a ||  free   ||     -3.045101637493 ||0.05940605706523
   1   ||         b ||  free   ||      15.23962450593 ||0.4485058952953
   2   ||         c ||  free   ||      101.7254376059 ||3.054649980729

a=-3.045\pm0.0594;b=15.2396\pm0.4485;c=101.725\pm3.0546

可见Minuit很“聪明”地找到了合理的参数值。事实上,上面的“实验值”是我根据函数f(x)=-3x^2+15x+100这个二次函数的函数值“制造”出来的(即,中心值稍有偏差,为方便起见将所有的误差都取为10)。在Minuit的参数初始值中,我将a,b,c三个参数的初始值和初始误差都设置成了1,但是Minuit还是很智能地找到了正确的值。

后记

最后,在这里大致说下写这篇文章的初衷:最近一段时间,我需要用计算得到的函数去拟合实验值,这个函数中含有五六个参数,后来了解到Minuit拟合包可以很好地处理这个问题,并且可以给出参数的误差。但是Minuit是用C++写的,而我的函数是用Mathematica推导得到的,这就急需将两者结合起来。虽然Mathematica有CForm[]函数,可以将对应的表达式转换成C语言的形式,但是函数中又含有积分,所以将两者结合,使其相互动态交换数据才是合适的办法。

一开始我也不知道该怎么弄,因为Mathematica还算比较熟悉,但是C++的话虽然之前零零散散看过一些,但好几年过去了,已然忘得一干二净。所以我就又将C++快速看了一下,再加上师弟Q. H在前期即编译器配置方面的帮助,及组内H. X之前做过相关的工作,有代码可以借鉴,所以很快就写出了测试代码,即上面的求体积和求积分的代码。然后再将该代码移植到Minuit程序拟合上,便得到了正确的答案。

刚开始我也在网上搜了一些相关的资料,但是很少,就算有也是说的云里雾里,有的根本就是全盘拷贝Mathematica官网上的说明资料。有的就算有部分代码,但也不全,还是干脆没法知道该怎么用。所以在我学会了之后就下决心要写这样的一篇文章,让初学者有C++和Mathematica的代码可看,看懂之后很容易就可以上手,并且根据我列出的参考资料,自己也可以进一步推广到自己的研究领域中去。上面的代码都是很基础的,当然要看懂上面的代码还是得具备C++和Mathematica的基本知识和编程经验。如果没有这方面的知识但又急需要用的话,我在这里悄悄告诉大家一个秘诀:还不赶紧滚回去看教材和资料?

最后的最后,文中的代码全部测试通过,但是从Visual Studio上拷贝过来之后又进行了编辑,没来得及检查,所以要是误删了分号或者编译提示错误,并且你找出了错误的话记得在评论中告知我,方便我修改。

明后天就要回家过年了,这篇文章是情急之中写出来的,所以不免有介绍不详细的地方或者错别字还没来得急更正,但是文中均给出了参考资料,所以看不懂的可以看参考资料中的更为详尽的内容。最后,祝愿大家鸡年大吉。

  • 本文严禁任何个人或网站转载
  • 如果本文给了你很大的帮助,帮你解了燃眉之急,那么请点赞和分享,或者你不吝啬的话可以给苦逼博士狗打赏点狗粮,本文已开启赞赏功能。
  • 喜欢作者本人的话”可以关注我的个人微博,个人微博可在我的知乎个人主页中的个人资料中找到链接。
编辑于 2017-12-03

文章被以下专栏收录

    本专栏旨在分享各种高质量的文章,包括在学习和研究过程中领悟出来的各种知识和技能,同时也欢迎其它领域的投稿。