11/17/2016

ExtendingMATLABwithNumericalRecipes

ExtendingM

withNumericalRecipes

[[Link],youcanprintthedocumentfromthisPDFversion.]

NumericalRecipescode,orforthatmatteranyotherC++code,caneasilybeinvokedfromwithinMATLAB,from
theconsole,[Link],C++executeshugelyfasterthan
MATLABmcode,andyoucanalsoaccessfeaturesinNumericalRecipes(orelsewhere)thatarenotavailablein
[Link]++,andcontrolthewhole
projectfromtheMATLABconsoleorcommandline.
Thisdocumentshowshowthisisdone,withexamplesoftwodifferentcodingmethodologies:(1)"API",using
theMATLABCAPIfunctioncalls,and(2)"NR3",usingthespecialNumericalRecipesincludefilenr3matlab.h.
Wethinkthat"NR3"isthebestchoice(it'ssimplertocode)butwealsoshow"API"examplessothatyoucan
understandtheunderlyingMATLABinterface.

Contents
ConventionsandOverview
"Hello,world!"UsingtheMatlabCAPI
DoJustOnce:IntroduceMatlabtoYourCompiler
TheNR3CodingMethod
"Hello,world!"Usingnr3matlab.h
AccessingMatlabVectorswithnr3matlab.h
LivingDangerously
LivingSafely
AccessingMatlabScalarswithnr3matlab.h
Matrices:"ThroughtheLookingGlass"(TTLG)
AccessingMatlabMatriceswithnr3matlab.h
WrapperFunctionsThatAccessMultipleMemberFunctionsinaClass
TheAPICodingMethod
TipsonMatlabCAPIProgramming
NumericalRecipesUsingtheMatlabCAPI
Appendix:UsingMicrosoftVisualStudio

ConventionsandOverview
MATLABmcode,ortextfromtheMATLABconsoleisshownonaredbackground:
% MATLAB code
[a b] = someroutine(c,d)
C++codeisshownonagreenbackground:
/* C++ code */
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
// don't worry yet about these weird types!
}
Acompiled,executablefilethatisinvokedbyaMATLABcommandiscalleda"mexfile".Mexfilesarewrittenas
sourcecodeoutsideoftheMATLABconsole,inatexteditorsuchasviorEmacsortheMATLABeditor,orinan
integrateddevelopmentenvironment(IDE)suchasMicrosoftVisualStudio.
[Link]
command"mex",[Link],ifnecessary,tomovethecompiledmexfileintoyour
MATLABworkingdirectoryorpath.
Amexfile,asonecompilationunit,[Link],thenameofthe
compiledmexfunctionistheMATLABfunctionnameusedtoaccessit.(Thatnametypicallydoesn'tevenappear
inthesourcecode,unlessyouhappentoincludeitinacomment!)
Ofcourse,withinyourmexfileyoucandoanynumberofC++functioncalls,memoryallocationsor
deallocations,inputoutput,andsoforth,beforereturningcontrol(eitherwithorwithoutreturnvalues)tothe
[Link]'s
workspace,independentofwhethertheyareinputargumentstoyourmexfunction.

Hello,world!UsingtheM
[Link]

CAPI
1/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

Hereissourcecodeforamexfunctionthatreceivesaninputarrayofarbitrarysizeandshapeandreturnsthe
[Link],yes,italsoprints"Hello,world!"totheMATLABconsole.
/* [Link] */
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
int i, n = mxGetNumberOfElements(prhs[0]);
double *indata = (double *)mxGetData(prhs[0]);
plhs[0] = mxCreateNumericMatrix(1,n,mxDOUBLE_CLASS,mxREAL);
double *outdata = (double *)mxGetData(plhs[0]);
for (i=0;i<n;i++) outdata[i] = indata[i]*indata[i];
printf("Hello, world!\n");
}
WecompileitfromtheMATLABconsole,
>> mex [Link]
>>
Anexampleofusingournewmexfunctionis
>> b = [1 2 3; 4 5 6]
b =
1
2
3
4
5
6
>> a = helloworld(b)
Hello, world!
a =
1
16
4
25
>>

DoJustOnce:IntroduceM

36

toYourCompiler

