levmar.c

A sparse variant of the Levenberg-Marquardt algorithm implemented by levmar has been applied to bund

/* ////////////////////////////////////////////////////////////////////////////////
// 
//  Matlab MEX file for the Levenberg - Marquardt minimization algorithm
//  Copyright (C) 2007  Manolis Lourakis (lourakis at ics forth gr)
//  Institute of Computer Science, Foundation for Research & Technology - Hellas
//  Heraklion, Crete, Greece.
//
//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//////////////////////////////////////////////////////////////////////////////// */

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>
#include <string.h>
#include <ctype.h>

#include <lm.h>

#include <mex.h>

/**
#define DEBUG
**/

#ifndef HAVE_LAPACK
#ifdef _MSC_VER
#pragma message("LAPACK not available, certain functionalities cannot be compiled!")
#else
#warning LAPACK not available, certain functionalities cannot be compiled
#endif /* _MSC_VER */
#endif /* HAVE_LAPACK */

#define __MAX__(A, B)     ((A)>=(B)? (A) : (B))

#define MIN_UNCONSTRAINED     0
#define MIN_CONSTRAINED_BC    1
#define MIN_CONSTRAINED_LEC   2
#define MIN_CONSTRAINED_BLEC  3

struct mexdata {
  /* matlab names of the fitting function & its Jacobian */
  char *fname, *jacname;

  /* binary flags specifying if input p0 is a row or column vector */
  int isrow_p0;

  /* rhs args to be passed to matlab. rhs[0] is reserved for
   * passing the parameter vector. If present, problem-specific
   * data are passed in rhs[1], rhs[2], etc
   */
  mxArray **rhs;
  int nrhs; /* >= 1 */
};

/* display printf-style error messages in matlab */
static void matlabFmtdErrMsgTxt(char *fmt, ...)
{
char  buf[256];
va_list args;

	va_start(args, fmt);
	vsprintf(buf, fmt, args);
	va_end(args);

  mexErrMsgTxt(buf);
}

/* display printf-style warning messages in matlab */
static void matlabFmtdWarnMsgTxt(char *fmt, ...)
{
char  buf[256];
va_list args;

	va_start(args, fmt);
	vsprintf(buf, fmt, args);
	va_end(args);

  mexWarnMsgTxt(buf);
}

static void func(double *p, double *hx, int m, int n, void *adata)
{
mxArray *lhs[1];
double *mp, *mx;
register int i;
struct mexdata *dat=(struct mexdata *)adata;
    
  /* prepare to call matlab */
  mp=mxGetPr(dat->rhs[0]);
  for(i=0; i<m; ++i)
    mp[i]=p[i];
    
  /* invoke matlab */
  mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->fname);

  /* copy back results & cleanup */
  mx=mxGetPr(lhs[0]);
  for(i=0; i<n; ++i)
    hx[i]=mx[i];

  /* delete the matrix created by matlab */
  mxDestroyArray(lhs[0]);
}

static void jacfunc(double *p, double *j, int m, int n, void *adata)
{
mxArray *lhs[1];
double *mp;
double *mj;
register int i, k;
struct mexdata *dat=(struct mexdata *)adata;
    
  /* prepare to call matlab */
  mp=mxGetPr(dat->rhs[0]);
  for(i=0; i<m; ++i)
    mp[i]=p[i];

  /* invoke matlab */
  mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->jacname);
    
  /* copy back results & cleanup. Note that the nxm Jacobian 
   * computed by matlab should be transposed so that
   * levmar gets it in row major, as expected
   */
  mj=mxGetPr(lhs[0]);
  for(i=0; i<n; ++i)
    for(k=0; k<m; ++k)
      j[i*m+k]=mj[i+k*n];

  /* delete the matrix created by matlab */
  mxDestroyArray(lhs[0]);
}

/* matlab matrices are in column-major, this routine converts them to row major for levmar */
static double *getTranspose(mxArray *Am)
{
int m, n;
double *At, *A;
register int i, j;

  m=mxGetM(Am);
  n=mxGetN(Am);
  A=mxGetPr(Am);
  At=mxMalloc(m*n*sizeof(double));

  for(i=0; i<m; i++)
    for(j=0; j<n; j++)
      At[i*n+j]=A[i+j*m];
  
  return At;
}

