Skip to content

Commit

Permalink
Fix matlab interface (#392)
Browse files Browse the repository at this point in the history 
* added a wrapper from mxArray to vector

* let travis wait for brew install

* using auto where possible in matlab interface

* Fix(matlab) Prevent memory leaks due to ExpData initialization errors.
  • Loading branch information
@paulstapor
paulstapor authored Aug 8, 2018
1 parent 1fc0374 commit 7153063
Show file tree
Hide file tree
Showing 3 changed files with 32 additions and 30 deletions.
6 changes: 3 additions & 3 deletions .travis.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ branches:

matrix:
fast_finish: true
include:
include:
- os: linux
dist: xenial
compiler: gcc
Expand Down Expand Up @@ -36,10 +36,10 @@ matrix:
- os: osx
osx_image: xcode9.3
compiler: clang
before_install:
before_install:
- brew update # without this homebrew can stumble over wrong ruby version
- travis_wait brew install gcc || brew link --overwrite gcc # fix linker warning regarding /usr/local/include/c++
- brew install hdf5 cppcheck swig doxygen ragel graphviz homebrew/cask/mactex
- travis_wait brew install hdf5 cppcheck swig doxygen ragel graphviz homebrew/cask/mactex
- brew upgrade python
after_success:
- cd $BASE_DIR # cd to base dir for correct relative path in deploy
Expand Down
2 changes: 1 addition & 1 deletion include/amici/interface_matlab.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ ReturnDataMatlab *setupReturnData(mxArray *plhs[], int nlhs);
* @param[in] prhs user input @type *mxArray
* @return edata: experimental data struct @type *ExpData
*/
ExpData *expDataFromMatlabCall(const mxArray *prhs[], const Model &model);
std::unique_ptr<ExpData> expDataFromMatlabCall(const mxArray *prhs[], const Model &model);

void amici_dgemv(AMICI_BLAS_LAYOUT layout, AMICI_BLAS_TRANSPOSE TransA,
const int M, const int N, const double alpha, const double *A,
Expand Down
54 changes: 28 additions & 26 deletions src/interface_matlab.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,14 @@ void amici_daxpy(int n, double alpha, const double *x, const int incx, double *y
FORTRAN_WRAPPER(daxpy)(&n_, &alpha, x, &incx_, y, &incy_);
}

/** conversion from mxArray to vector<realtype>
* @param array Matlab array to create vector from
* @param length Number of elements in array
* @return std::vector with data from array
*/
std::vector<realtype> mxArrayToVector(const mxArray *array, int length) {
return {mxGetPr(array), mxGetPr(array) + length};
}

/*!
* expDataFromMatlabCall parses the experimental data from the matlab call and
Expand All @@ -168,54 +176,45 @@ void amici_daxpy(int n, double alpha, const double *x, const int incx, double *y
* dimension checks
* @return edata pointer to experimental data object
*/
ExpData *expDataFromMatlabCall(const mxArray *prhs[],
std::unique_ptr<ExpData> expDataFromMatlabCall(const mxArray *prhs[],
Model const &model) {
if (!mxGetPr(prhs[RHS_DATA]))
return nullptr;

int nt_my = 0, ny_my = 0, nt_sigmay = 0,
ny_sigmay = 0; /* integers with problem dimensionality */
int ne_mz = 0, nz_mz = 0, ne_sigmaz = 0,
nz_sigmaz = 0; /* integers with problem dimensionality */
ExpData *edata = new ExpData(model);
if (edata->my.empty() && edata->mz.empty()) {
// if my and mz are both empty, no (or empty) data was provided
// in that case we simply return a nullptr for easier checking.
delete(edata);
return nullptr;
}

auto edata = std::unique_ptr<ExpData>(new ExpData(model));

// Y
if (mxArray *dataY = mxGetProperty(prhs[RHS_DATA], 0, "Y")) {
ny_my = (int)mxGetN(dataY);
auto ny_my = static_cast<int>(mxGetN(dataY));
if (ny_my != model.nytrue) {
throw AmiException("Number of observables in data matrix (%i) does "
"not match model ny (%i)",
ny_my, model.nytrue);
}
nt_my = (int)mxGetM(dataY);
auto nt_my = static_cast<int>(mxGetM(dataY));
if (nt_my != model.nt()) {
throw AmiException("Number of time-points in data matrix does (%i) "
"not match provided time vector (%i)",
nt_my, model.nt());
}
mxArray *dataYT;
mexCallMATLAB(1, &dataYT, 1, &dataY, "transpose");
edata->setObservedData(mxGetPr(dataYT));
auto observedData = mxArrayToVector(dataYT, ny_my * nt_my);
edata->setObservedData(observedData);

} else {
throw AmiException("Field Y not specified as field in data struct!");
}

// Sigma Y
if (mxArray *dataSigmaY = mxGetProperty(prhs[RHS_DATA], 0, "Sigma_Y")) {
ny_sigmay = (int)mxGetN(dataSigmaY);
auto ny_sigmay = static_cast<int>(mxGetN(dataSigmaY));
if (ny_sigmay != model.nytrue) {
throw AmiException("Number of observables in data-sigma matrix (%i) "
"does not match model ny (%i)",
ny_sigmay, model.nytrue);
}
nt_sigmay = (int)mxGetM(dataSigmaY);
auto nt_sigmay = static_cast<int>(mxGetM(dataSigmaY));
if (nt_sigmay != model.nt()) {
throw AmiException("Number of time-points in data-sigma matrix (%i) "
"does not match provided time vector (%i)",
Expand All @@ -224,49 +223,52 @@ ExpData *expDataFromMatlabCall(const mxArray *prhs[],

mxArray *dataSigmaYT;
mexCallMATLAB(1, &dataSigmaYT, 1, &dataSigmaY, "transpose");
edata->setObservedDataStdDev(mxGetPr(dataSigmaYT));
auto observedDataSigma = mxArrayToVector(dataSigmaYT, ny_sigmay * nt_sigmay);
edata->setObservedDataStdDev(observedDataSigma);
} else {
throw AmiException("Field Sigma_Y not specified as field in data struct!");
}

// Z
if (mxArray *dataZ = mxGetProperty(prhs[RHS_DATA], 0, "Z")) {
nz_mz = (int)mxGetN(dataZ);
auto nz_mz = static_cast<int>(mxGetN(dataZ));
if (nz_mz != model.nztrue) {
throw AmiException("Number of events in event matrix (%i) does not "
"match provided nz (%i)",
nz_mz, model.nztrue);
}
ne_mz = (int)mxGetM(dataZ);
auto ne_mz = static_cast<int>(mxGetM(dataZ));
if (ne_mz != model.nMaxEvent()) {
throw AmiException("Number of time-points in event matrix (%i) does "
"not match provided nmaxevent (%i)",
ne_mz, model.nMaxEvent());
}
mxArray *dataZT;
mexCallMATLAB(1, &dataZT, 1, &dataZ, "transpose");
edata->setObservedEvents(mxGetPr(dataZT));
auto observedEvents = mxArrayToVector(dataZT, nz_mz * ne_mz);
edata->setObservedEvents(observedEvents);
} else {
throw AmiException("Field Z not specified as field in data struct!");
}

// Sigma Z
if (mxArray *dataSigmaZ = mxGetProperty(prhs[RHS_DATA], 0, "Sigma_Z")) {
nz_sigmaz = (int)mxGetN(dataSigmaZ);
auto nz_sigmaz = static_cast<int>(mxGetN(dataSigmaZ));
if (nz_sigmaz != model.nztrue) {
throw AmiException("Number of events in event-sigma matrix (%i) does "
"not match provided nz (%i)",
nz_sigmaz, model.nztrue);
}
ne_sigmaz = (int)mxGetM(dataSigmaZ);
auto ne_sigmaz = static_cast<int>(mxGetM(dataSigmaZ));
if (ne_sigmaz != model.nMaxEvent()) {
throw AmiException("Number of time-points in event-sigma matrix (%i) "
"does not match provided nmaxevent (%i)",
ne_sigmaz, model.nMaxEvent());
}
mxArray *dataSigmaZT;
mexCallMATLAB(1, &dataSigmaZT, 1, &dataSigmaZ, "transpose");
edata->setObservedEventsStdDev(mxGetPr(dataSigmaZT));
auto observedEventsSigma = mxArrayToVector(dataSigmaZT, nz_sigmaz * ne_sigmaz);
edata->setObservedEventsStdDev(observedEventsSigma);
} else {
throw AmiException("Field Sigma_Z not specified as field in data struct!");

Expand Down Expand Up @@ -531,7 +533,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
std::unique_ptr<amici::ExpData> edata;
if (nrhs > amici::RHS_DATA && mxGetPr(prhs[amici::RHS_DATA])) {
try {
edata.reset(amici::expDataFromMatlabCall(prhs, *model));
edata = std::move(amici::expDataFromMatlabCall(prhs, *model));
} catch (amici::AmiException& ex) {
amici::errMsgIdAndTxt("AMICI:mex:setup","Failed to read experimental data:\n%s",ex.what());
}
Expand Down

0 comments on commit 7153063

Please sign in to comment.