YoumightwonderhowwegotthemexcommandtocompileC++source,sinceMATLABcomesonlywithabuilt
in,ratherbraindead,[Link]++compileralreadyonyourmachine
and,justonce,[Link],g++andIntelC++arelikelycompilersfor
WindowseitherMicrosoftVisualC++orIntelC++.
Tobindthe"mex"commandtoacompiler,dosomethinglikethefollowing:
>> mex -setup
Please choose your compiler for building external interface (MEX) files:
Would you like mex to locate installed compilers [y]/n? y
Select a compiler:
[1] Intel C++ 9.1 (with Microsoft Visual C++ 2005 linker)
in C:\Program Files\Intel\Compiler\C++\9.1
[2] Lcc-win32 C 2.4.1 in C:\PROGRA~1\MATLAB\R2007A~1\sys\lcc
[3] Microsoft Visual C++ 2005 in C:\Program Files\Microsoft Visual Studio 8
[4] Microsoft Visual C++ .NET 2003
in C:\Program Files\Microsoft Visual Studio .NET 2003
[5] Microsoft Visual C++ 6.0 in C:\Program Files\Microsoft Visual Studio
[0] None
Compiler: 3
Please verify your choices:
Compiler: Microsoft Visual C++ 2005
Location: C:\Program Files\Microsoft Visual Studio 8
Are these correct?([y]/n): y
Trying to update options file:
C:\Documents and Settings\...\MathWorks\MATLAB\R2007a\[Link]
From template: C:\PROGRA~1\MATLAB\R2007A~1\bin\win32\mexopts\[Link]

[Link]

2/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

Done . . .
>>
[Link].Options1or
3abovewouldbegoodchoices.Option2,thebraindeadCcompiler,wouldnotbe.

TheNR3CodingMethod
[Link]
[Link],[Link].
MexsourcefilesintheNR3codingmethodalwayshaveexactlytheframework
#include "nr3matlab.h"
/* you may put other stuff here */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
/* you must put something here */
}
ThevoidmexFunction(...)lineissomewhatequivalenttotheintmain(...)inordinaryC++
programming.(Youdon'tincludeamainroutineinmexfiles.)
Withinyourmexcode,youaregiventheintegernumberofarguments,nrhs,andtheexpectednumberof
[Link],we'llseehowtouse
theinputargumentsplhsandprhswiththeNR3codingmethod.
Althoughtheorderofpresentationmightatfirstseemodd,we'lldiscussvectors,thenscalars,andfinally
[Link],beforeanyofthese,wemustofcourserepeatthe"Hello,world!"example,nowinNR3coding
style:

Hello,world!Usingnr3matlab.h
UsingtheNR3codingmethod,our"Hello,world!"mexfunctionisrathermorereadablethanbefore.
/* [Link] */
#include "nr3matlab.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
MatDoub indata(prhs[0]);
Int i, j, k=0, m=[Link](), n=[Link]();
VecDoub outdata(m*n,plhs[0]);
for (i=0;i<m;i++) for (j=0;j<n;j++) outdata[k++] = SQR(indata[i][j]);
printf("Hello, world!\n");
}
Executionisnodifferentfrombefore.
>> mex
>> b =
b =
1
4
>> a =
Hello,
a =
1

[Link]
[1 2 3; 4 5 6]
2
3
5
6
helloworld2(b)
world!
16

AccessingM

25

36

Vectorswithnr3matlab.h

[Link]
[Link]
inanaturalway,withoutneedingtousetheAPImxfunctions.(Thatis,naturalforC++programmerswhohave
acopyofNumericalRecipes!)
Forreference,theadditionstothetemplatedclassNRvectorare:
template <class T> class NRvector {
/* ... */
NRvector(const mxArray *prhs);
// map Matlab rhs to vector
[Link]

3/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

NRvector(int n, mxArray* &plhs); // create Matlab lhs and map to vector

NRvector(const char *varname);
// map Matlab variable by name
void put(const char *varname);
// copy NRvector to a named Matlab variable
}
Youdon'treallyhavetounderstandtheabovedeclarations,ortheirdefinitions,[Link]
examples:
/* [Link] */
#include "nr3matlab.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
Int i,n;
// vector rhs arguments and lhs return values
VecDoub b(prhs[0]);
//
n = [Link]();
//
VecDoub c(n,plhs[0]); //
for (i=0;i<n;i++) c[i] =