/* check the supplied matlab function and its Jacobian. Returns 1 on error, 0 otherwise */
static int checkFuncAndJacobian(double *p, int  m, int n, int havejac, struct mexdata *dat)
{
mxArray *lhs[1];
register int i;
int ret=0;
double *mp;

  mexSetTrapFlag(1); /* handle errors in the MEX-file */

  mp=mxGetPr(dat->rhs[0]);
  for(i=0; i<m; ++i)
    mp[i]=p[i];

  /* attempt to call the supplied func */
  i=mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->fname);
  if(i){
    fprintf(stderr, "levmar: error calling '%s'.\n", dat->fname);
    ret=1;
  }
  else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || !(mxGetM(lhs[0])==1 || mxGetN(lhs[0])==1) ||
      __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0]))!=n){
    fprintf(stderr, "levmar: '%s' should produce a real vector with %d elements (got %d).\n",
                    dat->fname, n, __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0])));
    ret=1;
  }
  /* delete the matrix created by matlab */
  mxDestroyArray(lhs[0]);

  if(havejac){
    /* attempt to call the supplied jac  */
    i=mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->jacname);
    if(i){
      fprintf(stderr, "levmar: error calling '%s'.\n", dat->jacname);
      ret=1;
    }
    else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || mxGetM(lhs[0])!=n || mxGetN(lhs[0])!=m){
      fprintf(stderr, "levmar: '%s' should produce a real %dx%d matrix (got %dx%d).\n",
                      dat->jacname, n, m, mxGetM(lhs[0]), mxGetN(lhs[0]));
      ret=1;
    }
    else if(mxIsSparse(lhs[0])){
      fprintf(stderr, "levmar: '%s' should produce a real dense matrix (got a sparse one).\n", dat->jacname);
      ret=1;
    }
    /* delete the matrix created by matlab */
    mxDestroyArray(lhs[0]);
  }

  mexSetTrapFlag(0); /* on error terminate the MEX-file and return control to the MATLAB prompt */

  return ret;
}


/*
[ret, p, info, covar]=levmar_der (f, j, p0, x, itmax, opts, 'unc'                        ...)
[ret, p, info, covar]=levmar_bc  (f, j, p0, x, itmax, opts, 'bc',   lb, ub,              ...)
[ret, p, info, covar]=levmar_lec (f, j, p0, x, itmax, opts, 'lec',          A, b,        ...)
[ret, p, info, covar]=levmar_blec(f, j, p0, x, itmax, opts, 'blec', lb, ub, A, b, wghts, ...)
*/

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *Prhs[])
{
register int i;
register double *pdbl;
mxArray **prhs=(mxArray **)&Prhs[0], *At;
struct mexdata mdata;
int len, status;
double *p, *p0, *ret, *x;
int m, n, havejac, Arows, itmax, nopts, mintype, nextra;
double opts[LM_OPTS_SZ]={LM_INIT_MU, LM_STOP_THRESH, LM_STOP_THRESH, LM_STOP_THRESH, LM_DIFF_DELTA};
double info[LM_INFO_SZ];
double *lb=NULL, *ub=NULL, *A=NULL, *b=NULL, *wghts=NULL, *covar=NULL;

  /* parse input args; start by checking their number */
  if((nrhs<5))
    matlabFmtdErrMsgTxt("levmar: at least 5 input arguments required (got %d).", nrhs);
  if(nlhs>4)
    matlabFmtdErrMsgTxt("levmar: too many output arguments (max. 4, got %d).", nlhs);
  else if(nlhs<2)
    matlabFmtdErrMsgTxt("levmar: too few output arguments (min. 2, got %d).", nlhs);
    
  /* note that in order to accommodate optional args, prhs & nrhs are adjusted accordingly below */

  /** func **/
  /* first argument must be a string , i.e. a char row vector */
  if(mxIsChar(prhs[0])!=1)
    mexErrMsgTxt("levmar: first argument must be a string.");
  if(mxGetM(prhs[0])!=1)
    mexErrMsgTxt("levmar: first argument must be a string (i.e. char row vector).");
  /* store supplied name */
  len=mxGetN(prhs[0])+1;
  mdata.fname=mxCalloc(len, sizeof(char));
  status=mxGetString(prhs[0], mdata.fname, len);
  if(status!=0)
    mexErrMsgTxt("levmar: not enough space. String is truncated.");

  /** jac (optional) **/
  /* check whether second argument is a string */
  if(mxIsChar(prhs[1])==1){
    if(mxGetM(prhs[1])!=1)
      mexErrMsgTxt("levmar: second argument must be a string (i.e. row vector).");
    /* store supplied name */
    len=mxGetN(prhs[1])+1;
    mdata.jacname=mxCalloc(len, sizeof(char));
    status=mxGetString(prhs[1], mdata.jacname, len);
    if(status!=0)
      mexErrMsgTxt("levmar: not enough space. String is truncated.");
    havejac=1;

    ++prhs;
    --nrhs;
  }
  else{
    mdata.jacname=NULL;
    havejac=0;
  }

#ifdef DEBUG
  fflush(stderr);
  fprintf(stderr, "LEVMAR: %s analytic Jacobian\n", havejac? "with" : "no");
#endif /* DEBUG */

/* CHECK 
if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGetN(prhs[1])==1))
*/

  /** p0 **/