bind the first rhs argument

get its number of elements
create a lhs return value of same length
SQR(b[i]);

// get and put values for named variables in MATLAB space

VecDoub d("dd"); // bind the variable dd (throws error if not exist)
n = [Link]();
// get its number of elements
VecDoub e(2*n); // make a local vector twice as long
for (i=0;i<n;i++) e[2*i] = e[2*i+1] = d[i];
// duplicate each value
[Link]("ee");
// copy e to MATLAB variable ee
// advanced usage (OK, but dangerous! see section "Living Dangerously")
b[0] = 999.;
d[0] = 999.;

// modify a rhs argument in MATLAB space

// modify a named variable directly (not using put())

}
See,nomxfunctions!Except,ofcourse,mexFunctionitself,[Link]
exampleshowsonlyVecDoubexemplars,everythingwouldworkfineforothertypesofvectors,e.g.,[Link]
constructorscheckfortypeconsistencyasnecessary,andreportanerrorifyouhavedonesomethingtype
[Link]'sseeitwork:
>> bb = [1 2 3]
bb =
1
2
3
>> dd = [4 5 6]
dd =
4
5
6
>> cc = NRvectorDemo(bb)
cc =
1
4
9
>> ee
ee =
4
4
5
5
>> bb
bb =
999
2
3
>> dd
dd =
999
5
6
>>

[Link],flatteningittoavectorinthe
MATLABstorageorder("bycolumns"or"Fortranlike"or"reverseofodometer").So,forexample,

>> bb = [1 2; 3 4]
bb =
1
2
3
4
>> dd = [5 6; 7 8]
dd =
[Link]

4/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

5
7

>>
cc
>>
bb
>>
dd

6
8
cc = NRvectorDemo(bb)
=
1
9
4
16
bb
=
999
2
3
4
dd
=
999
6
7
8

>>
Noticethatthe"dangerous"operations,becausetheywritedirectlyintoMATLAB'sworkspace,don'tflattenthe
matricestovectors,sincethematricesstayinMATLABspace.

LivingDangerously
Officially,wehavetosay,"don'tdoit!".ButsincethereasonyouareprogammingC++withinMATLABis
probablyforperformanceonbigdatasets,youwillprobablywanttousetheoperationsthatwelabelaboveas
"dangerous".Weusethemallthetime.
WhatisdangerousabouttheseoperationsisnotthattheywillcrashMATLABorvoidyourwarranty(iftherewere
one,alwaysdoubtfulforsoftware).[Link]
youunderstandthese,youcanmanagethembutifyoudon't,youmightberathersurprised:
>> bb = [1 2 3]
bb =
1
2
3
>> bbsave = bb
bbsave =
1
2
3
>> dd = [4 5 6]
dd =
4
5
6
>> cc = NRvectorDemo(bb)
cc =
1
4
9
>> bb
bb =
999
2
3
>> bbsave
bbsave =
999
2
3
>>
KnowingthatbbwasgoingtobemodifiedbyNRvectorDemo,[Link]:itgot
[Link]'sclevermemorymanagementusesa"copyonmodify"[Link]
youmodifiedbbinsideMATLAB,itwouldhavecopiedbbsave'[Link]
a"dangerous"method,undetectablebyMATLAB,thesavedcopynevergotmade.
Moral:Ifyouaregoingtousenr3matlab.h'sabilitytomodifyvaluesdirectlyinMATLABspace,besurethatthere
[Link],youcanenforceuniquenessbytrickery,
suchas
bbsave = bb;
bb = bb .* 0.5; bb = bb .* 2;
butthisshouldalmostneverbenecessaryifyoujustkeepinmindthecopyhistoryofyourlargedatavariables.
So,shouldn'tthisscareyouawayfromlivingdangerously?[Link]
length108,andyouintendtoapplyaseriesofC++mexfunctionstosomespecificcomponentsofthevector.
[Link]
[Link]()method,[Link],thenthe
"dangerous"operationsarejustwhatyouneed.

LivingSafely
[Link]

5/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

IfyouworrythatyoumightaccidentallywriteintoMATLABspace,youshoulddeclarevectorsderivedfromrhs
[Link],[Link]
linesfromtheexampleaboveare:
const VecDoub b(prhs[0]);
// bind a rhs argument
const VecDoub d("dd");
// bind the variable dd
b[0] = 999.;
/* caught as compile time error */
d[0] = 999.;
/* caught as compile time error */
IfyoulaterdecidethatyoudowanttowriteintoMATLABspaceafterall,youcandosomethinglikethis:
VecDoub &bx = const_cast<VecDoub&>(b);
bx[0] = 999.;
// modifies b, no compile time error

AccessingM

Scalarswithnr3matlab.h

TurnnowtoMATLABscalars.
Sincescalarslikedoubleandintarefundamentaltypes,wecan'toverloadfunctionalityintotheirconstructors,
[Link],[Link]
mxScalararedefined.(Notethat,despitethenamingconvention,thisisaNumericalRecipes,nota
MathWorks,function.)Thesefunctionsletyougrabscalarvalues,eitherfromyourmexfunction'sarguments,
[Link]++variabletoanyofyourfunction'sreturn
values,andcopyscalarvaluesintonamedvariables(existingornew)inMATLAB'sworkspace.
Forreference,thefunctionsaredeclaredasfollows:
template
template
template
template

<class
<class
<class
<class

T>
T>
T>
T>

const T& mxScalar(const mxArray *prhs);

T& mxScalar(mxArray* &plhs);
const T& mxScalar(const char *varname);
void mxScalar(T val, const char *varname);

Easiertounderstandissamplecode.(NotethatwearenowprogramminginNR3style,usingtypedef'dtypes
[Link]'thavetodothisifyoudon'twantto.)
/* [Link] */
#include "nr3matlab.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
// normal usage of mxScalar:
// get a rhs argument by value
Doub a = mxScalar<Doub>(prhs[0]);
// bind a reference to a lhs returned value
Doub &b = mxScalar<Doub>(plhs[0]);
b = SQR(a); // set a return value, e.g.
// get a named variable in MATLAB space by value
Doub c = mxScalar<Doub>("cc");
// set a named variable in MATLAB space
Doub d = SQR(c); // e.g.
mxScalar<Doub>(d,"dd");
// abnormal usage (OK but dangerous! -- see "Living Dangerously")
// get a non-const reference to a rhs argument
Doub &e = const_cast<Doub&>( mxScalar<Doub>(prhs[1]) );
e = SQR(e); // overwrite its value in MATLAB space
// get a non-const reference to a variable by name
Doub &f = const_cast<Doub&>( mxScalar<Doub>("ff") );
f = SQR(f); // overwrite its value in MATLAB space
}
Theactionofthisfunctionistoreturnthesquareofitsfirstargumentandsetthenamedvariableddtobethe
[Link]"dangerous"operations:Itchangesthevalueofitssecond
[Link]

6/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

righthandsideargumenttoitssquare,anditoverwritesanonconstreferencetothenamedvariableffwith
[Link]'swatchitgo:
>> aa=2; bb=3; cc=4; dd=5; ee=6; ff=7;
>> [aa bb cc dd ee ff]
ans =
2
3
4
5
6
7
>> bb = mxScalarDemo(aa,ee)
bb =
4
>> [aa bb cc dd ee ff]
ans =
2
4
4
16
36
49
>>
Youcanverifythatthe"after"valuesarewhatweexpect.
AlthoughtheaboveexampleshowsonlyDoubscalars,everythingwouldworkfineforothertypes(e.g.,Int).The
mxScalarfuctionscheckfortypeconsistencyasnecessary,andreportanerrorifyouhavedonesomethingtype
inconsistent.
NoticethatthemxScalarfuctionsaretemplated,andalwaysrequireanexplicittemplateargumentlike
mxScalar<double>()ormxScalar<int>().ThisissothattypecheckingagainstMATLAB'stypecanbemade
automatic,withanerrorthrownifyoutrytodosomethingbad.
AlsonotethatmostMATLABscalarsaredouble,[Link]
youwilloftenfindyourselfimportingargumentslikethis:
int n = int(mxScalar<Doub>(prhs[0]));

Matrices:ThroughtheLookingGlass(TTLG)
Matricesarehandledprettymuchlikevectors,above,exceptforoneflyintheointment,thefactthatMATLAB
andC++[Link]
[Link],oryouwillberatherunhappy:yourmatriceswillsometimesbethe
transposeofwhatyouexpected!NRmatrixdatatypesincludeNR3typeslikeMatDoub,MatInt,etc.
Thegoalisperformance:Wewanttobeabletoworkwithpotentiallylargeamountsofdatafromboththe
[Link]
backandforth,notonlybecauseoftheactualcopying,butalsobecauseoftheoverheadinvolvedinrapidly
allocatinganddeallocatinglargechunksofmemory(thuspossiblydrivingMATLAB'smemorymanagerintopoor
performance).
Therefore,wewanttopointdirectlyintoMATLAB'smemory,bothtoreaddataandalso(withsomerestrictions)
[Link],wesawthatthisposednospecialproblemsexceptunderstandingthe
implicationsof"dangerous"operationsandavoidingthem,ornot,[Link]
additionalproblemofstorageorder:MATLABishardwiredtostorebycolumns("Fortranorder"or"reverseof
odometer"),whileC++ishardwiredtostorebyrows("Corder"or"odometerorder").
Weadopta"ThroughtheLookingGlass"(TTLG)approach:Inourinterface,whenaMATLABmatrixcrossesthe
interfaceintoC++,[Link](oranyotherC++matrix)issent
backacrosstheinterfacetoMATLAB,[Link],forexample,a
matrixthatmakesaroundtripintoC++andbackagain(withoutbeingmodifiedbyC++code)isunchanged.
[Link]
[Link]
thisonthefly:Whereyoumighthavewrittenmymatrix[i][j],youinsteadwritemymatrix[j][i],andso
[Link],suchasdoingsomethinginparalleltoalltheelementsofamatrix,thereisnothingto
[Link],whenyouabsolutelyneedtoprocessanuntransposedmatrix,youhavetotake
[Link].
ItisimportanttounderstandthateachsideoftheTTLGinterfaceisselfconsistentaccordingtoitsownusual
conventions,bothastosizeandastosubscriptnumbering(1basedinMATLABvs.0basedinC++).So,if
mymatrixis3x5inMATLAB(3rows,5columns),andthus5x3inC++(5rows,3columns),thenthefollowing
quantitiesarenumericallyequal:
size(mymatrix,1)=3=[Link]()
size(mymatrix,2)=5=[Link]()
mymatrix(1,1)=mymatrix[0][0]
mymatrix(2,4)=mymatrix[3][1]
[Link]

7/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

AccessingM

Matriceswithnr3matlab.h

Forreference,theadditionstothetemplatedclassNRmatrixare:
template <class T> class NRmatrix {
/* ... */
NRmatrix(const mxArray *prhs); // map Matlab rhs to matrix
NRmatrix(int n, int m, mxArray* &plhs); // create Matlab lhs and map to matrix
NRmatrix(const char *varname); // import Matlab variable by name
void put(const char *varname); // copy NRmatrix to a named Matlab variable
}
Easiertounderstandaresomeusageexamples:
/* [Link] */
#include "nr3matlab.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
Int i,j,m,n;
// matrix rhs arguments and lhs return values
MatDoub b(prhs[0]);
// bind a rhs argument
m = [Link]();
// get shape
n = [Link]();
MatDoub c(m,n,plhs[0]); // create a lhs return value of same shape
for (i=0;im;i++) for (j=0;jn;j++) c[i][j] = SQR(b[i][j]);
// get and put values for named variables in MATLAB space
MatDoub d("dd"); // bind the variable dd
m = [Link]();
// get shape
n = [Link]();
MatDoub e(m,2*n); // make a local vector with twice as many cols
// duplicate each value (duplicating cols in C++, rows in Matlab!)
for (i=0;im;i++) for (j=0;jn;j++) e[i][2*j] = e[i][2*j+1] = d[i][j];
[Link]("ee");
// copy e to MATLAB variable ee
// advanced usage (OK, but dangerous! see section "Living Dangerously")
MatDoub
MatDoub
b[0][1]
d[1][0]

&bx = const_castMatDoub&>(b);
&dx = const_castMatDoub&>(d);
= 999.;
// modify a rhs argument in MATLAB space
= 999.;
// modify a named variable with no copying

}
FromtheMATLABside,wecanusethemexfunctionNRmatrixDemolikethis.
>> mex [Link]
>> bb = [1 2 3; 4 5 6]
bb =
1
2
3
4
5
6
>> dd = [7 8 9; 10 11 12]
dd =
7
8
9
10
11
12
>> cc = NRmatrixDemo(bb)
cc =
1
4
9
16
25
36
>> ee
ee =
7
8
9
7
8
9
10
11
12
10
11
12
>> bb
bb =
1
2
3
999
5
6
[Link]

8/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

>> dd
dd =
7
10

999
11

9
12

Noticethatthemodifiedelementisb[0][1]ontheC++side,butb(2,1)ontheMATLABside,becauseofthe
TTLGapproach.

WrapperFunctionsThatAccessMultipleMemberFunctionsinaClass
Ifyouwantasinglemexfiletoperformmultiplefunctions(forexample,calldifferentmemberfunctionswithin
asingleclass),you'llneedtofigureoutontheC++sidewhatisintendedfromthenumberandtypeofright
[Link].
Forexample,asimplewrapperfortheNumericalRecipesclassGaumixmod(Gaussianmixturemodel)might
havethesepossiblecallsontheMATLABside:
>> gmm('construct',data,means)

% construct the model from data and means

>> loglike = gmm('step',nsteps)

% step the model and return log-likelihood

>> [mean sig] = gmm(k) % return the mean and covariance of the kth component
>> resp = gmm('response') % return the response matrix
>> gmm('delete') % delete the model
ThecorrespondingC++[Link]()andmxT<>()totestthetypeof
MATLABarguments.Thesefunctionsareinnr3matlab.h.
/* [Link] */
#include "nr3matlab.h"
#include "cholesky.h"
#include "gaumixmod.h"
Gaumixmod *gmm = NULL;
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
int i,j,nn,kk,mm;
if (gmm) {nn=gmm->nn; kk=gmm->kk; mm=gmm->mm;}
if (gmm && nrhs == 1 && mxT(prhs[0]) == mxT<Doub>()) {
// [mean sig] = gmm(k)
Int k = Int(mxScalar<Doub>(prhs[0]));
if (nlhs > 0) {
VecDoub mean(mm,plhs[0]);
for (i=0;i<mm;i++) mean[i] = gmm->means[k-1][i];
}
if (nlhs > 1) {
MatDoub sig(mm,mm,plhs[1]);
for (i=0;i<mm;i++) for (j=0;j<mm;j++) sig[i][j] = gmm->sig[k-1][i][j];
}
} else if (nrhs == 1 && mxScalar<char>(prhs[0]) == 'd') {
// gmm('delete')
delete gmm;
} else if (gmm && nrhs == 1 && mxScalar<char>(prhs[0]) == 'r') {
// gmm('response')
if (nlhs > 0) {
MatDoub resp(nn,kk,plhs[0]);
for (i=0;i<nn;i++) for (j=0;j<kk;j++) resp[i][j] = gmm->resp[i][j];
}
} else if (gmm && nrhs == 2 && mxT(prhs[1]) == mxT<Doub>()) {
// deltaloglike = gmm('step',nsteps)
Int nstep = Int(mxScalar<Doub>(prhs[1]));
Doub tmp;
for (i=0;i<nstep;i++) {
tmp = gmm->estep();
gmm->mstep();
}
if (nlhs > 0) {
Doub &deltaloglike = mxScalar<Doub>(plhs[0]);
deltaloglike = tmp;
}
[Link]

9/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

} else if (nrhs == 3 && mxT(prhs[0]) == mxT<char>()) {

// gmm('construct',data,means)
MatDoub data(prhs[1]), means(prhs[2]);
if ([Link]() != [Link]()) throw("wrong dims in gmm 1");
if ([Link]() >= [Link]()) throw("wrong dims in gmm 2");
if (gmm) delete gmm;
gmm = new Gaumixmod(data,means);
} else {
throw("bad call to gmm");
}
return;
}
Notethat,onceinstantiated,thepointer*[Link]'dneed
amorecomplicatedschemetoinstantiatemorethanoneGaumixmodobjectatatime.

TheAPICodingMethod
Youneedtoreadthefollowingsectionsonlyifyouwanttogobeyondthecapabilitiesprovidedbynr3matlab.h
anduseMATLABCAPIcallsdirectly.

TipsonM

CAPIProgramming

[Link]
pageisthechapteronexternalinterfaces,[Link],wewillmakeafew
additionalpointshere.
Mexsourcefilesalwayshaveexactlytheframework
#include "mex.h"
/* you may put other stuff here */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
/* you must put something here */
}
ThevoidmexFunction(...)lineissomewhatequivalenttotheintmain(...)inordinaryC++
programming.(Youdon'tincludeamainroutineinmexfiles.)
Withinyourmexcode,youaregiventhenumberofarguments,nrhs,theexpectednumberofreturned
valuesnlhs,whichmayofcoursebearrays,anarrayofpointerstotheinputstructures(mxArrays),and
anarrayofpointersthatyoupopulatewithoutputmxArrays,asshownintheexample.
AlthoughdatainMATLABis(asthenameimplies)organizedintosubscriptedmatricesandarrays,theC
APIalwayspresentsMATLABdataasaonedimensionalarray(above,indata).Theorderingisalwaysby
columns,notrows,or,moregenerally,withthelastsubscriptchangingleastrapidly("reverseof
odometer"or"Fortranorder").ItiseasytofindoutMATLAB'sviewoftheinternalmatrixarrangementby
APIcallslikemxGetM,mxGetN,mxGetNumberOfDimensions,mxGetDimensions,[Link],
onceyoudothis,itisuptoyoutolocatebysubscript(s)anydesiredelements.(Thisisbasicallywhatthe
nr3matlabinterfaceisdesignedtomakeeasier.)MATLAB'sFortranstorageorderinghasotherimplications
thatarediscussedinthesectionMatrices:"ThroughtheLookingGlass",above.
[Link]++
[Link],wetrustedtheusertosendusdoubledata,whichisMATLAB'sdefaulttype(evenforthings
thatyoumightotherwiseexpecttobeintegers).Thatis,ofcourse,[Link]
usedAPIfunctionslikemxIsDoubleormxGetClassIDtogetMATLAB'sopinionastothedatatype,and
thenabortedwithmexErrMsgTxtifitisnottheexpectedenumvaluemxDOUBLE_CLASS.
IfyouplantousetheAPIinterface,[Link]/extern/include
[Link],thisfiledefinestheclassnameslike
mxDOUBLE_CLASS.

NumericalRecipesUsingtheM

CAPI

AlthoughweprefertousetheNR3codingmethod,above,youcanperfectlywelluseNumericalRecipesThird
Edition(NR3)code,[Link],withtheMATLABCAPIinterface.
Forexample,hereisamexfilethatmakesNR3'sgeneralizedFermiDiracintegralroutineavailabletoMATLAB.
(TheNR3routinedoesafancyquadrature,usingtheDErule.)

[Link]

10/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

/* [Link] */
#include "mex.h"
#include
#include
#include
#include

"nr3.h" // see [Link]

"quadrature.h"
"derule.h"
"fermi.h"

Fermi *fermi = NULL;

// usage: f = fermidirac(k, eta, theta)
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
double k = mxGetScalar(prhs[0]);
double eta = mxGetScalar(prhs[1]);
double theta = mxGetScalar(prhs[2]);
double f;
if (fermi == NULL) fermi = new Fermi();
if (theta < 0. || k <= -1.) {
delete fermi;
mexErrMsgTxt("bad args to fermidirac (also kludge destructor)");
} else {
f = fermi->val(k,eta,theta);
}
plhs[0] = mxCreateDoubleScalar(f);
}
Compilethiswiththemexcommand,thenputitinyourMATLABworkingdirectory.
>> f = fermidirac(3/2, 0.7, 0.2)
f =
2.3622
NowyoucanalsodopowerfulMATLABthings,likemakingacontourplotofthelogarithmofthefunction.
% fermidirac demo
[eta theta] = meshgrid(-2:.2:2,0:.25:3.5);
f = arrayfun(@(e,t) fermidirac(5/2,e,t), eta, theta);
flog = log10(f);
[C,h] = contour(eta,theta,flog,-0.2:.1:1.7);
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
colormap cool
(seeMATLABhelpfor"contour"tounderstandthis)whichproduces

[Link]

11/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

Becausethisexamplepassesandreturnsscalarsonly,[Link]
theNR3codingmethod:mxScalarwouldreplacemxGetScalar,throwwouldreplacemexErrMsgTxt,and
mxCreateDoubleScalarwoulddisappear(sincefcouldbedeclareddirectlyasthereturnvalueplhs[0]).

Appendix:UsingMicrosoftVisualStudio
Inallsectionsabovewe'vedonethecompilationandlinkingfromtheMATLABconsole,usingthemexcommand.
WesawinthesectionDoJustOnce:IntroduceMatlabtoYourCompilerthatwecouldcausethemexcommand
[Link]
developmentenvironment(IDE)[Link]
MATLAB7.1+[Link].
Wefollowtheusefularticleby"abn9c".WesupposethatyourfunctioniscalledMyFunc(maketheobvious
changesbelowforwhateveritisactuallynamed).
Followthesesteps:
[Link]++[Link]
type,chooseDLL.
[Link]/AddExistingItem...[Link],locatedinyourMatlabinstall
directory(e.g.,atC:\ProgramFiles\MATLAB\R2007b\extern\include).You'llneedtoshowfilesoftype
"ResourceFiles"toseethisfile.
[Link]
LIBRARY MyFunc.mexw32
EXPORTS mexFunction
[Link]/New/File...,choosing"HeaderFile.h"
[Link],[Link],[Link]
viamenuProject/AddExistingItem...
NowgotothemenuProject/MyFuncPropertiesforthefollowingsteps:
[Link]:AllConfigurations
[Link]/General/ConfigurationType,[Link].
[Link]

12/13

11/17/2016

ExtendingMATLABwithNumericalRecipes

[Link]/C++/General/AdditionalIncludeDirectories,addtheMATLAB
extern\includedirectory(e.g.,C:\ProgramFiles\MATLAB\R2007b\extern\include).
[Link]/C++/Preprocessor/PreprocessorDefinitions,addMATLAB_MEX_FILE
[Link]/Linker/General/OutputFile,[Link]
change$(OutDir)tothefullpathofyourMATLABworkingdirectory.(Thiswillsaveyoufromhavingto
recopytheoutputfileeverytimeyourecompile.)
[Link]/Linker/General/AdditionalLibraryDependencies,add
extern\lib\win32\microsoft(e.g.,C:\ProgramFiles\MATLAB\R2007b\extern\lib\win32\microsoft).
[Link]/Linker/Input/[Link],[Link],and
[Link]'tneedpaths,justthenames.
[Link]/Linker/Input/[Link](thatis,[Link]
filethatyoualreadycreated).Youdon'tneedapath,justthefilename.
Youcannow"OK"outof"MyFuncPropertyPages".
[Link],youcannowerasewhateveris"helpfully"provided,andenteryourmex
[Link]
#include "stdafx.h"
#include "nr3matlab.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
printf("Hello, world!\n"); // remove this line after testing
}
[Link].
[Link]/BuildSolution,[Link],movethefile
[Link]:
>> MyFunc()
Hello, world!
IfyouexecuteMyFuncfromwithinMATLABbetweencompilations,youmightneedtodo
>> clear MyFunc
beforethelinkercanoverwriteanewfile,orbeforeyournewexecutableisrecognizedbyMATLAB.
MicrosoftVisualStudiodoesn'[Link]
arevariouscommercialtoolsforthis,sothatyoudon'thavetokeepgoingthroughalltheabovestepsforeach
newmexfunctionthatyouwrite.(WeuseonecalledCopyWiz.)

[Link]

13/13